Я хочу выполнить линейную регрессию в R с помощью lm()функции. Мои данные представляют собой годовой временной ряд с одним полем для года (22 года) и другим для штата (50 штатов). Я хочу подобрать регрессию для каждого состояния, чтобы в конце у меня был вектор ответов lm. Я могу представить, что выполняю цикл for для каждого состояния, а затем выполняю регрессию внутри цикла и добавляю результаты каждой регрессии в вектор. Однако это не кажется очень похожим на R. В SAS я бы сделал оператор «by», а в SQL - «group by». Как это сделать с помощью R?
r
regression
linear-regression
lm
JD Long
источник
источник
aggregateнеправильный ; ни то, ни другоеtapply.Ответы:
Вот один из способов использования
lme4пакета.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источник
Вот подход с использованием пакета plyr :
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)источник
layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/pageа затемl_ply(models, plot)вы получите также каждый из графиков остатков. Можно ли пометить каждый график группой (например, в данном случае «состояние»)?С 2009 года
dplyrбыл выпущен, который на самом деле обеспечивает очень хороший способ группировки, очень похожий на то, что делает 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Чтобы получить коэффициенты и Rsquared / p.value, можно использовать
broomпакет. Этот пакет обеспечивает: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)источник
rowwise(fitted_models) %>% tidy(model)чтобы пакет с метлой работал, но в остальном отличный ответ.d %>% group_by(state) %>% do(model = lm(response ~ year, data = .)) %>% rowwise() %>% tidy(model)d %>% group_by(state) %>% do(model=lm(response ~year, data = .)) %>% rowwise() %>% tidy(model) Error in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm) : Calling var(x) on a factor x is defunct. Use something like 'all(duplicated(x)[-1L])' to test for a constant vector. In addition: Warning messages: 1: Data frame tidiers are deprecated and will be removed in an upcoming release of broom. ...На мой взгляд, смешанная линейная модель лучше подходит для такого рода данных. В приведенном ниже коде фиксировано влияние общей тенденции. Случайные эффекты показывают, насколько тенденция для каждого отдельного состояния отличается от глобальной тенденции. Корреляционная структура учитывает временную автокорреляцию. Взгляните на Pinheiro & Bates (модели со смешанными эффектами в S и S-Plus).
library(nlme) lme(response ~ year, random = ~year|state, correlation = corAR1(~year))источник
Хорошее решение с использованием
data.tableбыло размещено здесь в CrossValidated @Zach. Я просто добавлю, что можно итеративно получить коэффициент регрессии 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а также все остальные выходные данные
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источник
Думаю, стоит добавить
purrr::mapподход к этой проблеме.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 = .)))См. Ответ @Paul Hiemstra для получения дополнительных идей по использованию
broomпакета с этими результатами.источник
## 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источник
Теперь я немного запоздал со своим ответом, но я искал аналогичную функциональность. Казалось бы, встроенная функция by в R также может легко группировать:
? by содержит следующий пример, который подходит для каждой группы и извлекает коэффициенты с помощью 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)источник
lm()Выше функция представляет собой простой пример. Кстати, я предполагаю, что ваша база данных имеет столбцы в следующей форме:год состояние var1 var2 y ...
На мой взгляд, вы можете использовать следующий код:
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)источник
Вопрос, похоже, заключается в том, как вызывать функции регрессии с формулами, которые изменяются внутри цикла.
Вот как вы можете это сделать (используя набор данных алмазов):
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() }источник