Regresión lineal y agrupar por en R

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 statement ‘por’ y en SQL haría un ‘grupo por’. ¿Cuál es la forma R de hacer esto?

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 

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) 

Desde 2009, se ha lanzado dplyr , que en realidad ofrece una manera 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:  # # state model # (fctr) (chr) # 1 CA  # 2 NY  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, uno puede usar el paquete de broom . Este paquete proporciona:

tres generics S3: ordenados, que resume los hallazgos estadísticos de un modelo como los coeficientes de una regresión; boost, 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) 

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

Una buena solución usando data.table fue publicada 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 del 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 
 ## 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 yx 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 > 

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:

? por 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) 

La función lm() arriba es un ejemplo simple. Por cierto, imagino que su base de datos tiene las columnas como en el siguiente formulario:

año estado var1 var2 y …

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

Creo que vale la pena agregar el enfoque purrr::map 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 el uso del paquete de broom con estos resultados.

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

Aquí es cómo puedes hacerlo en (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() } 
    Intereting Posts