15 模拟
15.1 已被证实
与统计研究不同，统计研究完全由至少有一点失实的内容构成，而统计学本身几乎完全是真实的。统计方法通常是通过数学证明来发展的。如果我们假设\(A\)、\(B\)和\(C\)成立，那么根据统计学法则，我们可以证明\(D\)是真实的。
例如，如果我们假设真实模型为 \(Y = \beta_0 + \beta_1X + \varepsilon\)，并且 \(X\) 与 \(\varepsilon\) 不相关1，那么我们可以从数学上证明，\(\beta_1\) 的普通最小二乘估计值 \(\hat{\beta}_1\) 平均而言将等于真实的 \(\beta_1\)。这就像证明存在无穷多个质数一样千真万确。与证明质数无穷的情况不同，我们可以（而且确实会）无休止地争论这些潜在的真实模型假设以及 \(X\) 与 \(\varepsilon\) 不相关的假设在我们的研究案例中是否适用。但我们都认同，如果这些假设适用，那么 \(\beta_1\) 的普通最小二乘估计值 \(\hat{\beta}_1\) 平均而言将等于真实的 \(\beta_1\)。这已经得到了证明。
在这本书里我所探讨的几乎所有方法和结论都像这样有数学证明作为支撑。我选择省略了实际的证明过程，但这些证明暗藏在字里行间2，也在那些引用文献里，我能想象你正快速略过这些引用呢。
为什么要费力对我们的统计方法进行数学证明呢？因为这些证明能让我们了解我们的估计方法是如何起作用的。普通最小二乘法背后的证明告诉我们它的抽样变异性是如何产生的，我们的估计要识别出因果效应必须满足哪些假设，以及如果这些假设不成立会发生什么等等。这极其重要！
所以我们有一个明显的问题。我还没教过你们如何进行计量经济学证明。很有可能你们不会去学习这些证明，或者至少不太可能擅长做这些证明，即便你们从事计量经济学研究。撰写证明实际上是开发方法的人要做的事，而不是使用这些方法的人要做的事。许多（或许是大多数）正在积极开展统计学研究的人，如果他们一开始就知道一些统计学证明的话，现在也已经忘得差不多了，只记得最基础的那些3。
但我们仍然必须了解我们的方法是如何运作的！这就需要用到模拟了。在本章的语境中，模拟指的是利用一个我们能够控制的随机过程（通常是一个数据生成过程）来生成数据的过程，我们可以用给定的方法对这些数据进行评估4。然后，我们一遍又一遍地重复这个过程。之后，我们查看所有重复操作后的结果。
为什么要这样做？通过使用我们控制的随机过程，我们可以选择真相。也就是说，我们不必担心或想知道哪些假设适用于我们的案例。我们可以让假设成立，也可以不成立，只需改变我们的数据生成过程。我们可以选择是否存在后门、碰撞，或者误差项是否相关，或者是否存在共同的支持。
然后，当我们看到得到的结果时，我们就能看到在那个特定的真实情况设定下，该方法会产生什么样的结果。由于我们 选定了 真实的因果效应，我们知道真实的因果效应是什么——我们的方法平均而言能得到真实结果吗？抽样误差有多大？与另一种方法相比，这个估计是更精确还是更不精确呢？如果我们确实得到了真实的因果效应，要是我们重新运行模拟，但让误差项相关，会发生什么情况呢？它还会得出真实的因果效应吗？
当应用一种新方法，或者在不确定某种已知方法是否可行的情境中应用该方法时，直接通过模拟来进行尝试是个很棒的主意。本章将为你提供开展模拟所需的工具。
15.1.1 剖析模拟
为了给出一个基本的例子，假设我们有一枚硬币，我们要抛掷\(100\)次，然后估计正面朝上的概率。 我们可能会使用什么估计方法？我们可能只是计算一下在我们的\(N=100\)的样本中正面朝上的比例。
模拟会是什么样子呢？
确定我们的随机过程/数据生成过程是什么。我们可以自主设定！假设我们选择一枚真正有\(75\%\)概率正面朝上的硬币。
使用随机数生成器来生成我们\(100\)次抛硬币的样本，其中每次抛硬币都有\(75\%\)的概率是正面朝上。
将我们的估计方法应用到这\(100\)次抛硬币的样本上。换句话说，计算该样本中正面朝上的比例。
将我们的估计值存储在某个地方。
多次重复步骤\(2-4\)，每次都抽取一个新的\(100\)次抛硬币的样本，将我们的估计方法应用到该样本上，并存储估计值（如果愿意，还可以存储估计过程中的其他一些信息）。
查看我们每次迭代得到的估计值的分布情况。
就是这样！对于这些估计值的分布，我们有很多事情可以做。首先，我们可以看看这个分布的均值。我们所设定的“真实情况”是硬币正面朝上的概率为\(75\%\)。所以，如果我们的估计值与\(75\%\)正面朝上这个结果相差甚远，那就表明我们的估计不太准确。
我们可以通过模拟来排除那些糟糕的估计量。例如，也许出于某种原因，我们不是用样本中正面朝上的比例来估计正面朝上的概率，而是用样本中正面朝上比例的平方来进行估计。你可以自己试试看——对每一次迭代都应用这个估计量，我们将无法得到正面朝上概率为\(75\%\)的估计结果。
我们也可以用这个方法来亲自验证优良估计量的价值。别盲目相信你那本陈旧乏味的统计学教材里所说的，即正面朝上的比例是对正面朝上概率的一个良好估计。自己动手试试！你可能早就知道教材是对的。但我发现，进行模拟也能让你切实体会到它们是对的，并且理解这个估计量为何有效，尤其是当你对其进行一番探究并尝试一些变体的时候5。
就我个人而言，当我想判断自己打算采用的方法是否合理时，模拟绝对是我首先会使用的工具。你有了一个新方法的点子，或者想对现有方法做些调整，又或者只是有了一种对特定变量进行转换或纳入的方式？试着对其进行模拟，看看你想做的事情是否”行得通”。这并不太难，而且肯定比不做验证就完成项目，然后在最后当别人读了你的成果并指出你的错误时，才发现自己的做法行不通要容易得多。
除了观察我们的估计量是否能准确得出真实值之外，我们还可以用模拟的方法来了解抽样变异性。当然，平均而言，我们得到正面朝上的比例是 \(75\%\)，但有多大比例的样本会得到正面朝上的比例为\(50\%\)呢？或者是\(90\%\)呢？又有多大比例的样本会发现正面朝上的比例与\(75\%\)在统计学上存在显著差异呢？
模拟是理解估计量抽样分布的重要工具。抽样分布（以及它为我们提供的我们非常渴望得到的标准误差）可能真的很难弄清楚，尤其是当分析包含多个步骤，或者使用多个估计量时。模拟有助于在无需撰写证明的情况下弄清楚抽样变异性的程度。
当然，要理解我们抛硬币估计量的标准误差，我们并不需要模拟。如果你找到合适的维基百科页面，答案会写得相当清楚。但要是我们先根据硬币的可观测特征，使用匹配法将这枚硬币与其他硬币关联起来，再应用一个平滑函数来处理过去观察到的异常硬币结果，然后剔除那些我们忘记记录硬币抛掷结果时产生的缺失观测值，这时标准误差是多少呢？呃……我迷糊了。我真的不想去为那个标准误差写一个证明。而且就算有人已经为我写好了，我也不确定自己能否找到并理解它6。即使在复杂的情况下，模拟也能帮助我们了解抽样变异性。
的确，即使是那些知道如何撰写证明的专业人士也会使用模拟来处理抽样变异问题——可参考本章中关于自助法标准误差和功效分析的部分。
最后，模拟对于理解一种估计方法在其假设被违背时的稳健性如何，以及其估计值对真实情况变化的敏感性如何至关重要。当然，我们知道，如果一个真实模型中针对某个变量存在二阶多项式关系，但我们在进行估计时采用的是不包含该多项式的线性回归，那么得出的结果将会是错误的。但这个错误会有多大呢？是小错误还是大错误？当该变量是一个极强的预测变量或者是一个较弱的预测变量时，错误会更大吗？当斜率为正或为负时，错误会更大吗，还是说这根本无关紧要呢？
在进行模拟时，对同一个模拟运行多种变体是非常常见的做法。对假设做一些小的调整——在一个模拟中设置一个更强的预测变量，在另一个模拟中设置一个较弱的预测变量。在误差项中加入相关性。大胆尝试各种情况！你现在有了一个”玩具”，而你现在的任务就是把它弄坏，这样你就能确切地知道它在何时以及如何出问题了。
15.1.2 我们为何进行模拟
综上所述，模拟是一个利用我们所控制的随机过程，将我们的估计量反复应用于不同样本的过程。这使我们能够完成一些非常重要的事情：
我们可以检验当把我们的估计方法应用于我们所选择的随机过程/数据生成过程时，是否能得出真实的效应。
我们可以了解我们的估计值的抽样分布情况。
我们可以观察一个估计量对某些假设的敏感程度。如果我们让这些假设不再成立，这个估计量是否就不再适用了呢？
现在，我们已经知道为什么模拟是有益的了7。现在我们需要搞清楚怎么去实现它
15.2 哎，我们必须编码了
没有别的办法。现代模拟就是要使用代码。而且你不只是要正确使用预先封装好的命令——实际上你还得自己编写代码8。
值得庆幸的是，模拟不仅相对简单，而且还是编写更复杂程序的一个不错的切入点。
一个典型的统计模拟只有几个部分：
生成随机数据的过程
一次估计
多次重复步骤\(1-2\)的过程
存储每次迭代估计结果的方法
查看各次迭代估计值分布情况的方法
让我们逐步讲解对这些内容进行编码的各个流程。
我们如何对数据生成过程进行编码呢？这是一个关键问题。通常，当你进行模拟时，其核心目的往往是观察某个估计量在特定的一组数据生成假设下的表现。所以我们需要知道如何生成满足那些假设的数据。
第一个问题是如何生成一般的随机数据9。通常，有生成根据各种分布的随机数据的函数（如\(第3章\)所述）。但我们通常可以用一些函数来解决：
- 随机均匀数据，用于生成必须位于特定范围内的变量。\(R\)中的 \(runif()\) 、\(Stata\)中的 \(runiform()\) 和 \(Python\)中的\(numpy.random.uniform()\)10。
- 正态随机数据，用于生成连续变量，并在许多情况下模仿现实世界大部分变量的分布。\(R\) 中的\(rnorm()\)、\(Stata\)中的\(normal()\) 和\(Python\)中的\(numpy.random.normal()\)。
- 分类随机数据，用于分类变量或二元变量。\(R\)中的\(sample()\)、\(Stata\)中的 \(runiformint()\)或\(Python\)中的\(numpy.random.randint()\)。后两者仅在您想要每个类别机会均等的情况下才有效。如果您想要类别的值为\(1\)且出现概率为\(p\)的二元数据，则\(Stata\)中有\(rbinomial(1,p)\)和\(Python\)中有 \(numpy.random.binomial(1,p)\)11。
你可以使用这些简单的工具来生成你喜欢的数据生成过程。例如，看看后面几段代码片段，我在其中生成了一个数据集，其中包含一个正态分布的\(\epsilon\)误差变量，我称之为\(epsilon\)12，一个概率为\(0.2\)时取值为\(1\)，否则为\(0\)的二元\(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 errors# If we want the random data to be the same every time, set a seed
set.seed(1000)
# tibble() creates a tibble; it's like a data.frame
# We must pick a number of observations, here 200
d <- tibble(eps = rnorm(200), # by default mean 0, sd 1
Y = runif(200), # by default from 0 to 1
X = sample(0:1, # sample from values 0 and 1
200, # get 200 observations
replace = TRUE, # sample with replacement
prob = c(.8, .2))) # Be 1 20% of the time
head(d)## # A tibble: 6 × 3
## eps Y X
## <dbl> <dbl> <int>
## 1 -0.446 0.821 0
## 2 -1.21 0.209 0
## 3 0.0411 0.767 0
## 4 0.639 0.671 0
## 5 -0.787 0.0463 0
## 6 -0.385 0.754 0使用这三个随机数生成器13，然后我们可以模拟因果图。毕竟，这是一本因果推论教科书。我们想知道我们的估计方法在面对不同的因果关系时会表现如何！
15.2.1 创建模拟数据
我们如何创建遵循因果图的数据？简单！如果我们使用变量\(A\)来创建变量\(B\)，那么\(A \rightarrow B\)。就这么简单。以下面的代码片段为例，它们修改了我们上面所做的，以添加一个后门混淆变量\(W\)（并将\(Y\)变成正态分布而不是均匀分布）。换句话说，这些代码片段是\(图15.1\)的表示。
如\(图15.1\)所示，我们知道\(W\)会导致\(X\)和\(Y\)的产生。因此，我们设定\(W\)会使二元变量\(X\)更有可能取值为\(1\)，但会降低连续变量\(Y\)的值。我们还知道，\(\epsilon\)是产生\(X\)的一个因素，因为它会导致\(X\)的出现；并且\(X\)是产生\(Y\)的一个因素，\(X\)对\(Y\)的真实影响是，\(X\)每增加一个单位会导致\(Y\)增加\(3\)个单位。最后，虽然图中未标明，但在\(Y\)的产生过程中，我们还有另一个随机因素\(\nu\)。如果真正指向\(Y\)的箭头只有来自\(X\)和\(W\)，那么一旦我们控制了这两个变量，\(Y\)就会完全由它们决定，不会有其他变化。你可以想象在\(图15.1\)中添加一个从ν指向\(Y\)的箭头。
但可以推测，还有其他的因素影响着\(Y\)，那就是\(\nu\)，如果我们用\(X\)和\(W\)对\(Y\)进行回归，它也是误差项。\(\nu\)在代码中的哪里呢？它就是创建\(Y\)时添加的随机正态数据。你还会注意到，\(\epsilon\)也从代码中消失了。但实际上它仍然存在——它是均匀随机数据！我们通过检查一些均匀随机数据\(\epsilon\)是否低于\(0.2\)加上\(W\)来创建\(X\)——\(W\)值越高，就越有可能低于\(0.2+W\)，因此，\(X\)就越有可能等于\(1\)。“分割”均匀随机变量的不同值有时可能比使用二项式随机数据更容易，特别是当我们希望不同观测值的概率不同时。
library(tidyverse)
# If we want the random data to be the same every time, set a seed
set.seed(1000)
# tibble() creates a tibble; it's like a data.frame. We must pick a number of observations, here 200
d <- tibble(W = runif(200, 0, .1)) %>% # only go from 0 to .1
mutate(X = runif(200) < .2 + W) %>%
mutate(Y = 3*X + W + rnorm(200)) # True effect of X on Y is 3
head(d)## # A tibble: 6 × 3
## W X Y
## <dbl> <lgl> <dbl>
## 1 0.0328 TRUE 3.95
## 2 0.0759 FALSE 0.806
## 3 0.0114 FALSE -1.67
## 4 0.0691 FALSE -0.164
## 5 0.0516 FALSE -0.608
## 6 0.00677 FALSE -0.0695我们可以更进一步。我们不仅想要随机创建这些数据，而且在反复运行模拟时，我们会想要这样做很多次。每当我们有想要运行很多次的代码时，把它放在一个独立的函数里通常是个好主意，这样我们就可以反复调用它。我们还会给这个函数一个参数，让我们能够选择不同的样本大小。
library(tidyverse)
# Make sure the seed goes OUTSIDE the function. It makes the random data the same every time,
# but we want DIFFERENT results each time, we run it (but the same set of different results, thus the seed)
set.seed(1000)
# Make a function with the function() function. The "N = 200" argument gives it an argument N that we'll use for sample size. The "=200" sets
# the default sample size to 200
create_data <- function(N = 200) {
d <- tibble(W = runif(N, 0, .1)) %>%
mutate(X = runif(N) < .2 + W) %>%
mutate(Y = 3*X + W + rnorm(N))
# Use return() to send our created data back
return(d)
}
# Run our function!
create_data(500)## # A tibble: 500 × 3
## W X Y
## <dbl> <lgl> <dbl>
## 1 0.0328 FALSE 0.724
## 2 0.0759 FALSE 0.465
## 3 0.0114 FALSE -0.0832
## 4 0.0691 FALSE -1.88
## 5 0.0516 FALSE -0.614
## 6 0.00677 FALSE 0.0189
## 7 0.0739 FALSE 0.495
## 8 0.0584 TRUE 2.34
## 9 0.0216 FALSE 0.450
## 10 0.0256 FALSE 0.00959
## # ℹ 490 more rows如果我们想要更复杂的东西呢？那我们就得更仔细地考虑如何创建我们的数据了！
在创建更复杂的随机数据时，一个常见的需求是模拟面板数据结构。我们如何为每个个体创建包含多条观测值的数据呢？我们要创建一个个体标识符，为该标识符赋予一些个体特征，然后将这些特征应用于整个数据的时间范围。在\(R\)和\(Python\)中，通过合并个体数据集和个体-时间数据集来实现这一点是最简便的；在\(Stata\)中，则可以使用\(egen\)命令来完成。
library(tidyverse)
set.seed(1000)
# N for number of individuals, T for time periods
create_panel_data <- function(N = 200, T = 10) {
# Create individual IDs with the crossing() function, which makes every combination of two vectors
# (if you want some to be incomplete, drop them later)
panel_data <- crossing(ID = 1:N, t = 1:T) %>%
# And individual/time-varying data (n() means "the number of rows in this data"):
mutate(W1 = runif(n(), 0, .1))
# Now an individual-specific characteristic
indiv_data <- tibble(ID = 1:N, W2 = rnorm(N))
# Join them
panel_data <- panel_data %>%
full_join(indiv_data, by = 'ID') %>%
# Create X, caused by W1 and W2
mutate(X = 2*W1 + 1.5*W2 + rnorm(n())) %>%
# And create Y. The true effect of X on Y is 3. But W1 and W2 have causal effects too
mutate(Y = 3*X + W1 - 2*W2 + rnorm(n()))
return(panel_data)
}
create_panel_data(100, 5)## # A tibble: 500 × 6
## ID t W1 W2 X Y
## <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 0.0328 0.594 0.550 1.61
## 2 1 2 0.0759 0.594 -0.397 -1.40
## 3 1 3 0.0114 0.594 -0.0922 -3.33
## 4 1 4 0.0691 0.594 1.73 2.59
## 5 1 5 0.0516 0.594 0.536 0.527
## 6 2 1 0.00677 -0.115 -0.0629 1.01
## 7 2 2 0.0739 -0.115 1.55 3.74
## 8 2 3 0.0584 -0.115 -0.799 -2.09
## 9 2 4 0.0216 -0.115 -0.868 -1.59
## 10 2 5 0.0256 -0.115 0.414 1.38
## # ℹ 490 more rows在生成随机数据时，另一个紧迫的问题是创建相关的误差项。如果我们担心误差项不独立所带来的影响，那么模拟是一种很好的尝试方法，前提是我们能够使误差项具有我们所期望的相关性。需要考虑的情况有无数种，但我将考虑两种：异方差性和聚类。
异方差性相对容易创建。每当你创建用于表示误差项的变量时，我们只需让该数据的方差根据观测值的特征而变化。
library(tidyverse)
set.seed(1000)
create_het_data <- function(N = 200) {
d <- tibble(X = runif(N)) %>%
# Make the standard deviation of the error term correlated with X. This is heteroskedasticity!
mutate(Y = 3*X + rnorm(N, sd = X*5))
return(d)
}
create_het_data(500)## # A tibble: 500 × 2
## X Y
## <dbl> <dbl>
## 1 0.328 1.96
## 2 0.759 1.84
## 3 0.114 0.672
## 4 0.691 5.23
## 5 0.516 0.619
## 6 0.0677 0.408
## 7 0.739 -0.0803
## 8 0.584 2.43
## 9 0.216 -1.30
## 10 0.256 2.12
## # ℹ 490 more rows如果我们绘制\(create\_het\_data\)的结果，就会得到\(图15.2\)。
那里的异方差性误差非常明显！在图的右侧，围绕回归线的波动比左侧大得多。这是异方差性的典型”扇形”形状。
library(ggplot2)
library(ggpubr)
d <- create_het_data(500)
ggplot(d, aes(x = X, y = Y)) +
geom_point(alpha = .4) +
geom_smooth(formula = 'y~x', method = 'lm', se = FALSE, color = 'black', linewidth = 1.5) +
theme_pubr() +
theme(text = element_text(family = 'sans', size =10)) 
ggsave('Simulation/heteroskedastic_data.pdf', width = 6, height =5,device=cairo_pdf)聚类分析要稍微难一些。首先，我们需要有可供聚类的对象。这通常是某种组标识变量，就像我们在本节前面讨论如何模拟面板结构时提到的那些变量一样。然后，我们可以通过创建或随机生成一个可由组内成员共享的单一“组效应”，再加上个体噪声，来创建聚类标准误。
换句话说，先生成一个个体层面的误差项，再生成一个个体/时间层面的误差项，并以某种方式将它们相加，从而得到聚类标准误。幸运的是，我们可以复用许多在创建面板数据时学到的技巧。
那我们的估计量呢？没有估计量，因果推断模拟就不完整。代码方面，我们懂的。就是——呃，在我们刚随机生成的数据上跑一下我们感兴趣的估计量。
假设我们在教科书中读到，聚类误差\(（clustered \space errors）\)只会影响线性回归的标准误\(（standrad \space errors）\)，而不会使估计量产生偏差。我们想亲自验证这一点。于是，我们取一些已经内置了聚类误差的数据，并对其估计线性回归。
此时的关键问题通常在于如何从回归结果中提取你想要的估计值（通常是一个系数），并将其返回并存储。如果你要获取的是一个线性回归系数：
在 \(R\) 中，可以使用 \(coef(m)[`varname']\) 来获取模型\(M\)中变量 \(vaename\) 的系数。
在 \(Stata\) 中，可以在运行模型后使用 \(\_b[varname]\) 来提取该系数。
在 \(Python\) 中，可以通过 \(m.params[`varname']\) 获取模型\(M\)中变量 \(varname\) 的系数14。
既然这又是我们之后会反复调用的功能，我们也将它单独封装成一个函数。而这个函数本身会调用我们刚刚创建的数据生成函数！我们将引用之前章节中编写的 \(create\_clus\_data\) 函数。
library(tidyverse)
set.seed(1000)
# N 代表个体数量，T 代表时间段数量
create_clus_data <- function(N = 200, T = 10) {
# 我们打算创建在个体标识（ID）层面聚集的误差。这样我们就可以按照制作面板数据的步骤来进行操作
panel_data <- crossing(ID = 1:N, t = 1:T) %>%
# 个体/随时间变化的数据
mutate(W = runif(n(), 0, .1))
# 1个特定个体的误差聚类
indiv_data <- tibble(ID = 1:N, C = rnorm(N))
# Join them
panel_data <- panel_data %>% full_join(indiv_data , by = 'ID') %>%
# 创建由 W1 和 W2 导致的 X
mutate(X = 2*W + rnorm(n())) %>%
# 误差项有两个组成部分：个体聚类（或个体群组）C，以及随个体和时间变化的元素
mutate(Y = 3*X + (C + rnorm(n())))
return(panel_data)
}
create_clus_data(100,5)## # A tibble: 500 × 6
## ID t W C X Y
## <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 0.0328 0.594 -0.341 0.681
## 2 1 2 0.0759 0.594 -1.29 -2.37
## 3 1 3 0.0114 0.594 -0.984 -4.23
## 4 1 4 0.0691 0.594 0.837 1.63
## 5 1 5 0.0516 0.594 -0.355 -0.416
## 6 2 1 0.00677 -0.115 0.109 1.18
## 7 2 2 0.0739 -0.115 1.72 3.84
## 8 2 3 0.0584 -0.115 -0.627 -1.97
## 9 2 4 0.0216 -0.115 -0.696 -1.44
## 10 2 5 0.0256 -0.115 0.586 1.52
## # ℹ 490 more rows现在我们必须进行多次迭代。我们不能只运行一次模拟——我们得到的任何结果都可能是偶然的。此外，我们通常想了解估计量所给出结果的 分布，要么是为了检查其离散程度，要么是为了看看其均值是否接近真实值。没有多个样本就无法得到抽样分布。你应该运行多少次模拟呢？并没有一个确切的数字，但应该进行很多次！你在测试代码以确保其能正常运行时，可以只进行几次迭代，但当你想要得到结果时，就应该增加迭代次数。根据经验法则，你希望结果有更多细节，或者估计中涉及的变量越多，你就应该增加模拟的次数——想查看均值的自助法分布？大约一千次应该就足够了。想了解相关系数或线性回归系数呢？几千次。想精确估计抽样分布的多个百分位数？也许要进行 数万 次。
迭代是计算机非常擅长的事情之一。在我们的三个代码示例中，我们都将以略有不同的方式使用”\(for \space 循环\)“。\(for \space 循环\)是计算机编程中的一种表达方式，意思是”针对这些值中的每一个，把同一件事情做很多次”。
library(tidyverse)
set.seed(1000)
# A function for estimation
est_model <- function(N, T) {
# Get our data. This uses create_clus_data from earlier
d <- create_clus_data(N, T)
# Run a model that should be unbiased if clustered errors themselves don't bias us!
m <- lm(Y ~ X + W, data = d)
# Get the efficient on X, which SHOULD be true value 3 on average
x_coef <- coef(m)['X']
return(x_coef)
}
# Run our model
est_model(200, 5)## X
## 3.050733例如，（通用语言/伪代码）中的 \(for \space 循环\)：\(for \space i \space in \space [1,4,5,10]: print \space i\) 会多次执行同一件事\(（print \space i）\)，针对我们列出的每个值\(（1、4、5、10）\)各执行一次。
首先，当 \(i\) 被设为 \(1\) 时执行该操作，从而打印出数字 \(1\)。接着，当 \(i\) 被设为 \(4\) 时执行该操作，打印出数字 \(4\)。然后打印 \(5\)，最后打印 \(10\)。 在模拟的情况下，我们的\(for \space 循环\)会是这样：\(for \space i \space in \space [1 \space to \space 100]: do \space my \space simulation\)。
我们甚至根本没有使用迭代数字\(i\)15，而是仅仅重复我们的模拟，在这种情况下是\(100\)次。
这种循环的实现方式在每种语言中略有不同。\(Python\)拥有最纯粹、“最常规”的循环，它看起来会和通用语言中的循环很相似。\(R\)理论上也可以实现类似的循环，但我建议使用 \(purrr\)包及其 \(map\) 函数，这将为我们后续省去一个步骤。对于\(Stata\)而言，循环本身会比较常规，但由于我们一直在创建新数据，再加上\(Stata\)一次只能使用一个数据集的限制，这就要求我们反复清空数据（此外，这还会在后续增加一个步骤）。
所有这些代码示例都会引用我们在前面代码块中定义的\(create\_clus\_data\)和\(est\_model\)函数。
现在需要一种存储结果的方法。我们已经有一个用于运行估计的函数，而这个函数本身会调用创建模拟数据的函数。我们需要一种方法来存储这些估计结果。
在 \(R\) 语言和 \(Python\) 语言中，我们采用的迭代方法已经足够了。\(R\) 语言里的 \(map()\) 系列函数，以及 \(Python\) 语言里的标准\(for\)循环结构，都能让我们得到涵盖所有迭代结果的向量。所以对于这两种语言而言，我们已经大功告成了。最后一个代码块里已经有你所需的指令了。
鉴于\(Stata\)一次只倾向于处理一个数据集，情况就有点棘手了。你如何能同时存储随机生成的数据和估计结果数据呢？我们可以使用 \(preserve/restore\)技巧——将结果存储在一个数据集中，使用\(preserve\)保存该数据集，这样\(Stata\)就会记住它，然后清空所有内容并运行我们的估计，将估计值保存为一个 \(local\)（局部宏），再使用 \(restore\) 恢复到我们保存结果的数据集，最后将那个 \(local\) 赋值到数据中以便记住该估计值16。
现在，我要说的是，有一种方法可以解决这种折腾，那就是使用\(Stata\)中的 \(simulate\) 命令。 使用\(simulate\)，您可以指定一个用于创建数据和估计模型的函数，然后将所有结果拉回来。 \(help \space simulate\) 应该包含您需要知道的一切。 我发现它不如显式编码灵活，而且工作量没有减少多少，所以我很少在本书中使用它17。 但您可能想看看它——您可能会发现自己更喜欢它。
library(tidyverse)
library(purrr)
set.seed(1000)
# Estimate our model 1000 times (from 1 to 1000)
estimates <- 1:1000 %>%
# Run the est_model function each time
purrr::map_dbl(function(x) est_model(N = 200, T = 5))
# There are many map functions in purrr. Since est_model outputs anumber (a "double"),
# we can use map_dbl and get a vector of estimates
mean(estimates)## [1] 3.002177sd(estimates)## [1] 0.04318719# A function for estimation
est_model_se <- function(N, T) {
# Get our data
# This uses the create_clus_data function from the previous section
d <- create_clus_data(N, T)
# Run a model that should be unbiased \# if clustered errors themselves don't bias us!
m <- lm(Y ~ X + W, data = d)
# Get the coefficient on X, which SHOULD be the true value 3 on average
x_coef <- coef(summary(m))['X', 'Std. Error']
return(x_coef)
}
set.seed(1000)
estimates <- 1:300 %>%
# Run the est_model function each time
map_dbl(function(x) est_model_se(N = 200, T= 5))
mean(estimates)## [1] 0.04463276最后，我们进行评估。现在我们可以检验抽样估计量的任何我们想了解的特征！手头有了所有随机生成样本的一系列结果后，比如说，我们可以查看抽样变异的均值是否就是我们所设定的真实效应！
让我们在这里检查一下。我们从数据生成函数中得知 \(X\) 对 \(Y\) 的真实效应是 \(3\)。在我们的模拟分析中，我们是否得到了一个接近 \(3\) 的平均值？在 R 中使用简单的\(mean()\)，在 \(Stata\) 中使用 \(summarize\) 或 \(mean\)，或者在 \(Python\) 中使用 \(numpy.mean()\)，得到的结果平均值为 \(3.00047\)。非常接近！看来我们的教科书是对的：单独的聚类误差不会偏离一个原本无偏的估计。
那么标准误差呢？我们所有随机样本的估计值分布应该近似于抽样分布。所以，我们结果的标准差就是估计值的标准误差。在 \(R\) 语言中使用 \(sd()\) 函数、在 \(Stata\) 软件中使用 \(summarize\) 命令，或者在 \(Python\) 中使用 \(numpy.std()\) 函数，得出的标准差为 \(0.045\)。
因此，一个合理的线性回归估计量应该给出\(0.045\)的标准误差。是这样吗？ 我们可以使用模拟来回答这个问题。 我们可以使用我们已经有的模拟代码，但将存储系数估计值的部分替换为存储标准误差。 这可以很容易地通过 \(R\) 语言中 \(broom\) 包的 \(tidy()\) 函数18、\(Stata\) 软件中的 \(\_se[X]\) 或 \(Python\) 中的 \(m.bse[`x']\) 来完成，其中 \(x\) 是你想要获得标准误差的变量名称。 当我们模拟几百次后，报告的标准误差的平均值为 \(0.044\)。 稍微小了一点。 如果我们改为估计异方差稳健标准误差，我们会得到 \(0.044\)。 没有更好！ 如果我们以最合适的方式处理标准误差，考虑到我们知道我们创建的聚类误差，并使用在组级别聚类的标准误差，我们会得到 \(0.045\)。 看来这效果很好！
library(fixest)
# A function for estimation
est_model_het <- function(N, T) {
# Get our data \# This uses the create_clus_data function from the previous section
d <- create_clus_data(N, T)
# Run a model that should be unbiased if clustered errors themselves don't bias us!
m <- feols(Y ~ X + W, data = d, se = 'hetero')
# Get the coefficient on X, which SHOULD be the true value 3 on average
x_coef <- broom::tidy(m)$std.error[2]
return(x_coef)
}
set.seed(1000)
estimates <- 1:300 %>%
# Run the est_model function each time
purrr::map_dbl(function(x) est_model_het(N = 200, T = 5))
mean(estimates)## [1] 0.04476346est_model_clus <- function(N, T) {
# Get our data
# This uses the create_clus_data function from the previous section
d <- create_clus_data(N, T)
# Run a model that should be unbiased if clustered errors themselves don't bias us!
m <- feols(Y ~ X + W, data = d, cluster = ~ID)
# Get the coefficient on X, which SHOULD be the true value 3 on average
x_coef <- broom::tidy(m)$std.error[2]
return(x_coef)
}
set.seed(1000)
estimates <- 1:300 %>%
# Run the est_model function each time
map_dbl(function(x) est_model_clus(N = 200,T = 5))
mean(estimates)## [1] 0.04464148library(fixest)
library(tidyverse)
library(progressr)
# 改进后的函数
est_model_clus <- function(N, T) {
d <- create_clus_data(N, T)
m <- feols(Y ~ X + W, data = d, cluster = ~ID)
return(tibble(coef = coef(m)["X"], se = se(m)["X"]))
}
# 带进度条的模拟
set.seed(1000)
with_progress({
results <- map_df(1:300, ~ est_model_clus(200, 5))
})
# 计算关键指标
metrics <- results %>%
summarise(
mean_coef = mean(coef),
mean_se = mean(se),
empirical_sd = sd(coef),
coverage = mean(coef - 1.96*se < 3 & coef + 1.96*se > 3),
se_ratio = mean_se / empirical_sd # 标准误估计偏差比例
)
print(metrics)## # A tibble: 1 × 5
## mean_coef mean_se empirical_sd coverage se_ratio
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3.00 0.0446 0.0452 0.957 0.98815.3 The Old Man and the C
正如我之前提到的19，在因果推断的背景下，模拟可以做的最强大的事情之一是让你知道一个给定的估计量是否有效，是否能平均提供真实的因果效应。将真实的因果效应融入你的数据生成函数中，然后将你感兴趣的估计量放入估计函数中，最后让它尽情发挥！如果估计值的均值接近真实值，并且分布不是太宽，那就没问题了。
这是一种绝佳的方式，可用来检验你所接触到的 任何 新方法或统计学论断。我们已经验证了：在存在聚类误差的情况下，普通最小二乘法平均而言是否仍然有效。但我们还能做更多的验证。如果对处理变量的影响具有相反的符号，那么一个后门路径能否抵消另一个后门路径的影响呢？当其中一个控制变量与结果变量存在非线性关系时，匹配法真的比回归法效果会更好吗？控制对撞因子 真的 会增加偏差吗？创建一个数据集，其中存在一条被对撞因子阻断的后门路径，然后看看在你的估计函数中控制或不控制这个对撞因子到底会出现什么情况20。
这使得模拟成为用于以下三件事的绝佳工具：\(（1）\)试用新的估计量；\(（2）\)比较不同的估计量；\(（3）\)了解为了使某个估计量失效我们需要破坏哪些条件。
我们无需撰写计量经济学证明就能完成所有这些事情。当然，我们得出的结果不会像证明那样具有广泛的普适性，但这是很好的第一步。此外，这让主动权掌握在我们自己手中！
我们可以想象有一位老人被困在某个荒岛上，不知怎的，他配备了一台笔记本电脑和他最喜欢的统计软件，但却无法查阅教科书，也找不到专家咨询。也许他叫阿古斯\(（Agus）\)。阿古斯有自己的统计学疑问——哪些方法可行？哪些不可行？也许他有一个关于新估计量的绝妙想法。他可以自己去尝试验证。我建议你也这么做。通过模拟来尝试各种方法，不仅是一种自己寻找答案的好方式，也是一种了解数据运行机制的好途径。
15.3.1 尝试新的估计量
教授们往往会向你展示那些效果良好的估计量。但这可能会给人一种错觉，以为这些估计量是凭空出现的，而且个个都能完美发挥作用。事实并非如此！总得有人先想出这些估计量，然后进行试验，并证明它们是有效的。
你完全有理由自己想出一个估计量并进行试验。当然，它很可能不会奏效。但我打赌你会有所收获。
咱们来试试自己搞一个估计器，或者试试\(Agus\)的估计器。\(Agus\)有个想法：你知道回归的标准误差会随着 \(X\) 的方差增大而减小吗？如果我们只选择 \(X\) 的极端值呢？那样 \(X\) 的方差就会变得非常大。也许那样我们最终会得到更小的标准误差。
所以这就是\(Agus\)的估计器：不是像往常一样把 \(Y\) 对 \(X\) 做回归，而是先筛选数据，只留下 \(X\) 的极端值（大概是顶部和底部各\(20\%\) ），然后用这些数据把 \(Y\) 对 \(X\) 做回归21。我们觉得这样也许能降低标准误差，不过还得确保这不会引入偏差。
阿古斯有许多这样的点子。他把它们写在棕榈叶上，然后用他那台老旧沙尘斑斑的电脑编程验证是否可行。对于行不通的那些，他会将棕榈叶轻放海面，轻轻一推任其漂向地平线，与之作别。这个估计量也会加入海中同僚的行列吗？
让我们亲自动手编码验证吧。你真的能自己编写估计算法吗？当然可以！只需编写让计算机执行你估计算法的代码。若不知如何编写这些代码，就上网搜索直至掌握为止
这段代码还有一处调整。我们不仅关注效应是否无偏，还希望了解标准误是否会降低。因此，可能需要同时收集系数及其标准误——需要返回两个值。该怎么做呢？看代码。
library(tidyverse)
library(purrr)
library(broom)
set.seed(1000)
# Data creation function. Let's make the function more flexible so we can choose our own true effect!
create_data <- function(N, true) {
d <- tibble(X = rnorm(N)) %>%
mutate(Y = true*X + rnorm(n()))
return(d)
}
# Estimation function. keep is the portion of data in each tail to keep. So .2 would keep the bottom and top 20% of X
est_model <- function(N, keep, true) {
d <- create_data(N, true)
# Agus' estimator!
m <- lm(Y~X, data = d %>% dplyr::filter(X <= quantile(X, keep) | X >= quantile(X, (1-keep))))
# Return coefficient and standard error as two elements of a list
ret <- tidy(m)
return(list('coef'= ret$estimate[2], 'se' = ret$std.error[2]))
}
# Run 1000 simulations. use map_df to stack all the results together in a data frame
results <- 1:1000 %>% map_df(function(x) est_model(N = 1000, keep = .2, true = 2))
mean(results$coef) ## [1] 1.99882sd(results$coef)## [1] 0.03373593mean(results$se)## [1] 0.03401751如果我们运行这个（程序或操作），我们会看到平均系数为 \(2.00\)（而真实值是 \(2\)），系数的标准差为 \(0.0332\)，系数的平均标准误差为 \(0.0339\)。看，我就说模拟系数的标准差与系数的标准误差是相符的吧！
所以平均值是 \(2.00\)，而真实值也是 \(2\)。没有偏差！相当不错，对吧？阿古斯应该庆祝一下了。当然，我们还没做完。我们不只是想要无偏差，我们还希望比常规估计方法有更低的标准误差。我们要怎么做才能实现这一点呢？
15.3.2 比较估计量
模拟能让我们做的一件事是比较不同方法的性能 。你可以拿完全相同的随机生成的数据，对其应用两种不同的方法，然后看看哪种方法表现更好。这个过程通常被称为“赛马比拼”。
当考虑用多种不同方法做同一件事时，“赛马比拼”法真的很有用。例如，在\(第14章\)关于匹配的内容中，我提到了在匹配过程中可以做出的一系列不同选择。以一种方式做出这些选择能够减少偏差，但以另一种方式做出则可以提高精度。放宽卡尺范围会纳入更多但质量更差的匹配项，这会增加偏差，但同时也会使抽样分布变窄，这样你就不太可能出现极大的误差。这是一个很难做出的抉择——是降低平均偏差，还是使抽样分布更窄呢？但如果结果恰好是你所实现的偏差减少量根本就不高，而在精度上的提升却非常大，那么做这个抉择可能就会容易一些！你怎么才能知道这种情况呢？可以尝试进行一场”赛马比拼”模拟。
当然，你所权衡的偏差量和精度会因你的数据情况而异。所以你可能想使用你自己的实际数据来进行模拟！在这种情况下，你并不知道”真实”值，但你可以检查平衡性来大致了解偏差情况。关于如何使用实际数据进行模拟，可参考下面的自助法部分内容22。
在阿古斯的估计量这个案例中，我们想知道它是否会比常规回归产生更低的标准误差。这就是核心想法，对吧？那么让我们来检验一下。这些代码示例都使用了上一组代码块中的
create_data 函数。
library(tidyverse)
library(purrr)
library(broom)
library(vtable)## Loading required package: kableExtra##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowsset.seed(1000)
# Estimation function. keep is the portion of data in each tail to keep. So .2 would keep the bottom and top 20% of X
est_model <- function(N, keep, true) {
d <- create_data(N, true)
# Regular estimator!
m1 <- lm(Y~X, data = d)
# Agus' estimator!
m2 <- lm(Y~X, data = d %>% filter(X <= quantile(X, keep) | X >= quantile(X, (1-keep))))
# Return coefficients as a list
return(list('coef_reg' = coef(m1)[2], 'coef_agus' = coef(m2)[2]))
}
# Run 1000 simulations. Use map_df to stack all the results together in a data frame
results <- 1:1000 %>% map_df(function(x) est_model(N = 1000, keep = .2, true = 2))
sumtable(results)| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| coef_reg | 1000 | 2 | 0.031 | 1.9 | 2 | 2 | 2.1 |
| coef_agus | 1000 | 2 | 0.034 | 1.9 | 2 | 2 | 2.1 |
尽管阿古斯和常规估计器得到的系数平均值都是\(2.00\)，但常规估计器的标准误差实际上更低。系数分布在所有模拟样本中的标准差，阿古斯的是 \(0.034\)，而常规估计器的是\(0.031\)23。因此，常规估计器表现更好。
15.3.3 打破常规
统计分析中充斥着各种各样的假设。说真的，这挺烦人的。不过呢，关键在于——为了确切了解我们所使用的估计量的性质，所有这些假设都必须成立。但也许有一些假设即便不成立，我们仍然能对那些估计量的性质有个相当不错的认识。
那么有多少假设不成立时我们仍然能勉强接受结果呢？模拟研究可以帮助我们弄清楚这一点。
例如，如果我们希望线性回归中自变量\(X\)的系数能无偏地估计\(X\)的因果效应，我们需要假设\(X\)与误差项不相关，也就是不存在开放的后门路径。
在社会科学领域，这几乎不可能是真的 。总会有一些我们无法测量或控制的小因素，而这些因素不巧地既与\(X\)相关，又与结果变量相关。那么我们是不是就永远都不用考虑做回归分析了呢？
没门儿！相反，我们应该问”这可能会是多大的问题呢？“要是这只是个小问题，嘿，即便没完全正确，但接近正确总好过差得老远！
通过创建一个非常灵活的数据生成函数，模拟可以用来探究这类问题。这个函数要足够灵活，以便尝试不同类型的假设违背情况或不同程度的假设违背。然后，使用不同的设置运行你的模拟，看看结果会糟糕到什么程度。如果即使数据违背了某个假设，结果仍然与真实情况相当接近，那么也许情况并没有那么糟糕24！
也许你可以将数据生成函数设计成带有一个名为 \(error\_dist\) 的参数。对其进行编程，使得当你将 \(error\_dist\) 设置为 “\(normal\)” 时，你会得到服从正态分布的误差项；设置为 “\(lognormal\)” 时，你会得到服从对数正态分布的误差项，以此类推，适用于各种分布。然后，你就可以模拟具有各种不同误差分布的线性回归，看看估计结果是否仍然理想。
在存在未解决的后门问题的应用中，我们如何运用模拟方法呢？我们可以对数据生成函数进行编程，为这个未解决的后门问题的影响强度设置参数。如果在我们的因果图中有 “\(X \leftarrow W \rightarrow Y\)” 这样的结构，那么 W 对 X 的影响是什么，\(W\) 对 \(Y\) 的影响又是什么呢？这两条箭头所代表的影响需要达到多强，才会使我们面临不可接受的高偏差水平呢？
一旦我们弄清楚这个问题需要严重到何种程度才会使我们面临不可接受的高偏差水平25，你就可以问问自己，在你所面临的情况中这个问题可能有多严重——如果问题的严重程度低于那个临界值，你或许还能接受（尽管在你的书面报告中，你肯定还是要提及这个问题，并且说明你认为情况尚可接受的原因）26。
让我们把这个实现为代码吧！我们要做的只是扩展我们的\(create\_data\)函数，使其包含一个变量\(W\)，该变量对我们的处理变量\(X\)和结果变量\(Y\)都有影响。然后，我们会设置 \(W \rightarrow X\) 和 \(W \rightarrow Y\) 关系的强度。实际上，我们还可以进行很多其他的调整 —— 比如 \(X \rightarrow Y\) 的强度、\(X\)和\(Y\)的标准差、添加其他控制变量、引入非线性关系、确定误差项的分布等等。但目前，让我们只使用这两个设置，把事情简化一些。
此外，我们还要把迭代/循环步骤也封装到一个函数里，因为我们后续要在不同的设置下多次进行这个迭代操作。
library(tidyverse)
library(purrr)
set.seed(1000)
# Have settings for strength of W -> X and for W -> Y
# These are relative to the standard deviation of the random components of X and Y, which are 1 each (rnorm() defaults to a standard deviation of 1)
create_data <- function(N, effectWX, effectWY) {
d <- tibble(W = rnorm(N)) %>%
mutate(X = effectWX*W + rnorm(N)) %>%
# True effect is 5
mutate(Y = 5*X + effectWY*W + rnorm(N))
return(d)
}
# Our estimation function
est_model <- function(N, effectWX, effectWY) {
d <- create_data(N, effectWX, effectWY)
# Biased estimator - no W control! But how bad is it?
m <- lm(Y~X, data = d)
return(coef(m)['X'])
}
# Iteration function! We'll add an option iters for number of iterations
iterate <- function(N, effectWX, effectWY, iters) {
results <- 1:iters %>%
map_dbl(function(x) {
# Let's add something that lets us keep track of how much progress we've made. Print every 100th iteration
if (x %% 100 == 0) {print(x)}
# Run our model and return the result
return(est_model(N, effectWX, effectWY))
})
# We want to know *how biased* it is, so compare to true-effect 5
return(mean(results) - 5)
}
# Now try different settings to see how bias changes!
# Here we'll use a small number of iterations (200) to speed things up, but in general bigger is better
iterate(2000, 0, 0, 200) # Should be unbiased## [1] 100
## [1] 200## [1] 0.0001086662iterate(2000, 0, 1, 200) # Should still be unbiased## [1] 100
## [1] 200## [1] 0.001805832iterate(2000, 1, 1, 200) # How much bias? ## [1] 100
## [1] 200## [1] 0.4988909iterate(2000, .1, .1, 200) # Now?## [1] 100
## [1] 200## [1] 0.009736692# Does it make a difference whether the effect is stronger on X or Y?
iterate(2000, .5, .1, 200)## [1] 100
## [1] 200## [1] 0.03934004iterate(2000, .1, .5, 200)## [1] 100
## [1] 200## [1] 0.05094061从这个模拟中我们能学到什么呢？首先，正如我们所预期的那样，当 \(effectWX\) 或 \(effectWY\) 设为\(0\)时，不存在偏差——估计值的均值减去真实值 \(5\) 几乎为\(0\)。我们还发现，如果 \(effectWX\) 和 \(effectWY\) 都为 \(1\)，偏差约为\(0.5\)——即误差项标准差的一半——并且大约是 \(effectWX\) 和 \(effectWY\) 都为 \(0.1\) 时偏差的 \(50\) 倍。偏差的大小与 \(effectWX\) 和 \(effectWY\) 的乘积有关吗？似乎有点这种关系……这是一个你可以通过更多模拟来探究的问题！
15.4 通过模拟进行功效分析
15.4.1 什么是功效分析？
在统计学领域，儿童电视节目传达的道理或多或少是正确的：只要你足够努力，一切皆有可能。同时，这个领域也或多或少印证了暴力电子游戏所蕴含的道理：如果你想解决一个极其棘手的问题，你需要投入大量的“火力”（此外，一台更快的计算机真的会大有帮助）。
功效分析27是试图弄清楚你所拥有的“火力”是否足够的过程。虽然有一些计算器和分析方法可用于进行功效分析，但功效分析通常通过模拟来完成，尤其是对于非实验性研究而言28。
为什么进行功效分析是个好主意呢？一旦我们确定了研究设计，有许多因素会使统计分析成为一个相当无力的工具，让我们更难发现真相。
- 真正微小的效应极难被发现——就像不借助超强力显微镜，你根本别想看到电子！（在公式 \(Y=\beta_1 X + \epsilon\) 中\(\beta_1\)₁的值极小）
- 大多数统计分析都围绕着观察数据的变异性展开。如果数据的变异性很小，我们就没什么可分析的了（\(X\) 的方差极小）。
- 如果除了我们所关注的效应之外，还有很多其他因素在起作用，那么就很难从噪声中提取出有用的信号（\(\epsilon\) 的方差很大）。
- 如果我们对所谓的”证据”设定了极高的标准，那么很多有价值的证据就会被忽略，这样一来，我们就需要更多的证据来弥补（如果你坚持要求 \(\beta_1\)的估计值 \(\hat \beta_1\) 的标准误极小）。
很巧的是，所有这些问题都可以通过增强”火力”来解决，这里所说的”火力”指的就是样本量。功效分析能帮助我们精确算出到底需要多大的”火力”。如果所需”火力”超出了我们的承受范围，那我们或许就得打道回府了。
如果你统计功效不佳且无法改善，你或许该考虑开展一项不同的研究。统计功效不佳带来的弊端可不仅仅是”我们做了研究却没得到显著结果——真倒霉”。统计功效不佳还会让你更有可能在发现统计效应时，出现效应被夸大或效应并非真实存在的情况！这是因为在统计功效不佳时，我们拒绝原假设基本上只是出于偶然。不妨这样想：如果你想在 \(95\%\) 的置信水平上拒绝“抛硬币正面和反面出现的概率各为 \(50\%\)”这一原假设，你需要看到一个在公平硬币的情况下出现概率低于\(5\%\)的结果。如果你只抛\(6\)次硬币（样本量极小，统计功效不佳），那么只有出现全是正面或全是反面的情况才会拒绝原假设。即便真实情况是正面出现的概率为\(55\%\)，即原假设确实不成立，但我们能拒绝“正面出现概率为\(50\%\)”这一原假设的唯一情况是，我们估计出正面出现的概率为\(0\%\)或\(100\%\)。而这两个估计值实际上都比我们刚刚拒绝的\(50\%\)更偏离真实情况！这就是为什么当你读到一项研究似乎发现了某一平常事物具有巨大影响时，该研究往往样本量极小（且无法重复验证）。
无论你进行的是哪种类型的研究，功效分析都是个不错的主意。不过，它在以下三种情况下尤其有用：
如果你要探寻的某种效应，你认为它可能并非当前所研究内容的核心或关键——它是一种微小的效应，或者处于一个有诸多其他因素在起作用的系统之中（与因变量 \(Y\) 的方差相比，自变量 \(X\) 的方差乘以回归系数 \(\beta_1\) 的值非常小），那么进行功效分析是个很棒的主意。要了解一个微小效应的有用信息，所需的样本量往往比你预期的大得多，现在了解这一点总比在你已经完成所有工作之后才发现要好。
如果你要探寻某种效应在不同组之间的差异情况（在回归分析中，如果你主要关注交互项），那么进行功效分析是个好主意。找出不同组在某种效应上的差异，比找出该效应本身需要更多的统计效力，你需要确保自己有足够的样本量。
如果你正在进行一项随机实验，那么功效分析是必不可少的——你实际上可以控制样本量，所以在确定样本量之前，你最好先弄清楚自己需要多大的样本量，以免样本量过小！
15.4.2 功效分析的作用
功效分析就像某种转盘子的马戏表演一样，需要在五件事情之间取得平衡。用 \(X\) 表示处理（干预），用 \(Y\) 表示结果，这五件事情分别是：
真实效应量（回归中的系数、相关性等）
处理（干预）的变异量（\(X\) 的方差，即 \(var(X)\) ）
\(Y\) 中的其他变异量（要么是用 \(X\) 解释 \(Y\) 后残差的变异，要么就是 \(Y\) 的方差，即 \(var(Y)\) ）
统计精度（估计值的标准误或统计功效，即真阳性率）
样本量
功效分析会固定其中四个因素，然后告诉你第五个因素的值会是多少。例如，它可能会说：“如果我们认为效应大小可能是\(A\)，且\(X\)的变异是\(B\)，与\(X\)无关的\(Y\)的变异是\(C\)，并且你希望统计精度达到\(D\)或更高，那么你至少需要\(E\)的样本量29。” 这告诉我们获得足够统计功效所需的最小样本量。
或者我们可以换个角度来看。“如果你打算采用样本量为\(A\)，并且变量\(X\)存在\(B\)程度的变异，同时与\(X\)无关的变量\(Y\)存在\(C\)程度的变异，而且你希望达到\(D\)或更高的统计精度，那么效应必须至少达到\(E\)这么大。”这就告诉了我们最小可检测效应30，也就是说，在给定你的样本量的情况下，你有合理机会能够测量到的最小效应。
你可以选择这五个要素中任意一个的最低可接受值。不过，人们最常选择的是最小样本量和最小可检测效应，紧随其后的是统计精度水平。
那个”统计精度”是怎么回事呢？具体来说，它指的是什么？通常，在功效分析中，统计精度是通过统计功效的目标水平来衡量的（这也就是“功效分析”这个名称的由来）。统计功效，也被称为真阳性率31，是假设检验中特有的一个概念。假设我们知道原假设是错误的。如果我们有一个规模无限大的样本，我们肯定会拒绝这个错误的原假设。然而，我们并没有无限大的样本。因此，由于抽样误差，我们有时会正确地拒绝原假设，而有时则会错误地未能拒绝它。我们正确拒绝错误原假设的比例就是真阳性率，也就是统计功效。如果我们在\(73\%\)的情况下拒绝了错误的原假设，这就相当于说我们有\(73\%\)的功效来拒绝该原假设3233。
统计功效也取决于你所进行的检验类型——如果你在\(95\%\)的置信水平下进行假设检验，相比于在\(99\%\)的置信水平下进行假设检验，你更有可能拒绝原假设（从而会有更高的统计功效）。你对假阳性结果越谨慎，就越有可能得到假阴性结果。所以这里存在一种权衡。
功效分析不一定非要从统计功效的角度来进行。事实上，你不一定非要从”拒绝原假设的能力”（这正是统计功效的核心所在）的角度去考虑问题。你可以以任何一种统计精度作为目标来进行功效分析，比如标准误差。已知\(A\)、\(B\)和\(C\)，你需要多大的样本量\(D\)才能使你的标准误差达到\(E\)或更小呢？
15.4.3 那些数字从何而来？
为了进行功效分析，你需要能够为那五个要素中的四个填入数值，这样功效分析才能告诉你第五个要素的值。我们如何知晓这些数值呢？
我们必须尽可能做出最佳猜测。我们可以利用以往的研究来大致了解效应可能有多大，或者变量\(X\)可能存在多大的变异性。如果没有以往的研究，那就尽你所能——思考可能的真实情况，做出有根据的猜测。在某种程度上，功效分析需要一些猜测。
其他方面并非基于猜测，而是基于标准。你希望达到多高的统计功效呢？统计功效越高，错过真实效应的可能性就越小。但这意味着你所需的样本量要大幅增加。人们通常在这里会采用经验法则。过去，达到\(80\%\)的统计功效是一个标准目标。如今，我更常看到人们将目标设定为\(90\%\)。
实际上，功效分析并非是一种极其精确的计算，它只是一种参考。它不需要做到丝毫不差。一定程度的猜测（尽管理想情况下应尽可能少）是必要的。毕竟，即使获得了所需的最小样本量，也不能保证你的分析就是可靠的，它只是让你在存在效应的情况下有相当大的机会得出结果。通常，人们会把功效分析的结果当作一个基准，然后确保实际样本量超过这个基准，他们的假设是自己在分析中可能过于保守了。所以，不必担心计算得是否精准，只要让分析中的数据足够接近实际情况，能发挥作用就行。
15.4.4 编写功效分析模拟代码
从本章的其他部分内容来看，进行功效分析模拟所需的大部分要素我们已经具备。“问题排查”部分的信息最为丰富，因为许多功效分析模拟在尝试不同操作时需要改变相当多的设置。
让我们来考虑一下我们已经讲过的三个要素：
数据生成。我们需要一个数据生成函数，它要能模拟出我们所认为的实际数据生成过程的样子。没错，又要回到绘制我们研究场景的因果图了！你已经完成这个步骤了，对吧？
你可以简化一些——也许不用考虑你计划测量和控制的 \(16\) 个变量以及 \(8\) 个你不打算处理的变量，你可以把它们全部合并为”控制变量”和”未测量变量”。但总体而言，你应该思考如何尽可能地生成与现实世界相似的数据。你还需要设法匹配一些情况，比如会有多大的差异——你预计处理变量会有多大的差异34？结果变量会有多大的差异？其他变量呢？你不可能做到完全精准，但要尽可能地接近，哪怕这只是一个有根据的猜测。
估计并存储结果。一旦我们获得了数据，就需要进行我们计划开展的分析。这和之前的操作方式完全一样——只需使用模拟数据进行你在研究中计划进行的分析即可！
唯一的不同之处在于，如果你打算使用统计功效作为统计精度的衡量标准。我们如何从分析中得出统计功效呢？
请记住，统计功效是具有统计学显著性的结果所占的比例（假设原假设为假，只要我们的数据生成过程存在非零效应，我们就知道原假设是假的）。因此，在模拟的每一次迭代中，我们需要提取一个二元变量，用以表明该效应是否具有统计学显著性，并将其输出。然后，在所有模拟中具有显著性的结果所占的比例就是你的统计功效。如果 \(45\%\) 的迭代产生了显著效应（且你的原假设为假），那么统计功效就是 \(45\%\)！
判断一个效应是否显著比提取系数或标准误差要稍微难一些。在 \(R\) 语言中，对于一个模型 \(m\)，你可以加载 \(broom\)包，然后执行\(d\leftarrow tidy(m)\) 。然后，如果 \(p\) 值低于 \(0.05\)（对于 \(95\%\) 的显著性检验），\(d\$p.value[2] < 0.05\) 将会返回 \(TRUE\)。这里的\([2]\)假设你感兴趣的系数是第二个（查看一次迭代的 \(tidy()\) 输出结果，以找到你感兴趣的系数所在的正确位置）。
在 \(Stata\) 中，你需要稍微深入了解一下支撑许多 \(Stata\) 命令的 \(Mata\) 矩阵代数。别担心，这不会太难。在进行回归分析之后，使用” \(matrix \space R = r(results)\)“语句提取结果矩阵。然后，使用” \(matrix \space p = R["pvalue", "X"]\)“语句获取\(p\)值。之后，你可以通过引用”\(p[1,1]\)“将该\(p\)值放置到你需要的地方。
值得庆幸的是，在\(Python\)里事情终于变得简单了。对于来自\(statmodels\)库的模型\(m\)，我们可以通过\(m.pvalues['X']\)获取变量\(X\)的系数的 \(p\) 值。可以将这个 \(p\) 值与（比如说）\(0.05\) 进行比较以进行显著性检验。如果该系数在\(95\%\)的置信水平下显著，那么\(m.pvalues['X']<0.05\)的结果为\(TRUE\)；反之则为\(False\)。
迭代。多次运行模拟的操作方式和往常一样，没什么变化！
不过，我们或许想在现有的迭代基础上再增加一层迭代35。比如说，我们想知道能达到\(90\%\)检验功效的最小可检测效应是多少。那我们要怎么找到这个值呢？我们可以尝试一系列不同的效应量，每次稍微增大一点，直到找到一个能产生 90%检验功效的效应量。而”尝试一系列不同的效应量”正是循环的用武之地。
同样地，如果我们想找到达到\(90\%\)检验功效所需的最小样本量，我们就得尝试一系列不同的样本量，直到检验功效达到\(90\%\)。
那我们就开始行动吧。我之前提到过，有一种情况检验功效会很低，就是当我们想了解两组之间效应的差异时，换句话说，就是回归中交互项的系数（其中一个参与交互的项是二分类变量）。
让我们来看看实际的检验功效到底有多低！假设我们想了解一项阅读培训项目对儿童阅读测试成绩的影响。此外，我们还想知道这种影响是否因性别而异。
在我们计划开展的研究中，我们知道测试成绩是特意设计成呈正态分布的（很多测试在设计时都会考虑这一点），其标准差为\(10\)36。
其他因素也会影响测试成绩，我们把很多这类因素打包成一个单一的变量，将其作为控制变量\(W\)纳入分析。控制变量\(W\)每变动一个标准差，测试成绩会提高\(4\)分。
实际上，这个阅读培训项目能使男孩的测试成绩提高\(2\)分，使女孩的测试成绩提高\(2.8\)分。这可是个相当大的差异——该项目对女孩的效果比对男孩的效果高出\(40\%\)！培训项目的分配不是随机的，控制变量\(W\)的值越高，接受培训的可能性就越低。
我们可以通过回归来估计培训对男孩和女孩的效果差异：
\[测试成绩 = \beta_0 + \beta_1 \times 培训 + \beta_2 \times 女孩 + \beta_3 \times 培训 \times 女孩 + \beta_4 \times W + \epsilon \quad (15.1)\]
在这个回归方程中，我们将重点关注 \(β^3\)
我们想知道为了有\(90\%\)的检验功效来检测男孩和女孩之间的效果差异所需的最小样本量。然后，假设我们不知道真实的效果差异，但知道可供抽样的儿童只有\(2000\)名，接着来询问在这个 \(2000\)个样本的情况下，我们能用\(90\%\)的检验功效检测到的最小可能效果差异是多少。让我们来编写代码。
library(generics)##
## Attaching package: 'generics'## The following object is masked from 'package:fixest':
##
## estfun## The following object is masked from 'package:lubridate':
##
## as.difftime## The following object is masked from 'package:dplyr':
##
## explain## The following objects are masked from 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, unionlibrary(broom)
library(vtable)
# power analysis
set.seed(1000)
# Follow the description in the text for data creation
# Since we want to get minimum sample size and minimum detectable effect, allow both sample size and effect to vary diff is the difference in effects between boys and girls
create_data <- function(N, effect, diff) {
d <- tibble(W = rnorm(N),girl = sample(0:1, N, replace = TRUE)) %>%
# A one-SD change in W makes treatment 10% more likely
mutate(Training = runif(N) + .1*W < .5) %>%
mutate(Test = effect*Training + diff*girl*Training + 4*W + rnorm(N, sd = 9))
return(d)
}
# Our estimation function
est_model <- function(N, effect, diff) {
d <- create_data(N, effect,diff)
# Our model
m <- lm(Test ~ girl*Training + W, data = d)
tidied <- tidy(m)
# By looking we can spot that the interaction term is in the 5th row
sig <- tidied$p.value[5] < .05
return(sig)
}
# Iteration function!
iterate <- function(N, effect, diff, iters) {
results <- 1:iters %>%
map_dbl(function(x) {
# To keep track of progress
if (floor(x/100) == x/100) {print(x)}
# Run our model and return the result
return(est_model(N, effect, diff))
})
# We want to know statistical power, # i.e. the proportion of significant results
return(mean(results))
}
# Let's find the minimum sample size
mss <- tibble(N = c(10000, 15000, 20000, 25000))
mss$power <- c(10000, 15000, 20000, 25000) %>%
# Before we had to do function(x) here, but now the argument we're passing is the first argument of iterate() so we don't need it
map_dbl(iterate, effect = 2, diff = .8, iter = 500) ## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500# Look for the first N with power above 90%
mss %>%
rename(Power = power) %>%
dftoHTML(out = 'viewer') ## NULLmss## # A tibble: 4 × 2
## N power
## <dbl> <dbl>
## 1 10000 0.602
## 2 15000 0.788
## 3 20000 0.882
## 4 25000 0.938library(generics)
library(vtable)
mde <- tibble(effect = c(.8, 1.6, 2.4, 3.2))
mde$power <- c(.8, 1.6,2.4, 3.2) %>%
map_dbl(function(x) iterate(N = 2000, effect = 2, diff =x, iter = 500)) ## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500mde %>%
rename(Power = power, Effect = effect) %>%
dftoHTML(out='viewer')## NULLmde## # A tibble: 4 × 2
## effect power
## <dbl> <dbl>
## 1 0.8 0.164
## 2 1.6 0.528
## 3 2.4 0.838
## 4 3.2 0.976从这些分析中我们能得到什么呢？通过最小样本量功效分析，我们得到了\(表15.1\)。从该表中我们可以看到，当样本量为 \(20000\) 时，功效为\(88.2\%\)；当样本量为 \(25000\) 时，功效为 \(93.8\%\)。所以，达到\(90\%\)功效的理想样本量肯定在这两个数值之间。我们可以尝试增加迭代次数，或者尝试\(20000\)到 \(25000\)之间选取不同的样本量数值，以获得更精确的答案。不过，“要使该交互项达到\(90\%\)的功效，需要略多于\(20000\)个观测值”，这样的信息本身已经相当不错了37。
表15.1：最小样本量功效分析模拟结果
等等，呃，两万个（样本量）？就为了那个庞大的交互项（你觉得针对所有人实施的一项干预措施，真正对某一个群体比对另一个群体有效\(40\%\)的情况会有多常见呢）？没错！“效果高出\(40\%\)” 远不如 “相对于测试标准差\(10\)而言为\(0.8\)” 重要。
如果我们只有 \(2000\)个观测值，对于真实交互效应为 \(0.8\) 的情况，我们的检验功效远达不到\(90\%\)，这并不奇怪。但对于我们 \(2000\) 的样本量，我们能检测到的最小可检测效应是多少呢？结果见\(表15.2\)。
在\(表15.2\)中，当效应值介于\(2.4\)到\(3.2\)之间的某个数值时，我们达到了\(90\%\)的检验功效。同样，我们可以尝试进行更多次迭代，并在\(2.4\)到\(3.2\)之间尝试更多的效应值，以获得更精确的结果，但这已经给了我们一个大致的范围。
如果我们只有\(2000\)个观测数据可用，为了能以\(90\%\)的检验功效检测出效应中的性别差异，我们需要女孩的效应达到男孩效应的约\(220\%\)（或者是男孩效应的\(40\%\)，但符号相反）。我们认为这是一个合理的假设吗？如果不合理，也许我们就不应该在这项研究中估计效应的性别差异，即便我们真的很想这么做。这就是功效分析的作用——知道什么时候仗已经输了，就该真的回家了（即适时放弃）38。
15.5 基于现有数据的模拟：自助法
说真的，这已经是我连续第三章都在讲自助法了？好吧，接受现实吧。自助法简直就是魔法。难道你不想让魔法走进你的生活吗39？
自助法\(（Bootstrapping）\)是指在有放回重抽样后多次重复进行分析的过程。有放回重抽样是指，如果你一开始有一个包含\(N\)个观测值的数据集，你通过从该数据集中随机抽取\(N\)个观测值来创建一个新的模拟数据集，但每次抽取一个观测值后，你都将其放回，这样你就有可能再次抽到它。如果你的原始数据集是\(1、2、3、4\)，那么一个自助法数据集可能是\(3、1、1、4\)，而另一个可能是\(2、4、1、4\)。
重采样过程模拟了实际数据中的抽样分布。平均而言，变量的均值将与您的数据均值一致，方差、变量间的相关性等统计特性也将保持一致。通过这种方式，我们可以直接使用真实存在的数据进行模拟。这种带替换的重采样过程取代了我们模拟中的\(create\_data\)函数。开展模拟研究最困难的部分之一就是如何生成接近现实世界的数据。如果我们能直接使用真实数据，那简直是天大的优势。
唯一的缺点在于：与自行生成数据的模拟不同，你无法获知真实的数据生成过程。因此，你无法通过自助法来——比方说——检验估计量是否存在偏差。虽然你能获得基于自助法模拟的估计值，但由于数据并非由你生成，你无从知晓真实参数值来进行比对，自然也就无法验证你的估计量是否平均而言准确无误！
这意味着，自助法作为模拟工具，其适用范围仅限于那些无需将结果与真实参数值进行对比的研究场景。通常情况下，这表现为自助法主要用于估计标准误——因为通过完整模拟重采样过程，该方法能自动复现变量间复杂的相互依存关系，使得真实抽样分布中的所有特殊形态都能自然渗透到自助法模拟的抽样分布之中。
不过，自助法的用途远不止于估计标准误！机器学习研究者将其广泛应用于各类场景，包括为模型生成大量训练数据集。您还可以通过自助法更全面地评估模型的变异性。具体而言，其典型工作流程包括：首先构建自助样本，接着估计模型参数，然后基于该模型生成预测值，如此循环往复。通过这种方式，您既能观察每个观测点的预测值如何随模型抽样变异而波动，也能汇总所有模拟预测结果来获得集成预测值4041。
自助法的实施步骤（参考\(第13章\)内容）如下：
数据准备：获取包含\(N\)个观测值的原始数据集
重抽样：从原始数据中有放回地随机抽取\(N\)个观测值（允许重复抽取），构成一个”自助样本”
统计量计算：基于该自助样本计算目标统计量
迭代重复：将步骤\(2-3\)42重复执行足够多次（通常\(\geq1000\)次）
分布分析：汇总所有迭代得到的统计量估计值形成经验分布。该分布的标准差即为统计量的标准误，其\(2.5\%\)和\(97.5\%\)分位数可构成\(95\%\)置信区间，以此类推。
\(步骤4\)中”足够多次”具体指代多少次？这需视情况而定：
代码测试阶段：执行数百次即可验证流程可行性（但不可采信其结果）
正式分析阶段：应在计算资源允许范围内最大化重复次数（通常需\(\geq 5000\)次）
重要研究建议：至少执行数千次——请预留充足计算时间！
自助法要获得准确的标准误估计，必须满足若干关键前提条件：
首先，样本容量需足够大（通常建议\(N\geq30\)）；
其次，数据分布应相对规范——若变量存在极端偏态，自助法的估计效果将大打折扣。
更重要的是，必须审慎处理观测值间的关联结构：例如面对面板数据（同一个体多次观测），应采用”聚类自助法”\(（cluster \space bootstrap）\)对个体而非单次观测进行重抽样（若机械地独立抽取\(2005年\)和\(2007年\)观测值，可能导致同一人在不同时点的数据被割裂）。当误差项存在聚类相关或自相关时，必须通过特定的重抽样策略（如块自助法）保持数据结构特征，否则会导致标准误低估。
15.5.1 编写自助法代码
为避免内容重复，本节将简明扼要。关于自助法标准误的具体代码实现，读者可参阅：\(第13章\)“回归分析”中内置的自助法标准误工具示例，以及\(第14章\)“匹配方法”中结合逆概率加权与回归调整的应用案例。
还剩下什么要做呢？在这一部分，我会向你展示如何做两件事：
一个关于如何手动执行自助法的通用演示。该示例将给出自助法标准误。对于\(\mathrm{R}\)和\(\mathrm{Stata}\)而言，借助大多数模型都具备的内置方法，计算自助法标准误会更轻松，但此处的代码能让你更清楚幕后的运算逻辑，还能轻松调整以实现标准误之外的功能（那些内置方法无法做到这一点），并且可用于处理没有内置自助法标准误计算方法的命令 。
在这次简短的回顾之后，我会深入讲讲自助法的一些变体
让我们手动执行一次自助法！过程非常简单。我们只需照搬在模拟代码中已经做过的所有操作。然后，把用于生成模拟数据的\(create\_data\)函数，替换为一个有放回地对数据重新抽样的函数 。
或者，我们可以省点事儿，不用专门为自助法流程编写独立函数 —— 通常这本来也就只需一行代码。要是我们有\(N\)个观测值，那就创建一个长度为\(N\)的向量，其中每个元素都是一个从\(1\)到\(N\)（在\(Python\)里则是从\(0\)到\(N-1\)）的随机整数。然后用这个向量来挑选要保留的观测值。在\(R\)和\(Python\)中，做法很简单：
# Load packages
library(tidyverse)
# Get our data set
data(mtcars)
# Resampling function
create_boot <- function(d) {
d <- d %>% slice_sample(n = nrow(d), replace = TRUE)
return(d)
# Or more "by-hand"
# index <- sample(1:nrow(d), nrow(d), replace = TRUE)
# d <- d[index ,]
}
# Get a bootstrap sample
create_boot(mtcars)## mpg cyl disp hp drat wt qsec vs am gear carb
## Pontiac Firebird...1 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Mazda RX4 Wag...2 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Pontiac Firebird...3 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Mazda RX4...4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Fiat X1-9...5 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Pontiac Firebird...8 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Pontiac Firebird...9 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Fiat 128...11 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Mazda RX4...14 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Camaro Z28...15 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2
## Fiat X1-9...18 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Mazda RX4 Wag...21 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Fiat 128...23 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Pontiac Firebird...25 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Camaro Z28...26 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Fiat 128...31 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3实现自助法其实非常简单！您可以直接将上述方法无缝嵌入常规模拟流程。不过在实际操作中，我们往往采取折中方案——手动编写核心估计方法并返回相应统计量，而将重抽样和迭代过程交由内置函数完成。三大统计语言\(（R/Python/Stata）\)均提供现成的自助法执行命令，其中\(Stata\)的相关功能已在\(第14章\)详述。
在所有这些场景中，我们首先需要定义一个估计函数。与先前代码示例不同，该函数不再调用数据生成函数来获取分析数据，而是通过以下协作机制运作：自助法函数将自动提供重抽样数据集作为输入，执行估计函数后捕获输出结果，最终汇总所有迭代的估计值。
# Load packages
library(tidyverse)
library(boot)
set.seed(1000)
# Get our data set
data(mtcars)
# Our estimation function. The first argument must be the data and the second the indices defining the bootstrap sample
est_model <- function(d, index) {
# Run our model
# Use index to pick the bootstrap sample
m <- lm(hp ~ mpg + wt, data = d[index ,])
# Return the estimate(s) we want to bootstrap
return(c(coef(m)['mpg'], coef(m)['wt']))
}
# Now the bootstrap call!
results <- boot(mtcars , est_model , 1000)
# See the coefficient estimates and their bootstrap standard errors
results##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = est_model, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -9.416845 -0.4034778 3.354603
## t2* -4.168012 -1.7372791 17.748946虽然操作复杂度并未降低，但使用这些专用自助法命令具有显著优势：它们通常提供多种高级选项，能实现比简单有放回抽样更复杂的重抽样方案，并能以比常规均值-标准差更精密的方式分析输出分布。
这些自助法命令最实用的优势在于支持偏差校正\(（bias \space correction）\)。此处的”偏差”并非指遗漏变量导致的因果效应估计偏差，而是特指数据分布形态对统计量估计的系统性影响。正如前文所述，自助法对偏态数据较为敏感——偏态分布会导致抽样分布同样产生偏斜，进而影响最终统计量及其标准误的估计精度。
针对此类偏差，最常用的校正方法是”\(BCa置信区间\)“\(（Bias-Corrected \space and \space Accelerated）\)。该方法通过量化以下两个关键指标来修正偏差：\(（1）\)自助法估计值低于原始估计值的比例（无偏分布下应为\(50\%\)）43；\(（2）\)自助分布的偏斜程度。基于这些指标计算得到的\(BCa置信区间\)，能够有效校正自助抽样分布的偏斜影响44，进而可用于假设检验。
在 \(R\) 语言中，将 \(boot()\) 函数的计算结果传递给 \(boot.ci(type='bca')\) 即可获取\(BCa置信区间\)。\(Stata\) 的实现方式为：在 \(bootstrap\) 命令中添加 \(bca\) 选项，运行完成后需额外执行 \(estata \space bootstrap，bca\) 命令（根据使用经验，该流程可能偶尔会出现异常）。\(Python\) 则需从 \(resample.bootstrap\)模块导入并调用 \(confidence\_interval()\) 而非 \(bootstrap\)。
15.5.2 自助法变体
正如几乎所有统计方法一样，自助法也顽固地拒绝成为单一概念，而是衍生出无数变体。我一直将演示的方法称为”自助法”而非更精确的”配对自助法\(（pairs \space bootstrap）\)“，这种做法肯定已经惹恼了某些人——虽然我实在不确定具体是谁45。
自助法的变体远比我在这里介绍的要多得多，但我将介绍其中几种，特别是那些针对非独立数据的变体。还记得我说过自助法只有在观测值互不相关时才有效吗？确实如此！至少对于基础版本是这样的。但如果我们遇到异方差性呢？或者自相关？又或者观测值是聚类在一起的？这并不是说自助法就失效了，而是意味着我们的自助程序必须考虑到观测值之间的关系。
集群/区块自助法\(（cluster/block \space bootstrap）\)似乎是个不错的入门选择——毕竟我之前已经有所提及。这种方法和常规自助法如出一辙，唯一的区别在于：我们不再对单个观测值进行重抽样，而是对观测值的集群（或称”区块”）进行重抽样。仅此而已！
当您遇到分层数据或面板数据时，就需要考虑使用集群自助法了46。分层数据指的是观测值以嵌套形式存在于不同组别中——比如学生考试成绩数据：学生们嵌套在班级里，班级嵌套在学校中，学校又归属于地区，地区最终从属于国家。
具体是否需要采用集群自助法，取决于您的分析需求（毕竟从某种角度看，所有数据都具有层级性）。但是，如果层级结构本身就是模型的一部分——比如模型中包含班级层面的特征变量，那么对单个观测值进行自助抽样就可能不太合理了。
并非所有层级数据集都需要集群自助法，但对于面板数据而言，这往往是理想选择。在面板数据中，我们会对同一个体（个人、企业等）进行多次观测。我们将在\(第16章\)对此展开更多讨论。在这种情况下，对单个观测值进行重抽样并无太大意义，因为这相当于将每个人分割开来，只对其部分观测值进行重抽样。因此，我们需要对个体及其所有观测值进行抽样，而非对每个观测值单独抽样。
集群/区块自助法另一个常见的应用场景是时间序列数据47。在时间序列中，我们会对同一变量进行多期观测，这通常意味着数据存在自相关性。若直接随机重抽样，会破坏时间维度上的相关结构。更棘手的是，当模型包含滞后项时，若恰巧未能抽到前一期数据，整个模型就会出问题！
在处理时间序列数据时，自助法的应用需要先将时间序列分割为连续的时间区块48。这些区块的划分通常由某种最优区块长度算法确定——比如对\(2000-2010年\)的每日股票指数数据，可能按月划分区块，或每六周为一个新区块，亦或按年划分，甚至采用不等长区块。关键在于确保每个区块内保持时间连续性，随后对这些完整区块（而非单个数据点）进行重抽样49。
在 \(Stata\) 中实施集群/区块自助法，只需在常规 \(bootstrap\) 命令后添加 \(cluster()\) 选项即可。若使用 \(R\) 语言，当您需要回归模型的\(bootstrap\)标准误时，可通过 \(sandwich \space 包\) 中的 \(vocvBS()\) 设置 \(cluster\) 参数（将 \(vcovBS()\) 结果传入 \(modelsummary\) 表格的 \(vcov\) 参数）。不过，\(R\) 和 \(Python\) 似乎都没有类似 \(boot包\) 那样适用于任意函数的通用实现方案。\(Python\)目前看来完全没有集群自助法的实现，但 \(arch包\) 的 \(arch.bootstrap\) 模块提供了多种区块自助法的变体。
所幸在 \(R\) 和 \(Python\) 中手动实现集群自助法并不复杂，只需三步：
\(（1）\)创建仅包含个体\(ID\)值的数据集；
\(（2）\)对该数据集执行常规自助重抽样；
\(（3）\)通过内连接（ \(R\) 中的 \(inner\_join()\) 或 \(Python\) 的 \(d.join(how='inner')\) 将抽样结果与原数据合并，这会自动复制每个个体对应的所有观测值。该方法同样适用于区块自助法，但前提是您需要自行定义区块边界——这通常不如算法划分来得精准。
# Load packages
library(tidyverse)
set.seed(1000)
# Example data
d <- tibble(ID = c(1,1,2,2,3,3), X = 1:6)
# Now, get our data set just of IDs
IDs <- tibble(ID = unique(d$ID))
# Our bootstrap resampling function
create_boot <- function(d, IDs) {
# Resample our ID data
bs_ID <- IDs %>% sample_n(size = n(), replace = TRUE)
print(bs_ID)
# And our full data
bs_d <- d %>% inner_join(bs_ID, by = 'ID', relationship = "many-to-many")
return(bs_d)
}
# Create a cluster bootstrap data set
create_boot(d, IDs)## # A tibble: 3 × 1
## ID
## <dbl>
## 1 3
## 2 2
## 3 3## # A tibble: 6 × 2
## ID X
## <dbl> <int>
## 1 2 3
## 2 2 4
## 3 3 5
## 4 3 5
## 5 3 6
## 6 3 6另一种广受欢迎的自助法实际运作方式大不相同。在残差重抽样中，你并非对整个观测值进行重抽样，而是对残差进行重抽样。这是怎么回事呢？
在残差重抽样流程下，需遵循以下步骤：
用原始数据开展分析。
借助分析得到预测结果值\(\hat{Y}\)，以及残差\(\hat{\epsilon}\)（实际结果与预测结果的差值，即\(\hat{\epsilon}=Y - \hat{Y}\) ）。
执行自助重抽样以获取重抽样后的残差\(r_b\) 。
把残差加到实际结果上，构建新的结果变量（即\(Y_b = Y + r_b\) ）。
用\(Y_b\)按常规方式估计模型。
存储你想要的任意结果，然后重复步骤 \(3-5\)。
这种方法的核心理念在于：模型中的预测变量始终保持不变，因此它们之间的任何依存关系都能被完整保留。但残差重抽样法存在一个明显缺陷——当误差项与预测变量存在任何相关性时，其表现就会欠佳。这一假设条件甚至比普通最小二乘法的要求更为严格（OLS仅要求误差项与预测变量均值无关）
在此类方法中，还有一种既拥有酷炫名称又能规避误差无关假设的流行变体——野自助法\(（wild \space bootstrap）\)。其优势在于：即便存在异方差性（且具体形式未知）时依然适用，同时对聚类数据也表现良好。虽然已有集群/区块自助法处理聚类数据，但当集群规模差异悬殊时（比如一个班级\(100\)人而另一个仅\(5\)人），传统方法会很快失效。此时基于集群的野自助法表现更优——在集群数量较少时（“较少”甚至可多达\(50\)个）同样适用50。
在不深入技术细节的前提下，我们简要说明野自助法的核心原理：它遵循残差重抽样的基本框架，但关键区别在于——并非直接重抽样残差，而是在每次生成自助数据时，始终保持各残差与原观测值的对应关系。
野自助法的核心操作在于：先通过转换函数处理残差，再乘以一个均值为\(0\)的随机变量（还需满足其他条件）——而非直接重抽样残差51。具体而言，对每个观测值生成随机数，乘以转换后的残差，由此构建新因变量\(Y_b\)，用其拟合模型并重复该过程。
需要特别说明的是，转换函数的选择和随机变量的分布类型存在多种方案——不同情境下各方法的性能表现各异，因此无法给出普适性建议。若准备采用野自助法，建议根据数据类型深入调研适用方案。如需了解技术细节，推荐参考\(Martin（2017）\)的权威指南。
Michael A. Martin. Wild Bootstrap, pages 1–6. John Wiley & Sons, Ltd., 2017.ISBN 9781118445112.野自助法何以有效看似令人费解，但其核心机理在于：通过随机变量与残差的乘积，能有效模拟抽样分布。更重要的是，这种方法完整保留了原始模型特性，使待估参数的抽样分布更具一致性。
在 \(R\) 语言中，可通过 \(sandwich包\) 的 \(vcovBS()\) 的 \(type\) 参数实现野自助法标准误计算。该参数支持多种随机乘数分布类型，默认‘\(wild\)’选项采用 \(Rademacher\)分布（即等概率取\(-1\)或\(1\)的离散分布）。\(Stata\) 用户可在回归命令后接 \(boottest包\) 中的 \(boottest\) 命令实现。目前 \(Python\) 尚未有成熟的实现方案。
以及其他一些较少被讨论的假设。↩︎
这对读者来说是幸事还是损失，这取决于你（或者也许取决于你的教授）。我这么做是有理由的，市面上已经有大量基于证明的计量经济学书籍供你研读。那些基于证明的书籍很有价值，但我不确定我们是否迫切需要再来一本这样的书。这涉及到边际价值递减、比较优势等等原理。要知道，我毕竟还是个经济学家。↩︎
这是否算个问题，因人而异。我认为做证明很有启发性，并且擅长进行计量经济学证明让我对数据的运行机制有了更好的理解。但为了能把某种方法当作工具运用自如，是否有必要了解该方法背后的证明过程呢？我对此表示怀疑。↩︎
然后再进行几次。↩︎
用模拟的方法向他人证明一个“优秀”估计量的价值更为困难。这是模拟方法更广泛弊端的一个方面——我们只能进行有限次数的模拟，因此我们只能表明该估计量在所选择的特定数据生成过程中是有效的。而理论证明能够更广泛、更具体地告诉我们某种方法在何时或何地可行或不可行。↩︎
或者，类似地，如果问题的复杂性在于数据而非估计量呢？每人有多个观测值会如何影响标准误差？或者一个组嵌套在另一个组内呢？又或者，如果不同组在不同时间接受处理，这会影响标准误差吗？↩︎
我向你保证确实如此。一旦你掌握了它，你就会一直进行模拟操作！它们可不只是有用，还让人上瘾。↩︎
一方面，很抱歉（可能给你带来了困扰）。另一方面，要是你还没学会做这件事（编写代码进行模拟），即便你这辈子再也不进行其他模拟了，花时间搞明白它也是非常非常值得的。在我写这段话的当下，已经是21世纪20年代了。计算机操控着你周围的一切，而且很可能是你日常经常接触的最强大的工具。难道你不想知道如何指挥这些工具做事吗？我不确定怎样才能有效地让你明白，在当今这个时代，哪怕你并不那么擅长编程，稍微学一学编程以及如何熟练使用计算机，都能实实在在地让你拥有一定的能力。这肯定很有吸引力，对吧？↩︎
如果你要较真的话，这一章讲的全是伪随机数据。编程语言中的大多数随机数生成器实际上是相当可预测的，毕竟，要知道，计算机本身就很容易预测。它们会从一个“种子”开始。将这个种子输入一个函数，生成一个新的数字，同时还会生成一个新的种子，用于生成下一个数字。就我们的目的而言，伪随机就足够了。甚至可能更好，因为这意味着我们可以自行设置种子，从而使我们的分析具有可重复性，也就是说，每个人运行我们的代码都应该得到相同的结果。↩︎
重要的是，几乎任何随机分布都可以通过使用随机均匀数据并通过反条件分布函数来生成。或者，对于分类数据，你可以这样做：“如果\(X>0.6\)，则\(A=2\)；如果\(X<=0.6 且 X>=0.25\)，则\(A=1\)；否则\(A=0\)”，其中\(X\)均匀分布在\(0\)到\(1\)之间。这等价于创建一个分类\(A\)，其中\(2\)的概率为\(0.4\times(1-0.6)\) ，\(1\) 的概率为\(0.35\times(0.6-0.25)\) ，\(0\)的概率为\(0.25\)。↩︎
想在\(Stata\)或\(Python\)中处理具有两个以上类别且概率不相等的情况吗？试试像前一条注释中那样对均匀分布数据进行分割操作。↩︎
是的，我们实际上可以看到误差项，通常是无法观察的！创造真相的美妙之处就在于此。这甚至可以让你做一些事情，比如，检查你模型中的残差与误差项的匹配程度。↩︎
以及其他相关内容。阅读文档以了解其他信息。在 R 语言中是使用
help(Distributions)，在 Stata 中是help random，在 Python 中是help(numpy.random)。我一直说要阅读文档，因为这真的是很好的建议！↩︎如果要提取其他内容，你可能需要更深入地查找。不过，通常在网上搜索“如何从LANGUAGE的ANALYSIS中获取STATISTIC”就能找到答案。↩︎
除了可能用于存储结果之外——我们之后会涉及这方面内容。↩︎
Stata 16 有同时使用多个数据集的方法，但我不想假设你安装了 Stata 16，而且我自己也从未学过这一功能，我猜想大多数 Stata 用户都在随心所欲地使用
preserve/restore命令。我们也可以尝试将结果存储在矩阵中，但如果你能避免使用 Stata 矩阵，你的日子会过得更舒心。↩︎哎呀，我从 2008 年左右就开始在 Stata 里运行模拟了，直到现在没有用它（
simulate命令）我也应付过来了。↩︎如果你愿意，你可以不用 broom 包，而是使用
coef(summary(m))['X', 'Std. Error']来获取标准误差。但 broom 包能让提取回归对象的很多其他部分变得容易得多，所以我建议你习惯使用它。↩︎这本书里的所有代码示例都没有使用 C 语言，但还望您见谅。↩︎
顺便说一下，生成包含对撞因子的数据是一种能让你直观地确信，对撞因子在未加控制的情况下的存在不会使你产生偏差的好方法。 例如，为了生成数据，一个带有路径\(X \rightarrow C \leftarrow Y\) 的对撞 \(C\)，你必须首先生成 \(X\)，然后生成带有 \(X\) 的任何影响的 \(Y\)，最后使用 \(X\) 和 \(Y\) 一起生成 \(C\)。 你甚至可以在生成 \(C\) 之前停止，仍然可以对 \(Y\) 进行关于 \(X\) 的回归。 仅仅创建变量肯定不会改变系数。↩︎
你大概已经预感到，这方法不会奏效。至少，如果它真管用，我可能早在这本书里提过了，对吧？但你知道为什么 它行不通吗？先坐下来琢磨一会儿。等你想出一个可能的解释……就可以用模拟来验证你的猜想！你会怎么设计一个模拟实验，来检验你的解释是否正确呢？↩︎
当然，你所权衡的偏差量和精度会因你的数据情况而异。所以你可能想使用你自己的实际数据来进行模拟！在这种情况下，你并不知道“真实”值，但你可以检查平衡性来大致了解偏差情况。关于如何使用实际数据进行模拟，可参考下面的自助法部分内容。↩︎
阿古斯的估计量并非一种改进，原因如下：首先，它舍弃了数据，从而减小了有效样本量。其次，虽然方差确实处于标准误差的分母位置，但普通最小二乘法残差却在分子位置！任何有损预测效果的操作，比如选取极端值数据，都会损害精度并增大标准误差。此外，选取极端值数据本身就更具噪声。↩︎
在进行这项工作时，请务必牢记，就像对待所有事情一样，某些假设的重要性可能会因具体情境而异。如果你在模拟中发现违背某个假设并没有那么严重，或许可以再次进行模拟，这一次依然违背相同的假设，但使用外观上截然不同的数据。如果结果仍然良好，那就是个好迹象。↩︎
这有点主观，并且要基于具体情况而定 —— 某种药物在导致致命副作用方面存在 0.1% 的偏差可能是非常严重的，但一个就业培训项目在增加收入方面存在 0.1% 的偏差可能就没那么严重。↩︎
这种整体方法是一种“敏感性分析”——检验结果对于假设违背的敏感程度。这也是一种部分识别方法，正如\(第11章\)所讨论的那样。这种用于思考剩余未闭合后门的敏感性分析方法被广泛使用，并且当你不确定是否存在未闭合后门时，它可用于确定一系列合理的回归系数。例如，可参考奇内利（Cinelli）和黑兹利特（Hazlett）2020年的研究。
Carlos Cinelli and Chad Hazlett. Making sense of sensitivity: Extending omitted variable bias. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(1):39– 67, 2020.
↩︎功效分析：一种将样本量、效应量和估计值的统计精度联系起来的分析方法，其目标通常是确定样本量是否足够大，以便在你所假设的效应量准确的情况下，有可能拒绝原假设。↩︎
功效分析计算器：必须针对你的估计和研究设计进行特定设置。在实验中这是可行的，因为一种简单的”我们将样本随机分为两组并测量一个结果”的设计涵盖了数千个案例。但对于观测数据，你要处理几十个变量，其中一些还未被测量，让计算器开发者提前猜到你的意图是不可行的。话虽如此，计算器对于某些因果推断设计也很适用，我在\(第20章\)讨论了一款我推荐用于断点回归设计的计算器。↩︎
在给定一定的效应大小、处理变异和与结果无关的变异的情况下，能产生特定统计精度的最小样本量。↩︎
在给定一定的样本量、处理变异和无关结果变异的情况下，能产生一定统计精度的最小效应量。↩︎
统计功效，即真阳性率。 当原假设实际上为假时，拒绝原假设的概率。↩︎
统计功效针对的是我们试图拒绝的原假设。 要拒绝一个完全荒谬的原假设相当容易，但要拒绝一个接近真实情况、但实际上并非真实的原假设，则需要强大得多的“火力”。↩︎
功效分析并不一定与假设检验相关。统计功效仅仅是衡量统计精度的一种方式。例如，我们可以转而询问，要使标准误差降至某个特定水平，我们需要多大的“效力”。↩︎
如果你有观测数据，也许你可以利用研究设定之前的旧数据来估计这些数据以及数据的其他特征。↩︎
谁来迭代这些迭代器？↩︎
我打算用 10 作为依据来创建测试分数，方法是将我们所有其他的决定因素相加，然后再加上一个标准差为 9 的正态分布误差。为什么是 9 而不是 10 呢？因为我们相加得到测试分数的其他组成部分会使标准差增大。这样做不会使标准差恰好达到 10，但会足够接近。↩︎
而且，如果所有这些（情况）都达到了 90%的功效（或者都没达到），那就表明我们需要尝试更小（更大）的样本量。这就是我确定此处要进行迭代的样本量范围的方法。↩︎
交互作用存在一些内在的统计问题，这些问题会降低检验功效，但真正的问题在于，实际上交互项的效应往往很小。治疗效果在不同组之间存在差异，但这种差异通常比主效应本身小得多。因此，无论我们需要多少观测数据来发现主效应，我们都需要更多的数据来发现那些微小的交互效应。用统计学方法来发现微小的效应是非常困难的！↩︎
我在这里与之交流的这个虚构人物真是个不知感恩的人。↩︎
这一方法被命名为”装袋法”\(（bagging）\)——这个术语我一方面觉得挺有意思，但另一方面听起来实在像某种\(1991\)年前后流行的、如今回想起来超级呆板的舞蹈风潮。↩︎
需要特别注意：若计划使用自助法分布进行假设检验，必须先将该分布的中心位置调整至未经过自助法的原始估计值，再将其与零假设分布进行比较。这种调整能确保假设检验的对比基准正确无误。↩︎
若数据或统计量服从极端非正态分布（如均值未定义的分布，参见\(第22章\)），则需在应用前对自助法分布的尾部进行截断处理。↩︎
除自助法外，研究者还常使用其近亲方法——“刀切法”（\(jackknife\)，与\(第13章\)提到的交叉验证同源）。该方法的实现流程是：对于包含\(N\)个观测值的样本，进行\(N\)次迭代计算，每次严格剔除一个观测值，从而确保每个观测点恰好被排除一次。在因果森林中，该方法主要用于：\(（1）\)评估单个观测点对估计结果的影响力；\(（2）\)为当前被排除的观测生成样本外预测值。↩︎
有趣的是，它们并非围绕估计值对称 —— 例如，若偏度导致过多自助法估计值低于（位于）非自助法估计值的左侧，那么我们对分布的左侧会有更多了解，此时\(95\%\)置信区间无需延伸过远；但在分布的右侧，我们则需要将区间延伸得更远一些。↩︎
关于自助法的研究文献可谓深不见底——究竟该采用哪种形式？搭配何种重抽样方法？选择怎样的随机生成过程？又该进行哪些调整？
这些讨论不仅艰深晦涩，还处处可见这样的警告：“本方法仅在特定情境下表现良好”。可惜入门门槛实在太高，若阁下执意要深入探索——那就祝君好运了！↩︎按照惯例，这种情况下我们称之为”集群自助法”而非”区块自助法”。↩︎
而这一次，按照惯例我们称之为”区块自助法”而非”集群自助法”——随你怎么理解吧↩︎
Dimitris N Politis and Halbert White. Automatic block-length selection for the dependent bootstrap. Econometric Reviews, 23(1):53–70, 2004.
↩︎另一种常用方法是”滑动区块自助法”（moving block bootstrap），其原理是采用滑动窗口进行重抽样，而非将数据分割为固定区块。↩︎
常规聚类标准误（clustered standard errors）依赖于渐近假设——即理论上当聚类数量趋近无穷大时才成立。实践中，对于聚类误差而言，50个聚类通常就被视为足够接近”无穷大”。但野聚类自助法（wild cluster bootstrap）无需这种假设，因此在聚类数量较少时表现更优。↩︎
均值为零的特性意味着：这些随机变量中会有许多负值——部分观测值的残差符号会发生翻转。↩︎