Lineare Regression und Gruppierung nach in R

Question

Ich möchte eine lineare Regression in R mit der lm()-Funktion durchführen. Meine Daten sind eine jährliche Zeitreihe mit einem Feld für das Jahr (22 Jahre) und einem anderen für den Bundesstaat (50 Bundesstaaten). Ich möchte eine Regression für jeden Bundesstaat anpassen, sodass ich am Ende einen Vektor von lm-Antworten habe. Ich kann mir vorstellen, für jede Zustandsschleife eine Schleife zu machen, dann die Regression innerhalb der Schleife durchzuführen und die Ergebnisse jeder Regression zu einem Vektor hinzuzufügen. Das scheint jedoch nicht sehr R-typisch zu sein. In SAS würde ich eine 'by'-Anweisung verwenden und in SQL würde ich ein 'group by' machen. Wie funktioniert das in R?

Answer 1

Seit 2009 ist dplyr veröffentlicht worden, was tatsächlich eine sehr schöne Möglichkeit bietet, diese Art von Gruppierung durchzuführen, die SAS sehr ähnelt.

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 = .))
# Quelle: Lokaler Datenrahmen [2 x 2]
# Gruppen: 
#
#    state   model
#   (fctr)   (chr)
# 1     CA 
# 2     NY 
fitted_models$model
# [[1]]
# 
# Aufruf:
# lm(formula = response ~ year, data = .)
# 
# Koeffizienten:
# (Intercept)         year  
#    -0.06354      0.02677  
#
#
# [[2]]
# 
# Aufruf:
# lm(formula = response ~ year, data = .)
# 
# Koeffizienten:
# (Intercept)         year  
#    -0.35136      0.09385  

Um die Koeffizienten und Rsquared/p-Wert abzurufen, kann man das broom Paket verwenden. Dieses Paket bietet:

drei S3 Generika: tidy, das die statistischen Ergebnisse eines Modells zusammenfasst, wie zum Beispiel Koeffizienten einer Regression; augment, das Spalten zu den Originaldaten hinzufügt, wie Vorhersagen, Residuen und Clusterzuordnungen; und glance, das eine einzeilige Zusammenfassung der Modellstatistiken liefert.

library(broom)
fitted_models %>% tidy(model)
# Quelle: Lokaler Datenrahmen [4 x 6]
# Gruppen: 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)
# Quelle: Lokaler Datenrahmen [2 x 12]
# Gruppen: 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
# Variablen nicht gezeigt: logLik (dbl), AIC (dbl), BIC (dbl), deviance (dbl),
#   df.residual (int)
fitted_models %>% augment(model)
# Quelle: Lokaler Datenrahmen [20 x 10]
# Gruppen: 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
# Variablen nicht gezeigt: .sigma (dbl), .cooksd (dbl), .std.resid (dbl)

Answer 2

Hier ist ein Ansatz unter Verwendung des plyr Pakets:

d <- data.frame(
  state = rep(c('NY', 'CA'), 10),
  year = rep(1:10, 2),
  response= rnorm(20)
)

library(plyr)
# Zerlegen Sie d nach Bundesland, passen Sie dann das angegebene Modell für jedes Stück an und geben Sie eine Liste zurück
models <- dlply(d, "state", function(df) 
  lm(response ~ year, data = df))

# Wenden Sie coef auf jedes Modell an und geben Sie ein Datenrahmen zurück
ldply(models, coef)

# Drucken Sie die Zusammenfassung jedes Modells
l_ply(models, summary, .print = TRUE)

Answer 3

Hier ist eine Möglichkeit, das lme4-Paket zu verwenden.

 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 4

Nach meiner Meinung ist ein gemischtes lineares Modell ein besseres Verfahren für diese Art von Daten. Der folgende Code gibt den Gesamttrend im Fixeffekt an. Die Zufallseffekte zeigen, wie sich der Trend für jeden einzelnen Staat vom globalen Trend unterscheidet. Die Korrelationsstruktur berücksichtigt die zeitliche Autokorrelation. Schauen Sie sich Pinheiro & Bates (Gemischte Effekte Modelle in S und S-Plus) an.

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

Answer 5

Ein schöne Lösung mit data.table wurde hier in CrossValidated von @Zach veröffentlicht. Ich würde nur hinzufügen, dass es auch möglich ist, iterativ den Regressionskoeffizienten r^2 zu erhalten:

## erstelle gefälschte Daten
    library(data.table)
    set.seed(1)
    dat <- data.table(x=runif(100), y=runif(100), grp=rep(1:2,50))

## berechne den Regressionskoeffizienten r^2
    dat[,summary(lm(y~x))$r.squared,by=grp]
       grp         V1
    1:   1 0.01465726
    2:   2 0.02256595

sowie alle anderen Ausgaben von 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