Regressão Linear e Agrupar por em R

Eu quero fazer uma regressão linear em R usando a function lm() . Meus dados são uma série temporal anual com um campo por ano (22 anos) e outro por estado (50 estados). Eu quero ajustar uma regressão para cada estado para que no final eu tenha um vetor de respostas lm. Eu posso imaginar fazendo loop para cada estado, em seguida, fazendo a regressão dentro do loop e adicionando os resultados de cada regressão a um vetor. Isso não parece muito semelhante a R, no entanto. No SAS eu faria uma instrução ‘by’ e no SQL eu faria um ‘group by’. Qual é o jeito de fazer isso?

   

Aqui está uma maneira de usar o pacote 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 

Aqui está uma abordagem usando o pacote 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, o dplyr foi lançado, o que realmente fornece uma maneira muito legal de fazer esse tipo de agrupamento, parecido com o que o SAS faz.

 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 os coeficientes e Rsquared / p.value, pode-se usar o pacote broom . Este pacote fornece:

três genéricos S3: arrumado, que resume as descobertas statistics de um modelo, como coeficientes de uma regressão; augment, que adiciona colunas aos dados originais, como previsões, resíduos e atribuições de cluster; e glance, que fornece um resumo de uma linha de statistics no nível do 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) 

Na minha opinião, é um modelo linear misto uma abordagem melhor para este tipo de dados. O código abaixo dado no efeito fixo a tendência geral. Os efeitos randoms indicam como a tendência para cada estado individual difere da tendência global. A estrutura de correlação leva em consideração a autocorrelação temporal. Dê uma olhada no Pinheiro & Bates (Modelos de Efeitos Mistos em S e S-Plus).

 library(nlme) lme(response ~ year, random = ~year|state, correlation = corAR1(~year)) 

Uma boa solução usando data.table foi postada aqui em CrossValidated by @Zach. Gostaria apenas de acrescentar que é possível obter iterativamente também o coeficiente de regressão 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 

bem como todos os outros resultados do 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 > 

Agora minha resposta chega um pouco atrasada, mas eu estava procurando por uma funcionalidade semelhante. Parece que a function interna ‘por’ em R também pode fazer o agrupamento facilmente:

by contém o exemplo a seguir, que se ajusta por grupo e extrai os coeficientes com solidez:

 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) 

A function lm() acima é um exemplo simples. A propósito, imagino que seu database tenha as colunas da seguinte forma:

ano estado var1 var2 y …

No meu ponto de vista, você pode usar o seguinte 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) 

Eu acho que vale a pena adicionar a abordagem purrr::map para 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 = .))) 

Veja a resposta do @Paul Hiemstra para mais idéias sobre o uso do pacote de broom com estes resultados.

A questão parece ser sobre como chamar funções de regressão com fórmulas que são modificadas dentro de um loop.

Aqui está como você pode fazer isso (usando dataset 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() }