Регрессия по подмножеству в R

Question

Я новичок в R и пытаюсь запустить линейную регрессию на нескольких подмножествах («случаев») данных в одном файле. У меня есть 50 разных случаев, поэтому я не хочу запускать 50 различных регрессий ... будьте добры, чтобы автоматизировать это. Я нашел и экспериментировал с методом ddply, но по какой-то причине он возвращает мне одинаковые коэффициенты для каждого случая. Код, я использую выглядит следующим образом:

ddply(MyData, "Case", function(x) coefficients(lm(Y~X1+X2+X3, MyData)))

Результаты я, опять же, одни и те же коэффициенты для каждого «Case». Любые идеи о том, как я могу улучшить свой код, чтобы регрессия выполнялась один раз для каждого случая и давала мне уникальные коэффициенты для каждого случая?

Answer 1

ddply передает данные. Кадров (от разделения входных данных.фрейма) к функции. Вы, наверное, хотите:

ddply(MyData, "Case", function(df) coefficients(lm(Y~X1+X2+X3, data=df))) 

(. Не тестировался, так как вы не обеспечивают воспроизводимый пример)

Вы прошли весь входной data.frame в lm для каждой группы ..

Answer 2

You может также сделать это, используя только базовые функции:

# load data 
data(warpbreaks) 

# fit linear model to each subset 
fits <- by(warpbreaks, warpbreaks[,"tension"], 
      function(x) lm(breaks ~ wool, data = x)) 

# Combine coefficients from each model 
do.call("rbind", lapply(fits, coef)) 

Answer 3

Хотя вы можете получить довольно далеко с (г) plyr, я получаю большую гибкость с простой цикл (как я не слишком конфи вмятина с делать(), когда речь идет, желая больше, чем просто выход коэффициентов(), например)

Начнем с загрузки данных и загрузки пакеты:

library(dplyr) 
library(broom) # get lm output (coefficients) as a dataframe 
library(reshape2) # Melt/DCAST for swapping between wide/long format 
data(warpbreaks) 

Далее создайте список групп, которые мы собираемся перебирать.

GROUPS <- unique(warpbreaks$tension) # Specify groups, one unique model per group. 

код, чтобы ответить на ваш вопрос, в данном случае, было бы что-то из рода:

for (i in 1:length(GROUPS)){ 
    CURRENT_GROUP <- GROUPS[i] 
    df <- filter(warpbreaks, tension==CURRENT_GROUP) # subset the dataframe 
    fit <- lm(breaks ~ wool, data = df) # Build a model 
    coeff <- tidy(fit) # Get a pretty data frame of the coefficients & p values 
    coeff <- coeff[,c(1,2,5)] # Extract P.Value & Estimate 

    # Rename (intercept) to INT for pretty column names 
    coeff[coeff$term=="(Intercept)", ]$term <- "INT" 

    # Make it into wide format with reshape2 package. 
    coeff <- coeff %>% melt(id.vars=c("term")) 

    # Defactor the resulting data.frame 
    coeff <- mutate(coeff, 
        variable=as.character(variable), 
        term=as.character(term)) 

    # Rename for prettier column names later 
    coeff[coeff$variable=="estimate", ]$variable <- "Beta" 
    coeff[coeff$variable=="p.value", ]$variable <- "P" 

    coeff <- dcast(coeff, .~term+variable)[,-1] 

    rsquared <- summary(fit)$r.squared 

    # Create a df out of what we just did. 
    row <- cbind(
    data.frame(
     group=CURRENT_GROUP, 
     rsquared=rsquared), 
    coeff 
) 

    # If first iteration, create data.frame -- otherwise: rowbind 
    if (i==1){ 
    RESULT_ROW = row 
    } 
    else{ 
    RESULT_ROW = rbind(RESULT_ROW, row) 
    } # End if. 
} # End for loop 

При написании внутренней для кода петли, просто проверить что-то путем моделирования цикла (первый запустите i < - 1, запустите внутренний код цикла шаг за шагом, затем запустите i < - 2 и т. д.). Полученный dataframe:

RESULT_ROW 
## group rsquared INT_Beta  INT_P woolB_Beta woolB_P 
## 1  L 0.26107601 44.55556 9.010046e-08 -16.333333 0.03023435 
## 2  M 0.07263228 24.00000 5.991177e-07 4.777778 0.27947884 
## 3  H 0.12666292 24.55556 9.234773e-08 -5.777778 0.14719571 

Я не говорю, что это лучше, чем ответы, описанные выше, я так же, как тот факт, что это дает мне гибкость, чтобы извлечь что-либо из любой модели (не только ая модель, которые имеют симпатичные функции коэффициентов()).