R 中的多元线性回归：示例教程

回归方法被广泛用于各行各业，以了解哪些变量会影响某一关注主题。

例如，经济学家可以用它来分析居民消费与国内生产总值（GDP）增长之间的关系。公共卫生官员可能希望基于个体的历史信息来理解其成本。在这两种情况下，关注点并非预测个体情境，而是把握总体关系的概貌。

在本文中，我们将先给出对回归的一般性理解。然后解释简单线性回归与多元线性回归有何不同，接着进入技术实现，并提供帮助您理解与解读回归结果的工具。

什么是多元线性回归？

在讨论多元线性回归之前，我们先了解什么是简单线性回归，因为多元线性回归只是简单线性回归的扩展。

简单线性回归 vs. 多元线性回归

简单线性回归旨在通过估计当单个自变量 X 变化一定量时，因变量 Y 将发生多大变化，来建模两者之间的关系。

• 自变量 X，也称预测变量，是用于进行预测的变量。

• 因变量 Y，也称响应变量，是我们试图预测的对象。

线性回归中的“线性”体现在我们尝试使用以下“线性”方程从 X 预测 Y。

Y = b0 + b1X 

• b0 是回归线的截距，对应于当 X 为零时的预测值。

• b1 是回归线的斜率。

那么，多元线性回归呢？

这就是在多个变量上使用线性回归，其方程为：

Y = b0 + b1X1 + b2X2 + b3X3 + … + bnXn + e

• Y 与 b0 与简单线性回归模型中的相同。

• b1X1 表示第一个自变量（X1）的回归系数（b1）。其余回归系数与变量亦作同样解读。

• e 是模型误差（残差），表示在估计 Y 时模型引入的变动量。

在多元回归的情形下，我们未必总能得到一条直线。不过，我们可以通过拟合更合适的模型来控制曲线的形状。

以下是多元线性回归为每个预测变量计算以寻找最佳拟合线的一些关键要素。

• 每个预测变量的估计系数。
• 模型的 R 平方，表示模型解释的方差比例，并以 0 到 100% 的便捷刻度来衡量模型与因变量 Y 之间关系的强度。
• 模型中每个回归系数的p 值与 t 统计量。这两个指标相互呼应，用于衡量某个系数在统计上的显著性。t 统计量越高（伴随接近 0 的 p 值），该预测变量越显著，意味着应将其保留在模型中。相反，非常低的 t 统计量（较高的 p 值）意味着应当剔除该预测变量。

多元线性回归的假设

构建多元线性回归模型时，一个重要方面是确保满足以下关键假设。

• 残差值服从正态分布。这可以通过正态概率图或直方图来检查。
• 因变量与自变量之间必须存在线性关系。可以通过散点图展示线性或曲线型关系来说明。
• 另一项假设是多重共线性，即自变量之间不应高度相关。多重共线性会使我们难以辨明哪些变量更好地解释因变量。该假设可通过计算所有自变量之间的皮尔逊双变量相关矩阵来验证。如果数据中不存在共线性，则所有值应小于 0.8。
• 同方差性假设认为，残差误差的方差在各个自变量取值范围内应大致相同。可通过将预测值与标准化残差作图来检查，看点是否在所有自变量取值范围内均匀分布。

在接下来的部分中，我们将覆盖其中的一些假设。

R 中多元线性回归的分步指南

本节我们将使用 R 编程语言深入实现一个多元线性回归模型。

我们将使用DataCamp Workspace 的客户流失数据集来估计客户价值。

我们所说的客户价值指什么？基本上，它衡量某个产品或服务对客户而言的价值，我们可以按如下方式计算：

Customer Value = Benefit — Cost. 其中 Benefit 和 Cost 分别表示某产品或服务的收益与成本。

如果公司能够为消费者提供更高的收益与更低的成本，或理想情况下两者兼具，则该数值更高。

此分析可帮助企业基于特定客户的价值，识别最有潜力的目标机会或下一步最佳行动。

快速数据探索

先对数据集做个快速概览，以便在拟合模型前应用合适的预处理。

churn_data = read_csv('data/customer_churn.csv', show_col_types = FALSE)
# Look at the first 6 observations
head(churn_data)
# Check the dimension
dim(churn_data)技术支持

数据前 6 行（作者制作的动画）

从上述结果可见，该数据集包含 3150 条观测与 14 列。

不过，基于问题陈述，我们不需要 churn 列，因为我们现在处理的是回归问题。

在拟合模型之前，为避免每次都要在变量名外加双引号，我们将列名中的空格替换为下划线，对列名做个预处理。

# Change the column names
names(churn_data) = gsub(" ", "_", names(churn_data))

head(churn_data)技术支持

列名转换后的前 6 行（作者制作的动画）

拟合多元线性回归模型

使用这份新格式化的数据，我们可以如下通过 R 中的 lm() 函数将其拟合到多元回归框架：

# Fit the multiple linear regression model
cust_value_model = lm(formula = Customer_Value ~ Call_Failure + 

                     Complaints + Subscription_Length + Charge_Amount +    

                     Seconds_of_Use +Frequency_of_use + Frequency_of_SMS +
                    Distinct_Called_Numbers + Age_Group + Tariff_Plan +
                    Status + Age,data = churn_data)技术支持

我们来理解一下刚才做了什么。

lm() 函数的格式如下：lm(formula = Y ~Sum(Xi), data = our_data)

• Y 为 Customer_Value 列，因为这是我们要估计的对象。

• Sum(X_i) 表示多元线性回归方程中的求和表达式。

• our_data 即 churn_data。

您可以通过我们的课程 Intermediate Regression in R 了解更多。

R 的替代选择是 Python：Intermediate Regression with statsmodels in Python。两者都能帮助您学习含多个解释变量的线性与逻辑回归。

多元线性回归模型告诉了我们什么？

既然我们已经建立了模型，下一步是检验假设并解读结果。为简明起见，我们不会覆盖所有方面。

模型残差的分布

可在 R 中使用 hist() 函数展示。

# Get the model residuals
model_residuals = cust_value_model$residuals

# Plot the result
hist(model_residuals)技术支持

模型残差分布（作者制图）

直方图看起来左偏，因此我们无法有足够信心断言其正态性。与其看直方图，不如观察残差的正态 Q-Q 图。如果满足正态性，点应沿一条直线分布。

# Plot the residuals
qqnorm(model_residuals)
# Plot the Q-Q line
qqline(model_residuals)技术支持

Q-Q 图与残差（作者制图）

从图中可以看到，只有部分残差位于直线上。因此我们可以认为该模型的残差并不服从正态分布。

多重共线性假设检验

可通过以下 R 代码完成。但我们需要先移除 Customer_Value 列。

# Install and load the ggcorrplot package
install.packages("ggcorrplot")
library(ggcorrplot)

# Remove the Customer Value column
reduced_data <- subset(churn_data, select = -Customer_Value)

# Compute correlation at 2 decimal places
corr_matrix = round(cor(reduced_data), 2)

# Compute and show the  result
ggcorrplot(corr_matrix, hc.order = TRUE, type = "lower",
          lab = TRUE)技术支持

数据的相关性结果（作者制图）

我们可以注意到两组较强的相关性，因为它们的值高于 0.8。

• Age 与 Age_Group：0.96

• Frequency_of_use 与 Seconds_of_Use：0.95

这一结果合乎情理，因为 Age_Group 由 Age 计算而来。同样地，总通话秒数（Seconds_of_Use）由通话总次数（Frequency_of_Use）推导而来。

在这种情况下，我们可以从数据集中去除 Age_Group 与 Seconds_of_Use。

我们尝试在没有这两个变量的情况下构建第二个模型。

second_model = lm(formula = Customer_Value ~ Call_Failure + Complaints +
                    Subscription_Length + Charge_Amount +
                    Frequency_of_use + Frequency_of_SMS +
                    Distinct_Called_Numbers + Tariff_Plan +
                    Status + Age,
                    data = churn_data)技术支持

第一模型（左）与第二模型（右）的 Q-Q 图

可以看到，消除数据中的多重共线性是有帮助的，因为在第二个模型中，相较于第一个模型，更多残差点落在直线上。

如何比较多元线性回归模型

现在我们面临一个问题：这两个多元线性回归模型哪个更好？

一种方法是对两个模型进行方差分析（ANOVA）检验。它检验原假设（H0）——我们先前移除的变量并不显著——相对于备择假设（H1）——这些变量是显著的。

如果新模型优于原模型，则我们无法拒绝 H0。反之，则意味着那些变量是显著的，因此我们拒绝 H0。

通用表达式为：anova(original_model, new_model)

# Anova test
anova(cust_value_model, second_model)技术支持

ANOVA 检验结果（作者制图）

从 ANOVA 结果可见，p 值（8.0893e-316）非常小（小于 0.05），因此我们拒绝原假设，意味着第二个模型并未优于第一个模型。

另一种查看模型中重要变量的方法是显著性检验。

若某变量的 p 值小于 0.05，则其具有显著性。该结果可以通过 summary() 函数生成。除了提供模型的相关信息，它还会给出调整后的 R 平方，用于在模型之间比较性能。

# Print the result of the model
summary(cust_value_model)技术支持

包含所有预测变量的原始模型摘要结果（作者制图）

表格中有两个关键部分：Residuals 与 Coefficients。Q-Q 图与 Residuals 部分传达了相同的信息。在 Coefficients 部分，Call_Failure、Complaints 和 Age_Group 的 p 值大于 0.05，被模型视为不显著。保留它们不会为模型带来额外价值。

将同样的分析应用于第二个模型，我们得到如下结果：

summary(second_model)技术支持

包含所有预测变量的第二个模型摘要结果（作者制图）

原始模型的调整 R 平方为 0.98，高于第二个模型的调整 R 平方（0.97）。这意味着包含所有预测变量的原始模型优于第二个模型。

改进多元线性回归模型

此分析的合乎逻辑的下一步是移除不显著变量并重新拟合模型，看看性能是否提升。

另一种高效选择相关预测变量的策略是使用赤池信息准则（AIC）。

它从包含全部特征开始，然后每次逐步剔除最差的预测变量，直到找到最佳模型。AIC 分数越小，模型越好。这可以通过 stepAIC() 函数实现。

结论

本教程涵盖了多元线性回归的主要方面，并探讨了构建稳健模型的一些策略。

我们希望本教程能帮助您掌握从数据中获取可操作洞见所需的技能。您可以通过不同的方法来改进这些模型，欢迎在我们的 Workspace 提供的源代码上尝试。