本篇是线性回归系列的第五篇推文,也是最后一篇。示例数据如下:
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
函数使用可以实现数据的嵌套化。
第二步:调用向量化操作函数
关于向量化操作函数,本号已经有所介绍,详见下列推文:
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最小,逐步回归结束,见结果的第二部分;输出最佳模型的结果,见结果的第三部分。
往期推荐阅读:
![](https://i-blog.csdnimg.cn/blog_migrate/30a8262ac0977e08406916b3f2abb2d1.jpeg)