stats | 线性回归（五）——亚组模型、加权模型和逐步回归

本篇是线性回归系列的第五篇推文，也是最后一篇。示例数据如下：

data <- mtcars

8 亚组模型

在进行模型分析时，为了研究某分类变量对结果是否具有异质性影响，常会进行亚组分析，即根据某变量的取值将样本分成若干份再分别建模，然后比较各亚组的模型结果。

lm函数的subset参数可以对data参数取子集，这一点和基础绘图系统的某些绘图函数如barplot、boxplot等类似。

示例数据的am变量有两个取值：

unique(data$am)

## [1] 1 0

选择示例数据中am变量取值为0的样本：

model <- lm(mpg ~ wt + qsec, data = data,
            subset = am == 0)
coef(summary(model))

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 11.2489412  6.7148019  1.675245 0.1133158633
## wt          -2.9962762  0.6635548 -4.515491 0.0003520832
## qsec         0.9454396  0.2945500  3.209776 0.0054642663

同样，也可以对am取值为1的样本进行建模。但当亚组较多时，这样逐个进行建模十分不便，下面介绍一种向量式的操作方法，以cyl变量作为分组变量。

• 第一步：数据嵌套

library(dplyr)
library(tidyr)

data %>%
  group_by(cyl) %>%
  nest() -> DATA
DATA

## # A tibble: 3 x 2
## # Groups:   cyl [3]
##     cyl data              
##   <dbl> <list>            
## 1     6 <tibble [7 x 10]> 
## 2     4 <tibble [11 x 10]>
## 3     8 <tibble [14 x 10]>

• tidyr::nest函数配合dplyr::group_by函数使用可以实现数据的嵌套化。

第二步：调用向量化操作函数

关于向量化操作函数，本号已经有所介绍，详见下列推文：

• base | 使用apply族函数进行向量化运算

• purrr | 使用map族函数进行向量化运算

library(purrr)
library(magrittr)

# 模型运行
DATA %<>%
  mutate(model = map(data, ~lm(mpg ~ wt + qsec, 
                               data = .x)))
# 查看结果
map(DATA$model, ~coef(summary(.x)))

## [[1]]
##              Estimate Std. Error   t value    Pr(>|t|)
## (Intercept) 25.461734  4.9799376  5.112862 0.006920375
## wt          -5.201906  2.6364265 -1.973090 0.119745053
## qsec         0.583864  0.5504117  1.060777 0.348594034
## 
## [[2]]
##               Estimate Std. Error   t value   Pr(>|t|)
## (Intercept) 24.8842725 12.4286868  2.002164 0.08024671
## wt          -7.5135758  2.3291008 -3.225956 0.01213033
## qsec         0.9903892  0.7884782  1.256077 0.24452757
## 
## [[3]]
##               Estimate Std. Error   t value    Pr(>|t|)
## (Intercept) 14.0209319  7.7550402  1.807977 0.098001952
## wt          -2.8137538  0.8457285 -3.327018 0.006746824
## qsec         0.7352592  0.5369923  1.369217 0.198237772

9 加权模型

加权模型是在模型中对不同样本给予不同的权重，目的有如下几种：

• 样本量加权：如样本来自不同地区，我们希望样本在各地区的分布是均匀的，即样本量与地区人口成正比，但实际样本分布可能并不完全满足这一要求，这时就可以使用地区人口数对模型进行加权；

• 影响量加权：如在时间序列研究中，我们假定影响效应与时间远近有关，这时就可以使用时间变量对模型进行加权。

lm函数中通过weights参数指定加权变量：

model <- lm(mpg ~ wt + qsec, data = data,
            weights = cyl)
coef(summary(model))

##              Estimate Std. Error    t value     Pr(>|t|)
## (Intercept) 17.095422  5.0826107   3.363512 2.177700e-03
## wt          -4.674319  0.4652289 -10.047352 5.928451e-11
## qsec         1.002560  0.2620742   3.825481 6.412188e-04

10 逐步回归

前篇已经介绍了模型评价和比较的几种指标。实际上，为了快速选择出“最佳”的模型表达式，还可以使用逐步回归。stats工具包中用于逐步回归的函数为step：

step(object, scope, scale = 0,
     direction = c("both", "backward", "forward"),
     trace = 1, keep = NULL, steps = 1000, k = 2, ...)

主要参数有：

• object：包含所有待选自变量的模型；

• direction：逐步回归的方法；默认为both（所有子集法），可选项有backward（后退法）和forward（前进法）。

step函数是以AIC为评价标准进行模型选择的。

model <- lm(mpg ~ wt + I(wt^2) + drat + qsec,
            data = data)
step(model)

# Start:  AIC=55.33
## mpg ~ wt + I(wt^2) + drat + qsec
## 
##           Df Sum of Sq    RSS    AIC
## - drat     1     1.374 133.31 53.663
## <none>                 131.94 55.332
## - I(wt^2)  1    51.581 183.52 63.891
## - qsec     1    71.529 203.47 67.193
## - wt       1   122.327 254.27 74.325
## 
## Step:  AIC=53.66
## mpg ~ wt + I(wt^2) + qsec
## 
##           Df Sum of Sq    RSS    AIC
## <none>                 133.31 53.663
## - I(wt^2)  1    62.149 195.46 63.908
## - qsec     1    70.431 203.75 65.236
## - wt       1   169.437 302.75 77.910
## 
## Call:
## lm(formula = mpg ~ wt + I(wt^2) + qsec, data = data)
## 
## Coefficients:
## (Intercept)           wt      I(wt^2)         qsec  
##     32.6418     -12.4331       1.0730       0.8599

此例共经过两步就选择出了最优的模型表达式，即mpg ~ wt + I(wt^2) + qsec。逐步回归的结果具体解读如下：

• 第一步，分别比较原始模型和去掉一个自变量的模型的AIC大小，<none>对应的行表示原始模型的结果，- drat对应的行表示去掉drat变量的模型结果，按照AIC从小到大排列，结果显示去除drat变量的模型AIC最小，见结果的第一部分；

• 第二步，以去除drat变量的模型为原始模型，重复上述步骤，结果显示原始模型的AIC最小，逐步回归结束，见结果的第二部分；

• 输出最佳模型的结果，见结果的第三部分。

---

往期推荐阅读：

• 《数据处理通识》专辑-base | 使用apply族函数进行向量化运算

• 《制表与可视化》专辑-ggplot2 | ggplot2作图语法入门

• 《数学模型》专辑-car | 线性回归（三）——残差分析和异常点检验

• 《地理计算与分析》专辑-spdep | 如何在R语言中计算空间自相关指数