Lineare Regression und Gruppierung nach in R.

92

Ich möchte mit der lm()Funktion eine lineare Regression in R durchführen . Meine Daten sind eine jährliche Zeitreihe mit einem Feld für das Jahr (22 Jahre) und einem anderen für den Staat (50 Staaten). Ich möchte für jeden Zustand eine Regression anpassen, damit ich am Ende einen Vektor von lm-Antworten habe. Ich kann mir vorstellen, für jeden Zustand eine for-Schleife auszuführen, dann die Regression innerhalb der Schleife durchzuführen und die Ergebnisse jeder Regression einem Vektor hinzuzufügen. Das scheint jedoch nicht sehr R-artig zu sein. In SAS würde ich eine 'by'-Anweisung machen und in SQL würde ich eine' group by 'machen. Was ist der R-Weg, um das zu tun?

JD Long
quelle
1
Ich möchte den Leuten nur sagen, dass es in R zwar viele Group-by-Funktionen gibt, aber nicht alle die richtige für die Group-by-Regression sind. Zum Beispiel aggregateist nicht richtig ; auch nichttapply .
18.

Antworten:

46

Hier ist eine Möglichkeit, das lme4Paket 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
ars
quelle
1
Gibt es eine Möglichkeit, R2 für diese beiden Modelle aufzulisten? Fügen Sie beispielsweise nach einem Jahr eine R2-Spalte hinzu. Addieren Sie auch den p-Wert für jeden der Koeffizienten?
ToToRo
3
@ToToRo hier finden Sie eine einfache Lösung (besser spät als nie): Your.df [, Zusammenfassung (lm (Y ~ X)) $ r.squared, by = Your.factor] wobei: Y, X und Your.factor sind deine Variablen. Bitte denken Sie daran, dass Your.df eine data.table-Klasse sein muss
FraNut
59

Hier ist ein Ansatz mit dem Plyr- Paket:

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)
Hadley
quelle
Angenommen, Sie haben eine zusätzliche unabhängige Variable hinzugefügt, die nicht in allen Bundesstaaten verfügbar war (z. B. Miles.of.ocean.shoreline), die von NA in Ihren Daten dargestellt wurde. Würde der lm-Anruf nicht scheitern? Wie könnte damit umgegangen werden?
MikeTP
Innerhalb der Funktion müssten Sie für diesen Fall testen und eine andere Formel verwenden
Hadley
Ist es möglich, jedem Aufruf in der Zusammenfassung den Namen der Untergruppe hinzuzufügen (letzter Schritt)?
Rote Beete
Wenn Sie laufen layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page und dann erhalten l_ply(models, plot)Sie auch jedes der Residuen-Diagramme. Ist es möglich, jedes der Diagramme mit der Gruppe zu kennzeichnen (z. B. "Zustand" in diesem Fall)?
Brian D
48

Seit 2009 dplyrwurde veröffentlicht, was tatsächlich eine sehr gute Möglichkeit bietet, diese Art der Gruppierung durchzuführen, die der von SAS sehr ähnlich ist.

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  

Um die Koeffizienten und Rsquared / p.value abzurufen, kann man das broomPaket verwenden. Dieses Paket bietet:

drei S3-Generika: ordentlich, das die statistischen Ergebnisse eines Modells wie die Koeffizienten einer Regression zusammenfasst; Augment, das den Originaldaten Spalten wie Vorhersagen, Residuen und Clusterzuweisungen hinzufügt; und Blick, der eine einzeilige Zusammenfassung der Statistiken auf Modellebene bietet.

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)
Paul Hiemstra
quelle
2
Ich musste tun rowwise(fitted_models) %>% tidy(model), um das Besenpaket zum Laufen zu bringen, aber ansonsten eine großartige Antwort.
Pedram
3
Funktioniert großartig ... kann dies alles tun, ohne die Pfeife zu verlassen:d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)
Holastello
23

Meiner Meinung nach ist ein gemischtes lineares Modell ein besserer Ansatz für diese Art von Daten. Der unten angegebene Code im festen Effekt gibt den Gesamttrend wieder. Die zufälligen Effekte zeigen an, 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 (Modelle mit gemischten Effekten in S und S-Plus) an.

library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
Thierry
quelle
3
Dies ist eine wirklich gute Antwort auf die allgemeine Statistik-Theorie, die mich über einige Dinge nachdenken lässt, die ich nicht berücksichtigt habe. Die Anwendung, die mich veranlasst hat, die Frage zu stellen, ist für diese Lösung nicht anwendbar, aber ich bin froh, dass Sie sie angesprochen haben. Danke dir.
JD Long
1
Es ist keine gute Idee, mit einem gemischten Modell zu beginnen. Woher wissen Sie, dass eine der Annahmen gerechtfertigt ist?
Hadley
7
Man sollte die Annahme durch Modellvalidierung (und Kenntnis der Daten) überprüfen. Übrigens können Sie die Annahme auch nicht für die einzelnen Filme garantieren. Sie müssten alle Modelle separat validieren.
Thierry
14

Eine nette Lösung mit data.tablewurde hier in CrossValidated von @Zach gepostet. Ich möchte nur hinzufügen, dass es möglich ist, iterativ auch den Regressionskoeffizienten r ^ 2 zu erhalten:

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

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
FraNut
quelle
8

Ich denke, es lohnt sich, den purrr::mapAnsatz für dieses Problem hinzuzufügen .

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

Weitere Ideen zur Verwendung des broomPakets mit diesen Ergebnissen finden Sie in der Antwort von @Paul Hiemstra .

ngm
quelle
Eine kleine Erweiterung für den Fall, dass Sie eine Spalte mit angepassten Werten oder Residuen wünschen: Wickeln Sie den Aufruf von lm () in einen Aufruf von resid () ein und leiten Sie dann alles in der letzten Zeile in einen Aufruf von unnest () um. Natürlich möchten Sie den Variablennamen von "Modell" in etwas Relevanteres ändern.
Randy
8
## 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
Eduardo Leoni
quelle
6

Ich jetzt meine Antwort kommt etwas spät, aber ich suchte nach einer ähnlichen Funktionalität. Es scheint, dass die eingebaute Funktion 'by' in R auch die Gruppierung einfach durchführen kann:

? by enthält das folgende Beispiel, das pro Gruppe passt und die Koeffizienten mit sapply extrahiert:

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)
Matthijs Cox
quelle
3

Die lm()obige Funktion ist ein einfaches Beispiel. Ich stelle mir übrigens vor, dass Ihre Datenbank die Spalten wie folgt enthält:

Jahreszustand var1 var2 y ...

Aus meiner Sicht können Sie den folgenden Code verwenden:

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)
Zack Mendes
quelle
0

Die Frage scheint zu sein, wie Regressionsfunktionen mit Formeln aufgerufen werden, die innerhalb einer Schleife geändert werden.

So können Sie dies tun (mithilfe des Diamanten-Datensatzes):

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()
  }
IVIM
quelle