【深度学习】图形模型基础(5)：线性回归模型第五部分：多变量线性回归模型

1.引言

当我们从基础的线性模型$y=a+bx+\text{error}$ 转向更复杂的模型$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\dots +\text{error}$ 时，我们面临了诸多挑战。这些挑战包括但不限于：

1. 预测变量的选择：确定哪些变量$x$应该被纳入模型中，这需要对数据进行深入分析和理解。
2. 系数的解释：理解每个系数${\beta }_i$ 的意义，以及它们如何共同影响模型的输出$y$。
3. 系数的相互作用：分析不同系数之间的相互作用，以及它们如何共同作用于模型的预测能力。
4. 新预测变量的构建：从现有的变量中创造新的预测变量，以捕捉数据中的离散性和非线性特征。

为了应对这些挑战，我们必须学会在引入新的预测变量时如何构建和理解模型。我们将通过一系列实例来探讨这些问题，这些实例将通过R代码、数据图表以及拟合模型的图形来展示。这些实例不仅有助于我们理解模型构建的过程，也有助于我们更好地把握模型的预测能力和局限性。通过这种方式，我们可以更深入地理解数据，更准确地预测结果，并更有效地应对模型中的复杂性。

2.在模型中增加变量

在模型中添加预测变量是一个复杂的过程，因为当多个预测变量存在时，回归系数的解释变得更加复杂。每个系数的解释部分依赖于模型中的其他变量。例如，系数 ${\beta }_k$ 表示在所有其他预测变量保持不变的情况下，两个个体在预测变量 $x_k$ 上相差一个单位时，结果变量 $y$ 的平均或预期差异。我们通过一系列例子来说明这一点，这些例子包括单一预测变量和多个预测变量的组合。

2.1.双变量预测模型

首先，我们考虑一个二元预测变量，即母亲是否高中毕业的指标，用1表示高中毕业，0表示没有。我们使用R语言来拟合这个模型，代码如下：

model_1 <- lm(kid_score ~ mom_hs, data = kidiq)
summary(model_1)

这个模型的公式可以表示为：

$$\text{kid\_score}=78+12\times \text{mom\_hs}+\text{error}\text{\hspace{1em}(5.1)}$$

其中，截距78表示母亲没有完成高中教育的儿童的平均测试分数。如果母亲完成了高中教育，儿童的平均测试分数预计会增加12分。

图5.1 展示了儿童的测试分数与母亲是否完成高中教育的指标之间的关系。图中叠加了回归线，这条线穿过由母亲教育水平定义的每个子群体的平均值。表示高中完成情况的指示变量已经进行了抖动处理；也就是说，每个$x$值上都加上了一个随机数，以防止这些点相互重叠。

2.2.连续单变量模型

接下来，我们尝试使用一个连续预测变量，即母亲的智商测试分数。拟合的模型如下：

model_2 <- lm(kid_score ~ mom_iq, data = kidiq)
summary(model_2)

模型公式为：

$$\text{kid\_score}=26+0.6\times \text{mom\_iq}+\text{error}\text{\hspace{1em}(5.2)}$$

这意味着，如果两个儿童的母亲智商相差1分，儿童的测试分数平均相差0.6分。如果智商相差10分，儿童的测试分数平均相差6分。

2.3.双变量模型的拟合

最后，我们考虑同时包括两个预测变量的线性回归模型。在R中，我们这样拟合模型：

model_3 <- lm(kid_score ~ mom_hs + mom_iq, data = kidiq)
summary(model_3)

Median MAD_SD
(Intercept) 25.7 5.9
mom_hs 6.0 2.4
mom_iq 0.6 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 18.2 0.6

2.4.双变量线性回归模型的理解

这个模型的公式为：

$$\text{kid\_score}=25.7+6\times \text{mom\_hs}+0.6\times \text{mom\_iq}+\text{error}\text{\hspace{1em}(5.3)}$$

图5.2 展示了儿童的测试分数与母亲的智商之间的关系，并叠加了回归线。线上的每个点可以被理解为，对于具有相应智商的母亲的孩子，预测的儿童测试分数，或者作为具有该智商的母亲的子群体儿童的平均分数。

在这个模型中，儿童的测试分数对母亲智商的回归斜率对于每个母亲的教育水平子群体是相同的。截距25.7在实际中没有太大意义，因为它假设母亲的智商为0，这在现实中是不存在的。在模型(5.3)中，$\text{mom\_hs}$的系数6表示，在智商相同的情况下，母亲是否完成高中教育对儿童测试分数的预期影响是6分。而$\text{mom\_iq}$的系数0.6表示，在母亲是否完成高中教育相同的情况下，母亲的智商每增加1分，儿童的测试分数预计增加0.6分。

图5.3 展示了儿童的测试分数与母亲的智商之间的关系。浅点代表那些母亲完成了高中教育的儿童，而深点代表那些母亲没有完成高中教育的儿童。图中叠加了基于儿童测试分数对母亲智商和母亲是否完成高中教育指标进行回归的线（深色线代表母亲没有完成高中教育的儿童，浅色线代表母亲完成了高中教育的儿童）。

3.多元线性回归系数的含义

在多元线性回归分析中，解释回归系数的过程并不总是直观的，因为我们通常假设在比较不同个体时，可以改变一个预测变量而保持其他变量不变。然而，在现实中，某些情况下这种假设并不成立。例如，如果模型中同时包含智商（IQ）和智商平方（IQ^{2）作为预测变量，那么在保持IQ}2不变的情况下单独改变IQ是没有意义的。同样，如果模型中包括母亲是否完成高中教育（mom_hs）、母亲的智商（mom_IQ）以及它们的交互项（mom_hs:mom_IQ），单独考虑任何一个变量而保持其他两个变量不变也是没有意义的。

我们对回归系数的解释通常有两种视角：预测性解释和反事实解释。

1. 预测性解释：这种解释关注的是，当比较两组在某一预测变量上相差1个单位但在所有其他预测变量上相同的个体时，结果变量平均上的差异。在线性模型中，这个系数代表了这两个项目之间预期的y值差异。这是我们之前描述过的解释方式。

2. 反事实解释：这种解释关注的是个体内部的变化，而不是个体之间的比较。在这里，系数表示的是在模型中保持所有其他预测变量不变，向相关预测变量增加1个单位所预期的y值变化。例如，如果我们说“将母亲的智商从100提高到101，预计会导致孩子的测试分数增加0.6分”，这就是一种反事实的解释。

尽管有时统计和回归的入门教材会警告不要使用反事实解释，但它们仍然会使用类似的表述，例如“母亲的智商每增加10分，孩子的测试分数就与增加6分有关”。然而，仅凭数据本身，回归分析只能告诉我们不同单位之间的比较，并不能告诉我们单位内部的变化。

因此，对回归系数最准确的解释应该是基于比较的，例如：“当比较两个母亲教育水平相同的孩子时，如果一个母亲的智商比另一个高出x分，那么预计这个孩子的测试分数平均会高出6x分。”或者，“如果两个项目i和j在预测变量k上相差x个单位，但在所有其他预测变量上相同，那么预测的y值差异$y_i-y_j$平均是 ${\beta }_k\times x$。”

这种表述可能显得笨拙，但它有助于解释为什么人们通常更喜欢使用简单的表述，如“$x_k$的变化1个单位导致或与y的变化β个单位有关”，尽管这些表述可能极具误导性。我们必须接受，回归分析虽然是一个强大的数据分析工具，但其解释有时可能相当复杂。我们将在后续的文章中进一步探讨在何种条件下可以对回归进行因果解释。

4.变量的交互作用

在模型(5.3)中，我们假定儿童测试分数对母亲智商的回归斜率在母亲是否完成高中教育的各个子群体中是相同的。然而，图5.3的数据观察表明，这些斜率实际上有显著的差异。为了解决这个问题，我们可以引入mom_hs和mom_iq之间的交互项，即这两个变量乘积的新预测变量，从而允许斜率在不同子群体中变化。在R中，我们可以通过以下代码来拟合并展示包含交互项的模型：

fit_4 <- lm(kid_score ~ mom_hs * mom_iq, data = kidiq)
summary(fit_4)

拟合后的模型可以表示为：

$$\text{kid\_score}=-11+51\times \text{mom\_hs}+1.1\times \text{mom\_iq}-0.5\times \text{mom\_hs}\times \text{mom\_iq}+\text{error}$$

这个模型的图形表示分别在图10.4a和图10.4b中展示。图10.4a为不同母亲教育水平子群体的回归线，而图10.4b则将x轴扩展到零，以展示截距，即线在y轴上穿过零点的位置。这突出了一个重要问题：当预测变量的可能值范围远离零时，截距没有直接的解释意义。

图5.4 (a) 展示了儿童测试分数对母亲智商的回归线，其中用不同的符号表示母亲完成了高中教育的儿童（浅圈）和母亲未完成高中教育的儿童（深点）。交互作用允许每个组别有不同的斜率，浅线和深线分别对应浅点和深点。 (b) 是相同的图表，但是水平和垂直轴都扩展到零，以显示截距。

在解释这个模型的系数时，我们需要谨慎。我们通过检视特定子群体内和跨特定子群体的平均或预测测试分数来从拟合模型中获得含义。一些系数仅对特定子群体可解释：

• 截距项代表母亲没有完成高中教育且智商为0的儿童的预测测试分数，但这是一个不具现实意义的假设场景。如果我们在将变量纳入回归预测之前先对输入变量进行中心化处理，截距将更具解释性。
• mom_hs的系数可以被理解为，在智商为0的情况下，未完成高中教育的母亲的儿童与完成高中教育的母亲的儿童预测测试分数之间的差异。这个系数的解释需要通过代入适当的数值并比较方程来实现。由于智商为0的母亲是不可想象的，这个系数的解释并不直观。
• mom_iq的系数可以被看作是未完成高中教育的母亲的儿童在智商上每增加1点时，平均测试分数的变化。这对应于图5.4中深色线的斜率。
• 交互项的系数表示了mom_iq斜率的差异，即完成了与未完成高中教育的母亲的儿童之间的比较，也就是图5.4中浅色线和深色线斜率之间的差异。

我们可以通过为完成和未完成高中教育的母亲的儿童分别查看回归线来理解模型：

对于未完成高中教育的母亲（$\text{mom\_hs}=0\}$）：
$$\text{kid\_score}=-11+1.1\times \text{mom\_iq}$$

对于完成高中教育的母亲（$\text{mom\_hs}=1\}$）：
$$\text{kid\_score}=40+0.6\times \text{mom\_iq}$$

这里，对于未完成高中教育的母亲的儿童，估计的斜率是1.1；对于完成高中教育的母亲的儿童，斜率是0.6。这些斜率是直接可解释的。然而，截距仍然面临一个问题，即仅在母亲的智商为0时才具有可解释性。

我们应该何时寻找交互作用？交互作用可能非常重要，我们通常首先在那些在没有交互作用时具有大系数的预测变量中寻找它们。例如，吸烟与癌症有很强的关联。在其他致癌物的流行病学研究中，调整吸烟作为一个未交互的预测变量和作为一个交互项是至关重要的，因为其他风险因素与癌症之间的关联强度可能取决于个体是否吸烟。我们通过家庭氡暴露的例子在图1.7中说明了这种交互作用：高水平的氡与更大的癌症可能性相关，但这种差异对于吸烟者来说比非吸烟者要大得多。包括交互作用是一种方式，允许模型以不同的方式适应数据的不同子集。

在存在交互作用的情况下解释回归系数时，如果通过将每个输入变量中心化到其均值或其他方便的参考点来预处理数据，通常可以更容易地解释模型。我们将在线性变换的背景下讨论这一点。

5.交互变量

在先前的章节中，我们探讨了利用指示变量在回归分析中进行比较的方法。现在，我们将此概念进一步应用于基于调查数据的模型拟合中，以预测身高并考虑体重和其他变量的影响。以下是来自数据集earnings的样本数据，其中包括身高、体重、性别、收入、种族、教育水平、步行量、运动量、是否吸烟、紧张程度、愤怒程度和年龄等变量。

为了从身高预测体重，我们首先建立了一个简单的线性回归模型：

fit_1 <- lm(weight ~ height, data = earnings)
summary(fit_1)

模型估计结果如下：

• 截距（Intercept）：(-172.9)
• 身高（height）的系数：(4.9)

辅助参数，即预测误差的标准差（sigma）为：(29.1)

得到的回归方程可以表示为：
[ \text{weight} = (-172.9) + 4.9 \times \text{height} ]

这个方程告诉我们，身高每增加一英寸，体重平均增加4.9磅。然而，对于身高为0英寸的人，模型预测的体重是-172.9磅，这在实际情况下是没有意义的。鉴于美国成年人的平均身高约为66英寸，我们可以重新表述模型对于平均身高人群的预测：

coefs_1 <- coef(fit_1)
predicted_1 <- coefs_1["(Intercept)"] + coefs_1["height"] * 66

或者，使用posterior_predict函数来预测特定身高的体重，并得到预测值的分布情况：

new_data <- data.frame(height = 66)
pred <- posterior_predict(fit_1, newdata = new_data)
cat("Predicted weight for a 66-inch-tall person is", round(mean(pred)),
    "pounds with a sd of", round(sd(pred)), "\n")

无论是通过直接使用模型系数进行计算，还是通过posterior_predict函数获取预测的分布，对于66英寸高的个体，我们都可以预测其体重，并且得到一个更加符合实际的预测结果。

5.1.中心化预测变量

为了增强模型预测结果的可解释性，我们对身高这一预测变量执行了中心化处理。中心化是一种常用的数据预处理方法，通过从原始数据中减去一个基准值（通常是均值）来实现。在我们的数据集earnings中，我们首先计算出所有个体的平均身高，然后从每个个体的身高中减去这个均值，得到中心化身高c_height：

earnings$c_height <- earnings$height - 66  # 假设66是身高的均值

接下来，我们使用这个中心化后的身高变量来拟合一个新的回归模型fit_2，以预测体重weight：

fit_2 <- stan_glm(weight ~ c_height, data = earnings)

这个模型的估计结果提供了以下信息：

• 截距（Intercept）：153.4，这表示在中心化身高为0时，模型预测的体重。
• 中心化身高（c_height）的系数：4.9，这表示中心化身高每增加一个单位，模型预测的体重将平均增加4.9磅。

此外，模型还提供了一个辅助参数，即预测误差的标准差sigma：

• sigma：29.1，这表示预测体重值围绕平均预测值的离散程度。

通过这种方式，中心化身高不仅使我们能够更直观地理解身高对体重影响的大小，还有助于我们更准确地进行模型预测。例如，如果我们知道某个人的身高是72英寸，而平均身高是66英寸，那么他们的中心化身高就是6英寸。根据模型，我们可以预测这个人的体重会比平均体重高出4.9 * 6磅。

5.2.二元回归分析

在构建更全面的回归模型时，我们可能会考虑将二元变量纳入分析，例如性别。在下面的文本中，我们将演示如何在回归模型中包含性别这一指示变量，并解释其对预测结果的影响。

我们首先通过引入性别变量male来扩展我们的模型，其中男性用1表示，女性用0表示。以下是在R语言中实现的代码：

fit_3 <- stan_glm(weight ~ c_height + male, data = earnings)

模型估计的结果如下：

• 截距：149.6
• c_height（以英寸为单位的中心化身高）的系数：3.9
• male（性别）的系数：12.0

辅助参数包括估计的标准差sigma，其值为28.8。

性别系数12.0的值表明，在控制了身高变量后，男性相较于女性，其体重平均高出12磅。这一发现提供了有关性别对体重影响的重要信息。

为了具体计算一个身高为70英寸的女性的预测体重，我们可以使用模型系数直接进行计算：

coefs_3 <- coef(fit_3)
predicted_female <- coefs_3["(Intercept)"] + coefs_3["c_height"] * 4.0 + coefs_3["male"] * 0

或者，使用posterior_predict函数来获取基于模型的预测分布，并计算平均值：

new_female <- data.frame(c_height = 4.0, male = 0)
pred_female <- posterior_predict(fit_3, newdata = new_female)
cat("Predicted weight for a 70-inch-tall woman is", round(mean(pred_female)),
    "pounds with a sd of", round(sd(pred_female)), "\n")

无论使用哪种方法，我们得到的女性预测体重都是165磅。对于相同身高的男性，相应的预测体重计算方法类似，但将性别变量male的值设为1：

new_male <- data.frame(c_height = 4.0, male = 1)
pred_male <- posterior_predict(fit_3, newdata = new_male)
cat("Predicted weight for a 70-inch-tall man is", round(mean(pred_male)),
    "pounds with a sd of", round(sd(pred_male)), "\n")
<