superponer - Regresión lineal y agrupar por en R
superponer graficas en r (10)
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?
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)
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)
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
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.
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: <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, uno puede usar el paquete de broom . Este paquete proporciona:
tres genéricos S3: ordenados, que resume los hallazgos estadísticos de un modelo como los coeficientes de una regresión; aumentar, 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))
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)
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()
}
Una buena solución usando data.table fue publicada here 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 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
>