Regresión lineal y agrupación en R

93

Quiero hacer una regresión lineal en R usando la lm()función. Mis datos son una serie de tiempo anual con un campo por año (22 años) y otro por estado (50 estados). Quiero ajustar una regresión para cada estado de modo que al final tenga un vector de respuestas lm. Puedo imaginarme haciendo el ciclo for para cada estado, luego haciendo la regresión dentro del ciclo y agregando los resultados de cada regresión a un vector. Sin embargo, eso no parece muy parecido a R. En SAS haría una declaración 'por' y en SQL haría un 'grupo por'. ¿Cuál es la forma R de hacer esto?

r regression linear-regression lm JD Long
fuente
1
Solo quiero decirle a la gente que aunque hay muchas funciones group-by en R, no todas son las adecuadas para la regresión group-by. Por ejemplo, aggregateno es correcto ; tampoco lo estapply .
李哲源

Respuestas:

48

Aquí hay una forma de usar el lme4paquete.

 library(lme4)
 d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)),
                 year=rep(1:10, 2),
                 response=c(rnorm(10), rnorm(10)))

 xyplot(response ~ year, groups=state, data=d, type='l')

 fits <- lmList(response ~ year | state, data=d)
 fits
#------------
Call: lmList(formula = response ~ year | state, data = d)
Coefficients:
   (Intercept)        year
CA -1.34420990  0.17139963
NY  0.00196176 -0.01852429

Degrees of freedom: 20 total; 16 residual
Residual standard error: 0.8201316
ars
fuente
2
¿Hay alguna forma de enumerar R2 para estos dos modelos? por ejemplo, agregue una columna R2 después del año. También agregue el valor p para cada uno de los coeficientes.
ToToRo
3
@ToToRo aquí puede encontrar una solución factible (más vale tarde que nunca): Your.df [, summary (lm (Y ~ X)) $ r.squared, by = Your.factor] donde: Y, X y Your.factor son tus variables. Tenga en cuenta que Your.df debe ser una clase
data.table
60

Aquí hay un enfoque usando el paquete plyr :

d <- data.frame(
  state = rep(c('NY', 'CA'), 10),
  year = rep(1:10, 2),
  response= rnorm(20)
)

library(plyr)
# Break up d by state, then fit the specified model to each piece and
# return a list
models <- dlply(d, "state", function(df) 
  lm(response ~ year, data = df))

# Apply coef to each model and return a data frame
ldply(models, coef)

# Print the summary of each model
l_ply(models, summary, .print = TRUE)
Hadley
fuente
Supongamos que agregó una variable independiente adicional que no estaba disponible en todos los estados (es decir, millas de costa oceánica) que estaba representada por NA en sus datos. ¿No fallaría la llamada? ¿Cómo podría solucionarse?
MikeTP
Dentro de la función, necesitaría probar ese caso y usar una fórmula diferente
hadley
¿Es posible agregar el nombre del subgrupo a cada llamada en el resumen (último paso)?
remolacha
si ejecuta layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page y luego l_ply(models, plot)obtiene cada uno de los gráficos de residuos también. ¿Es posible etiquetar cada una de las parcelas con el grupo (por ejemplo, "estado" en este caso)?
Brian D
51

Desde 2009, dplyrse ha lanzado, lo que en realidad proporciona una forma muy agradable de hacer este tipo de agrupación, muy similar a lo que hace SAS.

library(dplyr)

d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)),
                year=rep(1:10, 2),
                response=c(rnorm(10), rnorm(10)))
fitted_models = d %>% group_by(state) %>% do(model = lm(response ~ year, data = .))
# Source: local data frame [2 x 2]
# Groups: <by row>
#
#    state   model
#   (fctr)   (chr)
# 1     CA <S3:lm>
# 2     NY <S3:lm>
fitted_models$model
# [[1]]
# 
# Call:
# lm(formula = response ~ year, data = .)
# 
# Coefficients:
# (Intercept)         year  
#    -0.06354      0.02677  
#
#
# [[2]]
# 
# Call:
# lm(formula = response ~ year, data = .)
# 
# Coefficients:
# (Intercept)         year  
#    -0.35136      0.09385

Para recuperar los coeficientes y Rsquared / p.value, se puede usar el broompaquete. Este paquete proporciona:

tres genéricos S3: tidy, que resume los hallazgos estadísticos de un modelo, como los coeficientes de una regresión; aumento, que agrega columnas a los datos originales, como predicciones, residuos y asignaciones de grupos; y vistazo, que proporciona un resumen de una fila de estadísticas a nivel de modelo.

library(broom)
fitted_models %>% tidy(model)
# Source: local data frame [4 x 6]
# Groups: state [2]
# 
#    state        term    estimate  std.error  statistic   p.value
#   (fctr)       (chr)       (dbl)      (dbl)      (dbl)     (dbl)
# 1     CA (Intercept) -0.06354035 0.83863054 -0.0757668 0.9414651
# 2     CA        year  0.02677048 0.13515755  0.1980687 0.8479318
# 3     NY (Intercept) -0.35135766 0.60100314 -0.5846187 0.5749166
# 4     NY        year  0.09385309 0.09686043  0.9689519 0.3609470
fitted_models %>% glance(model)
# Source: local data frame [2 x 12]
# Groups: state [2]
# 
#    state   r.squared adj.r.squared     sigma statistic   p.value    df
#   (fctr)       (dbl)         (dbl)     (dbl)     (dbl)     (dbl) (int)
# 1     CA 0.004879969  -0.119510035 1.2276294 0.0392312 0.8479318     2
# 2     NY 0.105032068  -0.006838924 0.8797785 0.9388678 0.3609470     2
# Variables not shown: logLik (dbl), AIC (dbl), BIC (dbl), deviance (dbl),
#   df.residual (int)
fitted_models %>% augment(model)
# Source: local data frame [20 x 10]
# Groups: state [2]
# 
#     state   response  year      .fitted   .se.fit     .resid      .hat
#    (fctr)      (dbl) (int)        (dbl)     (dbl)      (dbl)     (dbl)
# 1      CA  0.4547765     1 -0.036769875 0.7215439  0.4915464 0.3454545
# 2      CA  0.1217003     2 -0.009999399 0.6119518  0.1316997 0.2484848
# 3      CA -0.6153836     3  0.016771076 0.5146646 -0.6321546 0.1757576
# 4      CA -0.9978060     4  0.043541551 0.4379605 -1.0413476 0.1272727
# 5      CA  2.1385614     5  0.070312027 0.3940486  2.0682494 0.1030303
# 6      CA -0.3924598     6  0.097082502 0.3940486 -0.4895423 0.1030303
# 7      CA -0.5918738     7  0.123852977 0.4379605 -0.7157268 0.1272727
# 8      CA  0.4671346     8  0.150623453 0.5146646  0.3165112 0.1757576
# 9      CA -1.4958726     9  0.177393928 0.6119518 -1.6732666 0.2484848
# 10     CA  1.7481956    10  0.204164404 0.7215439  1.5440312 0.3454545
# 11     NY -0.6285230     1 -0.257504572 0.5170932 -0.3710185 0.3454545
# 12     NY  1.0566099     2 -0.163651479 0.4385542  1.2202614 0.2484848
# 13     NY -0.5274693     3 -0.069798386 0.3688335 -0.4576709 0.1757576
# 14     NY  0.6097983     4  0.024054706 0.3138637  0.5857436 0.1272727
# 15     NY -1.5511940     5  0.117907799 0.2823942 -1.6691018 0.1030303
# 16     NY  0.7440243     6  0.211760892 0.2823942  0.5322634 0.1030303
# 17     NY  0.1054719     7  0.305613984 0.3138637 -0.2001421 0.1272727
# 18     NY  0.7513057     8  0.399467077 0.3688335  0.3518387 0.1757576
# 19     NY -0.1271655     9  0.493320170 0.4385542 -0.6204857 0.2484848
# 20     NY  1.2154852    10  0.587173262 0.5170932  0.6283119 0.3454545
# Variables not shown: .sigma (dbl), .cooksd (dbl), .std.resid (dbl)
Paul Hiemstra
fuente
2
Tuve que hacer rowwise(fitted_models) %>% tidy(model)para que el paquete de escobas funcionara, pero por lo demás, una gran respuesta.
pedram
3
Funciona muy bien ... puede hacer todo esto sin salir de la tubería:d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)
holastello
@pedram y @holastello, esto ya no funciona, al menos con R 3.6.1, broom_0.7.0, dplyr_0.8.3. d %>% group_by(state) %>% do(model=lm(response ~year, data = .)) %>% rowwise() %>% tidy(model) Error in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm) : Calling var(x) on a factor x is defunct. Use something like 'all(duplicated(x)[-1L])' to test for a constant vector. In addition: Warning messages: 1: Data frame tidiers are deprecated and will be removed in an upcoming release of broom. ...
Chris Nolte
23

En mi opinión, un modelo lineal mixto es un mejor enfoque para este tipo de datos. El código a continuación da en el efecto fijo la tendencia general. Los efectos aleatorios indican cómo la tendencia de cada estado individual difiere de la tendencia global. La estructura de correlación tiene en cuenta la autocorrelación temporal. Eche un vistazo a Pinheiro & Bates (modelos de efectos mixtos en S y S-Plus).

library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
Thierry
fuente
3
Esta es una muy buena respuesta a la teoría de estadísticas generales que me hace pensar en algunas cosas que no había considerado. La aplicación que me motivó a hacer la pregunta no sería aplicable a esta solución, pero me alegra que la hayas mencionado. Gracias.
JD Long
1
No es una buena idea comenzar con un modelo mixto, ¿cómo sabe que se justifica alguna de las suposiciones?
hadley
7
Se debe verificar la suposición mediante la validación del modelo (y el conocimiento de los datos). Por cierto, tampoco puede garantizar la suposición en las películas individuales. Tendría que validar todos los modelos por separado.
Thierry
14

Se data.tablepublicó una buena solución aquí en CrossValidated por @Zach. Solo agregaría que es posible obtener iterativamente también el coeficiente de regresión r ^ 2:

## make fake data
    library(data.table)
    set.seed(1)
    dat <- data.table(x=runif(100), y=runif(100), grp=rep(1:2,50))

##calculate the regression coefficient r^2
    dat[,summary(lm(y~x))$r.squared,by=grp]
       grp         V1
    1:   1 0.01465726
    2:   2 0.02256595

así como todos los demás resultados de summary(lm):

dat[,list(r2=summary(lm(y~x))$r.squared , f=summary(lm(y~x))$fstatistic[1] ),by=grp]
   grp         r2        f
1:   1 0.01465726 0.714014
2:   2 0.02256595 1.108173
FraNut
fuente
8

Creo que vale la pena agregar el purrr::mapenfoque a este problema.

library(tidyverse)

d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)),
                                 year=rep(1:10, 2),
                                 response=c(rnorm(10), rnorm(10)))

d %>% 
  group_by(state) %>% 
  nest() %>% 
  mutate(model = map(data, ~lm(response ~ year, data = .)))

Consulte la respuesta de @Paul Hiemstra para obtener más ideas sobre cómo usar el broompaquete con estos resultados.

ngm
fuente
Una pequeña extensión en caso de que desee una columna de valores ajustados o residuales: envuelva la llamada lm () en una llamada resid () y luego canalice todo en la última línea en una llamada unnest (). Por supuesto, querrá cambiar el nombre de la variable de "modelo" a algo más relevante.
randy
8 
## make fake data
 ngroups <- 2
 group <- 1:ngroups
 nobs <- 100
 dta <- data.frame(group=rep(group,each=nobs),y=rnorm(nobs*ngroups),x=runif(nobs*ngroups))
 head(dta)
#--------------------
  group          y         x
1     1  0.6482007 0.5429575
2     1 -0.4637118 0.7052843
3     1 -0.5129840 0.7312955
4     1 -0.6612649 0.9028034
5     1 -0.5197448 0.1661308
6     1  0.4240346 0.8944253
#------------ 
## function to extract the results of one model
 foo <- function(z) {
   ## coef and se in a data frame
   mr <- data.frame(coef(summary(lm(y~x,data=z))))
   ## put row names (predictors/indep variables)
   mr$predictor <- rownames(mr)
   mr
 }
 ## see that it works
 foo(subset(dta,group==1))
#=========
              Estimate Std..Error   t.value  Pr...t..   predictor
(Intercept)  0.2176477  0.1919140  1.134090 0.2595235 (Intercept)
x           -0.3669890  0.3321875 -1.104765 0.2719666           x
#----------
## one option: use command by
 res <- by(dta,dta$group,foo)
 res
#=========
dta$group: 1
              Estimate Std..Error   t.value  Pr...t..   predictor
(Intercept)  0.2176477  0.1919140  1.134090 0.2595235 (Intercept)
x           -0.3669890  0.3321875 -1.104765 0.2719666           x
------------------------------------------------------------ 
dta$group: 2
               Estimate Std..Error    t.value  Pr...t..   predictor
(Intercept) -0.04039422  0.1682335 -0.2401081 0.8107480 (Intercept)
x            0.06286456  0.3020321  0.2081387 0.8355526           x

## using package plyr is better
 library(plyr)
 res <- ddply(dta,"group",foo)
 res
#----------
  group    Estimate Std..Error    t.value  Pr...t..   predictor
1     1  0.21764767  0.1919140  1.1340897 0.2595235 (Intercept)
2     1 -0.36698898  0.3321875 -1.1047647 0.2719666           x
3     2 -0.04039422  0.1682335 -0.2401081 0.8107480 (Intercept)
4     2  0.06286456  0.3020321  0.2081387 0.8355526           x
Eduardo Leoni
fuente
6

Ahora mi respuesta llega un poco tarde, pero estaba buscando una funcionalidad similar. Parecería que la función incorporada 'por' en R también puede hacer la agrupación fácilmente:

? by contiene el siguiente ejemplo, que se ajusta por grupo y extrae los coeficientes con sapply:

require(stats)
## now suppose we want to extract the coefficients by group 
tmp <- with(warpbreaks,
            by(warpbreaks, tension,
               function(x) lm(breaks ~ wool, data = x)))
sapply(tmp, coef)
Matthijs Cox
fuente
3

La lm()función anterior es un ejemplo simple. Por cierto, imagino que su base de datos tiene las columnas de la siguiente forma:

año estado var1 var2 y ...

En mi punto de vista, puede utilizar el siguiente código:

require(base) 
library(base) 
attach(data) # data = your data base
             #state is your label for the states column
modell<-by(data, data$state, function(data) lm(y~I(1/var1)+I(1/var2)))
summary(modell)
Zack Mendes
fuente
0

La pregunta parece ser sobre cómo llamar a funciones de regresión con fórmulas que se modifican dentro de un ciclo.

Así es como puede hacerlo (usando el conjunto de datos de diamantes):

attach(ggplot2::diamonds)
strCols = names(ggplot2::diamonds)

formula <- list(); model <- list()
for (i in 1:1) {
  formula[[i]] = paste0(strCols[7], " ~ ", strCols[7+i])
  model[[i]] = glm(formula[[i]]) 

  #then you can plot the results or anything else ...
  png(filename = sprintf("diamonds_price=glm(%s).png", strCols[7+i]))
  par(mfrow = c(2, 2))      
  plot(model[[i]])
  dev.off()
  }
IVIM
fuente