Lineare Regression und group by in R

Möchte ich eine lineare regression in R mit der lm() Funktion. Meine Daten ist eine jährliche Zeitreihen mit einem Feld für das Jahr (22 Jahre) und eine weitere für den Staat (50 Staaten). Ich will fit eine regression für jeden Staat, so dass am Ende habe ich eine Vektor-lm-Antworten. Ich kann es mir vorstellen for-Schleife für jeden Staat, dann tut der regression innerhalb der Schleife und die addition der Ergebnisse der einzelnen regression zu einem Vektor. Das scheint nicht sehr R-wie, jedoch. In SAS würde ich tun, eine 'von' - Anweisung in SQL würde ich tun, eine 'group by'. Was ist der R Weg, dies zu tun?

InformationsquelleAutor der Frage JD Long | 2009-07-23

linear-regression lm r regression

Hier ist eine Möglichkeit, mit der lme4 Paket.

> 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

InformationsquelleAutor der Antwort ars

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)

InformationsquelleAutor der Antwort hadley

Seit 2009 dplyr veröffentlicht worden, die tatsächlich bietet eine sehr schöne Möglichkeit, das zu tun diese Art der Gruppierung sehr ähnlich, was SAS bedeutet.

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

Abrufen der Koeffizienten und Rsquared/p.Wert, man kann die broom Paket. Dieses Paket bietet:

drei S3 generics: ordentlich, in dem Sie beschreibt ein Modell ist
statistische Ergebnisse wie die Koeffizienten einer regression;
vermehren, die Spalten Hinzugefügt, um die ursprünglichen Daten wie
Vorhersagen, Residuen und die cluster-Zuordnungen; und der Blick, der
enthält eine Zeile Zusammenfassung der Modell-level-Statistiken.

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)

InformationsquelleAutor der Antwort Paul Hiemstra

22

Meiner Meinung nach ist eine gemischte lineare Modell ein besserer Ansatz für diese Art von Daten. Den code unten in den fixen Effekt der Allgemeine trend. Die zufälligen Effekte geben an, wie sich der trend für jeden einzelnen Zustand unterscheiden sich von dem globalen trend. Die Korrelation Struktur erfolgt die zeitliche Autokorrelation berücksichtigt. Haben Sie einen Blick auf Pinheiro & Bates (Mixed-Effects Models in S and S-Plus).
```
library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
```
InformationsquelleAutor der Antwort Thierry

Schöne Lösung mit data.table gepostet wurde hier in CrossValidated von @Zach.
Ich würde nur hinzufügen, dass es möglich ist, um iterativ auch der regressionskoeffizient 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

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

InformationsquelleAutor der Antwort FraNut

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

InformationsquelleAutor der Antwort Eduardo Leoni

4

Ich jetzt meine Antwort kommt ein bisschen spät, aber ich war auf der Suche für eine ähnliche Funktionalität. Es scheint, als würde die integrierte Funktion 'durch' in R können auch die Gruppierung leicht:

?durch enthält das folgende Beispiel, das passt pro Gruppe und extrahiert die Koeffizienten mit 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)
```
InformationsquelleAutor der Antwort Matthijs Cox
3

Den lm() obige Funktion ist ein einfaches Beispiel. Übrigens, ich kann mir vorstellen, dass Ihre Datenbank hat die Spalten wie in der folgenden form:

Jahr Zustand var1, var2 y...

In meiner Sicht, Sie können 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)
```
InformationsquelleAutor der Antwort Zack Mendes

Ich denke, es ist sinnvoll, fügen Sie die purrr::map Ansatz für dieses problem.

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

Siehe @Paul Hiemstra-Antwort für weitere Ideen zur Verwendung der broom Paket mit diesen Ergebnissen.

InformationsquelleAutor der Antwort ngm

Scheint die Frage zu sein, über wie nennen regression Funktionen mit Formeln, verändert sich innerhalb einer Schleife.

Hier ist, wie können Sie es tun (mit Diamanten dataset):

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

InformationsquelleAutor der Antwort IVIM

Schreibe einen Kommentar

Du musst angemeldet sein, um einen Kommentar abzugeben.