R中的线性回归和分组

【问题标题】：Linear Regression and group by in RR中的线性回归和分组
【发布时间】：2010-11-13 06:17:37
【问题描述】：

我想使用 lm() 函数在 R 中进行线性回归。我的数据是一个年度时间序列，其中一个字段代表年（22 年），另一个字段代表州（50 个州）。我想为每个状态拟合一个回归，以便最后我有一个 lm 响应向量。我可以想象为每个状态执行 for 循环，然后在循环内执行回归并将每个回归的结果添加到向量中。然而，这似乎不太像 R。在 SAS 中，我会做一个“by”语句，而在 SQL 中，我会做一个“group by”。这样做的 R 方式是什么？

【问题讨论】：

标签： r regression linear-regression lm


【解决方案1】：

这是使用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)

【讨论】：

  • 假设您添加了一个额外的自变量，该变量并非在所有州都可用（即miles.of.ocean.shoreline），在您的数据中由 NA 表示。 lm 调用不会失败吗？怎么处理？
  • 在函数内部，您需要针对这种情况进行测试并使用不同的公式
  • 是否可以将子组的名称添加到摘要中的每个呼叫（最后一步）？
  • 如果你运行layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page 然后l_ply(models, plot) 你也会得到每个残差图。是否可以用组标记每个图（例如，在这种情况下为“状态”）？
【解决方案2】：

自 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 泛型：tidy，它总结了一个模型的 统计结果，例如回归系数； 增加，它将列添加到原始数据中，例如 预测、残差和聚类分配；和一瞥，其中 提供模型级统计信息的单行摘要。

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) %&gt;% tidy(model) 才能让扫帚包工作，但除此之外，很好的答案。
  • 效果很好...可以在不离开管道的情况下完成这一切：d %&gt;% group_by(state) %&gt;% do(model = lm(response ~ year, data = .)) %&gt;% rowwise() %&gt;% tidy(model)
  • @pedram 和 @holastello，这不再有效，至少在 R 3.6.1、broom_0.7.0、dplyr_0.8.3 中是这样。 d %&gt;% group_by(state) %&gt;% do(model=lm(response ~year, data = .)) %&gt;% rowwise() %&gt;% 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. ...
  • 现在（dplyr 1.0.4，tidyverse 1.3.0），你可以这样做：library(broom); library(tidyverse) d %&gt;% nest(data = -state) %&gt;% mutate(model = map(data, ~lm(response~year, data = .)), tidied = map(model, tidy)) %&gt;% unnest(tidied)
【解决方案3】：

这是使用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

【讨论】：

  • 有没有办法同时列出这两种型号的 R2？例如一年后添加一个 R2 列。还要为每个系数添加 p 值？
  • @ToToRo 在这里你可以找到一个可行的解决方案（迟到总比没有好）： Your.df[,summary(lm(Y~X))$r.squared,by=Your.factor] 其中: Y, X 和 Your.factor 是你的变量。请记住，Your.df 必须是 data.table 类
【解决方案4】：

在我看来，混合线性模型是处理此类数据的更好方法。下面的代码给出了固定效果中的整体趋势。随机效应表明每个州的趋势与全球趋势有何不同。相关结构考虑了时间自相关。看看 Pinheiro & Bates（S 和 S-Plus 中的混合效果模型）。

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

【讨论】：

  • 这是一个非常好的一般统计理论答案，它让我想到了一些我没有考虑过的事情。导致我提出问题的应用程序不适用于此解决方案，但我很高兴你提出来。谢谢。
  • 从混合模型开始不是一个好主意 - 你怎么知道任何假设都是有根据的？
  • 应该通过模型验证（和数据知识）来检查假设。顺便说一句，您也不能保证对个人 lm 的假设。您必须分别验证所有模型。
【解决方案5】：

@Zach 在 CrossValidated 中发布了一个使用 data.table 的好解决方案 here。 我只想补充一点，也可以迭代地获得回归系数 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

【讨论】：

    【解决方案6】：

    我认为为这个问题添加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 = .)))

    有关使用 broom 包获得这些结果的更多想法，请参阅 @Paul Hiemstra 的回答。

    【讨论】：

    • 如果您想要一列拟合值或残差，请稍作扩展：将 lm() 调用包装在 resid() 调用中，然后将最后一行中的所有内容传递到 unnest() 调用中。当然，您可能希望将变量名称从“模型”更改为更相关的名称。
    【解决方案7】：

    现在我的答案来的有点晚，但我一直在寻找类似的功能。看起来 R 中的内置函数 'by' 也可以轻松进行分组：

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

    【讨论】：

      【解决方案8】： 
      ## 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

      【讨论】：

        【解决方案9】：

        上面的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)

        【讨论】：

          【解决方案10】：

          问题似乎是关于如何使用在循环内修改的公式调用回归函数。

          以下是您可以使用的方法（使用 diamonds 数据集）：

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

          【讨论】：

            猜你喜欢
            • 1970-01-01
            • 1970-01-01
            • 2014-09-21
            • 2022-01-16
            • 2015-01-02
            • 2016-03-04
            • 1970-01-01
            • 1970-01-01
            相关资源
            最近更新 更多
            热门标签
            Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode