Looping над комбинациями терминов модели регрессии

Question

Я бег регресса в виде

reg=lm(y ~ x1+x2+x3+z1,data=mydata) 

На месте последнего срока, z1, я хочу проходной набор различных переменных, z1 через z10 , выполняющий регрессию для каждого из них с последним термином. Например. во втором прогоне я хочу использовать

reg=lm(y ~ x1+x2+x3+z2,data=mydata) 

в 3-м периоде:

reg=lm(y ~ x1+x2+x3+z3,data=mydata) 

Как можно автоматизировать с помощью цикла через список г-переменных?

Answer 1

С этой фиктивными данными:

dat1 <- data.frame(y = rpois(100,5), 
x1 = runif(100), 
x2 = runif(100), 
x3 = runif(100), 
z1 = runif(100), 
z2 = runif(100) 
) 

Вы можете получить список из двух lm объектов таким образом:

lapply(dat1[5:6], function(x) lm(dat1$y ~ dat1$x1 + dat1$x2 + dat1$x3 + x)) 

Какие итерацию через эти две колонки и заменяет их в качестве аргументов в lm вызова ,

Как отмечает Алекс, предпочтительнее передавать имена по формуле, а не фактические столбцы данных, как я сделал здесь.

Answer 2

В то время как Сэм предоставил работы и является хорошим решением, я лично предпочел бы пойти немного по-другому. Его ответ уже принят, поэтому я просто размещаю это ради полноты.

dat1 <- data.frame(y = rpois(100, 5), 
        x1 = runif(100), 
        x2 = runif(100), 
        x3 = runif(100), 
        z1 = runif(100), 
        z2 = runif(100)) 

lapply(colnames(dat1)[5:6], 
     function(x, d) lm(as.formula(paste("y ~ x1 + x2 + x3", x, sep = " + ")), data = d), 
     d = dat1) 

Вместо цикла по фактическим колонкам кадра данных, это только петли над строкой имен. Это обеспечивает некоторые улучшения скорости, поскольку количество копий между итерациями меньше.

library(microbenchmark) 

microbenchmark({ lapply(<what I wrote above>) }) 
# Unit: milliseconds 
# expr 
# {lapply(colnames(dat1)[5:6], function(x, d) lm(as.formula(paste("y ~ x1 + x2 + x3", x, sep = "+")), data = d), d = dat1)} 
#  min  lq  mean median  uq  max neval 
# 4.014237 4.148117 4.323387 4.220189 4.281995 5.898811 100 

microbenchmark({ lapply(<other answer>) }) 
# Unit: milliseconds 
# expr 
# {lapply(dat1[, 5:6], function(x) lm(dat1$y ~ dat1$x1 + dat1$x2 + dat1$x3 + x))} 
#  min  lq  mean median  uq max neval 
# 4.391494 4.505056 5.186972 4.598301 4.698818 51.573 100 

Разница достаточно мал для этого игрушечного примера, но как число наблюдений и предикторов увеличиваются, разница, вероятно, станет более заметной.

Answer 3

Вот другой подход, использующий пакеты от семейства dplyr/tidyr. Это перестраивает данные в длинной форме, а затем использует group_by() из dplyr пакета вместо lapply():

library(dplyr) 
library(tidyr) 
library(magrittr) # for use_series() 
dat1 %>% 
    gather(varname, z, z1:z2) %>% # convert data to long form 
    group_by(varname) %>% 
    do(model = lm(y ~ x1 + x2 + x3 + z, data = .)) %>% 
    use_series(model) 

Это преобразует данные в длинном формате с использованием gather, где Z-значения занимают один и тот же столбец. use_series() из пакета magrittr вернуть список lm объектов вместо data.frame.Если вы загрузите broom пакет, вы можете извлечь коэффициенты модели в этом трубопроводе кода:

library(broom) 
dat1 %>% 
    gather(varname, z, z1:z2) %>% 
    group_by(varname) %>% 
    do(model = lm(y ~ x1 + x2 + x3 + z, data = .)) %>% 
    glance(model) # or tidy(model) 

Source: local data frame [2 x 12] 
Groups: varname 

    varname r.squared adj.r.squared sigma statistic p.value df logLik  AIC  BIC deviance df.residual 
1  z1 0.06606736 0.02674388 2.075924 1.680099 0.1609905 5 -212.3698 436.7396 452.3707 409.3987   95 
2  z2 0.06518852 0.02582804 2.076900 1.656192 0.1666479 5 -212.4168 436.8337 452.4647 409.7840   95 

данных:

dat1 <- data.frame(y = rpois(100, 5), x1 = runif(100), 
        x2 = runif(100), x3 = runif(100), 
        z1 = runif(100), z2 = runif(100)) 

Answer 4

В зависимости от того, что ваша конечная цель, это может быть много быстрее, чтобы соответствовать базовой модели, обновить его с add1 и извлечь F-тест/AIC вы хотите:

> basemodel <- lm(y~x1+x2+x3, dat1) 
> 
> add1(object=basemodel, grep("z\\d", names(dat1), value=TRUE), test="F") 
Single term additions 

Model: 
y ~ x1 + x2 + x3 
     Df Sum of Sq RSS AIC F value Pr(>F) 
<none>    477.34 164.31    
z1  1 0.0768 477.26 166.29 0.0153 0.9019 
z2  1 5.1937 472.15 165.21 1.0450 0.3093 

См. Также ?update для определения модели.