> (线性)回归:社会科学实证研究的主力方法 下文讨论的所有示例文件均可从 “数据> 管理” 页面加载。点击`examples`单选按钮并按 “加载(Load)”。 ### 功能说明 首先选择响应变量和一个或多个解释变量。如果模型中包含两个或多个解释变量,我们可能需要研究是否存在交互效应。当一个解释变量对响应变量的影响至少部分由另一个解释变量的水平决定时,交互效应就存在。例如,1 克拉与 2 克拉钻石的价格涨幅可能取决于钻石的净度等级。 在 “摘要(Summary)” 标签页中,我们可以通过在 “待检验变量(Variables to test)” 下拉菜单中选择变量,检验两个或多个变量是否共同对模型拟合有显著贡献。此功能对于检验因子(factor)类型变量的整体影响是否显著非常有用。 需要重新估计的额外输出: * 标准化(Standardize):如果解释变量的测量尺度不同,系数可能难以比较。通过在估计前对响应变量和解释变量进行标准化,我们可以看出哪些变量的影响更大。Radiant 中对数据的标准化方法是将响应变量Y替换为(Y−mean(Y))/(2×sd(Y)),并将所有解释变量X替换为(X−mean(X))/(2×sd(X))。详见Gelman 2008的讨论。 * 中心化(Center):将响应变量Y替换为Y−mean(Y),并将所有解释变量X替换为X−mean(X)。这在解释交互效应时可能有用。 * 逐步回归(Stepwise):一种数据挖掘方法,用于选择拟合效果最佳的模型。使用时需谨慎! * 稳健标准误(Robust standard errors):选择 “稳健(robust)” 后,系数估计值与普通最小二乘(OLS)相同,但标准误会进行调整,以解决(轻微的)异方差性和非正态性问题。 无需重新估计的额外输出: * RMSE:均方根误差(Root Mean Squared Error)和残差标准差(Residual Standard Deviation) * 平方和(Sum of Squares):响应变量的总方差分解为回归解释的方差和未解释的方差(即误差) * VIF:方差膨胀因子(Variance Inflation Factors)和 R 平方(Rsq),这些是解释变量间多重共线性的度量指标 * 置信区间(Confidence intervals):系数的置信区间 “预测(Predict)” 标签页允许计算回归模型的预测值。你可以选择基于数据集预测响应(即从 “预测输入(Prediction input)” 下拉菜单中选择`Data`)、基于命令预测(即从 “预测输入” 下拉菜单中选择`Command`),或结合两者(即从 “预测输入” 下拉菜单中选择`Data & Command`)。 如果选择`Command`,必须至少指定一个变量和值才能获得预测结果。如果未为模型中的每个变量指定值,则会使用均值或最频繁的水平。只能基于模型中的变量预测结果(例如,要预测 2 克拉钻石的`price`,`carat`必须是所选解释变量之一)。 * 要预测 1 克拉钻石的价格,输入`carat = 1`并按回车 * 要预测从 0.5 到 1 克拉、步长为 0.05 的钻石价格,输入`carat = seq(.5, 1, .05)`并按回车 * 要预测 1、2 或 3 克拉理想切工钻石的价格,输入`carat = 1:3, cut = "Ideal"`并按回车 生成所需预测后,可通过点击屏幕右上角的下载图标将其保存为 CSV 文件。要将预测结果添加到用于估计的数据集,点击 “存储(Store)” 按钮。 “绘图(Plot)” 标签页用于提供数据的基本可视化以及验证回归模型的诊断图。 ### 示例 1:目录销售数据 我们获取了一家通过邮购目录销售男女服装的公司的数据(数据集`catalog`)。该公司维护着关于过往和当前客户价值及特征的数据库。客户价值定义为过去一年对客户的总销售额(美元)。数据是从公司数据库中随机抽取的 200 名客户,包含以下 4 个变量: - `Sales` = 家庭过去一年的总销售额(美元) - `Income` = 家庭收入(千美元) - `HH.size` = 家庭规模(人数) - `Age` = 家庭户主年龄 该目录公司有意重新设计其客户关系管理(CRM)策略。我们将分步骤进行: 1. 使用去年的销售总额估计回归模型。响应变量:200 个家庭的销售总额;解释变量:家庭收入(以千美元计)、家庭规模和家庭户主年龄。要获取该数据集,前往 “数据> 管理”,从 “加载数据类型(Load data of type)” 下拉菜单中选择`examples`,按 “加载” 按钮,然后选择`catalog`数据集。 2. 解释每个估计系数。同时对模型整体进行统计评估。 3. 哪些解释变量是客户价值的显著预测因子(使用 95% 置信水平)? **答案:** 选择上述相关变量,按 “估计模型(Estimate model)” 按钮或按`CTRL-enter`(Mac 上为`CMD-enter`)。“模型> 线性回归(OLS)” 的输出如下:

F 检验的原假设和备择假设可表述为: - H0:所有回归系数均等于 0 - Ha:至少有一个回归系数不等于 0 F 统计量表明,回归模型整体解释了`Sales`的显著方差。计算得到的 F 统计量为 32.33,p 值极小(< 0.001)。模型解释的销售方差比例为 33.1%(见 R 平方)。 我们可以通过在 “待检验变量(Variables to test)” 框中选择`income`、`HH.size`和`Age`,复现所有回归输出中报告的标准 F 检验。相关输出如下: Regression 1 - F-test 注意,在本示例中,“模型 1(model 1)” 是不含解释变量的回归。正如预期,模型 1 的解释方差为 0。F 检验比较模型 1 和模型 2 的拟合优度,并根据两个模型估计系数数量的差异进行调整。使用的检验统计量如下。R22是模型 2 的解释方差,R12是模型 1 的解释方差。n等于数据中的行数,k2(k1)是模型 2(模型 1)中估计的系数数量。 $$ \begin{eqnarray} F & = & \frac{(R^2_2 - R^2_1)/(k_2 - k_1)}{(1 - R^2_2)/(n - k_2 - 1)} \\\\ & = & \frac{(0.331 - 0)/(3 - 0)}{(1 - 0.331)/(200 - 3 - 1)} \\\\ & = & 32.325 \end{eqnarray} $$ 我们可以使用与 F 统计量 32.325 相关的 p 值评估原假设。也可以使用概率计算器计算临界 F 统计量。从下方输出可知,该值为 2.651。由于相关 p 值 <0.001,且计算的 F 统计量大于临界值(32.325> 2.651),我们拒绝所有系数均等于 0 的原假设。

回归系数的解释如下: - 在模型中其他变量保持不变的情况下,收入每增加 1000 美元,我们预期销售额平均增加 1.7754 美元。 - 在模型中其他变量保持不变的情况下,家庭规模每增加 1 人,我们预期销售额平均增加 22.1218 美元。 - 在模型中其他变量保持不变的情况下,家庭户主年龄每增加 1 岁,我们预期销售额平均增加 0.45 美元。 对于每个解释变量,可提出以下原假设和备择假设: - H0:解释变量 x 相关的系数等于 0 - Ha:解释变量 x 相关的系数不等于 0 `Income`和`HH.size`的系数均显著(p 值 < 0.05),即我们可以拒绝这两个系数的 H0。`Age HH`的系数不显著(p 值 > 0.05),即我们无法拒绝`Age HH`的 H0。我们得出结论:家庭户主年龄的变化不会导致销售额的显著变化。 我们也可以使用 t 值评估系数的原假设和备择假设。由于`Income`和`HH.size`的计算 t 值**大于**临界 t 值,我们拒绝这两个效应的原假设。可通过 “基础(Basics)” 菜单中的概率计算器获取临界 t 值。对于自由度为 196 的 t 分布(见上文 “平方和(Sum of Squares)” 表中的 “误差(Errors)” 行),临界 t 值为 1.972。由于备择假设是 “双侧(two.sided)”,我们需在概率计算器中输入 0.025 和 0.975 作为上下概率界。


### 示例 2:理想的回归数据 数据集`ideal`包含模拟数据,非常适合演示回归数据及残差的理想状态。该数据有 1000 个观测值,涉及 4 个变量。`y`是响应变量,`x1`、`x2`和`x3`是解释变量。下方图表可作为真实世界数据回归的基准。我们将使用 “模型> 线性回归(OLS)” 进行分析。首先,前往 “绘图(Plots)” 标签页,选择`y`作为响应变量,`x1`、`x2`和`x3`作为解释变量。 `y`、`x2`和`x3`大致呈正态分布,而`x1`大致呈均匀分布。没有异常值或严重偏态分布的迹象。

在相关性图中,响应变量与解释变量之间、解释变量彼此之间存在明显关联。注意,在实验中,关注的 x 变量应具有零相关性。但在历史数据中,这几乎不可能。图表下三角部分的散点图显示,变量之间的关系(近似)线性。

响应变量`y`与每个解释变量的散点图证实了相关性图的结论。拟合在散点图中的线足够灵活,可捕捉任何非线性关系。但这些线非常平直,表明线性模型可能适用。

六个残差图组成的面板看起来非常理想,这符合我们对模拟数据的预期。响应变量的真实值与回归预测值形成带随机散点的直线,即随着响应变量实际值的增加,模型预测值也增加。残差(即响应变量实际值与回归预测值的差异)无模式,随机分布在水平线周围。任何模式都表明模型对数据某些部分的预测效果优于(或差于)其他部分。如果残差与行顺序图中存在模式,我们可能需关注自相关性。残差仍均匀分布在水平轴附近。注意,自相关性主要是时间序列数据中需要关注的问题。Q-Q 图显示整齐的对角直线,表明残差呈正态分布。残差的直方图和密度图(绿色)与理论正态分布密度(蓝色线)的对比也证实了这一结论。

我们要讨论的最后一个诊断是残差与解释变量(或预测因子)的一组图表。没有迹象表明存在趋势或异方差性。这些图表中的任何模式都值得关注。也没有异常值(即远离主要数据点集群的点)。

由于诊断结果良好,我们可以从回归中得出推论。首先,模型整体显著:F 统计量的 p 值小于 0.05,因此我们拒绝回归中三个变量斜率均为零的原假设。其次,每个变量都具有统计显著性。例如,`x1`的 t 统计量 p 值小于 0.05,因此当`x2`和`x3`也在模型中时(即 “保持模型中其他变量不变”),我们拒绝`x1`斜率为零的原假设。 `x1`和`x3`的增加与`y`的增加相关,而`x2`的增加与`y`的减少相关。由于是模拟数据,系数的具体解释不太重要。但散点图中`x3`的增加似乎与`y`的减少相关,这一差异如何解释?提示:考虑相关性图。


### 示例 3:线性回归还是对数 - 对数回归? 线性回归和对数 - 对数回归均常用于商业数据。在本示例中,我们将从数据和残差中寻找证据,判断哪种模型设定更适合现有数据。 数据集`diamonds`包含 3000 颗钻石的价格信息。“数据> 管理” 页面提供了更完整的数据和变量描述。选择`price`作为响应变量,`carat`和`clarity`作为解释变量。在查看回归参数估计前,前往 “绘图” 标签页查看数据和残差。下方是模型中变量的直方图。`Price`和`carat`似乎向右偏斜。注意,偏斜方向由 “尾部” 位置决定。

在相关性图中,响应变量与解释变量之间存在明显关联。`price`与`carat`的相关性非常高(即 0.93)。钻石的`carat`与`clarity`的相关性显著且为负。

响应变量`price`与解释变量的散点图不如示例 2 中的`ideal`数据整洁。拟合在散点图中的线足够灵活,可捕捉非线性关系。`carat`的线似乎有一定曲率,点围绕该线的分布并非随机。实际上,在价格和克拉数较高时,点似乎呈扇形展开。不同`clarity`水平下,`price`的变动似乎不大。即便有变动,钻石价格似乎也随净度增加而下降。这一令人惊讶的结果我们将在下文详细讨论。

六个残差图组成的面板表现欠佳。响应变量的真实值与回归预测值形成 S 形曲线。在实际值和预测值较高时,点围绕线的分布更分散,这与`price`对`carat`的散点图一致。残差(即实际数据与回归预测值的差异)呈现明显模式,并非随机分布在水平轴周围。残差与行顺序图非常平直,表明自相关性不是问题。最后,与示例 2 中的`ideal`数据不同,此处 Q-Q 图的点在右极端明显偏离直线,表明残差不呈正态分布。残差的直方图和密度图显示比正态分布更尖的形状,也证实了这一结论。

我们要讨论的最后一个诊断是残差与解释变量(或预测因子)的一组图表。在残差对克拉数的图中,残差从左到右呈扇形展开。`clarity`对残差的散点图显示,净度较低时存在强负值异常值,净度较高时存在强正值异常值。

由于诊断结果不佳,我们**不应**从该回归中得出推论。对数 - 对数设定可能更合适。可通过 “数据> 可视化” 标签页快速检查对数 - 对数模型的有效性。在 “散点(Scatter)” 图中选择`price`作为 Y 变量,`carat`作为 X 变量。勾选 “log X” 和 “log Y” 框,生成下方图表。对数价格与对数克拉数之间的关系接近线性,这正是我们所期望的!

我们将对`price`和`carat`均进行(自然)对数(ln)转换,重新运行分析,看对数 - 对数设定是否更适合现有数据。可在 “数据> 转换” 中完成此转换。选择变量`price`和`carat`。从 “转换类型(Transformation type)” 下拉菜单中选择 “转换(Transform)”,从 “应用函数(Apply function)” 下拉菜单中选择 “Ln(自然对数)”。确保点击 “存储(Store)” 按钮,将新变量添加到数据集。注意,不能对`clarity`应用对数转换,因为它是分类变量。 在 “模型> 线性回归(OLS)” 中,选择`price_ln`作为响应变量,`carat_ln`和`clarity`作为解释变量。在查看回归参数估计前,前往 “绘图” 标签页查看数据和残差。下方是模型中变量的直方图。注意,`price_ln`和`carat_ln`不再右偏,这是好迹象。

在相关性图中,响应变量与解释变量之间仍存在明显关联。`price_ln`与`carat_ln`的相关性极高(即 0.93)。钻石的`carat_ln`与`clarity`的相关性显著且为负。

响应变量`price_ln`与解释变量的散点图现在整洁得多。`price_ln`对`carat_ln`的散点图中的线(基本)平直。尽管点在线周围呈轻微块状分布,但散射基本随机。我们不再看到`price_ln`和`carat_ln`值较高时的扇形展开。不同`clarity`水平下,`price_ln`的变动似乎稍大。但钻石的`price_ln`仍随`clarity`增加而下降,这不符合预期。我们将在下文讨论这一结果。

六个残差图组成的面板比线性模型好得多。响应变量的真实值与回归预测值(几乎)形成直线。尽管在实际值和预测值较高及较低时,线可能仍略呈曲线。残差更接近围绕水平线的随机散射。残差与行顺序图仍非常平直,表明自相关性不是问题。最后,Q-Q 图显示整齐的对角直线,与示例 2 中的`ideal`数据相同,表明残差现在呈正态分布。残差的直方图和密度图也证实了这一结论。

我们要讨论的最后一个诊断是残差与解释变量(或预测因子)的一组图表。与线性模型相比,残差更接近围绕水平线的随机散射。尽管`carat_ln`值较低(较高)时,残差可能稍高(较低)。

由于诊断结果现在好得多,我们可以更有信心地从该回归中得出推论。回归结果见 “摘要” 标签页。注意,`clarity`变量有 7 个系数,而`carat_ln`只有 1 个。为什么?查看数据描述(“数据> 管理”)可知,净度是分类变量,水平从 IF(最差净度)到 I1(最佳净度)。在应用回归等数值分析工具前,分类变量必须转换为一组虚拟(或指示)变量。每个虚拟变量表示某颗钻石是否具有特定净度等级(=1)或不具有(=0)。有趣的是,要捕捉 8 级净度变量的所有信息,我们只需 7 个虚拟变量。注意,没有净度等级 I1 的虚拟变量,因为回归中实际上不需要它。当一颗钻石**不是**净度 SI2、SI1、VS2、VS1、VVS2、VVS1 或 IF 时,我们从数据中可知它一定是净度 I1。 F 统计量表明,回归模型整体解释了`price_ln`的显著方差。F 统计量非常大,p 值极小(< 0.001),因此我们可以拒绝所有回归系数均等于 0 的原假设。模型解释的`price_ln`方差比例为 96.6%。钻石价格似乎比钻石需求更容易预测。 F 检验的原假设和备择假设可表述为: H0:所有回归系数均等于 0 Ha:至少有一个回归系数不等于 0 回归系数的解释如下: - 在模型中其他变量保持不变的情况下,克拉数每增加 1%,我们预期钻石价格平均增加 1.809%。 - 在模型中其他变量保持不变的情况下,与净度 I1 的钻石相比,我们预期净度 SI2 的钻石价格平均高 100×(exp (.444)-1) = 55.89%。 - 在模型中其他变量保持不变的情况下,与净度 I1 的钻石相比,我们预期净度 SI1 的钻石价格平均高 100×(exp (.591)-1) = 80.58%。 - 在模型中其他变量保持不变的情况下,与净度 I1 的钻石相比,我们预期净度 IF 的钻石价格平均高 100×(exp (1.080)-1) = 194.47%。 净度各等级的系数表明,`clarity`提高会增加钻石价格。那为什么净度与(对数)价格的散点图显示价格随净度增加而下降?差异在于,回归中我们可以确定一个变量(如净度)变化的影响,同时保持模型中其他变量不变(即克拉数)。较大、较重的钻石比较小的钻石更可能有瑕疵,因此查看散点图时,我们实际看到的不仅是净度提高对价格的影响,还有与净度呈负相关的克拉数的影响。在回归中,我们可以比较**相同大小**(即保持克拉数不变)的钻石在不同净度水平下对(对数)价格的影响。如果模型中没有(对数)克拉数,由于遗漏变量偏差,净度的估计效应可能不正确。实际上,从`price_ln`对`clarity`的回归中,我们会得出数据中净度最高的钻石(IF)比净度最低的钻石(I1)价格低 59.22% 的结论。显然这不是合理结论。 对于每个解释变量,可提出以下原假设和备择假设: H0:解释变量 X 相关的系数等于 0 Ha:解释变量 X 相关的系数不等于 0 该回归中的所有系数均高度显著。

### 报告 > Rmd 通过点击屏幕左下角的图标或按键盘上的`ALT-enter`,向*报告 > Rmd*添加代码以(重新)创建分析。 如果已创建图表,可使用`ggplot2`命令或`patchwork`进行自定义。详见下方示例和*数据 > 可视化*。 ```r result <- regress(diamonds, rvar = "price", evar = c("carat", "clarity", "cut", "color")) summary(result) plot(result, plots = "scatter", custom = TRUE) %>% wrap_plots(plot_list, ncol = 2) + plot_annotation(title = "Scatter plots") ``` ### 视频教程 将以下完整命令复制粘贴到 RStudio 控制台(即左下角窗口),按回车即可获取 Radiant 教程系列中线性回归模块使用的所有材料:
usethis::use_course("https://www.dropbox.com/sh/s70cb6i0fin7qq4/AACje2BAivEKDx7WrLrPr5m9a?dl=1")
回归的数据探索与前置检查(一) - 本视频展示如何使用 Radiant 在运行线性回归前探索和可视化数据 - 主题列表: - 查看数据 - 可视化数据 回归结果解读与预测(二) - 本视频解释如何解读回归结果并从线性回归模型计算预测值 - 主题列表: - 解释系数(数值变量和分类变量) - 解释 R 平方和调整后 R 平方 - 解释 F 检验结果 - 从回归模型进行预测 处理分类变量(三) - 本视频展示如何在线性回归模型中处理分类变量 - 主题列表: - 在 Radiant 中查看基准类别 - 更改基准类别 向回归模型添加新变量(四) - 本视频演示如何检验添加新变量是否能得到解释力显著提高的更优模型 - 主题列表: - 在 Radiant 中设置添加新变量的假设检验 - 解释 F 检验结果 - 将此 F 检验与回归摘要中的默认 F 检验进行比较 线性回归验证(五) - 本视频演示如何验证线性回归模型 - 主题列表: - 线性性(散点图,与前置检查中的相同) - 正态性检验(正态 Q-Q 图) - 多重共线性(VIF) - 异方差性 对数 - 对数回归(六) - 本视频演示何时及如何运行对数 - 对数回归 - 主题列表: - 使用自然对数函数转换偏态分布数据 - 解释对数 - 对数回归中的系数 ### 技术说明 #### 线性模型的系数解释 为说明回归模型中系数的解释,我们从以下方程开始: $$ S_t = a + b P_t + c D_t + \epsilon_t $$ 其中St是t时刻的单位销售量,Pt是t时刻的价格(美元),Dt是指示产品在某周是否有陈列的虚拟变量,ϵt是误差项。 对于价格等连续变量,在保持模型中其他变量不变的情况下,我们可通过对销售方程求关于P的偏导数,确定 1 美元变化的影响。 $$ \frac{ \partial S_t }{ \partial P_t } = b $$ 因此,b是价格每变化 1 美元对销售量的边际效应。由于D等虚拟变量不是连续的,我们不能使用微分,确定边际效应的方法略有不同。比较D=1和D=0时的销售水平,我们得到: $$ a + b P_t + c \times 1 - a + b P_t + c \times 0 = c $$ 对于线性模型,c是产品有陈列时对销售量的边际效应。 #### 半对数模型的系数解释 为说明半对数回归模型中系数的解释,我们从以下方程开始: $$ ln S_t = a + b P_t + c D_t + \epsilon_t $$ 其中ln S_t是t时刻销售量的(自然)对数。对于价格等连续变量,在保持模型中其他变量不变的情况下,我们可通过对销售方程求关于P的偏导数,确定小幅度变化(如 100 美元产品变化 1 美元)的影响。对于方程左侧,我们可使用链式法则: $$ \frac {\partial ln S_t}{\partial P_t} = \frac{1}{S_t} \frac{\partial S_t}{\partial P_t} $$ 通俗地说,变量自然对数的导数是该变量的倒数乘以该变量的导数。从上述线性模型的讨论中,我们知道: $$ \frac{ \partial a + b P_t + c D_t}{ \partial P_t } = b $$ 结合这两个方程: $$ \frac {1}{S_t} \frac{\partial S_t}{\partial P_t} = b \; \text{or} \; \frac {\Delta S_t}{S_t} \approx b $$ 因此,价格每变化 1 美元会导致销售量变化100×b%。注意,此近似仅适用于解释变量的小幅变化,且可能受使用的尺度影响(如价格以美分、美元或千美元计)。下文针对虚拟变量的方法更通用,也可应用于连续变量。 由于D等虚拟变量不是连续的,我们不能使用微分,再次比较D=1和D=0时的销售水平,得到StΔSt。为使左侧为St而非lnSt,我们对两边取指数,得到St=ea+bPt+cDt。当Dt从 0 变为 1 时,St的百分比变化为: $$ \begin{aligned} \frac {\Delta S_t}{S_t} &\approx \frac{ e^{a + b P_t + c\times 1} - e^{a + b P_t + c \times 0} } {e^{a + b P_t + c \times 0} }\\ &= \frac{ e^{a + b P_t} e^c - e^{a + b P_t} }{ e^{a + b P_t} }\\ &= e^c - 1 \end{aligned} $$ 对于半对数模型,100×(exp(c)−1)是产品有陈列时销售量的百分比变化。类似地,价格每增加 10 美元,在保持其他变量不变的情况下,我们预期销售量变化100×(exp(b×10)−1)。 #### 对数 - 对数模型的系数解释 为说明对数 - 对数回归模型中系数的解释,我们从以下方程开始: $$ ln S_t = a + b ln P_t + \epsilon_t $$ 其中lnPt是t时刻价格的(自然)对数。为简化,忽略误差项,我们可通过对两边取指数,将模型改写为乘法形式: $$ \begin{aligned} S_t &= e^a + e^{b ln P_t}\\ S_t &= a^* P^b_t \end{aligned} $$ 其中a∗=ea。对于价格等连续变量,我们可通过对销售方程求关于Pt的偏导数,得到边际效应: $$ \begin{aligned} \frac{\partial S_t}{\partial P_t} &= b a^* P^{b-1}_t\\ &= b S_t P^{-1}_t\\ &= b \frac{S_t}{P_t} \end{aligned} $$ 弹性的一般公式是∂St/∂Pt*Pt/St。将此信息添加到上述方程中,我们看到对数 - 对数回归估计的系数b可直接解释为弹性: $$ \frac{\partial S_t}{\partial P_t} \frac{P_t}{S_t} = b \frac{S_t}{P_t} \frac{P_t}{S_t} = b $$ 因此,价格每变化 1% 会导致销售量变化b%。 ### R 函数 有关 Radiant 中用于估计线性回归模型的相关 R 函数概述,请参见*模型 > 线性回归(普通最小二乘法)*。 `regress`工具中使用的核心函数包括`stats`包中的`lm`,以及`car`包中的`vif`和`linearHypothesis`。