Lineare Regression und gruppieren nach in R

Ich möchte mit der function lm() eine lineare Regression in R durchführen. 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” -statement machen und in SQL würde ich eine “group by” machen. Wie kann man das machen?

    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 

    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) 

    Seit 2009 ist dplyr veröffentlicht worden, was tatsächlich eine sehr dplyr Möglichkeit bietet, diese Art von Gruppierung zu machen, ähnlich wie SAS es tut.

     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:  # # state model # (fctr) (chr) # 1 CA  # 2 NY  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 zu erhalten, kann man das Besenpaket verwenden. Dieses Paket bietet:

    drei S3-Generika: aufgeräumt, die die statistischen Ergebnisse eines Modells wie Koeffizienten einer Regression zusammenfassen; Augment, das den Originaldaten Spalten wie Vorhersagen, Residuen und Clusterzuordnungen hinzufügt; und glance, die eine einreihige 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) 

    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. Sehen Sie sich Pinheiro & Bates (Mixed-Effekt-Modelle in S und S-Plus) an.

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

    Eine nette Lösung mit data.table wurde hier in CrossValidated von @Zach gepostet. Ich würde 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 aus der 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 
     ## 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 yx 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 > 

    Ich komme jetzt etwas zu spät, aber ich habe nach einer ähnlichen function gesucht. Es scheint, dass die eingebaute function ‘by’ in R auch die Gruppierung leicht machen 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) 

    Die obige function lm() ist ein einfaches Beispiel. Übrigens stelle ich mir vor, dass Ihre database die Spalten wie folgt hat:

    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) 

    Ich denke, es lohnt sich, den purrr::map Ansatz zu diesem 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 = .))) 

    Siehe @Paul Hiemstras Antwort für weitere Ideen zur Verwendung des Besenpakets mit diesen Ergebnissen.

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

    Hier ist, wie Sie es tun können (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() }