JHU-数据科学-VII-笔记-全-

JHU 数据科学 VII 笔记（全）

001：回归简介 📈

在本节课中，我们将要学习回归分析的基本概念。回归是数据科学中最基础、最核心的工具之一，它用于探索变量之间的关系，并基于这种关系进行预测。我们将从回归的历史起源讲起，通过经典与现代的实例，理解回归分析能回答哪些问题，以及它在数据科学工作流程中的重要性。

回归的历史起源 🧬

上一节我们提到了回归的重要性，本节中我们来看看它的历史起源。回归的概念由弗朗西斯·高尔顿爵士提出，他同时发明了“回归”和“相关”这两个紧密相连的术语。

高尔顿的众多研究之一，是利用父母的身高来预测子女的身高。这个数据集因其历史意义，在本课程中将被多次探讨。我的同事兼课程讲师杰夫·利里指出，这个例子至今仍具有惊人的现实意义。

《美国人类遗传学杂志》上的一篇文章，将利用现代高通量生物信息学技术获得的遗传数据，与高尔顿在维多利亚时代收集的测量数据进行了比较。结果发现，现代方法并未显著优于高尔顿的古老方法。这既是对现代遗传学分析的一种有趣审视，也提醒我们弗朗西斯·高尔顿工作的超前性与重要性。

一个现代数据科学实例 🏀

了解了回归的历史背景后，本节中我们来看一个更贴近现代数据科学家日常工作的例子。这个例子来自“简明统计”博客，由本课程讲师杰夫·利里、罗杰·彭以及哈佛大学生物统计系的拉斐尔·伊里萨里共同撰写。

拉斐尔在一篇博文中分析了篮球运动员科比·布莱恩特的“独霸球权”行为。他绘制了科比出手投篮次数占比与球队单场净胜分差的关系图。

他对数据拟合了一条线性回归线，并在博文中提出了一些主张。例如，数据支持“如果科比减少独霸球权，湖人队会赢得更多比赛”这一说法。

更与本课程相关的是，他给出了一个具体的统计陈述：线性回归表明，科比出手占比每增加1%，会导致（湖人队）净胜分下降约1.6分。

并且，他给出了这个估计的标准误。如何生成并合理解读这样的陈述，如何遵循良好的统计实践（例如报告标准误），将是本课程的一大重点。

我们不仅要学习如何进行分析，也要思考分析本身是否合理。例如，这个分析是否完整？拉法是否考虑了所有相关因素？阅读那篇博文及其评论，会涉及许多与回归分析相关的主题。

回归分析能回答的问题 ❓

看过了现代应用，让我们回到高尔顿的经典例子，系统性地探讨回归分析能帮助我们回答哪些类型的问题。

以下是回归分析可以解决的核心问题：

  • 预测：利用父母身高预测子女身高。
  • 建立简约关系：寻找父母与子女身高之间简洁、易于描述的均值关系。回归擅长在保证解释力的前提下提供简约的模型，而许多机器学习方法虽然预测可能更精确，但模型往往复杂，不易解释和生成新知识。
  • 量化未解释的变异：讨论回归模型无法解释的变异，即所谓的残差变异。
  • 比较影响因素：量化如基因信息等因子，在父母身高之外对解释子女身高的额外贡献。这与之前提到的遗传学论文研究的问题相关。
  • 进行统计推断：将基于样本数据的结论，推广到更大总体所需的条件和假设。这是一个深刻的主题，称为统计推断。本课程将把推断工具应用到回归分析中。
  • 理解“回归”本质：回答高尔顿提出的一个有趣问题——“均值回归”或“回归平庸”。即：为什么非常高父母的孩子通常也高，但会比父母矮一点？为什么非常矮父母的孩子通常也矮，但会比父母高一点？这就是著名的“均值回归”思想。

总结 📝

本节课中我们一起学习了回归分析的简介。我们回顾了回归由弗朗西斯·高尔顿创立的历史，并通过一个现代篮球数据分析的实例，看到了回归在数据科学中的实际应用。最后，我们系统地列出了回归分析能够回答的各类核心问题，包括预测、建立简约模型、分析残差、比较变量贡献、进行统计推断以及理解“均值回归”现象。这些概念构成了我们深入学习回归模型的基础。

002：简单最小二乘法简介 📊

在本节课中，我们将学习回归分析的基础——最小二乘法。我们将通过弗朗西斯·高尔顿的历史数据，直观地理解如何寻找数据的“中心”，并以此作为预测的起点。

数据背景与准备

首先，我们来看弗朗西斯·高尔顿在1885年首次使用的数据。弗朗西斯·高尔顿是统计学史上一位非常有趣的人物，他的个性和贡献都颇具特色。

为了进行后续分析，我们需要在R中安装并使用 UsingR 包。这个包收录了《R语言入门统计学》一书中使用的多个数据集，非常适合我们的学习。

数据的边际分布

在探讨变量间关系之前，我们先分别查看父母身高和子女身高的边际分布。这意味着我们暂时忽略父母与子女之间的关联，单独分析每一组数据。

以下是使用 ggplot2 绘制的两个直方图：

  • 左侧直方图展示了子女的身高分布（单位：英寸，范围60-75英寸），Y轴表示每个身高区间内的子女数量。
  • 右侧蓝绿色直方图展示了父母的身高分布。这里对父母身高进行了性别校正（例如，将女性身高乘以1.08），并同样忽略了其与子女身高的关联。

我们将利用这些分布引入最小二乘法的概念，之后再探讨双变量间的关联。

寻找数据的“中心”

现在，我们暂时只关注子女身高，先不考虑用父母身高进行预测。我们想找到一个最佳预测值，在没有任何其他信息的情况下预测任意一个孩子的身高。这个最佳预测值很可能就是数据的“中心”。

那么，如何定义这个“中心”呢？

数学定义

我们将通过数学符号来精确定义。如果你对数学推导不感兴趣，可以稍作等待，我们随后会通过一个可视化的物理实验来直观解释。

我们用 Y_i 表示第 i 个孩子的身高，其中 i 从1到 N（在本数据集中 N=928）。我们定义的“中心”是一个值 μ，它使得所有观测值 Y_i 与 μ 之间的距离平方和最小。

用公式表示，即寻找 μ，使得以下值最小：
∑ (Y_i - μ)^2

这个定义与物理学中的质心概念相关。如果把之前直方图中的每个柱子看作有质量的物理实体，那么使这个系统平衡的支点就是质心，也就是我们通过最小二乘法找到的“中心”。



你可能已经猜到，这个“中心”就是数据的平均值。事实确实如此，我们稍后会进行证明。当然，如果你不想深究证明过程，也完全可以理解核心思想。


可视化实验

为了更直观地感受，我们使用RStudio的 manipulate 函数来做一个实验，手动寻找这个“质心”。

我已经加载了 manipulate 包和 Galton 数据集，并定义了一个绘制直方图的函数。这个函数接受一个参数 mu，即我们尝试寻找的中心值。

运行 manipulate 命令后，会出现一个滑块。上方会显示当前 mu 的值以及对应的均方误差。均方误差就是所有数据点与当前 mu 值距离平方的平均值（即总和除以N）。

以下是实验过程：

  • 当我移动滑块时，可以看到随着 mu 接近直方图的中心区域，均方误差逐渐减小。
  • 在某个位置（大约68.34英寸附近），均方误差达到最小。
  • 如果将滑块移动到远离中心的位置（比如很高或很低），均方误差会再次变大。

通过这个实验，你可以亲身体验到，使观测数据点与 mu 之间距离平方和最小的那个点，正是平衡这个直方图系统的“质心”。虽然我们知道答案就是身高的平均值，但实验能帮助我们直观地验证这一结论。

代码与结论

这里再次提供了实验的R代码。需要提醒的是，本专业所有课程的讲义都以R Markdown格式提供，你可以通过名为 index.Rmd 的文件获取每节课的所有代码。

最后，我们重申结论：最小二乘估计值就是样本的均值。下图展示了均值在数据分布中的位置，并附上了生成该图的 ggplot2 命令。


总结

本节课中，我们一起学习了最小二乘法的基本思想。我们从分析数据的边际分布开始，定义了通过最小化平方误差和来寻找数据“中心”的方法，并通过数学和可视化实验两种方式，验证了这个“中心”就是数据的平均值。这为后续学习如何利用一个变量（如父母身高）预测另一个变量（如子女身高）的回归模型奠定了重要基础。

003：证明样本均值是最小化平方和的最优解 📊

在本节课中，我们将学习一个重要的数学证明：为什么样本均值是使残差平方和最小的最优解。这个结论是回归模型和许多统计方法的基础。

上一节我们介绍了回归模型的基本概念，本节中我们来看看如何从数学上证明样本均值的最优性。

证明概述

我们想要证明，对于给定的观测数据点 ( Y_i )，样本均值 ( \bar{Y} ) 是使以下平方和函数最小的唯一值：

[
\sum_{i=1}^{n} (Y_i - \mu)^2
]

其中 ( \mu ) 是我们想要优化的参数。

证明步骤

以下是证明的核心步骤。

首先，我们在表达式中同时减去和加上样本均值 ( \bar{Y} )，这相当于加零，不会改变原式的值：

[
\sum (Y_i - \mu)^2 = \sum [(Y_i - \bar{Y}) + (\bar{Y} - \mu)]^2
]

接着，我们展开平方项：

[
\sum [(Y_i - \bar{Y})^2 + 2(Y_i - \bar{Y})(\bar{Y} - \mu) + (\bar{Y} - \mu)^2]
]

现在，让我们仔细观察交叉项 ( 2(Y_i - \bar{Y})(\bar{Y} - \mu) )。由于 ( (\bar{Y} - \mu) ) 不依赖于索引 ( i )，我们可以将其提到求和符号外面：

[
2(\bar{Y} - \mu) \sum (Y_i - \bar{Y})
]

根据样本均值的定义 ( \bar{Y} = \frac{1}{n} \sum Y_i )，我们可以推导出：

[
\sum (Y_i - \bar{Y}) = \sum Y_i - n\bar{Y} = n\bar{Y} - n\bar{Y} = 0
]

因此，整个交叉项为零。于是，原始平方和可以简化为：

[
\sum (Y_i - \mu)^2 = \sum (Y_i - \bar{Y})^2 + \sum (\bar{Y} - \mu)^2
]

请注意，( \sum (\bar{Y} - \mu)^2 ) 是一系列平方项的和，它总是大于或等于零。如果我们舍弃这一项，原式只会变小或保持不变：

[
\sum (Y_i - \mu)^2 \geq \sum (Y_i - \bar{Y})^2
]

这个不等式对任何 ( \mu ) 都成立，并且等号仅在 ( \mu = \bar{Y} ) 时取得。因此，样本均值 ( \bar{Y} ) 是使残差平方和最小的唯一解。

总结

本节课中我们一起学习了如何证明样本均值是最小化平方和的最优解。我们通过代数变换，展示了对于任何其他值 ( \mu )，其对应的平方和都大于或等于使用样本均值时的平方和。这个结论是理解最小二乘法原理的关键基础。

004：3_介绍性的数据示例

在本节课中，我们将通过一个具体的数据示例，学习回归分析的基本思想。我们将使用父母与子女身高的数据，探讨如何用一条直线来描述两者之间的关系，并介绍“最小二乘法”这一核心概念。

数据可视化与问题

到目前为止，我们尚未使用父母的身高数据。处理此类数据时，首先要做的是绘制子女身高与父母身高的散点图。这里我使用了 ggplot2 包。

这张图存在许多缺陷。一个主要缺陷是数据点过度重叠。在许多特定的坐标点上，存在多个不同的父母-子女身高配对（X，Y值）。

这里我提供了一张更好的图。😡 图中点的大小代表了在该特定（X，Y）位置上的父母-子女配对数量。颜色越浅（接近白色）代表频率接近40，而颜色越偏蓝或越深（相对而言）则代表频率向10或更低递减。因此，点的大小和颜色代表了该特定（X，Y）配对中父母-子女组合的频率。这是一张更好的图，因为它没有丢失这些信息。

我认为这张图还有另一个缺陷，例如，我没有在两个坐标轴上标注单位（英寸）。在坐标轴上标注单位是良好的实践，所以这也是此图的一个小缺点。如果你想查看代码，可以在R Markdown文件中找到。


回归直线的初步构想

现在，假设我们想用父母的身高来解释子女的身高，并且我们想用一条直线来实现。为了使问题简化，我们暂时强制这条直线通过原点。具体来说，它必须通过点（0，0），也就是说它的y轴截距为0。

那么，这条直线就是 y = βx，其中β是斜率。y = βx 定义了一条通过原点的直线。为了找到最佳直线，我们只需要找到最佳的斜率β。

以下是我们的做法：我们希望找到斜率β，使得观测数据点（Yi）与直线上拟合点（Xiβ）之间的垂直距离的平方和最小。我们将这个距离平方后求和。这直接类比于我们之前几页幻灯片中计算的最小二乘均值。这本质上就是以原点为支点，选择一条使所有点到该直线垂直距离的平方和最小的直线。

我们将使用RStudio的 manipulate 函数来实验并尝试找到这条直线。

通过原点的回归在解释某些问题时很有用，因为我们只有一个参数（斜率），而不是两个参数（斜率和截距）。但强制回归线通过点（0，0）通常是不可取的。一个简单的解决方法是：从父母身高和子女身高中分别减去各自的均值，这样（0，0）点就位于数据的中心。这将使这个解决方案更合理。我们稍后会讨论这与同时拟合斜率和截距的真实回归之间的关系。

核心概念图解

让我用一张图来说明其中一些概念。

这里有一个散点图，Y轴和X轴上都有一些数据。我想用我的X变量来预测Y变量。红色十字准线是我的数据点。我们正在做的“通过原点的回归”将点（0，0）视为支点，然后尝试找到最佳直线。

以这条线为例。对于每个观测点Yi，它会计算垂直距离。这个高度是Yi，这个长度是Xi，直线上的那个点是Xiβ。所以这里的距离就是 Yi - Xiβ。如果我们将其平方，就得到平方距离。

因此，通过原点的回归试图找到所有观测高度与拟合线之间垂直距离的平方和，并使这个总和最小。因为记住，通过原点的直线由一个简单方程 y = βx 定义，你只需要一个参数，即斜率。

现在，这条线看起来不会很好。如果你观察，这条线实际上应该穿过这里的某个地方，所以“回归到原点”似乎不太合理。我们通过对数据“中心化”（减去均值）所做的，就是将原点设置在数据的中心。通过减去均值，现在（0，0）点位于数据的中间。这基本上是重新定位坐标轴，使（0，0）点正好位于数据中央。现在，只考虑斜率来寻找穿过数据的回归线就显得合理多了。我们稍后会发现，这与同时拟合截距和斜率的解决方案是等价的。

让我们尝试用RStudio的 manipulate 函数来做这个。我认为因为这是接下来几节课的中心主题之一，我们会反复讨论这些要点，它们非常基础。

动手实验与R代码实现

我不展示运行代码的过程，因为你可以从R Markdown文件中获取代码。我希望在专业课程的这一阶段，你已经熟悉运行R代码。

这是我的图。在图上，我有一个回归线的特定值。请注意，子女身高的中心在这里是0，因为我减去了均值；父母身高的中心也是0。所以我有一个支点正好在中间的（0，0）点。

我还想指出，在这里我给出了斜率 β = 0.6 和均方误差。均方误差再次是Y数据点减去父母身高乘以这个因子0.6后的平方和。它计算每个点的垂直距离，将其平方后求和。同样，由于在任何（X，Y）组合处都有多个点，你可以考虑将这个平方距离乘以该特定组合的适当点数（因为过度绘图）。

在（0，0）点，我们的直线将围绕该点旋转，因为我们只拟合必须通过原点的直线，只考虑斜率。

让我们开始操作，看看均方误差如何随斜率变化。0.6看起来不错。让我们移动到一个不那么好的值。均方误差变得糟糕得多。现在把它一直调高，到了5.8。你可以看到均方误差在降低。当斜率降到大约0.6左右时，效果相当好。0.64对应5.0，然后0.62对应5.002，所以又上升了。看起来最佳斜率大约是0.64。

有趣的是，斜率为1并不好。斜率为1意味着你想预测子女身高，就直接用父母身高。但在这个案例中，我们必须乘以一个因子（0.6）才能得到比直接用父母身高更好的预测。这就是斜率为我们所做的。

这里我给出了 manipulate 代码，同样，所有代码都在R Markdown文件中。所以不要从幻灯片上重新输入，直接从R Markdown文件中获取实际文本。

既然我们已经用 manipulate 找到了斜率，我将展示如何在R中快速完成这个操作。函数 lm() 用于拟合线性模型。


这里是我使用 lm() 的代码。我的结果变量（y值）是中心化后的子女身高（减去了均值）。这里是中心化后的父母身高（减去了父母的均值）。-1 表示去掉截距，因为我们讨论的是通过原点的回归，不考虑Y轴截距。然后我指定要拟合的数据集是 Galton。它给出的斜率是 0.646。

我们将详细讨论这一点，并在课程的其余部分深入探讨。

最终拟合线与总结

在本讲的最后一张幻灯片中，我给出了拟合的直线。由于我们强制模型通过中心化子女身高和中心化父母身高的（0，0）点，它必须通过子女身高和父母身高的均值。我已经重新调整了绘图，使那个（0，0）点回到了原始位置。

所以，这里是我们斜率为 0.646 的直线。在接下来的课程中，我们将讨论如何得到这些结果、背后的动机以及我们可以用这条拟合线做的所有事情。我们可能会花接下来几节课的时间来讨论这个。

欢迎来到回归分析入门。实际上，我们在第一讲中已经涵盖了很多材料，我非常期待在接下来的几周里与大家一起学习。

本节课总结

在本节课中，我们一起学习了：

  1. 数据可视化的重要性：通过改进的散点图展示数据重叠问题。
  2. 回归的基本思想：用一条直线（y = βx）描述变量间关系。
  3. 最小二乘法原理：通过最小化垂直距离的平方和（∑(Yi - Xiβ)²）来寻找最佳拟合直线。
  4. 中心化的作用：通过减去均值将数据原点移至中心，使“通过原点的回归”更具解释性。
  5. R语言工具：使用 manipulate 包交互式探索斜率，并使用 lm() 函数快速拟合模型，得到斜率估计值 0.646。

我们为理解回归分析奠定了坚实的基础，接下来的课程将在此基础上深入展开。

005：符号与背景 📚

在本节课中，我们将学习回归分析所需的基本符号和背景知识。这些概念是理解后续回归模型的基础。

概述

本节将回顾统计学中的核心符号和定义，包括均值、方差、标准化以及协方差与相关系数。这些内容是回归分析的基石。

基本符号

在回归分析中，我们经常使用下标来表示有序的数据集合。

例如，我们可以用 x₁, x₂, ..., xₙ 来表示 n 个数据点。假设数据为 [1, 2, 5]，那么：

  • x₁ = 1
  • x₂ = 2
  • x₃ = 5
  • n = 3

字母本身没有特殊含义，我们同样可以用 y₁, y₂, ..., yₙ 来表示另一组数据。

一个重要的约定是：

  • 我们通常用希腊字母（如 μ）来表示未知的总体参数。
  • 我们用普通字母（如 x̄）来表示可以观测或计算的样本统计量。

样本均值

上一节我们介绍了数据索引的符号，本节中我们来看看如何计算样本均值。

样本均值是数据集中所有观测值的平均值，计算公式为：

x̄ = (1/n) * Σᵢ₌₁ⁿ xᵢ

其中，Σ 符号表示求和。例如，对于数据 [1, 2, 5]，其均值为 (1+2+5)/3 ≈ 2.67。

一个重要的性质是，如果我们从每个数据点中减去均值，得到的新数据集（称为中心化数据）其均值将为零。

令 x̃ᵢ = xᵢ - x̄，则新数据集 {x̃₁, x̃₂, ..., x̃ₙ} 的均值为 0。

此外，样本均值是“最小二乘”意义上的中心点，即它是使 Σᵢ₌₁ⁿ (xᵢ - μ)² 最小的 μ 值。

方差与标准差

在了解了数据的中心位置（均值）后，我们需要衡量数据的离散程度，这就是方差和标准差。

样本方差 s² 衡量了数据点围绕均值的平均平方偏差，计算公式为：

s² = [1/(n-1)] * Σᵢ₌₁ⁿ (xᵢ - x̄)²

样本标准差 s 是方差的平方根：

s = √s²

使用标准差的好处在于，它的单位与原始数据 x 的单位一致，而方差的单位是原始单位的平方。

类似于中心化，我们可以对数据进行缩放。如果将每个观测值除以标准差 s，得到的新数据集其标准差将为 1。

数据标准化

结合中心化和缩放，我们可以对数据进行标准化（或称为归一化）。

以下是标准化的步骤：

  1. 中心化：将每个数据点减去均值（xᵢ - x̄），得到均值为0的数据。
  2. 缩放：将中心化后的每个数据点除以标准差 s。

令 zᵢ = (xᵢ - x̄) / s，则新数据集 {z₁, z₂, ..., zₙ} 具有以下性质：

  • 经验均值为 0。
  • 经验标准差为 1。

标准化数据以零为中心，其数值表示该点距离均值有多少个标准差。例如，一个标准化后的值为 2，意味着该原始数据点比均值高出 2 个标准差。标准化使得不同尺度的数据集之间具有可比性。

协方差与相关系数

现在，我们进入回归分析最核心的概念之一：衡量两个变量之间关系的指标。

想象我们有两组一一对应的数据，例如每个人的体重指数 (x) 和血压值 (y)。我们可以有意义地绘制散点图来观察它们的关系。

样本协方差 Cov(x, y) 衡量了 x 和 y 共同变化的趋势，计算公式为：

Cov(x, y) = [1/(n-1)] * Σᵢ₌₁ⁿ (xᵢ - x̄)(yᵢ - ȳ)

然而，协方差的单位是 x 的单位乘以 y 的单位，不便于直接比较。因此我们引入相关系数 r，它是标准化后的协方差，是一个无单位的量：

r = Cov(x, y) / (sₓ * sᵧ)

其中 sₓ 和 sᵧ 分别是 x 和 y 的标准差。

相关系数的性质

关于相关系数，有以下重要事实：

  • 取值范围：相关系数 r 的值介于 -1 和 1 之间。
  • 极端情况：
    • 当 r = 1 时，所有数据点完全落在一条正斜率的直线上，表示完全正线性相关。
    • 当 r = -1 时，所有数据点完全落在一条负斜率的直线上，表示完全负线性相关。
  • 含义：相关系数衡量了 x 和 y 之间线性关系的强度和方向。
    • r 的绝对值越接近 1，线性关系越强。
    • r 接近 0，表示没有线性关系（但可能存在其他类型的关系）。
  • 核心地位：相关系数是线性回归分析中的核心度量。

总结

本节课中我们一起回顾了回归分析所需的基本符号和统计背景。我们学习了如何表示数据集合，计算样本均值和方差，以及通过中心化、缩放和标准化来调整数据。最后，我们深入探讨了衡量两个变量线性关系的核心工具——协方差和相关系数，并了解了相关系数的关键性质。掌握这些概念是理解后续线性回归模型原理和应用的基础。

006：线性最小二乘法

在本节课中，我们将要学习线性回归的核心估计方法——最小二乘法。我们将了解如何通过最小化误差平方和，找到数据的最佳拟合直线。

引言与数据背景

上一节我们介绍了回归模型的基本概念，本节中我们来看看如何具体计算回归线。

再次考虑高尔顿数据中父母身高与子女身高的散点图。图中圆圈的大小代表了特定（X，Y）组合出现的频率。我们希望利用父母身高来解释子女身高，并将通过线性回归来实现这一目标。

模型与符号定义

我们将使用上一讲中建立的符号体系。

  • 令 Y_i 表示第 i 个孩子的身高。
  • 令 X_i 表示第 i 个父母的身高。

我们希望找到一条最佳直线，其形式为：
子女身高 = 截距 + 斜率 × 父母身高

用公式表示为：
Y_i = β_0 + β_1 * X_i

其中，β_0（截距）和 β_1（斜率）是我们未知但希望求得的参数。

最小二乘准则

我们需要一个标准来定义“最佳”拟合线。一个著名的标准就是最小二乘准则。

其基本思想是：我们希望最小化数据点（即子女身高）与拟合线上对应点之间的垂直距离的平方和。

用公式可以写为：
最小化 Σ [ Y_i - (β_0 + β_1 * X_i) ]²

准则的直观理解

让我们通过图示来理解最小二乘准则的含义。

假设我们有一个数据集。最小二乘准则寻找最佳回归线的过程如下：

  1. 它为数据拟合一条直线。
  2. 对于每个数据点（例如点 (X_1， Y_1)），它计算该点的实际 Y 值（Y_1）与直线上对应 X 值的预测 Y 值（β_0 + β_1 * X_1）之间的垂直距离。
  3. 将这个距离（可能为正或负）进行平方，以消除符号影响。
  4. 对数据集中的所有点重复此过程，并将所有平方距离相加。

最终，最小二乘法通过调整 β_0 和 β_1 的值，来最小化这个总的平方误差和。每个数据点对误差的贡献是相等的。

最小二乘解

当我们执行这个模型拟合时，解是什么？

我们寻找的直线是 Y = β_0 + β_1 * X，通过一系列点 (X_i， Y_i) 进行拟合，其中 Y_i 是结果变量。我们在 β_0 和 β_1 上加上“帽子”符号（^）来表示它们的估计值。

解的具体形式如下：

  • 斜率估计值 β̂₁ 的计算公式为：
    β̂₁ = Cor(Y， X) * (S_y / S_x)
    其中 Cor(Y， X) 是 Y 与 X 的相关系数，S_y 是 Y 的标准差，S_x 是 X 的标准差。

  • 截距估计值 β̂₀ 的计算公式为：
    β̂₀ = Ȳ - β̂₁ * X̄
    其中 Ȳ 是 Y 的均值，X̄ 是 X 的均值。

解的性质与推论

了解这个结果后，我们可以推导出几个重要性质：

以下是关于最小二乘解的几个关键点：

  1. 单位：斜率 β̂₁ 的单位是 Y 的单位除以 X 的单位。这可以从公式中看出：相关系数是无单位的，S_y 带有 Y 的单位，S_x 带有 X 的单位。
  2. 穿过均值点：拟合直线总是穿过数据的中心点 (X̄， Ȳ)。从截距公式重新排列即可得到 Ȳ = β̂₀ + β̂₁ * X̄。
  3. 角色互换：如果我们将 Y 和 X 的角色互换，将 X 作为结果变量，Y 作为预测变量，那么得到的斜率将是 Cor(X， Y) * (S_x / S_y)。这与以 Y 为结果、X 为预测变量时得到的斜率不同。
  4. 数据中心化：如果先将数据中心化（即每个 X_i 减去 X̄，每个 Y_i 减去 Ȳ），然后强制回归线穿过原点（0， 0）进行拟合，得到的斜率是相同的。
  5. 数据标准化：如果先将数据标准化（即中心化并缩放，使标准差为1），那么斜率就等于相关系数 Cor(Y， X)，因为此时 S_y 和 S_x 都等于1。

总结

本节课中我们一起学习了线性回归的核心估计方法——最小二乘法。我们定义了通过最小化垂直距离平方和来寻找最佳拟合线的准则，并给出了斜率和截距估计值的具体计算公式。我们还讨论了该解的几个重要性质，包括其单位、必然穿过数据均值点，以及在数据经过中心化或标准化后的表现。理解最小二乘法是掌握线性回归建模的基础。

007：线性最小二乘法编码示例 📊

在本节课中，我们将通过一个具体的R语言编码示例，深入理解线性最小二乘法的计算过程。我们将使用高尔顿父母与子女身高数据集，手动计算回归系数，并与R内置的lm函数结果进行对比验证。

数据准备与可视化

首先，我们需要加载并查看将要分析的数据。这里使用的是高尔顿的父母身高与子女身高数据集。

以下是初始步骤，包括加载数据和绘制散点图。

# 绘制父母身高与子女身高的散点图
library(ggplot2)
# 假设数据已加载到名为 Galton 的数据框中
g <- ggplot(Galton, aes(x = parent, y = child)) + geom_point()
print(g)

上图展示了父母身高（x轴）与子女身高（y轴）之间的关系。

手动计算回归系数

上一节我们介绍了回归模型的基本概念，本节中我们来看看如何手动计算最小二乘估计的斜率和截距。我们将根据公式进行计算，并与R的lm函数结果进行比较。

定义变量以简化后续计算：

y <- Galton$child  # 结果变量：子女身高
x <- Galton$parent # 预测变量：父母身高

根据公式，斜率 β₁ 和截距 β₀ 的计算方法如下：

  • 斜率 β₁ 的公式为：β₁ = cor(y, x) * (sd(y) / sd(x))
  • 截距 β₀ 的公式为：β₀ = mean(y) - β₁ * mean(x)

以下是具体的计算代码：

# 手动计算斜率和截距
beta1 <- cor(y, x) * sd(y) / sd(x)
beta0 <- mean(y) - beta1 * mean(x)
c(beta0, beta1)

使用 lm 函数验证结果

为了验证我们手动计算的正确性，现在使用R内置的线性模型函数 lm 进行拟合。

以下是使用 lm 函数进行线性回归的代码：

# 使用lm函数拟合线性模型
fit <- lm(y ~ x)
coef(fit)

运行后，你会发现 coef(fit) 输出的截距和斜率与我们手动计算的 c(beta0, beta1) 结果完全一致（例如23.94和0.6463），这证实了我们的公式和计算是正确的。

反转X与Y的关系

一个有趣的现象是，如果我们交换自变量和因变量的角色，公式依然成立，但形式会发生变化。

现在，我们将子女身高作为预测变量（X），父母身高作为结果变量（Y）来建立模型。

以下是反转关系后的系数计算：

# 当x为子女身高，y为父母身高时
beta1_rev <- cor(y, x) * sd(x) / sd(y) # 注意公式中sd的位置交换了
beta0_rev <- mean(x) - beta1_rev * mean(y)
c(beta0_rev, beta1_rev)

# 使用lm函数验证
fit_rev <- lm(x ~ y)
coef(fit_rev)

对比 c(beta0_rev, beta1_rev) 和 coef(fit_rev)，两者结果依然相同，这证明了公式的普适性。



中心化数据与过原点回归

在讲座中提到的另一个要点是，如果先将Y和X中心化（减去各自的均值），那么过原点回归（即强制截距为0的回归）得到的斜率，与包含截距的普通线性回归斜率是相同的。

以下是验证这一点的计算步骤：

首先，对数据进行中心化处理：

yc <- y - mean(y) # 中心化的Y
xc <- x - mean(x) # 中心化的X

然后，计算过原点回归的斜率。其公式为：β₁_origin = sum(yc * xc) / sum(xc^2)

beta1_origin <- sum(yc * xc) / sum(xc^2)
beta1_origin

最后，与包含截距的线性回归模型的斜率进行比较：

# 普通线性回归的斜率（第二个系数）
coef(lm(y ~ x))[2]
# 使用中心化数据并强制截距为0的回归
coef(lm(yc ~ xc - 1))

你会发现，beta1_origin、coef(lm(y ~ x))[2] 和 coef(lm(yc ~ xc - 1)) 这三个值是完全相同的。代码中的 -1 即表示在 lm 公式中移除截距项。

标准化与相关系数

接下来，我们验证另一个重要结论：如果将Y和X标准化（均值为0，标准差为1），那么回归斜率就等于二者的相关系数。

以下是标准化的过程及验证：

# 标准化Y和X
yn <- (y - mean(y)) / sd(y)
xn <- (x - mean(x)) / sd(x)

# 验证标准差是否为1
sd(yn)
sd(xn)

# 比较相关系数与回归斜率
c(cor(y, x), cor(yn, xn), coef(lm(yn ~ xn))[2])

执行后，你会看到输出的三个数值完全相同。这证实了对于标准化后的数据，皮尔逊相关系数 cor(yn, xn) 就等于回归模型的斜率估计值。


在图形中添加回归线

最后，我们将学习如何在使用 ggplot2 绘制的散点图上添加回归线，这能更直观地展示变量间的关系。

以下是创建带回归线的增强散点图的代码。我们使用气泡图，其中点的大小和颜色代表该（父母身高，子女身高）组合出现的频数。

# 创建基础散点图（气泡图表示频数）
g <- ggplot(Galton, aes(x = parent, y = child, size = ..count.., color = ..count..)) +
  geom_point(stat = "sum")

# 添加线性回归平滑层
g + geom_smooth(method = "lm", formula = y ~ x)

geom_smooth(method = "lm") 层会自动拟合数据并绘制回归线。默认情况下，它还会添加回归线的置信区间带（图中灰色阴影区域），这很好地展示了估计的不确定性，我们将在后续课程中讨论其生成原理。

总结

本节课中我们一起学习了线性最小二乘法的实际编码应用。我们使用R语言手动计算了回归系数，并通过多种方式验证了其正确性，包括与 lm 函数对比、验证中心化与标准化后的性质。最后，我们还掌握了如何在 ggplot2 图形中直观地添加回归线。通过这些练习，你应当对最小二乘估计的计算和验证有了更扎实的理解。

008：线性回归斜率公式推导 🧮

在本节课中，我们将要学习线性回归斜率估计公式的数学推导。我们将证明，斜率估计值 β1_hat 等于 y 和 x 的相关系数乘以 y 的标准差再除以 x 的标准差。这个推导过程会涉及一些数学运算，但我们会尽量用简单直白的方式讲解，让初学者也能跟上思路。

从过原点回归开始

在进入完整的线性回归推导之前，我们先来看一个更简单的情况：过原点回归。这意味着我们拟合一条形式为 y = xβ 的直线，且这条直线强制通过原点（0,0）。

假设我们有 n 个数据点：(x1, y1), (x2, y2), ..., (xn, yn)。我们的目标是找到一个 β 值，使得以下最小二乘准则最小化：

∑ (yi - xi * β)^2，其中求和 ∑ 从 i=1 到 n。

我们设使这个表达式最小的解为 β_hat。为了找到它，我们对准则进行一些代数变换。

首先，我们在表达式中巧妙地加上一个 0：

∑ [ (yi - xi * β_hat) + (xi * β_hat - xi * β) ]^2

展开这个平方项，我们得到三项的和：

  1. ∑ (yi - xi * β_hat)^2
  2. -2 * ∑ [ (yi - xi * β_hat) * (xi * β_hat - xi * β) ]
  3. + ∑ (xi * β_hat - xi * β)^2

现在，关键的一步是观察第二项。如果我们选择的 β_hat 能使这一项等于 0，那么对于任何其他的 β 值，其对应的最小二乘准则值都会大于或等于 β_hat 对应的值。这就证明了 β_hat 是最小值点。

因此，我们令第二项为 0 并求解 β_hat：

∑ [ (yi - xi * β_hat) * xi * (β_hat - β) ] = 0

由于 (β_hat - β) 是一个不依赖于 i 的常数，我们可以将其从求和号中提出并消去。于是，方程简化为：

∑ (yi - xi * β_hat) * xi = 0

解这个方程，我们得到过原点回归的斜率估计公式：

β_hat = (∑ yi * xi) / (∑ xi^2)

一个特例：均值回归

为了验证这个公式，让我们考虑一个特例：假设所有的 xi 都等于 1。

此时，我们的模型 y = xβ 变成了 y = 1 * β，也就是 y = β。这被称为均值回归，因为我们试图用一个常数 β 来拟合所有的 y 值。

将 xi = 1 代入我们的公式：

β_hat = (∑ yi * 1) / (∑ 1^2) = (∑ yi) / n = y_bar

这个结果与我们之前学到的知识一致：在均值回归中，使平方和最小的估计值就是 y 的样本均值 y_bar。

推广到一般线性回归

上一节我们介绍了过原点回归，本节中我们来看看完整的一般线性回归模型。模型形式为：

y = β0 + β1 * x

我们的目标是最小化最小二乘准则：

∑ (yi - β0 - β1 * xi)^2

我们的推导将利用过原点回归的结论，使过程变得简单。

首先，假设 β1 是已知的固定值。那么，我们可以将准则重写为：

∑ [ (yi - β1 * xi) - β0 ]^2

令 y*_i = yi - β1 * xi，则上式变为 ∑ (y*_i - β0)^2。这正是一个均值回归问题！我们知道，使这个表达式最小的 β0 就是 y* 的均值：

β0_hat = (∑ y*_i) / n = (∑ (yi - β1 * xi)) / n = y_bar - β1 * x_bar

这个关系式非常重要，它意味着无论 β1 是多少，最优的 β0 都必须满足这个等式。这也表明，拟合的回归线必然穿过数据的中心点 (x_bar, y_bar)。

现在，我们将这个 β0 的表达式代回原始的最小二乘准则中：

∑ [ yi - (y_bar - β1 * x_bar) - β1 * xi ]^2

经过整理，我们可以将其重写为：

∑ [ (yi - y_bar) - β1 * (xi - x_bar) ]^2

现在，我们定义中心化的变量：

  • y_tilde_i = yi - y_bar
  • x_tilde_i = xi - x_bar

于是，我们的最小化问题变成了：

∑ (y_tilde_i - β1 * x_tilde_i)^2

看！这正好是一个以中心化数据进行的过原点回归问题。y_tilde 是结果变量，x_tilde 是预测变量。

根据我们之前推导的过原点回归公式，最优的 β1_hat 为：

β1_hat = (∑ y_tilde_i * x_tilde_i) / (∑ x_tilde_i^2)

将中心化变量的定义代入：

β1_hat = [ ∑ (yi - y_bar)(xi - x_bar) ] / [ ∑ (xi - x_bar)^2 ]

得到最终公式

现在，我们对分子和分母同时除以 (n-1)：

  • 分子除以 (n-1) 是 y 和 x 的样本协方差 Cov(y, x)。
  • 分母除以 (n-1) 是 x 的样本方差 Var(x)。

因此：

β1_hat = Cov(y, x) / Var(x)

我们知道，相关系数 Cor(y, x) 的定义是 Cov(y, x) / (Sd(y) * Sd(x))，其中 Sd 表示标准差。同时，方差 Var(x) = Sd(x)^2。

所以，我们可以将公式进一步改写：

β1_hat = [Cor(y, x) * Sd(y) * Sd(x)] / [Sd(x)^2] = Cor(y, x) * [Sd(y) / Sd(x)]

至此，我们完成了证明。一般线性回归的斜率估计公式为：

β1_hat = Cor(y, x) * (Sd(y) / Sd(x))

而截距的估计公式由之前的关系式给出：

β0_hat = y_bar - β1_hat * x_bar

总结

本节课中我们一起学习了线性回归核心公式的完整推导。

  1. 我们从过原点回归这个简化模型入手，推导出其斜率估计公式 β_hat = (∑ yi * xi) / (∑ xi^2)。
  2. 然后，我们将一般线性回归问题，通过中心化数据，巧妙地转化为一个过原点回归问题。
  3. 最终，我们证明了线性回归的斜率估计值 β1_hat 确实等于 y 和 x 的相关系数乘以 y 的标准差与 x 的标准差之比，即 β1_hat = Cor(y, x) * (Sd(y) / Sd(x))。
  4. 截距 β0_hat 则由中心点的关系确定：β0_hat = y_bar - β1_hat * x_bar。

这个推导揭示了相关系数、标准差与回归斜率之间的深刻联系，是理解线性回归模型本质的重要一步。

009：回归至均值 📈

概述

在本节课中，我们将学习“回归至均值”这一重要概念。回归至均值是回归分析发展史上的一个里程碑，由弗朗西斯·高尔顿发现。我们将探讨其含义、常见例子，并学习如何使用相关系数来量化这一现象。

什么是回归至均值？

上一节我们介绍了回归模型的基础。本节中，我们来看看一个有趣且普遍的现象——回归至均值。

回归至均值提出类似这样的问题：为什么非常高的父母所生的孩子也往往很高，但通常没有父母那么高？同样，为什么非常矮的父母所生的孩子也往往很矮，但通常没有父母那么矮？

这种现象也出现在其他领域。例如，表现最佳的运动员在第二年往往表现稍差，而表现最差的运动员在第二年往往表现稍好。股票市场也常被如此讨论，一些表现最佳的股票之后往往会下跌。

我们将讨论其原因，并区分这究竟是内在规律，还是回归至均值效应，这是一个重要的问题。

这些现象都是所谓的“回归至均值”的例子。

弗朗西斯·高尔顿的贡献

正如上一张幻灯片提到的，回归至均值由弗朗西斯·高尔顿在其著名论文《遗传身高中向平庸的回归》中提出。

我喜欢通过考虑100%回归至均值的极端情况来理解这个概念。

想象我模拟生成两对独立的标准正态随机变量，它们之间没有线性关系，是独立的。

在我生成的第一组标准正态随机变量中，如果我取最大值，那么它在第二组中的配对值更小的概率会很高。

这简单地表明，给定一个非常大的X值，Y小于X的概率会随着X趋向于非常大的值而增大。

换句话说，给定一个非常小的X值，Y大于X的概率会随着X趋向于更小的值而增大。

然而，在大多数情况下，是内在成分和噪声的混合。例如，如果我给班上每个学生两次非常难的测验。

表现最好的学生很可能确实更了解材料。但是，测验是不完美的工具，测验作为工具存在很多误差，因此存在一些噪声。

所以，表现最好的学生可能确实比其他人更了解材料一点，但也可能受益于一些噪声，因此在第二次测验中他们的表现可能会稍差一些，可能仍然很好，但会稍差。对于表现最差的学生也是如此。

思考体育讨论中有多少内容实际上只是在谈论回归至均值，这常常很有趣。一个棒球运动员在某一年打击率惊人，下一年往往会稍低一些。这只是回归至均值吗？如果能量化它就好了。

这就是弗朗西斯·高尔顿在首次处理回归至均值时用回归所做的事情。

如何量化回归至均值？

我们想讨论弗朗西斯·高尔顿如何利用回归的思想，特别是相关系数（我们知道它与线性回归密切相关）来量化回归至均值。我们将主要通过一张图来展示，但在展示R代码之前，我想先设置好这张图。

在这个案例中，我的X将是孩子的身高，Y将是父母的身高。我将使用一个数据集，其中父母是单亲，即父亲。这些数据都已标准化，Y和X的均值为0，方差为1，我们现在应该对此非常熟悉。

回想一下，如果我这样做了，那么我的回归线将通过点(0,0)。并且，无论孩子的身高是结果变量还是父母的身高是结果变量，斜率现在都简单地等于相关系数。

我还想提一下创建图表时的一个小技巧：如果X是结果变量，而你恰好将结果变量绘制在横轴上，那么你绘制的直线的斜率需要是1除以相关系数。这只是因为我将结果变量绘制在了X轴上，所以提醒你一下，如果你对代码感到有点困惑。

R代码实现与解读

以下是创建图表的R代码步骤：

# 加载数据
library(UsingR)
data(father.son)

# 标准化数据：减去均值，除以标准差
y <- (father.son$sheight - mean(father.son$sheight)) / sd(father.son$sheight)
x <- (father.son$fheight - mean(father.son$fheight)) / sd(father.son$fheight)

# 计算相关系数
rho <- cor(x, y)

# 加载绘图库
library(ggplot2)

# 创建基础散点图
g <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y))
g <- g + geom_point(colour = "salmon", alpha = 0.6, size = 2)
g <- g + xlim(-4, 4) + ylim(-4, 4) + coord_equal()

# 添加对角线（恒等线）
g <- g + geom_abline(intercept = 0, slope = 1, colour = "black")

# 添加坐标轴
g <- g + geom_hline(yintercept = 0)
g <- g + geom_vline(xintercept = 0)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/fdbce39ad029b7f2e3a5ddc13d658510_5.png)

# 添加两条回归线
# 1. 以父亲身高预测儿子身高 (y ~ x)
g <- g + geom_abline(intercept = 0, slope = rho, colour = "blue", size = 1.5)
# 2. 以儿子身高预测父亲身高 (x ~ y)，注意斜率是 1/rho
g <- g + geom_abline(intercept = 0, slope = 1/rho, colour = "darkgreen", size = 1.5)

print(g)

运行代码后，我们得到相关系数 rho 约为 0.5。这意味着父亲身高和儿子身高之间的相关系数是 0.5。

现在让我们解读这张图，并用它来描述回归至均值。

如果观测值完美地落在一条线上，那么在这张图上，它必须是恒等线，因为我们已经标准化了X和Y。提醒一下，父亲身高被绘制为X变量，儿子身高被绘制为Y变量。

所以，如果我们有一个父亲身高为2（即高于均值2个标准差），并且没有噪声，我们会预测儿子身高也为2（即高于均值2个标准差）。

然而，当我们查看数据时，存在很多噪声。所以现在的预测不是2，而是在回归线上（蓝色线）。由于我们的设置方式，这简单地等于2乘以斜率，即相关系数。现在的预测点在这里（蓝色线与X=2的垂线的交点），介于2和0之间。

事实上，它正好等于2乘以相关系数。这就是回归至均值。这种乘法收缩就是弗朗西斯·高尔顿测量回归至均值的方式。

这个相关系数向水平线收缩的程度，就给出了回归至均值的程度。

如果我们讨论没有噪声的情况，这条线将完美地落在恒等线上。考虑全是噪声的情况，父亲身高对儿子身高没有信息，那么相关系数为0，这条线就会在水平轴上。那么，如果你必须根据父亲身高来预测儿子身高，你最终总会得到0的预测。

这个现象就是回归至均值。当然，我们也可以将其翻转过来考虑。

如果我们把儿子身高作为预测变量，父亲身高作为结果变量，你可以在这张图上可视化这一点，因为我用这条线（深绿色线）展示了这个图。回归至均值就是它向纵轴收缩的程度。

这是一个非常聪明的想法，我认为这是导致我们现在所认为的现代回归的基本思想之一。然而，回归至均值的概念在统计学中仍然占有重要地位，例如，当你研究纵向数据时，思考回归至均值的概念非常重要。

总结

本节课中，我们一起学习了“回归至均值”的核心概念。我们了解到：

  • 回归至均值描述了极端观测值在后续测量中倾向于更接近平均值的现象。
  • 它由弗朗西斯·高尔顿在研究遗传身高时发现并命名。
  • 这种现象普遍存在于许多领域，如教育测试、体育表现和金融市场。
  • 我们可以通过相关系数 ρ 来量化回归至均值的程度。预测值会向均值“收缩”，收缩因子就是相关系数。
  • 在标准化数据中，回归线的斜率等于相关系数，这为我们可视化理解回归至均值提供了直观工具。

理解回归至均值有助于我们区分真实效应与统计波动，是数据分析中一个至关重要的概念。

010：统计线性回归模型

在本节课中，我们将学习如何为线性回归建立统计模型，并理解如何从数据中推断总体参数。我们将从最简单的概率模型开始，逐步探讨其核心假设、参数估计及统计推断的基本原理。

从数学过程到统计推断

上一节我们介绍了使用最小二乘法寻找最佳回归线的数学过程。然而，我们希望进行统计推断，即基于样本数据对总体进行概括。这正是本节“统计线性回归模型”要讨论的内容。

线性回归的概率模型

我们从一个简单的线性回归概率模型开始。该模型将响应变量表述为解释变量的线性函数，并加入随机误差项。

模型公式

模型的基本形式如下：

Yi = β₀ + β₁Xi + εi

其中：

  • Yi 是第 i 个观测的响应变量。
  • Xi 是第 i 个观测的已知解释变量。
  • β₀ 和 β₁ 是我们希望估计的未知总体参数（截距和斜率）。
  • εi 是随机误差项，我们假设它服从独立同分布的正态分布，即 εi ~ N(0, σ²)。

关于误差项的说明

独立误差项可以理解为许多未被纳入模型但共同影响响应变量的因素的累积效应，其总体行为可以近似为独立同分布的高斯（正态）误差。我们暂时搁置其复杂的解释问题，先关注如何进行回归的统计推断。

模型的期望与方差

在给定解释变量 Xi 的条件下，响应变量 Yi 的期望值和方差具有明确的统计含义：

  • 条件期望：E[Yi | Xi] = β₀ + β₁Xi。这恰好是回归线在该点的值。
  • 条件方差：Var(Yi | Xi) = σ²。这表示数据点围绕回归线的波动程度。

需要强调的是，这里的方差 σ² 是条件方差，即给定 X 后 Y 的方差。由于通过 X 解释了一部分 Y 的变异，这个值通常小于响应变量 Y 自身的总体方差。此外，期望值和方差目前都是我们想要了解的总体参数，后续我们将通过样本数据来估计它们。

本节总结

本节课我们一起学习了统计线性回归模型的基础。我们建立了包含随机误差项的概率模型 Yi = β₀ + β₁Xi + εi，并明确了模型的核心假设：误差项 εi 独立且服从正态分布 N(0, σ²)。在此基础上，我们探讨了响应变量的条件期望与条件方差的含义，为后续的参数估计和统计推断奠定了理论基础。

011：解释回归系数 📊

在本节课中，我们将学习如何解释线性回归模型中的系数。我们将重点关注截距和斜率的统计含义，并探讨如何通过平移或缩放自变量来改变系数的解释方式，使其更具实际意义。

建立统计框架

上一节我们介绍了线性回归的正式统计框架。本节中，我们将基于这个框架来解释回归系数。

回归模型的基本形式为：
[
Y = \beta_0 + \beta_1 X + \epsilon
]
其中，(\beta_0) 是截距，(\beta_1) 是斜率，(\epsilon) 是误差项。

解释截距

截距 (\beta_0) 是当自变量 (X = 0) 时，响应变量 (Y) 的期望值。公式表示为：
[
E[Y | X = 0] = \beta_0
]

然而，在许多研究中，自变量为零的情况可能没有实际意义。例如，如果 (X) 代表血压，那么血压为零的人群的响应值通常不是研究重点。

为了解决这个问题，我们可以对自变量进行平移。考虑将 (X) 减去一个常数 (A)：
[
Y = \beta_0 + \beta_1 (X - A) + \epsilon
]
这等价于：
[
Y = (\beta_0 - \beta_1 A) + \beta_1 X + \epsilon
]

我们定义一个新的截距 (\tilde{\beta_0} = \beta_0 - \beta_1 A)。模型变为：
[
Y = \tilde{\beta_0} + \beta_1 X + \epsilon
]

平移自变量不会改变斜率 (\beta_1)，但会改变截距。新的截距 (\tilde{\beta_0}) 表示当 (X = A) 时 (Y) 的期望值。

一个常见的做法是令 (A) 等于 (X) 的样本均值 (\bar{X})。这样，截距就被解释为当自变量取平均值时，响应变量的期望值。

解释斜率

斜率 (\beta_1) 表示自变量 (X) 每增加一个单位，响应变量 (Y) 的期望变化量。

我们可以通过以下推导来验证：
[
E[Y | X = x+1] - E[Y | X = x] = [\beta_0 + \beta_1 (x+1)] - [\beta_0 + \beta_1 x] = \beta_1
]

因此，(\beta_1) 确实代表了 (X) 每变化一个单位所引起的 (Y) 的期望变化。

缩放自变量的影响

如果我们将自变量 (X) 乘以一个正常数 (a)，会发生什么？

原始模型：
[
Y = \beta_0 + \beta_1 X + \epsilon
]

将 (X) 替换为 (aX)：
[
Y = \beta_0 + \beta_1 (aX) + \epsilon = \beta_0 + (\beta_1 a) X + \epsilon
]

为了保持模型形式不变，我们定义一个新的斜率 (\tilde{\beta_1} = \beta_1 a)。但更常见的解释是：如果我们放大了自变量（乘以 (a)），那么为了得到正确的 (Y) 预测值，斜率必须被缩小（除以 (a)）。实际上，模型可以重写为：
[
Y = \beta_0 + \left( \frac{\beta_1}{a} \right) (aX) + \epsilon
]
这显示了乘以 (a) 的 (X) 对应着除以 (a) 的系数。

一个具体例子

假设我们建立了一个模型，其中：

  • (X) 是身高，单位是米（m）。
  • (Y) 是体重，单位是千克（kg）。
  • 估计出的斜率 (\beta_1) 的单位是 千克/米（kg/m）。

如果我们将身高的单位从米转换为厘米（cm），我们知道 1 米 = 100 厘米。这相当于将自变量 (X) 乘以 100。

为了使模型在新单位下给出相同的预测，斜率系数必须进行相应调整。新的斜率 (\beta_1^{new}) 应为：
[
\beta_1^{new} = \frac{\beta_1}{100}
]
其单位变为 千克/厘米（kg/cm）。

通过这个例子，我们可以清楚地看到缩放自变量如何影响斜率的数值和单位解释。

总结

本节课中，我们一起学习了如何解释线性回归模型的系数：

  1. 截距 (\beta_0) 通常解释为当所有自变量为零时响应变量的期望值。通过平移自变量（例如，减去其均值），可以使截距的解释更具实际意义。
  2. 斜率 (\beta_1) 解释了自变量每增加一个单位，响应变量期望值的变化量。这是回归分析的核心解释。
  3. 对自变量进行平移会改变截距但不改变斜率。
  4. 对自变量进行缩放（乘以一个常数）会以相反的比例改变斜率的数值和单位，因此在报告和解释系数时，必须清楚说明自变量的测量单位。

理解这些概念对于正确理解和传达回归模型的结果至关重要。

012：使用线性回归进行预测 📈

在本节课中，我们将学习如何使用线性回归模型进行预测。我们将通过一个具体的例子，解释回归系数，并演示如何计算和解释预测值。

预测Y值

我们可以通过公式 Ŷ = β̂₀ + β̂₁X 来预测Y值，其中X是我们想要预测的自变量值。


如果我们代入观测到的X值，就会得到拟合值。这些拟合值是我们试图使其尽可能接近观测数据的值。

但这并不意味着我们只能在拟合值处进行预测。我们可以将任何我们想要的X值代入方程进行预测。

然而，如果我们代入的X值位于用于构建模型的数据范围内，预测结果将更为合理。

后续课程中，我们还将讨论如何使用预测区间来处理这种不确定性。目前，我们先讨论如何获得预测值。在接下来的几张幻灯片中，我们将通过一个具体例子来解读截距、解读斜率，并从特定的回归模型设置中生成预测。

代码示例：解读回归系数

现在，让我们通过一些代码来解读回归系数，并实时展示回归系数的运行过程。数据集来自UsingR包中的钻石数据集，数据包括钻石价格（新加坡元）和钻石重量（克拉）。



要获取数据，首先需要使用library(UsingR)加载UsingR包，然后加载diamond数据集。

我们还将使用ggplot2进行绘图。

以下是创建绘图的命令。我们将绘图赋值给变量g。数据集是diamond，横轴变量是carat（克拉），纵轴变量是price（价格）。我们为坐标轴添加标签：X轴标签为“质量（克拉）”，Y轴标签为“价格（新加坡元）”。

运行这些代码后，我们添加数据点（黑色背景，浅色半透明）和回归线。使用geom_smooth(method = lm)可以轻松添加回归线，该方法默认使用Y作为结果变量，X作为预测变量。我们将回归线颜色设为黑色。

运行代码后调用绘图，即可看到图形。X轴是质量，Y轴是价格，图中显示的直线就是拟合线，即最小化各点到直线垂直距离平方和的直线。

获取拟合线

lm函数是R中用于线性模型的程序，回归是它的一个特例。

在公式中，波浪号~左侧是结果变量price，右侧是预测变量carat。默认情况下，lm包含截距项。如果不需要截距，必须在模型中明确指定。

逗号后面指定我们使用的数据集为diamond数据框。我们将模型结果赋值给变量fit。

运行后查看输出：

运行模型后，输入fit会打印出系数β̂₀和β̂₁。若想获得更详细的输出，可以使用summary(fit)。在课程后期，你将能解读此摘要中的所有数字。

目前，我们只关注系数。使用coef(fit)可以获取系数向量，其中包含截距（标记为(Intercept)）和carat变量的回归斜率。

观察斜率3721，并尝试解读它。

这意味着，钻石质量每增加1克拉，其预期价格将增加3,721新加坡元。

截距-259是0克拉钻石的预期价格。这个值本身意义不大，因为我们不关心0克拉的钻石。

中心化X变量

现在，我们对X变量进行中心化处理，使截距处于更易解释的尺度上。

首先，将新模型赋值给fit2，以免覆盖原始的fit。在lm公式中，结果变量保持不变，预测变量变为carat - mean(carat)。在lm的公式语句中进行算术运算时，必须用I()函数括起来。数据集仍然是diamond。

运行代码后查看新系数：

斜率保持不变，仍是3721，但截距变为500。这意味着，500新加坡元是该数据集中平均大小钻石（约0.2克拉）的预期价格。

改变单位

1克拉的增加量实际上很大。如果我们将单位改为0.1克拉呢？

我们可以简单地将系数除以10。因此，我们预计钻石质量每增加0.1克拉，价格将增加372新加坡元。

让我们在R中演示这一点。现在将新模型赋值给fit3。在公式中，我们放入carat * 10，因此这个新变量的单位是0.1克拉。数据集仍是diamond。

运行后获取系数：


现在斜率是372，而不是3,721。

进行预测

假设有人带来了三颗新钻石，重量分别为0.16克拉、0.27克拉和0.34克拉。他们想知道这些钻石的估计价格。

你可以手动计算，即取截距加上斜率乘以这些新值。计算出的预测价格分别为336、745和1006新加坡元。根据散点图，这个拟合的线性回归模型似乎相当不错。

通常，你不想进行这么多编码，尤其是当有很多回归变量时，你需要一个更通用的方法。R中有一个通用函数叫predict，它可以处理多种模型拟合的输出（线性模型是其中之一）。

我们使用predict函数，并传入新的数据框（包含carat变量的新值）进行预测。这会给出与手动计算相同的答案。这种方法在我们有多个回归变量和更复杂设置时具有良好的扩展性。

如果你省略newdata参数，直接使用predict(fit)，它会预测在观测到的X值处的Y值，即给出Ŷ值。如果你想在新的X值处预测，必须提供newdata参数。


预测过程图示

我想简要说明预测是如何完成的。图中蓝色点是观测数据点。红色点是我们使用predict命令得到的拟合值，即所有观测X值在拟合线上对应的点。如果我们从观测数据点向拟合线画垂直线，垂足就是这些红点。

当我们预测一个新的X值时，我们是在横轴上找到这个点（例如0.16、0.27、0.34），向上画一条线到拟合回归线，然后水平移动到价格轴，得到预测的美元金额。

总结

在本节课中，我们一起学习了如何使用线性回归模型进行预测。我们通过钻石数据的例子，详细解读了回归系数（截距和斜率）的含义，并演示了如何对X变量进行中心化和单位转换以改善解释性。最后，我们介绍了使用predict函数进行预测的通用方法，并图示了预测的几何过程。

013：残差 📊

在本节课中，我们将要学习回归分析中的一个核心概念：残差。我们将了解残差的定义、性质及其在诊断模型拟合优度中的重要作用。

课程概述

本节课以钻石数据集为例，探讨如何利用钻石重量解释其价格的变化。我们将重点关注在考虑了预测变量（重量）的线性影响后，模型中剩余的、未被解释的变异，这部分变异就称为残差。

上一节我们介绍了回归线的基本概念，本节中我们来看看如何量化数据点与回归线之间的差异。

残差的定义

让我们回顾一下正在考虑的线性回归模型。在我们的例子中，结果变量是价格，解释变量是重量。

模型公式如下：
Y_i = β_0 + β_1 * X_i + ε_i

其中：

  • Y_i 是第 i 个观测的结果（钻石价格）。
  • X_i 是第 i 个观测的预测变量值（钻石重量）。
  • ε_i 是误差项，我们通常假设它服从均值为0、方差为 σ² 的正态分布。

对应于 X_i 的回归线上的点，我们称之为拟合值 Ŷ_i：
Ŷ_i = β̂_0 + β̂_1 * X_i

残差 e_i 定义为观测值 Y_i 与模型拟合值 Ŷ_i 之间的垂直距离：
e_i = Y_i - Ŷ_i

最小二乘法的目标正是最小化所有残差的平方和：∑ e_i²。

残差的性质与用途

理解残差的性质有助于我们更好地解释和使用它们。

以下是残差的一些关键性质：

  • 均值为零：如果模型中包含了截距项，那么残差的样本和（以及样本均值）为零。
  • 与预测变量正交：在包含截距项的线性模型中，残差与模型中任何一个预测变量的乘积和也为零。
  • 误差项的估计：残差可以被视为模型中随机误差项 ε_i 的一种估计。但需注意，随意增加无关的预测变量会人为地减小残差，因此需谨慎解释。

残差主要有两大用途：

  1. 诊断模型拟合：残差图是诊断模型假设是否成立、识别模型拟合不佳区域（如非线性、异方差性、异常值）的强大工具。
  2. 创建调整后的变量：残差可以看作是移除了预测变量 X 的线性影响后的结果变量 Y。例如，如果我们想分析已根据重量调整后的钻石价格（即所有钻石价格处于同一重量尺度上），就可以使用以价格为结果、重量为预测变量的回归模型所得到的残差进行后续分析。

残差变异与系统变异

区分以下两种变异非常重要：

  • 系统变异：由回归模型解释的变异，即预测变量 X 所能说明的结果 Y 的变化部分。
  • 残差变异：在通过线性方式考虑了预测变量 X 的影响后，结果 Y 中剩余的、未被解释的变异。它体现了数据点围绕回归线的分散程度。

正如在钻石数据集的散点图中所示，如果不考虑重量，价格本身变异很大。当我们用重量进行线性解释后，价格围绕回归线的变异（即残差变异）就小得多。

总结

本节课中我们一起学习了回归分析的核心概念——残差。

我们明确了残差是观测值与模型预测值之间的差值（e_i = Y_i - Ŷ_i），理解了其均值为零、与预测变量正交等重要性质。更重要的是，我们认识到残差是诊断模型拟合优度的关键工具，并能用于创建移除了特定变量线性影响的新变量。

下一节，我们将利用这些知识，深入探讨如何通过绘制和分析残差图来评估和改进我们的回归模型。

014：残差计算与可视化示例

在本节课中，我们将学习如何计算和可视化回归模型的残差。我们将使用钻石数据集作为示例，通过R语言演示残差的计算方法，并探讨如何通过残差图评估模型拟合的优劣。

🔢 计算残差

首先，我们加载钻石数据集并拟合一个简单的线性回归模型，其中价格（price）是因变量，克拉数（carat）是自变量。

# 加载数据并拟合模型
data(diamonds)
y <- diamonds$price
x <- diamonds$carat
n <- length(y)
fit <- lm(y ~ x)

拟合模型后，我们可以通过多种方式获取残差。最直接的方法是使用resid()函数。

# 方法一：使用resid()函数
e <- resid(fit)

另一种方法是手动计算残差，即观测值减去预测值。

# 方法二：手动计算残差
y_hat <- predict(fit)
e_manual <- y - y_hat

我们可以验证这两种方法得到的残差在数值精度上是相同的。

# 验证残差一致性
max(abs(e - e_manual))  # 结果应接近0

此外，我们还可以验证残差的两个重要性质：残差之和为零，以及残差与自变量的乘积之和为零。

# 验证残差性质
sum(e)               # 应接近0
sum(e * x)           # 应接近0

📊 可视化残差

为了更直观地理解残差，我们可以绘制散点图并在图上添加回归线和残差线段。

# 绘制散点图与回归线
plot(x, y, main = "钻石价格 vs 克拉数", xlab = "克拉数", ylab = "价格")
abline(fit, lwd = 2, col = "blue")

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/8ece94babfd1c8012331681ad0a16442_26.png)

# 添加残差线段
for (i in 1:n) {
  segments(x[i], y[i], x[i], y_hat[i], col = "red", lwd = 1)
}

然而，这种散点图在展示残差变异时可能不够清晰。为了更好地评估残差，我们绘制残差与自变量的关系图。

# 绘制残差图
plot(x, e, main = "残差 vs 克拉数", xlab = "克拉数", ylab = "残差")
abline(h = 0, lty = 2, col = "gray")

在残差图中，我们希望看到残差随机分布在水平线（y=0）上下，无明显模式。如果存在明显模式，可能表明模型拟合不足或存在异方差性。

🔍 残差图的诊断作用

残差图是诊断模型拟合问题的重要工具。以下通过两个模拟示例说明残差图如何揭示模型中的非线性或异方差性。

示例一：非线性模式

假设真实关系为 y = x + sin(x) + ε，但我们用线性模型 y ~ x 进行拟合。

# 生成模拟数据
set.seed(123)
x <- runif(100, -3, 3)
y <- x + sin(x) + rnorm(100, sd = 0.2)

# 拟合线性模型
fit_linear <- lm(y ~ x)

# 绘制残差图
plot(x, resid(fit_linear), main = "残差图（非线性示例）", xlab = "x", ylab = "残差")
abline(h = 0, lty = 2, col = "gray")

在残差图中，正弦函数的模式变得非常明显，这在线性模型中未被捕捉。

示例二：异方差性

假设真实关系为 y = x + ε，其中误差的标准差随x增大而增大。

# 生成异方差数据
set.seed(123)
x <- runif(100, 0, 3)
y <- x + rnorm(100, sd = 0.1 * x)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/8ece94babfd1c8012331681ad0a16442_36.png)

# 拟合线性模型
fit_hetero <- lm(y ~ x)

# 绘制残差图
plot(x, resid(fit_hetero), main = "残差图（异方差示例）", xlab = "x", ylab = "残差")
abline(h = 0, lty = 2, col = "gray")

在残差图中，可以清晰看到残差的变异性随x增加而增大，这表明存在异方差性。

💎 钻石数据集的残差分析

回到钻石数据集，我们拟合包含截距和斜率的模型，并比较仅含截距的模型与完整模型的残差。

# 拟合仅截距模型
fit_intercept <- lm(price ~ 1, data = diamonds)
e_intercept <- resid(fit_intercept)

# 拟合完整模型
fit_full <- lm(price ~ carat, data = diamonds)
e_full <- resid(fit_full)

# 创建比较数据框
resid_data <- data.frame(
  residual = c(e_intercept, e_full),
  model = factor(rep(c("仅截距", "截距+斜率"), each = nrow(diamonds)))
)

# 绘制残差比较图
library(ggplot2)
ggplot(resid_data, aes(x = model, y = residual, fill = model)) +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  labs(title = "残差分布比较", x = "模型", y = "残差") +
  theme_minimal()

左图展示了钻石价格围绕平均价格的变异，右图展示了围绕回归线的残差变异。通过引入克拉数作为预测变量，我们解释了大部分价格变异，剩余的残差变异较小。

📝 总结

本节课中，我们一起学习了如何计算和可视化回归模型的残差。我们通过钻石数据集演示了残差的计算方法，并探讨了残差图在诊断模型拟合问题中的重要作用。残差图能够有效揭示非线性、异方差性等模型缺陷，是回归分析中不可或缺的工具。在下一节中，我们将深入讨论R平方，进一步量化模型解释的变异比例。

015：14_残差方差 📊

在本节课中，我们将要学习回归分析中的一个核心概念——残差方差。我们将探讨它如何计算，它与总方差及回归方差的关系，以及它如何帮助我们理解模型拟合的度量指标R平方。

残差方差的概念

在结束关于R平方的讨论之前，我们快速讲解一下残差方差。

残差方差是数据点围绕回归线的变异程度。残差是观测结果与拟合回归线之间的垂直距离。

如果模型中包含截距项，残差之和必须为零，这意味着它们的均值为零。因此，若要计算残差的方差，只需计算其平方的平均值。

公式：残差方差的估计值为 sum(residuals^2) / n，其中 sigma^2 表示围绕回归线的真实总体方差。

自由度与分母调整

然而，在计算中，大多数人使用 n-2 作为分母，而不是 n。因此，它并非简单的平均平方残差。

对于大样本量 n，使用 1/(n-2) 和 1/n 的差异可以忽略不计，但对于小样本量，这种差异可能很重要。

这样考虑的原因是：如果模型包含截距项，残差之和为零，这是一个约束条件。如果你知道 n-1 个残差，就能推算出第 n 个。此外，如果模型中包含一个协变量（预测变量），这会对残差施加第二个约束条件，因此我们失去了两个自由度。每增加一个回归变量，就会增加一个约束，失去一个自由度。

本质上，这就像说你并非真正拥有 n 个独立的残差，而是 n-2 个。因为如果你知道 n-2 个残差，就能推算出最后两个。这就是分母使用 n-2 的原因。

在R中提取残差方差

上一节我们介绍了残差方差的理论计算，本节中我们来看看如何在R语言中实际操作并提取这个值。

以下是如何从你的 lm 拟合对象中提取残差标准差（即残差方差的平方根），并将其赋值给一个变量，以便在R程序中使用。

# 定义数据
y <- diamond$price
x <- diamond$carat
# 拟合回归模型
fit <- lm(y ~ x)

如果直接运行 summary(fit)，输出结果会包含回归模型的摘要信息，如截距、斜率估计值等，其中也会显示残差标准差的估计值。

然而，如果你想将其作为一个可以赋值的对象来使用，只需使用 $sigma 符号。

# 提取残差标准差并赋值
sigma_hat <- fit$sigma

在这个特定的钻石数据示例中，计算出的值约为 31.84。

我们可以通过手动计算来验证这个公式。

# 手动验证计算
sqrt(sum(resid(fit)^2) / (length(y) - 2))

这段代码依次执行以下操作：resid(fit) 提取残差，平方后求和，再除以 n-2 得到平均平方残差，最后开方得到标准差 31.84。这证实了提取方法的正确性。

方差分解：总方差、回归方差与残差方差

现在，让我们回到分析钻石价格变异性的图表。

总变异性是数据点围绕其均值（中心）的平均平方偏差。为简化，我们暂时忽略分母，只讨论平方偏差的和。

我们可以将回归变异性定义为总变异性中被回归线解释掉的那部分。即，取回归线上点的预测值，计算这些预测值围绕响应变量均值的变异程度。这就是回归方差。

剩余未被解释的变异性，即数据点围绕回归线的变异，就是残差方差。

一个重要的恒等式是：总方差等于回归方差加上残差方差。

公式：总平方和 (SST) = 回归平方和 (SSR) + 残差平方和 (SSE)

R平方的定义与性质

因为残差方差和回归方差相加等于总方差，我们可以定义一个量来表示模型所解释的总方差百分比。

只需将回归方差除以总方差即可。这个量被称为R平方。

对于我们的钻石示例，R平方表示钻石价格变异中，能被与质量（克拉）的线性回归关系所解释的百分比。

关于R平方，需要记住以下几点事实：

  • R平方是响应变量的变异中，能被与预测变量的线性关系所解释的百分比。
  • R平方的值必须在0和1之间。因为回归平方和与残差平方和都是正数，且相加等于总平方和，这迫使R平方落在此区间内。
  • 如果我们定义 R 为预测变量与结果之间的样本相关系数，那么R平方在数值上就等于这个相关系数的平方，即 R^2。
  • R平方可能误导对模型拟合度的判断。例如，如果数据有一定噪声，且删除了中间区域的许多数据点，可能会得到更高的R平方值。或者，如果随意向线性模型中添加回归变量，R平方会增加，而均方误差（平均残差方差）会减小。因此，在使用这些指标评估模型拟合时，必须考虑到这些因素。

R平方的局限性：安斯科姆四重奏

安斯科姆构造了一个特别鲜明的例子，展示了四组数据集。它们拥有相同的R平方值、X和Y的均值与方差，以及完全相同的线性回归关系。

然而，当你观察它们的散点图时，会发现拟合线在每种情况下的含义截然不同。

以下是该示例中的四组数据：

  1. 第一组是典型的、略带噪声的线性关系，符合我们对良好线性回归的常规认知。
  2. 第二组数据明显存在曲线趋势，表明模型中缺少了某个项（如二次项）来捕捉这种曲率。
  3. 第三组数据中存在一个明显的离群值。
  4. 第四组数据中，几乎所有数据都堆积在一个位置，只有一个点远离主体。这就像在第一组数据的基础上删除了中间大部分点后得到的结果。

在所有这四种情况下，你都得到了相同的R平方值。但显然，这个单一的汇总数字抛弃了从简单散点图中可以获得的大量重要信息。

总结

本节课中我们一起学习了残差方差的核心概念。我们了解到它是衡量数据围绕回归线变异程度的指标，其计算需要考虑自由度调整（使用 n-2 作为分母）。我们探讨了总方差如何分解为回归方差（模型解释的部分）和残差方差（模型未解释的部分），并由此引出了R平方的定义——模型解释的方差百分比。

同时，我们通过R代码演示了如何从拟合模型中提取残差标准差，并重点强调了R平方作为模型拟合度指标的局限性，特别是通过安斯科姆四重奏的例子，明白了不能仅依赖R平方来评估模型，必须结合数据可视化等其他方法进行综合判断。

016：回归中的推断 📊

在本节课中，我们将要学习如何在回归模型的背景下进行统计推断。我们将回顾假设检验和置信区间的基本概念，并将其应用于回归模型的参数估计中。

上一节我们介绍了回归模型的基本形式及其参数估计方法。本节中我们来看看如何对这些参数进行统计推断。

模型回顾 📝

我们的模型假设响应变量 Y_i 由以下公式描述：
Y_i = β_0 + β_1 * X_i + ε_i
其中，误差项 ε_i 服从均值为0、方差为 σ^2 的正态分布。

我们暂时假设真实的模型是已知的，这将作为本课程大部分内容的基础。

统计推断基础 🔍

我假设你已经学习过一些统计推断知识，因此了解什么是置信区间和假设检验。

让我们回顾一下统计推断中的几个核心概念。

以下是统计量的一般形式及其性质：

  • 统计量 (θ_hat - θ) / SE(θ_hat) 通常具有以下性质。
  • 该统计量通常服从正态分布。
  • 在有限样本中，如果用估计方差代替真实方差，该统计量通常服从 t 分布。
  • 我们将看到，在我们所做的假设下，这在回归背景下也成立。

这个统计量可用于对参数 θ 进行假设检验，检验其是否等于某个原假设值，或者小于、大于、不等于该值。

我们也可以用以下形式创建置信区间：
估计值 ± 相关分位数 * 标准误
其中，Q 是相关分位数。例如，如果 α 为 5%，我们想要 95% 的置信区间，我们就取 97.5% 的分位数。这个分位数可能来自正态分布或 t 分布。

在回归案例中，在误差项独立同分布且正态的假设下，推断将遵循与你在统计推断课程中看到的非常相似的路线。

回归斜率的方差 📉

回归斜率的方差公式信息量很大。

回归斜率 β_1_hat 的方差公式为：
Var(β_1_hat) = σ^2 / ∑(X_i - X_bar)^2

这个方差涉及两个因素：

  • 数据点围绕真实回归线的变异程度，即 σ^2。
  • 预测变量 X 的变异程度，即 ∑(X_i - X_bar)^2。

我认为分子部分，即数据点围绕回归线的变异程度，比较容易理解为什么它越小，我们对回归斜率的估计就越好。

然而，为什么我们希望预测变量 X 有更大的变异，以获得更低的回归斜率方差，这可能不那么直观。让我们画图说明。

如果我们的预测变量 X 都紧密地聚集在一起，那么显然我们无法很好地估计一条直线，直线可以很容易地围绕那团点云弯曲，得到等效的拟合。

另一方面，如果我们把 X 值分散开，我们就能得到一条拟合得更好的回归线，斜率的方差也会更低。

事实证明，你能得到的最小方差是将一半的观测值推到一端，另一半推到另一端。但这样做，你很大程度上是在假设这两端之间存在一条直线，因为你没有收集任何数据来评估中间部分的情况。

截距的方差与 t 分布 📐

这里我也给出了截距的方差公式。它可能信息量稍小，因为截距通常不如斜率那么受关注。

截距 β_0_hat 的方差公式为：
Var(β_0_hat) = (1/n + X_bar^2 / ∑(X_i - X_bar)^2) * σ^2

尽管如此，我们可以计算并写下这个公式。

在这两种情况下，我们都用其逻辑估计量 σ_hat^2 代替 σ^2，即 1/(n-2) * ∑ e_i^2，我们在上一讲中介绍过。

我们的斜率或截距估计值减去真实值，再除以其标准误，服从自由度为 n-2 的 t 分布，这或许并不令人惊讶。

因此，这可以用来创建置信区间和进行假设检验。我们将在接下来的几张幻灯片中看一些例子。

但首先，我要向你展示所有这些结果都完全成立，并且是 R 语言报告的内容。

本节课中我们一起学习了如何在回归分析中进行统计推断。我们回顾了用于参数假设检验和构建置信区间的 t 统计量，并探讨了回归系数方差的构成，特别是预测变量的变异性如何影响估计的精度。

017：手动计算与R验证 📊

在本节课中，我们将通过一个具体的编码示例，手动计算简单线性回归的各项统计量，并与R语言lm()函数的结果进行对比验证。我们将使用diamond数据集，演示如何从原始数据一步步计算出回归系数、残差、标准误、t统计量及p值，并最终构建出与R输出完全一致的系数表。

数据准备与基础计算

首先，我们加载所需的数据集和库，并定义关键变量。

# 加载数据
library(UsingR)
data(diamond)

# 定义变量
y <- diamond$price
x <- diamond$carat
n <- length(y)

接下来，我们根据公式手动计算回归系数。斜率 beta1 的计算公式为：

公式： beta1 = cor(y, x) * (sd(y) / sd(x))

截距 beta0 的计算公式为：

公式： beta0 = mean(y) - beta1 * mean(x)

# 计算斜率 beta1 和截距 beta0
beta1 <- cor(y, x) * sd(y) / sd(x)
beta0 <- mean(y) - beta1 * mean(x)

计算残差与回归标准差

上一节我们计算了回归线。本节中，我们来看看如何评估模型拟合的精度。首先计算残差，即观测值与预测值之间的差异。

公式： residual = y - (beta0 + beta1 * x)

# 计算残差
residuals <- y - (beta0 + beta1 * x)

随后，我们估计回归线周围的标准差 sigma。这衡量了数据点围绕回归线的离散程度。注意，这里除以 n-2 以进行无偏估计。

公式： sigma = sqrt( sum(residuals^2) / (n-2) )

# 估计回归标准差 sigma
sigma <- sqrt(sum(residuals^2) / (n - 2))

计算标准误与检验统计量

为了进行统计推断，我们需要计算回归系数的标准误。以下是计算过程。

首先，计算自变量 x 的离差平方和（SSxx）。

公式： SSxx = sum( (x - mean(x))^2 )

# 计算 x 的离差平方和
SSxx <- sum((x - mean(x))^2)

接着，分别计算截距 beta0 和斜率 beta1 的标准误。

公式：

  • SE(beta0) = sigma * sqrt( (1/n) + (mean(x)^2 / SSxx) )
  • SE(beta1) = sigma / sqrt(SSxx)
# 计算标准误
se_beta0 <- sigma * sqrt((1 / n) + (mean(x)^2 / SSxx))
se_beta1 <- sigma / sqrt(SSxx)

有了估计值和标准误，我们就可以构建用于检验系数是否为零的 t 统计量。

公式： t_statistic = estimate / standard_error

# 计算 t 统计量
t_beta0 <- beta0 / se_beta0
t_beta1 <- beta1 / se_beta1

最后，根据 t 统计量计算双侧检验的 p 值。

# 计算 p 值（双侧检验）
p_beta0 <- 2 * pt(abs(t_beta0), df = n - 2, lower.tail = FALSE)
p_beta1 <- 2 * pt(abs(t_beta1), df = n - 2, lower.tail = FALSE)

手动构建系数表

以下是手动汇总所有计算结果，并将其整理为与R标准输出格式一致的系数表。

# 手动创建系数表
coef_table <- matrix(c(beta0, beta1, se_beta0, se_beta1, t_beta0, t_beta1, p_beta0, p_beta1),
                     nrow = 2, byrow = FALSE)
rownames(coef_table) <- c("(Intercept)", "x")
colnames(coef_table) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")

# 打印手动计算的系数表
print(coef_table)

使用R的lm()函数验证结果

现在，我们使用R内置的lm()函数来拟合相同的线性模型，以验证我们手动计算的结果。

# 使用 lm() 函数拟合模型
fit <- lm(y ~ x)

# 打印模型摘要中的系数表
summary(fit)$coefficients

通过对比，你会发现手动计算的coef_table与summary(fit)$coefficients的输出在数值上完全一致。这验证了我们所使用公式的正确性。

计算置信区间 🔍

在验证了基础计算后，我们可以进一步计算回归系数的置信区间。这能为我们提供估计值的不确定性范围。

首先，从拟合模型fit中提取所需信息。

# 从模型拟合结果中提取系数表
fit_summary <- summary(fit)
coef_summary <- fit_summary$coefficients

# 获取自由度
df <- fit$df.residual

对于截距和斜率，其95%置信区间的计算公式为：

公式： estimate ± t_quantile * standard_error

其中，t_quantile 是自由度为 df 的 t 分布在97.5%分位点的值。

# 计算截距的95%置信区间
t_quantile <- qt(0.975, df = df)
ci_intercept <- coef_summary[1, 1] + c(-1, 1) * t_quantile * coef_summary[1, 2]

# 计算斜率的95%置信区间（针对1克拉的变化）
ci_slope_per_carat <- coef_summary[2, 1] + c(-1, 1) * t_quantile * coef_summary[2, 2]

在实际解释时，我们可能更关心钻石重量每增加0.1克拉对价格的影响。只需将上述斜率的置信区间除以10即可。

# 计算斜率针对0.1克拉变化的置信区间
ci_slope_per_point1_carat <- ci_slope_per_carat / 10

根据结果，我们可以有95%的把握认为，钻石重量每增加0.1克拉，其价格预计将增加356至389新加坡元。

总结

本节课中，我们一起学习了简单线性回归模型的核心计算过程。我们从定义变量开始，手动计算了回归系数、残差、回归标准差、标准误、t统计量和p值，并成功构建了与R语言lm()函数输出完全一致的系数表。最后，我们还演示了如何计算并解释回归系数的置信区间。这个过程不仅验证了理论公式的正确性，也加深了我们对回归模型统计推断底层逻辑的理解。

018：预测 📊

在本节课中，我们将要学习如何使用回归模型进行预测。预测是数据科学的核心概念之一，而回归模型是实现预测的基础且强大的工具。我们将探讨如何利用回归模型进行点预测，以及如何评估预测的不确定性，例如通过计算预测区间和置信区间。

预测的核心概念

上一节我们介绍了回归模型的基础。本节中，我们来看看如何利用这些模型进行预测。

回归和广义线性模型是执行预测的核心技术。它们通常能产生非常好的预测结果，模型简洁且易于解释。此外，一个额外的好处是，你可以在不进行任何数据重采样的前提下，获得关于预测的推断。这里的“推断”指的是，我们可以得到预测值，并围绕这些预测值计算置信区间，以评估预测的不确定性。

在回归模型中实现这一点非常容易，在广义线性模型中相对容易，但在一些更高级的机器学习算法中则可能相当困难。

如何进行预测

我们可能想要预测响应变量，例如特定克拉数钻石的价格，或者特定父母身高下孩子的身高。在这两种情况下，最直接的估计是取我们想要预测的预测变量值 ( x_0 )，乘以相关的估计斜率 ( \hat{\beta}_1 )，然后加上截距 ( \hat{\beta}_0 )。

公式如下：
[
\hat{y}_0 = \hat{\beta}_0 + \hat{\beta}_1 x_0
]

然而，如果我们想成为优秀的统计学家，就需要评估该预测中的不确定性。因此，拥有一个预测区间会很有帮助。

这里有一个细微的区别：预测回归线在特定点的值，与预测该点未来的 ( y ) 值，是两个不同的概念。这反映在预测方差的公式中。

预测方差公式

预测方差首先与数据点围绕回归线的变异性 ( \hat{\sigma} ) 有关。预测误差的公式如下：

对于预测回归线在 ( x_0 ) 处的值（置信区间）：
[
\text{Var}(\hat{y}_0) = \hat{\sigma}^2 \left[ \frac{1}{n} + \frac{(x_0 - \bar{x})2}{\sum_{i=1} (x_i - \bar{x})^2} \right]
]

对于预测 ( x_0 ) 处一个新的 ( y ) 值（预测区间）：
[
\text{Var}(\hat{y}_{\text{new}}) = \hat{\sigma}^2 \left[ 1 + \frac{1}{n} + \frac{(x_0 - \bar{x})2}{\sum_{i=1} (x_i - \bar{x})^2} \right]
]

请注意，预测新 ( y ) 值的公式中多了一个“1”。这使得预测新值的区间比预测回归线值的区间更宽。

理解公式中的关键项

以下是公式中各项的含义：

  • ( \hat{\sigma}^2 ): 残差方差，衡量数据点围绕回归线的离散程度。
  • ( 1/n ): 样本量 ( n ) 越大，估计越精确，此项贡献的方差越小。
  • ( (x_0 - \bar{x})^2 / \sum (x_i - \bar{x})^2 ): 此项衡量预测点 ( x_0 ) 距离 ( x ) 均值 ( \bar{x} ) 的远近，以及 ( x ) 变量自身的变异性。
    • 当 ( x_0 = \bar{x} ) 时，此项为0，预测方差最小。
    • ( x ) 的变异性越大（分母越大），此项贡献的方差越小。

在R中实现预测与绘图

让我们通过钻石数据集的例子，在R中生成预测区间和置信区间。

以下是关键的R代码步骤：

# 拟合线性模型
fit <- lm(price ~ mass, data = diamond)

# 创建一个密集的x值网格用于预测
newx <- data.frame(mass = seq(min(diamond$mass), max(diamond$mass), length = 100))

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/4624749abdf4bf52773b2de9883a6c3a_6.png)

# 预测回归线值及其置信区间
pred_conf <- predict(fit, newdata = newx, interval = "confidence")

# 预测新y值及其预测区间
pred_pred <- predict(fit, newdata = newx, interval = "prediction")

# 将结果合并到数据框中以便绘图
plot_data <- data.frame(newx, conf_lwr = pred_conf[, "lwr"], conf_upr = pred_conf[, "upr"],
                        pred_lwr = pred_pred[, "lwr"], pred_upr = pred_pred[, "upr"])

# 使用ggplot2绘图
library(ggplot2)
ggplot(plot_data, aes(x = mass)) +
  geom_ribbon(aes(ymin = pred_lwr, ymax = pred_upr, fill = "预测区间"), alpha = 0.3) +
  geom_ribbon(aes(ymin = conf_lwr, ymax = conf_upr, fill = "置信区间"), alpha = 0.5) +
  geom_line(aes(y = pred_conf[, "fit"]), color = "red", size = 1) + # 拟合线
  geom_point(data = diamond, aes(x = mass, y = price), alpha = 0.6) + # 原始数据点
  scale_fill_manual(values = c("预测区间" = "blue", "置信区间" = "salmon")) +
  labs(x = "质量（克拉）", y = "价格", fill = "区间类型") +
  theme_minimal()

运行上述代码将生成一张图。图中：

  • 蓝色区域 代表预测区间（用于预测新观测值）。
  • 浅橙色区域 代表置信区间（用于预测回归线本身）。
  • 红色线 是拟合的回归线。
  • 灰色点 是原始观测数据。

可以观察到，预测区间总是比置信区间宽。此外，两个区间在 ( x ) 的均值 ( \bar{x} ) 附近最窄，随着向数据范围的边缘移动而变宽。

两种区间的直观理解

为什么会有这样的区别呢？我们可以这样思考：

  • 置信区间（预测回归线）：如果我们能在所有 ( x ) 值处收集无限多的数据，我们就能极其精确地知道回归线在哪里。因此，随着数据量 ( n ) 的增加，这个区间会不断变窄，最终收敛到回归线本身。这体现了统计估计的精度随样本量增加而提高。
  • 预测区间（预测新观测值）：即使我们知道了真实的回归线（即真实的 ( \beta_0 ) 和 ( \beta_1 )），由于模型中存在固有的随机误差 ( \epsilon )，新的 ( y ) 值仍然会有波动。这种不确定性不会因为收集更多数据而消失。公式中的“+1”就代表了这种不可消除的残差变异。

此外，区间在数据云中心更窄，在两端更宽，这告诉我们：在预测变量取值接近其均值的地方进行预测，我们的信心更足。如果我们试图预测远离观测数据范围的值（外推），区间会变得非常宽，这提醒我们外推预测具有很高的不确定性。

总结

本节课中我们一起学习了回归预测的核心内容。

我们了解到，使用回归模型进行预测不仅包括计算点估计 ( \hat{y}_0 )，更重要的是评估预测的不确定性。我们区分了两种区间：

  1. 置信区间：用于估计回归线在特定 ( x ) 值处的取值范围，其宽度随样本量增加而减小。
  2. 预测区间：用于预测在特定 ( x ) 值处一个新观测值 ( y ) 的取值范围，它包含了模型固有的残差变异，因此总是更宽。

我们还学习了预测方差的公式，理解了其组成部分的含义，并通过R代码实践了如何计算和可视化这些区间。记住，在数据范围中心进行预测最可靠，而外推预测需格外谨慎。

掌握了这些，你就为使用回归模型进行可靠的预测和推断打下了坚实基础。接下来，我们将进入更复杂的主题，探讨包含多个预测变量的回归模型。

019：极速入门knitr 📘

在本节课中，我们将学习如何使用knitr来创建R Markdown文档，这是完成课程作业所需的重要工具。knitr能够将R代码、文本和图表整合到一个文档中，并生成HTML、PDF等多种格式的输出。

创建R Markdown文档

上一节我们介绍了knitr的基本概念，本节中我们来看看如何创建一个R Markdown文档。

首先，在RStudio中点击“File”菜单，然后选择“New File”，接着选择“R Markdown”。


此时，RStudio会弹出一个新文档，其中包含一个简单的R Markdown模板。你可以编辑标题，例如将其改为“这是一个测试knitr文档”。

理解代码块

文档中包含一个R代码块的示例。以下是R Markdown中代码块的基本结构：

```{r}
# 你的R代码放在这里


在代码块的开头，你可以使用大括号 `{}` 来设置各种选项。例如，`{r}` 表示这是一个R代码块。在大括号内，你可以添加多个选项来控制代码块的执行和输出。

以下是几个常用的knitr选项：

*   **eval**：控制是否执行代码。例如，`eval=FALSE` 表示只显示代码而不执行。
*   **echo**：控制是否在输出中显示代码。例如，`echo=FALSE` 表示只显示结果，不显示代码。
*   **results**：控制结果的显示方式。例如，`results='hide'` 会隐藏代码的输出结果。
*   **fig**：控制图形的输出选项，如 `fig.width` 和 `fig.height` 可以设置图形的尺寸。

## 保存与编译文档

编辑好文档后，需要将其保存为 `.Rmd` 文件，例如命名为 `test.Rmd`。

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_8.png)
![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_9.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_6.png)

保存后，点击RStudio中的“Knit”按钮（通常显示为“Knit HTML”或一个编织图案的图标），knitr就会开始编译文档。

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_8.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_9.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_11.png)

编译过程会执行文档中的所有R代码块，并将结果与文本内容合并，最终生成一个独立的输出文件（如HTML文档）。

## 查看输出结果

编译完成后，knitr会生成一个输出文件（例如 `test.html`）。你可以在RStudio的预览窗格中直接查看，也可以使用R命令在浏览器中打开它。

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/16ac7db9b1d327d77d04defbb501f788_11.png)

例如，在R控制台中输入以下命令：
```r
browseURL("test.html")

这将在你的默认浏览器中打开生成的HTML文档。

本节课中我们一起学习了knitr的基本使用方法。我们了解了如何创建R Markdown文档，如何编写和设置R代码块，以及如何保存、编译并查看最终的输出文档。掌握这些技能对于完成本课程的数据分析作业至关重要。

020：多变量回归第一部分 🧮

在本节课中，我们将要学习多变量回归。我们将探讨当存在多个潜在预测变量，而非单一预测变量时，如何在线性回归框架下进行分析。多变量回归的核心在于，在考察某个预测变量与响应变量关系的同时，能够“控制”或“调整”其他变量的影响。

为何需要多变量回归？🤔

在任何使用预测变量 X 预测响应变量 y 的场景中，如果你发现了一个显著的关系，但该预测变量并非随机分配给受试者或研究单元，那么总存在一个担忧：可能存在另一个已知或未知的变量，能够解释你所发现的这种关系。

让我们来看一个例子。假设你的朋友分析了一些数据，其中包含人们的各种健康信息和饮食信息。他说：“我发现了一个有趣的现象：薄荷糖的使用量与用力呼气容积之间存在显著的回归关系。”用力呼气容积是衡量肺功能的一个指标。

你可能会对此表示怀疑，因为从生物学角度看，薄荷糖（主要成分是糖）与肺功能之间似乎没有合理的直接联系。你可能会想，是否有其他变量可以解释这种关系？这里有两种可能的假设：

  1. 你的朋友在众多变量中反复测试，恰好找到了一个显著相关的变量，这可能是多重比较带来的偶然关联。
  2. 更可能的情况是，吸烟者往往更常使用薄荷糖，而吸烟本身与肺功能下降有明确的长期关联。因此，薄荷糖与肺功能的关系可能是通过吸烟这个中介变量产生的间接效应，而非直接效应。

那么，如何验证在排除吸烟影响后，薄荷糖是否仍有独立效应呢？一个直观的方法是：分别分析吸烟者群体和非吸烟者群体，看看在每个群体内部，肺功能是否因薄荷糖使用量不同而有差异。这样，我们就在“控制”了吸烟状态的情况下进行比较。

多变量回归正是以线性方式自动化实现这一过程的方法。它通过特定的假设，尝试在考察一个预测变量（如薄荷糖使用量）时，将其他变量（如吸烟状态）的影响“保持恒定”。本节课我们将解释其工作原理、调整方式以及局限性。

多变量回归的另一个用途：预测 📈

除了上述用于解释变量间关系的用途，多变量回归也是强大的预测工具。

例如，在著名的“遗产健康奖”竞赛中，目标是根据参保人过去的理赔历史，预测其未来一年的住院天数。这对于医院、保险公司和医疗服务提供者至关重要。

竞赛提供了包含海量预测变量的历史理赔数据。简单的单变量线性回归在这里作用有限，它适合快速考察单个变量，但要获得良好的预测效果，需要更复杂的模型。

这涉及到两个关键部分：

  1. 模型搜索：如何选择纳入模型的预测变量。
  2. 避免过拟合：如果向模型中放入过多变量（甚至是随机向量），仅凭变量数量就可能使残差为零，但这并不意味着模型好。因此，纳入无关变量或遗漏重要变量都会带来后果。

在数据科学专项课程的《实用机器学习》课程中，你将深入学习与预测相关的模型选择策略。而在本课程中，我们将更侧重于前一部分提到的场景：构建简约的模型，并深入解读线性模型系数的含义。

值得一提的是，在那场预测竞赛中，精心构建的多变量线性模型已经能取得非常接近最优结果的预测性能。更复杂的机器学习方法（如随机森林、提升算法等）通常只能带来微小的额外提升。这充分说明了，对于许多预测问题，构思良好的线性模型是一个强大且有效的起点。

多变量线性模型公式 📝

我们的多变量线性模型与简单线性模型非常相似，只是包含了更多的预测变量。

模型的基本形式如下：

公式：
y = β₀ + β₁X₁ + β₂X₂ + ... + βₚXₚ + ε

其中：

  • y 是响应变量（如肺功能、未来住院天数）。
  • β₀ 是截距项。
  • β₁, β₂, ..., βₚ 是各个预测变量 X₁, X₂, ..., Xₚ 的系数（可理解为“斜率”）。
  • ε 是误差项。

例如：

  • 在肺功能的例子中，X₁ 可能是薄荷糖使用量（二元变量），X₂ 可能是吸烟包年数。
  • 在保险预测的例子中，X₁ 可能是上一年的理赔次数，X₂ 可能是是否患有特定心脏疾病。

请注意，这里的“线性”指的是模型对于系数 β 是线性的。这是一个关键点。即使我们将某个预测变量进行变换（例如，使用 X² 而不是 X），只要系数本身不被平方或进行其他非线性运算，模型就仍然是线性模型。

代码示例（概念性）：

# 假设 lm 是线性模型拟合函数
model <- lm(y ~ X1 + X2 + X3, data = my_data)
# 模型公式 y = β0 + β1*X1 + β2*X2 + β3*X3 即体现了线性于系数的特性

最小二乘估计 ⚖️

上一节我们介绍了多变量线性模型的形式，本节我们来看看如何估计模型中的系数。最小二乘法的思想可以直接从简单线性回归推广到多变量情形。

我们的目标是最小化残差平方和。残差是观测值 y_i 与模型预测值 ŷ_i 之间的差异。模型预测值是所有预测变量与其系数乘积之和。

因此，我们需要最小化的目标函数（损失函数）是：

公式：
∑ (y_i - (β₀ + β₁X_{i1} + β₂X_{i2} + ... + βₚX_{ip}))²

对所有观测点 i 求和。

最小二乘法就是要找到一组系数 (β₀, β₁, ..., βₚ)，使得上述残差平方和达到最小。每个观测点的残差都以同等权重进入这个损失函数。

总结 🎯

本节课我们一起学习了多变量回归的基础知识。

我们首先探讨了其核心用途之一：在分析一个预测变量与结果的关系时，通过纳入其他变量来“控制”或“调整”混杂因素的影响，从而更准确地估计目标变量的独立效应，并以薄荷糖、吸烟和肺功能的关系为例进行了说明。

接着，我们了解了多变量回归的另一个重要用途——预测，它通过整合多个预测变量的信息来构建有效的预测模型。

然后，我们给出了多变量线性模型的一般公式 y = β₀ + β₁X₁ + ... + βₚXₚ，并强调了“线性”是指模型对于系数 β 是线性的这一关键点。

最后，我们介绍了用于估计模型系数的最小二乘法原理，其目标是最小化所有观测点的残差平方和。

在接下来的课程中，我们将继续深入多变量回归的细节，包括系数的解释、模型假设以及注意事项。

021：多变量回归第2部分 📊

在本节课中，我们将要学习多变量回归模型中回归系数的估计原理。我们将从一个简单的“过原点回归”模型出发，逐步推导出当模型包含两个或多个预测变量时，如何估计其系数。理解这一推导过程，有助于我们深入理解多变量回归“调整”其他变量影响的本质。

从过原点回归出发

上一节我们介绍了回归模型的基本概念，本节中我们来看看如何具体得到回归系数的估计值，或者说R等软件是如何计算出这些结果的。我们将进行一个推导，如果你对此不感兴趣，可以跳过这部分内容。

首先，回忆一下“过原点回归”模型。该模型强制回归线穿过原点，没有截距项。它只有一个预测变量X和一个结果变量Y。其斜率估计公式如下：

公式：
斜率估计值 = Σ(Xi * Yi) / Σ(Xi²)

这个公式我们在课程早期讨论过。

推导两个预测变量的回归估计

现在，让我们尝试利用上面的结果，来推导当模型包含两个预测变量时的最小二乘估计。如果你能理解这个过程，就能将其推广到三个、四个甚至更多预测变量的情况。

我们的目标是推导，或者说构建出估计量。我将展示一个比线性代数方法更直观的推导过程。如果你对线性代数推导感兴趣，我也有一些相关的视频资源。

考虑两个预测变量的情况。我们的最小二乘准则是最小化以下平方和：

公式：
Σ[ Yi - (X1i * β1) - (X2i * β2) ]²

例如，X1i可能代表截距项，X2i可能是我们感兴趣的协变量，比如血压。

以下是推导的关键技巧：

  1. 假设我们固定了β1的值。

  2. 定义一个新的变量：Ỹi = Yi - X1i * β1。

  3. 此时，最小二乘准则变为：Σ[ Ỹi - X2i * β2 ]²。

  4. 这正好是一个以X2i为唯一预测变量、Ỹi为结果变量的“过原点回归”问题。因此，我们知道β2的估计值必须满足：

    公式：
    β2估计值 = Σ(Ỹi * X2i) / Σ(X2i²)

  5. 注意，Ỹi是β1的函数。所以，我们可以将这个β2的表达式代回原始的最小二乘准则中。

  6. 这样，我们就得到了一个只包含β1的方程，并且它同样是一个关于β1的“过原点回归”问题。

最终推导出的有趣结果是：

  • 预测变量X1的回归系数β1，其估计值恰好等于：先将X2从Y和X1中线性剔除（即分别对Y和X1做关于X2的回归，取残差），然后对这两个残差做“过原点回归”所得到的系数。
  • 同理，预测变量X2的回归系数β2，其估计值等于：先将X1从Y和X2中线性剔除，然后对这两个残差做“过原点回归”所得到的系数。

这个观点非常重要，我们将反复强调。多变量回归中，一个变量的系数，是在剔除了所有其他变量对该变量以及结果变量的线性影响之后所得到的净效应。这解释了多变量回归在做什么：它通过调整其他变量的影响，来估计某个变量的独立贡献。

核心结论：系数的“调整”含义

为了确保你理解这个重要结论，我们在此重申一遍。

假设我们有两个协变量X1、X2和一个结果变量Y。我们拟合的模型包含这两个协变量及其系数β1和β2。

核心结论是：

  • 拟合得到的系数β1，等同于我们先将第二个变量X2从结果Y和第一个预测变量X1中线性剔除，然后对得到的残差进行“过原点回归”所获得的系数。
  • 这意味着，β1是针对X2的影响进行调整后的X1的效应，因为X2对Y和X1的线性影响已被移除。
  • 同样，系数β2也可以被类似地解释为针对X1调整后的效应。

这就是为什么多变量回归关系通常被认为是“调整了所有其他变量”后的关系。在下一部分，我们将讨论如何将此概念推广到包含多个预测变量的完整多变量回归中。

应用于简单线性回归的例子

让我们通过一个已知的例子来验证这个结论：简单线性回归，其中一个协变量是截距项，另一个是感兴趣的协变量。

假设模型为：Yi = β1 * X1i + β2 * X2i。这里我们假设X2是截距项（通常为常数1）。

根据我们的规则，β1的估计值是通过将X2从Y和X1中剔除后得到的。

  1. 将X2（截距）从Y中剔除：用Y对截距做回归，拟合值就是Y的均值Ȳ，残差为Yi - Ȳ，即Y的中心化形式。

  2. 将X2（截距）从X1中剔除：用X1对截距做回归，拟合值就是X1的均值X̄1，残差为X1i - X̄1，即X1的中心化形式。

  3. 然后，对这两个中心化后的变量进行“过原点回归”。β1的估计公式为：

    公式：
    β1估计值 = Σ[ (Xi - X̄) * (Yi - Ȳ) ] / Σ (Xi - X̄)²

这正是标准的简单线性回归斜率公式，也等于相关系数乘以Y与X的标准差之比。这验证了我们推导的正确性。

推广到更多变量

现在，让我们将这一思路扩展到两个以上的变量，从而深入理解多变量回归。

我们的最小二乘准则是最小化观测结果Yi与所有预测变量线性组合之间距离的平方和：

公式：
Σ[ Yi - (β1*X1i + β2*X2i + ... + βp*Xpi) ]²

  • 对于两个预测变量（一个截距，一个斜率），这相当于最小化数据点到回归线的垂直距离平方和。
  • 对于三个预测变量（一个截距，两个回归变量），拟合的将是一个平面，准则是最小化数据点到该平面的垂直距离平方和。
  • 当变量多于三个时，我们无法可视化，但思想是通用的：最小化Y到某个广义“超平面”的垂直距离。

我们所看到的核心思想是：多变量回归通过取残差的方式来“调整”其他变量。
例如，对于系数β1，我们的解释是：其余回归变量X2到Xp的线性影响已同时从Y和X1中被移除，然后我们研究这两个残差之间的关系。因此，在某种意义上，Y和X1都针对其余变量进行了调整。β1可以被视为在调整了与其他变量的线性关系后的“净效应”或“调整后效应”。

总结

本节课中我们一起学习了多变量回归系数估计的核心原理。我们从过原点回归出发，推导了两个预测变量情况下的估计过程，并得出了关键结论：多变量回归中每个变量的系数，都代表了在线性剔除所有其他变量对该变量和结果变量的影响后，该变量与结果之间的净关系。这个“调整”的思想是理解和解释多变量回归模型的基础。接下来，我们将通过计算机实验来验证这一原理，并在后续课程中通过具体例子学习如何解释多变量回归模型中的系数。

022：多元回归模型（续）📈

在本节课中，我们将继续探讨多元回归模型。我们将通过一个模拟演示来理解回归系数的估计过程，并详细解释多元回归系数的含义。此外，我们还将回顾线性模型的基本组成部分及其重要性。

模拟演示：回归系数的估计过程 🔍

上一节我们介绍了多元回归的基本概念，本节中我们来看看如何通过模拟来理解回归系数的估计。

我们进行一个模拟实验。生成100个观测值，并创建三个预测变量 x1、x2 和 x3，它们都服从标准正态分布。

响应变量 y 的生成公式如下：
y = 1 + x1 + x2 + x3 + ε
其中，所有理论上的总体回归系数均为1，ε 是随机误差项。

接下来，我们获取 y 在剔除 x2 和 x3 影响后的残差。lm 函数默认包含截距项。resid 函数用于提取残差。

同样地，我们对 x1 进行相同操作，即剔除 x2、x3 和截距项的影响。

然后，我们使用这两组残差进行过原点回归估计，得到的系数是 0.11053。

我们也可以用 lm 函数来验证这一点，将残差化的 y 对残差化的 x1 进行回归（需移除截距项），并提取系数。

该系数应与上述过原点回归的估计值完全相同。

需要指出的是，这个系数与将 y 对 x1、x2、x3 和截距项进行多元回归后得到的 x1 的系数是相同的。

这证明了，线性模型在估计某个预测变量的系数时，本质上是剔除了其他变量对该预测变量和响应变量的影响后，再进行估计。

多元回归系数的解释 📝

现在，我们来讲解如何解释多元回归模型的系数。下一讲我们将结合具体实例进行深入分析。

给定一组协变量取值 x1 到 xp，我们的回归预测值是各个 x 与其系数乘积之和，即线性模型公式。

现在，考虑第一个回归系数 β1 对应的预测变量 x1 增加一个单位，即取值从 x1 变为 x1 + 1。

那么，新的预测值公式中，x1 项变为 (x1 + 1)。

如果我们计算 x1 增加一个单位后的期望响应值与原始期望响应值之差，结果正好等于 β1。

请注意，在整个讨论中，其他变量 x2 到 xp 都保持固定不变。x1 增加一个单位（无论其单位是什么）。

因此，多元回归系数的解释是：在保持其他所有预测变量不变的情况下，该预测变量每增加一个单位，响应变量的期望变化量。这是解释多元回归关系的核心，即“保持其他变量恒定”或“调整了其他变量的影响”。

下一讲我们将通过多个具体例子来加深理解，但这就是其背后的数学原理。

线性模型的基本组成部分 📊

接下来，我们快速回顾线性模型的基本组成部分，它们与简单线性回归中的性质完全相同。

我们的线性模型是：
Y = β0 + β1*x1 + ... + βp*xp + ε
其中，我们假设误差项 ε 服从正态分布。

由此，我们得到一系列性质。这些性质与简单线性回归基本相同，只是现在我们引入了更多系数。

以下是线性模型的关键组成部分列表：

  • 拟合响应值：在任何特定预测变量值下的拟合响应值 ŷ，只需将拟合系数乘以对应的预测变量值并求和即可。这适用于新数据或观测数据。
  • 残差：观测值 y 与拟合值 ŷ 之差。我们的目标是最小化残差平方和。
  • 方差估计：残差的平均平方和。在简单线性回归中，我们除以 n-2；在多元回归中，我们除以 n-p。这是一个技术细节，因为知道了 n-p 个残差，由于线性约束，就能知道剩下的 p 个。你可以大致将其理解为残差的平均平方。
  • 新观测值的预测：对于一组新的预测变量值 x1 到 xp，将其代入线性模型，乘以各自系数后相加，就得到了预测的响应值。
  • 系数的标准误与检验：每个系数估计 β̂j 都有其标准误 σ(β̂j)。我们可以使用 t 检验 来检验某个特定系数是否为零，因为 (β̂j - βj) / σ(β̂j) 服从 t 分布。简单线性回归中的所有推断方法都延续到多元回归。
  • 预测区间：与简单线性回归一样，预测的响应值也有标准误，我们可以计算预测区间。

希望这些内容对你来说并不陌生。它们与简单线性回归几乎完全相同，只是现在我们有了更多的项。简单线性回归有截距和一个协变量两项，现在我们只是增加了更多的协变量。

线性模型对数据科学的重要性 ⚙️

在进入具体例子之前，我们想讨论一下线性模型对数据科学家的重要性。

在进行任何机器学习或复杂算法之前，线性模型应该是你的首选尝试。

以下是线性模型的主要优势：

  • 简洁且易于理解：它们以简洁且被深入理解的方式，描述了预测变量与响应变量之间的关系。
  • 强大的基础：尽管一些现代机器学习算法可能在线性假设等方面有所超越，但线性模型始终应该是你分析的起点。
  • 功能广泛：线性模型的能力可能超乎你的想象。例如，你可以将时间序列（如音频）分解为其谐波分量（离散傅里叶变换可视为线性模型的拟合）。你可以使用线性模型灵活地拟合相当复杂的函数和曲线。你可以将因子变量作为预测变量，此时方差分析（ANOVA）和协方差分析（ANCOVA）都是线性模型的特殊情况。你可以揭示响应变量与多个预测变量之间的复杂关系，并建立相当准确的预测模型。

总结 ✨

本节课中，我们一起学习了多元回归模型的核心内容。

我们通过一个模拟实验，演示了回归系数估计的本质是剔除其他变量的影响。我们明确了多元回归系数的解释：在保持其他变量不变的情况下，某个预测变量的单位变化所引起的响应变量期望变化。我们还系统回顾了线性模型的基本组成部分及其性质，这些是进行统计推断和预测的基础。最后，我们强调了线性模型在数据科学中的基础性和重要性，它应是分析工作的首要工具。

在接下来的课程中，我们将通过具体的实例来应用和巩固这些概念。

023：多变量回归示例（第一部分）🎯

在本节课中，我们将学习多变量回归分析。我们将通过一个著名的数据集——瑞士生育率数据，来探讨如何构建和解释包含多个预测变量的回归模型。我们将重点关注模型系数的解释、变量选择的影响，以及如何通过模拟理解“辛普森悖论”等现象。

数据介绍与探索 📊

我们使用的数据集来自一个著名的统计学教科书，即瑞士生育率数据。该数据集位于R的datasets包中。

require(datasets)
data(swiss)
?swiss

该数据集包含了1888年瑞士47个法语省份的标准化生育率指标和一些社会经济状况指标。变量包括：

  • Fertility：标准化生育率度量。
  • Agriculture：从事农业的男性百分比。
  • Examination：省份中在军事考试中获得特定分数的应征者百分比。
  • Education：接受过初等以上教育的人口百分比。
  • Catholic：天主教徒（相对于新教徒）的百分比。
  • Infant.Mortality：活产婴儿中未满一岁即夭折的百分比。

我们的目标是探究哪些因素可以解释各省的生育率差异。

在深入分析之前，我们先对数据进行基本的散点图探索。

library(GGally)
ggpairs(swiss)

生成的矩阵图中，对角线显示了每个变量的分布，下三角区域显示了任意两个变量之间的散点图及其平滑曲线，上三角区域则显示了对应的相关系数。例如，Fertility和Agriculture的散点图显示其相关系数为0.35。

构建多变量回归模型 🔧

上一节我们观察了数据的基本关系，本节中我们来看看如何构建一个包含所有预测变量的回归模型。

我们首先拟合一个模型，将Fertility作为因变量，所有其他变量作为自变量。

summary(lm(Fertility ~ ., data = swiss))

如果我们只关注系数表，可以运行：

fit <- lm(Fertility ~ ., data = swiss)
summary(fit)$coefficients

以Agriculture的系数-0.1721为例，其解释为：在保持Examination、Education、Catholic和Infant.Mortality不变的情况下，我们预计从事农业的男性比例每增加1%，标准化生育率将减少0.1721。

需要注意的是，Agriculture是以百分比（而非比例）度量的，这会影响系数的尺度。标准误0.0703衡量了该系数估计的精确度。t统计量-2.448是估计值除以标准误的结果，对应的p值0.018表明，在5%的显著性水平下，该系数是统计显著的。

对比单变量与多变量模型 🔄

上一节我们拟合了包含所有变量的模型，本节中我们来看看如果只使用Agriculture一个预测变量，结果会如何变化。

summary(lm(Fertility ~ Agriculture, data = swiss))$coefficients

有趣的现象出现了：在单变量模型中，Agriculture的系数约为0.1942，是正数。而在包含所有变量的多变量模型中，其系数变成了-0.1721，是负数。调整了其他变量（如教育水平）后，Agriculture对Fertility的影响方向发生了逆转。

这种现象常被称为“辛普森悖论”（或更准确地说是“辛普森表象”）。它并不真正矛盾，而是强调了在观察性研究中，忽略重要变量（混杂因素）可能导致对核心变量关系的错误解读。两个模型中，Agriculture的系数都是统计显著的，但科学解释却截然不同。

通过模拟理解系数逆转 🎲

为了更直观地理解系数为何会逆转，我们通过一个模拟示例来演示。

假设我们模拟100个数据点：

  • x2：代表时间，取值从1到100。
  • x1：代表财富，它随时间x2线性增长，并加上随机波动。
  • y：代表幸福感，其真实生成模型为：y = -1 * x1 + 1 * x2 + 随机噪声。即，幸福感与财富x1负相关，与时间x2正相关。
n <- 100
x2 <- 1:n
x1 <- .01 * x2 + rnorm(n, sd=.5)
y <- -x1 + x2 + rnorm(n, sd=.1)

首先，我们错误地只使用x1来预测y：

summary(lm(y ~ x1))$coefficients

结果会显示一个很大的正系数，这与真实的负关系完全相反。这是因为x1中包含了与y正相关的x2的信息。

然后，我们拟合正确的模型，同时包含x1和x2：

summary(lm(y ~ x1 + x2))$coefficients

此时，我们得到了接近真实值-1和+1的系数估计。

以下是理解这一过程的步骤：

  1. 绘制y与x1的散点图，并按x2着色。可以看到明显的正相关趋势，但颜色梯度揭示了x2的混杂作用。
  2. 计算残差：分别从y和x1中剔除x2的线性影响。
  3. 绘制y的残差与x1的残差的散点图。此时，x2的影响已被移除，图中清晰地显示出x1与y之间真实的负相关关系。

线性模型在多变量情况下的工作方式，正是通过从预测变量和响应变量中剔除其他变量的影响，来揭示“纯净”的关系。

回到瑞士数据与模型选择 🧐

经过模拟，我们对多变量回归有了更深的理解。现在回到瑞士数据集。

Agriculture的效应在加入其他变量（特别是Education和Examination）后发生逆转。值得注意的是，Education与Agriculture呈负相关（-0.64），而Education对Fertility有更强的负效应。此外，Education和Examination高度相关（0.7），它们测量的是类似的概念。

因此，如果有人仅通过Agriculture对Fertility的单变量回归就声称两者存在正向因果关系，这个结论是非常可疑的。即使在关联性层面，这个正向关联也很容易被其他合理变量的引入所打破甚至逆转。

最后，我们看看在模型中引入一个完全多余的变量会发生什么。这里“多余”指的是该变量是已有变量的线性组合。

swiss$z <- swiss$Agriculture + swiss$Education
summary(lm(Fertility ~ . + z, data = swiss))$coefficients

你会发现，新变量z的系数显示为NA。这是因为R无法估计一个完全由其他变量线性构成的变量的独立效应。如果你在拟合线性模型后看到NA系数，最可能的原因就是你引入了一个与其他变量存在完全或近似线性关系的变量。

总结 📝

本节课中我们一起学习了多变量回归分析的核心实践：

  1. 模型构建与解释：我们学习了如何拟合多变量线性回归模型，并正确解释在控制其他变量条件下某个预测变量的系数。
  2. 辛普森悖论：我们通过瑞士数据实例和模拟实验，深入理解了忽略重要混杂变量可能导致对核心变量关系的方向性误判。
  3. 变量选择的重要性：科学分析是一个动态过程，需要基于领域知识谨慎地选择纳入模型的变量，以得到可靠且可解释的结论。
  4. 冗余变量：我们了解到，在模型中引入与其他变量存在线性关系的冗余变量，会导致模型无法识别，系数估计为NA。

记住，回归分析不仅仅是运行代码，更需要对数据生成过程、变量间的潜在关系进行深入的思考和判断。

024：多元回归示例（第二部分）🎯

在本节课中，我们将学习线性回归模型如何处理因子变量。我们将通过一个具体的例子，演示如何将分类变量（如处理组/对照组、政治派别）作为回归预测因子，并解释其系数的含义。理解这一点对于正确解读模型结果至关重要。

线性回归模型的灵活性

你可能会惊讶地发现线性回归模型非常灵活。例如，你可以将因子变量作为回归预测因子，并实现诸如方差分析（ANOVA）这样的方法。如果你之前听说过方差分析，它实际上是线性模型的一个特例。

让我们通过一个例子来具体说明。假设我们有一个协变量 x，其取值为0或1。我们来看看当把这个变量放入线性回归模型时会发生什么。

模型公式：
y = β₀ + x * β₁ + ε

在这个模型中，y 是结果变量，β₀ 是截距，β₁ 是 x 的系数，ε 是误差项。这里，x 只取0或1。例如，x=0 可能代表对照组，x=1 代表处理组。

  • 对于处理组（x=1），其均值为 β₀ + β₁。
  • 对于对照组（x=0），其均值为 β₀。

拟合这个模型后，正如你所料：

  • 处理组的估计均值就是处理组样本的均值，即 β̂₀ + β̂₁。
  • 对照组的估计均值就是对照组样本的均值，即 β̂₀。

因此，β₁ 被解释为 x=1 的组（即处理组）相对于 x=0 的组（即对照组）的均值响应增加量（如果为负则是减少量）。这是一种将二水平因子变量作为线性回归变量的巧妙方法。它不仅通过拟合值告诉你两组的均值，还自动提供了比较两组的推断（例如T检验）。顺便提一下，β₁ 的T检验与假设方差齐性的两独立样本T检验完全相同。

扩展到多水平因子

我们可以将此方法扩展到两个以上的水平。例如，假设你有一个三水平的变量，比如你想比较某个结果与美国政治派别（民主党、共和党、独立人士）的关系。

你可以通过创建以下变量来实现：

  • 变量 x₁：如果是共和党则为1，否则为0。
  • 变量 x₂：如果是民主党则为1，否则为0。

我们暂时省略代表独立人士的变量 x₃（设为1代表独立人士，否则为0），因为它将是冗余的。

模型公式：
y = β₀ + x₁ * β₁ + x₂ * β₂ + ε

  • 如果一个人是共和党（x₁=1， x₂=0），其均值为 β₀ + β₁。
  • 如果一个人是民主党（x₁=0， x₂=1），其均值为 β₀ + β₂。
  • 如果一个人是独立人士（x₁=0， x₂=0），其均值仅为 β₀。

这就是我们不能包含第三项的原因。在我们这样设置变量的情况下，如果你不是共和党也不是民主党，那么你一定是独立人士。因此，加入第三个变量将是冗余的，不会提供任何新信息。我们有三个均值（共和党、民主党、独立人士）和三个参数（β₀， β₁， β₂）。如果添加一个额外参数，模型就会出问题。

系数解释与参照水平

如果我们观察这里的均值：

  • β₁ = 共和党均值 - 独立人士均值。所以 β₁ 比较了共和党与独立人士。
  • β₂ = 民主党均值 - 独立人士均值。所以 β₂ 比较了民主党与独立人士。
  • β₁ - β₂ 则比较了民主党与共和党。

通过省略代表独立人士的回归变量，截距 β₀ 就变成了独立人士的均值。所有其他系数都被解释为相对于独立人士的变化。β₁（共和党协变量前的系数）被解释为共和党与独立人士之间的差异，β₂（民主党协变量前的系数）被解释为民主党与独立人士之间的差异。

这都是省略独立人士这一回归项的结果。如果我们包含了独立人士的回归项而省略了共和党的，那么截距将代表共和党的均值，民主党前的系数将是民主党 vs 共和党，独立人士前的系数将是独立人士 vs 共和党。

R语言在默认情况下会自动为你处理这一点。如果你将一个变量作为因子纳入模型，它会自动选择其中一个水平作为参照水平。关键在于，在处理因子变量和线性模型时，你设定哪个水平作为参照水平会产生重大影响。这些系数的解释会因设置方式的不同而有很大差异。

R语言实例演示

现在，让我们通过一个R语言示例来观察如何处理因子变量，并看看R是如何对待它的。

首先，我们使用 InsectSprays 数据集，其中结果变量 count 是昆虫数量，预测因子 spray 是杀虫剂类型（因子变量，水平为A到F）。

我们先绘制数据以直观了解：

library(ggplot2)
ggplot(InsectSprays， aes(x = spray， y = count， fill = spray)) +
  geom_violin() +
  labs(title = "不同杀虫剂的昆虫数量分布"， x = "杀虫剂类型"， y = "昆虫数量")

从图中可以大致看出不同喷雾剂处理后的昆虫数量差异。

接下来，我们拟合一个线性模型，将 count 作为结果，spray 作为预测因子。

# 拟合模型
fit <- lm(count ~ spray， data = InsectSprays)
# 查看系数表
summary(fit)$coefficients

输出结果中，你会看到截距（Intercept），以及 sprayB， sprayC， sprayD， sprayE， sprayF 的系数，但 sprayA conspicuously 缺失了。这意味着所有系数都是相对于 sprayA 进行比较的。

  • 截距（14.5） 被解释为 sprayA 组的均值。
  • sprayB 的系数（0.833） 被解释为 sprayB 组均值与 sprayA 组均值的差异。
  • sprayC 的系数（-12.4） 被解释为 sprayC 组均值与 sprayA 组均值的差异。

如果你想得到 sprayB 的均值，就是截距加上 sprayB 的系数：14.5 + 0.833 = 15.333。

R在幕后自动创建了虚拟变量。我们可以手动创建这些变量来验证：

# 手动创建虚拟变量，以A为参照水平
fit_manual <- lm(count ~ I(spray == 'B') + I(spray == 'C') +
                   I(spray == 'D') + I(spray == 'E') + I(spray == 'F')，
                 data = InsectSprays)
summary(fit_manual)$coefficients

你会发现结果与让R自动处理因子变量时完全相同。

包含所有水平与省略截距

如果你在模型中包含了代表所有水平的变量（包括参照水平），会发生什么？

# 错误示例：包含所有水平（包括A）和截距
fit_redundant <- lm(count ~ I(spray == 'A') + I(spray == 'B') +
                      I(spray == 'C') + I(spray == 'D') +
                      I(spray == 'E') + I(spray == 'F')，
                    data = InsectSprays)
summary(fit_redundant)$coefficients

你会看到 sprayA 对应的系数显示为 NA。这是因为模型试图用7个参数（1个截距 + 6个水平）去拟合6个均值，造成了冗余，R会自动丢弃其中一个。

如果你希望每个系数直接代表各组的均值（而不是与参照组的差值），你可以通过省略截距来实现：

# 无截距模型，系数代表各组的均值
fit_no_intercept <- lm(count ~ spray - 1， data = InsectSprays)
summary(fit_no_intercept)$coefficients

现在，输出中包含了 sprayA 到 sprayF 所有水平的系数，它们分别等于各自组的样本均值。这个模型与包含截距的模型在数学上是等价的，信息量相同，唯一的区别在于系数的解释。

  • 有截距模型：截距是参照组（A）的均值，其他系数是各组与参照组的差异。P值检验的是这些差异是否显著不为零（即该组是否与A组不同）。
  • 无截距模型：每个系数是各组的均值。P值检验的是该组均值是否显著不为零。这通常不是我们关心的主要问题，我们更关心组间比较。

更改参照水平

默认情况下，R选择因子中字母顺序最低的水平作为参照水平（这里是A）。你可以使用 relevel 函数轻松更改参照水平。

# 将参照水平改为C
InsectSprays$spray2 <- relevel(InsectSprays$spray， ref = "C")
fit_new_ref <- lm(count ~ spray2， data = InsectSprays)
summary(fit_new_ref)$coefficients

现在，输出中截距代表 sprayC 的均值，sprayA 的系数代表A组与C组的差异，sprayB 的系数代表B组与C组的差异，依此类推。sprayC 本身不再出现在系数表中，因为它现在是参照水平。

重要总结与注意事项

让我们总结一下关于线性模型中因子变量的关键点：

  1. 默认行为：在R中，当你在 lm() 中包含一个因子变量时，R会自动包含一个截距，并将因子的第一个水平（按字母或因子顺序）设为参照水平。
  2. 系数解释：
    • 截距被解释为参照水平的均值。
    • 其他因子水平的系数被解释为该水平与参照水平均值的差异。
    • 该水平的实际均值 = 截距 + 该水平的系数。
  3. 假设检验：
    • 截距的检验是检验参照水平的均值是否为零。
    • 其他系数的检验是检验该组与参照组的均值差异是否显著。
  4. 无截距模型：如果省略截距，每个系数代表对应组的均值，但此时的假设检验（检验均值是否为零）通常不具实际意义。
  5. 更改参照水平：使用 relevel() 函数可以方便地更改参照水平，从而改变所有系数的比较基准。

关于本示例分析的注意事项：
本例中使用线性回归（假设正态分布）对计数数据建模可能并不理想，因为计数数据有下限（0），且方差可能不恒定（异方差性）。更合适的模型可能是泊松回归或负二项回归等广义线性模型（GLM），这将在课程后续部分涉及。目前，我们的均值估计可能是正确的，但基于正态假设的推断（如P值、置信区间）可能不准确。处理异方差性等问题需要更高级的统计方法。

本节课总结

在本节课中，我们一起深入探讨了在线性回归模型中如何处理和解释因子变量。我们学习了：

  • 如何将二水平及多水平分类变量纳入回归模型。
  • R语言如何通过创建虚拟变量并设置参照水平来自动处理因子变量。
  • 模型中有截距和无截距时，系数含义的根本区别。
  • 如何正确解释相对于参照水平的系数及其假设检验。
  • 使用 relevel() 函数更改参照水平以进行不同的组间比较。

理解这些概念对于避免错误解读回归结果至关重要，尤其是在涉及分类预测因子的研究中。

025：多变量回归示例第三部分 🧮

在本节课中，我们将学习多变量回归的一个重要示例，它构成了“协方差分析”（Ancova）主题的基础。我们将通过一个具体案例，演示如何拟合具有不同截距和斜率的多个回归线。

上一节我们介绍了多变量回归的基本概念，本节中我们来看看如何在实际数据中应用这些概念，特别是处理分类变量时的情况。

数据准备与探索

我们将使用之前见过的瑞士数据集。我们的目标是建立模型，将生育率（Fertility）作为结果变量，农业人口比例（Agriculture）作为预测变量。

首先，我们关注Catholic变量。该变量表示各省天主教徒的比例。通过绘制直方图可以发现，其分布呈现明显的双峰形态。

这是因为在数据所属的历史时期，大多数省份要么以天主教徒为主，要么以新教徒为主。

接下来，我们使用dplyr包创建一个二元的Catholic变量（CatholicBin）。如果一个省份的天主教徒比例超过50%，则该变量为1（多数为天主教），否则为0（多数为新教）。

# 创建二元分类变量
Swiss <- Swiss %>%
  mutate(CatholicBin = ifelse(Catholic > 50, 1, 0))

执行head(Swiss)命令后，可以看到数据集中新增了CatholicBin变量。


现在绘制数据图。图中展示了CatholicBin因子变量（0代表新教省份，1代表天主教省份）与农业人口比例的关系。图中存在一些潜在的异常观测值，但为了专注于当前的教学目标——拟合分组回归线，我们暂时忽略这些问题。后续课程将专门讨论异常值和影响点。


我们的目标是分别为天主教省份和新教省份拟合两条独立的回归线。

模型构建与解释

以下是我们可以考虑的三种模型。设：

  • Y 为生育率（Fertility）
  • X₁ 为农业人口比例（Agriculture）
  • X₂ 为二元宗教变量（CatholicBin，1代表天主教，0代表新教）

模型1：忽略分组（单一直线）

第一个模型完全忽略省份的宗教差异，对所有数据拟合一条直线。

公式：
E[Y | X₁, X₂] = β₀ + β₁X₁

这个模型假设宗教分组对生育率与农业人口的关系没有影响。

模型2：包含分组变量，无交互项（平行直线）

第二个模型在预测变量中加入了宗教分组X₂，但没有交互项。

公式：
E[Y | X₁, X₂] = β₀ + β₁X₁ + β₂X₂

现在，让我们分析这个公式在不同X₂取值下的含义：

  • 当 X₂ = 0（新教省份）：公式简化为 β₀ + β₁X₁
  • 当 X₂ = 1（天主教省份）：公式变为 (β₀ + β₂) + β₁X₁

因此，拟合这个模型相当于拟合了两条斜率相同（都是β₁）但截距不同的直线。新教省份的截距是β₀，天主教省份的截距是β₀ + β₂。

模型3：包含分组变量和交互项（不同斜率和截距的直线）

第三个模型不仅包含分组变量，还加入了X₁和X₂的交互项（X₁ * X₂）。

公式：
E[Y | X₁, X₂] = β₀ + β₁X₁ + β₂X₂ + β₃(X₁ * X₂)

同样，我们分析不同X₂取值下的情况：

  • 当 X₂ = 0（新教省份）：公式简化为 β₀ + β₁X₁
  • 当 X₂ = 1（天主教省份）：公式变为 (β₀ + β₂) + (β₁ + β₃)X₁

现在，我们拟合了两条斜率和截距都不同的直线。

  • 新教省份的直线是：截距 = β₀， 斜率 = β₁
  • 天主教省份的直线是：截距 = β₀ + β₂， 斜率 = β₁ + β₃

系数解释：

  • β₂ 代表了从新教省份切换到天主教省份时，截距的变化量。
  • β₃ 代表了从新教省份切换到天主教省份时，斜率的变化量。

代码实现

理论可能有些抽象，让我们通过实际代码来加深理解。

# 拟合模型1：忽略分组
fit1 <- lm(Fertility ~ Agriculture, data = Swiss)

# 拟合模型2：包含分组，无交互项（平行线）
fit2 <- lm(Fertility ~ Agriculture + CatholicBin, data = Swiss)

# 拟合模型3：包含分组和交互项（不同斜率的线）
fit3 <- lm(Fertility ~ Agriculture * CatholicBin, data = Swiss)
# 等价于 Fertility ~ Agriculture + CatholicBin + Agriculture:CatholicBin

# 查看模型3的摘要
summary(fit3)

通过运行summary(fit3)，我们可以查看β₀、β₁、β₂、β₃的具体估计值、标准误和显著性水平，从而判断宗教分组是否显著影响了生育率与农业人口关系的斜率和截距。

总结

本节课中我们一起学习了多变量回归中处理分组变量的核心方法。我们通过瑞士数据集示例，逐步构建了三个模型：

  1. 忽略分组变量的简单线性模型。
  2. 加入分组变量但无交互项的模型，该模型为不同组拟合了平行线（相同斜率，不同截距）。
  3. 加入分组变量和交互项的模型，该模型为不同组拟合了完全独立的线（不同斜率和不同截距）。

理解交互项（β₃(X₁ * X₂)）的含义至关重要，它量化了预测变量X₁（如农业人口比例）对结果Y（如生育率）的影响如何随着分组变量X₂（如宗教）的不同而发生变化。这是进行协方差分析（ANCOVA）和探索组间差异的基础。

026：多变量回归示例（第四部分）📊

在本节课中，我们将通过一个实际案例，学习如何在R语言中拟合包含分类变量的多元线性回归模型。我们将逐步构建模型，从最简单的回归线开始，逐步引入分类变量及其与连续变量的交互作用，并可视化不同模型的结果。

模型拟合与可视化

上一节我们讨论了如何将分类变量引入回归模型。本节中，我们来看看如何在R中实际操作，并绘制出相应的拟合线。

首先，我们拟合一个完全不考虑宗教因素的简单模型。该模型仅使用农业人口比例来预测生育率。

以下是拟合该模型的R代码：

# 拟合不考虑宗教的简单线性模型
fit <- lm(fertility ~ agriculture, data = swiss)
# 创建基础散点图并保存
g1 <- ggplot(swiss, aes(x = agriculture, y = fertility)) + geom_point()
# 获取截距和斜率
intercept <- coef(fit)[1]
slope <- coef(fit)[2]
# 将拟合线添加到图中
g1 <- g1 + geom_abline(intercept = intercept, slope = slope, color = "blue")
print(g1)
# 查看系数摘要
summary(fit)$coefficients

运行后，我们得到截距约为60，农业变量的斜率约为0.19。这条拟合线忽略了图中数据点的颜色差异。

拟合平行线模型

接下来，我们引入一个分类变量。我们创建一个二元变量CatholicBin，表示一个省份是否以天主教徒为主。然后，我们拟合一个模型，该模型为不同宗教的省份拟合具有不同截距但相同斜率（即平行线）的回归线。

以下是实现步骤：

# 创建二元分类变量
swiss$CatholicBin <- as.factor(1 * (swiss$Catholic > 50))
# 拟合包含分类变量的模型（平行线）
fit_parallel <- lm(fertility ~ agriculture + CatholicBin, data = swiss)
# 查看系数
summary(fit_parallel)$coefficients

在R中，将变量编码为0和1时，R会自动将其视为因子（分类变量）。然而，良好的习惯是显式使用as.factor()或factor()函数。这是因为，如果一个变量取值0、1、2（例如代表三种发色），若不声明为因子，R会将其视为连续变量，错误地认为“2是1的两倍”。

系数表显示，我们得到了一个截距、一个农业变量的斜率，以及一个代表天主教省份的截距变化量（约为7.88）。现在，我们可以将两条平行线添加到图中：

# 为新教徒省份绘制线（截距为coef[1]，斜率为coef[2]）
g1 <- g1 + geom_abline(intercept = coef(fit_parallel)[1],
                       slope = coef(fit_parallel)[2],
                       color = "darkred")
# 为天主教省份绘制线（截距为coef[1]+coef[3]，斜率仍为coef[2]）
g1 <- g1 + geom_abline(intercept = coef(fit_parallel)[1] + coef(fit_parallel)[3],
                       slope = coef(fit_parallel)[2],
                       color = "darkblue")
print(g1)

图中显示了两条斜率相同但截距不同的拟合线。

拟合带交互项的模型（非平行线）

最后，我们允许斜率也随宗教类别变化。在R公式中，使用星号*可以自动包含交互项以及所有主效应。

以下是具体操作：

# 拟合包含交互项的模型（允许斜率和截距都变化）
fit_interaction <- lm(fertility ~ agriculture * CatholicBin, data = swiss)
# 查看系数
summary(fit_interaction)$coefficients

系数表现在包含四项：

  1. 截距（新教徒省份）
  2. 农业斜率（新教徒省份）
  3. 天主教主效应（截距变化量）
  4. 农业与天主教的交互项（斜率变化量）

因此，对于天主教省份，其回归线的：

  • 截距 = coef[1] + coef[3]
  • 斜率 = coef[2] + coef[4]

现在，我们将这两条非平行线添加到图中：

# 为新教徒省份绘制线
g1 <- g1 + geom_abline(intercept = coef(fit_interaction)[1],
                       slope = coef(fit_interaction)[2],
                       color = "salmon", size = 1.2)
# 为天主教省份绘制线
g1 <- g1 + geom_abline(intercept = coef(fit_interaction)[1] + coef(fit_interaction)[3],
                       slope = coef(fit_interaction)[2] + coef(fit_interaction)[4],
                       color = "navy", size = 1.2)
print(g1)

从图中可以清晰地看到，现在两条线不再平行。根据系数（天主教相关项均为正数）可以判断，截距较高的蓝色线对应天主教省份，截距较低的鲑鱼色线对应新教省份。

总结与后续

本节课中，我们一起学习了多元线性回归中处理分类变量的完整流程：

  1. 拟合了基础模型。
  2. 引入了分类变量，拟合了平行线模型。
  3. 引入了交互项，拟合了非平行线模型，并进行了可视化比较。

通过观察最终的拟合图，我们可能会对少数远离拟合线的数据点（例如图中的两个蓝色点）产生疑问：它们是否对模型产生了过度影响？这引出了回归诊断中的重要话题。在下一讲中，我们将探讨残差分析和影响力诊断，学习如何评估数据点对模型拟合的影响，并判断模型假设是否合理。

027：调整变量的影响 🧮

在本节课中，我们将学习多元回归分析。我们将通过一系列示例，探讨在模型中调整一个变量后，另一个变量与结果之间观测到的关系会发生何种变化。理解这一点对于正确解读回归模型的结果至关重要。

概述 📋

多元回归模型允许我们同时考虑多个预测变量对结果的影响。然而，当我们引入新的变量进行调整时，原有变量系数的估计值和显著性可能会发生显著变化，甚至出现方向性的逆转。本节将通过可视化的模拟数据，展示几种典型的变化模式。

示例一：调整前后效应基本一致

首先，我们来看一个简单直接的情况。在这个模拟数据中，我们有一个二元分组变量（例如处理组和对照组）和一个连续变量X。

以下是用于拟合的模型公式：
Y = β₀ + β₁ * Group + β₂ * X + ε

这个模型会拟合出两条平行的直线。其中，β₁ 代表两组之间截距的差异，即处理效应；β₂ 代表X与Y关系的共同斜率。

在上图中，水平线展示了忽略X时，分组变量（红色与蓝色）的边际效应。当我们拟合上述包含X的模型后，得到的调整后处理效应（即β₁的估计值）与边际效应大致相同。此外，对于X的任意取值，我们都有足够多的红点和蓝点来进行直接比较，这种数据分布类似于随机化试验的结果。

示例二：调整后效应消失（辛普森悖论）

上一节我们看到了调整变量后效应基本不变的情况。本节中我们来看看一个更戏剧性的场景：调整一个变量后，原本显著的效应完全消失。

在这个例子中，如果忽略X（看水平线），红色处理组和蓝色对照组之间存在明显的边际差异。然而，当我们拟合包含X的模型并观察截距的变化时，处理效应变得微乎其微。

这种情况常发生在非随机化的观察性研究中。例如，X可能是胆固醇水平，而“处理”是是否服用降压药。胆固醇水平很可能决定了患者是否会开始服用该药。此时，调整X实质上是在调整导致接受处理的混杂因素。这个例子也凸显了观察性数据分析的挑战。

关于这个数据集，还有一个关键点：

  • 知识：在这个数据集中，X的值几乎可以完美预测个体属于哪个组（例如，X值小的是蓝组，X值大的是红组）。
  • 缺乏直接证据：不存在一个X的取值能让我们同时看到红组和蓝组的数据点，因此我们无法在控制X的条件下直接比较两组。我们严重依赖模型的假设来进行比较。

示例三：调整后效应方向逆转

我们已经看到了效应消失的例子。现在，让我们探讨一种更令人惊讶的现象：调整变量后，效应的方向发生了逆转。

在这个场景中，忽略X时，红色组的平均结果高于蓝色组（边际效应显著）。但是，当我们拟合模型并调整X后，估计出的处理效应显示蓝色组的结果反而高于红色组。调整后的估计值与未调整的估计值符号相反。

这种现象在统计学中常被称为辛普森悖论。它指的是，一个变量与结果的关系在考虑（或忽略）另一个变量时发生了方向性的逆转。从图中的数据分布可以直观理解其原因，这其实并非真正的“悖论”。

示例四：从无效应变为有效应

前面的例子展示了显著效应因调整而减弱或反转。相反的情况也可能发生：调整变量可能揭示出隐藏的效应。

在此例中，如果忽略X，红组和蓝组之间几乎没有边际差异。然而，当我们把X纳入模型进行调整后，却发现了显著的处理效应。

这再次说明，调整变量对结果的影响没有简单的规则。效应的显著性可能从有到无、从无到有、保持显著、保持不显著，或者符号翻转，所有排列组合都可能出现。

示例五：存在交互作用的情况

到目前为止，我们假设了处理组和对照组具有相同的斜率。现在，我们来看一个斜率不同的情况，这涉及到交互作用。

在这个数据集中，红线和蓝线明显不平行。为了拟合这样的数据，我们需要使用包含交互项的模型：
Y = β₀ + β₁ * Group + β₂ * X + β₃ * (Group * X) + ε

这个模型允许两组拥有不同的截距和不同的斜率。

这里有一个非常重要的概念：当模型中存在交互项时，主效应（β₁）本身不能单独被解释为处理效应。从图中可以清楚地看到，处理效应的大小和方向依赖于X的取值。在X较小时，蓝组结果更高；在X较大时，红组结果更高。因此，不存在一个统一的“处理效应”。这个例子展示了，如果忽略了变量间的交互作用（效应修饰），而错误地拟合了平行线模型，可能会得到完全误导性的结论。

连续变量的情况：三维视角与残差图

我们之前的讨论集中在二元处理变量和连续变量X上。但同样的原理完全适用于多个连续预测变量的情况。

假设我们有两个连续变量X1和X2，以及结果Y。在二维散点图上（Y vs X1），我们可能看不到明显的关系。

但是，如果我们能在三维空间中观察（Y, X1, X2），或者更实际地，使用残差图，关系就会显现。

以下是分析方法：

  1. 首先，分别用X2对Y和X1进行回归，得到残差。
    • Resid_Y = Y - (预测值 from Y ~ X2)
    • Resid_X1 = X1 - (预测值 from X1 ~ X2)
  2. 然后，绘制 Resid_Y 对 Resid_X1 的散点图。

这张残差图展示了在移除X2的线性影响后，Y与X1之间剩余的关系。图中清晰显示了一个强烈的正相关，这是在原始二维图中被X2掩盖了的效应。

这个例子的核心要点是：调整变量导致系数发生各种变化（反转、增强、减弱、显著性改变）的现象，并不仅限于二元变量的情况，在多个连续预测变量的回归中同样普遍存在。

如何选择正确的模型？ 🤔

在目睹了调整变量可能带来的巨大变化后，一个自然的问题是：究竟哪个模型才是“正确”的？

答案是：没有纯粹的统计答案。确定正确的模型需要将统计分析与专业领域知识紧密结合。

以下是构建模型时的关键考虑因素：

  • 科学逻辑：调整的变量在专业上是否有意义？例如，在血压模型中同时放入收缩压和舒张压，虽然它们高度相关，但让其中一个的效应消失并不具有科学洞察力。
  • 避免过度调整：如果调整的变量本身就是处理决策的中介或结果，可能会不恰当地“调整掉”真实的处理效应。
  • 团队协作：最好的做法是组建一个跨学科团队，包括领域专家、统计学家和数据分析师，共同讨论：
    • 不同调整策略下系数发生重大变化的原因。
    • 每种模型的利弊及其在科学上的合理性。
  • 自动化模型选择：机器学习课程中会提到的自动化算法（如逐步回归、LASSO）更侧重于预测精度，而非保证系数的可解释性。对于旨在理解变量作用的解释性建模，亲手进行细致的、基于知识的模型构建是不可替代的。

总结 📝

本节课我们一起学习了多元回归中调整变量的核心影响。我们通过多个模拟示例看到，在回归模型中引入新的协变量后，我们关心的变量系数可能会：

  1. 基本保持不变。
  2. 减弱甚至消失。
  3. 方向发生逆转（辛普森悖论）。
  4. 从隐藏变为显著。
  5. 在存在交互作用时，其含义发生根本变化。

这些变化在二元处理和连续预测变量，以及多个连续预测变量的情形中都会发生。最终，构建一个可信的、可解释的模型，依赖于统计工具与专业领域知识的深度融合与反复迭代。希望本课能帮助你未来在看到系数因调整而改变时，理解其背后可能的数据结构。

028：残差与诊断（第一部分）📊

在本节课中，我们将学习回归分析中的残差、诊断与验证。我们将从线性模型的基本概念出发，探讨如何利用残差评估模型拟合度，并介绍杠杆点和影响点的概念及其诊断方法。

线性模型回顾

上一节我们介绍了线性模型的基本形式。让我们先回顾一下线性模型的设定。

线性模型通常表示为：

[
Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_p X_{ip} + \epsilon_i
]

其中，( Y_i ) 是因变量，( X_{ij} ) 是自变量，( \beta_j ) 是模型系数，( \epsilon_i ) 是误差项。

在本讲座中，我们通常假设误差项 ( \epsilon_i ) 服从正态分布。请注意，部分结论依赖此假设，部分则不依赖。

残差的定义与估计

本节中，我们来看看残差是如何定义的，以及如何估计残差方差。

残差 ( e_i ) 定义为观测值 ( Y_i ) 与模型预测值 ( \hat{Y}_i ) 的差值：

[
e_i = Y_i - \hat{Y}_i
]

残差方差 ( \hat{\sigma}^2 ) 的估计量是残差平方和的平均值。为了获得无偏估计，我们使用 ( n - p ) 作为分母，而不是 ( n )：

[
\hat{\sigma}^2 = \frac{\sum_{i=1}^n e_i^2}{n - p}
]

其中，( n ) 是样本量，( p ) 是模型参数的数量（包括截距项）。

在R中生成诊断图

在R语言中，可以非常方便地获取残差图以及一系列衡量模型拟合优度与缺陷的图形。

以下是具体操作步骤：

  1. 加载数据集 swiss。
  2. 拟合一个包含所有自变量的线性模型。一个简便的方法是在公式的波浪符号 ~ 后使用点号 .，这表示纳入数据框中除因变量外的所有变量。
  3. 使用 par(mfrow = c(2, 2)) 命令将图形排列成2行2列，以便在一页中显示所有图形。
  4. 对拟合的模型对象直接使用 plot() 命令，R会自动生成一系列残差诊断图。

杠杆点与影响点

为了理解这些诊断图，我们需要区分两个重要概念：杠杆点和影响点（或异常点）。

考虑一个散点图，其中大部分数据在预测变量和响应变量之间呈现出清晰的线性关系。图中有四个用橙色标记的点，我们将逐一讨论。

第一个点位于左下角，处于数据云的中心。它既没有很高的杠杆值，也没有很大的影响力。

什么是杠杆值？

杠杆值衡量的是一个数据点因其预测变量值远离数据中心而潜在影响回归线位置的能力。

想象回归线必须穿过预测变量和响应变量的均值点，这个均值点就像一个支点。数据点距离这个支点越远，就像在杠杆上离支点越远的位置施加力，其“撬动”回归线的潜力就越大，即杠杆值越高。

因此，杠杆值的概念仅取决于数据点的预测变量值距离其中心有多远。

什么是影响力？

影响力则衡量一个数据点实际对回归线估计产生的影响大小。一个点可以有高杠杆值，但如果它本身紧密遵循由其他数据点确立的回归关系，那么它可能选择不“行使”这种杠杆力，从而对最终模型的影响很小。

例如，图中右上角的点，虽然预测变量值很大（高杠杆值），但它恰好落在由其他数据点暗示的回归线上。如果移除这个点，回归线的变化很小。因此，这个点具有高杠杆值但低影响力。

高杠杆值且高影响力的点

现在，看看右下角的点。这个点不仅预测变量值远离中心（高杠杆值），而且其响应变量值也严重偏离由其他数据点暗示的趋势。如果将其纳入回归分析，它会对回归线的斜率和位置产生巨大影响。因此，这个点同时具有高杠杆值和高影响力。

低杠杆值但可能是异常值的点

最后，观察左上角的点。它的预测变量值处于数据中心（低杠杆值），但其响应变量值却严重偏离整体趋势。虽然它是一个异常值，但由于杠杆值低，它对回归线整体形状（斜率）的影响，远小于那个高杠杆值且高影响力的右下角点。它主要影响的是截距估计。

在本讲座中，我们将详细探讨杠杆值和影响力的概念，并学习如何使用数据来诊断它们。

总结

本节课我们一起学习了回归诊断的基础知识。我们回顾了线性模型和残差的定义，演示了如何在R中快速生成诊断图形。更重要的是，我们区分了杠杆点（由预测变量的极端值引起）和影响点（实际对模型参数估计产生较大影响的点）这两个关键概念。理解这两者的区别对于正确解释诊断图、识别潜在问题数据点至关重要。在接下来的课程中，我们将深入探讨具体的诊断统计量和图形。

029：残差与诊断（第二部分）🔍

在本节课中，我们将深入探讨回归诊断中的核心概念：离群值和高杠杆点。我们将学习如何更精确地定义和识别它们，并介绍一系列实用的诊断工具，帮助你评估数据点对模型的影响。

概述

上一节我们介绍了残差的基本概念。本节中，我们将重点关注如何利用残差及其他诊断指标，系统性地识别和处理回归模型中的离群值与高杠杆点。我们将使这些术语的定义更加精确，并学习如何运用R语言中的诊断工具。

离群值与杠杆点：更精确的定义

简单地将一个点称为“离群值”是非常模糊的。在本讲中，我们将尝试使这些术语更加精确。

离群值可能由真实过程或虚假过程导致。如果是虚假过程（如数据录入错误），我们通常希望将其从拟合模型中移除；如果它是真实过程的一部分，我们则不应简单地忽略它。

我们之前通过图示讨论过，离群值可能符合回归关系，也可能不符合。有多种方法可以探查离群值，以评估它们对模型的实际影响或潜在影响。

诊断指标清单

R语言中的 ?influence.measures 命令提供了一个相当完整的默认诊断指标清单。以下是你可以使用的一系列潜在影响力度量指标。

学生化残差与标准化残差

学生化残差和标准化残差是残差的不同变体。残差的一个问题是它们具有单位特异性。例如，如果因变量y以磅为单位测量，那么残差也以磅为单位。这使得它们在不同设置下不一定具有可比性，你无法简单地设定一个阈值（比如“残差大于4即为离群值”）来判断。

为了更有效地判断离群值，我们希望对残差进行标准化。学生化残差和标准化残差就是尝试将残差除以一个适当的标准误度量。它们之间的细微区别在于计算标准误时是否剔除了当前数据点（即内部标准化与外部标准化），这会影响标准化残差是否遵循t分布。

我个人认为区分内部与外部标准化的残差帮助不大。但我喜欢使用标准化残差，因为它们处于一个更好的尺度上。不过，我倾向于将它们与数据的整体分布云图进行比较，而不是设定严格的阈值。

帽子值

帽子值是另一个有用的诊断指标，它衡量的是杠杆作用（潜在影响力）。回顾之前的图示：如果一个点远离其他X值的聚集区，它就具有高杠杆值。但这并不表示它一定会对模型产生实际影响。杠杆值仅仅衡量该点在X空间中的偏远程度。在单变量回归中这很容易可视化，但在多变量回归中则很困难，此时帽子值就变得非常有用。

帽子值的一个巨大用途是发现数据录入错误。你可以检查数据集中每一行（对应模型中的一个数据点）的帽子值。如果某些帽子值异常大，不仅意味着它们可能影响模型，还提示你应该检查该行数据是否存在录入错误。

影响力度量

接下来我们看看如何度量实际影响力，而不仅仅是潜在影响力。几乎所有的影响力度量都遵循以下思路：

  1. 剔除该数据点。
  2. 重新拟合模型。
  3. 比较该点存在与否时，模型某个方面的变化。

以下是具体的度量指标：

  • DFFITS：对于每个数据点，衡量在包含与不包含该点时，其对应X值处的拟合值发生了多大变化。
  • DFBETAS：专门查看斜率系数。衡量在包含与不包含该特定数据点时，每个回归系数发生了多大变化。对于一个线性回归模型（包含截距和一个斜率项），DFBETAS将是一个 2 x n（数据点数量）的矩阵。
  • 库克距离：衡量系数整体变化的综合指标，在某种意义上汇总了DFBETAS的信息。

PRESS残差

rstandard 函数返回普通残差，而有一种巧妙的方法可以获取所谓的“预测残差”（PRESS残差）。其思想同样是：剔除一个数据点，然后查看在该X值处预测的Y值（基于剔除该点后拟合的模型）与实际Y值之间的差异。这些PRESS残差的平方和类似于留一法交叉验证误差。

在线性回归中，计算这些指标实际上不需要重新拟合模型，可以利用线性代数关系快速计算。例如，将普通残差除以 (1 - 帽子值) 即可得到PRESS残差。

PRESS残差在检测离群值和确定影响力方面不如其他指标具体，但它更多地用于评估模型拟合优度等。

诊断工具的使用建议

在进入计算机实验之前，我先给出一些关于如何使用这些工具的建议。

首先，我不太喜欢简单的阈值规则（例如“超过某个值就宣布为离群值”）。也许对于大型数据集和自动化系统，你需要采用这种方法。但如果你正在深入进行数据分析，最好不要用简单的阈值规则来宣布离群值，而应采取更全面的方法。

  • 观察残差集合相对于其他残差的分布。
  • 观察帽子值集合相对于其他帽子值的分布。
  • 一些诊断指标没有有意义的绝对尺度，它们可能依赖于样本量等因素。

所有这些度量的一个重要方面是，它们以不同的方式探查数据，以寻找影响力和判断是否为离群值。你应该像医生诊断病情一样看待它们——通过多种方式探查你的数据，以发现可能存在的问题。

常用诊断图

以下是几种常用的诊断图：

残差 vs. 拟合值图：这是最常做的残差图。你需要在图中寻找任何系统性模式。如果存在系统性模式，则表明模型存在某种失拟。

残差Q-Q图：用于检验残差的正态性假设。

杠杆值图：用于识别具有高杠杆值的点。你应该专门检查这些数据点所在的行，它们可能是数据录入错误。

影响力度量图：这些图能直接回答核心问题：如果移除这个数据点，会对拟合模型、系数等产生多大影响？

总结

本节课中，我们一起学习了如何更精确地定义回归模型中的离群值和高杠杆点。我们介绍了一系列诊断工具，包括学生化残差、帽子值、DFFITS、DFBETAS、库克距离和PRESS残差，并讨论了如何在实际分析中综合运用这些工具和诊断图。通过接下来的示例练习，你将掌握如何将这些方法融入自己的数据分析流程中。

030：残差与诊断（第三部分）🔍

在本节课中，我们将通过模拟实验来理解几种回归诊断指标的工作原理，并学习如何解读R语言默认生成的诊断图。我们将重点关注杠杆值、DFBETA和残差图，这些工具能帮助我们识别数据中的异常点、高影响力点以及模型可能存在的问题。

通过模拟实验理解诊断指标

上一节我们介绍了杠杆值和DFBETA等概念，本节中我们来看看如何通过模拟数据来直观感受这些指标的作用。

实验一：不相关数据中的离群点

首先，我们观察一个由大量不相关的数据点构成的云团，并在其中加入一个明显的离群点 (10, 10)。

以下是生成和绘制该数据的R代码：

# 生成独立的正态随机数对
x <- rnorm(100)
y <- rnorm(100)
# 添加离群点 (10, 10)
x <- c(x, 10)
y <- c(y, 10)
# 绘制散点图
plot(x, y)

在这个例子中，仅仅因为这个离群点的存在，数据会显示出很强的相关性估计，否则相关性应估计为零。

现在，让我们查看几个诊断指标的值。

以下是相关诊断指标的计算结果：

  • DFBETA值：离群点 (10, 10) 的DFBETA值比其他所有点高出几个数量级。
  • 杠杆值：该点的杠杆值也远大于其他点。杠杆值介于0和1之间，这个点的值明显异常。

这些指标都能清晰地识别出这个离群点。

实验二：回归线上的离群点

接下来，我们看另一个例子。这里的数据基本沿一条直线分布，但同样加入了一个远离数据云团的点。不过，这个点恰好落在了回归线上。

以下是相关诊断指标的计算结果：

  • DFBETA值：这个离群点的DFBETA值仍然较大，表明它对拟合有一定影响，但不像第一个例子中那样具有压倒性的差异。
  • 杠杆值：该点的杠杆值比大多数其他点大十倍左右。原因在于，这个点的X值超出了其他数据的范围，因此具有高杠杆值。但由于它符合回归关系，所以对系数估计（DFBETA）的影响相对较小。

这个例子说明，高杠杆点不一定对模型参数有巨大影响。

残差图的重要性：Anscombe数据集示例

为了理解我们为何要绘制残差图，让我们看看Anscombe构造的著名数据集。这个例子完美展示了残差图如何放大模型潜在的问题。

我们可以从其网站下载该数据。其两两散点图除了显示大量重叠外，没有展示出任何有趣的信息。

以下是拟合线性模型和查看P值的R代码：

# 假设数据已加载为 `anscombe`
fit <- lm(y1 ~ x1 + x2 + x3 + x4, data = anscombe)
summary(fit)

查看模型摘要会发现，每个系数的P值都高度显著。如果不进一步检查，我们可能会认为模型拟合得很好。

然而，问题在于残差图。在多元回归中，我们无法像简单线性回归那样将残差对单个X作图，通常选择绘制残差与拟合值的图。

以下是绘制残差与拟合值图的R代码：

plot(fitted(fit), residuals(fit), xlab = "Fitted Values (Y_hat)", ylab = "Residuals")

在这个特定实例中，残差图揭示出一个非常巧妙的图案。如果不看残差图，我们完全无法发现这个清晰的系统模式，它表明了模型的严重缺陷。回归分析的目标是建模系统性的部分，而将无法解释的部分视为噪声。但像图中这样明显的系统性图案，本应被模型捕捉。

Anscombe构造这个数据集的目的，正是为了说明残差图的作用：它们能放大并揭示模型拟合不佳的方面。

解读R语言的默认诊断图

现在让我们回到瑞士数据集，解读R语言默认生成的一系列诊断图。课程开始时我们展示过这些图，现在来看看如何解释它们。

R的 plot() 函数对线性模型对象会生成四张诊断图。

以下是四张默认诊断图的解读：

  1. 残差与拟合值图：此图与我们刚才用Anscombe数据绘制的图类似。目标是寻找任何系统性模式。例如，如果残差随拟合值增大而扩散或收缩，可能暗示异方差性。在此数据集中，此图看起来没有太大问题。
  2. 正态Q-Q图：此图专门用于评估误差项的正态性假设。点越接近对角线，表明残差越接近正态分布。
  3. 标准化残差位置-尺度图：此图绘制标准化残差（使其在不同实验间更具可比性，类似于t统计量）的平方根与拟合值的关系。它同样用于检查异方差性或其他模式，只是换了一种尺度。
  4. 残差与杠杆图：此图将（标准化）残差与杠杆值对应起来。目标是检查高杠杆点是否伴有异常大或小的残差。例如，一个点可能杠杆值很高但残差很小（因为它强拉了回归线穿过自己），也可能杠杆值和残差都很高。

在本次分析的瑞士数据集中，这些诊断图看起来都没有严重问题。理想情况下，诊断图应具备交互功能，例如将鼠标悬停在点上可以显示该点的详细信息（如观测编号），这有助于进一步调查。我们将在后续的“数据产品”课程中讨论如何创建此类交互图。

总结

本节课中，我们一起学习了如何利用模拟实验理解杠杆值和DFBETA等诊断指标。我们通过Anscombe数据集的经典例子，深刻认识到绘制并检查残差图对于发现模型系统性缺陷至关重要。最后，我们系统解读了R语言默认提供的四张回归诊断图，包括残差与拟合值图、Q-Q图、尺度-位置图以及残差与杠杆图。掌握这些诊断工具，并将其纳入未来的回归分析流程中，将极大地帮助您评估模型假设、识别异常点，并最终建立更稳健的回归模型。

031：模型选择（第一部分）🔍

在本节课中，我们将要学习在回归背景下，当模型包含多个变量时如何进行模型拟合。我们将重点讨论模型选择的基本原则，并理解在构建解释性模型时，简约性和可解释性的重要性。

上一节我们介绍了回归模型的基础，本节中我们来看看当模型包含多个变量时需要注意的事项。

首先，需要区分多元回归与数据科学专业中预测课程所涉及的内容。在预测任务中，我们通常不太关心模型的可解释性，因此更能容忍复杂的模型、大量的交互项，以及使用自动化搜索算法和规则来针对特定损失函数寻找最佳预测模型。

相比之下，在许多回归分析的应用场景中，我们希望模型简单且易于解释。这意味着我们将高度重视简约性原则，即寻找能够解释我们所观察现象的最简单模型。这个原则类似于奥卡姆剃刀定律：在其他条件相同的情况下，最简单的解决方案很可能是正确的。在这里，我们倾向于选择简单的模型而非复杂的模型。

这种偏好会产生广泛的影响。例如，如果一个回归模型中的某个变量只能解释极少量的额外变异，但严重损害了模型的可解释性，我们通常会选择省略这个变量，即使它确实能解释一小部分正面的变异性。

为了帮助我们思考模型构建，可以引用Scott Zeger的一个观点：将模型视为观察数据的透镜。从这个角度看，模型本身没有对错之分。任何能让你从数据中学到有价值且真实信息的模型，就是“正确”的模型。透镜的比喻很贴切：有时你需要显微镜来观察超精细的细节，有时则需要望远镜来获取宏观层面的信息。合适的模型就像合适的透镜，能突出显示那些能揭示真实情况的特征。

在特定情况下，例如当我们关注某种处理效应时，模型可能存在无数种出错的方式，导致我们得出错误结论。在本讲座中，我们将探讨包含不必要回归变量或遗漏重要回归变量时会发生什么。

Donald Rumsfeld曾提出一个关于认知的分类，虽然备受争议，但对理解回归模型构建很有启发：

  • 已知的已知：我们知道并且已经记录下来的回归变量，应该直接纳入模型。
  • 已知的未知：我们知道应该纳入回归模型，但不幸的是我们没有收集到的变量。
  • 未知的未知：我们应该纳入回归模型，但我们既没有收集到，也不知道它们是否重要的变量。

这个分类有助于我们在构建回归模型时进行系统的思考。

以下是关于如何处理这些不同类型变量的一些通用规则：

遗漏重要变量（已知的未知或未知的未知）
这通常会导致偏差。如果你感兴趣的某个回归变量（如处理效应）与一个被遗漏的变量相关，那么你对感兴趣变量的估计就可能不准确。防止这种情况的一种方法是尝试随机化。如果被遗漏的变量与你感兴趣的变量不相关，那么是否包含它通常影响不大。这就是为什么A/B测试和临床试验如此有效。然而，如果存在太多混杂因素，即使随机化也可能无济于事。

包含不必要变量
包含不重要的变量通常不会引入偏差。例如，如果你生成一堆与任何事物无关的标准正态随机数并放入回归模型，它们不会导致偏差。然而，这确实会增大标准误。向回归模型中引入不必要的变量会膨胀标准误，我们可以实际测量其影响程度。

此外，像R平方这样的模型统计量，并不一定是评判包含不同回归变量的模型优劣的好方法。原因是，随着你放入更多回归变量，R平方保证会上升，无论这些变量是否有用。同样，在回归中纳入更多回归变量时，残差平方和或均方误差也保证会下降。

以下是一个示例，我模拟了一组数据，并逐步将生成的随机正态向量作为回归变量加入模型，结果显示R平方持续上升，尽管加入的只是随机噪声。

# 模拟示例代码示意
set.seed(123)
n <- 100
y <- rnorm(n)
# 逐步加入随机预测变量并计算R平方
rsq <- numeric(20)
for (i in 1:20) {
  x <- matrix(rnorm(n * i), ncol = i)
  fit <- lm(y ~ x)
  rsq[i] <- summary(fit)$r.squared
}
plot(1:20, rsq, type = "b", xlab = "Number of Random Predictors", ylab = "R-squared")

本节课中我们一起学习了多元回归模型选择的第一部分核心思想。我们明确了在解释性建模中追求简约性和可解释性的目标，并将模型比作观察数据的透镜。我们探讨了Donald Rumsfeld的认知分类在回归变量选择中的应用，并分析了遗漏重要变量会导致偏差，而包含不必要变量则会增大标准误的后果。最后，我们指出R平方等指标在模型比较中的局限性，因为它们会随着变量增加而必然提升。在接下来的课程中，我们将继续探讨更具体的模型选择方法和准则。

032：模型选择（第二部分）🔍

在本节课中，我们将深入探讨回归模型选择的另一个关键方面：包含不必要变量所带来的影响，即方差膨胀。我们还将学习如何量化这种影响，并介绍一种实用的模型选择策略。

上一节我们讨论了遗漏重要变量会导致估计偏差。本节中，我们来看看相反的情况：如果在模型中包含了不重要的变量，会发生什么。

方差膨胀：一个模拟演示 📈

为了理解包含不必要变量的影响，我们通过一个模拟实验来说明。以下是核心概念：

  • 理想情况：如果加入的变量与我们关心的变量不相关，则方差膨胀很小。
  • 糟糕情况：如果加入的变量与我们关心的变量高度相关，则会导致严重的方差膨胀。

以下是模拟的核心代码逻辑：

# 模拟设置：样本量100，重复1000次
n <- 100
nSim <- 1000

# 生成三个独立的解释变量
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)

# 生成因变量，它只依赖于x1和随机误差
y <- x1 + rnorm(n)

# 拟合三个模型，并观察x1系数估计的标准误变化
# 模型1: 仅包含 x1
# 模型2: 包含 x1 和 x2
# 模型3: 包含 x1 和 x3

当 x2 和 x3 与 x1 独立时，包含它们对 x1 系数估计的标准误影响微乎其微。然而，如果我们让 x2 和 x3 依赖于 x1（例如 x2 = x1 + rnorm(n)），重复上述模拟，会发现 x1 系数的标准误显著增大，有时甚至膨胀数倍。

结论：包含与核心变量高度相关的冗余变量，会严重增加其系数估计的不确定性（方差）。例如，在模型中同时放入舒张压和收缩压（两者高度相关），若不必要，会不当地增大舒张压效应的标准误。

方差膨胀因子

我们无法直接获知真实的误差方差，但可以通过比较来评估方差膨胀的程度。一个常用的指标是方差膨胀因子。

方差膨胀因子 衡量了当前模型中某个系数的方差，与所有其他变量均与该变量不相关的理想情况下方差的比值。其公式基于线性代数推导，但在实践中，我们通常使用统计软件直接计算。

以下是使用R语言计算VIF的步骤：

# 加载必要的包
library(car)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_32.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_34.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_35.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_37.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_38.png)

# 使用瑞士生育率数据，拟合包含所有变量的模型
data(swiss)
fit <- lm(Fertility ~ ., data = swiss)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_40.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_41.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_42.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_43.png)

# 计算方差膨胀因子
vif(fit)
# 更直观地查看对标准误的膨胀影响（取平方根）
sqrt(vif(fit))

对于VIF结果的解读：

  • VIF = 1：表示该变量与其他变量完全不相关，是理想情况。
  • VIF > 1：表示存在相关性。例如，VIF=4意味着标准误是理想情况下的2倍（sqrt(4)=2）。
  • VIF值很大（如>5或10）：表明存在严重的多重共线性，需要关注。

在瑞士数据中，Examination和Education的VIF很高，因为它们测量的是相似的概念（教育水平），彼此高度相关。而Infant.Mortality的VIF接近1，说明它与其他变量基本无关。

对残差方差估计的影响

除了回归系数，模型选择还会影响另一个重要参数——残差方差σ²的估计。

以下是其影响规律：

  1. 模型欠拟合（遗漏重要变量）：残差方差估计是有偏的。因为被遗漏变量的系统性影响被错误地归入了随机误差。
  2. 模型正确拟合或过拟合（包含不必要变量）：残差方差估计是无偏的。
    • 但是，如果包含了不必要变量，方差估计本身的变异性（方差）会增大，使得估计更不可靠。

核心规律总结：对于参数估计（无论是回归系数β还是残差方差σ²），遗漏变量导致偏差，而包含冗余变量导致方差增大/估计精度下降。这再次凸显了良好实验设计（如随机化）的重要性，因为它能平衡已知和未知的混杂因素，减少对复杂模型调整的依赖。

自动化模型选择与嵌套模型

当变量很多，特别是考虑交互项和多项式时，模型选择的空间会爆炸式增长。完全的自动化搜索（如基于AIC、BIC）更属于机器学习的范畴。对于解释性回归，我们通常依赖分析者、数据和科学背景之间的交互。

有一种自动化策略在特定场景下非常有用：嵌套模型比较。

当我们特别关心某一个核心变量（如处理效应）时，可以构建一系列模型：

  1. 仅包含核心变量。
  2. 在模型1基础上，加入一个可能重要的协变量。
  3. 在模型2基础上，再加入另一个协变量。
  4. ...

每个后续模型都完全包含了前一个模型的变量，这就是“嵌套”关系。我们可以使用似然比检验（F检验） 来评估新增的一组变量是否提供了显著额外的解释力。

以下是使用R进行嵌套模型比较的示例：

# 使用瑞士数据，假设我们关注Agriculture的影响
fit1 <- lm(Fertility ~ Agriculture, data = swiss)
fit2 <- lm(Fertility ~ Agriculture + Examination + Education, data = swiss)
fit3 <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data = swiss)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_115.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_116.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_117.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_118.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_119.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_120.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/849e472ba79624ded97b40e4795fbbf3_121.png)

# 使用anova函数进行嵌套模型比较
anova(fit1, fit2, fit3)

anova函数的结果会列出模型的残差平方和、自由度，并计算F统计量和p值，从而判断每次新增的变量集是否显著改善了模型。

注意：此方法仅适用于严格的嵌套模型。对于非嵌套模型的比较，需要用到信息准则（AIC/BIC）等其他方法。

总结与思考 🎯

本节课我们一起学习了模型选择的几个核心要点：

  1. 方差膨胀：包含与核心变量高度相关的冗余变量，会显著增大其系数估计的标准误，降低估计精度。方差膨胀因子是诊断这一问题的有用工具。
  2. 全面影响：模型选择不仅影响回归系数，也影响残差方差的估计——遗漏变量导致偏差，冗余变量导致方差增大。
  3. 实用策略：当分析聚焦于一个核心变量时，嵌套模型比较是一种结构清晰、统计严谨的模型选择策略。
  4. 设计优先：良好的实验设计（如随机化、区组设计、交叉设计）能从源头上减少混杂，简化模型选择的负担。

最后，请记住统计学家乔治·博克斯的名言：“所有的模型都是错的，但有些是有用的。” 我们的目标不是找到“正确”的模型，而是构建一个能够帮助我们理解数据、揭示有用真相的“有用”透镜。

本节课到此结束，我们下节课再见！

033：模型选择（第三部分）🔍

在本节课中，我们将要学习模型选择对回归模型中另一个关键参数——残差方差的影响，并探讨自动化模型选择、实验设计的作用，以及一种实用的嵌套模型检验方法。

残差方差的影响 📊

上一节我们讨论了包含或遗漏变量对回归系数的影响。本节中，我们来看看这些操作对回归模型中另一个估计参数——残差方差 σ² 的影响。

在假设误差独立同分布且线性模型形式正确的前提下，我们可以从数学上描述遗漏必要变量或包含不必要变量的影响。

其规律与我们之前讨论的系数影响类似：

  • 如果模型欠拟合（即遗漏了重要变量），那么方差估计是有偏的。这是因为我们将那些本可以由被遗漏的协变量系统解释的变异，错误地归因于随机误差。
  • 如果模型正确拟合或过拟合（即包含了所有正确项或额外的不必要项），那么方差估计是无偏的。然而，如果我们包含了不必要变量，方差估计值本身的方差会增大，导致估计的可靠性下降。

简而言之，对于方差估计：遗漏变量导致偏差，包含不必要变量导致估计可靠性降低。

关于自动化模型选择 🤖

现在让我们来谈谈广义的自动化模型选择。需要指出的是，自动化模型选择如今更多属于机器学习的范畴。即使对于相对简单的线性回归模型，当你开始考虑交互项、多项式项（如预测变量的平方）时，需要搜索的模型项空间也会迅速爆炸式增长。

如果你有大量预测变量并希望缩减这个空间，可以使用主成分分析等因子分析技术来降低协变量空间的维度。然而，这些方法也有代价：你得到的主成分或因子可能不如原始数据易于解释。这部分内容在多变量分析或机器学习课程中会有更深入的探讨。

在本课程中，我们主要考虑预测变量数量相对较少的情况，模型选择将是一个结合分析师、数据和科学背景的高度交互过程。

实验设计的作用 🧪

我想提到的另一点是，良好的实验设计通常可以消除大量模型讨论的必要性。我们之前讨论过随机化如何能有效避免许多问题，例如使我们感兴趣的变量与那些我们不感兴趣甚至未知的干扰变量无关。

然而，设计还有其他方面可以达到相同目的。例如：

  • 分层与区组内随机化：经典的例子是田间作物实验。如果你要试验不同的种子，你可能会根据田地的不同区域进行“区组”划分，然后将不同种子随机分配到这些区组中。
  • 交叉设计：在生物统计学中非常常见。例如，比较两种阿司匹林时，不是让两组不同的人分别服用一种，而是让同一组人在不同时期先后服用两种阿司匹林（通常会随机化服药顺序），这样每个人都可以作为自身的对照，从而控制个体间的固有差异。

实验设计是一个相当广泛的领域。我想强调的更广泛的要点是：深思熟虑的良好实验设计，通常可以消除在单纯观察性数据建模过程中必须考虑的许多主要问题。

一个实用的方法：嵌套模型检验 🔗

最后，我想介绍一种我非常喜欢且觉得非常有用的自动化模型搜索技术：查看嵌套模型。

我通常对某个特定变量（例如某种处理效应）非常感兴趣，但又担心处理组在其他变量（如年龄、性别）上不平衡，这些变量可能会影响结果。这时，嵌套模型分析就很有用。

以下是进行嵌套模型检验的思路：

  1. 拟合一个只包含核心变量（如“处理”）的模型。
  2. 拟合一个包含核心变量和一个潜在混杂变量（如“年龄”）的模型。
  3. 拟合一个包含核心变量和多个潜在混杂变量（如“年龄”和“性别”）的模型。
    以此类推。

嵌套模型是指每一个后续模型都完全包含了前一个模型的所有项。这为检验每个新增的变量集是否必要提供了一种非常简便的方法。

以下是如何在R语言中进行嵌套模型检验的示例代码。我们使用swiss数据集，假设我们关注的核心变量是Agriculture（农业）。

# 拟合三个嵌套线性模型
fit1 <- lm(Fertility ~ Agriculture, data = swiss)
fit3 <- lm(Fertility ~ Agriculture + Examination + Education, data = swiss)
fit5 <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data = swiss)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_111.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_112.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_113.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_114.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_115.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_116.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_117.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/jhu-ds-7/img/443ef12ddaca7b7d4661a06bab6aa2bc_119.png)

# 使用anova函数进行嵌套模型比较
anova(fit1, fit3, fit5)

anova函数会列出模型，并给出残差平方和、自由度等。通过计算F统计量和对应的P值，它可以检验新增的变量集是否提供了显著的额外解释力。

例如，输出可能会显示，在模型1（仅有Agriculture）基础上加入Examination和Education是必要的；同样，在模型3基础上再加入Catholic和Infant.Mortality也是必要的。

要点：当你感兴趣的分析思路自然形成一系列逐步复杂的嵌套模型时（尤其当你专注于某个特定变量时），嵌套模型检验是一个合理且有用的方法。但这种方法仅适用于模型是嵌套的情况。对于非嵌套模型的比较，你需要使用信息准则等其他方法，这通常在我们的预测课程中涉及。

总结 📝

本节课中我们一起学习了：

  1. 模型选择对残差方差的影响：遗漏重要变量会使方差估计有偏，而包含不必要变量则会降低估计的精确性。
  2. 自动化模型选择的复杂性：在变量较多时，模型空间巨大，需借助其他技术或依赖于具体情境的交互式选择。
  3. 实验设计的重要性：良好的随机化、分层、交叉设计等可以从源头减少对复杂模型调整的依赖。
  4. 嵌套模型检验：当分析围绕一个核心变量展开，并逐步控制其他变量时，这是一种检验变量必要性的有效且直观的方法。

最后，我们必须意识到，模型几乎总是不完美的。正如统计学家乔治·博克斯的名言：“所有的模型都是错的，但有些是有用的。” 我们的目标不是找到一个完全正确的模型，而是找到一个能帮助我们从中学习数据真相的有用透镜。

希望本讲座为你提供了一些可用的模型选择技巧，并加深了对模型构建后果的理解。我们下节课再见。

034：广义线性模型

在本节课中，我们将学习广义线性模型。我们将从之前讨论的线性可加响应模型，转向更广泛的广义线性模型类别。这些模型能够处理线性模型可能难以应对的一些情况，但可能会牺牲线性模型所具有的一些数学上的简洁性。

在开始讨论广义线性模型之前，我们先回顾一下线性模型，了解其优点和局限性。

🔍 线性模型的回顾与局限

线性模型是应用最广泛的统计技术之一。然而，它并非没有局限性。可加性假设在某些情况下难以应用。

例如，对于二元数据，将其视为可加的会显得很奇怪，其误差结构在这种情况下会变得难以解释。

此外，如果数据是严格为正的，那么使用可加的高斯误差会带来问题，因为这允许预测值为负的概率存在。虽然有时这并非大问题，但如果正态分布对负值赋予了很高的正概率，而你的响应变量已知必须为正，那么这对模型来说就是有问题的。

🔄 数据转换的替代方案

为了解决上述问题，一种常见的做法是对数据进行转换。

对于严格为正的结果变量，一个常见的转换是取其自然对数。自然对数可能是可解释性最强的转换之一。当然，也存在许多其他转换方法，试图使数据更接近正态分布。

对于二项分布数据，人们常采用所谓的反正弦平方根变换，但这通常会严重破坏模型系数的可解释性。

相比之下，对数变换在很多情况下反而能使系数更具可解释性，我们将单独讨论这一点。

🎯 采用广义线性模型的理由

除了对结果变量进行转换或用线性模型近似外，采用广义线性模型还有另一个原因：直接在数据记录的尺度上建立模型，无需转换，这能更好地尊重数据的已知特性。

如果我们有二元数据，一个真正尊重数据二元性且无需我们进行转换的模型，会显得更加自然和合理。

此外，最常用的自然对数变换不适用于负值或零值。虽然有一些修正方法，但这会损害该变换的一些优良特性，特别是其良好的可解释性。

📜 广义线性模型的起源与构成

广义线性模型源于1972年Nelder和Wedderburn的一篇著名论文。

一个广义线性模型包含三个组成部分：

  1. 随机成分：描述随机性的分布必须来自一个特定的分布族，称为指数族。这是一个庞大的分布族，包括正态分布、二项分布和泊松分布等。
  2. 系统成分：即所谓的线性预测器。这是我们建模的部分，在线性模型中我们已经非常熟悉。它由回归协变量和系数构成。
  3. 连接函数：连接函数将指数族分布的重要均值与线性预测器联系起来。

因此，我们需要三样东西：一个属于指数族的分布、作为线性预测器的系统成分，以及一个连接两者的连接函数。

📐 示例一：线性模型作为特例

让我们从一个最熟悉的例子开始：线性模型。

在这个例子中：

  • 随机成分：我们假设响应变量 Y_i 服从正态分布，均值为 μ_i，方差为 σ^2。正态分布属于指数族。
  • 系统成分：线性预测器 η_i 定义为协变量 X_i 与其系数 β 的线性组合。
  • 连接函数：我们使用恒等连接函数，它规定正态分布的均值 μ_i 恰好等于线性预测器 η_i。

这最终导出了我们一直讨论的线性模型：Y_i = Σ X_i β + ε_i。我们只是以另一种方式写出了线性模型：我们没有说误差项服从正态分布，而是说 Y 本身服从正态分布（这是前一种设定的结果），并单独指定了线性预测器。

虽然用可加误差的求和形式写出来似乎更简单，但当我们转向泊松或二项分布等不同设定时，这样做的原因就会变得更加明显。

🪙 示例二：逻辑回归（二项分布）

接下来，我们看看广义线性模型中可能最有用的变体之一：逻辑回归。

在这个设定中，我们的数据是0或1，即二元数据。因此，假设它们来自正态分布没有太大意义。自然地，描述抛硬币或0/1结果的唯一可用分布是伯努利分布。

  • 随机成分：我们假设结果变量 Y_i 服从伯努利分布，其成功概率（即期望值）为 μ_i。我们将数据建模为一连串的抛硬币结果，只是每次正面朝上的概率 μ_i 可能不同，且不一定等于0.5。
  • 系统成分：线性预测器 η_i 保持不变，仍然是协变量与其系数的线性组合。
  • 连接函数：在这种情况下，最著名和最常用的是逻辑连接函数或对数几率。我们通过对几率取自然对数，将均值（正面概率）与线性预测器联系起来。几率定义为 μ / (1 - μ)。

这意味着，我们将分布的均值（成功概率）进行转换，然后在这个转换后的尺度上，协变量及其系数（我们线性模型的标准部分）才会出现。

请注意，我们转换的是分布的均值，而不是数据 Y 本身。这是广义线性模型的一个巧妙之处。我们假设数据来自一个具有特定成功概率的“硬币”，我们以特定方式转换这个概率，使其与我们的协变量和系数相关联。

我们可以通过逆逻辑函数（或称Sigmoid函数）从对数几率 η 返回到均值 μ。这样，我们就可以写出似然函数，并通过最大化这个似然函数来获得参数估计。

🔢 示例三：泊松回归

另一个重要例子是泊松回归。

  • 随机成分：假设 Y_i 服从泊松分布，均值为 μ_i，且 μ_i > 0。泊松分布非常适合对计数数据建模，这些计数是正的且无上限的（不像二项分布计数受试验次数限制）。
  • 系统成分：线性预测器 η_i 与之前相同。
  • 连接函数：泊松回归中最常用的连接函数是对数连接函数。我们通过对均值 μ 取自然对数，将其与线性预测器 η 联系起来。同样，我们是对分布的均值取对数，而不是对数据本身取对数。反过来，我们可以通过对 η 取指数来回到 μ。

同样，我们可以写出似然函数，并通过最大化似然函数来获得参数估计。

⚙️ 模型拟合与推断

广义线性模型的参数估计通过最大化似然函数获得。其求解过程与最小二乘法有相似之处，但形式上更复杂一些，涉及权重和方差项。不过，对于大多数应用者来说，具体的拟合细节是透明的，我们更关心的是模型的解释。

与线性模型不同，广义线性模型的方程需要通过迭代算法求解。这意味着有时程序可能会失败（例如，在逻辑回归中数据全为0或全为1时）。除此之外，大部分分析对我们来说应该是熟悉的。

  • 预测响应值：使用估计的系数 β_hat 乘以回归变量，得到预测响应值。但请注意，对于逻辑回归，这是在对数几率尺度上；对于泊松回归，这是在对数尺度上。如果需要回到原始数据的尺度，需要进行相应的逆转换（例如，逻辑回归用逆逻辑函数，泊松回归取指数）。
  • 系数解释：系数的解释方式与线性回归类似，表示在保持其他变量不变的情况下，该变量每单位变化所引起的期望响应的变化。唯一的区别是，这种解释是在线性预测器的尺度上进行的（逻辑回归是对数几率尺度，泊松回归是对数均值尺度）。虽然解释稍复杂，但我们获得了在数据原始尺度上建模的好处，且无需转换结果变量。
  • 统计推断：我们失去了线性模型中那种基于t分布的简洁的封闭形式推断。但广义线性模型的推断在很大程度上是类似的，统计学家已经推导出检验系数显著性（如P值）所需的近似分布。需要注意的是，这些结果基于大样本渐近理论，因此需要较大的样本量。如果样本量很小，P值等可能不够可靠。

⚠️ 方差结构与过度离散

线性模型有一个常数方差的假设。然而，对于广义线性模型，情况有所不同：

  • 伯努利分布：一个硬币翻转的方差是 p(1-p)，即 μ(1-μ)。由于 μ 随观测点 i 变化，方差也随之变化。
  • 泊松分布：泊松分布的方差等于其均值 μ，因此方差也随 i 不同而不同。

这是一个建模假设，可以进行检查。例如，如果你的数据不符合这种均值-方差关系（例如，泊松数据的观测方差远大于其均值），就会产生问题，这被称为过度离散。

为了解决这个问题，可以使用更灵活的方差模型，即使这会失去广义线性模型的一些严格假设。在R等统计软件中，这通常通过“拟”家族来实现（例如 quasipoisson, quasibinomial）。它们提供了更灵活的方差结构，是处理过度离散的标准选项。

📝 总结

本节课我们一起学习了广义线性模型的核心概念。我们了解到，GLM通过三个核心部分——属于指数族的随机成分、线性预测器构成的系统成分以及连接两者的连接函数——扩展了线性模型，使其能够直接处理非正态分布的响应变量（如二元数据、计数数据）。

我们重点回顾了线性模型的局限性，探讨了采用GLM而非简单数据转换的理由，并通过线性模型、逻辑回归和泊松回归三个具体例子，详细阐述了GLM的构成与应用。我们还简要介绍了GLM的模型拟合、推断、预测解释以及需要注意的方差结构问题。

在接下来的课程中，我们将深入探讨两个最重要的具体案例：用于二元结果的逻辑回归和用于计数数据的泊松回归。

035：逻辑回归第1部分 📊

在本节课中，我们将要学习针对二分类数据的广义线性模型，即逻辑回归。我们将重点讨论当结果变量是0或1的随机变量时，如何建立模型。课程会以一个具体的例子——预测巴尔的摩乌鸦队（一支美式橄榄球队）的胜负——来贯穿始终，帮助你理解核心概念。

概述：什么是二分类数据？

在许多分析场景中，结果变量只能取两个值。例如，在生物统计学中，研究结束时病人可能“存活”或“死亡”；在体育分析中，比赛结果是“胜利”或“失败”。这类数据可以被视为一系列“抛硬币”实验，但每次实验的“成功”概率可能不同，并依赖于一系列协变量。

如果一系列0和1的结果是独立的，且成功概率恒定，那么总成功次数就服从二项分布。当协变量为常数时，二项分布数据可以通过二元逻辑回归来处理。本节课我们将通过讲解二元数据的情况，来覆盖二元和二项两种情况。

案例引入：乌鸦队胜负数据

作为教学案例，我们使用由课程讲师Jeff Leek整理的数据集。该数据集记录了乌鸦队的比赛信息，我们将探究乌鸦队的得分在多大程度上能预测其比赛胜负。

以下是数据中包含的关键变量：

  • ravenWinNumber: 乌鸦队获胜记为1，否则为0。
  • ravenScore: 乌鸦队得分。
  • opponentScore: 对手得分。

在本例中，我们暂时只关注乌鸦队得分 (ravenScore) 这一个预测变量。你可以想象，如果你是球队管理层，会希望用包括得分、对手得分、主客场等多种协变量来预测胜负。通过这个简单的例子理解逻辑回归原理后，你可以很容易地将知识扩展到多变量的情况。

为什么不用线性回归？

上一节我们介绍了案例数据，本节中我们来看看为什么标准的线性回归模型可能不适用于此类数据。

如果我们试图用线性回归直接建模乌鸦队获胜（记为1）的概率，即建立如下模型：
获胜概率 = β0 + β1 * 得分 + 误差
这个模型存在明显问题。线性回归通常假设误差服从正态分布，但我们的结果变量只有0和1两个值，该假设显然不成立。此外，线性模型预测出的“概率”值可能小于0或大于1，这不符合概率的定义。

虽然偶尔用线性模型快速查看趋势并非不可，但更严谨、更合适的方法是使用逻辑回归模型。

核心概念：几率、对数几率与Logit函数

既然不能直接对概率建模，逻辑回归选择对几率进行建模。几率定义为事件发生的概率与不发生的概率之比。

公式：
几率 = p / (1 - p)
其中 p 代表事件发生的概率。

几率可以转换回概率：
公式：
p = 几率 / (1 + 几率)

我们更常用的是几率的自然对数，称为对数几率或Logit。
公式：
Logit(p) = log( p / (1 - p) )

逻辑回归的巧妙之处在于：它假设对数几率与预测变量之间存在线性关系，而不是直接假设概率与预测变量线性相关。

逻辑回归模型

现在，让我们正式建立逻辑回归模型。我们不再像线性回归那样直接建模期望值，而是通过一个连接函数来建模。

对于二分类结果，其期望值就是概率 p。逻辑回归模型设定为：
公式：
p = exp(β0 + β1*x) / [1 + exp(β0 + β1*x)]
或者等价地，写成对数几率的形式：
公式：
log( p / (1 - p) ) = β0 + β1*x

这个模型可以理解为：我们将数据视为一系列成功概率不同的“抛硬币”实验，每个实验的成功概率 p 通过Logit函数与预测变量 x 线性相关。从对数几率反向计算概率的函数（即 exp(a)/[1+exp(a)]）通常被称为Logistic函数或Sigmoid函数。

模型解释

逻辑回归系数的解释非常直观。所有提到的对数均指自然对数（ln）。

  • 截距 β0：当预测变量 x 为0时，结果的对数几率。在乌鸦队的例子中，β0 是乌鸦队得0分时获胜的对数几率。exp(β0)/[1+exp(β0)] 则是模型预测的得0分时的获胜概率。当然，现实中得0分不可能获胜，这说明了模型在数据范围外的推断需要谨慎。
  • 斜率 β1：它表示预测变量 x 每增加一个单位，结果的对数几率的变化量。
    公式推导：
    Logit(p | x+1) - Logit(p | x) = [β0 + β1*(x+1)] - [β0 + β1*x] = β1
  • 指数化斜率 exp(β1)：这是优势比。它表示预测变量 x 每增加一个单位，结果的几率变为原来的多少倍。这是更常用的解释方式。

如果模型中有多个预测变量（如同时考虑得分和主客场），exp(β1) 的解释则是在保持其他变量不变的情况下，得分每增加一分，获胜的优势比。这与多元线性回归中“控制其他变量后”的解释思路一脉相承。

几率的起源与趣味解读

对数几率之所以有用，是因为对数能将加法关系转化为乘法（比率）关系，而比率是人们易于理解的概念。几率的概念源于“公平博弈”的思想。

想象一个赌局：抛一枚成功概率为 p 的硬币。如果正面朝上，你赢 X 元；如果反面朝上，你输 Y 元。为了让游戏公平（期望收益为0），需要满足：
Y / X = p / (1 - p)
这正是几率。如果设 X = 1，那么 Y 就是你为赢得1美元而愿意支付的钱数。这就是赛马等博彩中“赔率”的由来。例如，“10比1”的赔率意味着如果你下注1美元，赢了就能拿回10美元。

一个有趣的现象是，赛马场通过动态调整赔率来平衡投注资金，确保无论哪匹马赢，庄家几乎不亏钱（他们通过收取手续费盈利）。长期来看，公众投注形成的聚合赔率，竟然能相当准确地反映赛马的真实胜率，这体现了“群体智慧”和市场效率。

总结与下节预告

本节课中我们一起学习了逻辑回归的基础知识。我们了解到，对于二分类数据，直接使用线性回归存在模型假设不满足、预测值可能超出合理范围等问题。逻辑回归通过建模对数几率与预测变量的线性关系，巧妙地解决了这些问题。模型的系数，特别是经过指数化后的优势比，具有清晰直观的解释。

下一讲，我们将通过计算机实验，直观地展示如何拟合和可视化逻辑回归曲线。敬请期待。

036：逻辑回归第2部分 📈

在本节课中，我们将学习逻辑回归曲线，并探讨R语言是如何拟合数据的。我们将通过可视化工具来理解逻辑回归模型的参数如何影响曲线形状，以及模型如何通过最大似然估计来寻找最佳拟合。

逻辑回归曲线与参数调整 🔧

上一节我们介绍了逻辑回归的基本概念，本节中我们来看看逻辑回归曲线的具体形态以及其参数的意义。

逻辑回归模型的核心函数是一个S形曲线，其公式如下：

公式： P(Y=1|X) = e^(β₀ + β₁X) / (1 + e^(β₀ + β₁X))

其中，P(Y=1|X) 表示在给定预测变量 X 时，结果为1（或“成功”）的概率。β₀ 和 β₁ 是模型的参数。

我们可以使用R语言的 manipulate 包来动态观察这些参数如何改变曲线的形状。

  • β₁（斜率参数）的作用：当 β₁ 的绝对值增大时，S形曲线会变得更陡峭。如果 β₁ 为正值，曲线从左下向右上增长；如果 β₁ 为负值，曲线则从左上向右下下降（即形状翻转）。
  • β₀（截距参数）的作用：改变 β₀ 的值会使整个S形曲线沿着X轴左右平移，但不改变其基本形状。

R语言如何拟合逻辑回归曲线 🤖

理解了曲线形态后，我们来看看R语言在幕后是如何工作的。

想象我们有一组数据，其中结果变量 Y 是二元的（0或1）。在X值较小的区域，观测到的结果大多是0；在X值较大的区域，观测到的结果大多是1。当然，中间区域会存在一些混合。

R语言拟合逻辑回归的目标是：在所有可能的S形曲线（即所有 β₀ 和 β₁ 的组合）中，找到一条最能匹配数据背后概率规律的曲线。

以下是R语言实现这一目标的过程：

  1. 尝试不同曲线：系统通过调整 β₀ 和 β₁ 这两个“滑块”，生成无数条不同的逻辑回归曲线。
  2. 评估拟合优度：对于每一条尝试的曲线，计算其能“解释”当前观测数据的可能性。
  3. 选择最佳曲线：采用最大似然估计的原则，最终选择那个能使观测数据出现“可能性”最大的参数组合，其对应的曲线就是最佳拟合曲线。

这个过程就像让这些数据点固定不动，然后我们调整曲线的形状和位置，直到找到最贴合数据点分布趋势的那条S形路径。

总结 📝

本节课中我们一起学习了逻辑回归的核心机制。

我们首先通过可视化演示，理解了模型参数 β₀ 和 β₁ 如何控制逻辑回归S形曲线的位置与陡峭程度。接着，我们探讨了R语言拟合逻辑回归模型的本质：它通过最大似然估计方法，系统地搜索所有可能的参数组合，以找到那条能最优描述二元结果与预测变量之间概率关系的曲线。

最终得到的拟合函数，即上述概率公式，就是将线性预测值 (β₀ + β₁X) 转换回我们更容易理解的概率尺度。

037：逻辑回归第3部分 📊

在本节课中，我们将学习如何在实际数据上拟合逻辑回归模型，解读模型输出，并理解优势比与相对风险等核心概念。

上一节我们介绍了逻辑回归曲线的拟合原理，本节中我们来看看如何将模型应用于实际数据并解读结果。

拟合乌鸦队数据

我们使用R语言的glm函数来拟合乌鸦队（Ravens）数据。该函数的使用方式与lm函数类似，结果变量位于波浪号~左侧，预测变量位于右侧。关键区别在于需要指定family = binomial参数，这告诉glm函数我们处理的是0/1二分类数据。如果是计数数据或真正的二项分布数据，则还需要提供样本量。

默认情况下，对于二项分布和二元情况，glm会假设我们使用逻辑连接函数，这符合我们的需求。

# 拟合逻辑回归模型
model <- glm(win ~ score, data = ravens_data, family = binomial)

解读模型输出

使用summary()函数查看模型摘要，其输出格式与线性模型lm类似。

以下是模型输出的核心部分：

  • 系数：在Logit尺度上，截距约为-1.68，score（得分）的系数约为0.1。对于得分变量，我们关注其系数在Logit尺度上是否接近0，或在指数化后（优势比尺度）是否接近1。
  • 标准误、z值和p值：这些统计量的解读方式与线性模型类似，但需要认识到它们的推导方法不同。

预测概率与拟合曲线

模型给出的预测响应值被转换回概率尺度。具体过程是：将得分值乘以对应的系数估计值，加上估计的截距，然后计算 e^(该值) / (1 + e^(该值))，从而得到预测概率。

公式：P(win) = exp(β0 + β1 * score) / [1 + exp(β0 + β1 * score)]

我们得到的S型曲线只是完整曲线的一部分，因为观测数据范围有限。完整的逻辑曲线从0到1，而我们的拟合值只显示了其中一段（例如，从0到0.4）。

优势比及其解释

如果我们对乌鸦队得分的系数进行指数化，截距变为0.1864，得分系数变为1.1125。这意味着乌鸦队每多得1分，获胜的优势（而非直接概率）增加约11.25%。

我们可以使用confint函数轻松获取系数的置信区间。为了在更易解释的优势比尺度上查看，通常会对置信区间的上下限进行指数化。

# 获取系数置信区间并指数化
exp(confint(model))

在我们的例子中，得分优势比的置信区间包含1（例如0.99到1.3）。这表明，尽管我们知道得分是获胜的关键，但该系数在统计上并不显著。

模型比较与ANOVA

anova函数在逻辑回归中的工作方式与在线性模型中相同。你可以输入拟合的模型对象，也可以比较多个嵌套模型，anova会提供一系列序贯检验。

以下是anova函数的一些应用场景：

  • 当模型中只有一个预测变量（如本例的score）时，它会进行单自由度的检验。
  • 当模型包含因子变量时特别有用。summary输出的系数表会独立检验因子的每个水平，而anova可以检验整个因子变量（所有水平）是否应该作为一个整体放入或移出模型。

优势比与相对风险

在解释优势比时，请记住：

  • 优势比是概率的函数，但其本身不是概率。
  • 优势比为1表示无差异（对应Logit尺度上的0）。
  • 通常，优势比低于0.5或高于2被认为有较强的效应，但这高度依赖于具体研究背景。

相对风险是另一个常被提及的指标，它是两个概率的比值。许多人更喜欢相对风险，因为它基于更直观的概率概念。

然而，与优势不同，相对风险会施加一些相当严格的模型约束，因此基于相对风险的回归并不像逻辑回归那样普遍。一个经典的结论是：当研究事件的概率很小时，相对风险近似等于优势比。因此，有时人们会拟合优势比模型，但将其解释为相对风险。

本节课中我们一起学习了如何拟合和解读逻辑回归模型，包括使用glm函数、理解模型输出中的系数和p值、计算预测概率和优势比，以及使用anova进行模型比较。我们还讨论了优势比与相对风险的区别与联系。

在接下来的课程中，我们将探讨泊松回归，用于处理计数型数据。

038：泊松回归第1部分

在本节课中，我们将要学习泊松回归模型。泊松回归是广义线性模型（GLM）的一种，专门用于处理计数数据和比率数据。我们将从理解泊松分布的基本性质开始，探讨为何它适用于计数建模，并学习如何将线性回归的思路应用到计数数据上。

泊松分布简介

上一节我们介绍了广义线性模型的基本框架，本节中我们来看看专门用于计数的泊松模型。计数数据在应用中非常常见，例如呼叫中心的来电数量、流感病例数等。这些计数在理论上是无界的，或者说其理论上界（如世界总人口）相对于我们观察到的计数值来说非常大。

此外，数据也可能以比率或比例的形式出现，例如测试通过率，或单位时间内的事件发生率（如核泵故障率）。在公共卫生领域，一个常见的比率是发病率，即每单位人-时间风险中新发病例的数量。所有这些计数和比率数据都可以用泊松广义线性模型来处理。

泊松分布是建模计数和比率的有效工具。比率是单位监测时间内的计数，例如发病率和网络流量等常使用泊松分布建模。

泊松分布的一个常见用途是近似二项分布概率，当成功概率非常小且试验次数n非常大时。另一个我喜欢的应用是列联表数据建模，例如统计不同发色和不同眼色的人数交叉表。泊松模型为这类数据提供了非常优雅的建模框架。

以下是泊松分布的概率质量函数。单位时间的计数率为 λ，总时间为 t。如果 X 服从均值为 tλ 的泊松分布，那么其期望值 E[X] = tλ。因此，我们对比率的自然估计是计数除以总时间，即 X/t。一个有用的性质是，这个比率估计的期望值 E[X/t] 恰好等于我们想要估计的比率 λ。

泊松分布的方差等于其均值，即 Var(X) = tλ。这是模型的一个假设，我们可以进行检验，并在其不成立时寻求解决方案。

另一个有趣的性质是，当泊松分布的均值变大时，它会趋近于正态分布。这可以通过模拟来展示：随着泊松分布均值的增大，其分布形状越来越接近正态分布。

我们可以通过数学推导或模拟来验证泊松分布的均值与方差相等这一性质。

一个实例：网站流量分析

让我们来看一个杰夫·利克的网站流量实例。网站地址是 www.biostat.jhsph.edu/~jleek。在这个例子中，泊松分布的均值被解释为每天的网站点击量，即我们的时间单位 t = 1 天。如果我们想将估计的比率解释为每小时点击量，则需要设 t = 1/24。

以下是数据预览。我们将日期从标准字符格式转换为儒略日（自1970年1月1日以来的天数），这样它就是一个简单的计数，便于分析。

数据集中包含日期、访问次数、来自“简单统计”博客的访问次数以及儒略日。

下图绘制了数据集，X轴是儒略日，Y轴是访问次数。

在上一讲中，我们讨论了用线性回归直接建模计数数据（或二值数据）的缺陷。对于计数数据，直接使用线性模型存在一些问题。然而，正如我们之前看到的，当计数的均值非常大时，这种担忧会减少，因为泊松分布会趋近于正态分布。如果计数值非常大，使用线性模型的反对意见就会少很多。

让我们建立符号：以访问次数 NH 作为结果变量，儒略日 JD 作为预测变量，构建一个线性回归模型。我们可以绘制拟合线，但可以看到存在一些曲线关系，或许应该加入二次项。这是我们的第一种方法，虽然可能不是最佳选择，且线性模型的解释性不佳。在接下来的部分，我们将探讨如何调整线性模型以获得更好的解释性。对于网站点击量这类数据，我们通常希望从相对尺度来思考，而线性模型是在线性加性尺度上处理的。

对数变换模型

我们首先可以尝试对结果变量取自然对数。模型如下：

log(NH) = β0 + β1 * JD + ε

让我们讨论一下对数变换的作用。一个随机变量对数的期望值的指数，我称之为总体几何均值。其原因是，样本的几何均值是数据乘积的 n 次方根。对这个几何均值取对数，就得到了对数数据的算术均值。因此，几何均值就是对数据取对数后的算术均值的指数。

当我们在线性回归中对结果变量取自然对数时，指数化的回归系数就可以用几何均值来解释。例如：

  • exp(β0) 是第0天（起点日）的估计几何均值点击量。需要重申的是，这里的截距意义不大，因为1970年1月1日并非我们关心的日期。为了使截距可解释，我们应该从所有日期中减去数据集中的最早日期，这样截距就代表了数据集第一天的几何均值点击量。
  • exp(β1) 是估计的每天几何均值点击量的相对增加（或减少）量。

这里存在一个问题：如果计数为零，则无法取对数。常见的处理方法是加一个常数，通常是加1，即对 log(结果 + 1) 进行建模。

如果我们对 log(NH + 1) 与儒略日进行线性回归，会得到截距（如前所述，此处意义不大）和斜率。在指数尺度上，斜率约为1.002。这意味着我们的模型估计每天网站流量增加约0.2%。这是一个很好的解释。如果加入其他协变量，这个0.2%就是在保持其他协变量不变的情况下，每天的估计增加量。

本节课中我们一起学习了泊松回归的基础知识。我们了解了泊松分布适用于计数和比率数据，探讨了直接使用线性回归和对数变换线性回归处理计数数据的思路及其解释。在下一部分，我们将正式引入泊松回归模型，它为我们处理计数数据提供了更自然、更强大的框架。

039：泊松回归第二部分 📊

在本节课中，我们将深入探讨泊松回归模型，重点比较其与线性回归的区别，学习如何处理计数数据中的速率和比例问题，并了解模型诊断与改进方法。

线性回归 vs. 泊松回归 🔄

上一节我们介绍了泊松回归的基本概念。本节中，我们来看看它与线性回归的核心区别。

在广义线性模型中，我们并非对结果变量本身进行变换，而是对其均值进行变换。

  • 在线性模型中，结果变量等于线性部分加上误差项。其公式可表示为：E(Y) = β₀ + β₁X。
  • 在泊松对数线性模型中，是结果变量期望值的对数等于线性部分。其公式为：log[E(Y)] = β₀ + β₁X。

我们可以通过对等式两边取指数来逆转这个过程，得到：E(Y) = exp(β₀ + β₁X)。

主要区别在于，我们假设数据服从泊松分布，其均值具有 exp(β₀ + β₁X) 的形式。这改变了系数的解释方式，为我们观察到的结果提供了更可信的分布，并且我们得到了相对性的解释（因为一切都在对数尺度上）。

对数变换与系数解释 📈

对结果变量进行对数变换通常是处理正数数据的一个好方法，它能使系数在对数尺度上保持甚至增强可解释性。

当我们把系数转换回原始尺度时，所看到的是乘性差异。

以下是理解系数的方法：
如果我们查看模型 E(Y) = exp(β₀ + β₁ * JulianDate)，那么下一天（JulianDate + 1）的期望均值为 exp[β₀ + β₁ * (JulianDate + 1)]。
用后者除以前者，我们得到 exp(β₁)。

因此，我们的斜率系数 exp(β₁) 可以解释为：预测变量每增加一个单位，结果变量均值相对增加或减少的倍数。在对数尺度上，我们关注系数是否接近零；在指数化后，我们关注其是否接近一。在多变量背景下，exp(β₁) 表示在保持其他预测变量不变的情况下，结果变量的相对变化。

模型拟合与诊断 🔍

将拟合的泊松回归模型叠加到数据上，可以看到它虽然与线性模型拟合得相当接近，但具有我们想要的曲线形态。一个更简单、系数更少的模型似乎能更好地拟合数据。

一个常见的担忧是泊松分布的方差必须等于均值。因此，随着均值上升，方差也应上升。然而，如果我们绘制拟合值与残差的散点图，可以清楚地看到对于较低的均值，方差反而更高。这是一个问题。

以下是处理方差异质性的一些方法：

  • 拟泊松模型：该模型将方差视为均值的常数倍，而非严格等于均值。
  • 模型无关的标准误：例如使用“三明治”方差估计量（源自广义估计方程）。这是一个更高级但非常重要的应用主题。

主要要点是进行残差图分析以判断模型假设是否成立。可以尝试拟泊松模型，这在R中很容易实现。如果问题依然存在（如本例），则需要探索其他解决方案。

在本例中，使用模型无关的置信区间与常规置信区间相比，差异并不大。对于像本例这样的小系数，取指数后结果大约为1.002，意味着每天估计约有0.2%的增长。

处理速率与比例数据 ⚖️

当计数数据有一个“偏移量”时，我们需要处理速率或比例。例如，核泵故障次数应随监控时间变长而增加，流感病例数应随城市人口增多而增加。

我们希望解释的不是结果变量的期望值，而是结果变量除以某个相对项（如监控时间、风险人时等）后的期望值。

假设我们想将比例 E(Y)/T 建模为 exp(β₀ + β₁X)。通过对两边取对数并重新整理，我们得到：log[E(Y)] = β₀ + β₁X + log(T)。

在R中，只需将 log(T) 作为一个没有系数的偏移量项加入线性模型即可。简单的方法是在 glm 函数中使用 offset = log(visits + 1) 参数（加1是为了避免对零取对数）。默认的 family = poisson 即假定为对数连接函数，这正是我们需要的。

通过比较调整偏移量前后的拟合速率，可以直观看到模型如何根据总访问量来调整来自特定来源的访问量比例。

零膨胀与其他问题 ⚠️

泊松数据中另一个常见问题是零值出现得过于频繁，这被称为零膨胀。泊松模型在试图拟合其他计数值时，可能无法很好地处理过多的零值。

在本例早期存在大量零值，且具有时间成分，这使得建模更具挑战性。有专门的R包（如 pscl）可以帮助处理零膨胀泊松数据。

总结 🎯

本节课我们一起学习了泊松回归的核心内容。我们比较了泊松回归与线性回归在模型设定和解释上的根本区别，理解了如何通过对数连接函数处理计数数据。我们探讨了模型诊断的重要性，特别是对方差等于均值这一假设的检验，并介绍了拟泊松模型和稳健标准误等解决方案。最后，我们学习了如何在泊松回归中通过引入偏移量来处理速率或比例数据，并简要提及了零膨胀问题及其处理思路。

希望你能将二元逻辑回归、线性回归及多元回归中的技能，直接应用到泊松回归的实例中。下一讲将是本系列的最后一课，我们将介绍一些使用线性模型的实用技巧。

040：回归模型扩展应用 🎯

在本节课中，我们将探讨回归模型的扩展应用，展示如何利用线性模型拟合复杂函数，并介绍一个有趣的音乐和弦识别案例。这些例子旨在激发大家继续深入学习回归模型和广义线性模型的兴趣。

概述 📋

回归模型是数据分析中极其重要的工具。我们已经学习了如何通过添加多项式项来拟合非线性关系。本节将介绍两种更强大的扩展方法：回归样条和傅里叶基，它们能帮助我们拟合更复杂的函数模式。

使用回归样条拟合复杂函数 📈

上一节我们介绍了通过添加多项式项来拟合曲线。然而，对于更复杂的形状（如S型曲线或正弦波），我们需要更灵活的方法。

回归样条的基本思想

我们的目标是拟合模型：y = f(x) + ε，其中 f(x) 是一个复杂函数。回归样条提供了一种解决方案。

以下是构建回归样条模型的基本步骤：

  1. 定义模型：模型形式为：
    y = β₀ + β₁x + Σ [ (x - κᵢ)₊ * γᵢ ]
    其中，(x - κ)₊ 是一个函数，如果 (x - κ) 为正则返回该值，否则返回0。κᵢ 被称为节点。

  1. 理解节点：可以将节点想象成在一根棍子上选定的“折点”。在这些点处，函数是连续的，但可能不光滑（存在“拐角”）。这些数学上的“折点”允许我们拟合复杂的形状。

  1. 一个简单例子：我们首先生成一些带有噪声的正弦波数据。然后，在数据范围内均匀选择20个节点，并构建包含截距项、线性项和所有节点项的预测矩阵。最后，拟合一个普通的线性模型。结果能够较好地追踪正弦波趋势，但在节点处曲线显得不够光滑。

改进：光滑样条

为了使拟合曲线在节点处也光滑（即一阶导数连续），我们可以对模型进行一个简单的改进：添加平方项。

改进后的模型为：
y = β₀ + β₁x + β₂x² + Σ [ (x - κᵢ)₊² * γᵢ ]

通过将节点项平方，我们确保了函数在节点处不仅连续，而且拥有连续的一阶导数。再次用相同的数据和节点拟合模型，得到的曲线就变得非常光滑了。

关于回归样条的注意事项

  • 基函数：上述节点项的集合构成了一组基函数，它们是构建更复杂函数的“积木”。回归样条只是众多基函数中的一种。
  • 节点选择：一个关键挑战是如何选择节点的数量和位置。节点太少可能导致欠拟合，太多则可能导致过拟合。
  • 现代解决方案：现代方法通常放置大量节点，然后通过添加正则化项来惩罚节点项的系数大小，从而自动控制模型复杂度，这属于更高级的主题。

通过这个例子，我们展示了仅用普通回归就能拟合相当复杂的函数，并且这种方法可以自然地推广到广义线性模型中。

应用案例：使用傅里叶基识别音乐和弦 🎵

接下来，我们看一个更有趣的例子：如何使用线性模型识别一段音频中的和弦。

问题设定

一个和弦由多个音符同时演奏组成。例如，C大调和弦包含C、E、G三个音符。每个音符对应一个特定的频率。如果我们录制了一段C大调和弦的数字化音频，能否让R语言分析出它是由哪几个音符构成的？

建模过程

  1. 生成数据：我们首先获取钢琴上中音C到高音C一个八度内8个音符的标准频率。然后，为C、E、G三个音符分别生成对应频率的正弦波，并将它们相加，合成一段持续2秒的C大调和弦数字信号。
  2. 构建傅里叶基：我们创建一个基，它包含这8个可能音符频率对应的正弦函数。
  3. 拟合线性模型：将和弦信号作为因变量，8个正弦函数作为自变量（去除截距项），拟合一个线性模型。
  4. 解读结果：绘制模型的系数图。结果显示，对应C、E、G三个频率的系数远大于其他频率的系数。因此，模型成功地识别出了构成该和弦的音符。

与离散傅里叶变换的联系

我们的模型只使用了正弦项。实际上，完整的分析还应包括余弦项。离散傅里叶变换 所做的正是这件事：它用所有可能的正弦和余弦项（数量等于采样点数量）来拟合一个完全饱和的线性模型。

在信号处理中，工程师们经常计算信号的离散傅里叶变换并绘制其频谱图，这本质上就是在进行一个包含大量正弦和余弦回归项的线性模型拟合。之所以使用快速傅里叶变换算法而不是直接求解线性模型，是因为前者计算效率极高。

这个案例表明，线性模型的能力可能超乎我们的直觉想象，能够解决诸如音频分析这类看似不相关的问题。

总结与展望 🚀

本节课我们一起学习了回归模型的两种强大扩展。

  • 我们介绍了回归样条，它通过引入节点和特定的基函数，使我们能够用线性模型拟合非常复杂的非线性关系。
  • 我们探讨了傅里叶基的应用，并通过一个识别音乐和弦的生动案例，展示了线性模型在信号处理等领域的实用性。

这两个例子都凸显了回归模型技术的强大与灵活。它们不仅功能丰富，而且实现起来代码量并不大。

后续学习建议

如果你希望在此基础上继续深入，建议从以下两个方向拓展：

  1. 深入学习广义线性模型：本课程仅初步涉及，GLM家族庞大，适用于更多类型的数据。
  2. 学习处理相关数据和纵向数据：本课程假设误差项独立。但在许多实际场景中（如重复测量、家族数据），数据点之间存在相关性。处理这种相关性是线性模型扩展到更广阔领域（如多水平模型、纵向数据分析）的重要部分。

感谢学习本课程！希望你能运用所学知识，并在回归建模的道路上继续探索。

posted @ 2026-03-26 08:55  布客飞龙I  阅读(24)  评论(0)    收藏  举报