Regresión lineal y agrupar por en R

Question

Quiero hacer una regresión lineal en R usando la función lm(). Mis datos son series temporales anuales con un campo por año (22 años) y otro por estado (50 estados). Quiero ajustar una regresión para cada estado para que al final tenga un vector de respuestas de lm. Puedo imaginar hacer un ciclo para cada estado, luego hacer la regresión dentro del ciclo y agregar los resultados de cada regresión a un vector. Eso no parece muy parecido a R, sin embargo. 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?

Answer 1

Aquí hay una forma de usar el paquete lme4.

> 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 

Answer 2

## 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 
> 

Answer 3

En mi opinión, es un modelo lineal mixto 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 para 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)) 

Answer 4

He aquí un método que utiliza 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) 

Answer 5

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

años var1 var2 Estado Y ...

En mi punto de vista, se 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) 

Answer 6

Una buena solución usando data.table se publicó here en CrossValidated by @Zach. que acababa de añadir que es posible obtener de forma iterativa 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 salida 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 

Answer 7

Desde el año 2009, tiene dplyr ha sido lanzado, que en realidad ofrece una manera muy agradable de hacer este tipo de agrupamiento, muy parecido a lo que SAS hace.

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/Valor PD, se puede usar el paquete broom. Este paquete proporciona:

tres genéricos S3: ordenado, donde se resuman resultados estadísticos de un modelo como coeficientes de una regresión; aumento, que agrega columnas a los datos originales, como predicciones, residuos y asignaciones de clúster; y vistazo, que proporciona un resumen de una sola 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) 

Answer 8

ahora mi respuesta llega un poco tarde, pero yo estaba buscando una funcionalidad similar.Al parecer, la función integrada 'por' en R también puede hacer el agrupar fácilmente:?

mediante 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) 

Answer 9

El la pregunta parece ser acerca de cómo llamar a las funciones de regresión con fórmulas que se modifican dentro de un bucle.

Aquí es cómo puede hacerlo en (el uso de diamantes conjunto de datos):

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() 
    }