Régression linéaire et groupe par dans R

Question

Je wan pour faire une régression linéaire dans R à l'aide de l' lm() fonction. Mes données est un temps annuel de la série avec un champ pour l'année (22 ans) et un autre pour l'état (50 membres). Je veux ajustement d'une droite de régression pour chaque état, de sorte qu'à la fin j'ai un vecteur de réponses lm. Je peux imaginer faire de boucle pour chaque état ensuite de faire de la régression à l'intérieur de la boucle, et en additionnant les résultats de chaque régression à un vecteur. Cela ne semble pas très R-comme, cependant. Dans le SAS, je ferais une 'par' énoncé, et en SQL, je voudrais faire un 'group by'. Quel est le R moyen de faire cela?

Answer 1

Voici une approche utilisant le paquet 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 2

Voici un moyen d'utiliser le package 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 3

À mon avis, un modèle linéaire mixte est une meilleure approche pour ce type de données. Le code ci-dessous donne dans l'effet fixe la tendance globale. Les effets aléatoires indiquent en quoi la tendance de chaque État diffère de la tendance globale. La structure de corrélation prend en compte l'autocorrélation temporelle. Jetez un coup d’œil à Pinheiro & Bates (Modèles à effets mixtes en S et S-Plus).

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

Answer 4

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

La fonction lm() ci-dessus est un exemple simple. Au fait, j'imagine que votre base de données a les colonnes comme dans le formulaire suivant:

année état var1 var2 y ...

De mon point de vue, vous pouvez utiliser le code suivant:

 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)