Lineare Regression und Gruppe von in R

Question

Ich möchte eine lineare Regression in R unter Verwendung der lm() Funktion. Meine Daten sind eine jährliche Zeitreihe mit einem Feld für das Jahr (22 Jahre) und ein anderes für den Staat (50 Staaten). Ich möchte eine Regression für jeden Zustand anpassen, so dass ich am Ende einen Vektor von lm Antworten habe. Ich kann mir vorstellen, für jeden Zustand eine for-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-ähnlich zu sein. In SAS würde ich eine "by" -Anweisung machen und in SQL würde ich eine "group by" machen. Wie kann man das machen?

Answer 1

Hier ist eine Möglichkeit mit dem 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 

Answer 2

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

Meiner Meinung nach ist ein gemischtes lineares Modell ein besserer Ansatz für diese Art von Daten. Der unten angegebene Code wirkt sich auf den Gesamttrend aus. 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. Werfen Sie einen Blick auf Pinheiro & Bates (Mixed-Effekte-Modelle in S und S-Plus).

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

Answer 4

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) 

Answer 5

Die lm() Funktion oben ist ein einfaches Beispiel. By the way, ich stelle mir vor, dass die Datenbank die Spalten wie in der folgenden Form besitzt:

Jahre Zustand var1 var2 y ...

In 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) 

Answer 6

Eine nette Lösung mit data.table wurde in CrossValidated von @Zach here gepostet. Ich würde nur hinzufügen, dass es möglich ist, iterativ zu erhalten auch die Regressionskoeffizienten 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 andere Ausgabe 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 

Answer 7

Seit 2009 dplyr hat wurde veröffentlicht, die tatsächlich eine sehr schöne Möglichkeit bietet, diese Art von Gruppierung zu machen, ähnlich wie SAS.

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 

auf die Koeffizienten und Rsquared/p.value abzurufen, die eine broom Paket verwenden kann. Dieses Paket bietet:

drei S3-Generika: ordentlich, die statistische Ergebnisse eines Modells wie Koeffizienten einer Regression zusammenfasst; Augment, das Spalten zu den ursprünglichen Daten wie Vorhersagen, Residuen und Clusterzuordnungen hinzufügt; und Blick, die bietet eine einreihige Zusammenfassung der Modellebene 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) 

Answer 8

I kommt meine Antwort jetzt ein bisschen spät, aber ich war für eine ähnliche Funktionalität suchen.Es würde die eingebaute Funktion ‚durch‘ in R scheinen kann die Gruppierung leicht auch tun:

durch enthält das folgende Beispiel, die pro Gruppe passt 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) 

Answer 9

Die Frage scheint zu sein, wie man Regressionsfunktionen mit Formeln bezeichnet, die innerhalb einer Schleife modifiziert werden.

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

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