13 回归
13.1 回归的基础
回归是我们拟合直线以解释变异的最常见方式。用一个变量来解释另一个变量的变异，是我们迄今为止在整本书中所做工作的关键部分！
在识别因果效应时，回归是通过控制其他变量来估计两个变量之间关系的最常用方法，这使您能够通过那些控制变量关闭后门路径1。自然，这与我们正在做的事情相关。我们不仅将在本章中进一步讨论它，而且在本书第二部分的其他章节中描述的许多方法本身也是基于回归的。
我们已经在第4章介绍了回归的基础知识。请回顾一下那一章的内容。
以下是该章节的一些关键点：
我们可以用一个变量\(X\)的值来预测另一个变量\(Y\)的值。这被称为用\(X\)解释\(Y\)2。
虽然实现这一目标的方法有很多，但其中一种是通过拟合一条直线或曲线来描述变量间的关系。例如，方程\(Y=\beta_0+\beta_1X\)采用普通最小二乘法（标准线性回归）估计这条直线时，会选择使残差平方和最小的直线——即计算预测误差的平方并求和。线性回归能给出\(X\)和\(Y\)之间关系的最佳线性近似，其近似质量部分取决于真实模型的线性程度。
优点：高效利用数据变异；直线或曲线易于解释
缺点：会忽略部分有意义的变异；若选择的形状与真实关系不符，结果将产生偏差
变量前的系数\(\beta_1\)可解释为斜率，表示\(X\)每增加一个单位，\(Y\)预计相应变化\(\beta_1\)个单位。
当仅有一个预测变量时，斜率估计值等于\(X\)和\(Y\)的协方差除以\(X\)的方差。若存在多个预测变量，计算方法类似，但需额外考虑不同预测变量之间的相关性。
将观测值的预测变量代入回归方程后，可得到预测值\(\hat Y\)，即回归方程解释的\(Y\)部分。\(Y\)与\(\hat Y\)的差值称为”残差”，代表未被解释的部分。
若在回归方程中新增变量\(Y=\beta_0+\beta_1X+\beta_2Z\)，每个变量的系数将通过剔除其他变量解释的部分后剩余的变异来估计。因此，\(\beta_1\)的估计值不再反映 \(Y\)与\(X\)之间的最佳拟合直线，而是反映\(Z\)无法解释的\(Y\)部分 与\(Z\)无法解释的\(X\)部分 之间的关系。这一过程即为”控制\(Z\)的影响”。
如果我们认为\(Y\)和\(X\)之间的关系无法用直线很好地解释，可以改用曲线拟合。普通最小二乘法\(OLS\)能够处理诸如\(Y=\beta_0+\beta_1X+\beta_2 X^2\)这类”参数线性”的模型（注意参数\(\beta_1\)和\(\beta_2\)只是简单地与变量相乘后相加），也可以选择非线性回归方法，如概率单位回归\(probit\)、逻辑回归\(logit\)等众多其他选项。
以上内容基本解释了回归分析的工作原理。只要回归模型与总体模型相符（即我们选择的函数形态能准确描述真实关系），模型通常就能取得不错的效果。不过，还有许多重要内容需要进一步探讨。我们还需要了解哪些内容呢？
接下来我们将补充介绍普通最小二乘法\(OLS\)中尚未深入讲解的几个关键部分：
误差项：模型中的随机扰动因素
抽样变异：样本差异对估计结果的影响
\(OLS\)的统计特性：估计量的统计性质
回归结果的解读：如何正确理解输出结果
二元变量和转换变量的系数解释：特殊变量的参数含义
需要掌握的内容相当丰富，让我们开始逐一探讨。
13.1.1 误差项
回归分析的精妙之处，往往隐藏在那些”未能触及”的领域。此前我将直线拟合描述为：通过确定合适的截距项\(\beta_0\)和斜率系数\(\beta_1\)，最终构建出线性方程\(Y = \beta_0 + \beta_1 X\)。
但在实际数据中，这条拟合直线显然存在不足。它几乎不可能完美预测每一个观测值，更不用说所有数据点了。拟合直线与实际观测值之间必然存在差异——这正是我所说的”未能触及”之处。我们可以将这个差异作为”误差项”\(\varepsilon\)纳入方程，于是完整的回归方程变为：
\[Y = \beta_0 + \beta_1 X + \varepsilon \quad (13.1)\]
这一差异具有两个专业术语：残差 \(residual\) 指代我们使用拟合直线得到的预测值与实际观测值之间的差异；而误差 \(error\) 则表示真实最佳拟合直线与实际值之间的偏差。为何需要区分这两个概念？根源在于抽样变异——由于我们仅拥有有限样本数据，基于当前数据得到的最佳拟合直线，与若能获取总体数据时得到的最佳拟合直线必然存在差异。在实际分析中，我们只能观测到残差，但必须始终牢记误差的概念3。正如\(第3章\)所述，我们的研究目标始终是推断总体特征，而理解总体层面的误差将帮助我们揭示数据背后的真实规律。
残差：模型预测值与实际观测值之间的差异。
误差：实际观测值与”在拥有无限观测数据时所能得到的理想预测值”之间的差异。
\(图13.1\)清晰地展示了残差与误差的区别。图中两条直线分别表示：\(1）\)最初用于生成数据的真实模型；\(2）\)基于随机样本通过普通最小二乘法\(OLS\)估计得到的最佳拟合线。对于每一个数据点：
- 该点到\(OLS\)拟合线的垂直距离称为残差（因为\(OLS\)拟合线代表我们的预测值）
- 该点到真实模型的垂直距离称为误差（即使我们已知\(Y\)与\(X\)的真实关系模型，由于影响\(Y\)的其他因素存在，我们仍然无法完美预测该数据点）
library(tidyverse)## ── 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 3.5.1 ✔ 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 errorslibrary(vtable)## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(Cairo)
library(extrafont)## Registering fonts with Rlibrary(modelsummary)
library(Cairo)
library(lubridate)
library(ggpubr)##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:cowplot':
##
## get_legendlibrary(jtools)# Residual example
set.seed(2897)
# add y and yreal columns
tb <- tibble(x = runif(30)) %>% mutate(y = .6*x + rnorm(30)) %>% mutate(yreal = .6*x) %>% arrange(x)
ggplot(tb, aes(x = x, y = y)) +
geom_point(size=1) +
geom_smooth(formula='y~x', method='lm', se=FALSE, color='red') +
geom_line(aes(y=yreal)) +
geom_segment(aes(x=x, xend=x, y=y, yend=predict(lm(y~x, data=tb))[20]),
data=tb[20,],
size=1,
linetype='dashed') +
geom_segment(aes(x = x, xend = x, y = y, yend = .6*x),
data = tb[22,],
size = 1,
linetype = 'dashed') +
scale_x_continuous(limits = c(0, 1)) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,label = 'Residual',x = .62, y = 0) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,label = 'Error',x = .81, y = -.08) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,label = 'OLS',x = .95, y = .83) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,label = 'True\nModel',x = .95, y = .47) +
theme_pubr() +
labs(x = "X",y = "Y")+
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave('Regression/residual.pdf', width = 6, height = 5,device=cairo_pdf)误差项本质上包含什么？它源自何处？实际上，误差项\(\epsilon\)包含了所有影响\(Y\)但未被模型包含的因素。在方程\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \epsilon\)中，\(Y\)作为因变量位于等式左侧，其值由右侧各项决定。因此，任何影响\(Y\)但不属于\(X\)或\(Z\)的因素，都必须通过\(\epsilon\)来体现。
\(\epsilon\)既包含可观测变量也包含不可观测变量。例如，若真实模型如\(图13.2\)所示，而我们使用\(OLS\)模型\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \epsilon\)进行估计，那么\(\epsilon\)就由\(A\)、\(B\)和\(C\)共同组成。从图中可以明确看出\(A\)、\(B\)和\(C\)都是\(Y\)的成因，但它们并未被纳入当前模型中。
13.1.2 回归假设与抽样变异
关于误差项的思考将引出我们需要对其做出的第一个假设。关于误差项我们将提出若干假设，我会在需要时逐步介绍。使用回归分析的人往往对这些假设极为关注。参加任何研究研讨会，你都会听到无数关于误差项及其假设是否成立的讨论。
第一个要讨论的假设将本书前后两部分联系起来，这就是外生性假设。如果我们希望\(OLS\)对\(\beta_1\)的估计量在期望意义上等于总体参数\(\beta_1\)，就必须满足\(X\)与\(\epsilon\)不相关这一条件4。
外生性假设\(（Exogeneity \space assumption）\)：在回归分析中，该假设认为模型中的变量（或可能仅指处理变量）与误差项\(\varepsilon\)不存在相关性。
这一概念听起来应当很熟悉，因为它本质上重申了识别因果效应的条件！在\(图13.2\)这样的因果图中，所有影响\(Y\)的变量最终都会以某种形式进入回归方程——不论是作为模型中的解释变量，还是被纳入误差项\(\varepsilon\)。若某变量直接出现在回归方程中，它就能阻断所有通过该变量的因果路径。通过向回归模型添加变量，我们实现了”控制该变量”或”将其作为控制变量”的效果。因此，对于回归方程\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \varepsilon\)而言，路径\(X \leftarrow Z \to Y\)就被阻断了。
但如果某个因素仍留在误差项\(\varepsilon\)中且与\(X\)相关，就意味着我们未能阻断这条路径。由于\(A\)包含在\(\varepsilon\)中，且从因果图可知\(A \to X\)，那么我们可以表述为：\(（a）\)我们未能阻断\(X \leftarrow A \to Y\)的后门路径，因此\(X\)的效应未被识别； 或 \(（b）\)\(X\)与\(\varepsilon\)相关，因此具有”内生性”而非外生性，导致\(X\)的效应未被识别5。这两种表述本质上说的是同一件事，只是采用了不同领域的术语6。
同样地，我们可以用两种术语体系说明不需要担心\(C\)是否在模型中。虽然\(C\)存在于误差项\(\varepsilon\)中（所有能预测\(Y\)但未被纳入模型的因素都属于误差项），但我们可以表述为：\(（a）\)不存在通过\(C\)的从\(X\)到\(Y\)的后门路径； 或\(（b）\)\(X\)与\(C\)无关，因此\(C\)不会导致\(X\)具有内生性。
当存在内生性问题时，我们对\(\beta_1\)的估计在期望意义上将不等于总体参数值。这种情况——当估计量在期望意义上给出错误答案时——我们称之为估计偏差。具体而言，这种由于与\(X\)相关的变量留在误差项中导致的偏差，被称为遗漏变量偏差，因为它源于我们在回归方程中遗漏了重要变量。本书先前会表述为”未能阻断通过\(A\)的后门路径”，而用当前术语则表述为”\(A\)导致遗漏变量偏差，但\(C\)不会”。
偏差\(（Bias）\)：无论样本量多大，在期望意义上给出错误答案的估计量被称为有偏估计。
遗漏变量偏差\(（Omitted \space variable \space bias）\)：当与处理变量相关的因素被遗漏在误差项\(\varepsilon\)中而未纳入模型时，所导致的回归估计偏差。
这就是我们对误差项的第一个假设——外生性假设，即误差项\(\varepsilon\)与我们想要估计因果效应的变量都不相关。判断该假设是否成立的方法，正是本书前半部分的核心内容——如何识别因果效应。
当然，这只是关于误差项的众多假设中的第一个。其他许多假设实际上都与\(OLS\)估计量的抽样变异特性相关。
与\(第3章\)讨论的均值类似，回归系数也是估计量。尽管存在真实的总体模型，但由于抽样变异，我们获得的估计值会因样本不同而波动。
幸运的是，我们能够很好地理解这种抽样变异的特征。我们知道可以将观测值视为来自理论分布的抽样，也明白诸如均值等统计量同样可被视为来自理论分布。就均值而言，其服从正态分布（参见\(第3章\)）。如果从总体中抽取大量样本并计算每个样本的均值，这些样本均值的分布将遵循正态分布。据此，我们可以利用获得的样本均值来推断总体均值。
回归系数同样服从正态分布，且我们能够确定该正态分布的均值和标准差——当然，这需要我们对误差项做出若干额外假设。
\(OLS\)系数服从的正态分布具有何种特征7？
在单预测变量的回归模型\(Y = \beta_0 + \beta_1 X + \epsilon\)中，\(OLS\)估计量\(\hat{\beta}_1\)服从均值为\(\beta_1\)、标准差为\(\sqrt{\sigma^2/(var(X)n)}\)的正态分布。其中\(n\)为观测数，\(\sigma\)为误差项\(\epsilon\)的标准差，\(var(X)\)为\(X\)的方差8。不过抽样分布的标准差通常被称为标准误，本文也采用此称谓。
对于多变量模型\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \epsilon\)，其数学表达需引入矩阵代数，但核心思想相同。OLS估计量\(\hat{\beta}_1\)和\(\hat{\beta}_2\)服从联合正态分布，其均值分别为总体参数\(\beta_1\)和\(\beta_2\)，标准差为\(\sigma^2(A'A)^{-1}/n\)对角元素的平方根（矩阵\(A\)包含\(X\)和\(Z\)两列）。可将\((A'A)^{-1}\)理解为”除以\(X\)和\(Z\)的方差及协方差”，正如单变量公式中”除以\(X\)的方差”。
既然无法直接观测误差项\(\varepsilon\)，我们如何确定其标准差\(\sigma\)？实际上我们无法精确获知\(\sigma\)，因此对抽样分布的认识将基于对\(\sigma\)的最佳估计。我们使用可观测的残差替代不可观测的误差，并通过残差的标准差来估计\(\sigma\)的值。
标准误\(（Standard \space error）\)：抽样分布的标准差。
减小\(OLS\)估计量抽样变异的方法仅涉及标准误公式中的三个要素：\(（1）\)降低误差项标准差\(\sigma\)，即提升模型预测精度；\(（2）\)选择高变异度的解释变量\(X\)——\(X\)的剧烈波动更易捕捉\(Y\)的共变规律；\(（3）\)扩大样本量\(n\)。
关于误差项的其他假设有何作用？这些假设实际上是确保前述标准误结论成立的必要条件。我们暂时搁置这个潜在问题，留待后续”你的标准误可能存在问题”章节讨论。当前关于标准误应用的所有结论仍然有效，仅需在计算方法上稍作调整——这正是该章节将要探讨的内容。
13.1.3 普通最小二乘法中的假设检验
好了，我们为何需要了解\(OLS\)系数的分布特性呢？这与研究任何理论分布的目的一致——通过理论分布，可以基于观测数据推断某些总体参数值的可能性。而且，由于\(OLS\)估计量\(\hat{\beta}_1\)的理论分布以总体参数\(\beta_1\)为中心，这意味着我们可以利用样本估计值\(\hat{\beta}_1\)来判定某些总体参数的合理性，例如”\(X\)对\(Y\)的影响系数为\(13.2\)的可能性极低”。
总结当前分析流程（对\(第3章\)方法稍作调整）：
- 选定理论分布——具体为以某个\(\beta_1\)值为中心的正态分布
- 基于观测数据用\(OLS\)估计\(\beta_1\)，得到\(\hat{\beta}_1\)
- 利用理论分布评估：若总体真值为第一步设定值，获得\(\hat{\beta}_1\)的极端程度
- 若概率极低，则初始设定值很可能错误9
该流程可进一步延伸至假设检验概念。假设检验将上述四步规范化，形成判断是否”拒绝”理论分布（从而拒绝初始设定的\(\beta_1\)值）的决策框架10。
在假设检验中，我们首先设定”原假设”——即待检验的初始参数值\(\beta_1\)。实践中通常设为零（尽管并非必须），因此简化为\(\beta_1 = 0\)。接着，选定显著性水平\(\alpha\)。若基于\(\beta_1 = 0\)的理论分布，计算得出获得当前估计结果（或更极端情况）的概率小于\(\alpha\)，则认为原假设不成立。通常取\(\alpha = 0.05\)，即当\(p<=5\%\)时拒绝原假设。拒绝\(\beta_1 = 0\)后能得出什么结论？并非确认估计值准确——仍有众多非零值未被排除。我们仅能判定参数很可能不为零。这一结论具有价值：确认变量间存在某种统计关联。
为演示该过程，我们生成\(200\)个随机观测值，真实模型为\(Y = 3 + 0.2X + \epsilon\)（其中\(\epsilon \sim N(0,1)\)）。然后，假装我们不知道真实的\(\beta_1 = 0.2\) ，通过\(OLS\)回归\(Y = \beta_0 + \beta_1 X + \epsilon\)进行估计。
得到的初始估计值为\(\hat{\beta}_1 = 0.142\)，其标准误\(se(\beta_1) = 0.077\)。因此检验所依据的理论分布为均值为\(0\)、标准差为\(0.077\)的正态分布。
在该正态分布\(N(0, 0.077^2)\)下，观测到的估计值\(\hat{\beta}_1 = 0.142\)位于第\(96.7\)百分位。这意味着，当真实参数\(\beta_1 = 0\)时，这意味着与\(0\)的距离达到\(0.142\)（或更远）的情况，发生的概率为\(2 \times (100\% - 96.7\%) = 6.6\%\)11。若设定显著性水平\(\alpha = 0.05\)，由于\(6.6\%>5\%\) ，我们无法拒绝原假设\(\beta_1 = 0\)——尽管已知真实模型设定为\(\beta_1 = 0.2\)。
# Simulated model
set.seed(31564)
tb <- tibble(X = rnorm(200)) %>% mutate(Y = 3 + .2*X + rnorm(200))
summary(lm(Y~X, data=tb))##
## Call:
## lm(formula = Y ~ X, data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6509 -0.6735 -0.0111 0.6796 3.3546
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.97538 0.07565 39.331 <2e-16 ***
## X 0.14169 0.07680 1.845 0.0666 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.068 on 198 degrees of freedom
## Multiple R-squared: 0.0169, Adjusted R-squared: 0.01193
## F-statistic: 3.403 on 1 and 198 DF, p-value: 0.06656est <- coef(lm(Y~X, data=tb))[2] %>% unname()
se <- summary(lm(Y~X,data = tb))$coefficients[2,2]
# gives the distribution function
pnorm(est, sd=se)## [1] 0.9674659# gives the density
densdata <- tibble(x=-1000:1000/2000) %>% mutate(y=dnorm(x, sd=se))
ggplot() +
geom_ribbon(data=densdata %>% dplyr::filter(x <= -est),
mapping=aes(ymax=y,ymin=0,x=x),
alpha = .5,
color = 'gray') +
geom_ribbon(data=densdata %>% dplyr::filter(x >= est),
mapping=aes(ymax = y,ymin=0,x=x),
alpha = .5,
color = 'gray') +
geom_vline(aes(xintercept=c(-est, est)), linetype='dashed') +
# density
geom_line(data=densdata, mapping=aes(x=x, y=y)) +
scale_x_continuous(
breaks = c(-.4, -est, 0, est, .4),
labels = function(x) scales::number(x, accuracy=.001)
) +
labs(x='X', y='Density') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank())
ggsave('Regression/estimate_vs_dist.pdf', width = 6, height = 5,device=cairo_pdf)
我们可以从\(图13.4\)的另一个角度，看看这是如何运作的。这次，我生成了一堆随机数据样本，展示抽样分布与给定的原假设值对比情况。我神奇地让自己知道\(\beta_1 = 0.2\) 是真实值，用于进行检验。我多次从同一个真实模型\(Y = 3 + 0.2X + \varepsilon\) 生成随机数据，用\(OLS\)估计\(\hat{\beta}_1\) ，并针对\(\beta_1 = 0.2\) 的原假设进行检验。有时，即便原假设确实为真，我还是会拒绝原假设，但这种情况不常发生。实际上，由于\(\alpha = 0.05\) 的决策规则，我恰好有\(5\%\)的可能会拒绝原假设。抽样变异就会导致这样的结果！这是对真实情况的拒绝，我们拒绝真实原假设的频率就是 “假阳性” 率，或者名字很糟糕的 “第一类错误率” 12。
ests <- c()
set.seed(31564)
for (i in 1:1000) {
tb <- tibble(X=rnorm(200)) %>% mutate(Y=3 + .2*X + rnorm(200))
summary(lm(Y~X, data=tb))
ests[i] <- coef(lm(Y~X,data = tb))[2] %>% unname()
}
tb <- tibble(x = density(ests)$x, y=density(ests)$y)
rejl <- .2 - 1.96*se
rejr <- .2 + 1.96*se
p1 <- ggplot() +
geom_ribbon(data=tb %>% dplyr::filter(x <= rejl),
mapping=aes(ymax = y,ymin=0, x=x),
alpha = .5,
color = 'gray') +
geom_ribbon(data=tb %>% dplyr::filter(x >= rejr),
mapping=aes(ymax = y,ymin=0, x=x),
alpha = .5,
color = 'gray') +
geom_line(data=tb, mapping=aes(x = x, y = y)) +
geom_vline(aes(xintercept = .2), linetype = 'dashed') +
annotate(geom = 'label', hjust = .5, x = .2, y = 2, label = 'True\nNull\nValue', family = 'sans', size = 10/.pt) +
annotate(geom = 'text', hjust = .5, x = -.02, y = 1.25, label = 'Estimates\nRejecting\nthe Null', family = 'sans', size = 10/.pt) +
annotate(geom = 'text', hjust = .5, x = .43, y = 1.25, label = 'Estimates\nRejecting\nthe Null', family = 'sans', size = 10/.pt) +
labs(x = 'X',y = 'Density') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank())
rejl <- 0 - 1.96*se
rejr <- 0 + 1.96*se
p2 <- ggplot() +
geom_ribbon(data=tb %>% dplyr::filter(x >= rejr),
mapping=aes(ymax = y,ymin=0, x=x),
alpha = .5,
color = 'gray') +
geom_line(data=tb, mapping=aes(x = x, y = y)) +
geom_vline(aes(xintercept=0), linetype = 'dashed') +
annotate(geom = 'label', hjust = .5, x = 0, y = 2, label = 'False\nNull\nValue', family = 'sans', size = 10/.pt) +
annotate(geom = 'text', hjust = .5, x = .4, y = 1.3, label = 'Estimates\nRejecting\nthe Null', family = 'sans', size = 10/.pt) +
labs(x = 'X',y = 'Density') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank())
plot_grid(p1,p2)
ggsave('Regression/true_false_null.pdf', width = 8, height = 4,device=cairo_pdf)然后，我用同样的这些估计值，针对\(\beta_1 = 0\) 的原假设进行检验。有时，即便原假设是错误的，我也没能拒绝原假设。抽样变异也会导致这种情况！我们没能拒绝错误原假设的情况，就是 “假阴性” 率，或者 “二类错误率” 。
这些图表所展示的是，假设检验其实不是为了确定真实值，也不是要得出绝对的、非真即假的结果。它是为了表明某些理论分布是不太可能的。然后我们可以说，如果这些分布太不可能了，我们就可以不用考虑它们了。
在实际情况中，应用于大多数情境下的回归时，我们讨论的是正态分布的回归系数，以及\(\alpha = 0.05\)的假阳性率（采用双侧检验 ），检验的原假设是系数为\(0\)。所以，要是我们的估计值（或者比\(0\)离得更远的值 ）出现的概率小于\(5\%\)，我们就称之为 “统计显著”，并认为关系是非零的。
概率的计算，源于计算均值为\(0\)、标准差为我们估计出的\(se(\hat{\beta}_1)\)的正态分布的百分位数。如果我们的估计值\(\hat{\beta}_1\)低于第\(2.5\%\)百分位数，或者高于第\(97.5\)百分位数，那么在\(\alpha = 0.05\)时，就是统计显著的。
我们也可以查看得到的精确百分位数并进行评估。比如，如果我们发现对\(\hat{\beta}_1\)的估计处于第\(1.3\)百分位数，那我们可以把它翻倍到\(2.6\)，然后说我们有\(2.6\%\)的概率距离原假设这么远或者更远。这个概率被称为 “\(p\)值” 。\(p\)值越低，仅因抽样变异就得到距离原假设这么远或更远结果的可能性就越小。要是你回看一下\(图13.3\)，整个阴影区域就是\(p\)值，因为它是距离原假设和实际估计值一样远或者更远的区域。如果\(p\)值低于\(\alpha\)，那就具有统计显著性。
\(p\)值\(（p-value）\)：在选定原假设的理论抽样分布下，得到与原假设距离一样远（或更远 ）的估计结果的概率。
利用正态分布的一些已知性质，我们可以走个捷径，不用费劲去计算百分位数。在这种情况下，\(t\) 统计量是系数除以它的标准误，即\(\hat{\beta}_1 / se(\hat{\beta}_1)\) 。当我们把估计值按其标准误缩放时，就把我们带到了标准差为\(1\)的 “标准正态” 理论分布上。这样一来，如果 \(t\) 统计量低于\(-1.96\)（第\(2.5\)百分位数 ）或者高于\(1.96\)（第\(97.5\)百分位数 ），那么在\(95\%\)的置信水平下，就是统计显著的，且存在非零关系。
运用这种解读方式，能让我们读懂回归表格，我们一会儿就会讲到。但在这之前，先就 “统计显著性” 这整个概念说句题外话，或者说一番恳求。
我教过很多学统计方法的学生，也和很多学过统计方法的人交流过，学生所学和我所教内容之间，最大的差异就在于对统计显著性的理解。我觉得这是因为统计显著性很有迷惑性、很诱人13 。它强大、诱人、看似简单，但也很容易被误用，而且极易让人误入歧途、做出不当行为。所以，请牢记以下几点，并不厌其烦地反复念叨，直到它们融入你的思维：
一个估计值不具有统计显著性，不代表它是错的。它只是不具有统计显著性而已。
永远、永远、永远、永远、永远、永远、永远别陷入这样的想法：“哎呀，我的结果不显著，我最好改改分析，弄出个显著的来。” 这很糟糕 14。我很确定教授们总说别这么做，但学生们不知怎的就记成相反的了15。
除了 “真实关系是否非零” 之外，还有很多因素会影响显著性 —— 抽样变异、样本量，当然还有你做分析的方式、对\(\alpha\)的选择等等。显著性检验不是对结果的最终定论16。
统计显著性仅能说明某个特定的原假设值是否不太可能。它无法说明你发现的效应是否有意义。换句话说，“统计显著” 和 “有意义” 不是一回事17 。一个结果显示你的处理让智商提高了 0.000000001 分，不管在统计上是否显著，都不是有意义的效应。
总体而言，请记住：显著性检验的意义在于，不仅要考量估计值本身，还要考量这些估计值的精确性。要是你得到了一个超棒的结果，但标准误极大，那很有可能这只是个巧合，可能与众多真实关系都能契合。这正是显著性检验想要促使人们进行的思考。不过，还有其他方式能让你留意到精确性。你可以直接去思考标准误本身；问问自己，不管估计值与给定的原假设值距离多远，这个估计值是否精确。你可以问问合理的原假设值范围是多少（构建一个置信区间 ），而非聚焦于某一个特定值。你可以完全采用贝叶斯方法，去做那些狂热贝叶斯支持者会做的事儿。显著性检验只是其中一种方法，和其他方法一样，有其优缺点。
13.1.4 回归表格与模型拟合统计量
到目前为止，我们对普通最小二乘法\(OLS\)估计得出的直线，以及其系数的统计性质，已经有了不错的理解。但我们该如何解读整个模型呢？怎样才能一次性理解整个估计出的回归呢？为此，我们可以借助回归结果最常用的呈现方式：回归表格。
我们将使用餐厅和食品检查的数据进行一些回归分析。我们可能会好奇，连锁餐厅是否比门店较少（或仅一家门店 ）的餐厅在卫生检查中获得更好的成绩18 。在本章中，我们会一直使用这些数据。数据的一些基本描述性统计量见\(表13.1\)。我们有检查分数（满分\(100\)分 ）、检查年份，以及连锁餐厅的门店数量的数据。
# 加载必要的库
library(dplyr) # 用于数据处理
library(lubridate) # 用于日期处理
library(stargazer) # 用于生成统计表格（st()函数通常来自此包）##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer# 读取数据
res <- read_csv('Regression/restaurant_inspections.csv',show_col_types = FALSE)
# 数据处理与表格1的回归准备
tab1 <- res %>%
group_by(business_name) %>%
mutate(NumberofLocations = n()) %>% # 计算每个餐厅的门店数量
mutate(Year = year(mdy(inspection_date))) %>% # 从检查日期中提取年份
select(inspection_score, Year, NumberofLocations) %>% # 选择需要的变量
filter(!is.na(inspection_score)) %>% # 移除检查分数为NA的记录
rename(
`Inspection Score` = inspection_score,
`Number of Locations` = NumberofLocations,
`Year of Inspection` = Year
)## Adding missing grouping variables: `business_name`# 生成描述性统计表格
st(tab1,out = 'viewer', title = 'Summary Statistics for Restaurant Inspection Data')## NULL
现在我将运行两个回归分析。第一个回归仅将卫生检查分数对连锁店的门店数量进行回归：
\[ InspectionScore = \beta_0 + \beta_1NumberOfLocations + \epsilon \quad (13.2) \]
第二个回归加入检查年份作为控制变量：
\[ InspectionScore = \beta_0 + \beta_1NumberOfLocations + \beta_2YearOfInspections + \epsilon \quad (13.3) \]
然后，我在\(表13.2\)中展示了两个回归的估计结果 。
m1 <- lm(`Inspection Score`~`Number of Locations`, data = tab1)
m2 <- lm(`Inspection Score`~`Number of Locations` + `Year of Inspection`, data = tab1)
knitr::opts_current$set(label = 'statisticaladjustment-restaurantregs')## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.msummary(list('Inspection Score'= m1,'Inspection Score' = m2),
stars = TRUE,
gof_omit = 'AIC|BIC|Lik',
output = 'Regression/restaurant_regs.tex')## Warning: To compile a LaTeX document with this table, the following commands must be placed in the document preamble:
##
## \usepackage{tabularray}
## \usepackage{float}
## \usepackage{graphicx}
## \usepackage{codehigh}
## \usepackage[normalem]{ulem}
## \UseTblrLibrary{booktabs}
## \UseTblrLibrary{siunitx}
## \newcommand{\tinytableTabularrayUnderline}[1]{\underline{#1}}
## \newcommand{\tinytableTabularrayStrikeout}[1]{\sout{#1}}
## \NewTableCommand{\tinytableDefineColor}[3]{\definecolor{#1}{#2}{#3}}
##
## To disable `siunitx` and prevent `modelsummary` from wrapping numeric entries in `\num{}`, call:
##
## options("modelsummary_format_numeric_latex" = "plain")
## This warning appears once per session.knitr::opts_current$set(label = 'statisticaladjustment-restaurantregs2')## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.msummary(list('Inspection Score' = m1,'Inspection Score' = m2),
stars = TRUE,
gof_omit = 'AIC|BIC|Lik',
output = 'markdown')| Inspection Score | Inspection Score | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 94.866*** | 225.333*** |
| (0.046) | (12.411) | |
| Number of Locations | -0.019*** | -0.019*** |
| (0.000) | (0.000) | |
| Year of Inspection | -0.065*** | |
| (0.006) | ||
| Num.Obs. | 27178 | 27178 |
| R2 | 0.065 | 0.068 |
| R2 Adj. | 0.065 | 0.068 |
| F | 1876.705 | 997.386 |
| RMSE | 6.05 | 6.04 |
在\(表13.2\)中我们能看到什么呢？每一列代表一个不同的回归。第一列结果展示的是方程\(（13.2）\)回归的结果，第二列结果展示的是方程\(（13.3）\)回归的结果。
我们首先注意到的是，每个用于预测检查分数\(（Inspection \space score）\)的变量，都有自己对应的两行数据。“截距\(（Intercept）\)” 也有两行 —— 这是\(\hat{\beta}_0\)，在有些表格里也被称作 “常数项\(（Constant）\)”。第一行显示的是系数估计值。我们对 “门店数量\(（Number \space of \space Location）\)” 变量系数\(\beta_1\)的估计值是\(\hat{\beta}_1 = -0.0019\) 。在两个回归中，这个估计值是一样的 —— 加入 “年份\(（year）\)” 作为控制变量后，估计值的变化小到可以忽略不计。还有一些星号 —— 我一会儿再讲这些星号。
在系数估计值下方，括号里是我们对精确性的度量 。具体来说，在这个表格里，这些是系数的标准误19。本次回归中\(se(\hat{\beta}_1)\)是\(0.0004\)—— 在我看来，精度相当高。衡量系数估计值精确性的方法有好几种。在这里放标准误是最常见的，但有时你也会看到置信区间之类的内容。在有些领域，使用 \(t\) 统计量（系数除以标准误 ）比使用标准误更常见。理想情况下，表格会说明用的是哪种，但也不总是这样 20。
最后，我们看到那些星号，叫做 “显著性星号”。它们能让你一眼看出系数与原假设值\(0\)在统计上是否显著不同。严格来说，这些并非确切的显著性检验。相反，它们是\(p\)值的一种体现 21。\(p\)值越低，星号就越多。
每个星号数量对应的\(p\)值临界值因领域而异，但应像此处一样在表格注释中说明。在许多社会科学领域，也是本书采用的标准是：\(*\)表示\(p\)值小于\(0.1\)\(（10\%）\)，\(**\)表示小于\(0.05\)\(（5\%）\)，\(***\)表示小于\(0.01\%\)\(（1\%）\)22 。所以，要是我们在表格中看到 “门店数量\(（Number \space of \space Location）\)” 系数旁有\(***\)，就像这里的情况，那就意味着，如果我们选定\(\alpha = 0.01\)（或更高 ），我们会拒绝\(\beta_1 = 0\)的原假设。要是看到\(**\)，那在\(\alpha = 0.05\)或更高时，我们会发现统计显著性23。基本上，它衡量的是对于这个系数，在哪些\(\alpha\)值下你会发现显著性 。这些星号能让你一眼看出哪些系数在统计上与\(0\)显著不同，这里所有系数都是如此。
接着看表格下方，有一堆不是系数或精确性度量的内容。具体呈现哪些统计量，不同表格会有差异 24，但总体而言，这些要么是对所进行分析的描述，要么是模型质量的度量。对于前者，常见的是把估计的额外细节列在这儿（比如任何标准误调整 ），而且几乎每个表格都会告诉你用于估计回归的 “观测值数量\(Observations\)”（有时也叫 “\(N\)” ）。
对于后者（衡量模型质量 ），有无数种不同的方法。这里可能纳入了两种最常见的 ——\(R^2\)和调整后\(R^2\)\(（adjusted \space R^2）\)。它们是对模型所解释的因变量方差占比的度量。第一个模型的\(R^2\)是\(0.065\)，这意味着检查分数变异的\(6.5\%\)能由门店数量来预测。要是我们用门店数量去预测检查分数，然后减去预测值，剩下的残差变量就只具有原变量\(（100-0.065）93.5\%\)的方差。调整后\(R^2\)思路类似，但它会针对模型中使用的变量数量进行调整，所以它只统计那些超出仅添加一个随机变量到模型中所能解释的方差部分25 。
最后，我们还有 “F统计量\(（F-statistic）\)”。这是用于进行假设检验的统计量。具体来说，它的原假设是模型中所有系数（截距 / 常数项除外 ）同时为零，然后检验在该原假设下，你的结果有多不可能出现。对于任何还算像样的模型来说，这个统计量不显著的情况都极为罕见，所以我通常会忽略这个统计量。
这个特定表格中没有，但在许多标准回归表格输出样式中会出现的，是残差标准误\(（residual \space standard \space error）\)，有时也叫均方根误差\(（root \space mean \space squared \space error, RMSE）\)。简单来说，它是我们基于残差的标准差，对误差项标准差的估计。我们用\(OLS\)模型得到预测的检查分数值，用实际值减去这些预测值得到残差。然后，计算这些残差的标准差，并针对 “自由度”—— 数据中的观测数量减去模型中的系数数量 —— 做一个小调整。这个数值越大，模型预测中的平均误差就越大。
我们能用所有这些模型质量度量做些什么呢？快速看一眼就行，但总体而言，别太在意这些。它们通常是衡量你的因变量被\(OLS\)模型预测得有多好的指标。但如果你在读这本书，你可能不太关心预测。你感兴趣的是识别因果效应，以及很好地估计特定系数，而非总体上预测因变量。
要是你的\(R^2\)或调整后\(R^2\)值很低，或者残差标准误很高，这告诉你的是，除了你建模的部分，你的因变量还有很多其他影响因素。这值得担心吗？可能吧。如果你原以为自己画出了一个能很好解释因变量、涵盖了所有导致它的因素的因果图，并且基于此选择了模型和要控制的变量…… 结果却得到了一个很小的\(R^2\)，那你最初认为自己真正理解了因变量的数据生成过程的假设可能是错的。但如果你不在乎因变量的大多数成因，并且很确定已经把识别处理效应所需的变量纳入模型了，那么\(R^2\)就没什么重要性。
和统计显著性类似，我发现\(R^2\)值是另一个学生们往往会执着于的指标，让他们很纠结 。但虽然它可以作为不错的诊断工具，绝对没必要执着于此26。肯定别为了最大化\(R^2\)来设计模型 —— 要构建能识别你想要识别的效应、回答你的研究问题的正确模型，而不是构建预测因变量的正确模型。即便你关心预测，\(R^2\)在构建预测模型时也有很多缺陷，不过这就是另一本书的内容了。
了解了回归表格中回归各部分内容的位置后，我们该如何解读结果呢？咱们再看看\(表13.3\)中的回归表格 。
解读\(OLS\)回归结果时，要牢记一个关键表述：“…… 的一个单位变化” 。回归系数的作用是估计两个变量间的线性关系，并以斜率形式呈现结果，也就是预测变量变化一个单位时，因变量的变化情况27 。
要是想精确解读，对变量\(X\)的系数\(\beta_1\)的解读是 “在控制模型中其他变量的情况下，\(X\)变化一个单位，与\(Y\)变化\(\beta_1\)个单位呈线性关联” 。要是想更精确，还能说 “如果两个观测值在模型中其他变量上取值相同，但其中一个的\(X\)值高一个单位，那么\(X\)值高一个单位的观测值，其\(Y\)的平均值会高\(\beta_1\)个单位” 。
咱们从第一行结果开始看，这列只有一个变量，“门店数量\(Number \space of \space Locations\)” 的系数是\(-0.019\) 。先简单想想这里的单位。“门店数量” 指连锁餐厅的门店数，因变量 “检查分数\(Inspection \space Score\)” 是检查员给出的分数，满分\(100\)。
由于模型中没有其他控制变量，这个\(-0.019\)的意思是 “连锁餐厅门店数量每增加一个单位，在\(0-100\)分的评分体系中，检查员评分平均线性减少\(0.019\)分” 。或者说，“对比两家餐厅，若其中一家所属连锁的门店数比另一家多一个，那么前者的检查分数平均会低\(0.019\)分” 。
注意，不管我怎么表述，都刻意避免说 “门店数量增加\(1\)个单位，会使检查员评分降低\(0.019\)分” 。这种说法会暗示我估计出了因果效应。只有依据本书第一部分所学内容，认为我们识别出了门店数量对检查分数的因果效应时，才该用这种表述。
那常数项\(94.866\)怎么理解呢？它是当所有预测变量都为\(0\)时，对因变量的预测值。想想看，因变量的条件均值就是我们把合适的值代入\(OLS\)模型后得到的结果。要是把所有变量都代入\(0\)，其他项都消掉了，就只剩下常数项。所以，对于一家门店数为\(0\)的餐厅，我们预测其检查分数是\(94.866\)。当然，现实中不可能有门店数为\(0\)的餐厅，所以这个结果没什么实际意义 28。
接着看第二行。这里引入了一个控制变量 —— 检查年份。它似乎没怎么改变 “门店数量” 的系数，但会改变解读方式。
现在，我们得把 “对每个系数的解读，是基于我们在模型中‘控制’其他变量这一前提” 这一情况考虑进来 29。我们可以用几种方式表述。
我们可以说 “在控制检查年份的情况下，连锁餐厅门店数量每增加一个单位，在\(0-100\)分的评分体系中，检查员评分平均线性减少\(0.019\)分” 。除了 “\(controlling \space for（控制）\)”，我们也可以说 “\(adjusting \space for （调整）\)”，“\(adjusting \space linearly \space for （线性调整）\)”，甚至 “\(conditioning \space on （以......为条件）\)”30。
我们可以这样说：“对比同一年份的两次检查，若一家餐厅所属连锁的门店数比另一家多一个，那么前者的检查分数平均会低\(0.019\)分。” 毕竟，这就是控制变量的意义所在。我们想要进行同类对比，关闭后门路径，所以要剔除 “门店数量”和”检查分数” 关系中由 “检查年份” 驱动的部分。思路是，通过纳入 “年份” 作为控制变量，我们去除了由年份解释的部分，之后就可以把剩余估计值当作是在对比两次实际年份相同的检查情况。此时，年份已无变异 —— 我们把年份固定住了31 。正因如此，你也常听到这样的解读：“在保持检查年份不变的情况下，门店数量每增加\(1\)个单位，检查分数平均降低\(0.019\)分。”
这是一份解读回归结果的简单指南，重点放在语义表述上。不过，训练自己用这些表述方式解释回归，真的能帮你更好地解读结果。当你开始在更复杂的情境中进行回归分析时（我们会在本章详细讨论 ），同样的直觉也适用。
13.1.5 回归方程中的下标
我遗漏了回归表达式的一个常见特征。在写出回归方程时，人们常会对变量使用下标（比如，\(X_{\text{小文本下标}}\) ）。我们的系数已经有下标了（\(\beta_1\) 里的\(1\)等等 ），但下标也会出现在变量本身里32。在本书中，我大多选择省略这些下标 ，不过当你阅读论文或撰写自己的论文时，熟悉这个概念很重要。
一个回归可能会表示为：
\[ Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i \quad (13.4) \]
这里的 \(i\) 告诉我们数据是按什么索引变化的。换句话说，单个观测值是按什么区分的呢？这里我们用了\(i\)，它一般表示 “个体”，比如，一个人、一家公司或一个国家，具体依情境而定。在这个回归中，\(Y\)和\(X\)因个体不同而有差异。
或者，我们可能会看到：
\[ Y_i = \beta_0 + \beta_1 X_{t-1} + \varepsilon \quad (13.5) \]
这里的t是 “时间段” 的简写。这描述的是一种回归，其中每个观测值对应不同的时间段。\(X_{t - 1}\)表示我们把时间段 \(t\) 的\(Y\)，与前一个时间段（\(t - 1\) ）的\(X\)关联起来。
当下标用于存在多个可能的变化维度时，能发挥最大作用。考虑以下示例：
\[ Y_{it} = \beta_g + \beta_t + \beta_1 X_{it} + \beta_2 W_i + \varepsilon_{it} \quad (13.6) \]
这描述的是一种回归，其中\(Y\)和\(X\)在一组个体\(i\)和时间\(t\)（面板数据集 ）上变化。每个时间段\(t\)有不同的截距\(\beta_t\)，每个\(g\)值（那\(g\)是什么呢？我们得去看研究人员对其回归的描述，但它可能代表个体\(i\)所分组别 ）也有不同截距。我们还看到一个控制变量\(W_i\)。这里的\(i\)下标（且没有时间下标\(t\)）表明\(W\)仅在个体间变化，不随时间改变。或许\(W\)是像出生地这样的变量，个体间不同，但不随时间变化。
13.1.6 将因果图转化为回归
好了，我们已经讲了回归的基本原理。而且，整本书的前半部分都在讲如何用因果图识别因果效应。
人们很容易认为自己自然而然就知道如何把两者联系起来，但我发现学生们往往难以实现这种跨越。别误会，这不是特别难的跨越，但也不是我们能想当然的。
\(图13.5\)中的图表可能会有帮助。在该图表中，我们要思考如何处理变量\(A\)。它是结果变量，还是其他某个变量\(Y\)是结果变量呢？它是处理变量，还是其他某个变量\(X\)是处理变量呢？如果它两者都不是，那么在图表中，以及在回归里（倘若有其位置的话 ），它该处于什么位置呢？
如\(图13.5\)所示，每个变量在回归中可能扮演几种角色之一。它可以是我们的结果变量 —— 我们对这些很了解。结果变量会成为回归中的因变量，即\(Y = \beta_0 + \beta_1 X\)里的\(Y\)部分。从对因果图的讨论中，我们也了解处理变量，它们会成为\(\beta_1 X\)部分。
那其他变量呢？我们的因果图会告诉我们，为了关闭后门路径需要调整的变量（以及那些不应被控制，以免因控制对撞变量而关闭前门路径或开启后门路径的变量 ）。任何应作为控制变量纳入的，都应该…… 纳入作为控制变量。这就是\(Y = \beta_0 + \beta_1 X + \beta_2 A\)里的\(+\beta_2 A\)部分，其中\(X\)是处理变量，\(Y\)是结果变量。
有两样东西我们还没讲到，它们会出现在回归中，但在因果图里不常涉及。首先，要是变量\(A\)不需要被控制，但控制它也不会破坏识别，那该怎么办？这种情况下，我们可以把它作为控制变量纳入，而且或许我们也想这么做 —— 只要\(A\)与\(Y\)相关，加入它就能解释\(Y\)中更多的变异，减少误差项的方差。而由于像\(\hat{\beta}_1\)这样的回归系数的标准误依赖于误差项的方差，标准误就会缩小。
最后，图表中问到 “处理效应是否会因\(A\)的不同取值而有差异” 的那部分该怎么处理呢？这就要用到交互项了。我们会在本章后面讨论交互项。
13.1.7 编写回归代码
在继续深入之前，我们真该实际运行一次回归了。这似乎是个不错的想法。
以下代码块会在\(R\)、\(Stata\) 和 \(Python\)中复现\(表13.2\)和\(表13.3\)。代码会加载餐厅检查数据，计算门店数量，将检查分数对门店数量进行回归，再做另一个纳入年份作为控制变量的回归，然后输出回归表格到文件。
需要注意的是，虽然每种语言中执行回归的代码都相当标准，但生成回归表格的过程通常依赖可下载的软件包。这类软件包有很多。我会使用\(R\)的\(modelsummary\)包（执行 \(install.packages(“modelsummary”)\) 安装）、\(Stata\) 的 \(estout\)（执行 \(ssc \space install \space estout\) 安装 ），以及 \(Python\) 中处理多个回归展示时用到的 \(stargazer\)（执行 \(pip \space install \space stargazer\) 安装 ）。不过，也有一些不错的替代选择，比如 \(R\) 中 \(jtools\) 包的 \(export\_summ\)，或者 \(Stata\) 里的 \(outreg2\)、\(regsave\)，以及 \(Stata \space 17\) 专属的\(collect\)33 。
library(tidyverse); library(modelsummary)
res <- causaldata::restaurant_inspections
res <- res %>%
# Create NumberofLocations
group_by(business_name) %>%
mutate(NumberofLocations = n())
# Perform the first, one-predictor regression
# use the lm() function, with ~ telling us what the dependent variable varies over
m1 <- lm(inspection_score ~ NumberofLocations, data = res)
# Now add year as a control
# Just use + to add more terms to the regression
m2 <- lm(inspection_score ~ NumberofLocations + Year, data = res)
# Give msummary a list() of the models we want in our table and save to the file "regression_table.html"
# (see help(msummary) for other options)
msummary(list(m1, m2), stars=TRUE,output= 'markdown')| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 95.093*** | 247.627*** |
| (0.048) | (12.368) | |
| NumberofLocations | -0.037*** | -0.038*** |
| (0.001) | (0.001) | |
| Year | -0.076*** | |
| (0.006) | ||
| Num.Obs. | 27178 | 27178 |
| R2 | 0.075 | 0.080 |
| R2 Adj. | 0.075 | 0.080 |
| AIC | 174690.2 | 174540.5 |
| BIC | 174714.9 | 174573.4 |
| Log.Lik. | -87342.111 | -87266.266 |
| F | 2198.532 | 1181.426 |
| RMSE | 6.02 | 6.00 |
# Default significance stars are +/*/**/*** .1/.05/.01/.001. Social science
# standard */**/*** .1/.05/.01 can be restored with
msummary(list(m1, m2), stars=c('*' = .1, '**' = .05, '***' = .01),output= 'markdown')| (1) | (2) | |
|---|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | ||
| (Intercept) | 95.093*** | 247.627*** |
| (0.048) | (12.368) | |
| NumberofLocations | -0.037*** | -0.038*** |
| (0.001) | (0.001) | |
| Year | -0.076*** | |
| (0.006) | ||
| Num.Obs. | 27178 | 27178 |
| R2 | 0.075 | 0.080 |
| R2 Adj. | 0.075 | 0.080 |
| AIC | 174690.2 | 174540.5 |
| BIC | 174714.9 | 174573.4 |
| Log.Lik. | -87342.111 | -87266.266 |
| F | 2198.532 | 1181.426 |
| RMSE | 6.02 | 6.00 |
import pandas as pd
import statsmodels.formula.api as smf
from stargazer.stargazer import Stargazer
from causaldata import restaurant_inspections
res = restaurant_inspections.load_pandas().data
# Perform the first, one-predictor regression
# use the sm.ols() function, with ~ telling us what the dependent variable varies over
m1 = smf.ols(formula = 'inspection_score ~ NumberofLocations', data = res).fit()
# Now add year as a control
# Just use + to add more terms to the regression
m2 = smf.ols(formula = 'inspection_score ~ NumberofLocations + Year', data = res).fit()
# Open a file to write to
f = open('Regression/regression_table.html', 'w')
# Give Stargazer a list of the models we want in our table and save to file
regtable = Stargazer([m1, m2])
f.write(regtable.render_html())## 1743f.close()* ssc install causaldata if you haven't yet!
causaldata restaurant_inspections.dta, use clear download
* Create NumberofLocations
by business_name, sort: g NumberofLocations = _N
* Perform the first, one-predictor regression
regress inspection_score NumberofLocations
* Store our results for estout
estimates store m1
* Now add year as a control
regress inspection_score NumberofLocations year
estimates store m2
* use esttab to create a regression table
* note esttab defaults to t-statistics, so we use
* the se option to put standard errors instead
* we'll save the result as regression_table.html
* (see help esttab for other options)
esttab m1 m2 using regression_table.html, se replace13.2 回归的进阶应用
回归是一种工具 —— 一种灵活性很高的工具。而我们其实才学会一种使用它的方式！本节甚至不会改变我们使用它的方式，只是会改变纳入其中的变量。这本身就会开启一个全新的世界。
到目前为止，我们讨论的是\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \varepsilon\)形式的回归。我们仍会围绕这个展开。但一直以来，我们纳入回归的都是连续变量本身。有没有其他情况呢？
首先，我们可能没有连续变量。相反，我们可能有离散变量，最常见的是二元变量 —— 真或假，而非具体数值。如果\(Z\)是 “此人有金色头发”，该如何解读\(\beta_2\)呢？或者，如果 “头发颜色” 不是我们连续测量的变量，而是有 “黑色”“棕色”“金色”“红色” 等类别，我们该如何控制这个变量呢？
其次，我们可能不会直接使用变量，而是先对其进行转换（transform） 。早在\(第4章\)我们就讨论过，有些关系用直线可能无法很好地描述。有时我们需要曲线。具体该怎么做呢？
第三，关系本身可能会受其他变量影响。比如，在\(第4章\)我们讨论过 \(Emily \space Oster\) 的一项研究，研究问题是 “服用维生素\(E\)和健康结果之间的关系，在医生推荐服用维生素\(E\)的时期，是否比不推荐的时期更强？” 我们可以用\(\text{Outcomes} = \beta_0 + \beta_1 \text{Vitamin}E + \varepsilon\)来建模服用维生素\(E\)和健康结果之间的关系。但该如何建模这种关系随时间的变化呢？我们需要交互项\(（interaction）\) 。
所以，这三件事 —— 全都与我们在普通最小二乘法模型中用作预测变量的变量类型有关 —— 将是我们要着手探讨的内容起点。一旦讲完这些，我们会接着探讨如何处理标准误，以及当因变量也不是连续变量时该怎么办 。
13.2.1 处理离散变量
二元变量在社会科学中简直无处不在。你接受处理了吗？你是左利手还是右利手？你是男性还是女性34？你是天主教徒吗？你结婚了吗？
二元变量对因果分析尤为重要，因为我们往往感兴趣的很多原因，本质上就是二元的。你接受处理了吗？
任何时候，只要涉及 “是” 或 “否” 的情况（也就是当我们进行某种定性描述时 ），我们就是在处理二元变量。这些变量可以像普通变量一样纳入回归模型。所以，我们仍可以用\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \varepsilon\)的形式，只不过现在\(X\)或\(Z\)（或两者 ）可能只能取两个值：\(0\)（“否”/ 假 ）或 \(1\)（“是”/ 真 ）35 。如果二元变量是控制变量，我们可以像平常那样看待它 —— 我们只是在关闭通过该变量的后门路径。但要是我们对这个二元变量的效应感兴趣，该如何解读它的系数呢？
简而言之，它是取值为真和取值为假时，因变量的差异 36。比如，如果我们进行\(\text{Sales} = \beta_0 + \beta_1 \text{Winter} + \varepsilon\)这样的回归，其中\(\text{Winter}\)是二元变量（冬季时为\(1\)，非冬季时为\(0\) ），那么\(\beta_1\)表示冬季的销售额平均比非冬季高多少 37。
当我们估计这个模型时，如果得到\(\hat{\beta}_1 = -5\)，那就意味着，平均而言，冬季的销售额比非冬季低\(5\)。看看\(表13.4\)中的模拟数据。注意到 “冬季\(（Winter）\)” 的系数\(（-5）\)，刚好是 “非冬季\(（Not \space Winter）\)” 的均值\(（15）\)和 “冬季\(（Winter）\)” 的均值\(（10）\)之间的差值。还要注意，截距项的系数\(（15）\)是所有变量为\(0\)时的预期均值 —— 当\(\text{Winter}\)为\(0\)时，我们处于非冬季，非冬季的平均销售额是\(15\)。
有一点很重要，需要注意：我们只纳入了 ” 是 / 否” 问题的其中一方。模型是\(\text{Sales} = \beta_0 + \beta_1 \text{Winter} + \varepsilon\)，而非\(\text{Sales} = \beta_0 + \beta_1 \text{Winter} + \beta_2 \text{NotWinter} + \varepsilon\) 。
为啥这么做呢？首先，想想看要是那样（纳入双方）该怎么解读。我们想让\(\beta_1\)用来对比 “冬季\(（Winter）\)” 和 “非冬季\(（Not \space Winter）\)”，可\(\beta_2\)是对比 “非冬季” 和 “冬季”。那…… 它们之间有啥区别呢？解读起来会非常混乱。其次，这是为了满足\(OLS\)的另一个假设 ——无完全多重共线性 。也就是说，你不能用模型中的其他变量，对模型里的任何一个变量做出完美的线性预测。在这里，\(\text{Winter} + \text{NotWinter} = 1\) 。不过，我们的模型里没有\(1\)这个变量，所以没问题，是吧？错！实际上，一个全为\(1\)的变量一直潜藏在我们的模型中 —— 它正和常数项相乘呢！
我们为啥不能有那样的线性组合呢？因为\(OLS\)没办法搞清楚该如何估计参数。假设冬季的平均销售额是\(10\)，非冬季是\(15\)。要是你的回归是\(\text{Sales} = \beta_0 + \beta_1 \text{Winter} + \beta_2 \text{NotWinter} + \varepsilon\) ，你可以用\(\beta_0 = 15\)、\(\beta_1 = -5\)、\(\beta_2 = 0\)得到完全相同的预测结果。也可以用\(\beta_0 = 10\)、\(\beta_1 = 0\)、\(\beta_2 = 5\) 。或者\(\beta_0 = 3\)、\(\beta_1 = 7\)、\(\beta_2 = 12\) 。还有无数其他组合方式！虽然这些组合都能 “起作用”，但 OLS 没办法从这么多选项里选出一个估计值。它没法再给你一个最佳估计了。所以，我们得去掉\(\text{NotWinter}\)，限制它的选项，实际上就是迫使它选择\(\beta_0 = 15\)、\(\beta_1 = -5\)、\(\beta_2 = 0\)这种形式 。
我们可以稍微拓展一下对二元变量的直观理解，让自己有能力在模型中纳入分类变量\(（categorical \space variables）\) 。二元变量就是简单的是 / 否。但分类变量可以有任意数量的类别。这类变量包括像 “你住在哪个国家？” 或者 “你出生于哪个十年？” 这样的问题。你没法用是或否来回答这些问题，但它们是离散且互斥的类别 38。
我们可以通过给每个类别分配专属的二元变量 ，来处理分类变量。所以，对于 “你住在哪个国家？” 这个问题，不用一个变量来表示，而是用 “你住在法国吗？” 一个变量、“你住在冈比亚吗？” 另一个变量、“你住在新西兰吗？” 又一个变量，依此类推。
\[ Income = \beta_0 + \beta_1 France + \beta_2 Gambia + \beta_3 NewZealand + \dots \]
如果我们把分类变量当作控制变量纳入模型，通过为每个类别设置这些二元变量，就可以说，我们关闭了像本例中 “你住在哪个国家” 这类变量所涉及的后门路径。
要是我们想解读这些二元类别中某一个的系数，该怎么做呢？这比二元变量的情况要稍微复杂一点。就像我们不能在模型中同时纳入二元变量的正反两面（模型里只能放 “冬季\(（Winter）\)”，不能同时放 “冬季” 和 “非冬季\(（NotWinter）\)” ）一样，我们也不能把分类变量的所有类别都纳入模型。得去掉其中一个类别。被去掉的那个类别就成了 “参照类别\(（reference \space category）\)” 。所有系数都是相对于该类别 的差值。
之前，对于二元变量，系数表示的是 “是” 和 “否” 之间的差异。现在，对于有参照类别的分类变量，系数表示的是 “本类别” 和 “参照类别” 之间的差异。
比如，要是我们从回归模型里去掉 “法国\(（France）\)”39 ，法国就成了参照类别。接着，我们估计\(\text{Income} = \beta_0 + \beta_1 \text{Gambia} + \beta_2 \text{NewZealand} + \dots\) ，得到\(\hat{\beta}_1 = 0.5\)，\(\hat{\beta}_2 = 3\) 。这就意味着，冈比亚\(（Gambia）\)的平均收入比法国的平均收入高\(0.5\)，新西兰\(（New \space Zealand）\)的平均收入比法国的平均收入高\(3\)。这两种解读都是相对于参照类别而言的，而非相互对比。
我们可以从\(表13.5\)的模拟数据中看到这一点，该表只包含法国、冈比亚和新西兰的数据。注意，冈比亚的平均收入\(（30.5）\)比法国的平均收入\(（30）\)高 \(0.5\) ；当法国作为参照类别时，冈比亚的系数是 \(0.5\) 。同样，新西兰的平均收入\(（33）\)比法国高 \(3\) ，系数就是 \(3\) 。截距项是所有预测变量都为\(0\)时的平均收入 —— 也就是既不是冈比亚，也不是新西兰，那就只能是法国。法国的平均收入是\(30\)，所以截距项是\(30\) 。
这种参照类别机制意味着，类别上的系数（以及它们的显著性 ）本身相对没什么独立意义。它们只有在相互对比时才有意义。要是你改变参照类别，系数会完全改变。要是我们把新西兰设为参照类别，冈比亚的系数就会从 \(0.5\) 变成 \(-2.5\) 40。
要是你想知道一个分类变量整体是否有显著效应，不能看单个系数。相反，得看所有类别系数 。这会以 “\(联合F检验（joint \space test）\)” 的形式呈现，也就是把包含该分类变量作为预测因子的模型，和去掉该分类变量的模型的预测能力做对比。
在 \(R\) 中，\(car\)包的 \(linearHypothesis\) 命令可以执行联合\(F\)检验。你得把分类变量各系数的名称完整列出来（花点功夫用 \(paste\) 和 \(unique\) 自己就能做出来，或者用\(matchCoefs()\)函数 ）。\(fixest\)包的 \(wald()\) 用起来更简单，能轻松对变量名进行模式匹配，不过这也需要用 \(fixest\) 的 \(feols()\) 回归函数，而非\(lm()\) 。在 \(Stata\) 中，\(test \space i.cat\)会对\(cat\)变量的所有类别进行联合显著性检验。在 \(Python\) 中，要是你用 \(smf.ols().fit()\) 创建了回归结果对象（假设叫\(ml\)），可以用 \(ml.f_test()\)做联合\(F\)检验。你得看看文档，了解如何创建一个矩阵或字符串，指明要检验的因子变量。
13.2.2 多项式
有时候，一条直线是不够的。\(OLS\)假设因变量和每个预测变量之间的关系可以用一条直线（线性）来描述。如果情况并非如此，\(OLS\)的表现就会很差。比如，看看\(图13.6\)。真实的关系显然遵循右侧画出的曲线。要是我们试图用一条直线来估计这种关系，得到的结果拟合效果会很差，就像左侧的直线所示。
res <- read_csv('Regression/restaurant_inspections.csv', show_col_types = FALSE)
#table 1 regressions
tab1 <- res %>%
group_by(business_name) %>%
mutate(NumberofLocations = n()) %>%
mutate(Year = year(mdy(inspection_date))) %>%
select(inspection_score, Year, NumberofLocations) %>%
dplyr::filter(!is.na(inspection_score)) %>%
rename(`Inspection Score` = inspection_score, `Number of Locations` = NumberofLocations, `Year of Inspection` = Year)## Adding missing grouping variables: `business_name`m1 <- lm(`Inspection Score`~`Number of Locations` + I(`Number of Locations`^2), data = tab1)
knitr::opts_current$set(label = 'statisticaladjustment-restaurantpoly')## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.msummary(m1,
output = 'markdown',
coef_map = c('(Intercept)' = 'Constant',
'`Number of Locations`' = 'Number of Locations',
'I(`Number of Locations`^2)' = 'Number of Locations Squared'
),
title = 'A Regression with a Second-Order Polynomial using Restaurant Inspection Data',
stars = TRUE,
gof_omit = 'AIC|BIC|Lik')| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| Constant | 97.518*** |
| (0.059) | |
| Num.Obs. | 27178 |
| R2 | 0.194 |
| R2 Adj. | 0.194 |
| F | 3279.008 |
| RMSE | 5.62 |
不过，对于这类关系，我们仍可以使用普通最小二乘法。只是需要在回归方程中小心地对曲线形状进行建模。
我们主要有两种方法可以做到这一点：要么添加多项式项\(（polynomial \space terms）\)，要么对数据进行转换\(（transform）\) 。
多项式\(（polynomial）\)指的是在一个方程中，同一个变量以自身形式，以及自身的幂次形式出现。比如，\(\beta_1 X + \beta_2 X^2 + \beta_3 X^3\) 就是一个 “三阶多项式”，因为它包含\(X\)的一次方（\(X\)）、二次方（\(X^2\) ）和三次方（\(X^3\) ） 4142。
通过添加这些多项式项，我们能够拟合出非直线的曲线。实际上，添加足够多的项，几乎可以模拟出任何形状 43。这一点在\(图13.7\)中可以看到。直接观察数据点，显然这是一种非线性关系。在数据点上方，我们有线性模型（\(Y = \beta_0 + \beta_1 X\) ）拟合出的最佳直线，也就是那条直线；还有二阶多项式模型拟合出的弯曲虚线（\(Y = \beta_0 + \beta_1 X + \beta_2 X^2\) ）；以及三阶多项式模型拟合出的弯曲实线（\(Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 X^3\) ）
这些曲线显然能更好地拟合这种关系，而我们要做的仅仅是添加那些多项式项而已！
# Polynomial example
set.seed(1000)
tb <- tibble(X = runif(100)) %>% mutate(Y = X + 4*X^2 - 5*X^3 + .3*rnorm(100))
ggplot(tb, aes(x = X, y = Y)) +
geom_point() +
geom_smooth(formula = 'y ~ x', method='lm', color = 'black', se = FALSE) +
geom_smooth(formula = y ~ x + I(x^2), color = 'gray', method='lm', se = FALSE, linetype = 'dashed') +
geom_smooth(formula = y ~ x + I(x^2) + I(x^3), method='lm', color = 'red', se = FALSE) +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks.y = element_blank())
ggsave('Regression/linear_square.pdf', width = 6, height = 4,device=cairo_pdf)这就是添加多项式项的原因。添加多项式项，你就能更好地用一条线去拟合非线性关系。但这就给我们留下了两个紧迫的问题：\(（1）\)一旦加入多项式项，该如何解读回归模型；\(（2）\)为啥不一直添加一大堆多项式项呢？
该如何解读含多项式的模型？ 恐怕对于这个问题，我们得用一点微积分知识。拿我们之前的三次模型来说：
\[ Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 X^3 \quad (13.7) \]
我们通常对像\(\hat{\beta}_1\)这样的回归系数估计值的解读是 “在其他条件不变的情况下，\(X\)增加一个单位，\(Y\)会增加\(\hat{\beta}_1\)个单位” 。但这就有问题了！我们没法 “让其他条件不变”，因为要改变\(X\)，就必然会改变\(X^2\)和\(X^3\) 。所以，对\(X\)效应的解读必须考虑其每个多项式项的系数。当回归模型中纳入了\(X\)的多项式时，\(X\)各项的单个系数本身意义不大，必须结合起来解读。
我们该如何把它们结合起来解读呢？让我们试着回到 “在其他条件不变的情况下，\(X\)增加一个单位，\(Y\)会有（???）单位的变化” 这样的解读方式。那（???）是什么呢？嗯，要弄清楚一个变量如何随另一个变量变化，我们通常可以借助什么工具？导数！如果我们对\(Y\)关于\(X\)求导，就能知道\(X\)变化一个单位与\(Y\)有怎样的关联 44。
\[ \frac{\partial Y}{\partial X} = \beta_1 + 2\beta_2 X + 3\beta_3 X^2 \quad (13.8) \]
\(X\)变化一个单位，\(Y\)会有 \((\beta_1 + 2\beta_2 X + 3\beta_3 X^2)\) 单位的变化45 。注意，这种关系会随着 \(X\) 取值的不同而变化。如果 \(X = 1\)，那么效应就是 \(\beta_1 + 2\beta_2 + 3\beta_3\) 。但如果 \(X = 2\)，效应就是 \(\beta_1 + 4\beta_2 + 12\beta_3\) 。\(X\) 对 \(Y\) 的效应取决于当前 \(X\) 的取值。这和我们预期的一致。回顾\(图13.7\) ，随着 \(X\) 增大，一开始 \(Y\) 也会增大，但之后 \(Y\) 停止增长并开始下降。当我们沿着 \(X\) 轴移动时，\(X\) 变化一个单位所产生的效应理应发生变化。
咱们用一个更具体的例子来说明。回到\(表13.2\) 的回归，但添加一个二次项，如\(表13.6\)所示。
在这里，我们看到线性项的系数是\(-0.0802\)，二次项的系数是\(0.0001\) 。这意味着，连锁餐厅门店数量每增加一家，健康检查分数的变化量为\(\beta_1 + 2\beta_2 \times \text{门店数量}\)，也就是\(-0.0802 + 0.0002 \times \text{门店数量}\) 。所以，对于只有\(1\)家门店的连锁，新增第\(2\)家门店时，健康分数会降低\(-0.0802 + 0.0002×1 = -0.0800\) 。但对于有\(1000\)家门店的连锁，新增第\(1001\)家门店时，健康分数会增加\(-0.0802 + 0.0002×1000 = 0.1198\) 。这里不存在单一的效应，它取决于我们处于\(X\)轴的哪个位置。如果忽略二次项，就会忽视这一点。
多项式似乎很有用！ 那为啥不一直添加一大堆多项式项呢？既然要这么做，为啥只停留在二次项和三次项呢？为啥不给所有模型都加上四次、五次、六次或七次多项式项呢？
有几个合理的理由。一方面，模型会变得更难解读。当然了，如果模型是合适的，我们还是会想用它。
另一方面，高阶多项式项往往没什么实际作用。以\(图13.7\)为例，二阶多项式模型和三阶多项式模型拟合出的曲线几乎完全一样46 。所以在很多情况下，添加越来越多的多项式项，会让模型更复杂，却不会实际改善拟合效果。
还有一个问题是，添加多项式项会给数据带来压力，可能导致”过拟合”。过于灵活的模型会为了拟合数据中的噪声，扭曲成奇怪的形状，结果得到更差的模型。比如，想象一下，用\(100\)阶多项式去拟合一个只有\(100\)个观测值的数据集。我们能完美预测每个点，但显然，这样的曲线没法告诉我们什么有用信息。
一般来说，添加更多多项式项会让模型对数据中的微小变化变得非常敏感，还会使模型在观测数据的边缘附近做出奇怪的预测。以\(图13.8\)为例，它以\(图13.7\)为基础，添加了一个十次多项式的回归拟合。其形状看起来并没有比二次多项式更好地拟合数据。而且，我们看到奇怪的上下波动，似乎是在追逐单个观测值，而非整体关系。最糟糕的是最右侧，在最后一刻，曲线突然大幅上升，试图拟合数据最右侧的单个点。多项式中这些额外的次数不只是没必要，还在让模型变得更差。
ggplot(tb, aes(x = X, y = Y)) +
geom_point() +
geom_smooth(method='lm', formula = y ~ x, color = 'black', se = FALSE) +
geom_smooth(method='lm', formula = y ~ x + I(x^2), color = 'black', se = FALSE, linetype = 'dashed') +
geom_smooth(method='lm', formula = y ~ x + I(x^2) + I(x^3), color = 'black', se = FALSE) +
geom_smooth(method='lm', formula = y ~ poly(x, degree = 10), color = 'black', se = FALSE) +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks.y = element_blank())
ggsave('Regression/ten_degree.pdf', width = 6, height = 4,device=cairo_pdf)所以，我们不想要 “太多” 多项式项。但多少才是恰到好处呢？理想情况下，你希望多项式项的数量刚好能对变量间关系的实际真实潜在形状进行建模。但我们实际上并不知道真实形状是什么，所以这没什么用。
要解决该纳入多少个多项式项（或者是否纳入多项式项 ）的问题，一个不错的方法是可视化\(（graphically）\) 。当观测值不太多时，这种方法效果最好。你可以画个散点图，看看数据的形状。然后，加上回归线，看看它是否能解释你看到的形状。要是不行，就试着添加另一个多项式项，直到能解释为止。
你也可以通过可视化分析残差。先尝试一个回归模型，计算残差，然后将残差在\(y\)轴、自变量在\(x\)轴作图。要是有足够多的多项式项，\(X\)和残差之间 不应该 有任何明显的关系。这在\(图13.9\)中能看到。\(图13.8\)中的同一组数据，用线性\(OLS\)和二阶多项式\(OLS\)进行估计，然后将这些回归的残差对\(X\)作图。线性回归后，左侧残差和\(X\)之间的关系仍有明显形状。说明我们的直线不够弯曲。但右侧，用了二阶多项式后，残差只是围绕\(0\)的一堆噪声。这就够了，我们可以停止（添加多项式项 ）了！
p1 <- ggplot(tb, aes(x=X, y = resid(lm(Y~X, data=tb)))) +
geom_point(na.rm=TRUE) +
geom_hline(aes(yintercept=0), linetype = 'dashed', color = 'black')+
scale_y_continuous(limits = c(-.75,.75)) +
labs(y = 'Residual',title = 'Residuals from Linear OLS') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks.y = element_blank())
p2 <- ggplot(tb, aes(x = X, y = resid(lm(Y~X + I(X^2), data = tb)))) +
geom_point(na.rm=TRUE) +
geom_hline(aes(yintercept=0), linetype = 'dashed', color = 'black')+
scale_y_continuous(limits = c(-.75,.75)) +
labs(y = 'Residual', title = 'Residuals from Two-Order Polynomial') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family="sans"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks.y = element_blank())
plot_grid(p1,p2)
ggsave('Regression/linear_sq_resid.pdf', width = 8, height = 4,device=cairo_pdf)除此之外，作为一条很好的经验法则，除非你有非常充分的理由，否则 几乎从来不要用 到三次以上的多项式（\(cubic\)，即三次多项式 ）。到那个时候（用高阶多项式时 ），你就得开始担心，项数太多引发的统计问题会成为麻烦。用多项式（来解决非线性问题 ）这种 “疗法”，可能比非线性本身这个 “疾病” 更糟糕。
也有一些更自动化的方法来选择多项式的项数。你很容易就能找到有人建议：尝试不同数量的多项式项，然后去掉那些不显著的高阶项。换句话说，先试四次多项式\(（quartic）\)，要是\(X^4\)的系数不显著，就去掉它。然后试三次多项式，要是\(X^3\)项不显著…… 依此类推。我不推荐 这种方法，因为此时你是基于统计显著性检验来对模型做决策，而统计显著性检验 并非 用于干这个事儿的47 。
一种稍好的方法是使用像 最小绝对收缩和选择算子\(（LASSO, Least \space Absolute \space Shrinkage \space and \space Selection \space Operation ）\)这样的、有理论依据的模型选择算法。关于\(LASSO\)的具体内容，得等到本章后面 “惩罚回归\(（Penalized \space Regression）\)” 部分再讲。目前，你可以把它当成是普通最小二乘法的一个变体，它不仅试图最小化误差平方和，还在最小化的过程中对模型里的系数进行压缩。和因为不显著就去掉多项式项类似，它采取的思路是：那些对预测提升不大的额外多项式项，可以去掉。不过，它这么做的方式，是 专门 为判断一个额外的项是否 “值得纳入” 而设计的，而非像显著性检验那样，只是偶然地朝着类似目标去尝试。
13.2.3 变量转换
多项式项并非拟合非线性形状的唯一方式 。很多时候，直接对变量本身进行转换\(（transform）\)更有意义。
通常，这意味着在使用变量前，先让它通过一个函数。比如，你可以不进行\(Y = \beta_0 + \beta_1 X + \varepsilon\)这样的回归，而是进行：
\[ Y = \beta_0 + \beta_1 \ln(X) + \varepsilon \quad (13.10) \]这里的\(\ln()\)是自然对数（以\(e\)为底）48。
我们为啥要这么做呢？在使用变量前对其进行转换，主要有两个原因，其中一些我们在\(第3章\)已经讲过49 。
第一个原因 可能是为了让变量具备更契合回归分析的统计性质。比如，你可能想减少偏态\(（skew）\)。如果一个变量有很多低值观测，却只有少量极高值观测，那它就是”右偏”的50。收入就是偏态变量的典型例子。大多数人的收入接近中位数，但像首席执行官和电影明星这类超高收入者，挣的钱比中位数多得多51 。
偏态数据容易产生异常值\(（outlier）\)。异常值指的是取值与其他观测差异极大的观测。比如，要是有\(100\)个人，每年收入都差不多\(4\)万美元，而有一个人年收入\(100\)万美元，那这个人就是异常值52。
异常值\(（outlier）\)：异常值指的是取值与其他观测值相比，要么大得多、要么小得多的观测数据 。
异常值观测本身并非内在就有问题。有时，人们会把大的异常值当作数据记录有误的迹象，觉得数十万美元的薪资更可能是输入错误（本应是十万美元 ），但数据往往是没问题的，只是个异常值而已。然而，即便数据准确，异常值的存在也会让回归模型在统计上变得不稳定。因为这些观测值偏离程度大，往往会产生大的误差，所以仅少数几个观测值就能对回归线产生极大的拉扯作用，带来过度的影响。
使用一种能把大的观测值 “压缩” 得更贴近其他值的转换方法，就能减少这种偏态，让回归模型表现得更好一点。这在\(图13.10\)中可以看到。左侧是原始数据，有一个单独的大异常值远远高于其他数据。我对这些数据进行了两次\(OLS\)估计，一次剔除异常值，一次纳入异常值。估计出的\(OLS\)斜率变化极大！
# Outliers
set.seed(1000)
df <- data.frame(X = c(runif(30), 6)) %>% mutate(Y = 1.5*X + 2 + .5*rnorm(31))
df[31,2] <- df[31,2] +4
lin1 <- lm(Y~X, data = df[1:30,])
lin2 <- lm(Y~X, data = df)
log1 <- lm(log(Y)~log(X), dat = df[1:30,])
log2 <- lm(log(Y)~log(X), dat = df)
p1 <- ggplot(df, aes(x=X, y=Y)) +
geom_point() +
geom_smooth(formula = y~x, color='black', method='lm', se=FALSE) +
geom_line(aes(x=X, y=predict(lin1, newdata = df)), color = 'black', linetype = 'dashed') +
annotate(geom = 'text', label = 'OLS fit\nwithout outlier', hjust = .5, x = 5, y = 5, size = 10/.pt, family = 'sans') +
annotate(geom = 'text', label = 'OLS fit\nwith outlier', hjust = .5, x = 4, y = 14, size = 10/.pt, family = 'sans') +
labs(title = 'Linear Regression Without Transformation') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks.y = element_blank())
p2 <- ggplot(df, aes(x = log(X), y = log(Y))) +
geom_point() +
geom_smooth(formula = y~x, color = 'black', method = 'lm', se = FALSE) +
geom_line(aes(x = log(X), y = predict(log1, newdata = df)), color = 'black', linetype = 'dashed') +
annotate(geom = 'text', label = 'OLS fit\nwithout outlier', hjust = .5, x = 1, y = 1, size = 10/.pt, family = 'sans') +
annotate(geom = 'text', label = 'OLS fit\nwith outlier', hjust = .5, x = .25, y = 1.75, size = 10/.pt, family = 'sans') +
labs(title = 'Linear Regression With Transformation') +
theme_pubr() +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks.y = element_blank())
plot_grid(p1,p2)
ggsave('Regression/log_outliers.pdf', width = 8, height = 4,device=cairo_pdf)在右侧图中，我对\(X\)和\(Y\)都取了对数。你能立刻看到，异常值更贴近其他数据了（在图上表现为其他数据 “展开” 了，因为我们可以在不丢失异常值的情况下，更近距离地聚焦观察它们 ）。观测值在高端的间距被缩小了。而且，纳入或剔除异常值，对回归线的外观产生的影响也更微弱了。
要记住，异常值对\(OLS\)直线的斜率产生影响是正常的。实际上，本就该如此 53。所有观测值都该对回归线有某种影响。左侧图的问题在于，单个观测值的影响太大 了，而通过取对数，我们可以稍微缓解这种情况。
转换变量的第二个原因是尝试让变量之间呈现线性关系。记住，\(OLS\)天然适合拟合直线。只有对变量进行转换（无论是用多项式，还是其他转换方式 ），它才能拟合曲线。
这取决于你认为自己在建模的是何种关系。\(X\)增加\(1\)个单位，是否对应\(Y\)增加\(\beta_1\)个单位？如果是，那你处理的是线性关系，无需转换。但要是\(X\)增加\(1\)个单位，对应\(Y\)以某个百分比增加；或者\(X\)以某个百分比增加，对应\(Y\)有一定量的增加；又或者是类似的情况，该怎么办呢？
比如，要是我们用\(10\)年前的投资利率，对当前财富进行回归，真实关系是\(\text{Wealth} = \text{InitialWealth} \times e^{10 \times \text{InterestRate}}\) 。利率增加的效应是乘法层面的 ，而非线性的。所以我们需要进行转换，让关系 “线性化”。这种情况下，对数又能派上用场了。如果我们对两边取对数，得到：
\[ ln(Wealth) = ln(InitialWealth) + 10 \times InterestRate \quad (13.11) \]
这就是线性关系了54。\(OLS\)处理这种关系会很轻松。对数有助于把这类乘法关系转化为加法线性关系，因为\(log(XY) = log(X) + log(Y)\) 。
所以，要是你有一个理论上的关系，你得思考这两个变量是如何关联的 。这种关系应该是线性的吗？还是说，你得进行某种转换，才能让关系呈现出线性？
那有哪些类型的转换呢？
我们已经知道转换的原因了，但具体能进行哪些转换呢？
可选的转换有很多，我不打算逐一讲解。不过，有几种在社会科学研究中经常出现。
第一种，不出所料，是对数转换 。我们已经讲过这个。
第二、第三和第四种转换，源于对数转换存在的问题。具体来说，对数无法处理\(0\)值，因为\(log(0)\)无定义。那要是我们既想要对数处理偏态的良好特性，数据中又有 \(0\)值，该怎么办呢？
有三种常用的解决办法。对数转换之后，第二种转换是，在取对数前给变量加 1 ，即\(log(x+1)\)。这种方法能解决问题，但也有点凑合，还会让对数原本的一些优点打折扣。别深究细节，我 不推荐 这么做。
第三种转换是平方根转换 。和对数类似，对变量取平方根能 “压缩” 大的正异常值。而且，它能处理\(0\)值，因为\(\sqrt{0} = 0\) 。不过，虽然平方根能抑制大的观测值，但效果远不如对数。比如，从\(4\)到\(10000\)，平方根会乘以\(50\)（\(\sqrt{10000} / \sqrt{4} = 100 / 2 = 50\) ），而对数只乘以约\(6.6\)（\(\log(10000) / \log(4) \approx 9.2 / 1.4 \approx 6.6\) ）。另外，你还会失去对数特有的、很有用的 “百分比变化” 解读方式，我们马上会讲到这一点。
第四种是反双曲正弦变换\(（inverse \space hyperbolic \space sine \space transformation）\) ，也叫 “\(asinh\)”，等于\(\ln(x + \sqrt{x^2 + 1})\) 。这种变换有个好特性，既接受\(0\)值（\(\text{asinh}(0) = 0\) ），也接受负数 55。而且，对于足够大的值（比如大于\(10\)），它的结果和普通对数非常相似。所以，我们喜欢对数的那些优点，反双曲正弦变换也能提供。对于有\(0\)值、且我们的唯一目标是减轻偏态影响的偏态数据，反双曲正弦变换通常是我的推荐。
不过，如果目标是从变量中得到某种理论上的 “百分比变化” 解读，那么任何为处理\(0\)值而临时想的办法（本质上，\(0\)值和百分比变化不太搭 ），都会得出奇怪的结果。有一些方法可以减轻这个问题，比如通过缩放让变量的均值变大（比如用美元而非千美元 ）（Bellemare and Wichman 2020 ） 。你或许也可以用天然就能处理\(0\)值的方法，而不是用\(\text{asinh}(x)\)或\(\ln(1 + x)\) ，比如泊松回归，或者 “两部分” 模型 —— 把 “是\(0\)还是大于\(0\)”，以及 “大于\(0\)的值是多少” 当作两个独立问题来处理（Mullahy and Norton 2024 ） 。
我们要讲的第五种变换，它其实不是处理偏态或让关系线性化的，而是处理异常值的。缩尾（Winsorizing） 在金融学中特别流行，就是拿数据来，强行把两端 “挤一挤”。你只需找出那些离中心足够远的值，把它们向中心方向调整。要对数据进行\(X\%\)的缩尾，就把所有高于第X百分位数的观测值，替换成第X百分位数的值。
比如，要是我们有\(100\)个观测值，取值从\(1\)到\(100\)，对顶部和底部各\(5\%\)进行缩尾，那就意味着让\(6-95\)号观测值保持不变，然后把\(1,2,3,4,5\)号观测值替换成\(6\)，把 \(96,97,98,99,100\)号观测值替换成\(95\)。
这是处理异常值的一种非常简单粗暴的方法，但往往有效。它让大部分数据保持原样，要是你认为真实关系确实是线性的，但又不想让模型受异常值过度影响，这就很好。
我们要讲的第六种也是最后一种变换，完全不会改善模型的统计性质 。它不处理偏态、异常值，也不用于线性化。这就是标准化\(（standardizing）\)变换 —— 用变量减去其均值，再除以标准差，即\((X - \text{mean}(X)) / \text{sd}(X)\) 56。
既然它对模型统计性质没帮助，为啥还要这么做呢？因为它能让模型更易解读。我们知道，在模型\(Y = \beta_0 + \beta_1 X + \varepsilon\)中，对\(\beta_1\)的解读基于 “\(X\)增加\(1\)个单位”。但这意味着我们得时刻记住所有变量的尺度，才知道\(1\)个单位的增加有多大。但如果先进行标准化，\(\beta_1\)的解读就变成基于 “\(X\)增加\(1\)个标准差”，某些情况下，这样更容易理解。
对\(Y\)进行标准化也是如此。假设\(Y\)是学生考试分数，\(X\)是每天看电视的时间。如果\(\hat{\beta}_1 = -5\)，那么每天看电视时间增加\(1\)个单位，考试分数会降低\(5\)分。这是大的影响，还是小的影响呢？我们得对考试分数本身有更多了解才行。但如果先对\(Y\)标准化，\(\hat{\beta}_1 = -0.5\)，那就意味着每天多看\(1\)小时电视，考试分数会降低\(0.5\)个标准差 —— 这影响就挺大的！现在我们能立刻解读结果了。当然，这只对那些我们对其标准差的感知，比对单位的感知更清晰的变量有效。
所有这些变换的一个共同点是它们都是连续变换\(（continuous \space transformations）\) 。原始变量的微小变化，会转化为变换后版本的微小不同变化。不过，我们可能会认为，对某个变量来说，相关的变换是不连续的\(（discontinuous）\) 。比如，假设我们想研究年龄对酒驾可能性的影响。在美国，法定饮酒年龄是\(21\)岁。所以，思考年龄对酒驾影响的相关方式，可能不是 “年龄增加\(1\)岁的影响”，而是 “从\(21\)岁以下到\(21\)岁以上的影响”。从\(20\)岁\(11\)个月到\(21\)岁 \(0\)个月，不是平滑的变化，而是有一个急剧的转变。
我们该如何将其付诸实践呢？嗯，基于我们认为事物会急剧变化的相关节点，我们可以创建一组分类指示变量。如果我们想在这里做到完全灵活，我们可以为年龄的每一个可能取值创建一个二元指示变量（这被称为 “饱和\(（staurating）\)” 模型 ）57。这会让酒驾和年龄之间的关系呈现出几乎任意的形状，让不连续性自行显现。但另一方面，我们的数据中可能不会有很多人正好是，比如说，\(22\)岁\(3\)个月大，所以我们的估计不会很精确。
我们可以通过提前确定我们认为的相关不连续性，来获得更强的分析效力。我们可能会说 “满\(21\)岁是相关的不连续性，所以我要创建一个分类指示变量，区分年龄在 \(21\)岁以下和以上，然后把这个变量纳入我的模型”。或者，我们可以把这个二元指示变量和一个交互项（本章后面会讲到 ）结合起来，说 “年龄对酒驾有线性影响，但在\(21\)岁时会有一个不连续的跳跃”。如果你的变量存在特殊的、有意义的取值，那么直接对其建模真的有助于准确呈现数据！
认真思考变量中的不连续性很重要，如果你知道存在不连续性却不加以考虑，可能会损害你的模型。以 “扎堆\(（bunching）\)” 现象为例。“扎堆” 指的是，即便你的变量是连续的，但存在一个特别常见的取值。通常这是一个 “临界值”。比如，你每周工作多少小时？合理的取值可能是\(12,13,42,16.25\)等等。但在给定的数据集中，你可能会发现很多人正好工作\(40\)小时（或者在你们国家，不管多少小时算全职 ），还有很多没有工作的人，工作时长正好是\(0\)小时。
这种 “扎堆” 可能会隐藏 “后门（\(back \space doors\)，可理解为影响模型的潜在因素 ）”！工作\(0\)小时的人和工作\(1\)小时的人之间，可能存在本质上的、不连续且无法观测的差异，而这种差异在工作\(1\)小时和工作\(2\)小时的人之间是不存在的。所以，你的数据中存在大量工作时长正好为\(0\)小时的人，可能意味着存在一些 “后门”，而这些 “后门” 可能是敞开的 58！
在存在 “扎堆” 现象的情况下，有一个相当简单的变换方法，可用于检验你是否真的需要担心这些 “敞开的后门”。只需采用你原本要用的任意变换，再添加一个二元指示变量，用于表示是否恰好处于 “扎堆点”（比如，在这个例子中，就是 “工作时长恰好为\(0\)小时” 的指示变量 ）。要是该系数与\(0\)相差较大（即你能拒绝系数为\(0\)的假设 ），那么这种 “扎堆” 现象很可能让一些 “后门” 处于敞开状态（Caetano et al. 2021 ） 。
就这些了 —— 这就是我们要讲的所有变换。世界上还有许许多多其他变换，但以我的经验，这些是社会科学研究中最常见的。说实话，我不确定自己除了这七种变换或者多项式变换外，还有没有用过别的变换。
话虽如此，仅仅列出这些变换，并不意味着我们对它们的探讨已至终点。尤其是对数变换 —— 我似乎一直讲个没完 —— 会让我们的结果解读起来更棘手。那该如何解读呢？
一开始看似没问题。接着你仔细想想，就会觉得这是个大难题。然后你对它了解得更多，又会发现它根本不是问题。生活就是这样 59。
咱们先看对\(X\)进行对数变换的回归（之后再讲对\(Y\)的变换 ）：
\[ Y = \beta_0 + \beta_1 ln(X) + \varepsilon \quad (13.12) \]
一开始看起来怎样？没问题呀！我们可以把\(\hat{\beta}_1\)解读为 “\(\ln(X)\)变化 1 个单位，\(Y\)会有\(\hat{\beta}_1\)个单位的变化” 。很好。
现在仔细想想 ——\(\ln(X)\)变化 \(1\)个单位，到底是什么意思呢？这和\(X\)变化 \(1\)个单位完全不是一回事。
咱们先看对\(X\)进行对数变换的回归（之后再讲对\(Y\)的变换 ）。
\[ Y = \beta_0 + \beta_1 ln(X) + \varepsilon \quad (13.12) \] 一开始看起来咋样？没问题呀！我们可以把\(\hat{\beta}_1\)解读成 “\(\ln(X)\)变化 1 个单位，Y会有\(\hat{\beta}_1\)个单位的变化” 。很好。
现在仔细琢磨一下 ——\(\ln(X)\)变化 1 个单位，到底是 什么意思呢？这和X变化 1 个单位完全不是一回事。
咱们深入了解一下，之后就又会觉得简单了。\(\ln(X)\)的线性变化，和\(X\)的 比例变化 相关。好好想想，我们知道\(\ln(AB) = \ln(A) + \ln(B)\) 。所以，如果从\(\ln(X)\)开始，给它 加1 ，就会得到\(1 + \ln(X) = \ln(e) + \ln(X) = \ln(eX)\) 60。给\(\ln(X)\)加上某个数，等价于给\(X\)乘以某个数。
具体来说，不管我们给\(\ln(X)\)加上什么数，把这个要加的数叫\(c\)，实际上就是把\(X\)乘以\(e^c\) 。现在我们可以利用一个实用的小结论：如果\(c\)接近 0，那么\(e^c\)接近\(1 + c\) 。为啥实用呢？因为这意味着\(\ln(X)\)增加\(0.05\)，大致等价于把\(X\)乘以\(1.05\)，换句话说，\(\ln(X)\)增加\(0.05\)，大致等价于\(X\)增加\(5\%\) 61。现在你知道我之前为啥说本质上成比例的关系了吧62。
咱们把这个直观认识用起来，这样不管对数在方程里的哪个位置，我们都能解读。
首先，回到对\(X\)进行对数变换的原始方程：\(Y = \beta_0 + \beta_1 \ln(X) + \varepsilon\) 。
这里咱们用\(X\)变化\(10\%\)，而非\(100\%\)，这样更贴近实际（对大多数社会科学研究而言，\(10\%\) 的变化比 \(100\%\) 的变化常见得多，而且\(1 + c \approx e^c\)这个近似也更准 ）。\(X\)增加 \(10\%\)，即\(X\)变为原来的\(1.1\)倍，对应\(\ln(X)\)的变化是\(\ln(1.1) \approx 0.1\)个单位，这会让\(Y\)变化\(0.1 \times \beta_1\)个单位。所以，我们可以这样解读\(\beta_1\)：\(X\)增加\(10\%\)，与\(Y\)变化\(0.1 \times \beta_1\)个单位相关。
接下来，试试把对数用到\(Y\)上：\(\ln(Y) = \beta_0 + \beta_1 X + \varepsilon\) 。现在，\(X\)变化 1 个单位，与\(\ln(Y)\)变化\(\beta_1\)个单位相关，这意味着\(Y\)变化\(\beta_1 \times 100\)（只要\(\beta_1\)足够小，能让我们用近似的话 ）。
最后，把对数用到两边：\(\ln(Y) = \beta_0 + \beta_1 \ln(X) + \varepsilon\) 。现在，\(X\)增加\(10\%\)，对应\(\ln(X)\)变化约\(0.1\)个单位，这会让\(\ln(Y)\)变化\(0.1 \times \beta_1\)个单位，转化为\(Y\)变化\(0.1 \times \beta_1 \times 100\)，也就是\(10 \times \beta_1\) 63。
就这些！你要是愿意，可以记住这三种解读方式，但更妙的办法是每次都梳理逻辑，这样你就知道背后的原理。另外，作为本书附赠内容，思考任何变换时，不是所有变换都有百分比解读，但你总能思考 “变换后的变量变化\(1\)个单位，对应原始变量是什么意思”，并据此得出解读。
13.2.4 交互项
多项式和变换能让\(Y\)与\(X\)之间的关系更灵活，这很棒。我们可以拟合出各种曲线，让\(Y\)和\(X\)的关系随\(X\)取值变化。但仍有遗漏：要是\(Y\)与\(X\)的关系，不是基于\(X\)自身取值，而是基于另一个变量\(Z\)的取值 ，该咋整呢？
比如，汽油价格和个人驾车里程数有啥关系？对有车的人来说，关系可能很强且为负；对没车的人，关系可能接近\(0\)；对有车但常骑车出行的人，关系可能很弱。汽油价格和驾车里程的关系，因 是否拥有汽车 而异。我们或许想让这种关系可变化，以便恰当建模，也可能想探究不同群体间关系是否有差异 64。
但用多项式或变换，无法体现这种关系的变化。这时，我们需要用 交互项\(（interaction \space terms）\) 。
交互项，就是把两个不同变量相乘，将乘积纳入回归，比如\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \beta_3 XZ + \varepsilon\) 。
交互项\(（interaction \space terms）\)：两个变量相乘的乘积，纳入回归模型，通常和未相乘的变量版本一起纳入。
构建模型时，操作上就这么简单。要是你觉得\(Y\)与\(X\)的关系因\(Z\)取值不同而有差异，就把\(X \times Z\)和\(Z\)本身都纳入模型65 。
交互项的真正难点在于解读 。解读起来有点烧脑，挺恼人的；在标准社会科学研究设计里，没走两步就会碰到交互项是核心的研究，躲都躲不开。
解读含交互项模型的问题，实则包含两个小问题。第一个是：当模型中\(X\)和其他变量存在交互时，\(X\)的效应是怎样的？第二个是：我该如何解读交互项本身？
这两个问题的直观理解，都要回归到 “微积分” 这个老朋友。要是你不懂微积分，也能跟上这个直观思路；要是你懂微积分，嘿，微积分超酷的，对吧？咱也算 “微积分同好会” 成员啦 66。
先解决第一个问题：当模型里\(X\)和其他变量有交互时，\(X\)的效应是啥？要解决这个，写出回归方程，然后对\(X\)求导就行。搞定！毕竟，\(Y\)对\(X\)的导数，就是\(X\)变化时\(Y\)的变化情况，这和回归系数的解读是一回事。要是你不懂微积分，这导数也很好求。只要\(X\)没有多项式或变换项，就把所有含\(X\)的项找出来，然后把\(X\)去掉，就成。要是\(X\)有多项式或变换项…… 网上有导数计算器，能帮你算。
所以，要是我们的回归方程是：\[ Y = \beta_0 + \beta_1 X + \beta_2 Z + \beta_3 XZ + \varepsilon \quad (13.13) \]那么，我们求导可得：
\[ \frac{\partial Y}{\partial X} = \beta_1 + \beta_3 Z \quad (13.14) \]
这就是\(X\)对\(Y\)的效应 67。
当然，这还没完，因为我们算出的\(X\)的效应里包含\(Z\)。这意味着，\(X\)没有单一的效应 68。相反，我们得问，在\(Z\)取特定值 时，\(X\)的效应是多少。我们可以说，当\(Z = 0\)时，\(X\)的效应就是\(\beta_1\)；当\(Z = 16\)时，效应是\(\beta_1 + 16\beta_3\)；其他取值以此类推。
这也意味着，单个系数本身没法告诉我们太多信息，就像多项式回归里那样。\(\beta_1\)不再代表 “\(X\)的效应” 了。所以，要是你关注效应大小，就得考虑整个\(X\)的效应表达式。只看\(\beta_1\)，能得到的信息非常少（而且，\(\beta_1\)的显著性，能告诉你的信息也非常少 ）。
这就是当\(X\)涉及交互项时，我们如何得到它的效应。那另一个问题呢？交互项本身该怎么解读？
从形式上看，同样用求导的方法也能解决 —— 在回归方程\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \beta_3 XZ\)中，交互项\(\beta_3\)是\(Y\)对\(X\)和\(Z\)的交叉导数\(（cross-derivative）\)，也就是\(\frac{\partial^2 Y}{\partial X \partial Z}\) 。不过这理解起来有点难，所以我们可以退一步，用更直观、不那么数学的方式来思考。
解读交互项，有个不错的思路：在回归方程\(Y = \beta_0 + \beta_1 Z\)里，\(Z\)的系数能告诉我们，\(Z\)增加\(1\)个单位时，对\(Y\)的预测值怎么变，对吧？但咱们已经知道，在方程\(Y = \beta_0 + \beta_1 X + \beta_2 Z + \beta_3 XZ\)中，求\(X\)的效应，要么求导，要么把含\(X\)的项找出来、去掉\(X\)，得到\(\frac{\partial Y}{\partial X} = \beta_1 + \beta_3 Z\) 。这看着眼熟，是吧？现在，在这个方程里，\(Z\)的系数\(\beta_3\)，能告诉我们：\(Z\)增加\(1\)个单位时，对\(\frac{\partial Y}{\partial X}\)（也就是\(X\)对\(Y\)的效应 ）的预测值怎么变。换句话说，\(\beta_3\)表示，\(Z\)增加\(1\)个单位时，\(X\)对\(Y\)的效应会增强多少 69。
通常，参与交互的变量是二元变量（比如\(0-1\)变量 ）。这不会改变解读方式，但能让计算更简单。假设我们有回归方程\(Y = \beta_0 + \beta_1 X + \beta_2 \text{Child} + \beta_3 X \times \text{Child}\) ，其中\(\text{Child}\)是二元变量，孩子取\(1\)，非孩子取\(0\)。从这个方程能得到啥？和之前一样，\(X\)对\(Y\)的效应是\(\beta_1 + \beta_3 \text{Child}\) ，而\(\beta_3\)表示，\(\text{Child}\)增加\(1\)个单位时，\(X\)对\(Y\)的效应会增强多少70。
但\(\text{Child}\)只能从\(0\)变到\(1\)—— 它就这两个取值。所以，要是\(X\)对\(Y\)的效应是\(\beta_1 + \beta_3 \text{Child}\) ，那就意味着：当\(\text{Child} = 0\)（非孩子 ）时，\(X\)对\(Y\)的效应是\(\beta_1\)（代入\(\text{Child} = 0\)，得到\(\beta_1 + \beta_3 \times 0 = \beta_1\) ）；当\(\text{Child} = 1\)（孩子 ）时，效应是\(\beta_1 + \beta_3 \times 1 = \beta_1 + \beta_3\) 。
这两种效应的差值 —— 孩子组的\(\beta_1 + \beta_3\) 与非孩子组的\(\beta_1\)之差 —— 就是\(\beta_3\) 。这代表孩子和非孩子之间，\(X\)对\(Y\)的效应差异。换个说法，当\(\text{Child}\)从 \(0\)（非孩子 ）增加到 \(1\)（孩子 ）时，\(X\)对\(Y\)的效应增强了多少，就是\(\beta_3\) 。和之前的解读一致。要是我们发现\(\beta_3\)的系数有实际意义（或者你更愿意说 “统计显著” ），那就可以说 “\(X\)对\(Y\)的效应在孩子和非孩子之间存在差异”。
res <- read_csv('Regression/restaurant_inspections.csv', show_col_types = FALSE)
# interaction table regressions
inx <- res %>%
group_by(business_name) %>%
mutate(NumberofLocations = n(),
Weekend = lubridate::wday(lubridate::mdy(inspection_date)) %in% c(1,7)) %>%
mutate(Year = year(mdy(inspection_date))) %>%
select(inspection_score, Year, NumberofLocations, Weekend) %>%
dplyr::filter(!is.na(inspection_score)) %>%
rename(`Inspection Score` = inspection_score,
`Number of Locations` = NumberofLocations,
`Year of Inspection` = Year)## Adding missing grouping variables: `business_name`m1 <- lm(`Inspection Score`~`Number of Locations` + `Year of Inspection` + Weekend, data = inx)
m2 <- lm(`Inspection Score`~`Number of Locations` + `Year of Inspection` + Weekend + `Number of Locations`:`Year of Inspection`, data = inx)
m3 <- lm(`Inspection Score`~`Number of Locations` + `Year of Inspection` + Weekend + `Number of Locations`:Weekend, data = inx)
knitr::opts_current$set(label = 'statisticaladjustment-restaurantregsinx')## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.msummary(list('Inspection Score' = m1,'Inspection Score' = m2, 'Inspection Score' = m3),
coef_map = c('`Number of Locations`' = 'Number of Locations',
'WeekendTRUE' = 'Weekend',
'`Year of Inspection`' = 'Year of Inspection',
'`Number of Locations`:`Year of Inspection`' = 'Number of Locations x Year of Inspection',
'`Number of Locations`:WeekendTRUE' = 'Number of Locations x Weekend'),
stars = TRUE,
gof_omit = 'AIC|BIC|Lik',
output = 'markdown')| Inspection Score | Inspection Score | Inspection Score | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| Weekend | 1.432*** | 1.459*** | 1.759*** |
| (0.419) | (0.418) | (0.488) | |
| Num.Obs. | 27178 | 27178 | 27178 |
| R2 | 0.069 | 0.072 | 0.069 |
| R2 Adj. | 0.069 | 0.072 | 0.069 |
| F | 669.086 | 525.621 | 502.255 |
| RMSE | 6.04 | 6.03 | 6.04 |
# For code examples
res %>%
group_by(business_name) %>%
mutate(NumberofLocations = n(),
Weekend = lubridate::wday(lubridate::mdy(inspection_date)) %in% c(1,7)
) %>%
mutate(Year = year(mdy(inspection_date))) %>%
select(inspection_score, Year, NumberofLocations, Weekend) %>%
dplyr::filter(!is.na(inspection_score)) %>%
write_csv('Regression/restaurant_data.csv')## Adding missing grouping variables: `business_name`在\(表13.7\)中，我们能看到本章前面\(表13.2\)的延续内容，不过这次纳入了一些交互项。第二列里有”门店数量\(（NumberOfLocations）\)” 和 “检查年份\(（YearofInspection）\)” 的交互项；第三列中，有 “门店数量” 和”周末\(（Weekend）\)“的交互项，”周末”是个二元变量，用来表示检查是否在周末进行。
看第二列，“门店数量” 的效应是\(-1.466\)（即 “门店数量” 的系数 ）加上\(0.001\times检查年份\)（交互项的系数 ）。咋解读这个呢？单独看\(-1.466\)说明不了啥 —— 这代表的是，当我们代入 “检查年份\(= 0\)”（也就是公元\(0\)年！）时71，“门店数量” 增加\(1\)个单位和 “检查分数” 之间的关系，远远超出了我们的数据范围，反映不出实际预测情况 。但我们可以代入数据里有的值，比如 “检查年份\(= 2008\)”，这样就能说，\(2008\)年时，“门店数量” 增加\(1\)个单位和 “检查分数” 的关系是\(-1.466 + 0.001×2008 = 0.542\)，是正的 ！实际上，在数据涵盖的所有年份里，这个效应都是正的，尽管有个很大的负数\(-1.466\)。别被单独的系数骗了，解读时一定要结合相关交互项。我们还能说，每过一年，“门店数量” 增加\(1\)个单位和 “检查分数” 的关系，会更正向地变化\(0.001\)。
在第三列中，“门店数量\(（NumberOfLocations）\)” 的效应是\(-0.019 - 0.010×周末（Weekend）\) 。由此，我们可以说：在非周末时，“门店数量” 增加\(1\)个单位和 “检查分数” 的关联是\(-0.019\)；在周末时，是\(-0.019 - 0.010 = -0.029\) 72。我们还能说，周末时，这种关系比工作日时更负面 \(0.010\) 73。
使用交互项时要牢记的几点
第一，得非常仔细地思考 为啥 要纳入某个交互项，因为要是盲目 “钓鱼”（指无明确理论驱动地尝试各种交互项），容易得到侥幸的结果。回到\(第1章\)想想 —— 要是任由数据 “随意发挥”，我们会得到没多大意义的结果。做研究、分析数据，得有特定的设计，要有 强有力的理论依据 ，能说明为啥预期\(X\)的效应会随\(Z\)的不同取值而变化，之后再添加\(X×Z\)交互项。“我也不知道，可能效应不一样吧”，这可不是个好理由。
无论有没有交互项，带着设计和意图去分析数据，总归是好的。为啥讲到交互项时，我要这么强调呢？首先，因为你 能尝试的交互项太多了，很多看着还挺诱人。比如，工作培训的效应在不同性别间有差异吗？不同种族间？不同年龄段间？确实，有可能。我们能给自己编个好听的故事，解释为啥这些变量和工作培训做交互会有意义。但要是把所有交互都试一遍，即便实际啥都没有，也很可能很快就发现 “显著的交互项”。所以，盲目尝试交互项，容易导致 假阳性结果 。
当发现 “无效应” 时，人们尤其难抗拒 “啥都试试” 的冲动。比如，也许我们没发现工作培训有啥效应，但说不定这效应只在女性中存在！那就试试性别交互项，看看即便\(\beta_1\)不显著，\(\beta_1 + \beta_3\)会不会很大、很显著。或者，效应只在太平洋岛民中存在！试试种族交互项。又或者，只在艾奥瓦州存在，那就试试州交互项。然而，因为可供我们研究的子群体几乎无穷多 （尤其是深入挖掘时 —— 比如，说不定效应只在艾奥瓦州的女性太平洋岛民中存在 ），要是把每个子群体都检查一遍，几乎肯定能发现一些 “效应”。总的来说，对于 “全样本中根本不存在、但在某些子群体中存在” 的效应，你得高度怀疑。当然，确实有很多效应，实际上只在特定子群体中存在。但盲目用一堆交互项 “钓鱼”，会产生 大量假阳性结果 。所以，只在有 强有力理论依据 、预期效应仅存在于某子群体或在群体间有差异时，才用交互项；而且，在检验数据前，得先相信自己的理论 “故事”。
第二，即便我们坚信可能存在效应差异，交互项本身也很 “噪声”（\(noisy\)，指估计值不稳定 ）。本质上，交互项是在 找两个 “有噪声” 事物的差异 —— 两个群体的效应，而这两个效应本身都有抽样变异74。两个 “有噪声” 事物的差异，会更 “噪声”75。要精确估计交互项，比精确估计全样本效应，需要多得多的观测值 。即便差异大到一个群体的效应是另一个的两倍，估计交互项所需的样本量，也可能是主效应的\(16\)倍（Gelman 2018 ） ，这样才能有足够的统计效力。
为啥 “噪声” 是个问题？因为 “噪声” 的估计值，更常出现 “大得惊人” 的效应。抽样变异大，估计值的波动范围就宽。所以，当你真的发现一个 “大到值得关注” 的交互项时，很可能只是抽样变异，即便结果在统计上显著！
话虽如此，要是你有充分理由纳入交互项，那交互项仍然非常值得使用。而且，对于因果推断来说，学习交互项至关重要。实际上，本书后半部分的很多研究设计，都是围绕交互项展开的。事实证明，通过探究处理变量的效应在不同条件下的差异（比如，在理应存在因果效应的条件，和仅存在后门路径效应的条件之间的差异 ），能在识别因果关系上发挥很大作用。所以，花些时间练习解读和使用交互项吧，肯定会有回报的。
13.2.5 编写多项式和交互项的代码
在这两个代码块中，我们将继续使用本章一直在处理的餐厅检查数据。
首先，我们要做一个包含多项式的回归。然后，针对该变量的某个给定取值，计算 “门店数量（NumberofLocations）” 的效应，同时算出该取值下效应的标准误。
# Load packages and data
library(tidyverse); library(modelsummary); library(marginaleffects);library(causaldata)
df <- causaldata::restaurant_inspections
# Use I() to add calculations of variables
# Including squares
m1 <- lm(inspection_score ~ NumberofLocations + I(NumberofLocations^2) + Year, data = df)
# Print the results to a table
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 362.834*** |
| (11.601) | |
| NumberofLocations | -0.084*** |
| (0.001) | |
| I(NumberofLocations^2) | 0.000*** |
| (0.000) | |
| Year | -0.132*** |
| (0.006) | |
| Num.Obs. | 27178 |
| R2 | 0.210 |
| R2 Adj. | 0.210 |
| AIC | 170415.0 |
| BIC | 170456.1 |
| Log.Lik. | -85202.504 |
| F | 2402.353 |
| RMSE | 5.56 |
# Use avg_slopes to calculate the effect of locations at 100 locations.
# variables is the variable we want the effect of and newdata can be used with datagrid to get the effect at a particular value
avg_slopes(m1, variables = 'NumberofLocations', newdata = datagrid(NumberofLocations = 100))##
## Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## -0.0597 0.000706 -84.5 <0.001 Inf -0.0611 -0.0583
##
## Term: NumberofLocations
## Type: response
## Comparison: dY/dX# The Estimate is the AME (average marginal effect), which we want. The standard error is reported, too.import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# Use I() to insert calculations of your variables
# and ** to square
m1 = smf.ols(formula = '''inspection_score ~ NumberofLocations + I(NumberofLocations**2) + Year''', data = df).fit()
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.210 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.210 |
| Method: | Least Squares | F-statistic: | 2402. |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.00 |
| Time: | 10:57:59 | Log-Likelihood: | -85203. |
| No. Observations: | 27178 | AIC: | 1.704e+05 |
| Df Residuals: | 27174 | BIC: | 1.704e+05 |
| Df Model: | 3 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 362.8342 | 11.601 | 31.276 | 0.000 | 340.096 | 385.573 |
| NumberofLocations | -0.0844 | 0.001 | -82.888 | 0.000 | -0.086 | -0.082 |
| I(NumberofLocations ** 2) | 0.0001 | 1.77e-06 | 69.684 | 0.000 | 0.000 | 0.000 |
| Year | -0.1319 | 0.006 | -22.870 | 0.000 | -0.143 | -0.121 |
| Omnibus: | 2526.441 | Durbin-Watson: | 1.943 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 3411.908 |
| Skew: | -0.775 | Prob(JB): | 0.00 |
| Kurtosis: | 3.782 | Cond. No. | 1.70e+07 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.7e+07. This might indicate that there are
strong multicollinearity or other numerical problems.
# Use t_test to test linear combinations of coefficients be sure to test them against 0 to get the appropriate result
# coef is the estimate here, std err the SE.
# We know this is the right equation to use because we know the derivative - we need to figure that out first
m1.t_test('NumberofLocations + 2*I(NumberofLocations ** 2)*100 = 0')## <class 'statsmodels.stats.contrast.ContrastResults'>
## Test for Constraints
## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## c0 -0.0597 0.001 -84.505 0.000 -0.061 -0.058
## ==============================================================================causaldata restaurant_inspections.dta, use clear download
* ## interacts variables
* Interact a variable with itself to get a square!
* use c. to let Stata know it's continuous
reg inspection_score c.numberoflocations##c.numberoflocations year
* look at the results in a nicer table
esttab
* Use margins with the dydx() option to get the effect of numberoflocations
* and use at() to specify that you want it when there are 100 locations
margins, dydx(numberoflocations) at(numberoflocations = 100)
* The first term is the effect, second is the SE接下来，我们会纳入一个交互项。同样地，我们会在交互变量的特定取值下，计算门店数量（\(Locations\)，结合前文推测是 \(NumberOfLocations\) 这类代表门店数量的变量 ）的效应（记住，交互项本身就已经体现了效应之间的差异，我们无需通过单独的命令来获取这一差异 ） 。
# Load packages
library(tidyverse); library(modelsummary); library(marginaleffects); library(causaldata)
df <- causaldata::restaurant_inspections
# Use * to include two variables independently plus their interaction
# (: is interaction-only, we rarely use it)
m1 <- lm(inspection_score ~ NumberofLocations*Weekend +Year, data = df)
# Print the results to a table
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 225.126*** |
| (12.415) | |
| NumberofLocations | -0.019*** |
| (0.000) | |
| WeekendTRUE | 1.759*** |
| (0.488) | |
| Year | -0.065*** |
| (0.006) | |
| NumberofLocations × WeekendTRUE | -0.010 |
| (0.008) | |
| Num.Obs. | 27178 |
| R2 | 0.069 |
| R2 Adj. | 0.069 |
| AIC | 174871.9 |
| BIC | 174921.2 |
| Log.Lik. | -87429.963 |
| F | 502.255 |
| RMSE | 6.04 |
# Use avg_slopes to calculate the effect of locations on the weekend
# variables is the variable we want the effect of and newdata can be used with datagrid to get the effect at a particular value
avg_slopes(m1,variables = 'NumberofLocations', newdata = datagrid(Weekend = TRUE))##
## Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## -0.0289 0.00751 -3.85 <0.001 13.1 -0.0437 -0.0142
##
## Term: NumberofLocations
## Type: response
## Comparison: dY/dX# The Estimate is the AME (average marginal effect), which we want. The standard error is reported, too.import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# Use * to include two variables independently plus their interaction
# (: is interaction-only, we rarely use it)
m1 = smf.ols(formula = '''inspection_score ~ NumberofLocations*Weekend + Year''', data = df).fit()
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.069 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.069 |
| Method: | Least Squares | F-statistic: | 502.3 |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.00 |
| Time: | 10:58:00 | Log-Likelihood: | -87430. |
| No. Observations: | 27178 | AIC: | 1.749e+05 |
| Df Residuals: | 27173 | BIC: | 1.749e+05 |
| Df Model: | 4 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 225.1260 | 12.415 | 18.134 | 0.000 | 200.793 | 249.460 |
| Weekend[T.True] | 1.7592 | 0.488 | 3.606 | 0.000 | 0.803 | 2.715 |
| NumberofLocations | -0.0191 | 0.000 | -43.759 | 0.000 | -0.020 | -0.018 |
| NumberofLocations:Weekend[T.True] | -0.0098 | 0.008 | -1.307 | 0.191 | -0.025 | 0.005 |
| Year | -0.0648 | 0.006 | -10.494 | 0.000 | -0.077 | -0.053 |
| Omnibus: | 2829.404 | Durbin-Watson: | 1.931 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 3861.648 |
| Skew: | -0.849 | Prob(JB): | 0.00 |
| Kurtosis: | 3.725 | Cond. No. | 6.82e+05 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.82e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
# Use t_test to test linear combinations of coefficients be sure to test them against 0 to get the appropriate result
# coef is the estimate here, std err the SE.
# We got the proper coefficient names from reading the reg table
m1.t_test('NumberofLocations + NumberofLocations:Weekend[T.True] = 0')## <class 'statsmodels.stats.contrast.ContrastResults'>
## Test for Constraints
## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## c0 -0.0289 0.008 -3.851 0.000 -0.044 -0.014
## ==============================================================================causaldata restaurant_inspections.dta, use clear download
* ## interacts variables
* use c. to let Stata know locations is continuous
* and i. to let it know weekend is binary
reg inspection_score c.numberoflocations##i.weekend year
* look at the results in a nicer table
esttab
* Use margins with the dydx() option to get the effect of numberoflocations
* and use at() to specify that you want it when it's the weekend
margins, dydx(numberoflocations) at(weekend = 1)
* The first term is the effect, second is the SE13.2.6 非线性回归
到目前为止，尽管我们允许预测变量呈现出各种形式—二元的、有偏态的、与因变量呈非线性关系的—但我们对结果变量的限制更为严格。我们一直假定，我们可以使用预测变量的线性函数来预测因变量的条件均值。换句话说，我们可以通过获取预测变量，将每个预测变量乘以其系数，然后将所有结果相加，来预测因变量76。
然而，这并非总是可行。要用线性函数进行预测，事物的表现得像，嗯，一条直线。它得随着预测变量的变化而平滑变化（所以，比方说，它不能是一个离散值，一下跳上去，而非逐步攀升 ）。而且，它得能向任意方向无限延伸，也就是能取任何值 —— 要是你估计的模型是\(\text{NumberofChildren} = 4.5 - 0.1\text{IncomeinThousands}\) ，那么对于收入为\(50,000\)美元（即\(\text{IncomeinThousands} = 50\) ）的人，预测的孩子数量是\(4.5 - 0.1×50 = -0.5\) 个。这显然说不通 。
这就引出了非线性回归。这是一个内容广泛的研究领域，包含各种不依赖预测变量和结果变量之间线性关系的回归方法。这些方法通常会根据你所处理的因变量的具体类型进行调整。所以，在本章中（很遗憾 ），我们只能涵盖用于处理二元因变量的概率单位\(（probit）\)模型和对数单位\(（logit）\)模型。不过，还有用于计数变量的泊松\(（possion）\)回归、用于分类数据的多项对数\(（multinomial \space logit）\)和多项概率单位\(（multinomial \space probit）\)模型、用于删失数据的托宾\(（tobit）\)模型等等77。
即便对这些特殊的因变量使用线性回归并无意义，仍有人会这么做。实际上，即便因变量不太合适，使用普通最小二乘法也有一些合理的理由。所有这些非线性回归模型都有其自身的假设，而且通常没有解析解78，这会使得它们的估计值稳定性稍差。它们在处理诸如对类别众多的分类变量进行控制这类情况时可能会遇到困难。此外，当对二元结果变量使用\(OLS\)时，会导致误差项存在异方差性（我们会在几节之后讨论这个问题 ），而在对数单位\(（logit）\)和概率单位\(（probit）\)模型中，异方差性的存在实际上会使系数出现偏差，因此需要使用专门的异方差稳健型估计量。在二元因变量的情境下，不管怎样使用\(OLS\)的做法被称为“线性概率模型”，而且确实有相当多的应用 。
话虽如此，虽然在某些情况下，非线性回归这种”治疗方法”比线性模型用错的”病症”更糟，但一般来说并非如此，而且你会希望使用合适的模型。到底为啥这么说呢？正如我之前提到的，我会详细讲清楚我们为啥要这么做、它的原理，特别是在二元因变量的应用方面。
广义线性模型\(（GLMs）\)是一种将我们已掌握的回归知识，拓展到更广泛情境的方法，在这些情境中我们需要一定的非线性。广义线性模型不是进行非线性回归的唯一方式，但概率单位\(（probit）\)模型和对数单位\(（logit）\)模型常用这种方法，而且它也适用于其他一些情境。一旦我们学过了回归，理解广义线性模型就相当容易了。
广义线性模型：一种将线性指标通过特定函数转换生成预测值的线性模型。
广义线性模型\(（GLM）\)是这么说的：拿那个回归模型，把它代入一个函数，以此来进行预测79。我们把这个函数\(F()\)称作 “链接函数”，把输入到这个函数里的部分称作 “指数”。就这么简单！所以，我们的回归方程是：
\[Y = F(\beta_0 + \beta_1 X) \quad (13.15)\]
注意，函数\(F()\)里的\(\beta_0 + \beta_1 X\)这个指数，和我们之前用到的回归方程完全一样80。还要注意，要是我们让\(F(x) = x\)，那我们就又回到了\(Y = \beta_0 + \beta_1 X\)，基本上就是在做普通最小二乘法了。
广义线性模型\(（GLM）\)是如何让我们进行非线性回归的呢？只要我们选取一个能模拟\(Y\)所具备特性的函数，就能确保不会做出不切实际的预测。
对于二元变量，我们要预测\(Y\)等于\(1\)的概率。换句话说，\(F(\beta_0 + \beta_1 X)\)得给我们一个概率。概率应在\(0\)到\(1\)之间。那么，对于与二元\(Y\)适配的链接函数，我们要找什么样的呢？
它得能接收从\(-\infty\)到\(\infty\)的任意值（因为\(\beta_0 + \beta_1 X\)能取从\(-\infty\)到\(\infty\)的任意值 ）
它只能输出 \(0\) 到 \(1\)之间的值
只要输入（\(\beta_0 + \beta_1 X\) ）增加，输出也得增加
有无数个函数能满足这些条件，但最常用的有两个：\(logit\)链接函数和 \(probit\)链接函数 。
\(logit\)链接函数为：
\[ F(x) = \frac{e^x}{1+e^x} \quad (13.16) \]
\(probit\)链接函数是\(\Phi(x)\)，其中\(\Phi\)是标准正态分布的 “累积分布函数”，即\(\Phi(x)\)会给出\(X\)在标准正态分布（均值为\(0\)、标准差为\(1\)）中所处的百分位数。
\(probit\)和\(logit\)这两种链接函数，产生的回归预测结果几乎完全相同。通常无需太纠结用哪种81。但关键在于，二者都不会做出\(0\)到\(1\)范围之外的预测，而且在估计系数时，会考虑到这些边界，因此更适合用于这类非连续的因变量。
“和平、爱与线性”有啥问题？为啥不就用普通最小二乘法呢？我们真正感兴趣的是\(X\)的增加如何导致\(Y\)的变化。按理说，把识别问题解决好才是真正重要的步骤，而且只要还能得到那些效应，我们不一定在意预测出的\(Y\)值是多少，对吧？
这种观点存在两个问题。第一，任何时候，当你估计的模型不是正确的模型时，奇怪的结果可能会以意想不到的方式出现，而且是那种很难提前预测到的方式82。第二，没有妥善考虑非线性，可能会搞砸你对\(X\)如何影响\(Y\)的估计。
我们可以在\(图13.11\)中看到实际情况。这里有一堆观测值，在\(X\)的不同取值下，\(Y\)要么是 \(0\)，要么是 \(1\)。在那张图上，我们有拟合的\(OLS\)模型，还有拟合的\(logit\)模型（\(probit\)模型看起来会非常相似 ）。
## Probit/logit graph
set.seed(200)
tib <- tibble(X = runif(250)*8) %>% mutate(Y = ((X - 4)/4 + rnorm(250)*.5 > 0)*1)
ggplot(tib, aes(x = X, y = Y)) +
geom_point() +
geom_smooth(formula = y ~ x, method = 'lm', se = FALSE, color = 'black', linetype ='dashed') +
geom_smooth(formula = y ~ x, method = 'glm', method.args = list(family = binomial(link='logit')), se = FALSE, color = 'black') +
theme_pubr() +
scale_y_continuous(breaks = c(0,1)) +
annotate(geom = 'text', x = 1.75, y = .35, color = 'black', label = 'OLS Prediction',size = 10/.pt, family = 'sans') +
annotate(geom = 'text', x = 3.75, y = .35, color = 'black', label = 'Logit',size = 10/.pt, family = 'sans') +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"))
ggsave('Regression/ols_logit.pdf', width = 6, height = 5,device=cairo_pdf)我们看到了什么？正如之前提到的，我们看到\(OLS\)的预测值超出了\(0\)到\(1\)的范围。但除此之外，\(OLS\)模型和\(logit\)模型的斜率不同。在图表的右侧或左侧附近，这种情况尤为明显。这意味着什么呢？这意味着我们之前考虑无论如何都使用\(OLS\)的理由 —— 获取\(X\)对\(Y\)的影响，而不关心现实的预测 —— 无论如何都行不通了。
但为什么会这样呢？嗯，想想看 —— 如果\(X\)对\(Y\)有正向影响，但你已经预测\(Y\)有 \(0.999\) 的发生概率，那么\(X\)增加一个单位会有什么作用呢？几乎没什么作用，对吧？\(Y\)能增加的空间实在不大了。但如果当前对\(Y\)的预测是\(0.5\)，那么就有更多的增长空间，所以影响可能会更大。它应该更大！\(OLS\)无法考虑到这一点83。
普通最小二乘法无法解释\(X\)对\(Y\)的影响应如何产生差异。这与我们讨论多项式时影响产生差异的方式不同。曲线性不应像多项式那样基于\(X\)的值，而应基于\(Y\)当前的预测值。当\(Y\)的样本均值处于中间位置 —— 接近 \(0.5\)左右时，使用\(OLS\)并非大问题，尽管你可能还是想使用概率单位\(（probit）\)或对数单位\(（logit）\)模型。但如果\(Y\)的均值接近 \(0\) 或 \(1\)，你就会面临大问题，绝对需要一个非线性模型 。
解读广义线性模型\(（GLM）\)，会比解读普通最小二乘法\(（OLS）\)模型，稍微难一点儿。从一种非常狭隘、而且没什么用的意义上讲，解读方式完全一样。也就是说，每个系数仍然意味着”该变量变化一个单位，与指数（发生一个系数大小的变化）相关联”。不管我们说的是线性回归，还是非线性回归，这一点都成立。然而，不同之处在于，对于\(OLS\)而言，指数函数就是一个回归方程，能给出因变量的条件均值，而且我们知道如何解读因变量的变化。在广义线性模型中，指数函数是要代入链接函数的那个部分。所以，它关乎的是 “代入链接函数的那个量的变化单位”，这要难理解得多。
呈现概率单位\(（probit）\)模型、对数单位\(（logit）\)模型，以及更一般的非线性回归的结果时，常用边际效应84。边际效应基本上能把结果转换回线性模型的表述方式。在\(logit\)模型中，\(X\)的系数表示\(X\)的变化对指数的影响，而\(X\)的边际效应表示\(X\)的变化对因变量取值为\(1\)的概率的影响。
边际效应：在\(probit\)或\(logit\)模型中，解释变量\(X\)每增加一个单位，对二分类因变量\(Y\)取值为\(1\)/“真”的概率所产生的影响/关联。
从数学角度看，边际效应很直观，不过我们得再回到微积分知识。\(X\)对\(Y = 1\)的概率的边际效应，就是\(Pr(Y = 1)\)对\(X\)的导数，对于这些二元因变量模型，有一个简洁的解：
\[\frac{\partial Pr(Y = 1)}{\partial X} = \beta_1 Pr(Y = 1)(1 - Pr(Y = 1)) \quad (13.17)\]
所以，\(X\)的边际效应等于\(X\)的系数\(\beta_1\)乘以\(Pr(Y = 1)\)的预测概率（由模型得出，且以所有预测变量值为条件 ），再乘以 \(1\) 减去该概率。要注意，当\(Pr(Y = 1)\)接近\(0\)或 \(1\)时，\(Pr(Y = 1)(1 - Pr(Y = 1))\)很小，当\(Pr(Y = 1)=0.5\)时最大，所以我们延续\(图13.11\)的结论，即图 “中间” 的斜率最陡，两侧最平缓。
不过，有个问题。问题在于，对于给定变量，不存在单一的边际效应。我刚让大家回顾\(图13.11\)，那就继续看它！随着\(X\)增加，\(Y = 1\)的\(logit\)预测概率，可能只增加一点点（在图的左侧或右侧 ），也可能增加很多（在中间 ）。实际上，数据中的每个观测值，基于其自身的\(Pr(Y = 1)\)，都有自己的边际效应，这与边际效应的\(\beta_1 Pr(Y = 1)(1 - Pr(Y = 1))\)公式相符。如果模型预测你有\(90\%\)的概率\(Y = 1\)，你的边际效应就是\(\beta_1×0.9×0.1\)；如果预测你有\(40\%\)的概率，你的边际效应就是\(\beta_1×0.4×0.6\)。那我们该用什么值作为 “边际效应” 呢85？
对于这个问题，有两种常见的处理方法，其中一种我会大力推荐，另一种则不会。第一种（我不推荐的）是 “均值处的边际效应”。均值处的边际效应，首先取所有预测变量的均值。要是我们只研究\(X\)，那就取\(X\)的均值。然后，用这个均值去预测\(Pr(Y = 1)\)。最后，用这个预测出的\(Pr(Y = 1)\)来计算边际效应，即\(\beta_1 Pr(Y = 1)(1 - Pr(Y = 1))\)。这算出的是预测变量完全处于均值的观测值的边际效应。我为啥不推荐呢？因为它计算的不是针对特定个体的边际效应。说到底，这个边际效应是针对谁的呢86？而且，均值处的边际效应会独立地取每个变量的均值，但实际上这些变量很可能是相关的。要是比如有个人，\(X\)接近均值，而\(Z\)几乎总是接近最小值，这时取所有变量的均值就很奇怪。
平均边际效应：样本中所有观测个体边际效应的算术平均值。
我推荐的是平均边际效应。平均边际效应会计算每个个体观测值的边际效应，用该观测值的预测变量算出\(Pr(Y = 1)\)，再算出\(\beta_1 Pr(Y = 1)(1 - Pr(Y = 1))\)。接着，把每个个体观测值的边际效应算出来后，取平均值得到平均边际效应。就这么简单！这样你就得到了样本中边际效应的均值，能让你了解具有代表性的边际效应情况 87。
举个简单例子，假设我们估计了\(logit\)模型\(Pr(Y = 1) = logit(.1 +.05X -.03Z)\)。为简单起见，假设样本中只有两个观测值。第一个观测值，\(X = 1\)，\(Z = 2\)；第二个观测值，\(X = 0\)，\(Z = 20\)。
第一个观测值的预测值为：
\(Pr(Y = 1) = logit(.1 +.05(1) -.03(2)) = logit(.99) = \frac{e^{.99}}{1 + e^{.99}} \approx.729\)
类似地，第二个观测值的预测值为：
\(Pr(Y = 1) = logit(.1 + 0 -.03(20)) = logit(-.5) = \frac{e^{-.5}}{1 + e^{-.5}} \approx.622\)
那么，\(X\)的边际效应，因为\(X\)的系数是\(0.05\)，对第一个观测值而言，就是\(.05(.729)(1 -.729) =.010\)；对第二个观测值而言，是\(.05(.622)(1 -.622) =.012\)。将它们平均，得到平均边际效应为\((.010 +.012)/2 =.011\)
咱们来编写一些概率单位\(（probit）\)和对数单位\(（logit）\)模型的代码！在下面的代码中，我会展示如何估计一个\(logit\)模型，以及如何估计这些模型的平均边际效应。那\(probit\)模型呢？很方便的是，下面所有的代码示例都可用于估计\(probit\)模型，只需把看到的”\(logit\)“一词替换成”\(probit\)“就行。
在每种情况下，我们要看看随着时间推移，周末检查是变得更普遍了还是没那么普遍了，依旧用我们的餐厅检查数据。
对于\(R\)代码，你经常会在网上看到指南说用\(mfx\)包中的\(logitmfx\)函数，而非\(marginaleffects\)函数。这么做没问题，但要注意，它默认计算的是均值处的边际效应，而非平均边际效应。
# Load packages and data
library(tidyverse); library(modelsummary); library(marginaleffects)
df <- causaldata::restaurant_inspections
# Use the glm() function with
# the family = binomial(link = 'logit') option
m1 <- glm(Weekend ~ Year, data = df,family = binomial(link = 'logit'))
# See the results
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 44.236* |
| (23.500) | |
| Year | -0.024** |
| (0.012) | |
| Num.Obs. | 27178 |
| AIC | 2460.5 |
| BIC | 2476.9 |
| Log.Lik. | -1228.236 |
| F | 4.362 |
| RMSE | 0.09 |
# Get the average marginal effect of year
avg_slopes(m1, variables = 'Year')##
## Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## -0.000187 9.06e-05 -2.07 0.0389 4.7 -0.000365 -9.58e-06
##
## Term: Year
## Type: response
## Comparison: dY/dXimport statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# smf.logit wants the dependent variable to be numeric
df['Weekend'] = 1*df['Weekend']
# Use smf.logit to run logit
m1 = smf.logit(formula = "Weekend ~ Year", data = df).fit()## Optimization terminated successfully.
## Current function value: 0.045192
## Iterations 10m1.summary()| Dep. Variable: | Weekend | No. Observations: | 27178 |
|---|---|---|---|
| Model: | Logit | Df Residuals: | 27176 |
| Method: | MLE | Df Model: | 1 |
| Date: | 周六, 23 8月 2025 | Pseudo R-squ.: | 0.001782 |
| Time: | 10:58:01 | Log-Likelihood: | -1228.2 |
| converged: | True | LL-Null: | -1230.4 |
| Covariance Type: | nonrobust | LLR p-value: | 0.03625 |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 44.2360 | 23.502 | 1.882 | 0.060 | -1.828 | 90.300 |
| Year | -0.0244 | 0.012 | -2.088 | 0.037 | -0.047 | -0.002 |
# And get marginal effects
m1.get_margeff().summary()| Dep. Variable: | Weekend |
|---|---|
| Method: | dydx |
| At: | overall |
| dy/dx | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Year | -0.0002 | 9.05e-05 | -2.068 | 0.039 | -0.000 | -9.79e-06 |
causaldata restaurant_inspections.dta, use clear download
* Use the logit command to regress
logit weekend year
* and the margins command
margins, dydx(year)13.3 你的标准误很可能是错的
13.3.1 你的标准误为何很可能是错的
到目前为止，我们在回归分析中所做的所有工作，都聚焦于正确设定模型、妥善处理非线性关系、解读系数，当然，还有识别效应。但要是你没有对估计值的精确性（即标准误，\(SE\)）进行估计，那么估计出的模型对你而言意义不大88。
你应该还记得本章前面的内容，在估计模型\(Y = \beta_0 + \beta_1 X + \varepsilon\)时，我们得到了\(OLS\)系数估计值\(\hat{\beta}_1\)，以及标准误\(se(\hat{\beta}_1)\)的估计值（我们依据\(var(X)\)、\(n\)和对\(\sigma\)的估计来计算 ）。那么，对于\(\hat{\beta}_1\)的抽样分布，我们最合理的推测是，它服从均值为\(\hat{\beta}_1\)、标准差为\(se(\hat{\beta}_1)\)的正态分布。这就是标准误的作用 —— 它能让我们了解抽样分布，从而确定诸如合理估计值的范围、置信区间、统计显著性，以及排除某些不可能的理论 / 总体分布等信息。
这就涉及到我之前提到过、并说过之后会再回来讨论的误差项假设了。我们极其渴望知道抽样分布是什么样的。我们费心计算标准误，是因为想了解该抽样分布的标准差。但如果某些假设不成立，那么标准误就无法描述系数抽样分布的标准差。
我们必须回到的第一个假设是，我们得假定误差项\(\varepsilon\)本身必须服从正态分布89。关于误差项正态性的这一假设，让我们能够从数学上证明\(OLS\)系数服从正态分布。
所幸，这个假设并非至关重要。即便 \(\varepsilon\) 与正态分布相差甚远，\(OLS\)估计量的抽样分布往往也很接近正态。虽不完美，但一般无需对此大动肝火。话虽如此，的确存在一些怪异的误差项分布，会把事情搞砸。不过那是另一回事了。
第二个假设是误差项独立同分布。也就是说，我们需要假定误差项的理论分布，与其他观测的误差项以及同一观测的其他变量不相关（独立），且每个观测的误差项分布相同（同分布）。无论查看哪个数据点，误差项可能因观测而异，但始终来自同一分布。
这到底是什么意思呢？从假设不成立的情况去理解更容易。一种不成立的情况是自相关，即误差项之间以某种方式相关。处理跨多个时间段的数据（时间自相关）或地理上聚类的数据（空间自相关）时，这种情况常出现。
自相关：即误差项之间相互关联
以时间自相关为例，假设你在对美国失业率增长率按年做回归。经济往往会经历持续数年的起伏周期，如\(图13.12\)所示。想想这里的误差 —— 往往会连续出现几个正误差，接着是几个负误差。所以分布不可能独立，因为上一时期的误差会预示本时期的误差。
# GDP graph
data("economics")
ggplot(economics %>% dplyr::filter(date <= as.Date('2008-01-01')), aes(x=date, y=unemploy/pop)) +
geom_point(size=1) +
geom_smooth(formula='y~x', method='lm', se=FALSE, color='black') +
annotate(geom='text', x=as.Date('1970-01-01'), y=.03, label='OLS Fit', family='sans', size=10/.pt) +
theme_pubr() +
expand_limits(x = as.Date('2010-02-01')) +
labs(x='Month',y='US Unemployment Rate') +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"))
ggsave('Regression/unemployment.pdf', width = 6, height = 5,device=cairo_pdf)“独立同分布”假设不成立的另一种常见情况是存在异方差性。好长的一个词！这可真是教科书里的经典内容了。
异方差性：即误差项的方差与回归模型中的变量相关
异方差性指的是误差项分布的方差与模型中的变量相关。比如，假设你要回归分析某人的\(Instagram\)粉丝数量与其每日发帖所花时间的关系。你可能会发现，正如\(图13.13\)所示，在几乎不花时间发帖的人群中，粉丝数量的差异很小。但在花大量时间发帖的人群中，粉丝数量的差异极大。花少量时间意味着粉丝数量差异小，相应地，误差项的差异也小。但花大量时间发帖时，误差项的差异就很大。误差分布并非相同，因为分布的方差会随着模型中变量的值而变化。
# Instagram
set.seed(5403)
ig <- tibble(`Hours per Day`=runif(100, min=0, max=10)) %>%
mutate(`Followers`=runif(100, min=0, max=3*`Hours per Day`))
ggplot(ig, aes(x = `Hours per Day`, y = `Followers`)) +
geom_point(size=1) +
geom_smooth(formula='y~x', method='lm', se=FALSE, color='black') +
annotate(geom='text', x=2.5, y=5.5, label='OLS Fit', family='sans', size=13/.pt) +
theme_pubr() +
labs(x='Hours per Day on Instagram', y='Followers (Thousands)') +
theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family= "sans"))
ggsave('Regression/instagram.pdf', width = 6, height = 5,device=cairo_pdf)当 “独立同分布” 假设不成立时，我们又会得到一个还算便利的结果。这些不成立的假设，实际上并不会使系数的分布偏离正态。然而，它们确实会导致我们对标准误的估计（即\(\sqrt{\sigma^2 / (var(X)n)}\) ）出现错误。这意味着，我们必须想办法弄清楚该抽样分布的标准差究竟是多少。要进行这样的计算，就需要以某种方式考虑自相关或异方差性的影响。
13.3.2 修正你的标准误
在本节中，我要讲讲几种常见的标准误修正方法，当我们认为 “独立同分布” 假设不成立时，就可以用这些方法90。
首先，我们来考虑异方差性问题。异方差性指的是误差项\(\varepsilon\)的方差会随着\(X\)（或者说模型中任何一个预测变量）取值的不同而变化。在估计标准误时，异方差性会带来这样的问题：在方差较高的区域，不同样本中\(Y\)的条件均值波动会更大。不管这种情况出现在\(X\)取值较高、较低，还是中等的区域，都会使\(\hat{\beta}_1\)在不同样本间的变化更大（增大抽样分布的标准差 ），以适应这些变化的取值。但常规的\(OLS\)标准误只是对整个样本估计了相同的方差，这样就会低估实际的变化程度。
为了直观呈现这种情况，来看一个高度简化的例子，见\(图13.14\) 。该图包含两个数据集，在整个样本中，它们的总体误差项方差水平相同。左侧的数据集存在异方差性。想象从左侧的簇中选一个观测值，再从右侧的簇中选一个，然后在两者间画一条线。这就是样本量为\(2\)时我们的\(OLS\)估计。现在想象不断重复这个过程91。由于右侧的\(Y\)值相对固定，左侧的\(Y\)值变动很大，斜率在不同样本间变动也大 —— 是高方差的簇导致了这种情况，而样本整体的残差方差估计不会考虑到这一点。再看右侧的数据集，斜率在不同样本间仍会变动，但变动程度由每个簇内的方差决定，而这些簇内方差是相同的。所以，当我们用样本整体的残差方差去估计这种变动时，就能正确估计 。
# Heteroskedasticity demonstratoin
set.seed(200)
tibA <- tibble(X = c(runif(10),runif(10, 2, 3))) %>%
mutate(Y = c(rnorm(10, 0, 3),rnorm(10)))
set.seed(200)
tibB <- tibble(X = c(runif(10),runif(10, 2, 3))) %>%
mutate(Y = rnorm(20,0,sd(tibA$Y)))
p1 <- ggplot(tibA, aes(x = X, y = Y)) +
geom_point() +
theme_pubr()+
theme(text = element_text(size = 13, family="sans"),
axis.title.x = element_text(size = 13, family="sans"),
axis.title.y = element_text(size = 13, family= "sans"))
p2 <- ggplot(tibB, aes(x = X, y = Y)) +
geom_point() +
theme_pubr()+
scale_y_continuous(limits = c(min(tibA$Y),max(tibA$Y))) +
theme(text = element_text(size = 13, family="sans"),
axis.title.x = element_text(size = 13, family="sans"),
axis.title.y = element_text(size = 13, family= "sans"))
plot_grid(p1,p2)
ggsave('Regression/heterosked.pdf', width = 8, height = 4,device=cairo_pdf)面对这个问题，我们能做些什么呢？嗯，我们可以让标准误的估计量针对误差方差的情况 “量身定制”。不再使用\(\sqrt{\sigma^2 / (var(X)n)}\)，我们可以采用 “三明治估计量” ，简而言之，它能让用于计算\(var(X)\)的各个\(X\)值被放大、缩小，甚至与其他观测值的\(X\)值相乘92。
比较常用的一种异方差稳健三明治估计方法是休伯 - 怀特\(（Huber-White）\)法。在休伯 - 怀特法中，对单个\(X\)值进行缩放的方式是：取每个观测值的残差，将其平方，然后用这个平方残差作为权重来进行缩放。实际上，这样做会在计算方差时，给残差大的观测值赋予更高的权重。
实际上，有几种不同的估计量，它们以极为相似的方式、通过细微调整来针对异方差性调整标准误。一般来说，你在这里要找的是 “异方差稳健标准误”，或者简称 “稳健标准误”，其中大多数是三明治估计量，休伯 - 怀特标准误是经典应用93。不过，一种名为\(HC3\)的休伯 - 怀特替代方法，通常被认为是更理想的改进，也是理想的默认选项（Long 和 Ervin，2000 ），因为它使用 “刀切法” 来估计三明治估计量中的\(\Sigma\) 94。\(HC3\)在大样本中的表现大致相同，但在小样本中表现要好得多（MacKinnon 和 White，1985 ）。除了休伯 - 怀特和\(HC3\)，还有其他版本，可用于不同场景，但如今\(HC3\)被视为一个很好的起点。
接下来，我们可以探讨不满足独立性的误差。换句话说，误差之间是相关的。为什么相关误差会要求我们改变计算标准误的方式呢？相关误差会改变估计量的抽样分布，就如同相关数据会改变均值的抽样分布一样。假设你在收集人群样本并计算他们的身高。要是你每次从全球随机抽样选人，这会给身高均值的抽样分布带来特定的标准差。但要是你每次随机抽取家庭并计算平均身高，这些值就会是相关的 —— 你的身高与父母身高的相似度，很可能高于与陌生人身高的相似度。相较于随机抽取所有人的方法，我们更有可能抽到全是高个子的样本，从而得到很高的平均身高；或者抽到全是矮个子的样本，得到很矮的平均身高。所以均值的波动更大，抽样分布的标准差也更大。必须进行一些调整。
在确定要做何种调整时，我们有点犯难 —— 得弄清楚误差之间是如何相关的。毕竟，说 “我的误差是相关的” 没什么意义，不可能所有误差都相关。相反，我们会说它们依据某种特征呈现相关性。误差可能存在的每种不同相关方式，都需要不同的修正办法。
误差相关的一种常见情况是跨时间相关，也就是我们之前讨论过的基于时间的自相关。存在自相关时，我们可以使用异方差和自相关一致\(（HAC）\)的标准误95。计算\(HAC\)标准误的一种常用方法是纽维 - 韦斯特\(（Newey-West）\)估计量，它是三明治估计量的另一种形式。这种估计量从异方差稳健标准误估计入手，然后进行调整96。调整方法是，首先选定滞后阶数（即你预期误差存在相关性的时间段数量 ）。接着，计算在每个滞后阶数下观测到的自相关程度，将它们相加，且较短长度的滞后在求和中权重更大，以此得到调整因子97。
跨时间的自相关并非自相关的唯一形式，跨地理区域的相关性也会产生重大影响。虽然我不会在此深入探讨这个主题，但 康利\(（Conley）\)空间标准误 是三明治估计量的另一种形式，它在计算标准误时会考虑地理相邻区域间的自相关情况。
误差相关的另一种常见情形是存在分层结构。例如，假设你在研究给学生配备笔记本电脑对他们考试成绩的影响。实际上，某个班级的学生不只是配备或不配备笔记本电脑的问题，他们还在同一个班级，由同一位老师授课。基于是否配备笔记本电脑对他们考试成绩进行预测时，由于面临相似的环境，他们的成绩可能都会同时高于或低于预测值，也就是说，他们的误差会相互关联。我们会说这些误差在班级内存在聚类。
所幸，在这些情况下，我们可以应用聚类标准误。最常见的聚类标准误是\(（Liang-Zeger）\)标准误，它同样是三明治估计量的一种形式。使用聚类标准误时，你要明确指定一个分组，比如班级。这样一来，不再是用三明治估计量去放大或缩小单个\(X\)值，而是让同一组内的\(X\)值相互作用9899。
使用聚类标准误可以解释每个分组内误差之间的任何形式的相关性100。这种灵活性很让人欣慰，无需再对误差项做相关性为零的限制性假设。当然，统计精确性源于假设。要是你什么都不假设，数据也没法告诉你任何信息。所以，在过于宽泛的层面聚类（比如，要是你做的笔记本电脑与考试成绩的分析里只有两所学校，却在学校层面聚类，或者干脆在整个样本层面聚类 —— 即所有人的误差都相互关联 ），就等同于说 “呃…… 我也不知道咋回事，估计啥情况都有可能”，这样你会得到非常大的标准误，可能低估数据实际能说明的问题。所以，到一定程度时，你得确定这种相关性可以延伸到什么范围。
那该在什么层面聚类呢？首先，要是你有很强的理论依据，知道误差应该在哪些地方聚类，比如在班级内，那就这么做。除此之外，另一种常用方法是在处理层面聚类。比如，在笔记本电脑的例子中，如果笔记本电脑分配给了某些班级而非其他班级，你就在班级层面聚类。但如果是分配给了某些学校而非其他学校，你就在学校层面聚类。即便如此，也只有当处理效应在不同聚类间存在差异时，才需要聚类（Abadie 2017 ）。不过，这些都不是死板的规则，聚类决策需要你对所研究的特定领域有大量了解。看看有没有其他研究你这个主题的人，了解他们是如何聚类的。也可以看看第 \(固定效应\) 章 “专业人士怎么做” 部分对聚类的进一步讨论。
不过，关于\(Liang-Zeger\) 聚类标准误，有个重要点得记住：它们在聚类数量很多（比如超过\(50\)个）时效果很好。要是聚类数量比这少，它们仍比常规标准误有改进，但远达不到你期望的精确程度101。对于聚类数量较少的情况，你或许可以改用野生聚类自助法标准误，它有点像是我们正在讨论的聚类标准误和我下一节要讲的自助法标准误的混合体。关于野生聚类自助法标准误，可在\(第15章\)了解更多内容。
所以，难道全都是三明治估计量的天下了？不！还有另一种估计标准误的方法，应用非常广泛，那就是自助法标准误。自助法在统计学里近乎神奇。我们用标准误是想干嘛呢？它们的意义何在？我们是想衡量抽样分布，对吧？想知道如果在大量不同样本中估计同一个统计量，会出现什么情况。
所以自助法的思路是：“那我们就在大量不同样本中估计同一个统计量，看看这些样本间的分布是怎样的。嗒哒！这就是抽样分布。” 哦，嘿，有道理呀。我们一开始怎么就没想到呢？
但等等，怎么能在大量不同样本中估计统计量呢？我们只有一个样本呀！这就是自助法的巧妙之处。它利用我们现有的这一个样本，但对其进行有放回抽样。自助法的流程如下102：
从包含\(N\)个观测值的数据集开始。
从该数据集中随机抽取\(N\)个观测值，允许同一观测值被多次抽取103。这就是一个 “自助样本”。
用这个自助样本估计我们感兴趣的统计量。
重复步骤\(1-3\)很多次。
查看我们多次重复步骤\(1-3\)得到的估计值的分布。该分布的标准差就是估计值的标准误。分布的第\(2.5\)和第\(97.5\)百分位数给出\(95\%\)置信区间，依此类推。
这真的有效吗？咱们举个非常简单的例子。我们的数据里有四个观测值：\(1, 2, 3, 4\) 我们要估计均值。样本均值显然是\((1 + 2 + 3 + 4)÷4 = 2.5\)。数据的样本方差是\([(-1.5)^2 + (-0.5)^2 + (0.5)^2 + (1.5)^2]÷(4 - 1) = 1.67\)。所以，用常规方法计算的标准误是\(\sqrt{1.67÷4} = 0.645\)。
好啦，现在试试自助法。我们得到第一个自助样本：\(1,1,4,4\)—— 记住，同一观测值能被重复抽样，这就意味着原始数据里的有些观测值会多次出现，有些则根本不出现。这个自助样本的均值是\((1 + 1 + 4 + 4)÷4 = 2.5\)。我们再次抽样，得到 \(1,3,4,4\)，均值是\(3\)。再抽一次！下一个样本是\(1,2,4,4\)，均值是\(2.75\)。接着抽！\(1,2,2,3\)，均值是 \(2\)。到目前为止，我们抽取的自助样本的均值的平均值是\((2.5 + 3 + 2.75 + 2)÷4 = 2.56\)，标准差是\(0.427\)。
\(0.427\)和我们知道的正确值 \(0.645\) 不完全一样，但我们还远没完成。进行自助法时，一般就算不抽几千次，也得抽几百次。要记住，我们这是在用实际的重复抽样来模拟理论分布。我们得更接近理论分布所基于的理论上无限多的观测值。
咱们把规模搞大些。别用\(1-4\)的观测值了，换成\(1-1000\)，而且计算\(2000\)个自助样本，而非\(4\)个。现在，样本均值是\(500.5\)，用常规方法算出的均值标准误是 \(9.13\)。自助分布的均值也是\(500.5\)，\(2000\)个自助样本的标准差是\(9.06\)。比\(9.13\)略小，但基本和我们预期的一致104。这表明，在这个特定例子中，那些基于理论分布的假设对精确性的估计可能有点悲观了105。
这就是计算自助法标准误的基本流程。确实，计算量有点大。但反过来讲，它有时能指出传统估计值是过于保守（或不够保守 ），还能让你探究抽样分布中任何你感兴趣的部分 —— 中位数、不同百分位数、偏度等等，不只是标准差 —— 而且在估计非常规情况的标准误时，它非常有用。要是我们有个估计量，不确定它是否服从正态分布，或任何已知分布，咋办？没问题！用自助法就行。
当然，自助法也有自身的一系列假设和不足 —— 它或许很神奇，但就算真有魔法，也没有免费午餐这回事。自助法在小样本情况下表现不太好106，对于 “极值” 分布（高度偏态的数据 ）表现也不佳，而且要是数据中的观测值相互关联（比如在面板数据中，多个观测值代表同一主体在不同时间的情况 ），会存在一些至少需要解决的难题。再比如，标准自助法在存在自相关时表现不好。但对于样本量合适、数据是行为良好的独立同分布的估计，自助法会是很棒的选择，尤其是在常规标准误难以计算的情境中。另外，\(第15章\)讨论了一些自助法的替代版本，解决了其中部分问题。
咱们把这些标准误调整方法付诸实践。在接下来的每个代码示例中，我会先展示如何实现异方差稳健标准误、异方差和自相关一致\(（HAC）\)标准误，然后是聚类标准误。
首先是稳健标准误。我们这里继续使用餐厅检查数据。
关于稳健标准误，有个重要点要了解，尤其是在编写代码实现它们时，实际上有很多不同的计算方法。每种方法都基于同一思路进行不同的细微调整。有些版本在小样本下表现更好；有些在特定类型的异方差情况下表现更佳。以下所有方法都能调用多种不同的稳健估计量（通常命名为\(HC0\)、\(HC1\)、\(HC2\)、\(HC3\)之类 ）。但不同编程语言会选择不同的默认值107。比如，\(Stata\) 中那个广为人知的简便稳健选项，遗憾的是默认用\(HC1\)，所以我们得用\(vce(hc3)\)来获取理应作为默认的\(HC3\)。大多数\(R\)语言方法默认用\(HC3\)。一定要查看文档，确认你用的是哪种，还要想想自己应该用哪种。
# Load packages
library(tidyverse)
# For approach 1
library(AER); library(sandwich)## Loading required package: car## Loading required package: carData##
## Attaching package: 'carData'## The following object is masked from 'package:causaldata':
##
## Mroz##
## Attaching package: 'car'## 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##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Loading required package: survival# For approach 2
library(modelsummary)
# For approach 3
library(fixest)
df <- causaldata::restaurant_inspections
# In R there are many ways to get at robust SEs
# For the first two methods we'll look at, we need to estimate
# the regression on its own
# (most types of regressions will work, not just lm()!)
m1 <- lm(inspection_score ~ Year + Weekend, data = df)
# First, the classic way, sending our regression
# object to AER::coeftest(), specifying the kind
# of library(sandwich) estimator we want
# vcovHC for heteroskedasticity-consistent
coeftest(m1, vcov = vcovHC(m1))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 185.3800329 12.1498220 15.2578 < 2.2e-16 ***
## Year -0.0456403 0.0060438 -7.5516 4.434e-14 ***
## WeekendTRUE 2.0571662 0.3529064 5.8292 5.632e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Second, we can take that same regression object
# and, do the robust SEs inside of our msummary()
msummary(m1, vcov = 'robust',
stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 185.380*** |
| (12.150) | |
| Year | -0.046*** |
| (0.006) | |
| WeekendTRUE | 2.057*** |
| (0.353) | |
| Num.Obs. | 27178 |
| R2 | 0.003 |
| R2 Adj. | 0.003 |
| AIC | 176731.8 |
| BIC | 176764.7 |
| Log.Lik. | -88361.913 |
| F | 46.723 |
| RMSE | 6.25 |
| Std.Errors | HC3 |
# Third, we can skip all that and use a regression
# function with robust SEs built in, like
# fixest::feols()
feols(inspection_score ~ Year + Weekend, data = df, se = 'hetero')## OLS estimation, Dep. Var.: inspection_score
## Observations: 27,178
## Standard-errors: Heteroskedasticity-robust
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 185.380033 12.148904 15.25899 < 2.2e-16 ***
## Year -0.045640 0.006043 -7.55219 4.4144e-14 ***
## WeekendTRUE 2.057166 0.351245 5.85678 4.7735e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 6.24818 Adj. R2: 0.00267# (this object can be sent to modelsummary
# or summary as normal)import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# We can simply add cov_type = 'HC3' to our fit!
m1 = smf.ols(formula = 'inspection_score ~ Year + Weekend',
data = df).fit(cov_type = 'HC3')
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.003 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.003 |
| Method: | Least Squares | F-statistic: | 46.72 |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 5.54e-21 |
| Time: | 10:58:03 | Log-Likelihood: | -88362. |
| No. Observations: | 27178 | AIC: | 1.767e+05 |
| Df Residuals: | 27175 | BIC: | 1.768e+05 |
| Df Model: | 2 | ||
| Covariance Type: | HC3 |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 185.3800 | 12.150 | 15.258 | 0.000 | 161.567 | 209.193 |
| Weekend[T.True] | 2.0572 | 0.353 | 5.829 | 0.000 | 1.365 | 2.749 |
| Year | -0.0456 | 0.006 | -7.552 | 0.000 | -0.057 | -0.034 |
| Omnibus: | 3522.402 | Durbin-Watson: | 1.921 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 5062.647 |
| Skew: | -1.002 | Prob(JB): | 0.00 |
| Kurtosis: | 3.673 | Cond. No. | 6.79e+05 |
Notes:
[1] Standard Errors are heteroscedasticity robust (HC3)
[2] The condition number is large, 6.79e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
causaldata restaurant_inspections.dta, use clear download
* You commonly see people do ", robust" in Stata, but for the
* reasons mentioned above you want vce(hc3) instead.
reg inspection_score weekend year, vce(hc3)接下来，我们要用纽维 - 韦斯特（Newey - West）标准误处理自相关问题。要记住，在每种情况下，你都得选定最大滞后阶数（或者让系统为你选定一个 ）。本书不太聚焦时间序列内容，所以我不会深入探讨。但基本来说，要仔细研究数据的自相关结构，以此确定滞后阶数设为多少合适。
因为这些方法是用于时间序列数据的，所以我们首先得把数据调整成时间序列形式。
# Load packages
library(tidyverse)
# For approach 1
library(AER); library(sandwich)
# For approach 2
library(fixest)
df <- causaldata::restaurant_inspections
# Turn into a time series
df <- df %>%
group_by(Year) %>%
summarize(Avg_Score = mean(inspection_score), Pct_Weekend = mean(Weekend)) %>%
# Only the years without gaps
filter(Year <= 2009) %>%
# feols requires an index to do this
mutate(index = 1)
# Perform our time series regression
m1 <- lm(Avg_Score ~ Pct_Weekend, data = df)
# Use the same coeftest method from robust SEs, but use the NeweyWest matrix. We can set maximum lags ourselves,
# or it will pick one using an automatic procedure (read the docs!)
coeftest(m1, vcov = NeweyWest)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.7736 1.9631 47.767 4.081e-11 ***
## Pct_Weekend 42.7314 48.6663 0.878 0.4055
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Or just use feols from fixest
# Note that fixest applies a small-sample correction and uses a "classic" form of HAC that does not apply "pre-whiting",
# so results are different
feols(Avg_Score ~ Pct_Weekend, data = df,panel.id = ~ index + Year, vcov = 'NW')## OLS estimation, Dep. Var.: Avg_Score
## Observations: 10
## Standard-errors: Newey-West (L=1)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.7736 1.06529 88.02613 1.5970e-14 ***
## Pct_Weekend 42.7314 27.21926 1.56989 1.5089e-01
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.79959 Adj. R2: -0.050873import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# Get our data into a single time series!
df = df.groupby('Year').agg({'inspection_score': 'mean','Weekend':'mean'})
# Only use the years without a gap
df = df.query('Year <= 2009')
# We can simply add cov_type = 'HAC' to our fit, with maxlags specified
m1 = smf.ols(formula = 'inspection_score ~ Weekend',data = df).fit(cov_type = 'HAC', cov_kwds={'maxlags':1})
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.066 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | -0.051 |
| Method: | Least Squares | F-statistic: | 3.081 |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.117 |
| Time: | 10:58:03 | Log-Likelihood: | -20.065 |
| No. Observations: | 10 | AIC: | 44.13 |
| Df Residuals: | 8 | BIC: | 44.74 |
| Df Model: | 1 | ||
| Covariance Type: | HAC |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 93.7736 | 0.953 | 98.416 | 0.000 | 91.906 | 95.641 |
| Weekend | 42.7314 | 24.346 | 1.755 | 0.079 | -4.985 | 90.448 |
| Omnibus: | 2.274 | Durbin-Watson: | 0.500 |
|---|---|---|---|
| Prob(Omnibus): | 0.321 | Jarque-Bera (JB): | 1.042 |
| Skew: | 0.392 | Prob(JB): | 0.594 |
| Kurtosis: | 1.626 | Cond. No. | 89.4 |
Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 1 lags and without small sample correction
# Note that Python uses a "classic" form of HAC that does not apply "pre-whiting", so results are often different from other languagescausaldata restaurant_inspections.dta, use clear download
* Turn our data into a single time series
collapse (mean) weekend inspection_score, by(year)
* Only use the run of years without gaps
keep if year <= 2009
* Tell Stata we have a time series dataset
tsset year
* Use "newey" instead of "regress"
* We must set maximum lags with lag()
newey inspection_score weekend, lag(3)对于最后一种三明治估计量，我们可以使用聚类标准误。在数据中，其实并没有特别合理的分组来进行聚类。但由于这只是技术演示，我们可以随意操作。我们将按检查是否在周末进行聚类。这个分组只有两个取值，聚类数量很少，可能会得到奇怪的结果。
# Load packages
library(tidyverse)
# For approach 1
library(AER); library(sandwich)
# For approach 2
library(modelsummary)
# For approach 3
library(fixest)
df <- causaldata::restaurant_inspections
# In R there are many ways to get at clustered SEs. For the first
# two methods we'll look at, we need to estimate the regression on its own. (most commands will work, not just lm()!)
m1 <- lm(inspection_score ~ Year + Weekend, data = df)
# First, the classic way, sending our regression object to AER::coeftest(), specifying the kind of library(sandwich) estimator.
# We want vcovCL for clustering.
# Note the ~ before Weekend; this is technically a formula
coeftest(m1, vcov = vcovCL(m1, ~Weekend))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 185.3800329 3.2643454 56.789 < 2.2e-16 ***
## Year -0.0456403 0.0016238 -28.108 < 2.2e-16 ***
## WeekendTRUE 2.0571662 0.0014011 1468.257 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Second, we can pass it to msummary which will cluster if we send a formula with the cluster variable to vcov
msummary(m1, vcov = ~Weekend,stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 185.380*** |
| (3.264) | |
| Year | -0.046*** |
| (0.002) | |
| WeekendTRUE | 2.057*** |
| (0.001) | |
| Num.Obs. | 27178 |
| R2 | 0.003 |
| R2 Adj. | 0.003 |
| AIC | 176731.8 |
| BIC | 176764.7 |
| Log.Lik. | -88361.913 |
| RMSE | 6.25 |
| Std.Errors | by: Weekend |
# Third, we can use a regression function with clustered SEs built in, like fixest::feols(). Don't forget the ~ before Weekend.
feols(inspection_score ~ Year + Weekend, cluster = ~Weekend, data = df)## OLS estimation, Dep. Var.: inspection_score
## Observations: 27,178
## Standard-errors: Clustered (Weekend)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 185.380033 3.264345 56.7893 0.01120904 *
## Year -0.045640 0.001624 -28.1076 0.02263987 *
## WeekendTRUE 2.057166 0.001401 1468.2537 0.00043359 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 6.24818 Adj. R2: 0.00267import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# We can simply add cov_type = 'cluster' to our fit, and specify the groups
m1 = smf.ols(formula = 'inspection_score ~ Year',data = df).fit(cov_type = 'cluster',cov_kwds={'groups': df['Weekend']})
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.002 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.002 |
| Method: | Least Squares | F-statistic: | 2810. |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.0120 |
| Time: | 10:58:03 | Log-Likelihood: | -88373. |
| No. Observations: | 27178 | AIC: | 1.768e+05 |
| Df Residuals: | 27176 | BIC: | 1.768e+05 |
| Df Model: | 1 | ||
| Covariance Type: | cluster |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 186.1690 | 1.714 | 108.619 | 0.000 | 182.810 | 189.528 |
| Year | -0.0460 | 0.001 | -53.008 | 0.000 | -0.048 | -0.044 |
| Omnibus: | 3528.301 | Durbin-Watson: | 1.921 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 5073.965 |
| Skew: | -1.003 | Prob(JB): | 0.00 |
| Kurtosis: | 3.674 | Cond. No. | 6.79e+05 |
Notes:
[1] Standard Errors are robust to cluster correlation (cluster)
[2] The condition number is large, 6.79e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
causaldata restaurant_inspections.dta, use clear download
* for basically any regression command,
* just add ", vce(cluster)"
reg inspection_score weekend year, vce(cluster weekend)最后，来说说自助法标准误。这里需要做一个必要的选择，即生成多少个自助样本。那么…… 该生成多少个呢？其实没有严格固定的规则。对于大多数应用场景，几千个样本应该就可以了。那我们就这么做。下面的代码仅展示了一种简单的、假设每个观测值相互独立的自助法。但一般来说，如果你查阅相关文档或\(第15章\)，就会发现像聚类自助法（某些观测值会被同时纳入或排除 ）这类方法也并不复杂。
# Load packages and data
library(tidyverse); library(modelsummary); library(sandwich)
df <- causaldata::restaurant_inspections
# Let's run our standard model from before
m <- lm(inspection_score ~ Year + Weekend, data = df)
# And use the vcov argument of msummary with vcovBS from the sandwich package to get boostrap SEs with 2000 samples
msummary(m, vcov = function(x) vcovBS(x, R = 2000),stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 185.380*** |
| (12.280) | |
| Year | -0.046*** |
| (0.006) | |
| WeekendTRUE | 2.057*** |
| (0.346) | |
| Num.Obs. | 27178 |
| R2 | 0.003 |
| R2 Adj. | 0.003 |
| AIC | 176731.8 |
| BIC | 176764.7 |
| Log.Lik. | -88361.913 |
| RMSE | 6.25 |
| Std.Errors | Custom |
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# There are bootstrap estimators in Python but generally not designed for OLS, So we'll make our own. Also see the Simulation chapter.
def our_reg(DF):
# Resample with replacement
resampler = np.random.randint(0,len(DF),len(DF))
DF = DF.iloc[resampler]
# Run our estimate
m = smf.ols(formula = 'inspection_score ~ Year',data = DF).fit()
# Get all coefficients
return(dict(m.params))
# Run the function 2000 times and store results
results = pd.DataFrame([our_reg(df) for i in range(0,2000)])
# Mean and standard deviation are estimate and SE
results.describe()## Intercept Year
## count 2000.000000 2000.000000
## mean 186.038522 -0.045960
## std 12.316787 0.006126
## min 144.732242 -0.069386
## 25% 177.768489 -0.050014
## 50% 185.928616 -0.045910
## 75% 194.228833 -0.041830
## max 233.185451 -0.025420causaldata restaurant_inspections.dta, use clear download
* for most commands, just add ", vce(boot, reps(number-of-samples))"
reg inspection_score weekend year, vce(boot, reps(2000))13.4 回归的其他注意事项
讲了这么多，我们居然还没讲完回归？这一章得有多长啊108？唉，和任何统计方法一样，回归有无数细微的角落和缝隙，得去探究这种方法对你的数据是否真的适用。而且，像回归这样应用广泛的方法，花时间去探究我们力所能及的方面是值得的。
在最后这一部分，我们会看看在运用回归时，还需要牢记的另外四个技巧或注意事项。
13.4.1 样本权重
我们之前讲利用样本关系来了解总体关系时，都假定样本是从总体中随机抽取的。所以，要是我们有一个比如关于巴西人的样本，那么每个巴西人进入我们样本的可能性是均等的，不管他们住在里约热内卢还是乡村，不管是富人还是穷人，不管受教育与否等等。
由于诸多原因，这种情况很少见。有些人比其他人更容易被调查或采集数据。有些抽样方法甚至都不尝试对所有人进行随机抽样，而是基于便利性、地理位置或其他诸多特征来抽样。
出现这种情况时，显然，结果更能代表那些更有可能被抽样的人群。我们可以用样本权重来解决这个问题，即给样本中的每个观测值赋予某种重要性度量，让有些观测值被视为比其他观测值更重要（从而在估计中更有影响力 ）。
在二元回归应用里，斜率估计值是\(Cov(X, Y)/Var(X)\)。不使用样本权重时，我们计算协方差的方法是，对每个观测值，用\((X - \bar{X})\)乘以\((Y - \bar{Y})\)，然后对所有观测值的这个乘积取平均。使用样本权重\(W\)时，操作基本一样，但对每个观测值，先把\((X - \bar{X})(Y - \bar{Y})\)乘以该观测值的权重\(W\)，之后再对所有结果取平均109 。权重\(W\)越大，在计算协方差、方差等时，该观测值的影响就越大，在最终估计中也就更能体现其代表性。
样本权重：计算样本统计量时，应用于每个观测值的一种缩放因子 。
不出所料，整个过程被称为加权最小二乘法。
加权最小二乘法：运用样本权重进行的普通最小二乘法。
加权的理由不止一个。我刚提到的问题 —— 某些人比其他人更可能被抽样 —— 在社会科学中是个很常见的问题。实际上，这种情况太普遍了，很多调查，当然还有大多数大型且执行良好的调查，比如多数国家的人口普查，都会为你准备好样本权重。调查规划者很清楚他们的抽样是如何进行的，所以对于不同的个体、家庭、企业或其他研究对象被纳入样本的可能性，他们心里有数。
这些权重有时是通过对比数据中具有某一特定属性的人群比例与总体中的相应比例来确定的。如果总体中有\(500\)名男性和\(500\)名女性（各占\(50\%\)），而我们的\(100\)人样本中有\(70\)名男性和\(30\)名女性，那么我们可以说每个男性被纳入样本的概率是\(70/500=14\%\)，每个女性是\(30/500=6\%\)。这些概率会为我们确定样本权重提供依据。
然后，我们会用这些概率的倒数来对每个个体进行加权。这样，男性会得到\(1/0.14=7.143\)的权重，女性会得到\(1/0.06=16.667\)的权重。这样一来，那些原本不太可能被抽样的人群在分析中被赋予了更高的权重，让我们还原出总体中的人口结构比例，使我们的估计结果能更好地代表总体。
样本概率倒数权重：一种样本权重，其中每个观测值的权重为其被纳入样本的概率的倒数 。
比如，要是我们仅用样本去估计男性所占比例，会估计出总体中男性占\(70\%\)。显然错了！但如果进行加权，计算结果会是\((7.143×70 + 16.667×0)/(7.143×70 + 16.667×30) = 50\%\) 110。完美！估计回归系数时也是同样的道理。不加权的话，系数会因抽样过程产生偏差；加权后，能更好地反映总体中的关系。
样本权重是有效的。要是你的数据带有样本权重（而且这些权重估计得当，专业调查里的权重通常会是这样 ），仔细研读其含义并加以运用。一定要选用契合你分析场景的权重 —— 调查数据针对不同的观测层级，往往有不同的权重。若要推广到个体层面，就用个人权重；若要推广到家庭层面，就用家庭权重，以此类推。
不过，加权还有其他原因。要是你使用的是聚合数据，一个重要原因就来了。比如，你可能要分析发放笔记本电脑对考试分数的影响，但没用单个学生及其分数的数据，而是用了班级平均分111。
直接分析班级平均数据会存在一些问题。有些班级规模比其他班级大。规模大的班级代表的学生数量更多，所以要是你想得出 “一台笔记本电脑使学生分数提高\(\hat{\beta}\)” 这样的结果，聚合分析就很难处理，因为你的结果实际是基于班级平均层面的。其次，不管你用的是什么平均数据，来自规模大的班级的数据噪声会更小，因为样本量越大，均值的抽样分布标准差越小。所以，你数据里的某些均值会比其他均值估计得更准确，但普通最小二乘法会把它们同等看待。
这是两个不同问题，有两种不同（但类似）的解决办法。
第一个问题及解决办法很简单。要是一个班级代表\(30\)人，另一个代表\(40\)人，咋办？容易！就把第一个班级的数据算\(30\)次，第二个班级算\(40\)次。这些是 “频率权重”，权重就是变量所代表的观测数量。这样能让结果回到个体、非聚合层面来解释112。虽然不如真的有非聚合数据好（因为会丢失个体间的差异 ），但它让解释更容易，还能体现规模差异的影响。
频率权重：在聚合数据中，样本权重等于所代表的观测数量，用于让结果能在非聚合层面进行解释 。
第二个问题和均值方差的差异有关，稍复杂点，但也没复杂多少。针对这个问题，我们用 “逆方差加权”。可以计算每个聚合观测均值的方差，然后用方差的倒数为每个观测加权。
样本均值的方差是\(\sigma^2/N\)，这里\(N\)是样本中的观测数，也就是聚合到该观测里的观测数量。其倒数是\(N/\sigma^2\)。假设\(\sigma^2\)在样本中恒定，那这部分在加权时就不用考虑了（因为对所有人都适用，权重是相对的 ），最后权重就是\(N\)。每个聚合观测的权重就是聚合到它里面的观测数量。
等等…… 频率权重也是按用于聚合的观测数量来加权的呀。那频率权重和逆方差权重有啥区别呢？实际上，在估计系数时两者没啥区别 —— 频率权重和逆方差权重会得出完全相同的系数。不过，在计算标准误时就有区别了。频率权重的处理方式，就好像每个聚合后的观测值真的是一组相互独立、完全相同的单独观测值的集合。所以，在计算标准误里的\(\sigma/\text{var}(X)\)时，它用的是未聚合的样本量来计算\(\sigma\)。而逆方差权重呢，会把数据认作是一组均值，所以用聚合后的样本量来计算\(\sigma\)。这样一来，频率权重的标准误就总会比逆方差权重的标准误小113。
逆方差权重除了用于聚合数据，还有很多其他应用场景。在很多情况下，我们知道有些观测值比其他观测值的方差大得多。比如，处理变量本身对有些观测值来说，可能比其他观测值的变异性大很多，这些观测值在估计中会被赋予更大的权重，而我们可以通过取方差的倒数来构造权重，纠正这种情况。问题就解决了。这种方法在固定效应模型中的一种应用，你可以在\(第16章\)看到。
逆方差权重出现的另一个地方是元分析（meta-analysis），不过我们没篇幅展开讲了。在元分析中，你会研究一堆关于同一主题的不同研究，尝试估计出整体效应，把这些研究的结果聚合起来。有些研究的估计值标准误大，有些则小。逆方差权重能让估计更精确的效应得到更大的权重。
咱们用样本权重来估计一些回归模型吧。大多数命令都允许你直接给它们指定一个权重变量。确定用哪种加权方式是你的活儿，你只要把最终的权重传进去就行。所以，要是你想用抽样概率倒数权重，最好自己算出那个倒数！
这也意味着，要是某种加权方法需要你调整标准误，它不会自动发生。比如，你想做频率权重，那你要么自己调整标准误，要么完全不用权重，直接把每个观测值复制若干次（次数等于权重），然后不加权重跑回归。在\(R\)里，加载\(tidyverse\)包后，你可以把数据框\(df\)传给\(slice(rep(1:nrow(df)), weight)\) 来实现。
\(Python\)里要麻烦点 —— 你可以用\(expanded_df = df.loc[numcopies]\)来做，但得自己算出\(numcopies\)（即你想复制每一行的次数）。\(itertools\)包可能在这方面有帮助。
\(Stata\) 在这方面是个例外，因为它有好几种权重类型\(（aweight, fweight, iweight, pweight）\)，处理方式不同，而且不是所有命令都接受所有类型的权重。一般来说，\(pweight\)用于抽样概率倒数权重，\(aweight\)用于方差加权，\(fweight\)用于频率加权，这些会帮你妥善处理标准误调整114。用\(Stata\)的权重之前，看看 \(Stata\) 的\(help \space weight\)文档。
对于这些例子，咱们要用餐厅检查数据的聚合版本做些逆方差加权回归。首先，咱们把数据聚合到不同的连锁餐厅层面，统计一下要平均的检查次数。具体来说，看看平均检查分数是否和连锁餐厅首次出现在样本中的年份有关。然后，因为样本均值的方差与样本中观测数量的倒数成正比，所以咱们要用观测数量的倒数的倒数（也就是观测数量）来加权。
# Load packages and data
library(tidyverse); library(modelsummary)
df <- causaldata::restaurant_inspections
# Turn into one observation per chain
df <- df %>%
group_by(business_name) %>%
summarize(Num_Inspections = n(),
Avg_Score = mean(inspection_score),
First_Year = min(Year))
# Add the weights argument to lm
m1 <- lm(Avg_Score ~ First_Year, data = df,weights = Num_Inspections)
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 178.197*** |
| (37.448) | |
| First_Year | -0.042** |
| (0.019) | |
| Num.Obs. | 1618 |
| R2 | 0.003 |
| R2 Adj. | 0.003 |
| AIC | 10386.8 |
| BIC | 10403.0 |
| Log.Lik. | -5190.416 |
| F | 5.098 |
| RMSE | 4.60 |
import statsmodels.formula.api as smf
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# Aggregate the data
df['Num_Inspections'] = 1
df = (df.groupby('business_name').agg({'inspection_score': 'mean','Year': 'min','Num_Inspections': 'sum'}))
# Here we call a special WLS function
m1 = smf.wls(formula = 'inspection_score ~ Year', weights = df['Num_Inspections'], data = df).fit()
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.003 |
|---|---|---|---|
| Model: | WLS | Adj. R-squared: | 0.003 |
| Method: | Least Squares | F-statistic: | 5.098 |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.0241 |
| Time: | 10:58:28 | Log-Likelihood: | -5190.4 |
| No. Observations: | 1618 | AIC: | 1.038e+04 |
| Df Residuals: | 1616 | BIC: | 1.040e+04 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 178.1965 | 37.448 | 4.758 | 0.000 | 104.744 | 251.649 |
| Year | -0.0422 | 0.019 | -2.258 | 0.024 | -0.079 | -0.006 |
| Omnibus: | 850.037 | Durbin-Watson: | 1.682 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 7498.460 |
| Skew: | -2.301 | Prob(JB): | 0.00 |
| Kurtosis: | 12.489 | Cond. No. | 6.59e+05 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.59e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
causaldata restaurant_inspections.dta, use clear download
* Aggregate the data
g Num_Inspections = 1
collapse (mean) inspection_score (min) year ///
(sum) Num_Inspections, by(business_name)
* Since we want inverse-variance weights,
* we use aweights, or aw for short
reg inspection_score year [aw = Num_Inspections]当然，这些代码示例只是展示如何用权重运行单个回归。要是你处理的调查数据有复杂的权重结构，而且需要反复使用，你可能希望编程语言能明确识别这一点并提供帮助。在\(R\)中，你可以用\(survey\)包实现。在\(Stata\)中，你可以通过\(help \space svy\)和\(help \space svyset\)查看一系列\(svy\)命令及前缀 。
13.4.2 共线性
我之前讲分类预测变量时，提到过完全多重共线性。完全多重共线性是指模型中某些变量的线性组合能完美预测模型里的另一个变量。出现这种情况时，模型无法估计。要是模型是\(Sales = \beta_0 + \beta_1 Winter + \beta_2 NotWinter\)，那么\(\beta_0 = 15\)、\(\beta_1 = -5\)、\(\beta_2 = 0\)产生的预测（和残差 ），与\(\beta_0 = 10\)、\(\beta_1 = 0\)、\(\beta_2 = 5\)，以及无数其他参数估计组合的结果完全一样。普通最小二乘法没法在这些组合间做选择，也就没法估计模型115。
但要是没有完全多重共线性，可预测变量之间高度相关，会怎样呢？毕竟，换个角度看，\(Winter\)和\(NotWinter\)就是两个相关系数为\(-1\)的变量。要是有两个变量相关系数是\(0.99\)或\(-0.99\)（而非1或\(-1\) ），甚至\(0.9\)或\(-0.9\)呢？这种情况我们称为高度共线性（但非完全共线性 ）。
共线性：模型中的一个预测变量可由模型中其他变量的线性组合很好地预测 。
除了分类变量的例子，当有多个变量实际测量的是同一事物时，也常出现高度共线性。比如，你想用 “情商” 分数预测某人是否会升职。你给大家做一个情商测试，测试分两部分，题目基本一样，只有细微差别。然后，你把测试的第一部分和第二部分作为模型里的单独预测变量，希望从整体上了解情商如何预测升职。这两部分的分数很可能高度相关116。
高度共线性往往会让估计值变得非常敏感、波动大，导致抽样分布范围很宽，还会使标准误增大。
共线性为何会增大标准误？咱们接着用\(Winter/Sales\)的例子来说。把\(Winter\)变量换成\(WinterExceptJan1\)（\(1\)月除外的冬季 ）。现在，预测变量是\(WinterExceptJan1\)（冬季里除\(1\)月\(1\)日外的日子取\(1\)）和\(NotWinter\)（非冬季日子取\(1\)）。\(1\)月\(1\)日对这两个变量而言都是\(0\)。这两个变量的相关性会非常接近\(-1\)。除\(1\)月\(1\)日外，对于所有观测，\(WinterExceptJan1 + NotWinter = 1\)都成立。所以，我们接近完美线性组合，但又不完全是。
在这种情况下，对于模型\(Sales = \beta_0 + \beta_1 WinterExceptJan1 + \beta_2 NotWinter\)，\(\beta_0 = 15\)、\(\beta_1 = -5\)、\(\beta_2 = 0\)不会和\(\beta_0 = 10\)、\(\beta_1 = 0\)、\(\beta_2 = 5\)产生完全相同的预测（和残差 ），但会非常接近。除\(1\)月\(1\)日外，其他日子的预测结果完全一样。此时，和完全多重共线性不同，我们能选出一组 “最佳” 估计值。只需把\(\beta_0\)设为 \(1\)月\(1\)日的\(Sales\)值，\(\beta_1\)设为冬季其他日子与\(1\)月\(1\)日\(Sales\)的平均差值，\(\beta_2\)设为非冬季日子与\(1\)月\(1\)日\(Sales\)的平均差值就行。
然而，我们选定的那些估计值完全取决于\(1\)月\(1\)日的数据。要是 \(1\)月\(1\)日的销售额是\(5\)，那我们会得到\(\beta_0 = 5\)、\(\beta_1 = 5\)、\(\beta_2 = 10\) 。要是\(1\)月\(1\)日销售额碰巧是\(15\)，所有系数就会剧烈变化，变成\(\beta_0 = 15\)、\(\beta_1 = -5\)、\(\beta_2 = 0\) 。单个观测值改变数值，每个系数的变动量差不多就等于这个观测值的变动量！这样一来，\(\beta_1\)的波动程度和单个销售额观测值的波动程度一样。\(\beta_1\)的标准误本质上就是销售额的标准差，而不是像通常那样随着样本量增大快速变小117 。
那面对共线性，我们能做些什么呢？和往常一样，第一步是仔细思考我们的模型。甚至在考虑计算出的相关性是否真的很高之前，问问自己 “我为什么要把这些预测变量纳入模型？” 记住，给模型添加另一个预测变量，意味着在估计其他变量的效应时，要对这个变量进行控制。你真的想控制它吗？
当你有同一概念的多个度量指标时，一个诱人的想法是把它们都纳入模型，这样它们能以某种方式相互补充，或者累加出效应。但实际上，这只会迫使每个指标去体现其与其他指标无关的那部分变异的效应。
回到情商测试两部分的例子，为什么我们要在控制第二部分的情况下，了解第一部分的效应呢？“在控制考试第二部分的情况下，第一部分每多一分，会增加……” 这到底是什么意思呢？第一部分的系数代表的是，去除第二部分所能解释的关系后，第一部分分数和升职概率之间的关系，但实际上，第二部分几乎能解释全部关系。这样一来，剩下的变异就很少，会让\(\text{var}(X)\)很小，\(\sigma^2/\text{var}(X)\)很大，从而夸大标准误。不管怎样，当相关性真的那么强时，你可能并没有关闭任何混淆路径，只是把你想要的处理效应给控制没了。
在这种有同一事物多个度量指标的情况下，一种在纳入所有预测变量的同时避免共线性的好方法，是使用降维方法。这类方法会将同一事物的多个度量指标提炼出它们的共同部分，捕捉你认为它们所代表的潜在概念（而不是像把它们都纳入回归那样，关注非共同部分 ）。
降维是一个很大的研究领域，主要源自心理测量学，我在这儿就不试图充分展开讲了。不过，像潜在因子分析和主成分分析这类方法，会接收多个变量，生成一小组因子，这些因子代表变量不同的共同方面。它们要么寻找彼此无关的成分（这样就完全不会产生共线性 ），要么专注于提取那些可能相互关联的潜在概念。不出所料，机器学习研究者经常得一次性处理成千上万的变量，近年来也为降维做出了很多贡献。想了解综述可看 Jia 等人（2022 年 ）的研究 。
那在其他情境中，共线性问题并非源于同一事物的多个度量指标，而是不同事物的多个度量指标恰好高度相关，我们该怎么处理呢？
比如，你可能想研究工资上涨对企业盈利能力的影响，同时控制在巴黎和布鲁塞尔（企业在这两个市场频繁交易 ）的股市指数。巴黎和布鲁塞尔的股价可能高度相关，但它们不是一回事。我们能想到要同时控制这两个指数，但也许相关性太高，会产生共线性问题。是该去掉其中一个，还是用降维方法？多高的相关性算 “太高” 呢？
你可以去检查所有预测变量之间的相关性，找出高度相关的情况，但这没法涵盖那种需要多个变量才能对另一个变量形成非常接近的线性预测的情况。
一种用于检查是否存在过度共线性的常用工具（反正在线性回归中是这样 ）118，它能一次性考虑所有变量，这就是方差膨胀因子\(（VIF）\)。模型中每个变量的\(VIFVIF\)都不同。具体而言，对于变量\(j\)，其\(VIF_j = 1/(1 - R_j^2)\)，其中\(R_j^2\)是将变量\(j\)对模型中所有其他变量进行回归得到的\(R^2\) 。
大致的经验法则是，要是一个变量的\(VIF\)超过\(10\)（也就是，把该变量对其他预测变量做回归得到的\(R^2\)超过\(0.9\)），那你可能得考虑把它从模型里去掉，或者和与之强相关的其他变量一起做降维处理。也有人把\(VIF\)超过\(5\)作为临界值。总体而言，在统计学里，用硬性临界值作为决策规则通常不太合理。相反，把高 V\(VIF\)当作仔细审视模型的契机，思考去掉某些变量是否合理。
在\(R\)中，你可以用\(car\)包的\(vif(m1)\) 函数计算模型中每个变量的\(VIF\)，其中\(m1\)是你的模型。在\(Stata\)中，先运行回归，然后用\(estata \space vif\)命令。
\(Python\)里，你可以从\(statsmodels.stats.outliers\_influence\)导入\(variance\_inflation\_factor\)函数。不过，它的用法和\(R\)、\(Stata\)版本有点不同，因为它是对预测变量的数据框操作，而非对模型操作，而且一次只能计算一个变量的\(VIF\)，所以要是你想算所有变量的，必须循环处理。要是\(x\)是仅包含预测变量的数据框，你可以用\([variance\_inflation\_factor(x.values, i) \space for \space i \space in \space range(x.shape[1])]\)计算每个变量的\(VIF\)。
13.4.3 测量误差
有个看起来很合理的假设，就是数据里的变量就是数据里呈现的变量。但要是…… 它们并非如此呢？
我的意思是，你在数据中记录的值往往并不完全精确，这种情况相当常见。而且，数据里的变量常常只是你想要研究的变量的替代指标，而非真实的变量本身。
不精确是生活的自然组成部分。没有哪种测量是完全精确的。要是你有一个身高数据集，数据里的 “\(5英尺6英寸\)” 和 “\(6英尺1英寸\)”，很可能隐藏着更精确的值，比如 “5.5132161 英尺” 和 “6.0781133 英尺”，只是你无法获取这些更精确的值。要是唯一的问题是你把数值四舍五入到了合理的精确程度，那可能不算大问题。但要是测量误差更严重 —— 比如你把两个差异很大的身高都四舍五入成\(6\)英尺，或者那个 “\(5英尺6英寸\)” 实际是 “\(5英尺8英寸\)”，但有人记录错了，或者量身高时尺子用错了，那么数据里的数值就和真实数值不匹配，我们的估计也就无法捕捉到数据中实际存在的数值了。
替代指标的使用也很常见。替代指标是用来替代另一个变量的变量，因为你没有真正想要的那个变量。要是你想了解给学生配备笔记本电脑对未来考试成绩的影响，你可能觉得数学能力是一个混淆因素，所以想控制学生的数学能力。你实际上看不到他们的数学能力，但能获取他们数学测试的分数，这虽不完全等同于数学能力，但应该密切相关。测试分数是数学能力的替代指标，可并非一回事。要是两个学生能力相同但测试分数不同，模型会把他们当作不同个体处理，可实际上他们是一样的，这样模型就会出错。
替代变量：应与模型中实际想要的变量高度相关的变量，因无法获取实际想要变量的测量值，故用其替代 。
不精确和替代指标都属于测量误差，也叫变量误差（errors-in-variables ）。意思是，模型里的变量代表你想要的真实潜在值，但带有一些误差。\(X = X^* + \varepsilon \quad (13.18)\) 这里，\(X\)是你数据里的变量，\(X^*\)是潜在的未观测变量，\(\varepsilon\)是误差项119。
测量误差：因测量错误、不精确或使用替代变量，你所使用的变量无法精确代表你实际想要的变量 。
潜在值 / 潜在变量：是模型的一部分，但并非数据的一部分的值或变量。一般来说，只有当你试图用其他变量来表示该潜在变量时，才会使用这个术语
我们为何要关注测量误差呢？不难想象测量误差会如何让我们的估计偏离正轨。在极端情况下，比如用某个完全不相关的变量当作替代指标。要是我想了解笔记本电脑对考试成绩的影响，却不用 “有笔记本电脑” 这个变量，而用 “头发是棕色”，那我们没理由认为这样能得出有笔记本电脑带来的影响。
但在没那么极端的情况下呢？主要有两类测量误差，它们会产生不同影响。
最容易处理的测量误差是经典测量误差，即误差\(\varepsilon\)与潜在变量\(X^*\)无关的测量误差。出现这种情况时，测得的变量\(X\)就是 “真实潜在变量加上一些噪声” 。
经典测量误差有明确影响：在回归\(Y = \beta_0 + \beta_1 X\)中，若\(X\)存在经典测量误差，那么估计值\(\hat{\beta}_1\)会比用无误差的\(X^*\)（即运行\(Y = \beta_0 + \beta_1 X^*\) ）得到的真实\(\beta_1\)更接近\(0\) 120。它会 “向零收缩”，也叫 “衰减” 。
为啥会这样呢？\(X^*\)对\(Y\)的效应是\(\beta_1\)，对吧？每当\(X^*\)变动，我们预期\(Y\)也会变动。但要是我们观测的是\(X\)（\(X^*\)带噪声的版本 ），就看不到\(X^*\)所有的变动了。我们会看到\(Y\)因\(X^*\)变动（这些变动在\(X\)中看不到 ）而波动。而且，\(Y\)中与\(X\)无关的变动越多，\(X\)和\(YY\)看起来就越不相关。所以估计值会向零收缩。
经典测量误差也有一些明确的应对办法，或者说，它给了一种不 “修正” 它的理由。最常见的是，人们会直接用普通最小二乘法，不做实际修正，然后把估计值当作 “下限”。他们知道真实效应至少和得到的估计值一样远离零，只是不知道大多少。
你可能注意到，我只讨论了\(X\)的测量误差。那\(YY\)呢？事实证明，\(Y\)的经典测量误差不算大问题。当然，它会降低模型的\(R^2\)（因为\(Y\)中出现了更多无法预测的噪声 ），但至少平均来看，不会影响系数。这些系数本来就已经适应\(Y\)中无法预测的噪声了，现在只是噪声更多了一点。这些噪声会被归入回归模型的误差项。
那什么是非经典测量误差呢？这是指误差与真实值相关的情况，会让问题复杂得多。为啥测量误差会是非经典的呢？比如，若\(X^*\)是二元变量，\(X\)也是（只能取\(0\)或\(1\)）。若\(X^* = 1\)，误差只能是\(0\)或\(-1\)；若\(X^* = 0\)，误差只能是\(0\)或\(1\)。\(X^*\)的值决定了误差的可能取值，这就使误差和真实值相关了。
非经典测量误差在连续变量中也会出现。比如，你有一堆企业主的税务数据，要衡量收入。你知道，有些人会在报税时虚报收入，好少缴税。以现金交易为主的企业，比主要用信用卡或支票收款的企业，虚报收入容易得多。而且，现金交易型企业的平均收入往往也比其他企业低。所以，我们会预期，低收入企业的负向测量误差（收入少报 ）比高收入企业多。此时，误差和真实值相关。
再举个例子，假设有项研究，变量是人们报告的每日运动量。经常锻炼的人，可能愿意如实跟调查人员说；而很少锻炼的人，可能会把数字往高了报，让自己看起来更好。
非经典测量误差比经典测量误差更难预测，因为要是不了解误差和\(X^*\)具体的关联方式，我们就不知道\(\beta_1\)会向零收缩、远离零膨胀，还是向上偏、向下偏，甚至完全不偏！啥情况都有可能。而且，和经典测量误差不同，非经典测量误差不管出现在\(Y\)还是\(X\)中，都是问题。
所以，要是你觉得存在测量误差，得好好想想，这误差有没有可能和潜在变量相关。要是相关，那你面临的问题比想象中更棘手。
要是我们存在测量误差，还想修正它，该咋做呢？要是你掌握了一些关于测量误差本身的信息，或许能做些调整，让结果更接近真实的\(\beta_1\)。我就不深入讲这些方法啦，你得去别处找资料。测量误差属于那种一深入研究就会涉及很深内容的课题。我能找到的每本聚焦该主题的教材，都比你现在读的这本高深多了，但像韦恩・富勒（Wayne Fuller ）2009 年写的《测量误差模型》（Measurement Error Models ） 这类书，看看或许还是值得的。
不过，咱别深入探讨（或者说，先把你之后能谷歌搜索的关键词列一列 ），能做些啥呢？要是你知道回归误差方差和测量误差方差的比率，可进行德明回归\(（Deming \space regression）\)或总体最小二乘法\(（Total \space Least \space Squares）\)，以此修正向零收缩的问题。
还有其他办法你可以去了解，包括广义矩方法\(（GMM）\)。另一种常用方法是工具变量法，不过和\(第19章\)里工具变量用于识别因果效应不同，这儿是用来调整测量误差的。要是你有同一潜在变量的两种测量方式，且两种测量的误差互不相关（比如，有两个不完美的温度计，在相隔一英里的地方记录温度，你想知道该区域的温度 ），你可以用一种测量作为工具变量去预测另一种，然后在模型里用预测值。这样能分离出两种测量的共同部分，这部分会比单独一种测量更接近真实潜在值。
这还没涉及专门为二元预测变量测量误差设计的那整类方法呢！这类方法里，很多都会尝试估计\(0\)变\(1\)、\(1\)变\(0\)的概率，再把这些纳入回归模型估计中。
当然啦，我列的这些方法都和线性模型里的测量误差有关。要是你用的是非线性模型，像\(logit\)模型、\(probit\)模型或者其他模型，咋办呢？嗯…… 你得去找专门为那些方法设计的测量误差修正变体。这些方法是存在的。通常它们是我已经讲过的模型的变体，比如广义矩方法\(（GMM）\)或者工具变量法。祝你探索愉快！
13.4.4 惩罚回归
本书前半部分一直在反复强调：依据理论构建模型。运用你所掌握的知识，极为谨慎地判定哪些控制变量是必要的，哪些应被排除，诸如此类。
然而，这一过程常常会给我们留下太多候选控制变量，这种情况并不罕见。要是有三十个、五十个乃至一千个潜在控制变量，都有可能帮我们阻断混淆路径，那该如何是好？即便我们有充分理由把它们全都纳入模型，实际这么做的话，也会陷入统计乱局，极难解释，还会陷入共线性噩梦。此外，数据很可能会过拟合—— 有这么多预测变量，模型会非常贴合样本，但对噪声很敏感（部分原因和共线性类似 ），这样就无法很好地契合真实模型，应用到新样本时预测效果也不佳。
我们或许得舍弃一些控制变量。但舍弃哪些呢？理想情况下，我们得有一套基于原则的模型选择流程，帮我们做抉择。
于是就有了惩罚回归。惩罚回归是源自机器学习领域的工具。机器学习方法总体而言，相较于因果推断，更关注样本外预测 121，但我们能在这儿用上它们的一个技巧。还记得普通最小二乘法是通过最小化残差平方和来选取系数的吧？我们可以把这个目标表示为：
\[ argmin_\beta \{\sum(Y-\hat{Y})^2\} \quad (13.19) \]
换句话说，就是选取能使\(\sum (Y - \hat{Y})^2\)（即残差\((Y - \hat{Y})\)的平方和 ）最小的\(\beta\)（参数 ）。要注意，\(\beta\)会出现在\(\hat{Y}\)里，\(\hat{Y}\)的表达式是\(\hat{\beta}_0 + \hat{\beta}_1 X\)（再加上其他可能有的系数和预测变量 ）。
惩罚回归的思路是：“行，最小化残差平方和是挺好。但我希望你同时达成两个目标。一方面，尽量把残差平方和变小，另一方面，也得让这个关于\(\beta\)系数的函数变小。” 此时，目标就变成了：
\[ argmin_\beta \{\sum(Y - \hat Y)^2 + \lambda F(\beta) \} \quad (13.20) \]
其中，\(F(\cdot)\)是一个函数，它接收所有系数122，以某种方式把它们汇总起来，通常是随着\(\beta\)的绝对值增大，这个函数值也变大 ；\(\lambda\)是惩罚参数 —— 很小的\(\lambda\)意味着 “别太在意惩罚项，就尽量把残差平方和最小化就行”，而大的\(\lambda\)意味着 “把\(\beta\)系数变小很重要，比起最小化残差，更得操心这个” 123。
你可以为\(F(\cdot)\)选很多函数，不过最常用的是 “范数” 函数，具体来说，\(L1\)范数会给出 “\(LASSO\)回归”，\(L2\)范数会给出 “岭回归”。\(L1\)范数函数就是把\(\beta\)的绝对值加起来（\(\sum |\beta|\)，其中\(|\cdot|\)是绝对值函数 ），\(L2\)范数函数则是把\(\beta\)的平方加起来（\(\sum \beta^2\) ）124 。
引入这个惩罚项，会促使估计量选择更接近\(0\)的\(\beta\)值，以此避免惩罚。为啥这样做可取呢？因为没那么激进的预测，在拟合当前数据时会有所保留，但作为交换，你很可能对未来数据有更好的预测，因为你的估计不会被样本里的噪声过多干扰。
但这些方法不只是统一收缩系数；有点像用系数 “预算”，而且花得很明智。只有当系数真的有助于降低残差平方和时，才会让系数变大。\(L2\)范数（岭回归 ）最终会得到收缩后的系数。而\(L1\)范数\(（LASSO）\)版本，也叫 “最小绝对收缩和选择算子”，不只是收缩系数，还常常把很多系数变成恰好为\(0\)，实际上就是把那些变量完全从模型里剔除了 —— 这就是为啥它是 “选择算子”。
这样，我们就能把\(LASSO\)用到协变量过多的问题上。它会找出哪些变量可以剔除，同时对残差平方和的损害最小。然后，还会进一步收缩保留下来的系数。
但要是我们想识别因果效应，\(LASSO\)会不会给我们一个糟糕的模型呢？\(LASSO\)可能在两方面让我们陷入麻烦。第一，有些协变量即便对预测帮助不大，但为了阻断混淆路径，也需要纳入模型。这很好解决 —— 就告诉\(LASSO\)别剔除特定的必要协变量。或者，在它建议剔除之后…… 别照做，再把这些变量加回来。
\(LASSO\)可能让我们陷入麻烦的第二个方面是，我们的目标和机器学习领域人们追求的（也是\(LASSO\)设计初衷 ）很不一样，我们不是要得到最佳的样本外预测，而是要尽可能精准估计因果效应。可\(LASSO\)呢，明确地说，甚至都没尝试去得到\(\beta\)的最佳估计 —— 它收缩系数的操作，并不会让系数成为更准确的估计，只是让样本外预测更好。实际上，\(LASSO\)通常会给出有偏的系数估计。不过这也很好解决。要是你只用\(LASSO\)做变量选择，那就选。之后，用它选中的变量去预测处理变量或结果变量，再用这些变量重新跑常规的无惩罚\(OLS\)（贝洛尼、切尔诺朱科夫和汉森，2014 ）。
我们实际该如何执行\(LASSO\)呢？运行\(LASSO\)时，有几个实际要点得牢记。
这里要指出的一个重要点是，在估计\(LASSO\)模型前，对所有变量进行标准化很关键。也就是，对每个变量，减去其均值，再除以标准差。
为啥要这么做呢？因为惩罚项\(\sum |\beta|\)只关心系数的大小。而系数大小对变量的度量尺度很敏感！比如，对比模型\(Y = \beta_0 + \beta_1 X\)和\(Y = \beta_0 + \beta_1 \text{Halfof}X\)（其中\(\text{Halfof}X = X/2\) ），要是在第一个模型里估计出\(\hat{\beta}_1 = 3\)，在第二个模型里\(\hat{\beta}_1 = 6\)的话，得到的预测结果完全一样。\(\text{Halfof}X\)里的 “除以 2”，刚好被\(\beta_1\)翻倍抵消了。常规\(OLS\)里这不是问题，因为我们的解释会相应调整。\(\text{Halfof}X\)变动\(1\)个单位，对应\(X\)变动\(2\)个单位。所以\(\text{Halfof}X\)变动\(1\)个单位的影响，是\(X\)变动\(1\)个单位的两倍，系数是\(6\)而非\(3\)，也说得通。
变量尺度越小，系数越大，这就意味着，仅仅把\(X\)除以\(2\)（本应对模型无实质影响 ），会让\(X\)的系数受到的惩罚翻倍，也让\(LASSO\)更有理由收缩或剔除它。
我们不想这样，所以会先标准化，确保所有变量处于同一尺度。有些\(LASSO\)命令会自动做标准化（或有选项可开启 ），其他的则得你自己动手。
第二点要牢记，\(LASSO\)能让你大胆尝试，看看哪些变量需要保留。所以，运行\(LASSO\) 时，给几乎所有变量都加上一堆多项式项和交互项，这种做法挺常见。要是多项式项不重要，就会被剔除。
当然，这么做也不能太离谱。要是你真把所有变量两两交互，\(LASSO\)又选了些奇怪的交互项，你能理解模型吗？要是纳入一个多项式，结果它只让保留三次项，却把一次项、二次项扔了，你到头来还是得把低次项加回去。
第三，我们跳过了\(LASSO\)的一个重要步骤 —— 选定\(\lambda\)值125，它决定对系数施加多大惩罚 。你可以按需设定，选大值能剔除很多变量，选小值则只剔除少数。通常，人们会用交叉验证来选这个值。
交叉验证是把数据随机分成多块。选一个\(\lambda\)值，然后依次去掉每一块，用剩下的数据预测被去掉块的值，相当于做样本外预测。把去掉的块加回来，再去掉下一块，重新估计模型、预测被去掉块的值。重复此过程，直到所有块都被去掉过一次，评估所选\(\lambda\)下样本外预测的效果。接着换个\(\lambda\)值，重复整套流程。不断重复，最终选出让样本外预测效果最好的\(\lambda\)。
话说回来，咱来跑些\(LASSO\)模型吧126 。要是你运行代码，得到的剔除变量结果可能略有不同。为啥？因为\(LASSO\)估计过程存在随机性。要想结果完全一致，得设定随机种子 127。
# Load packages and data
library(tidyverse); library(modelsummary); library(glmnet)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 4.1-9df <- causaldata::restaurant_inspections
# We don't want business_name as a predictor
df <- df %>% select(-business_name)
# the LASSO command we'll use takes a matrix of predictors rather than a regression formula. So! Create a matrix with
# all variables and interactions (other than the dependent variable) -1 to omit the intercept
X <- model.matrix(lm(inspection_score ~ (.)^2 - 1, data = df))
# Add squared terms of numeric variables
numeric.var.names <- names(df)[2:3]
X <- cbind(X,as.matrix(df[,numeric.var.names]^2))
colnames(X)[8:9] <- paste(numeric.var.names,'squared')
# Create a matrix for our dependent variable too
Y <- as.matrix(df$inspection_score)
# Standardize all variables
X <- scale(X); Y <- scale(Y)
# We use cross-validation to select our lambda family = 'gaussian' does OLS (as opposed to, say, logit)
# nfolds = 20 does cross-validation with 20 "chunks" alpha = 1 does LASSO (2 would be ridge)
cv.lasso <- cv.glmnet(X, Y, family = "gaussian", nfolds = 20, alpha = 1)
# I'll pick the lambda that maximizes out-of-sample prediction which is cv.lasso$lambda.min and then run the actual LASSO model
lasso.model <- glmnet(X, Y, family = "gaussian", alpha = 1, lambda = cv.lasso$lambda.min)
# coefficients are shown in the beta element - . means LASSO dropped it
lasso.model$beta## 9 x 1 sparse Matrix of class "dgCMatrix"
## s0
## Year -0.0496933820
## NumberofLocations -1.1281884176
## WeekendFALSE -0.0049024323
## WeekendTRUE .
## Year:NumberofLocations .
## Year:WeekendTRUE 0.0002087634
## NumberofLocations:WeekendTRUE -0.0055607783
## Year squared -0.0746412335
## NumberofLocations squared 0.9453048415# The only one it dropped is Year:NumberofLocations
# Now we can run OLS without that dropped Year:NumberofLocations
m1 <- lm(inspection_score ~ (.)^2 - Year:NumberofLocations + I(Year^2) + I(NumberofLocations^2), data = df)
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 83262.170*** |
| (5446.779) | |
| Year | -82.613*** |
| (5.419) | |
| NumberofLocations | -0.082*** |
| (0.001) | |
| WeekendTRUE | -105.326 |
| (119.606) | |
| I(Year^2) | 0.021*** |
| (0.001) | |
| I(NumberofLocations^2) | 0.000*** |
| (0.000) | |
| Year × WeekendTRUE | 0.053 |
| (0.059) | |
| NumberofLocations × WeekendTRUE | -0.008 |
| (0.007) | |
| Num.Obs. | 27178 |
| R2 | 0.216 |
| R2 Adj. | 0.216 |
| AIC | 170190.9 |
| BIC | 170264.8 |
| Log.Lik. | -85086.453 |
| RMSE | 5.54 |
import pandas as pd
import numpy as np
import statsmodels.formula.api as sm
from sklearn.linear_model import LassoCV
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import Pipeline
from causaldata import restaurant_inspections
df = restaurant_inspections.load_pandas().data
# Create a matrix with predictors, interactions, and squares
X = df[['NumberofLocations','Year','Weekend']]
transform_pipeline = Pipeline(
steps=[
# Standardize the data
('scaler', StandardScaler()),
# and add polynomials
('poly', PolynomialFeatures(2, interaction_only = False,include_bias = False))]
)
X_all = transform_pipeline.fit_transform(X)
# Use LassoCV to pick a lambda with 20 chunks.
reg = LassoCV(cv = 20, max_iter = 10000,).fit(X_all, df['inspection_score'])
reg.coef_## array([-5.56235824e+00, -7.15979525e-01, 0.00000000e+00, 8.41441424e-01,
## 2.01207410e-01, -1.73728423e-02, 7.15930183e-01, 3.43634122e-02,
## 1.22796760e-04])# Looks like Weekend gets dropped. So we can redo OLS without it
m1 = smf.ols(formula = '''
inspection_score ~ NumberofLocations + I(NumberofLocations**2) + Year +
I(Year**2) + Weekend:Year + NumberofLocations:Weekend + NumberofLocations:Year''', data = df).fit()
m1.summary()| Dep. Variable: | inspection_score | R-squared: | 0.217 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.217 |
| Method: | Least Squares | F-statistic: | 1077. |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.00 |
| Time: | 10:58:32 | Log-Likelihood: | -85070. |
| No. Observations: | 27178 | AIC: | 1.702e+05 |
| Df Residuals: | 27170 | BIC: | 1.702e+05 |
| Df Model: | 7 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 8.274e+04 | 5438.956 | 15.212 | 0.000 | 7.21e+04 | 9.34e+04 |
| NumberofLocations | -0.8998 | 0.142 | -6.349 | 0.000 | -1.178 | -0.622 |
| NumberofLocations:Weekend[T.True] | -0.0066 | 0.007 | -0.956 | 0.339 | -0.020 | 0.007 |
| I(NumberofLocations ** 2) | 0.0001 | 1.79e-06 | 66.684 | 0.000 | 0.000 | 0.000 |
| Year | -82.0679 | 5.412 | -15.165 | 0.000 | -92.675 | -71.461 |
| Weekend[T.True]:Year | 0.0001 | 0.000 | 0.559 | 0.576 | -0.000 | 0.001 |
| I(Year ** 2) | 0.0204 | 0.001 | 15.136 | 0.000 | 0.018 | 0.023 |
| NumberofLocations:Year | 0.0004 | 7.05e-05 | 5.772 | 0.000 | 0.000 | 0.001 |
| Omnibus: | 2528.009 | Durbin-Watson: | 1.942 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 3426.234 |
| Skew: | -0.773 | Prob(JB): | 0.00 |
| Kurtosis: | 3.799 | Cond. No. | 6.55e+11 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.55e+11. This might indicate that there are
strong multicollinearity or other numerical problems.
causaldata restaurant_inspections.dta, use clear download
* NOTE: this code relies on functions introduced in Stata 16.
* If your Stata is older, instead do ssc install lars and use
* the lars function (the syntax will be a bit different)
* Standardize all variables
foreach var of varlist inspection_score numberoflocations ///
year weekend {
summ `var'
replace `var' = (`var' - r(mean))/r(sd)
}
* Use the lasso command to run LASSO while including all interactions
* and squares, using sel(cv) to select lambda using cross-validation
* "linear" gives us OLS, but logit or probit would work
lasso linear inspection_score c.numberoflocations##c.year##c.weekend c.numberoflocations#c.numberoflocations c.year#c.year, sel(cv)
* get list of included coefficients
lassocoef
* Looks like the only one not included is numberoflocations*weekend, so we can rerun OLS without that.
* first, reload the non-standardized values
causaldata restaurant_inspections.dta, use clear
reg inspection_score c.numberoflocations##c.year c.year##c.weekend c.numberoflocations#c.year#c.weekend ///
c.numberoflocations#c.numberoflocations c.year#c.year要是你阅读三种语言的代码，可能会注意到，不同语言的结果有差异。尽管给的是完全相同的数据（\(Stata\) 版本除外，它还额外给了三个变量间的交互项 ），但每种语言的\(LASSO\)估计在筛选剔除变量时都不一样。在统计软件里，这并不像你想的那么罕见 —— 不同语言的结果时常不同（麦卡洛和维诺德，2003 ） ，尤其是像\(LASSO\)这类尚未 “解决” 的方法，得在一堆不同系数里搜索可行解，不同软件包的搜索方法可能大相径庭。
而且，默认选项也可能差异很大。就算底层估计代码一样，不同语言基于估计搜索从哪开始、如何做交叉验证、交叉验证时检查哪些\(\lambda\)值这类问题，仍可能产生不同结果。不管用哪种语言，一定要读文档，确保清楚自己让软件估计的是什么，最好在多种语言里试试估计，看看是否一致。
说我们“调整”这些变量，而不是控制它们，可能更准确，因为我们实际上并不像在受控实验中那样对任何东西拥有控制权。但“控制”是最常见的说法。↩︎
这是一种纯粹的统计学意义上的“解释”，并不说明\(Y\)发生的原因（除非我们已经确定了\(X\)对\(Y\)的因果效应），而仅表示用\(X\)预测 YY 的能力。↩︎
还记得第3章的内容吗？↩︎
有趣的是，这并不要求\(X\)与\(\varepsilon\)完全无关，只需不存在线性相关即可——二者之间可能存在某种关系，只要不是线性关系。例如，若\(X\)与\(\varepsilon\)的关系呈现V形或U形，它们显然存在关联，但描述这种关系的直线将是水平的（斜率为零），因为下降部分会被上升部分抵消。实际上，两个变量存在关联却完全无线性相关的情况十分罕见，因此我们通常无法借此忽略\(X\)与\(\varepsilon\)之间的关联性。↩︎
“内生性”与”外生性”术语源自生物学词汇”内源”与”外源”，大致分别表示”来自系统内部”和”来自系统外部”。内生变量由系统内其他因素（即影响结果的其他决定因素）所决定。↩︎
其他领域对此还有不同表述方式。也可以说\(A\)是”未观测/未测量的混杂因素”，或者说\(X\)与\(Y\)存在”伪相关”。↩︎
以下计算将基于前文所述的外生性假设成立为前提↩︎
既然无法直接观测误差项\(\varepsilon\)，我们如何确定其标准差\(\sigma\)？实际上我们无法精确获知\(\sigma\)，因此对抽样分布的认识将基于对\(\sigma\)的最佳估计。我们使用可观测的残差替代不可观测的误差，并通过残差的标准差来估计\(\sigma\)的值。↩︎
若采用贝叶斯视角，则该初始值需要大幅修正。↩︎
假设检验虽存在诸多争议（如依赖人为设定的显著性阈值——\(p = 0.049\)拒绝原假设而\(p = 0.051\)则不拒绝；以及将”拒绝明显错误假设”等同于”发现真实效应”的逻辑缺陷），但由于其在社会科学研究的普遍性，研究者必须理解其运作机制。需注意的是，社会科学中真实效应恰好为零的情况极为罕见，因此只要样本量足够大，我们总能拒绝原假设。建议：可将显著性检验作为学术交流工具，但不应过度依赖——切勿基于显著性选择模型（本书前半部分已阐述更优的模型选择方法）。↩︎
注意这里将百分位数乘以\(2\)——因为”相同或更大差异”包含分布的两侧。这是标准的”双尾检验”。若进行”单尾检验”（仅关注同方向的差异），则不需要乘以\(2\)。↩︎
为什么说它是假 “阳性” 呢？因为我们在某种程度上，会把假设检验中拒绝原假设当作一种 “成功”。当然，这很愚蠢 —— 我们是想找到真实值，而不只是到处去拒绝原假设。一次合理估计后的不拒绝，比一次糟糕分析后的拒绝原假设要有价值得多。但人类心理中似乎存在某种因素，让我们表现得好像并非如此。↩︎
要是《指环王》讲的是统计学，显著性检验就是那枚至尊魔戒。更别提还有那么多人想把它扔进火山里了。当然啦，如果《指环王》真的是关于统计学的，他们也就不会拍那些电影了。↩︎
为啥不行呢？首先，可能真的没什么效应。不显著不代表错误呀！而且，为了追求显著性而改变分析，会让你未来的显著性检验变得不正确、毫无意义 —— 这些检验都假定你不会这么做。你的结果会变成垃圾，哪怕它们看起来不错，你还能自欺欺人地觉得它们挺好。↩︎
公平地说，有时候研究人员也会这么干。这叫 “\(p\)值篡改”，就是不断改变统计分析，直到弄出一个统计显著的结果。这很糟糕。这是件坏事，因为这么做本身就不好，而且会让结果变糟。太糟了。↩︎
话虽如此，猜测一个不显著的结果，如果样本量变大（比如大得多的样本 ）是否会变得显著，也不是完全没用。↩︎
尽管包括我在内的研究人员，常常会把 “显著” 当作 “统计显著” 的简称。↩︎
此数据来自 Kaggle 平台上的 Louis - Ashley Camus 。↩︎
有些回归表格会把精确性度量放在系数右侧，而非下方，这样每个回归就会有两列（系数列和精确性度量列 ）。↩︎
如果没说明，通常就是标准误；如果是单个数值，几乎肯定是标准误或者\(t\)统计量。你可以通过看显著性星号来区分这两者 —— 如果与系数相比，精确性数值很大却给出显著性，那很可能是\(t\)统计量；如果与系数相比，小数值给出显著性，那很可能是标准误。↩︎
记住，p值是得到与原假设值距离一样远（或更远 ）结果的概率。如果p值低于\(\alpha\)，就具有统计显著性。↩︎
在许多领域，且在大多数 R 命令默认设置中，所有（星号对应临界值 ）都整体后移一步。\(*\)表示小于 0.05，\(**\)表示小于 0.01，\(***\)表示小于 0.001。这种情况下，通常还会有另一种指示符，比如\(+\)表示 0.1（的临界值 ）。↩︎
这有点本末倒置 —— 我们通常想在开始前选定\(\alpha\)值。知道在更高\(\alpha\)值下的显著性检验有啥意义呢？我们之前说过，那不算数。要是我们的\(\alpha = 0.05\)，那么\(***\)也不会给我们额外信息 —— 某事物要么显著，要么不显著。不存在 “更显著” 一说。当然啦，我描述的是一种理想化的显著性检验情况，人们实际上并不会这么做。↩︎
这里常见的其他内容可能是一两个信息准则（比如 AIC 或 BIC，分别是赤池信息准则和贝叶斯信息准则 ），或者平方和。我不会在这里讲这些。↩︎
没错，往模型里加一个随机变量确实会解释更多变异。实际上，往模型里加任何变量，都会让\(R^2\)上升一点，哪怕这个变量毫无意义。↩︎
我不懂这话啥意思（作者调侃，用于形容学生因执着\(R^2\)而纠结的状态 ）。↩︎
有时回归结果会以 “标准化系数” 呈现，这时表述就变成 “…… 的一个标准差变化” 。不过你也可以简单把它想成是变量经过缩放、标准差为 1 时的 “一个单位变化” 。↩︎
要是你的预测变量数据的相关范围离 0 很远，常数项有时会出现奇怪的值。这不算问题 —— 别尝试对数据范围外的情况做预测就行。一般来说，我们不用太操心常数项。↩︎
也就是，关闭那些经过这些变量的后门路径。↩︎
还记得第 4 章讲的 “条件均值” 内容不？我这儿说的 “conditioning（以…… 为条件 ）” 就是那个意思。↩︎
或者说，至少我们让它的线性预测值保持恒定了。要是 “检查年份” 与 “检查分数” 或 “门店数量” 的关系更适合用曲线描述，那么在回归方程中把 “检查年份” 作为线性预测变量纳入，无法完全控制这种关系带来的影响。↩︎
我对下标没啥意见！但我觉得在我们当前的知识层面，下标带来的干扰大于帮助。↩︎
我知道
estout可能是制作 Stata 回归表格更好用的包，所以我在这里用它。我听说regsave也不错，但就我个人而言，我会一直用到死那天 。你们甭想从我这儿拿走它晦涩的语法，还有它比outreg更新滞后这让人挠头的事实（除非我死了 ）。↩︎有些人指出，这（某些情况 ）不应被当作二元变量处理。希望这些人知道他们给回归分析增添了多大难度。要为研究人员想想啊。↩︎
有时你会看到人们用 - 1 和 1，或者 1 和 2，而非 0 和 1。不过这通常只会造成混淆，而且在社会科学中相当少见。↩︎
为啥会这样呢？普通最小二乘法试图拟合一条使残差平方和最小的直线。但\(X\)轴上的数据只能取两个值：\(0\) 和 \(1\)。正因如此，\(OLS\)直线只能产生两个实际预测值：一个对应\(X\)轴为 \(0\) 时，一个对应\(X\)轴为 \(1\) 时。每种情况下能做出的最佳预测，就是使残差平方和最小的值，也就是以二元变量取 \(0\) 或 \(1\) 为条件时结果的均值。然后，直线的斜率就是当\(X\)轴变量增加 \(1\)（即从 \(0\) 变为 \(1\) ）时，预测值的变化量。↩︎
而且，如果冬季的销售额平均更低，\(\beta_1\)就会是负数。↩︎
有时这些类别可能会有重叠 —— 比如有些人确实住在多个国家。但希望我们能避免这种情况，或者这种情况极少，我们可以直接忽略，因为它们会让问题复杂很多。↩︎
这并不意味着我们要把法国的所有观测数据都去掉。只是说，我们要从回归模型里移除 “法国” 这个二元变量。↩︎
你能看出为啥吗？（提示：原本冈比亚比法国高 0.5，新西兰比法国高 3；若新西兰为参照，冈比亚与新西兰的差值就是 0.5 - 3 = -2.5 ）↩︎
这些也分别被称为线性项、二次项和三次项，四次方的则称为四次项（quartic terms ）。↩︎
注意，这里包含了直到三次方的所有项。你几乎永远都不想遗漏任何项，比如像\(\beta_1 X^2+ \beta_2 X^3\)种情况，它有二次项和三次项，但没有一次项↩︎
任何函数关系都必须满足一个基本前提——我们需要确保对于每一个\(X\)的取值，都能获得一个有意义的\(Y\)的条件期望值。↩︎
当然，这在线性模型中也适用。在 \(Y = \beta_0 + \beta_1 X\) 里，Y 对 X 求导就是 \(\beta_1\) 。这也就是我们对线性模型中 \(\beta_1\) 的标准解读。↩︎
就算你不懂微积分，只要面对的是多项式，求这样的导数也很简单。对每一项，用该项的指数去乘它，然后把指数减 1（要记住 \(X = X^1\) ，指数减 1 后就得到 \(X^0 = 1\) ）。所以 \(\beta_3 X^3\) 乘以指数 3 就得到 \(3\beta_3 X^3\) ，然后指数减 1 变成 \(3\beta_3 X^2\) 。对每一项都重复这个操作就行。↩︎
我可以告诉你，因为这数据是我自己生成的，真实的潜在模型其实确实有一个三次项。尽管如此，回归中的三次项几乎没起到什么作用。↩︎
需谨记，“统计不显著”并不实际意味着”该系数为零”，而是表明”我们尚未掌握足够证据证明其非零”。若将其从模型中剔除，则等同于强制设为零值——这种做法可能是错误的。↩︎
为何称其”自然”？因为只要涉及复合增长（即利息产生新利息，或动物繁衍后代后，后代继续繁衍），数字\(e\)就会如圆周率\(\pi\)之于圆和曲线的计算那样，必然出现在你的运算过程中。↩︎
“希望你的变量服从正态分布” ，不是转换变量的合理理由。有一种常见的误解是，回归中使用的变量需要服从正态分布。回归的标准误计算，关注的是误差项（error term） 是否正态分布（即便如此，只要误差表现良好，也不是超级关键 ），但对于变量本身，分布如何并不重要。↩︎
左偏指的是有很多高值，而低值极少。不过在社会科学研究中，左偏比右偏少见得多。↩︎
想要了解更多关于有偏态变量的内容，回到第3章查阅即可 。↩︎
异常值（Outliers）也可能出现在无偏态（non - skewed ） 的数据中 。↩︎
要是异常值没这影响，我们大可省掉转换这一堆事，直接把所有异常值都删掉就行（但实际不能这么做，所以说明异常值有影响是合理且必要的 ）。↩︎
\(ln()\) 指的是 “自然对数”，也就是以\(e\)为底的对数。在实证社会科学中，你几乎不会看到其他类型的对数。↩︎
不过，把它用在同时有正值和负值的变量上，可能达不到你想要的效果。↩︎
有些人把这个叫作 “归一化（normalizing）”，而非 “标准化（standardizing）”，但这有误导性，因为听起来像是在强行让数据服从正态分布。标准化和产生正态分布毫无关系，与普遍的误解相反。↩︎
或者，你可以把数据 “分箱（bin）”，创建像 “0 - 10 岁”“11 - 20 岁” 之类的类别。不过分箱要谨慎 —— 如果箱子太宽，而且箱子之间的分界点不是基于有意义的真实不连续性来选的，那这种方法其实也没那么灵活。↩︎
天啊！↩︎
要注意，本节所有的复杂性都源于我们试图解读有实际意义的大变化 ，比如问 “X变化 10% 会有什么影响” ，或者 “X变化 1 个单位，Y会有百分之几的变化” 。大多数时候，应用研究者就是这么做的。你也可以转而关注无穷小的变化 （或者，如果你是经济学家，就看弹性或半弹性 ），这种情况下，你可以直接求导，或者用近似方法，这样计算不会给你带来太多麻烦。当然，之后你还得把效应大小当作导数或斜率来解读，要是你不适应这种方式，把效应大小放到具体情境中解读就会更棘手。比如，你还能自信地说 “X变化 1 个单位，与…… 相关” 吗？你觉得这样更容易吗？取决于你自己！↩︎
记住，\(\ln()\)是以\(e\)为底的对数，所以\(\ln(e) = 1\) 。↩︎
这和X增加 5.271%完全等价 ，但研究者通常会简化，直接说\(\ln(X)\)增加c，对应X增加\(c \times 100\)% 。你可能想更精确些，不过对于大的变化，需要用\(e^c\) 。比如，\(\ln(X)\)变化 1 个单位，对应X变化\(e^1 - 1 \approx 1.718\)，也就是 71.8% 的增长，不是 100%。差得远呢！但对于小的百分比变化，这在社会科学中大多够用，近似效果也不错。一旦超过 0.1/10%，就开始有点难搞了。↩︎
你可以换对数的底数来算这个（比如，用以 1.1 为底的对数\(\log_{1.1}(X)\) ），这样\(\log_{1.1}(X)\)变化 1 个单位，就对应X变化 10%。\(\ln(X)\)变化 1 个单位，对应X变化无需近似，你只需选好百分比增长的大小就行。这种方法没普及开，但要是我能说了算，会把这个脚注内容放到正文中。要是你读的是本书第二版，就知道我没赢（编辑注：你现在读的就是第二版，我没赢 ）。↩︎
重要的是，这个过程也依赖小\(c\)近似，只能用于小的百分比变化，和其他近似一样。↩︎
而且，有大量研究就是围绕 “关系是否随其他变量取值变化” 来设计的。我们在第 4 章提到的 Oster 研究，就属于这类论文。↩︎
为啥还要单独纳入参与交互的变量（即Z ）呢？因为不然的话，\(X \times Z\)的系数，不仅体现交互作用，还会混进Z的直接效应，搞不清Y与X的关系随Z如何变化，以及Y与Z的效应。不行！所以，我们单独纳入Z，把Z的直接关系转移到另一个项，从交互项里剥离出来。↩︎
作为一个，自认为 “清晰讲解概念” 技能强到敢写一本好几百页教材的人，我最近意识到，我都不记得 “不会求导、忘了不用微积分工具（这些工具能轻松解决问题，而不是用更费劲的基础方法 ），怎么解一道基础数学题” 是啥滋味了。要是你看到一本署名是我的微积分教材，放心，那要么很烂，要么是代笔，要么又烂又是代笔 。↩︎
again，要是你不懂微积分，注意看我们是怎么处理含\(X\)的项的 —\(\beta_1 X\)和\(\beta_3 XZ\) ，把\(X\)去掉，就得到\(\beta_1\)和\(\beta_3 Z\) 。↩︎
本就该如此。毕竟，搞交互项的核心 ，就是让\(X\)的效应随\(Z\)变化呀。↩︎
反过来也成立 —— 它也表示，\(X\)增加 1 个单位时，\(Z\)对\(Y\)的效应会增强多少。要是觉得非得让这两件事数值一样（确实有点奇怪 ），但就就是导数的性质导致的。↩︎
要是你想狠狠锻炼一下 “解读肌肉”，可以试试解读三向交互，比如\(Y = \beta_0 + \beta_1 W + \beta_2 X + \beta_3 Z + \beta_4 WX + \beta_5 WZ + \beta_6 XZ + \beta_7 WXZ\) 。逻辑是一样的，就是多了一步。\(\beta_7\)表示，Z增加 1 个单位时，W增加 1 个单位对 “X增加 1 个单位带来的Y的变化” 的增强程度。↩︎
早说过，有交互项时，单个系数没那么大意义，得把所有情况结合起来看。↩︎
为啥第二列里效应一直是正的，第三列却是负的？是不是很奇怪？确实奇怪！这种奇怪的发现，会促使我们更深入地探究模型和数据，看看究竟发生了啥。↩︎
不过，这个差异看着不算大，而且在统计上也不显著。↩︎
这个解释用了二元交互，但同样的逻辑也适用于连续交互。↩︎
多数情况下如此。要是两个 “噪声” 在特定方式下相关，做差可能抵消部分噪声。↩︎
即便该线性函数有时包含对单个预测变量的变换或多项式，它们在方程中仍以线性形式呈现。\(Y = \beta_0 + \beta_1 X + \beta_2 X^2\) 是\(X\)和 \(X^2\) 的线性函数 。↩︎
此外还存在诸如非线性最小二乘法等方法，用于处理预测变量中的非线性关系，例如当模型形式为\(Y = \beta_0 \beta_1^X\)时的估计问题。↩︎
“解析解”指的是，你可以推导出一个函数，告诉你如何找到合适的系数。在普通最小二乘法中，它是通过微积分最小化误差平方和得到的方程。代入你的数据，它会给出精确正确的答案。很多这类非线性回归方法做不到这一点，所以，它们不是代入数据就能得到系数，而基本上得尝试一大堆不同的系数，直到找到效果还不错的系数。虽然它们实际操作时，比我描述的要聪明得多，但也并不总是能很好地发挥作用。↩︎
要是你忘了，函数是一种数学对象，接收若干输入，输出一个结果。所以，要是我们的函数是\(F(x) = x + 2\)，那么\(F(5)\)就是\(5 + 2 = 7\) 。↩︎
而且，当我们添加其他控制变量等情况时，它会继续和我们早已熟悉的那些模型一样 。↩︎
要是你打算通过代数运算来证明模型的某些性质，或者推导其边际效应，\(logit\)（模型）用起来要容易得多。在计算方面，估计\(logit\)模型通常也更简便。但一般情况下不用操心这个，至少在开始估计多项式模型之前不用。↩︎
这一点适用于任何类型的模型设定错误。运行非线性回归却没有纳入合适的多项式项或控制变量？这也会得出奇怪的结果。↩︎
即便在普通最小二乘法中使用多项式，虽然这会引入一条曲线，甚至可能看起来有点像\(logit\)曲线，但也无法以一致的方式很好地完成这项工作 。↩︎
严格来讲，普通最小二乘法没有 “链接函数”，尽管线性广义线性模型有。但我觉得这仍是一种清晰的解释方式，所以…… 就这么讲吧。↩︎
关键在于，样本中这种效应的变化，和我们在\(第10章\)探讨的、或者谈及交互项时的处理效应异质性，不是一回事。边际效应随\(P(Y = 1)\)变化，仅捕捉了效应变化中很狭隘的一个方面 —— 只是与在\(Y\)预测值的最大值或最小值附近 “没多少增长空间” 相关的部分。并非像 “某种处理对男性比对女性更有效” 这类情况。处理效应异质性那套内容，在概率单位\(（probit）\)/ 对数单位\(（logit）\)模型中也适用，但你得像在普通最小二乘法里那样，明确去处理它 。↩︎
统计学里从不缺那种实际上不代表任何特定个体的代表性数值 —— 每次我们用均值时，就会出现这种情况。但在这个例子里，还有其他问题，而且我们有更好的选择。“没人真的有\(2.3\)个孩子” 这类问题，还不足以让我们彻底放弃使用均值，毕竟均值在其他方面有无可估量的价值，但它足以让我们放弃均值处的边际效应，因为有更优的替代方法。↩︎
当然，没人要求你非得算均值 —— 对于效应分布，你想咋处理就咋处理。算中位数、标准差，啥都行。或者看看整个分布。但一般来说，人们在这儿会关注均值。↩︎
我想要强调的是，标准误真的不只是用来计算统计显著性的！它们能提供关于估计值抽样分布的信息，这对于运用观测数据了解理论 / 总体分布及关系的诸多不同方法而言，都很有价值。↩︎
一种常见的误解是模型中的变量必须服从正态分布，但这是错误的；我们并不需要这样。↩︎
在深入探讨之前，我想指出，标准误出错的方式远不止这些，修正它们的方式也不止这些。德尔塔法、克林斯基 - 罗布法…… 有超多修正方法要学。在这里，我仅向你们展示标准误出错方式的皮毛。要是你决定继续做研究，你得做好准备，未来会在标准误上犯下奇妙的错误，那些错误是你如↩︎
要是你更擅长直观学习，就拿支铅笔在手里 —— 每只手的高度代表样本中\(Y\)的条件均值，你握铅笔的角度就是\(\hat{\beta}_1\)。现在让双手各自独立上下移动。有时这会改变角度 —— 要是左手大幅上下移动，但右手不动，角度会有很大变化；要是双手一起上下移动，角度变化就小。这是常规的抽样变动。现在试着只让一只手上下移动，角度变化会快得多！这就是异方差性的一种形式。↩︎
它被称为 “三明治估计量”，是因为在该估计量的矩阵代数形式\((X'X)^{-1}(X'\Sigma X)(X'X)^{-1}\)中，起到缩放作用（且在后续能让观测值相关联 ）的矩阵\(\Sigma\)被 “夹” 在一堆\(X\)中间了。↩︎
回归分析的一些使用者，尤其是经济学家，往往假定异方差始终是个问题，所以默认使用异方差稳健标准误。我没法说他们错了 。
Long, J. Scott, and Laurie H. Ervin. 2000. “Using Heteroscedasticity Consistent Standard Errors in the Linear Regression Model.” The American Statistician 54 (3): 217–24.↩︎
在 “刀切法（jackknife）” 中，模型会被多次估计，每次去掉一个观测值，通过观察去掉的单个观测值对估计值的影响程度来计算方差。
MacKinnon, James G., and Halbert White. 1985. “Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties.” Journal of Econometrics 29 (3): 305–25.↩︎
自相关一致标准误只是时间序列数据从业者众多工具中的一种。这类工具在本书中几乎完全没提到。虽然时间序列分析在本书的一些地方会出现，比如\(第17章\)的事件研究，但要深入学习时间序列分析，你最好去看更传统的计量经济学教材。↩︎
为什么从稳健标准误入手呢？为什么不直接修正自相关，不管异方差性呢？因为本质上，自相关往往也会引入异方差性。要是误差是相关的，那么当一个观测值的误差小时，附近其他观测值的误差往往也小，这样就会产生低误差方差。但有时一个观测值的误差大，附近其他观测值的误差往往也大，就会产生高方差。有时方差低，有时方差高，这就是异方差性呀！↩︎
更详细地说，我们将去均值的时间序列与误差项相乘，计算其与滞后项的相关性。把这些相关性加总，给一期滞后的权重略低于 1，然后随着滞后阶数增加，权重递减。将这个总和乘以 2，再加 1，这就是调整因子。↩︎
给矩阵代数爱好者：还记得几节之前夹在所有\(X\)中间的那个\(\Sigma\)矩阵不？对于聚类标准误来说，这个矩阵是分块对角矩阵，每个分组在其分块内有不受限制的非零值，矩阵其余部分为\(0\)。↩︎
在模型中纳入组指示变量作为控制变量，和在该组层面进行聚类，二者有啥区别呢？首先，你可以两者都做，人们也经常这么做。其次，在模型里纳入指示变量，意味着组归属是结果的一个重要预测因素，还可能是一个重要的混淆因素；而聚类则意味着组归属与模型的良好预测能力相关。↩︎
它们甚至可用于解释面板数据中任何形式的自相关，在面板数据里，你有多个时间序列都在一个数据集中 —— 只需按时间序列进行聚类就行 。↩︎
就算你有很多聚类，可要是处理（干预）是按聚类分组分配的，而且你只有几个接受处理的聚类，也会出现这种情况（指前面提到的聚类标准误不够精确的情况 ） 。↩︎
我在这里讲的是自助法标准误，但自助法的概念可不只限于标准误，你可能会在其他地方看到它。它在机器学习领域肯定很流行。↩︎
根据数据的结构 —— 比如，在面板数据中，同一个体被多次测量，有多个观测行 —— 完全随机抽样可能就不合适了。例如，“聚类自助法” 会随机抽取数据块，可能是个体或组，而非单个观测值。↩︎
每次进行自助抽样时，由于过程中存在随机性，这个确切数值会略有不同。要是你做自助抽样且希望结果可重复，一定要在软件中进行抽样前设置随机种子 。↩︎
要记住，就实际数据而言，那些理论分布都只是近似情况 ↩︎
要是你好奇我为啥在上面的例子里把样本从\(4\)个换成\(1000\)个，原因就是这个（指前面提到的自助法在小样本下表现不好 ） 。↩︎
当试图在一种语言里复现用另一种语言得到的结果时，这会造成很多困惑 。↩︎
要是你看这章看烦了，想想我写这章的感受吧。这一章字数快超\(3\)万了 —— 都赶上中篇小说篇幅了！我的手指都酸了。但我得接着写。你们肯定会弄懂线性建模的。哎哟，等等，停一下。我的无名指掉啦 。↩︎
通常，权重会被标准化，让整个样本中所有权重的平均值为\(1\)。这样能避免估计结果出现问题，比如，不会算出一个整体上更大的方差，而只是算出一个能更好代表某些观测值的方差。而且，严格来说，这里的情况是每个变量的每个取值都要乘以\(\sqrt{w}\)。然后，当你把两个变量相乘时，比如计算协方差时的\(XY\)，或者计算方差时的\(XX\)，就会变成权重\(\sqrt{w^2} = w\) 。↩︎
那个分母（即\(7.143×70 + 16.667×30\) ）是样本中所有权重的总和，其作用是对权重进行缩放，使它们的平均值为\(1\) 。↩︎
在很多实际情况中，你可能只能获取班级平均分的数据。要是你有个体数据，通常你会想用个体数据 。↩︎
除非组内每个观测值确实完全相同↩︎
所以，在使用频率权重之前，一定要确保你的数据符合 “一组相互独立、完全相同的观测值” 这一情况。这样的例子可以是一个具有二元结果的随机实验，经过聚合的四个观测值分别是 “未处理，\(Y = 0\)”“未处理，\(Y = 1\)”“处理，\(Y = 0\)” 以及 “处理，\(Y = 1\)”，这些就真正代表了相互独立观测值的集合 。↩︎
不过，你也可以用 expand weight（扩展权重），采用和 R、Python 里类似的复制每个观测值的方法 。↩︎
实际上，回归软件会对模型进行估计，但只会通过剔除其中一个变量，从而消除完全多重共线性 。↩︎
要是它们（两部分分数）不相关，那这个测试很可能不咋地 。↩︎
这种特定结果基于这样一个非常特殊的例子：两类（变量）仅排除了一个观测值，但相同的基本原理适用于任何高共线性情况 。↩︎
反正，在线性回归中是这样 。↩︎
这个方程把误差项表示为加到真实值上，但它也完全可以以其他方式起作用，比如与真实值相乘 。↩︎
在含控制变量的模型中，情况也是如此 。↩︎
但也并非总是如此。也可查看第 22 章的部分内容 。↩︎
而且，这是一个仅关于\(\beta\)系数的函数 ——\(X\)不会在这里出现，\(Y\)也不会 。↩︎
我在此处介绍的是线性 / 普通最小二乘法（OLS）形式的惩罚回归，但你也完全可以对其他任何回归方法进行惩罚，比如逻辑回归（logit）或概率单位回归（probit） 。↩︎
你也可以对 L1 范数和 L2 范数加权后相加，把两者结合起来。这被称为 “弹性网回归” 。↩︎
有时，它不叫\(\lambda\)（兰姆达 ），而叫\(\alpha\)（阿尔法 ） 。↩︎
这段代码也可在 https://lost - stats.github.io/ 的 LASSO 页面找到 。↩︎
那我为啥不在这儿这么做呢？是为了演示不设随机种子会让代码无法复现。把种子设好呀！↩︎