ddply с функцией lm () - PullRequest
Новая
       35

ddply с функцией lm ()

15 голосов
/ 23 сентября 2011

Привет, ребята, как я могу использовать функцию ddply для линейной модели: 

x1 <- c(1:10, 1:10)
x2 <- c(1:5, 1:5, 1:5, 1:5)
x3 <- c(rep(1,5), rep(2,5), rep(1,5), rep(2,5))

set.seed(123)
y <- rnorm(20, 10, 3)
mydf <- data.frame(x1, x2, x3, y)

require(plyr)
ddply(mydf, mydf$x3, .fun = lm(mydf$y ~ mydf$X1 + mydf$x2))

Генерирует эту ошибку:

Ошибка в model.frame.default (формула = mydf $ y~ mydf $ X1 + mydf $ x2, drop.unused.levels = TRUE): недопустимый тип (NULL) для переменной 'mydf $ X1'

Благодарим вас за помощь.

Ответы [ 3 ]

28 голосов
/ 23 сентября 2011

Вот что вам нужно сделать. 

mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)

моды - это список двух объектов, содержащих результаты регрессии.Вы можете извлечь то, что вам нужно из модов.Например, если вы хотите извлечь коэффициенты, вы можете написать 

coefs = ldply(mods, coef)

Это даст вам 

  x3 (Intercept)         x1 x2
1  1    11.71015 -0.3193146 NA
2  2    21.83969 -1.4677690 NA

РЕДАКТИРОВАТЬ.Если вы хотите ANOVA, то вы можете просто сделать 

ldply(mods, anova)

  x3 Df    Sum Sq   Mean Sq   F value     Pr(>F)
1  1  1  2.039237  2.039237 0.4450663 0.52345980
2  1  8 36.654982  4.581873        NA         NA
3  2  1 43.086916 43.086916 4.4273907 0.06849533
4  2  8 77.855187  9.731898        NA         NA
12 голосов
/ 23 сентября 2011

То, что объяснил Рамнатх, совершенно верно. Но я немного уточню.

ddply ожидает кадр данных и затем возвращает кадр данных. Функция lm() принимает фрейм данных в качестве входных данных, но возвращает объект линейной модели в ответ. Вы можете увидеть это, посмотрев документы для lm через ?lm:

Значение

lm возвращает объект класса "lm" или для нескольких ответов класса с ("млм", "лм").

Так что вы не можете просто поместить объекты lm во фрейм данных. Вы можете либо принудительно настроить вывод lm во фрейм данных, либо вы можете поместить объекты lm в список вместо фрейма данных.

Итак, чтобы проиллюстрировать оба варианта:

Вот как можно поместить объекты lm в список (очень похоже на то, что проиллюстрировал Рамнатх): 

outlist <- dlply(mydf, "x3", function(df)  lm(y ~ x1 + x2, data=df))

С другой стороны, если вы хотите извлечь только коэффициенты, вы можете создать функцию, которая запускает регрессию, а затем возвращает только коэффициенты в виде фрейма данных, например: 

myLm <- function( formula, df ){
  lmList <- lm(formula, data=df)
  lmOut <- data.frame(t(lmList$coefficients))
  names(lmOut) <- c("intercept","x1coef","x2coef")
  return(lmOut)
}

outDf <- ddply(mydf, "x3", function(df)  myLm(y ~ x1 + x2, df))
1 голос
/ 23 сентября 2011

Используйте это 

mods <- dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
coefs <- llply(mods, coef)

$`1`
(Intercept)          x1          x2 
 11.7101519  -0.3193146          NA 

$`2`
(Intercept)          x1          x2 
  21.839687   -1.467769          NA 



anovas <- llply(mods, anova)

$`1`
Analysis of Variance Table

Response: y
      Df Sum Sq Mean Sq F value Pr(>F)
x1         1  2.039  2.0392  0.4451 0.5235
Residuals  8 36.655  4.5819               

$`2`
Analysis of Variance Table

Response: y
      Df Sum Sq Mean Sq F value Pr(>F)  
x1         1 43.087  43.087  4.4274 0.0685 .
Residuals  8 77.855   9.732                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
...