本篇是线性回归系列的第五篇推文，也是最后一篇。示例数据如下：
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语言中计算空间自相关指数
