Die lineare Regression und die Gruppe, die durch R in
https://stackoverflow.com/questions/1169539
-
19-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich möchte eine lineare Regression in R tun, um die lm() Funktion. Meine Daten ist eine jährliche Zeitreihe mit einem Feld für Jahr (22 Jahre) und ein anderer für den Zustand (50 Staaten). Ich möchte eine Regression für jeden Zustand passen, so dass am Ende ich einen Vektor von lm Antworten haben. Ich kann für jeden Zustand dann die Regression in der Schleife und Addieren der Ergebnisse der einzelnen Regression auf einen Vektor zu tun für Schleife vorstellen. Das scheint nicht sehr R-like, aber. In SAS würde ich eine ‚durch‘ Erklärung und in SQL würde ich eine ‚Gruppe von‘ tun. Was ist die R Art und Weise, dies zu tun?
Lösung
Hier ist ein Weg, 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
Andere Tipps
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 dplyr wurde veröffentlicht, die tatsächlich eine sehr schöne Art und Weise, diese Art von Gruppierung zu tun, sehr ähnlich, was SAS hat.
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 abrufen und Rsquared / p.value, kann man das broom Paket verwenden. Dieses Paket enthält:
drei S3 Generika: ordentlich, die ein Modell der zusammenfasst statistische Ergebnisse wie Koeffizienten einer Regression; ergänzen, welche Spalten mit den Originaldaten wie fügt Vorhersagen, Residuen und Clusterzuweisungen; und Blick, die bietet eine einzeilige Zusammenfassung der Modell-Ebene 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)
Meiner Meinung nach ist ein gemischtes lineares Modell ein besserer Ansatz für diese Art von Daten. Der folgende Code in dem festen Effekt der allgemeine Trend gegeben. Die zufälligen Effekte zeigen, wie der Trend für jeden einzelnen Staat vom globalen Trend abweichen. Die Korrelationsstruktur nimmt die zeitliche Autokorrelation berücksichtigt. Hier finden Sie aktuelle Pinheiro & Bates (Mixed Effects Modelle in S und S-Plus).
library(nlme)
lme(response ~ year, random = ~year|state, correlation = corAR1(~year))
Eine schöne Lösung data.table mit gepostet hier in CrossValidated durch @Zach.
Ich würde nur hinzufügen, dass es möglich ist, iterativ auch die 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 Ausgang 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
## 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
ich jetzt meine Antwort kommt ein bisschen spät, aber ich war auf der Suche für eine ähnliche Funktionalität. Es würde die eingebaute Funktion ‚durch‘ in R scheint auch tun können, die Gruppierung leicht:
nach 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)
Ich denke, es lohnt sich, die purrr::map Ansatz 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 = .)))
Siehe @ Paul Hiemstra Antwort für weitere Ideen auf das broom Paket mit diesen Ergebnissen.
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:
Jahr 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)
Die Frage ist über scheint, wie Regressionsfunktionen mit Formeln nennen, die in einer Schleife verändert 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()
}
