Линейная регрессия и группа по R

Question

Я хочу сделать линейную регрессию в R, используя функцию lm(). Мои данные - это ежегодный временной ряд с одним полем в течение года (22 года), а другой для штата (50 штатов). Я хочу соответствовать регрессии для каждого состояния, так что в конце у меня есть вектор ответов lm. Я могу представить, как сделать цикл для каждого состояния, а затем выполнить регрессию внутри цикла и добавить результаты каждой регрессии к вектору. Однако это не похоже на R-like. В SAS я бы сделал оператор «by», а в SQL я бы сделал «группу». Каков способ R сделать это?

Answer 1

В одну сторону с помощью пакета 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 

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

На мой взгляд, смешанная линейная модель - лучший подход для данных такого типа. Приведенный ниже код содержит фиксированный эффект общей тенденции. Случайные эффекты показывают, как тренд для каждого отдельного государства отличается от глобальной тенденции. Корреляционная структура учитывает временную автокорреляцию. Посмотрите на Pinheiro & Bates (Модели смешанных эффектов в S и S-Plus).

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

Answer 4

Вот подход с использованием 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) 

Answer 5

Вышеприведенной lm() функцией является простым примером. Кстати, я думаю, что ваша база данных столбцов как в следующем виде:

год состояние var1 var2 у ...

С моей точки зрения, вы можете использовать следующий код:

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

Было найдено хорошее решение с использованием data.tablehere в CrossValidated by @Zach. Я бы только добавить, что можно получить итеративно также коэффициент регрессии г^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 

Answer 7

С 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 пакет. Этот пакет обеспечивает:

три S3 дженерик: опрятные, которые обобщаются статистических выводов МОДЕЛИ, такие как коэффициенты регрессии; дополнение, которое добавляет столбцы к исходным данным, таким как предсказания, остатки и назначения кластеров; и взгляд, который обеспечивает однострочное резюме статистики на уровне модели.

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

Я теперь мой ответ приходит немного поздно, но я искал аналогичной функциональностью.Казалось бы, встроенную функцию «на» в R также может сделать группировку легко:

по содержит следующий пример, который подходит для каждой группы и извлекает коэффициенты с 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

вопрос заключается в том, как вызвать функции регрессии с формулами, которые модифицируются внутри цикла.

Вот как вы можете сделать это в (с использованием алмазов набора данных):

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