数字本无实形，亦无情感。可你若逼得太紧，就连数字也会显露边界。
— 莫斯・戴夫（美国知名说唱歌手、演员，本名 Dante Terrell Smith）
基础概率理论
在实践中，因果推断依赖从极简到极复杂的各类统计模型，而构建这类模型需要基础的概率理论知识，因此我们从核心定义开始。随机过程是指可多次重复、且每次结果都可能不同的过程。样本空间是随机过程所有可能结果的集合。我们将随机过程分为离散型与连续型：离散过程的结果为整数（或可数的有限 / 无限值），而连续过程的结果可以是任意数值（包括小数、分数）。
独立事件有两种定义，第一种是逻辑独立性 ：两个事件同时发生，但没有任何逻辑依据表明它们会互相影响。若错误地认为两者存在因果关联，就会陷入后此谬误（post hocer go propter hoc） —— 拉丁语意为 “在此之后，因此之故”。这种谬误的本质是：事件发生的时间先后顺序，不足以证明前者是后者的原因。
独立事件的第二种定义是统计独立性 。我们将借助有放回抽样与无放回抽样的案例来解释该定义。 以一副充分洗匀的 52 张扑克牌为例：从中抽取第一张牌，恰好是 A 的概率是多少？ 样本空间（即该随机过程所有可能结果的集合）共有 52 种结果。 在这 52 种结果里，我们只关注抽到 A 的情况。整副牌有 4 张 A，因此：
\[\frac{4}{52} = 0.077\]
假设第一张抽出的牌是 A。现在我们再提出同样的问题：在剩余牌中继续抽一张，第二张也恰好是 A 的概率是多少？ 这时概率不再是 \(\dfrac{1}{13}\) ，因为我们没有采用有放回抽样，而是无放回抽样。 因此新的概率为： \[\frac{3}{51} \approx 0.059\] 在无放回抽样条件下，第一张牌为 A与第一张为 A 时第二张也为 A这两个事件，并不是独立事件。 要使两个事件相互独立，就必须把抽出的 A 放回牌堆，重新洗牌。 因此事件 \(A\) 、\(B\) 相互独立的充要条件为： \[P(A\mid B) = P(A)\] 独立事件举例：先用一枚骰子掷出 3 点，再用另一枚骰子掷出 5 点。 这两个事件相互独立，因此无论第一颗骰子掷出几点，第二颗掷出 5 点的概率恒为 \(0.17\) 。 但如果我们想求解需要多个事件依次发生才能实现某一结果的概率，该如何计算？
举个例子：以克利夫兰骑士队争夺 NBA 总冠军为例。 2016 年 NBA 季后赛为七场四胜制，金州勇士大比分 3–1 领先。 勇士若输掉这轮系列赛，需要满足什么条件？骑士必须连胜三场。 这种情况下，计算概率需要将所有边际概率相乘，公式记作： \(\Pr(\cdot)^n\) 其中 \(\Pr(\cdot)\) 为单个事件发生的边际概率，\(n\) 为该事件连续发生的次数。 若骑士每场比赛获胜的无条件概率为 \(0.5\) ，且各场比赛相互独立， 那么骑士从 3–1 落后完成逆转的概率，为每场获胜概率的连乘积： \[\Pr(W,W,W) = 0.5^3 = 0.125\] 再看一个例子便于理解： 在德州扑克中，每位玩家会拿到两张暗牌。当拿到两张同点数牌时，称为手握 “口袋对子”。 那么拿到口袋 A的概率是多少？ \[\frac{4}{52} \times \frac{3}{51} \approx 0.0045 = 0.45%\] 没错，概率仅为 \(0.45%\) 。 我们把上述规律整理为通用形式： 对于独立事件，计算联合概率的规则是各边际概率相乘： \[\Pr(A,B) = \Pr(A)\Pr(B)\] 其中 \(\Pr(A,B)\) 表示事件 \(A\) 与事件 \(B\) 同时发生的联合概率， \(\Pr(A)\) 表示事件 \(A\) 单独发生的边际概率。
接下来我们看一个稍难一点的应用问题： 掷两枚六面骰子，掷出点数和为 7的概率是多少？该概率与掷出点数和为 3的概率是否相等？ 我们通过对比两种情况的概率来解答，借助表格帮你直观理解。
首先，先看两枚六面骰子掷出点数和为 7 的所有组合方式。 掷两枚骰子，总共有 36 种可能结果（\(6^2=36\) ）。 如表 2 所示，仅用两枚骰子凑出点数和为 7，共有 6 种不同组合。 因此掷出点数和为 7 的概率为 \(6/36 = 16.67%\) 。 再来看两枚六面骰子掷出点数和为 3 的所有组合方式。 如表 3 所示，掷两枚骰子凑出点数和为 3，仅有 2 种组合。 因此掷出点数和为 3 的概率为 \(2/36 = 5.56%\) 。 由此可得结论：掷出点数和为 7 与掷出点数和为 3 的概率并不相等。
事件与条件概率
在讲解概率的三种表示方法 之前，我们先引入几个新的术语和概念：事件 与条件概率 。
设事件为A，另一个事件为B。对于任意两个事件，一共存在四种可能的情况：
A 且 B ：事件 A 和事件 B 同时发生；
非 A 且 B ：事件 A 不发生，事件 B 发生；
A 且非 B ：事件 A 发生，事件 B 不发生；
非 A 且非 B ：事件 A、事件 B 均不发生。
概率树 我们来思考一个场景：你正在考取驾照。假设要拿到驾照，你必须通过笔试 和路考 两项考试。但如果笔试不及格，你就没有资格参加路考。我们可以用一棵概率树 来表示这两个事件的关系。
概率树直观易懂、便于解读²。 首先，我们可以看到笔试通过的概率为 0.75，笔试不通过的概率为 0.25。 其次，从任意一个节点分出的所有分支，其对应概率之和为 1.0。所有联合概率的总和也为 1.0。这一规律称为全概率公式，它等于事件 A 与所有\(B_n\) 事件同时发生的联合概率之和： \[\Pr(A) = \sum_n \Pr(A \cap B_n)\] 在驾照的概率树中，我们也能看到条件概率的概念。例如，“笔试已通过的前提下，路考不及格” 的概率，可表示为： \(\Pr({Fail} \mid {Pass}) = 0.45\)
维恩图与集合 表示多个事件同时发生的第二种方式是使用维恩图。维恩图由约翰・维恩于 1880 年首次提出。它们常被用于教授基础集合论，也可用来表示概率论与统计学中的集合关系。下面的示例将包含两个集合：集合 A 和集合 B。
德克萨斯大学的橄榄球主教练整个赛季都与体育主管和校董会关系紧张。经历了几个平庸的赛季后，他在学校的未来岌岌可危。如果长角牛队没能打进一场大型碗赛，他很可能不会被续约；但如果打进了，他大概率会留任。我们就以这位教练的处境为例，来探讨基础集合论。
在此之前，我们先回顾一下相关术语：A 和 B 为事件，U 是包含 A、B 两个子集的全集。设事件 A 为“长角牛队获得大型碗赛的邀请”，事件 B 为“教练获得续约”。已知 \(P(A)=0.6\) ，\(P(B)=0.8\) ，且 A、B 同时发生的概率 \(P(A,B)=0.5\) 。
注意：\(A \cup \sim{A} = U\) ，其中 \(\sim{A}\) 是事件 A 的补集，即全集中所有不属于 A 的部分。同理，事件 B 的补集 \(\sim B\) 也满足 \(B \cup \sim B = U\) 。因此： \[A \cup \sim A = B \cup \sim B\]
我们可以据此改写以下定义： \[A = B \cup \sim B - \sim A\] \[B = A \cup \sim A - \sim B\]
当我们想描述“事件A或事件B可能发生”的事件集合时，记作：\(A \cup B\) ，读作“A并B”，表示这个新集合包含A中的所有元素与B中的所有元素。只要元素属于集合A或集合B中的任意一个，就属于这个新的并集。 当我们想描述“事件A和事件B同时发生”的事件集合（即交集）时，记作：\(A \cap B\) ，读作“A交B”。这个新集合只包含同时属于A和B的元素，也就是说，只有同时在A和B中的元素才会被加入这个新集合。
现在我们来仔细看一个与集合A相关的关系： \[A = A \cup B + A \cup \sim B\] 这句话的意思是：集合A可以分为两部分。第一部分是A与B同时发生的所有情况；第二部分是A发生但B不发生的情况，也就是式子中的\(A \cup \sim B\) ，这部分覆盖了集合A中剩下的所有元素。
类似的推理方式可以帮助你理解下面的表达式： \[A \cap B = A \cup \sim B + \sim A \cup B + A \cup B\] 要得到A交B的集合，我们需要三个部分：A发生但B不发生的集合、B发生但A不发生的集合，以及它们的交集。把这三部分加起来，就得到了\(A \cap B\) 。
接下来只需简单计算就能求出所有缺失的值。回顾一下：A代表球队打进季后赛，\(P(A)=0.6\) ；B代表教练获得续约，\(P(B)=0.8\) ；同时，\(P(A,B)=0.5\) ，即A和B同时发生的概率。于是我们可以得到：
\[A = A \cap B + A \cap \sim B\]
\[A \cap \sim B = A - A \cap B\]
\[P(A, \sim B) = P(A) - P(A, B)\]
\[P(A, \sim B) = 0.6 - 0.5\]
\[P(A, \sim B) = 0.1\]
在处理集合问题时，要理解概率是通过计算子集（例如 \(A \cap B\) ）在集合（例如 \(A\) ）中所占的比例得到的，这一点非常重要。当我们计算事件 \(A \cap B\) 发生的概率时，它是相对于全集 \(U\) 而言的。但如果我们问：“在事件 \(A\) 发生的前提下，\(A \cap B\) 所占的比例是多少？”，那么我们就需要进行如下计算：
\[? = A \cap B \div A\]
\[? = 0.5 \div 0.6\]
\[? \approx 0.83\]
我特意没有在左侧写明定义，是为了先聚焦计算本身。现在我们来明确一下计算的含义：在事件A已经发生的前提下，事件B也发生的概率是多少？这就是：
\[P(B|A) = \frac{P(A,B)}{P(A)} = \frac{0.5}{0.6} \approx 0.83\]
\[P(A|B) = \frac{P(A,B)}{P(B)} = \frac{0.5}{0.8} \approx 0.63\]
注意，这些条件概率在维恩图中并不直观。我们本质上是在问：某个子集（如 \(P(A)\) ）中，由交集（如 \(P(A,B)\) ）贡献的部分占比是多少。这种推导方式，正是条件概率概念的定义方式。
蒙提霍尔问题
我们来看另一个例子——蒙提霍尔问题。这个问题很有趣，因为大多数人都觉得它反直觉，甚至曾难倒不少数学家和统计学家。但贝叶斯定理能把答案解释得非常清楚，也正因如此，贝叶斯定理本身也曾引发过争议（McGrayne, 2012）。
假设有三扇关闭的门：门1（\(D_1\) ）、门2（\(D_2\) ）和门3（\(D_3\) ）。其中一扇门后藏着一百万美元奖金，另外两扇门后各有一只山羊。主持人蒙提·霍尔会让参赛者先选一扇门。参赛者选好门后，在他打开这扇门之前，蒙提会打开另外两扇门中的一扇，露出一只山羊，然后问参赛者：“你要不要换另一扇门？”
面对蒙提的提议，人们的普遍反应是觉得没必要换门，因为奖金在剩下两扇门后的概率是均等的，都是50%。但直觉告诉我们，这显然不是正确答案——蒙提打开那扇门的行为，其实已经传递了信息。
我们用概率符号来形式化描述这个问题。假设你最初选择了门1（\(D_1\) ）。此时，门1后藏有奖金的概率是\(P(D_1 = 奖金) = 1/3\) ，我们称之为事件\(A_1\) 。游戏开始时，每扇门后有奖金的概率都是1/3，因此\(P(A_1) = 1/3\) 。根据全概率公式，门1后没有奖金的概率\(P(\sim A_1) = 2/3\) 。接着，假设蒙提打开了门2（\(D_2\) ），露出一只山羊，然后问你：“要不要换成门3？”
我们需要计算门3后藏有奖金的概率，并将门1的概率进行对比。我们把“蒙提打开门2”称为事件\(B\) ，把“奖金在第\(i\) 扇门后”称为事件\(A_i\) 。我们可以用贝叶斯分解来解答这个问题，最终要计算的是：在蒙提打开了门2（事件\(B\) ）的条件下，门1后藏有奖金（事件\(A_1\) ）的概率，这是一个条件概率问题。根据公式2.11，我们可以写出：
\[P(A_1|B) = \frac{P(B|A_1)P(A_1)}{P(B|A_1)P(A_1) + P(B|A_2)P(A_2) + P(B|A_3)P(A_3)} \tag{2.12}\]
公式右侧主要有两类概率：一是奖金在某扇门后的边缘概率\(P(A_i)\) ；二是在奖金在门\(i\) 后的条件下，蒙提打开门2的条件概率\(P(B|A_i)\) 。
在没有额外信息的情况下，奖金在门\(i\) 后的先验概率（也叫无条件概率）是\(1/3\) 。
条件概率\(P(B|A_i)\) 则需要仔细分析： 1. \(P(B|A_1)\) ：如果奖金在门1后，蒙提打开门2的概率是多少？此时他可以打开门2或门3，因此\(P(B|A_1) = 0.5\) 。 2. \(P(B|A_2)\) ：如果奖金在门2后，蒙提打开门2的概率是多少？蒙提不会打开有奖金的门，因此这个概率是\(0\) 。 3. \(P(B|A_3)\) ：如果奖金在门3后，蒙提打开门2的概率是多少？他不能打开你选的门1，也不能打开有奖金的门3，因此只能打开门2，概率为\(1\) 。
将这些值代入公式2.12： \[P(A_1|B) = \frac{\frac{1}{2} \cdot \frac{1}{3}}{\frac{1}{2} \cdot \frac{1}{3} + 0 \cdot \frac{1}{3} + 1.0 \cdot \frac{1}{3}} = \frac{\frac{1}{6}}{\frac{1}{6} + \frac{2}{6}} = \frac{1}{3}\]
结果显示，你最初选择的门1后有奖金的概率仍然是\(1/3\) 。那门3的概率呢？我们用同样的方法计算：
\[P(A_3|B) = \frac{P(B|A_3)P(A_3)}{P(B|A_3)P(A_3) + P(B|A_2)P(A_2) + P(B|A_1)P(A_1)} = \frac{1.0 \cdot \frac{1}{3}}{1.0 \cdot \frac{1}{3} + 0 \cdot \frac{1}{3} + \frac{1}{2} \cdot \frac{1}{3}} = \frac{2}{3}\]
有趣的是，你对最初选择的门的概率信念没有改变，但对另一扇门的信念却发生了变化。门3的先验概率\(P(A_3)=1/3\) ，经过更新后变成了后验概率\(P(A_3|B)=2/3\) 。这意味着，蒙提打开门2的行为，让你获得了新信息，从而修正了对奖金位置的判断。
正如脚注中提到的，即使是聪明人，基于贝叶斯定理的推导也常常令人惊讶——这或许是因为我们缺乏将信息整合到概率中的系统性方法。贝叶斯定理为我们提供了一种逻辑、准确的方式来完成这一点。此外，贝叶斯定理还开启了一种从“结果”推断“原因”的全新思维方式。本书大部分内容都是从已知原因推断结果，而贝叶斯定理提醒我们，也可以从已知结果出发，形成对原因的合理推断。
求和符号
用于因果推断的工具建立在概率论的基础之上。我们会经常用到统计学中的数学工具与概念，比如期望和概率。线性回归模型是本书中最常用的工具之一，但在深入学习之前，我们需要先掌握一些基础符号。我们从求和符号开始：希腊字母大写Σ代表求和运算符。设 \(x_1, x_2, ..., x_n\) 为一组数，我们可以用求和符号简洁地表示它们的和：
\[\sum_{i=1}^n x_i \equiv x_1 + x_2 + ... + x_n\]
字母 \(i\) 被称为求和下标 ，也可以用 \(j\) 或 \(k\) 等其他字母作为下标。下标变量代表随机变量 \(x\) 的一个特定取值。数字 \(1\) 和 \(n\) 分别是求和的下限 和上限 。表达式 \(\sum_{i=1}^n x_i\) 的含义是：“对 \(i\) 从 \(1\) 到 \(n\) 对应的所有 \(x_i\) 求和”。下面的例子可以帮助理解：
\[\sum_{i=6}^9 x_i = x_6 + x_7 + x_8 + x_9\]
求和符号有三个基本性质：
常数规则 对于任意常数 \(c\) ： \[\sum_{i=1}^n c = nc \tag{2.13}\] 例如： \[\sum_{i=1}^3 5 = (5 + 5 + 5) = 3 \cdot 5 = 15\]
数乘规则 \[\sum_{i=1}^n cx_i = c\sum_{i=1}^n x_i \tag{2.14}\] 例如： \[\begin{align*}
\sum_{i=1}^3 5x_i &= 5x_1 + 5x_2 + 5x_3 \\
&= 5(x_1 + x_2 + x_3) \\
&= 5\sum_{i=1}^3 x_i
\end{align*}\]
线性性质 结合前两个性质，可得：对于任意常数 \(a\) 和 \(b\) ： \[\sum_{i=1}^n (ax_i + by_i) = a\sum_{i=1}^n x_i + b\sum_{i=1}^n y_i\]
在结束对求和符号的介绍前，我们还需了解几个常见误区 ：
比值的和不等于和的比值： \[\sum_i \frac{x_i}{y_i} \neq \frac{\sum_{i=1}^n x_i}{\sum_{i=1}^n y_i}\]
变量平方的和不等于和的平方： \[\sum_{i=1}^n x_i^2 \neq \left( \sum_{i=1}^n x_i \right)^2\]
我们可以用求和符号进行多种计算，其中一些会在本书中反复出现。例如，计算平均值：
\[\begin{align*}
\bar{x} &= \frac{1}{n}\sum_{i=1}^n x_i \\
&= \frac{x_1 + x_2 + \dots + x_n}{n}
\end{align*} \tag{2.15}\]
其中 \(\bar{x}\) 是随机变量 \(x_i\) 的平均值（均值）。另一个常见计算是随机变量对其均值的离差，离差的和恒为 \(0\) 。
\[\sum_{i=1}^n (x_i - \bar{x}) = 0 \tag{2.16}\]
你可以在\(表5\) 中看到这一点。
考虑两个序列 \(\{y_1, y_2, ..., y_n\}\) 和 \(\{x_1, x_2, ..., x_n\}\) 。我们可以对 \(x\) 和 \(y\) 的所有可能取值进行双重求和。例如，当 \(n=m=2\) 时，\(\sum_{i=1}^2 \sum_{j=1}^2 x_i y_j\) 等于 \(x_1y_1 + x_1y_2 + x_2y_1 + x_2y_2\) 。推导过程如下： \[\begin{align*}
x_1y_1 + x_1y_2 + x_2y_1 + x_2y_2 &= x_1(y_1 + y_2) + x_2(y_1 + y_2) \\
&= \sum_{i=1}^2 x_i(y_1 + y_2) \\
&= \sum_{i=1}^2 x_i \left( \sum_{j=1}^2 y_j \right) \\
&= \sum_{i=1}^2 \left( \sum_{j=1}^2 x_i y_j \right) \\
&= \sum_{i=1}^2 \sum_{j=1}^2 x_i y_j
\end{align*}\]
本书中一个非常有用的结论是： \[\sum_{i=1}^n (x_i - \bar{x})^2 = \sum_{i=1}^n x_i^2 - n(\bar{x})^2 \tag{2.17}\]
以下是详细的分步证明（为便于阅读，后续步骤省略了求和下标）： \[\begin{align*}
\sum_{i=1}^n (x_i - \bar{x})^2 &= \sum_{i=1}^n (x_i^2 - 2x_i\bar{x} + \bar{x}^2) \\
&= \sum x_i^2 - 2\bar{x}\sum x_i + n\bar{x}^2 \\
&= \sum x_i^2 - 2\frac{1}{n}\sum x_i \sum x_i + n\bar{x}^2 \\
&= \sum x_i^2 + n\bar{x}^2 - \frac{2}{n}\left( \sum x_i \right)^2 \\
&= \sum x_i^2 + n\left( \frac{1}{n}\sum x_i \right)^2 - 2n\left( \frac{1}{n}\sum x_i \right)^2 \\
&= \sum x_i^2 - n\left( \frac{1}{n}\sum x_i \right)^2 \\
&= \sum x_i^2 - n\bar{x}^2
\end{align*}\]
该结论的更一般形式为： \[\begin{align*}
\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) &= \sum_{i=1}^n x_i(y_i - \bar{y}) \\
&= \sum_{i=1}^n (x_i - \bar{x})y_i \\
&= \sum_{i=1}^n x_i y_i - n(\bar{x}\bar{y})
\end{align*} \tag{2.18}\]
或者写成： \[\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) = \sum_{i=1}^n x_i(y_i - \bar{y}) = \sum_{i=1}^n (x_i - \bar{x})y_i = \sum_{i=1}^n x_i y_i - n(\bar{x}\bar{y}) \tag{2.19}\]
期望值
随机变量的期望值（也称为期望，有时也称总体均值），是该变量所有可能取值的加权平均值，权重为总体中各取值发生的概率。假设随机变量\(X\) 的取值为\(x_1,x_2,\dots,x_k\) ，对应的概率分别为\(f(x_1),f(x_2),\dots,f(x_k)\) ，则其期望值定义为： \[\begin{align*}
E(X) &= x_1f(x_1)+x_2f(x_2)+\dots+x_kf(x_k) \\
&= \sum_{j=1}^k x_j f(x_j)
\end{align*} \tag{2.20}\]
我们来看一个数值例子：若\(X\) 的取值为\(-1,0,2\) ，对应概率分别为\(0.3,0.3,0.4\) ，则其期望值为： \[\begin{align*}
E(X) &= (-1)(0.3)+(0)(0.3)+(2)(0.4) \\
&= 0.5
\end{align*}\]
我们也可以计算该变量函数的期望，例如\(X^2\) 。\(X^2\) 的取值为\(1,0,4\) ，对应概率仍为\(0.3,0.3,0.4\) ，因此： \[\begin{align*}
E(X^2) &= (-1)^2(0.3)+(0)^2(0.3)+(2)^2(0.4) \\
&= 1.9
\end{align*}\]
期望值的基本性质：
对任意常数\(c\) ，有\(E(c)=c\) ；
对任意常数\(a,b\) ，有\(E(aX+b)=aE(X)+b\) ；
对常数\(a_1,\dots,a_n\) 和随机变量\(X_1,\dots,X_n\) ，有： \[ E(a_1X_1+\dots+a_nX_n)=a_1E(X_1)+\dots+a_nE(X_n)
\] 用求和符号表示为： \[ E\left(\sum_{i=1}^n a_iX_i\right) = \sum_{i=1}^n a_iE(X_i)
\] 当\(a_i=1\) 时，有： \[ E\left(\sum_{i=1}^n X_i\right) = \sum_{i=1}^n E(X_i)
\]
此外，对任意两个随机变量\(W\) 和\(H\) ，期望算子还有以下性质：
\(E(aW+b)=aE(W)+b\) （对任意常数\(a,b\) ）
\(E(W+H)=E(W)+E(H)\)
\(E(W-E(W))=0\)
方差
方差的定义：随机变量\(W\) 的总体方差为： \[V(W)=\sigma^2=E\left[(W-E(W))^2\right]\] 可推导得： \[V(W)=E(W^2)-E(W)^2 \tag{2.21}\]
在样本数据中，方差的估计式为： \[\widehat{S}^2=(n-1)^{-1}\sum_{i=1}^n(x_i-\bar{x})^2\] 除以\(n-1\) 是为了对均值估计进行自由度校正；在大样本中，该校正对\(S^2\) 的实际影响可忽略。
方差的性质： 1. 对任意常数\(a,b\) ，有\(V(aX+b)=a^2V(X)\) ； 2. 常数的方差为\(0\) ，即对任意常数\(c\) ，\(V(c)=0\) ； 3. 两个随机变量和的方差： \[ V(X+Y)=V(X)+V(Y)+2\left(E(XY)-E(X)E(Y)\right) \tag{2.22}
\] 若\(X\) 与\(Y\) 相互独立，则\(E(XY)=E(X)E(Y)\) ，此时\(V(X+Y)=V(X)+V(Y)\) 。
协方差
方程2.22中的最后一项被称为协方差 。协方差衡量两个随机变量之间的线性相关程度，我们用算子\(C(X,Y)\) 表示： - \(C(X,Y) > 0\) 表示两个变量同向变化 - \(C(X,Y) < 0\) 表示两个变量反向变化
因此，方程2.22可改写为： \[V(X + Y) = V(X) + V(Y) + 2C(X,Y)\]
需要注意：协方差为0，并不代表两个变量之间没有关系——它们可能存在非线性关系。协方差的定义为： \[C(X,Y) = E(XY) - E(X)E(Y) \tag{2.23}\]
如前所述，若\(X\) 与\(Y\) 相互独立，则总体中\(C(X,Y)=0\) 。两个线性函数的协方差满足： \[C(a_1 + b_1X, a_2 + b_2Y) = b_1b_2C(X,Y)\] 其中常数\(a_1\) 和\(a_2\) 的影响会被抵消，因为常数的均值为其本身，差值为0。
协方差的数值大小难以直接解读，因此我们引入相关系数 ： 令\(W = \frac{X-E(X)}{\sqrt{V(X)}}\) ，\(Z = \frac{Y-E(Y)}{\sqrt{V(Y)}}\) ，则相关系数定义为： \[\text{Corr}(W,Z) = \frac{C(X,Y)}{\sqrt{V(X)V(Y)}} \tag{2.24}\]
相关系数的取值范围为\([-1,1]\) ： - 正相关表示变量同向变化，负相关表示反向变化 - 系数越接近1或-1，线性关系越强
总体模型
我们从横截面分析开始。假设我们能从感兴趣的总体中抽取一个随机样本，存在两个变量\(x\) 和\(y\) ，我们想了解\(y\) 如何随\(x\) 的变化而变化。
这里会立刻出现三个问题： 1. 如果\(y\) 还受\(x\) 以外的其他因素影响，我们该如何处理？ 2. 连接这两个变量的函数形式是什么？ 3. 如果我们关心\(x\) 对\(y\) 的因果效应，如何将它与单纯的相关性区分开？
我们先从一个具体模型开始： \[y = \beta_0 + \beta_1 x + u \tag{2.25}\]
该模型被假定在总体中成立，方程2.25定义了一个线性双变量回归模型 。在关注因果效应的模型中，方程左侧的项通常被视为“结果”，右侧的项被视为“原因”。
方程2.25通过引入一个称为误差项 的随机变量\(u\) ，明确允许其他因素影响\(y\) 。它还假设\(y\) 与\(x\) 呈线性关系，以此明确了模型的函数形式。 - 系数\(\beta_0\) 被称为截距参数 - 系数\(\beta_1\) 被称为斜率参数
这些参数描述了总体特征，我们在实证工作中的目标是估计它们的值。我们永远无法直接观测到这些参数，因为它们不是数据（本书会反复强调这一点），但我们可以利用数据和假设来估计它们。要做到这一点，我们需要可靠的假设，以便用数据准确估计这些参数。在这个简单的回归框架中，所有决定\(y\) 的未观测变量都被包含在误差项\(u\) 中。
首先，我们做一个不影响一般性的简化假设：令总体中\(u\) 的期望值为0，即： \[E(u) = 0 \tag{2.26}\] 其中\(E(\cdot)\) 是前面讨论过的期望算子。将\(u\) 标准化为0不会产生任何影响，因为截距项\(\beta_0\) 总能提供这种灵活性。如果\(u\) 的平均值不为0（例如为\(\alpha_0\) ），我们可以调整截距项，将模型重写为： \[y = (\beta_0 + \alpha_0) + \beta_1 x + (u - \alpha_0)\] 其中新的误差项为\(u-\alpha_0\) ，新的截距项为\(\beta_0 + \alpha_0\) 。调整截距项并不会改变斜率参数 \(\beta_1\) 。
均值独立性
均值独立性是一个与统计学基础内容相契合的假设，它针对总体中由 \(x\) 的取值所划分的每个“切片”，定义了误差项的均值： \[E(u|x) = E(u) \quad \text{对所有 } x \text{ 成立} \tag{2.27}\] 其中 \(E(u|x)\) 表示“给定 \(x\) 时 \(u\) 的期望值”。如果方程 2.27 成立，我们就说 \(u\) 与 \(x\) 均值独立。
举个例子：假设我们要估计教育对工资的影响，其中 \(u\) 代表未观测到的能力。均值独立性要求： \[E(\text{ability} \mid x=8) = E(\text{ability} \mid x=12) = E(\text{ability} \mid x=16)\] 也就是说，八年级教育、十二年级教育和大学教育对应的总体子群体中，平均能力水平是相同的。但在现实中，人们会根据自身未观测到的技能和特质来决定教育投入，因此这个假设很可能不成立。
不过，如果我们愿意接受这一假设，将 \(E(u|x)=E(u)\) （均值独立性）与 \(E(u)=0\) （标准化假设）结合，就能得到新的假设： \[E(u|x) = 0 \quad \text{对所有 } x \text{ 成立} \tag{2.28}\]
方程 2.28 被称为零条件均值假设 ，是回归模型中的关键识别假设。由于条件期望是线性算子，\(E(u|x)=0\) 意味着： \[E(y|x) = \beta_0 + \beta_1 x\] 这表明总体回归函数是 \(x\) 的线性函数，也就是 Angrist 和 Pischke（2009）所说的条件期望函数。这一关系是理解参数 \(\beta_1\) 作为因果参数的关键。
普通最小二乘法（OLS）
给定 \(x\) 和 \(y\) 的数据，如何估计总体参数 \(\beta_0\) 和 \(\beta_1\) ？设 \((x_i, y_i)\) （\(i=1,2,\dots,n\) ）是来自总体的随机样本，将任意观测值代入总体方程： \[y_i = \beta_0 + \beta_1 x_i + u_i\] 其中 \(i\) 表示特定观测值，我们能观测到 \(y_i\) 和 \(x_i\) ，但无法观测到 \(u_i\) 。利用前面讨论的两个总体约束： \[\begin{align*}
E(u) &= 0 \\
E(u|x) &= 0
\end{align*}\] 可以得到 \(\beta_0\) 和 \(\beta_1\) 的估计方程。第二个约束意味着误差项的均值不随 \(x\) 的不同而变化，这种独立性假设意味着 \(E(xu)=0\) 、\(E(u)=0\) 且 \(C(x,u)=0\) 。\(C(x,u)=0\) 表明 \(x\) 和 \(u\) 不相关。将 \(u = y - \beta_0 - \beta_1 x\) 代入，得到： \[\begin{align*}
E(y - \beta_0 - \beta_1 x) &= 0 \\
E\left(x(y - \beta_0 - \beta_1 x)\right) &= 0
\end{align*}\] 这两个总体条件决定了 \(\beta_0\) 和 \(\beta_1\) ，而样本对应的条件为： \[\frac{1}{n}\sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = 0 \tag{2.29}\]
\[\frac{1}{n}\sum_{i=1}^n \left(x_i(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)\right) = 0 \tag{2.30}\] 其中 \(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) 是从数据中得到的估计值，这是两个关于未知参数 \(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) 的线性方程。
利用求和算子的性质，从方程 2.29 推导： \[\begin{align*}
\frac{1}{n}\sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) &= \frac{1}{n}\sum_{i=1}^n y_i - \hat{\beta}_0 - \hat{\beta}_1 \left(\frac{1}{n}\sum_{i=1}^n x_i\right) \\
&= \bar{y} - \hat{\beta}_0 - \hat{\beta}_1 \bar{x}
\end{align*}\] 由于方程 2.29 等于 0，因此 \(\bar{y} = \hat{\beta}_0 + \hat{\beta}_1 \bar{x}\) ，可得截距的表达式： \[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\]
将 \(\hat{\beta}_0\) 代入第二个方程，整理可得： \[\sum_{i=1}^n x_i(y_i - \bar{y}) = \hat{\beta}_1 \left[\sum_{i=1}^n x_i(x_i - \bar{x})\right]\]
\[\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) = \hat{\beta}_1 \left[\sum_{i=1}^n (x_i - \bar{x})^2\right]\] 若 \(\sum_{i=1}^n (x_i - \bar{x})^2 \neq 0\) ，则： \[\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} = \frac{\text{样本协方差}(x_i,y_i)}{\text{样本方差}(x_i)} \tag{2.31}\]
\(\hat{\beta}_1\) 被称为普通最小二乘（OLS）斜率估计量，当 \(x_i\) 的样本方差不为 0 时可计算。\(x\) 的变异是识别其对 \(y\) 影响的关键，若样本中所有人的教育年限都相同，则无法估计斜率。
计算出 \(\hat{\beta}_1\) 后，可计算截距估计量： \[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\] 这是 OLS 截距估计量，它是样本均值的线性组合。
对于任意候选估计值 \(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) ，每个 \(i\) 的拟合值为： \[\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i\] 残差定义为： \[\hat{u}_i = y_i - \hat{y}_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i\] 注意残差与误差项的区别：残差是基于拟合值与实际值的预测误差，可通过数据计算；而误差项 \(u\) 是未观测到的总体扰动项，永远不会出现在数据集中。
残差平方和定义为： \[\sum_{i=1}^n \hat{u}_i^2 = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)^2\] OLS 的核心思想是通过选择 \(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) 最小化残差平方和，微积分可证明其解与矩估计结果一致。
最终，OLS 回归方程为： \[\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x \tag{2.32}\]
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 4.0.0 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
set.seed (1 )<- tibble (x = rnorm (10000 ),u = rnorm (10000 ),y = 5.5 * x + 12 * u)
# A tibble: 10,000 × 3
x u y
<dbl> <dbl> <dbl>
1 -0.626 -0.804 -13.1
2 0.184 -1.06 -11.7
3 -0.836 -1.04 -17.0
4 1.60 -1.19 -5.45
5 0.330 -0.500 -4.19
6 -0.820 -0.525 -10.8
7 0.487 -0.302 -0.948
8 0.738 0.472 9.72
9 0.576 -0.248 0.186
10 -0.305 1.26 13.4
# ℹ 9,990 more rows
<- tb %>% lm (y ~ x, .) %>% print ()
Call:
lm(formula = y ~ x, data = .)
Coefficients:
(Intercept) x
-0.04991 5.55690
(Intercept) x
-0.04990882 5.55690164
<- tb %>% mutate (yhat1 = predict (lm (y ~ x, .)),yhat2 = 0.0732608 + 5.685033 * x, uhat1 = residuals (lm (y ~ x, .)),uhat2 = y - yhat2summary (tb[- 1 :- 3 ])
yhat1 yhat2 uhat1 uhat2
Min. :-20.45096 Min. :-20.7982 Min. :-51.5275 Min. :-51.5247
1st Qu.: -3.79189 1st Qu.: -3.7550 1st Qu.: -8.1520 1st Qu.: -8.2751
Median : -0.13842 Median : -0.0173 Median : -0.1727 Median : -0.3147
Mean : -0.08624 Mean : 0.0361 Mean : 0.0000 Mean : -0.1223
3rd Qu.: 3.71578 3rd Qu.: 3.9258 3rd Qu.: 7.9778 3rd Qu.: 7.8506
Max. : 21.12342 Max. : 21.7348 Max. : 44.7176 Max. : 44.4416
%>% lm (y ~ x, .) %>% ggplot (aes (x= x, y= y)) + ggtitle ("OLS Regression Line" ) + geom_point (size = 0.05 , color = "black" , alpha = 0.5 ) + geom_smooth (method = lm, color = "black" ) + annotate ("text" , x = - 1.5 , y = 30 , color = "red" , label = paste ("Intercept = " , - 0.0732608 )) + annotate ("text" , x = 1.5 , y = - 30 , color = "blue" , label = paste ("Slope =" , 5.685033 ))
Warning: `fortify(<lm>)` was deprecated in ggplot2 3.6.0.
ℹ Please use `broom::augment(<lm>)` instead.
ℹ The deprecated feature was likely used in the ggplot2 package.
Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
`geom_smooth()` using formula = 'y ~ x'
我们来看一下回归的输出结果： 首先，如果你对数据进行汇总，会发现拟合值可以通过两种方式得到：一种是使用 Stata 的 `predict` 命令，另一种是手动使用 `generate` 命令计算。为了让读者更好地理解，我两种方法都演示了一遍。 其次，我们来看一下数据，图3中展示了估计出的系数：截距项和x的斜率。可以看到，这两个估计系数都非常接近数据生成过程中预设的真实参数值。
一旦得到估计系数和OLS回归线，我们就可以对任意（合理的）\(x\) 值预测结果变量\(y\) 的值。将特定的\(x\) 值代入方程，就能立即计算出\(y\) 的预测值，预测会存在一定误差。OLS方法的优势在于它能最小化线性函数的预测误差，对于所有线性估计器而言，它是预测\(y\) 的最优估计，因为它将预测误差降到了最低。也就是说，任何估计方法都会存在预测误差，但OLS是“误差最小的那一个”。
需要注意，截距项是当\(x=0\) 时\(y\) 的预测值，在这个样本中，截距项为-0.0750109。斜率项则描述了\(x\) 的合理变化对\(y\) 的预测影响： \[\Delta\hat{y} = \hat{\beta}_1 \Delta x\] 如果\(\Delta x=1\) （即\(x\) 增加1个单位），在这个数值例子中，\(\hat{y}\) 的变化量等于斜率系数5.598296。
计算出\(\hat{\beta}_0\) 和\(\hat{\beta}_1\) 后，将每个\(x_i\) 代入方程，即可得到OLS拟合值（\(i=1,\dots,n\) ）： \[\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i\]
OLS残差的计算公式为： \[\hat{u}_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i\]
大部分残差不为0（即数据点不在回归线上），这一点可以从图3中看出：残差有正有负。正残差表示回归线（预测值）低估了\(y_i\) 的真实值；负残差则表示回归线高估了真实值。
拟合值记为\(\hat{y}_i\) ，残差记为\(\hat{u}_i = y_i - \hat{y}_i\) 。残差与拟合值的散点图呈球状分布，表明两者不相关（图4）。这是OLS方法的固有性质——最小二乘法生成的残差与拟合值不相关，这是由算法本身决定的，并非什么特殊的“魔法”。
OLS的代数性质
回忆一下我们是如何得到 \(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) 的：当模型包含截距项时，有： \[\sum_{i=1}^n \left(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i\right) = 0\]
根据定义，OLS残差的和恒为0： \[\sum_{i=1}^n \hat{u}_i = 0 \tag{2.33}\]
耳听为虚，眼见为实。让我们一起验证这一点：请将以下命令原封不动地输入Stata中运行。
library (tidyverse)set.seed (1 )<- tibble (x = 9 * rnorm (10 ),u = 36 * rnorm (10 ),y = 3 + 2 * x + u,yhat = predict (lm (y ~ x)),uhat = residuals (lm (y ~ x))
x u y yhat
Min. :-7.521 Min. :-79.729 Min. :-48.014 Min. :13.49
1st Qu.:-4.916 1st Qu.: -1.359 1st Qu.: -7.242 1st Qu.:14.09
Median : 2.309 Median : 17.707 Median : 19.612 Median :14.27
Mean : 1.190 Mean : 8.958 Mean : 14.338 Mean :14.34
3rd Qu.: 4.983 3rd Qu.: 32.875 3rd Qu.: 45.343 3rd Qu.:14.73
Max. :14.358 Max. : 54.424 Max. : 50.268 Max. :14.90
uhat
Min. :-61.502
1st Qu.:-21.943
Median : 5.161
Mean : 0.000
3rd Qu.: 30.739
Max. : 36.282
x u y yhat uhat
1.189825e+01 8.958419e+01 1.433807e+02 1.433807e+02 1.776357e-15
输出结果可以总结如下表（表6）所示。
注意误差项 \(u\) 、拟合值 \(\hat{y}\) 和残差 \(\hat{u}\) 三列的区别。将这十行数据求和，会发现误差项和 \(y\) 的拟合值的和都不为0，但残差的和为0。正如我们所说，这是OLS的代数性质之一——估计系数的选择本身就保证了残差之和为0。
根据定义，\(y_i = \hat{y}_i + \hat{u}_i\) （这一点也可以从表6中看出），我们对等式两边取样本均值： \[\frac{1}{n}\sum_{i=1}^n y_i = \frac{1}{n}\sum_{i=1}^n \hat{y}_i + \frac{1}{n}\sum_{i=1}^n \hat{u}_i\] 由于残差之和为0，因此 \(\bar{y} = \overline{\hat{y}}\) 。同理，我们在估计过程中也得到： \[\sum_{i=1}^n x_i\left(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i\right) = 0\]
解释变量与残差之间的样本协方差（因此样本相关系数也）始终为0（参见表6）。
\[\sum_{i=1}^n x_i \hat{u}_i = 0\] 由于拟合值 \(\hat{y}_i\) 是 \(x_i\) 的线性函数，因此拟合值与残差也不相关（参见表6）： \[\sum_{i=1}^n \hat{y}_i \hat{u}_i = 0\] 这两个性质是由OLS的定义自动成立的，换句话说，\(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) 的选择本身就保证了这一点。
第三个性质是：当 \(x\) 取其样本均值 \(\bar{x}\) 时，模型预测的 \(y\) 值也等于其样本均值 \(\bar{y}\) ，即点 \((\bar{x},\bar{y})\) 一定在OLS回归线上： \[\bar{y} = \hat{\beta}_0 + \hat{\beta}_1 \bar{x}\]
拟合优度
对每个观测值，我们都有： \[y_i = \hat{y}_i + \hat{u}_i\]
定义总平方和（SST）、解释平方和（SSE）与残差平方和（SSR）： \[\text{SST} = \sum_{i=1}^n (y_i - \bar{y})^2 \tag{2.34}\]
\[\text{SSE} = \sum_{i=1}^n (\hat{y}_i - \bar{y})^2 \tag{2.35}\]
\[\text{SSR} = \sum_{i=1}^n \hat{u}_i^2 \tag{2.36}\] 将它们除以 \(n-1\) 后，分别对应 \(y_i\) 、\(\hat{y}_i\) 和 \(\hat{u}_i\) 的样本方差。对SST做如下变换： \[\begin{align*}
\text{SST} &= \sum_{i=1}^n (y_i - \bar{y})^2 \\
&= \sum_{i=1}^n \left[(y_i - \hat{y}_i) - (\hat{y}_i - \bar{y})\right]^2 \\
&= \sum_{i=1}^n \left[\hat{u}_i - (\hat{y}_i - \bar{y})\right]^2
\end{align*}\] 由于拟合值与残差不相关，因此有： \[\text{SST} = \text{SSE} + \text{SSR}\]
假设 \(\text{SST} > 0\) ，我们定义 \(y_i\) 的总变异中被解释变量 \(x_i\) （即OLS回归线）解释的比例为： \[R^2 = \frac{\text{SSE}}{\text{SST}} = 1 - \frac{\text{SSR}}{\text{SST}}\] 这就是回归的\(R^2\) （决定系数），它也等于 \(y_i\) 与 \(\hat{y}_i\) 相关系数的平方。因此 \(0 \le R^2 \le 1\) ： - \(R^2=0\) 表示 \(y_i\) 与 \(x_i\) 之间不存在线性关系； - \(R^2=1\) 表示存在完美的线性关系（例如 \(y_i = x_i + 2\) ）； - \(R^2\) 越大，说明 \(y_i\) 越接近落在OLS回归线上。
不过，在以估计因果效应为目标的研究中，不必过度关注 \(R^2\) 。它虽然是一个有用的描述性指标，但本身并不代表因果性。如果你关心的是某个变量的因果效应，那么 \(R^2\) 无关紧要，此时关键在于零条件均值假设（方程2.28）。
OLS的期望值
到目前为止，我们一直基于总体模型来介绍简单回归，但分析完全是代数层面的，仅基于样本数据。因此，无论模型如何，对样本应用OLS后，残差的平均值始终为0。接下来，我们要进入更深入的统计性质分析，研究OLS估计量的统计性质——这需要基于总体模型，并假设数据是随机抽样得到的。
数理统计学关注的核心问题之一是：不同样本下估计量的表现如何？比如，如果我们反复抽样并计算估计值，平均来看，结果会是“正确”的吗？我们需要求出OLS估计量的期望值（即所有可能随机样本下估计结果的平均值），判断其平均结果是否等于真实参数。这就引出了无偏性 的概念——这是所有估计量都希望具备的优良性质。
\[E(\hat{\beta}) = \beta \tag{2.37}\]
我们的目标是估计 \(\beta_1\) ，即描述 \(y\) 与 \(x\) 关系的总体斜率参数。\(\hat{\beta}_1\) 是基于特定样本得到的估计值，不同样本会得到不同的估计结果。无偏性的含义是：如果我们能从总体中抽取无限多个样本，每次计算一个估计值，那么这些估计值的平均值将等于真实的 \(\beta_1\) 。
OLS估计量要满足无偏性，需要以下几个假设： 1. 参数线性 ：总体模型为 \(y = \beta_0 + \beta_1 x + u\) ，其中 \(\beta_0\) 和 \(\beta_1\) 是未知的总体参数，\(x\) 和 \(u\) 是由数据生成过程产生的随机变量，因此 \(y\) 作为 \(x\) 和 \(u\) 的函数，同样是随机变量。 2. 随机抽样 ：样本 \(\{(x_i, y_i)\}\) 是从总体中随机抽取的，因此每个观测值都满足 \(y_i = \beta_0 + \beta_1 x_i + u_i\) （注意 \(u_i\) 是观测 \(i\) 的未观测误差，不是残差）。 3. 解释变量存在样本变异 ：\(x_i\) 的取值不完全相同，即 \(x\) 的样本方差不为0。如果 \(x_i\) 全部相同（为常数），则OLS估计无法计算（相当于分母为0）。 4. 零条件均值假设 ：总体中，误差项的条件均值为0，即 \(E(u|x) = E(u) = 0\) 。这是因果推断中最关键的假设，也是证明OLS无偏性的核心前提。
下面我们来证明 \(\hat{\beta}_1\) 是 \(\beta_1\) 的无偏估计：
步骤1：写出 \(\hat{\beta}_1\) 的公式 \[\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})y_i}{\sum_{i=1}^n (x_i - \bar{x})^2}\] 令 \(SST_x = \sum_{i=1}^n (x_i - \bar{x})^2\) ，则公式可简化为： \[\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})y_i}{SST_x}\]
步骤2：代入总体模型 \(y_i = \beta_0 + \beta_1 x_i + u_i\) \[\begin{align*}
\sum_{i=1}^n (x_i - \bar{x})y_i &= \sum_{i=1}^n (x_i - \bar{x})(\beta_0 + \beta_1 x_i + u_i) \\
&= \beta_0 \sum_{i=1}^n (x_i - \bar{x}) + \beta_1 \sum_{i=1}^n (x_i - \bar{x})x_i + \sum_{i=1}^n (x_i - \bar{x})u_i \\
&= 0 + \beta_1 SST_x + \sum_{i=1}^n (x_i - \bar{x})u_i
\end{align*}\] 因此： \[\hat{\beta}_1 = \beta_1 + \frac{\sum_{i=1}^n (x_i - \bar{x})u_i}{SST_x}\]
步骤3：求 \(E(\hat{\beta}_1)\) 令 \(w_i = \frac{x_i - \bar{x}}{SST_x}\) ，则 \(\hat{\beta}_1 = \beta_1 + \sum_{i=1}^n w_i u_i\) 。在随机抽样和零条件均值假设下，\(E(u_i|x_1,\dots,x_n)=0\) ，因此： \[\begin{align*}
E(\hat{\beta}_1) &= E\left(\beta_1 + \sum_{i=1}^n w_i u_i\right) \\
&= \beta_1 + \sum_{i=1}^n w_i E(u_i) \\
&= \beta_1 + 0 \\
&= \beta_1
\end{align*}\]
同样，截距项也满足 \(E(\hat{\beta}_0) = \beta_0\) 。
为了更直观地理解这一点，我们可以进行蒙特卡洛模拟。设定总体模型： \[y = 3 + 2x + u \tag{2.38}\] 其中 \(x \sim N(0,9)\) ，\(u \sim N(0,36)\) ，且 \(x\) 与 \(u\) 相互独立。重复抽样1000次并估计OLS，最终 \(\hat{\beta}_1\) 的平均值将趋近于真实值2。
# ============================== # 蒙特卡洛模拟：无偏性 vs 一致性 # ============================== set.seed (123 ) # 固定随机种子 <- 10 # 真实总体均值（我们要估计的参数） <- 1000 # 重复抽样 1000 次 <- c (10 , 50 , 200 , 1000 , 5000 ) # 从小到大的样本量 # 存储结果 <- data.frame ()for (n in n_seq) {for (b in 1 : B) {# 生成数据：正态分布 N(10, 4) <- rnorm (n, mean = true_mu, sd = 2 ) # 估计量 1：样本均值（无偏 + 一致） <- mean (x) # 估计量 2：有偏但一致（偏差随 n 消失） <- mean (x) + 10 / n # 估计量 3：无偏但不一致（永远波动） <- x[1 ] # 保存 <- rbind (result, data.frame (n = n,无偏一致 = est1,有偏一致 = est2,无偏不一致 = est3))# ============================== # 看结果：每个样本量下的平均估计值 # ============================== library (dplyr)%>% group_by (n) %>% summarise (真实值 = true_mu,无偏一致_均值 = mean (` 无偏一致 ` ),有偏一致_均值 = mean (` 有偏一致 ` ),无偏不一致_均值 = mean (` 无偏不一致 ` )) %>% print (digits= 4 )
# A tibble: 5 × 5
n 真实值 无偏一致_均值 有偏一致_均值 无偏不一致_均值
<dbl> <dbl> <dbl> <dbl> <dbl>
1 10 10 10.0 11.0 10.0
2 50 10 10.0 10.2 9.93
3 200 10 10.0 10.1 9.98
4 1000 10 10.0 10.0 9.92
5 5000 10 10.0 10.0 9.94
# 画图看收敛性 library (ggplot2)%>% :: pivot_longer (cols= c (` 无偏一致 ` ,` 有偏一致 ` ,` 无偏不一致 ` )) %>% ggplot (aes (x= factor (n), y= value)) + geom_boxplot () + geom_hline (yintercept= true_mu, color= "red" , lwd= 1 ) + facet_wrap (~ name) + labs (title= "无偏性 vs 一致性" , x= "样本量 n" , y= "估计值" ) + theme_minimal ()
# 加载 tidyverse 包（包含 tibble、dplyr、ggplot2 等工具） library (tidyverse)# 循环执行 1000 次回归模拟 <- lapply (1 : 1000 ,function (x) tibble (# 生成自变量 x：服从 N(0,9)，即均值为0、标准差为3 x = 9 * rnorm (10000 ),# 生成误差项 u：服从 N(0,36)，即均值为0、标准差为6 u = 36 * rnorm (10000 ),# 生成因变量 y：真实模型为 y = 3 + 2x + u y = 3 + 2 * x + u%>% # 对每次生成的数据拟合线性回归模型 y ~ x lm (y ~ x, .)# 提取每次回归的系数（截距与斜率），并转为 tibble 格式 as_tibble (t (sapply (lm, coef))) %>% # 对系数进行统计汇总（查看均值、标准差等） summary (x)
(Intercept) x
Min. :1.789 Min. :1.880
1st Qu.:2.777 1st Qu.:1.975
Median :3.013 Median :2.001
Mean :3.013 Mean :2.001
3rd Qu.:3.266 3rd Qu.:2.026
Max. :4.110 Max. :2.132
# 再次提取系数，绘制直方图观察斜率估计值的分布 as_tibble (t (sapply (lm, coef))) %>% ggplot () + geom_histogram (aes (x), binwidth = 0.01 )
表7给出了1000次重复抽样下，估计值\(\hat{\beta}_1\) 的平均值。你的结果与这里的差异仅源于模拟中的随机性，但整体趋势应与文中展示的一致。尽管每次样本的斜率估计值各不相同，但所有样本的\(\hat{\beta}_1\) 平均值为1.998317，非常接近方程2.38中设定的真实值2。该估计量的标准差为0.0398413，也与回归结果中记录的标准误非常接近。由此可见，估计值本身是重复抽样下系数的平均值，而标准误则是这些重复估计值的标准差。这些系数估计值的分布可在图5中查看。
问题在于，我们永远无法知道自己拿到的是哪一种样本：是“几乎正好等于2”的样本，还是“与2相差较大”的样本？我们永远无法判断自己的估计结果是否接近总体真实值，只能希望当前的样本是“典型样本”，能得到接近真实值的斜率估计值。
但需要明确：无偏性是估计方法的性质，而非单次估计结果的性质 。例如，假设我们估计出教育回报率为8.2%，很容易误以为“8.2%本身是教育回报率的无偏估计”，但严格来说，这是错误的。如果我们假设误差项\(u\) 与教育水平无关，那么得到 \(\hat{\beta}_1=0.082\) 的这个估计方法（规则）是无偏的，而非这一个0.082的数值本身具有无偏性。
迭代期望定律
条件期望函数（CEF）是指在协变量 \(x\) 固定时，结果变量 \(y\) 的均值。
我们来仔细研究这个函数。先明确一下符号表示和基本概念：如前所述，条件期望函数记作 \(E(y_i | x_i)\) 。它是 \(x_i\) 的函数，由于 \(x_i\) 是随机变量，因此条件期望函数本身也是随机的——尽管我们有时会用到 \(x_i\) 的特定取值，例如 \(E(y_i | x_i=8年教育)\) 或 \(E(y_i | x_i=女性)\) 。当存在处理变量时，条件期望函数有两个取值：\(E(y_i | d_i=0)\) 和 \(E(y_i | d_i=1)\) ，但这些只是特殊情况。
条件期望函数的一个重要补充是迭代期望定律（LIE） 。该定律指出，无条件期望可以表示为条件期望的无条件平均。换句话说，\(E(y_i) = E\left\{E(y_i | x_i)\right\}\) 。这个概念很简单：如果你想知道某个随机变量 \(y\) 的无条件期望，只需计算其关于协变量 \(x\) 的所有条件期望的加权和即可。
举个例子：假设女生的平均绩点是3.5，男生的平均绩点是3.2，男女各占总人口的一半。那么： \[\begin{align*}
E[\text{GPA}] &= E\left\{E(\text{GPA}_i | \text{Gender}_i)\right\} \\
&= (0.5 \times 3.5) + (0.5 \times 3.2) \\
&= 3.35
\end{align*}\] 你可能一直在不知不觉中使用迭代期望定律。它的证明也不复杂。设 \(x_i\) 和 \(y_i\) 为连续分布，联合密度函数为 \(f_{xy}(u,t)\) ，\(x=u\) 时 \(y\) 的条件分布为 \(f_y(t | x_i=u)\) ，边缘密度函数分别为 \(g_y(t)\) 和 \(g_x(u)\) 。
\[\begin{align*}
E\left\{E(y | x)\right\} &= \int E(y | x=u)g_x(u)du \\
&= \int \left[\int t f_{y|x}(t | x=u)dt\right] g_x(u)du \\
&= \int t \left[\int f_{y|x}(t | x=u)g_x(u)du\right] dt \\
&= \int t \left[\int f_{xy}(t,u)du\right] dt \\
&= \int t g_y(t) dt \\
&= E(y)
\end{align*}\] 这个证明过程非常清晰：第一行使用了期望的定义；第二行使用了条件期望的定义；第三行交换了积分顺序；第四行使用了联合密度的定义；第五行进行了替换；第六行将联合密度在 \(x\) 的支撑域上积分，得到了 \(y\) 的边缘密度。因此，迭代期望定律可表述为：\(E(y_i) = E\left\{E(y_i | x_i)\right\}\) 。
条件期望函数的分解性质
我们要讨论的条件期望函数的第一个性质是分解性质。迭代期望定律的强大之处在于，它能将随机变量分解为两部分：条件期望函数部分和一个具有特殊性质的残差部分。条件期望函数的分解性质指出： \[y_i = E(y_i | x_i) + \varepsilon_i\] 其中： (i) \(\varepsilon_i\) 与 \(x_i\) 均值独立，即： \[E(\varepsilon_i | x_i) = 0\] (ii) \(\varepsilon_i\) 与 \(x_i\) 的任何函数都不相关。
该定理表明，任何随机变量 \(y_i\) 都可以分解为由 \(x_i\) 解释的部分（条件期望函数）和剩余部分，且该剩余部分与 \(x_i\) 的任何函数正交。我们先证明第(i)部分。注意 \(\varepsilon_i = y_i - E(y_i | x_i)\) ，代入下式： \[\begin{align*}
E(\varepsilon_i | x_i) &= E\left(y_i - E(y_i | x_i) | x_i\right) \\
&= E(y_i | x_i) - E(y_i | x_i) \\
&= 0
\end{align*}\]
定理的第二部分指出，\(\varepsilon_i\) 与 \(x_i\) 的任何函数都不相关。设 \(h(x_i)\) 为 \(x_i\) 的任意函数，则 \(E\left(h(x_i)\varepsilon_i\right) = E\left\{h(x_i)E(\varepsilon_i | x_i)\right\}\) 。由均值独立性可知，内部条件期望项等于零。
条件期望函数的预测性质
条件期望函数的第二个性质是预测性质 ：\(E(y_i | x_i) = \arg\min_{m(x_i)} E\left[(y - m(x_i))^2\right]\) ，其中 \(m(x_i)\) 是 \(x_i\) 的任意函数。换句话说，条件期望函数是在给定 \(x_i\) 时，对 \(y_i\) 的最小均方误差预测 。
我们可以通过在等式右侧加上 \(E(y_i | x_i) - E(y_i | x_i) = 0\) 来展开： \[\left[y_i - m(x_i)\right]^2 = \left[(y_i - E[y_i | x_i]) + (E(y_i | x_i) - m(x_i))\right]^2\] 用更简洁的符号表示为： \[(a - b + b - c)^2\] 将式子展开、整理，并代回原表达式即可得到证明。
对 \(m(x_i)\) 求最小化
对上述函数关于 \(m(x_i)\) 求最小值时，注意第一项 \((y_i - E(y_i | x_i))^2\) 与 \(m(x_i)\) 无关，因此不影响求导结果。第二、三项与 \(m(x_i)\) 有关。令 \(2(E(y_i | x_i) - m(x_i)) = h(x_i)\) ，\(\varepsilon_i = y_i - E(y_i | x_i)\) ，则式子变为： \[\arg\min \varepsilon_i^2 + h(x_i)\varepsilon_i + \left[E(y_i | x_i) + m(x_i)\right]^2\] 对该式求导并令导数为零，得到 \(h'(x_i)\varepsilon_i\) ，根据条件期望函数的分解性质，该项的期望值为零。
方差分析（ANOVA）定理
条件期望函数的最后一个性质是方差分析（ANOVA）定理 ：随机变量的无条件方差，等于条件期望的方差加上条件方差的期望，即： \[V(y_i) = V\left[E(y_i | x_i)\right] + E\left[V(y_i | x_i)\right]\] 其中 \(V\) 表示方差，\(V(y_i | x_i)\) 为条件方差。
线性条件期望函数定理
你现在应该已经知道，在实证研究中使用最小二乘法非常普遍。这是因为回归分析有多种理论支撑。我们已经讨论过其中一种——在误差项满足特定假设时，OLS估计量是无偏的。
我想提出一些略有不同的观点。安格里斯特和皮施克（2009）认为，即使真实的条件期望函数（CEF）本身不是线性的，线性回归依然可能有用，因为它是对条件期望函数的良好近似。让我们来深入拆解一下这个观点。
安格里斯特和皮施克（2009）给出了几个使用回归分析的理由，其中线性条件期望函数定理 可能是最容易理解的。我们假设条件期望函数本身就是线性的，那么会怎样？如果条件期望函数是线性的，那么线性条件期望函数定理指出，总体回归就等于这个线性条件期望函数。既然条件期望函数是线性的，而总体回归又等于它，那么我们当然应该用总体回归来估计条件期望函数。
如果你需要证明这一点，这里给出一个简单的推导：如果 \(E(y_i | x_i)\) 是线性的，那么它可以表示为 \(E(y_i | x_i) = x'\hat{\beta}\) 。根据分解性质，我们有： \[E\left(x\left(y - E(y | x)\right)\right) = E\left(x\left(y - x'\hat{\beta}\right)\right) = 0\] 解这个方程，我们得到 \(\hat{\beta} = \beta\) ，因此 \(E(y | x) = x'\beta\) 。
最佳线性预测定理
还有几个其他的线性定理值得一提。根据条件期望函数的预测性质，条件期望函数是在所有函数中，给定 \(x\) 预测 \(y\) 的最小均方误差预测。因此，总体回归函数是我们在所有线性函数 中能得到的最佳预测。
回归条件期望函数定理
我还想介绍回归的另一个性质：函数 \(X\beta\) 提供了对条件期望函数的最小均方误差线性近似，即： \[\beta = \arg\min_{b} E\left[\left\{E(y_i | x_i) - x_i'b\right\}^2\right]\]
让我们退后一步，把握大局。我讲这些是为了向你说明回归的价值：即使条件期望函数本身不是线性的，线性回归也有其合理性。既然我们无法确定条件期望函数是否为线性，那么这至少是一个值得考虑的理由。回归本质上就是一个将数据转化为估计值的“黑箱”，而我想说明的是，即使在条件不佳的情况下，这个“黑箱”也能产出有价值的结果。接下来，我们再通过另一个著名的“回归解剖定理”来深入了解这个“黑箱”。
回归解剖定理
除了之前讨论的条件期望函数（CEF）和回归定理，我们现在来剖析回归本身。回归解剖定理基于弗里希和沃（1933）以及洛弗尔（1963）的早期研究。我认为通过具体例子和数据可视化来理解这个定理会更直观。它能帮助我们解释多元线性回归模型中各个系数的含义。
假设我们关心家庭规模对劳动供给的因果效应，设定回归模型： \[Y_i = \beta_0 + \beta_1X_i + u_i\] 其中 \(Y\) 为劳动供给，\(X\) 为家庭规模。
如果家庭规模是真正随机的，那么孩子数量与未观测误差项不相关。这意味着，将劳动供给对家庭规模回归得到的估计值 \(\hat{\beta}_1\) ，可以解释为家庭规模对劳动供给的因果效应。我们可以在散点图中绘制所有数据点，回归线的斜率就是数据云的最佳线性拟合。在孩子数量随机的情况下，这个斜率也代表了家庭规模对劳动供给的平均因果效应。
但在现实中，家庭规模几乎不可能是随机的——人们会根据多种因素（如收入、职业规划等）主动选择生育子女的数量。此时，之前的模型假设 \(E(u|X)=E(u)=0\) 很可能不成立。
假设我们认为家庭规模在控制了种族和年龄后是条件随机的，那么模型可写为： \[Y_i = \beta_0 + \beta_1X_i + \gamma_1R_i + \gamma_2A_i + u_i\] 其中 \(Y\) 为劳动供给，\(X\) 为孩子数量，\(R\) 为种族，\(A\) 为年龄，\(u\) 为总体误差项。
要估计家庭规模对劳动供给的平均因果效应，我们需要两个条件：一是包含所有变量的样本数据，二是在给定种族和年龄时，孩子数量是条件随机的。
那么如何解释 \(\hat{\beta}_1\) 呢？回归解剖定理既解释了这个系数估计的含义，也让我们能在二维空间中可视化这个高维问题。
设多元回归模型为： \[y_i = \beta_0 + \beta_1x_{1i} + \dots + \beta_kx_{ki} + e_i\]
构造辅助回归，将 \(x_{1i}\) 对其余自变量回归： \[x_{1i} = \gamma_0 + \gamma_{k-1}x_{k-1,i} + \dots + \gamma_Kx_{Ki} + f_i\] 令 \(\tilde{x}_{1i} = x_{1i} - \hat{x}_{1i}\) 为辅助回归的残差，则参数 \(\beta_1\) 可表示为： \[\beta_1 = \frac{\text{Cov}(y_i, \tilde{x}_i)}{\text{Var}(\tilde{x}_i)}\] 这个公式表明，回归系数估计是一个经过缩放的协方差，缩放因子为辅助回归残差的方差。
证明过程： 首先，利用 \(E[\tilde{x}_{ki}] = E[x_{ki}] - E[\hat{x}_{ki}] = 0\) ，将 \(y_i\) 和残差代入协方差公式： \[\beta_k = \frac{\text{Cov}(\beta_0 + \beta_1x_{1i} + \dots + \beta_Kx_{Ki} + e_i, \tilde{x}_{ki})}{\text{Var}(\tilde{x}_{ki})}\] 由于构造可知 \(E[f_i] = 0\) ，且 \(f_i\) 与其他自变量正交，交叉项均为0。又因为误差项 \(e_i\) 与所有自变量不相关，故 \(E[e_if_i] = 0\) 。
将 \(x_{ki} = E[x_{ki}|X_{-k}] + \tilde{x}_{ki}\) 代入，可得： \[E[\beta_kx_{ki}\tilde{x}_{ki}] = \beta_k \text{Var}(\tilde{x}_{ki})\] 因此： \[\text{Cov}(y_i, \tilde{x}_{ki}) = \beta_k \text{Var}(\tilde{x}_{ki})\] 证明完成。
为了更直观地理解，我们将使用 Stata 的汽车数据集进行演示。
library (tidyverse)library (haven)<- :: read_dta ("data/auto.dta" ) %>% mutate (length = length - mean (length))<- lm (price ~ length, auto)<- lm (price ~ length + weight + headroom + mpg, auto)<- lm (length ~ weight + headroom + mpg, auto)<- %>% mutate (length_resid = residuals (lm_aux))<- lm (price ~ length_resid, auto)<- lm1$ coefficients<- lm2_alt$ coefficients<- lm2$ residuals<- tibble (price = coef_lm2_alt[1 ] + coef_lm1[2 ]* auto$ length_resid,length_resid = auto$ length_resid)<- tibble (price = coef_lm2_alt[1 ] + coef_lm2_alt[2 ]* auto$ length_resid,length_resid = auto$ length_resid)%>% ggplot (aes (x= length_resid, y = price)) + geom_point () + geom_smooth (data = y_multi, color = "blue" ) + geom_smooth (data = y_single, color = "red" )
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
我们来看一下表8中的回归结果，以及我将其称为“短双变量回归”和“长多元回归”中斜率参数的直观可视化。
在价格对车长的短回归中，车长的系数为57.20：车长每增加1英寸，价格就会上涨57美元，这在图6中表现为一条向上倾斜的虚线，其斜率正是57.20。
在回归右侧加入更多变量，即所谓的“控制”这些变量，对你来说很快就会成为自然而然的事。但在这个回归解剖分析中，我希望能让你对回归中“控制变量”的行为有一个不同的理解。
首先，注意在控制了其他变量后，车长的系数发生了符号变化，并且绝对值也变大了。此时车长的效应为-94.5。这表明车长的影响被其他几个变量混淆了，在控制了这些变量后，更长的车实际上反而更便宜。你可以在图6中看到这种情况，多元回归的斜率是负的。
那么这张图里到底发生了什么？首先，它把维度（变量）数量从四个压缩到了只有两个。 它是通过我们之前介绍的回归解剖过程实现的。基本思路是：先运行辅助回归，用它得到残差，然后再按公式 \(\frac{\text{cov}(y_i, \tilde{x}_i)}{\text{var}(\tilde{x}_i)}\) 计算斜率系数。这样我们就能画出辅助回归残差与结果变量观测值的散点图，并在图中画出回归线（图6）。这是一种非常直观的方式，用来展示多元回归中两个变量之间的多维相关关系。注意，图中的黑色实线是负斜率，而双变量回归的斜率是正的。回归解剖定理表明，这两个估计量——一个是多元OLS，另一个是价格与残差的双变量回归——是完全等价的。
OLS估计量的方差
以上内容大致总结了线性回归的核心要点。在零条件均值假设下，我们可以推断出，从样本回归中得到系数的方法是无偏的。这一点很重要，因为它为我们相信回归结果提供了合理依据。但现在，我们需要进一步拓展这一认知，以捕捉抽样过程本身固有的不确定性。这种不确定性通常被称为推断 ，接下来我们将进行讨论。
回想之前的模拟：我们从总体中重复抽样，估计了1000次回归系数，并绘制了这1000个估计值的直方图。系数的均值约为1.998，非常接近数据生成过程设定的真实效应值2；标准差约为0.04。这意味着，在对总体进行重复抽样时，我们会得到不同的估计值，但这些估计值的平均值接近真实效应，且其分布的标准差为0.04。在接下来的学习中，重复抽样下的分布离散度概念是最值得牢记的要点。
在之前讨论的四个假设下，OLS估计量是无偏的。但这些假设本身并不足以告诉我们估计量自身的方差信息。这些假设帮助我们确信，估计的系数平均而言等于总体参数，但要深入探讨估计量的方差，我们需要衡量其在重复抽样分布中的离散程度。正如之前提到的，这引出了方差，并最终推导出标准差。
在四个假设下，我们可以推导出OLS估计量的方差。为简化计算，我们引入第五个假设：同方差性（或常数方差）假设 。该假设规定，总体误差项 \(u\) 在任何解释变量 \(x\) 取值下的方差都是相同的，即： \[V(u|x) = \sigma^2\]
当我第一次学习这部分内容时，对 \(\sigma^2\) 的理解总是很吃力，部分原因是我文科背景，对随机变量的离散程度缺乏直观认识。后来我发现，可以把 \(\sigma^2\) 理解为一个正数（比如2或8），它衡量的是误差项本身的离散程度，即误差项在解释变量条件下的方差是一个有限的正数，代表了除 \(x\) 之外所有影响 \(y\) 的因素的方差。由于零条件均值假设的存在，同方差性假设也可以表示为： \[E(u^2|x) = \sigma^2 = E(u^2)\]
基于假设1、4和5，我们可以得到： \[E(y|x) = \beta_0 + \beta_1x\] \[V(y|x) = \sigma^2\]
也就是说，\(y\) 的均值会随 \(x\) 变化，但在同方差假设下，\(y\) 的方差不随 \(x\) 变化。不过，常数方差假设并非总能成立，需要根据具体情况判断。
定理：OLS的抽样方差
在假设1和2下，我们可以得到： \[V(\hat{\beta}_1|x) = \frac{\sigma^2}{\sum_{i=1}^n (x_i - \bar{x})^2} = \frac{\sigma^2}{SST_x}\] \[V(\hat{\beta}_0|x) = \frac{\sigma^2 \left(\frac{1}{n}\sum_{i=1}^n x_i^2\right)}{SST_x}\]
推导过程如下：我们可以将 \(\hat{\beta}_1\) 表示为： \[\hat{\beta}_1 = \beta_1 + \sum_{i=1}^n w_i u_i\] 其中 \(w_i = \frac{(x_i - \bar{x})}{SST_x}\) 。由于 \(\beta_1\) 是常数，不影响方差，我们可以利用“不相关随机变量和的方差等于方差之和”的性质，以及 \(u_i\) 之间相互独立、不相关的特点，得到： \[V(\hat{\beta}_1|x) = Var\left(\sum_{i=1}^n w_i u_i | x\right) = \sum_{i=1}^n w_i^2 Var(u_i|x) = \sigma^2 \sum_{i=1}^n w_i^2\]
又因为： \[\sum_{i=1}^n w_i^2 = \frac{1}{SST_x}\]
因此最终得到： \[V(\hat{\beta}_1) = \frac{\sigma^2}{SST_x}\]
需要注意的是，这个公式依赖于同方差假设。如果该假设不成立，公式则无效。同时，同方差性是推导该标准方差公式的前提，但不是证明OLS无偏性的前提，后者仅需前四个假设即可。
通常我们更关注 \(\beta_1\) ，影响其方差的因素有两个： - 分子 ：误差方差 \(\sigma^2\) 越大，估计量的方差越大。这意味着，当 \(y\) 和 \(x\) 的关系中“噪音”越多（即 \(u\) 的变异性越大），我们就越难准确估计 \(\beta_1\) 。 - 分母 ：解释变量 \(x\) 的样本变异性 \(SST_x\) 越大，估计量的方差越小。随着样本量 \(n\) 的增加，\(SST_x \approx n\sigma_x^2\) ，因此 \(V(\hat{\beta}_1)\) 会以 \(\frac{1}{n}\) 的速度递减，这也是大样本能降低估计量抽样方差的原因。
\(\hat{\beta}_1\) 的标准差是其方差的平方根： \[sd(\hat{\beta}_1) = \frac{\sigma}{\sqrt{SST_x}}\] 这一指标是置信区间和假设检验中的关键度量。
估计误差方差
在方差公式 \(V(\hat{\beta}_1) = \frac{\sigma^2}{SST_x}\) 中，\(SST_x\) 可以通过样本数据计算，但我们需要估计 \(\sigma^2\) 。由于 \(\sigma^2 = E(u^2)\) ，如果我们能观测到误差项 \(u_i\) ，其无偏估计量将是样本均值： \[\frac{1}{n}\sum_{i=1}^n u_i^2\]
但 \(u_i\) 是不可观测的，因此我们用OLS残差 \(\hat{u}_i\) 来代替： \[u_i = y_i - \beta_0 - \beta_1x_i\] \[\hat{u}_i = y_i - \hat{\beta}_0 - \hat{\beta}_1x_i\]
由于 \(\hat{\beta}_0\) 和 \(\hat{\beta}_1\) 只是总体参数的估计值，因此 \(\hat{u}_i\) 与 \(u_i\) 几乎不可能完全相等，除非是巧合。因此，用 \(\frac{1}{n}\sum_{i=1}^n \hat{u}_i^2\) 来估计 \(\sigma^2\) 是有偏的，因为它没有考虑残差受到的两个约束： \[\sum_{i=1}^n \hat{u}_i = 0\] \[\sum_{i=1}^n x_i \hat{u}_i = 0\]
未观测的误差项不存在这样的约束，因此无偏估计量需要进行自由度调整：残差的自由度为 \(n-2\) （而非 \(n\) ），因此： \[\hat{\sigma}^2 = \frac{1}{n-2}SSR\]
在假设1到5下，该估计量满足 \(E(\hat{\sigma}^2) = \sigma^2\) 。大多数统计软件包的回归结果中，会输出回归标准误： \[\hat{\sigma} = \sqrt{\hat{\sigma}^2} = \sqrt{\frac{SSR}{n-2}}\] 这是对总体误差标准差 \(sd(u)\) 的估计，在Stata中也被称为“均方根误差”（Root Mean Squared Error）。
有了 \(\hat{\sigma}\) ，我们就可以计算 \(\hat{\beta}_1\) 和 \(\hat{\beta}_0\) 的标准差估计值，即标准误： \[se(\hat{\beta}_1) = \frac{\hat{\sigma}}{\sqrt{SST_x}}\] 这些标准误通常会与系数估计值一同出现在回归结果中，报告格式一般为系数下方括号内的数值。
稳健标准误
解释变量 \(x\) 的不同取值下，误差项的方差完全相同，这在现实中可能吗？简单来说，这几乎是不现实的。异方差性才是常态，而非例外，因此我们通常更应该怀疑同方差假设，而不是轻易接受它。可以直接认为误差项永远存在异方差性，然后继续寻找解决方案。
不过这也并非完全是坏消息：在重复抽样下，回归估计量的无偏性并不依赖于误差项的方差假设。前四个假设，尤其是零条件均值假设，保证了系数在重复抽样下的中心趋势等于真实参数（在本书中即因果参数）。问题出在系数的离散程度上。如果不存在同方差性，OLS估计就不再具有最小均方误差的性质，这意味着常规的标准误估计是有偏的。用抽样的比喻来说，系数的分布范围可能比我们想象的更宽。幸运的是，这个问题有解决办法。异方差条件下，方差公式为： \[\text{Var}(\hat{\beta}_1) = \frac{\sum_{i=1}^n (x_i - \bar{x})^2 \sigma_i^2}{SST_x^2}\]
注意式中 \(\sigma_i^2\) 的下标 \(i\) ，这说明方差不是常数。当 \(\sigma_i^2 = \sigma^2\) 对所有 \(i\) 成立时，该公式就退化为常见的 \(\frac{\sigma^2}{SST_x}\) 。但当这一条件不成立时，就会出现异方差误差 。适用于任何形式异方差（包括同方差）的有效估计量为： \[\text{Var}(\hat{\beta}_1) = \frac{\sum_{i=1}^n (x_i - \bar{x})^2 \hat{u}_i^2}{SST_x^2}\] 这可以在OLS回归后直接用数据计算得到。该方法由Friedhelm Eicker、Peter J. Huber和Halbert White提出（White [1980]），解决异方差问题的这一方法有多个名称，最常见的是“稳健”标准误。
聚类稳健标准误
在推断问题上，人们有时会质疑你构建标准误的方式，以此给你制造困扰。但异方差并不是推断过程中唯一需要担心的问题。有些现象并非单独影响每个观测值，而是会影响一组包含多个个体的观测，并且对组内的个体产生共同影响。
例如，你想估计班级规模对学生成绩的影响，但班级中存在一些不可观测因素（如教师）会对所有学生产生同等影响。如果我们可以假设这些不可观测因素在班级之间是独立的，但班级内学生的不可观测因素是相关的，那么我们就需要对标准误进行聚类处理。在进入具体例子之前，我先用一个模拟来展示这个问题。
作为模拟的基准，我们先模拟一组非聚类数据，并对其进行最小二乘估计，这有助于我们清晰理解数据存在聚类时，最小二乘法会出现什么问题。
#- Analysis of Clustered Data #- Courtesy of Dr. Yuki Yanai, #- http://yukiyanai.github.io/teaching/rm1/contents/R/clustered-data-analysis.html library ('arm' )
Loading required package: MASS
Attaching package: 'MASS'
The following object is masked from 'package:dplyr':
select
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Loading required package: lme4
arm (Version 1.15-2, built: 2026-3-22)
Working directory is E:/mixtape-master
Attaching package: 'mvtnorm'
The following object is masked from 'package:arm':
standardize
library ('lme4' )library ('multiwayvcov' )library ('clusterSEs' )
Loading required package: AER
Loading required package: car
Loading required package: carData
The following object is masked from 'package:arm':
logit
The following object is masked from 'package:dplyr':
recode
The following object is masked from 'package:purrr':
some
Loading required package: lmtest
Loading required package: zoo
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Loading required package: sandwich
Loading required package: survival
Loading required package: Formula
Loading required package: plm
The following objects are masked from 'package:dplyr':
between, lag, lead
When using this package, cite:
Justin Esarey and Andrew Menger (2019).
"Practical and Effective Approaches to Dealing with Clustered Data."
Political Science Research and Methods 7(3): 541-549.
URL: https://doi.org/10.1017/psrm.2017.42.
library ('ggplot2' )library ('dplyr' )library ('haven' )
<- function (param = c (.1 , .5 ), n = 1000 , n_cluster = 50 , rho = .5 ) {# Function to generate clustered data # Required package: mvtnorm # individual level <- matrix (c (1 , 0 , 0 , 1 - rho), ncol = 2 )<- rmvnorm (n = n, sigma = Sigma_i)# cluster level <- rep (1 : n_cluster, each = n / n_cluster)<- matrix (c (1 , 0 , 0 , rho), ncol = 2 )<- rmvnorm (n = n_cluster, sigma = Sigma_cl)# predictor var consists of individual- and cluster-level components <- values_i[ , 1 ] + rep (values_cl[ , 1 ], each = n / n_cluster)# error consists of individual- and cluster-level components <- values_i[ , 2 ] + rep (values_cl[ , 2 ], each = n / n_cluster)# data generating process <- param[1 ] + param[2 ]* x + error<- data.frame (x, y, cluster = cluster_name)return (df)
# Simulate a dataset with clusters and fit OLS # Calculate cluster-robust SE when cluster_robust = TRUE <- function (param = c (.1 , .5 ), n = 1000 , n_cluster = 50 ,rho = .5 , cluster_robust = FALSE ) {# Required packages: mvtnorm, multiwayvcov <- gen_cluster (param = param, n = n , n_cluster = n_cluster, rho = rho)<- lm (y ~ x, data = df)<- coef (fit)[2 ]if (! cluster_robust) {<- vcov (fit)<- sqrt (diag (Sigma)[2 ])<- confint (fit)[2 , ]else { # cluster-robust SE <- cluster.vcov (fit, ~ cluster)<- sqrt (diag (Sigma)[2 ])<- qt (.025 , df = n - 2 , lower.tail = FALSE )<- b1 - t_critical* se<- b1 + t_critical* se<- c (lower, upper)return (c (b1, se, b1_ci95))
# Function to iterate the simulation. A data frame is returned. <- function (n_sims = 1000 , param = c (.1 , .5 ), n = 1000 ,n_cluster = 50 , rho = .5 , cluster_robust = FALSE ) {# Required packages: mvtnorm, multiwayvcov, dplyr <- replicate (n_sims, cluster_sim (param = param, n = n, rho = rho,n_cluster = n_cluster,cluster_robust = cluster_robust))<- as.data.frame (t (df))names (df) <- c ('b1' , 'se_b1' , 'ci95_lower' , 'ci95_upper' )<- df %>% mutate (id = 1 : n (),param_caught = ci95_lower <= param[2 ] & ci95_upper >= param[2 ])return (df)
# Distribution of the estimator and confidence intervals <- c (.4 , 0 ) # beta1 = 0: no effect of x on y <- run_cluster_sim (n_sims = 10000 , param = sim_params, rho = 0 )<- ggplot (sim_nocluster, aes (b1)) + geom_histogram (color = 'black' ) + geom_vline (xintercept = sim_params[2 ], color = 'red' )print (hist_nocluster)
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
<- ggplot (sample_n (sim_nocluster, 100 ),aes (x = reorder (id, b1), y = b1, ymin = ci95_lower, ymax = ci95_upper,color = param_caught)) + geom_hline (yintercept = sim_params[2 ], linetype = 'dashed' ) + geom_pointrange () + labs (x = 'sim ID' , y = 'b1' , title = 'Randomly Chosen 100 95% CIs' ) + scale_color_discrete (name = 'True param value' , labels = c ('missed' , 'hit' )) + coord_flip ()print (ci95_nocluster)
%>% summarize (type1_error = 1 - sum (param_caught)/ n ())
如图7所示，最小二乘估计量是以其真实总体参数为中心的。
在显著性水平设为5%的情况下，我们的模拟中错误拒绝原假设 \(H_0: \beta_1=0\) 的概率应为约5%。现在我们来检查置信区间：从图8中可以看到，约95%的95%置信区间包含了 \(\beta_1\) 的真实值（即0）。换句话说，这意味着我们错误拒绝原假设的概率约为5%。
但当我们对聚类数据 使用最小二乘法时，情况会如何变化呢？为了回答这个问题，我们重新模拟数据，使同一聚类内的观测不再是独立的随机抽取。
library ('arm' )library ('mvtnorm' )library ('lme4' )library ('multiwayvcov' )library ('clusterSEs' )library ('ggplot2' )library ('dplyr' )library ('haven' )#Data with clusters <- c (.4 , 0 ) # beta1 = 0: no effect of x on y <- run_cluster_sim (n_sims = 10000 , param = sim_params)<- hist_nocluster %+% sim_cluster_ols
Warning: <ggplot> %+% x was deprecated in ggplot2 4.0.0.
ℹ Please use <ggplot> + x instead.
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
如图 9 所示，当数据存在聚类相关性时，最小二乘估计值的分布比数据独立时更窄。不过，为了更清楚地看到这一点，我们再来看一下置信区间。
library ('arm' )library ('mvtnorm' )library ('lme4' )library ('multiwayvcov' )library ('clusterSEs' )library ('ggplot2' )library ('dplyr' )library ('haven' )#Confidence interval <- ci95_nocluster %+% sample_n (sim_cluster_ols, 100 )print (ci95_cluster_ols)%>% summarize (type1_error = 1 - sum (param_caught)/ n ())
图 10 展示了最小二乘估计得到的 95% 置信区间的分布。可以看到，当数据存在聚类时，有更多的估计结果错误地拒绝了原假设。这是因为在聚类数据下，估计量的标准差被低估了，导致我们错误拒绝原假设的概率过高。那么，我们该如何解决这个问题呢？
library ('arm' )library ('mvtnorm' )library ('lme4' )library ('multiwayvcov' )library ('clusterSEs' )library ('ggplot2' )library ('dplyr' )library ('haven' )#clustered robust <- c (.4 , 0 ) # beta1 = 0: no effect of x on y <- run_cluster_sim (n_sims = 10000 , param = sim_params, cluster_robust = TRUE )<- hist_nocluster %+% sim_cluster_olsprint (hist_cluster_robust)
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
#Confidence Intervals <- ci95_nocluster %+% sample_n (sim_cluster_robust, 100 )print (ci95_cluster_robust)%>% summarize (type1_error = 1 - sum (param_caught)/ n ())
在这个案例中，我们在回归命令里加入了 , cluster(cluster_ID) 这个语法。在深入解释这个语法的作用之前，我们先来看一下置信区间发生了怎样的变化。图11展示了95%置信区间的分布，其中颜色最深的区域依然代表那些错误拒绝原假设的估计结果。
当观测值的误差项在组内存在相关性时，使用这种方法对模型进行最小二乘估计，能让第一类错误率显著降低，回到我们预期的水平。
这就引出了一个很自然的问题：对估计量方差的调整，到底做了什么，能让一类错误率大幅下降？不管它做了什么，显然是有效的！我们用一个例子来深入理解这种调整。考虑以下模型： \[y_{ig} = x_{ig}'\beta + u_{ig}, \quad g=1, \dots, G\] 并且 \[E[u_{ig}u_{jg}']\] 当 \(g \neq g'\) 时，该值为0；当 \(g = g'\) 时，该值为 \(\sigma_{(ij)g}\) 。 我们先按组（cluster）堆叠数据。
\[y_g = x_g'\beta + u_g\] OLS估计量仍然是 \(\hat{\beta} = (X'X)^{-1}X'Y\) 。堆叠数据本身并不会改变估计量本身，但会改变其方差： \[V(\beta) = E\left[(X'X)^{-1}X'\Omega X(X'X)^{-1}\right]\] 基于此，我们可以写出聚类数据的方差-协方差矩阵： \[\hat{V}(\hat{\beta}) = (X'X)^{-1}\left[\sum_{i=1}^G x_g' \hat{u}_g \hat{u}_g'\right] (X'X)^{-1}\]
在实证研究中，对聚类数据进行调整非常普遍，因为聚类数据本身就很常见。在面板数据或双重差分等重复截面设计中，这一点尤为关键。此外，在实验设计中，这也很重要，因为处理的聚合层级往往高于微观数据本身。
在现实世界中，绝不能假设误差项是从同一分布中独立抽取的。你必须了解变量的构建方式，才能为计算标准误选择正确的误差结构。如果有班级规模这类聚合变量，就需要在班级层面聚类；如果某项政策在州一级实施，就需要在州一级聚类。文献中还有更复杂的误差结构，例如多向聚类（Cameron et al., 2011）。
不过，将样本作为标准误基础的概念本身也在发生变化。如今，研究人员越来越少使用随机样本，更多地使用包含总体本身的行政数据，因此抽样不确定性的概念也变得模糊。例如，Manski and Pepper (2018) 指出：“当把州或县作为观测单位时，随机抽样假设并不自然。” 因此，尽管超总体的比喻有助于扩展这些经典的不确定性概念，但数字化行政数据集的普及，促使计量经济学家和统计学家以其他方式思考不确定性。
当以州或县为观测单位时，随机抽样假设并不自然。因此，尽管超总体的比喻有助于扩展这些经典的不确定性概念，但数字化行政数据集的普及，促使计量经济学家和统计学家以其他方式思考不确定性。
Abadie等人（2020）的新研究探讨了基于抽样的标准误概念，在因果推断情境下可能并非思考不确定性的正确方式，他们称之为基于设计的不确定性 。这项研究为接下来两章的内容埋下伏笔，因为它直接提到了反事实的概念。基于设计的不确定性反映了我们不知道如果干预措施不同，反事实中会出现什么值。Abadie等人（2020）推导了基于设计的不确定性的标准误，而非基于抽样的不确定性。结果发现，这些标准误通常更小。
现在，我们来探讨实证研究中使用的这些基本因果概念，并尝试构建工具，理解反事实与因果关系如何协同作用。
注释翻译
Notes 1. 使用一枚六面骰子掷出5点的概率是 \(1/6 = 0.167\) 。 2. 集合符号 \(\cup\) 表示“并集”，指两个事件同时发生。 3. 为什么它们不同？因为0.83是条件概率 \(P(B|A)\) 的近似值，其精确值是 \(0.833...\) （无限循环）。 4. 这里有个很讽刺的故事：有人把蒙提霍尔问题抛给了专栏作家玛丽莲·沃斯·莎凡特。沃斯·莎凡特的智商极高，人们经常给她寄谜题想难住她。她没有用贝叶斯分解，仅凭逻辑就答对了。然而，她的专栏激怒了很多人，批评者写信居高临下地解释她错了，但事实证明错的是他们。 5. 关于回归的更全面综述，请参见伍德里奇（2010）和伍德里奇（2015）。我是站在巨人的肩膀上。 6. 全概率定律要求所有边缘概率之和为1。 7. 只要有可能，我都尽量用“帽子符号（^）”来表示估计统计量，因此用 \(\hat{S}^2\) 而不是 \(S^2\) 。不过，用 \(S^2\) 表示样本方差可能更常见。 8. 这不一定是因果语言。我们一般是指两个随机变量以某种可测量的方式系统地共同变动。 9. 注意，条件期望经过线性函数后，常数项保持不变，而常数乘以 \(x\) 的项也保持不变。这是因为期望算子的第一个性质：\(E[X|X]=X\) 。在零条件均值假设下，这使我们得到 \(E[u|X]=0\) 。 10. 见方程2.23。 11. 注意，我们除以的是 \(n\) ，而不是 \(n-1\) 。换句话说，在用样本计算均值时，没有自由度校正。只有在计算更高阶矩时，才会进行自由度校正。 12. 回忆前面的公式： \[\begin{align*}
\sum_{i=1}^n (x_i-\bar{x})(y_i-\bar{y}) &= \sum_{i=1}^n x_i(y_i-\bar{y}) \\
&= \sum_{i=1}^n (x_i-\bar{x})y_i \\
&= \sum_{i=1}^n x_iy_i - n(\bar{x}\bar{y})
\end{align*}\] 13. 即使 \(u\) 和 \(x\) 相互独立，它也不会恰好为0。可以理解为，\(u\) 和 \(x\) 在总体中是独立的，但在样本中不是。这是由于抽样误差导致样本特征与总体属性略有不同。 14. 使用表6中的Stata代码，你可以自己证明这些代数性质。我鼓励你通过创建等于这些项乘积的新变量，并像处理其他变量一样进行合并，来亲自验证这些代数性质总是成立。 15. 回忆之前关于自由度校正的讨论。 16. 本节是对传统计量经济学教学的回顾。我们讨论它是为了内容的完整性。传统上，计量经济学家通过无偏性和一致性等概念来论证因果关系。 17. 说过我们会经常用到这个结论。 18. 我发现回归中出现这么多 \(\frac{\text{cov}}{\text{var}}\) 项很有意思。它们不断出现，要留意它们。 19. 我在一个数据样本上运行得到的标准误是0.0403616。 20. 强烈建议感兴趣的读者研究安格里斯特和皮施克（2009），他们对迭代期望定律（LIE）有精彩的论述。 21. 我们来举个具体的例子来证明这一点。设 \(h(x_i) = a + \gamma x_i\) 。然后求联合期望 \(E(h(x_i)\varepsilon_i) = E((a+\gamma x_i)\varepsilon_i)\) 。再求条件期望 \(E(a|x_i) + E(\gamma|x_i)E(x_i|x_i)E(\varepsilon|x_i) = a + x_iE(\varepsilon_i|x_i) = 0\) ，条件期望传递后即可得证。 22. 证明参见安格里斯特和皮施克（2009）。 23. 弗里希-沃-洛维尔定理的一个有用证明可在洛维尔（2008）中找到。 24. 虽然随机生孩子听起来很有趣，但我还是建议你在想要孩子的时候再要。可以联系当地的高中健康老师，了解更多能合理控制孩子数量的方法。 25. 这几乎肯定是个荒谬的假设，但请先跟我一起假设下去。 26. 存在一个 \(\sigma\) 的无偏估计量，但它很繁琐，经济学里几乎没人用。参见霍尔茨曼（1950）。 27. 没人再引用怀特（1980）了，就像用微积分时没人引用莱布尼茨或牛顿一样。艾克、休伯和怀特提出的这个方法太有价值了，以至于被吸收进了统计工具包，与原始论文分离开来。 28. 感谢本·奇德米（Ben Chidmi），他帮我在Stata中编写了这个模拟。 29. 通常，我们会诉诸超总体的概念，即我们观测到的总体本身只是某个“超”总体的一次抽样。
