14 匹配
14.1 关闭后门的另一种方法
回归分析受到了很多关注，但关闭后门的方法不止这一种。在本书的前半部分，我对此非常坚持：选择一个在变量\(W\)中没有变化的样本，就可以关闭\(W\)依赖的任何开放的后门1。这听起来很简单，而且确实是一个直观的想法。
” 我们选择了两个看起来很可比的群体，并将它们进行了比较 “，这比” 我们做出了一些线性假设并最小化了误差平方和，以便调整这些后门变量 ” 更容易向你的老板解释。
2但一旦你真的开始尝试进行匹配，很快就会意识到，在挑选可比组这类操作上，存在一堆关于具体怎么做的问题。这并不显而易见！而且，做匹配的方法不同，结果可能大相径庭。肯定有聪明的人已经思考这个问题有一阵子了吧？没错，确实有。
匹配\(（Matching）\)是通过构建在一组匹配变量上具有相似性的对照组，来阻断处理变量和结果变量之间混淆路径的过程。这种方法通常适用于二元的处理 / 未处理变量，因此你要挑选一个与恰好接受”处理”的组非常相似的”对照”组34。
举个最基本的匹配案例：假设我们想评估职业培训项目对获取优质工作的影响。在符合项目资格的失业人群中，男女比例原本各占\(50\%\)，但项目宣传却主要面向男性，导致实际参与者中男性占\(80\%\)、女性仅\(20\%\)。鉴于劳动力市场表现存在性别差异，此时就形成了明显的后门路径：\(就业结果 \leftarrow 性别 \rightarrow 职业培训项目\)。
匹配方法的操作方式是：从未接受培训群体中筛选出一个同样符合\(80\%\)男性/\(20\%\)女性比例的控制组，与原本就呈现\(80-20\)性别比例的处理组进行对比。此时，处理组和对照组之间的性别差异被完全消除——“\(性别 \rightarrow 职业培训项目\)”的路径就此切断，后门关闭，我们就能准确识别出”\(职业培训 \rightarrow 就业结果\)“这一我们真正关心的因果效应5。
匹配与回归是阻断后门路径的两种不同方法。尽管二者都要求所选取的控制/匹配变量足以阻断所有后门路径，但它们依赖的其他假设却存在本质差异（并无优劣之分）。举例而言：
匹配方法可轻松规避线性假设——这一假设在回归分析中却至关重要（如回归章节所述，我们不得不依赖或费力规避该假设）
匹配在估计目标处理效应时更具灵活性（如\(第10章\)所示），能更自由地获取特定类型的平均处理效应
需要强调的是：两种方法不存在绝对优劣，实际应用中可组合使用，而且不同假设也能形成互补优势（本章后续将展开讨论）
唉，如果匹配法纯粹是个更差的选择，我们大可直接忽略它；若它明显更优，之前的回归章节也大可不写。但现实偏偏是两者各有千秋，所以我们最好都掌握。当然，决定采用匹配法时还有大量细节需要厘清——这也正是本章要探讨的核心内容6。
14.2 加权平均值
首先，我们究竟在尝试做什么来匹配？正如我所说，我们希望使我们的处理组和对照组具有可比性。但这究竟意味着什么呢？
匹配方法会为每个观测值生成一组权重（记作\(W\)），这些权重的设计目的就是让实验组和对照组具备可比性。随后，当我们需要估计处理效应时，可以分别计算处理组和对照组结果的加权平均值，再进行比较7。
正如前几章所述，加权平均的计算方法是：将每个观测值乘以其权重后相加，再除以权重之和。以数值\(1，2，3，4\)为例——常规的非加权均值其实也是一种加权均值，只不过每个值的权重均为\(1\)。具体计算如下：\(1*1=1， 2*1=2， 3*1=3， 4*1=4\)；将这些结果相加：\(1+2+3+4=10\)；最后用总和除以权重之和\(1+1+1+1=4\)，得到\(10/4=2.5\)8。
现在我们用不同的权重再算一次。假设权重分别为\(0.5，2，1.5和1\)。计算过程相同：首先，每个数值乘以对应权重：\(1*0.5=0.5，2*2=4， 3*1.5=4.5，4*1=4\)；然后相加：\(0.5+4+4.5+4=13\)；最后除以权重总和 \(（0.5+2+1.5+1=5）\)，得到加权平均数 \(13/5=2.6\)。
我们可以将\(Y\)的加权平均值表示为方程：
换句话说，将每个值\(Y\)乘以其权重\(W\)，把它们全部加起来（\(\sum\)），然后除以所有权重的总和\(\sum w\) 。你可能注意到，要是所有权重都是\(w = 1\) ，这就和计算一个普通的平均值一回事 。
那么这些权重从何而来？其实存在多种不同的匹配方法，每种方法生成权重的路径各不相同。但它们都基于一组”匹配变量”，并通过这些变量构建权重体系，从而阻断这些变量所对应的后门路径。其核心思想是：通过权重调整，最终使实验组和对照组在匹配变量上不再存在任何差异。
在上一节中，我举了一个男女比例\(80-20\)的例子。现在让我们更精确地剖析这个案例。假设我们有一个接受过职业培训的实验组，组内包含\(80\)名男性和\(20\)名女性。在这\(80\)名男性中，\(60\)人最终成功就业，\(20\)人未就业；女性成员中，\(12\)人就业，\(8\)人未就业。现在来看由\(500\)名男性和\(500\)名女性组成的对照组。其中男性成员有\(350\)人就业、\(150\)人未就业；女性成员则有\(275\)人就业、\(225\)人未就业。
如果直接比较原始数据：接受职业培训的群体中，\((60+12)/100=72\%\)的人成功就业；而对照组中\((350+275)/1000=62.5\%\)的人就业，得出处理效应为\(9.5\)个百分点——这个结果看起来不错！但这里存在一个通过性别变量的后门路径：不同性别的平均就业率存在差异，而性别又与是否接受职业培训相关。
有很多不同的匹配方法，每种方法可能会生成不同的权重，但有一种方法可能会生成如下这样的权重：
给每个接受处理的个体权重为 \(1\)
给所有未接受处理的男性权重为 \(80/500=0.16\)
给所有未接受处理的女性权重为 \(20/500=0.04\)
用这些权重，咱们看看是否消除了组间性别的差异。处理组依旧会是\(80\%\)男性 —— 给所有接受处理的个体相同权重，不会改变这一侧的情况。那未接受处理的个体呢？要是我们计算男性比例，男性赋值为\(1\)，女性赋值为\(0\)，他们的加权均值是\((0.16×500 + 0.04×0)/(0.16×500 + 0.04×500) = 80\%\) 。完美！9
现在，我们已经平衡了处理组和对照组之间的性别差异，那么处理效应是多少呢？处理组中最终找到工作的人依旧占\(72\%\)—— 同样，这部分没变化。但在未处理组中，比例是\((0.16×350 + 0.04×275)/(0.16×500 + 0.04×500) = 67\%\)10，处理效应就是\(0.72 - 0.67 = 5\)个百分点。虽说不如我们之前得到的\(9.5\)个百分点那么好，但也还不错。之前的部分收益是因为存在通过性别的混淆路径（后门路径 ）。
14.3 概念匹配：单一匹配变量
匹配方法有多种实现方式，具体包括哪些呢？为便于演示，我们将从单一变量匹配入手分析。虽然实际应用中很少仅匹配单个变量，但这种方法更易于理解，能帮助我们厘清核心概念——至于多变量匹配的重要细节，我将在后续章节专门阐述。
通过具体案例能更直观地理解。我们以信用卡违约数据为例：采用\(2005年\)台湾地区名\(30000\)持卡人的月度账单及\(4-9月\)还款记录（按时还款、延迟还款等。本研究聚焦信用问题的”粘性”效应，即分析 \(4月\) 逾期还款\(LateApril\)对 \(9月\) 逾期还款\(LateSpet\)的影响。虽然存在诸多未控制的混杂因素，但关键要控制\(4月\)账单金额\(BillApril\)这一后门变量——该变量既影响\(4月\)也可能影响\(9月\)的还款行为。
Moshe Lichman. UCI machine learning repository,2013.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(extrafont)## Registering fonts with Rlibrary(modelsummary)
library(ggpubr)
library(gghighlight)
library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:ggpubr':
##
## get_legend
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(Matching)## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## ##
## ## Matching (Version 4.10-15, Build Date: 2024-10-14)
## ## See https://www.jsekhon.com for additional documentation.
## ## Please cite software as:
## ## Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching
## ## Software with Automated Balance Optimization: The Matching package for R.''
## ## Journal of Statistical Software, 42(7): 1-52.
## ##library(vtable)## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(haven)
library(conflicted)
library(dplyr)uci <- read_csv('Matching/UCI_Credit_Card.csv', show_col_types = FALSE) %>%
mutate(LateSept = PAY_0 > 0,
LateApril = PAY_6 > 0,
BillApril = BILL_AMT6/1000) %>%
dplyr::select(LateSept,LateApril,BillApril, AGE)
head(uci)## # A tibble: 6 × 4
## LateSept LateApril BillApril AGE
## <lgl> <lgl> <dbl> <dbl>
## 1 TRUE FALSE 0 24
## 2 FALSE TRUE 3.26 26
## 3 FALSE FALSE 15.5 34
## 4 FALSE FALSE 29.5 37
## 5 FALSE FALSE 19.1 57
## 6 FALSE FALSE 20.0 37在进行单变量匹配时，我们需要做出以下关键选择：
匹配标准如何设定？
采用样本选择匹配法还是构建加权匹配样本？
若选择样本匹配，匹配数量如何确定？
若构建加权样本，权重随距离衰减的规则如何设置？
可接受的最差匹配质量阈值是多少？
\(（1）\)匹配标准如何设定？在进行匹配时，我们通常需要从对照组中筛选与处理组相似的观测样本。但”相似性”的具体定义是什么？要明确匹配目标，就必须首先确立相似性的判定标准11。
确定匹配标准主要有两种方法：距离匹配\(（distance \space matching）\)与倾向得分匹配\(（propensity \space score \space matching）\)。距离匹配法的核心原则是：“当观测对象的匹配变量值相近时，即视为相似”。该方法通过最小化处理组与对照组协变量之间的差异程度来实现匹配。
距离匹配的合理性直观明确：通过强制确保处理组与对照组在匹配变量上高度相似，从而有效阻断后门路径。如\(图14.1\)所示（该图示与经典后门路径模型一致），这正是距离匹配的典型应用场景——当我们能够测量后门路径上的关键变量时，通过变量匹配即可实现因果识别。
以信用卡违约数据为例，我们选取\(1\)个处理组样本（\(4月\)逾期还款）：第\(10305\)行观测对象。该持卡人\(4月\)账单金额为\(91630\)新台币，其\(4月\)还款状态为逾期\(（LateApril=true）\)，但\(9月\)还款正常\(（LateSept=false）\)。
uci[10305,]## # A tibble: 1 × 4
## LateSept LateApril BillApril AGE
## <lgl> <lgl> <dbl> <dbl>
## 1 FALSE TRUE 91.6 23由于我们的匹配变量是\(BillApril\)，我们要找未接受处理、且\(BillApril\)值与新台币\(91630\)元非常接近的匹配观测值。一个\(BillApril\)为新台币\(113023\)元（距离为\(\vert91,630 - 113,023\vert = 21,393\) ）或为\(0\)元（距离为\(\vert91,630 - 0\vert = 91,630\) ）的对照组，不会是这个接受处理观测值的理想匹配。另一方面，\(BillApril\)值为\(91613\)元（距离\(\vert91,630 - 91,613\vert = 17\)）的个体，可能就是非常好的匹配，因为他们之间的距离很小。
要是我们为第\(10305\)行选一个单独的匹配对照观测值12，会选第\(27281\)行，它\(4月\)没有逾期（所以未接受处理 ），且\(BillApril\)为新台币\(91690\)元（距离\(\vert91,630 - 91,609\vert = 21\)），在未接受处理的观测值里，这是数据中最接近新台币\(91630\)元的。我们能在\(图14.2\)里直观看到这个匹配。我们有来自第 \(10305\)行的接受处理观测值，你能在图中心附近看到它。匹配的是在\(X\)轴（我们的匹配变量 ）上最接近它的未接受处理观测值。第\(27281\)行的观测值就在左边，非常近。再往左稍远一点的是第\(27719\)行的观测值，我们很快会说到它。
set.seed(1000)
sub <- uci %>%
mutate(id = row_number()) %>%
mutate(BillApril = BillApril*1000) %>%
# the distance is within 200
dplyr::filter(abs(BillApril - uci$BillApril[10305]*1000) < 200) %>%
# convert the LateSept variable to numerical form and add random perturbations
mutate(yjit = jitter(as.numeric(LateSept), amount = .5)) %>%
# standardize the BillApril variable
mutate(BA_std = (BillApril - mean(BillApril))/sd(BillApril)) %>%
mutate(blank = '')
refba <- sub %>% dplyr::filter(id == 10305) %>% pull(BA_std)
# obtain weights using the Epanechnikov function
epan <- function(x, r = 1) {
kernel <- (3/4) * (1 - r * x^2)
pmax(kernel, 0)
}
sub <- sub %>% mutate(wt = epan(BA_std - refba, r = 1))
head(sub)## # A tibble: 6 × 9
## LateSept LateApril BillApril AGE id yjit BA_std blank wt
## <lgl> <lgl> <dbl> <dbl> <int> <dbl> <dbl> <chr> <dbl>
## 1 TRUE FALSE 91829 37 2097 0.828 2.03 "" 0
## 2 FALSE FALSE 91459 31 2111 0.259 -1.20 "" 0
## 3 FALSE FALSE 91753 26 2947 -0.386 1.37 "" 0
## 4 FALSE FALSE 91825 43 3877 0.191 2.00 "" 0
## 5 FALSE TRUE 91783 39 4049 0.0164 1.63 "" 0
## 6 FALSE TRUE 91534 33 10098 -0.432 -0.544 "" 0.222ggplot(sub) +
geom_point(aes(x = BillApril, y = yjit, shape = LateApril)) +
geom_hline(aes(yintercept = .5), linetype = 'dashed') +
geom_text(
aes(x = BillApril, y = yjit, label = scales::number(id, big.mark = ',')),
data = sub %>% dplyr::filter(id %in% c(10305, 27281, 27719)),
family = 'sans',
size = 10/.pt,
hjust = -.2) +
labs(x = 'Bill in April',shape = 'Late in April') +
scale_shape_manual(values = c(1,19), labels = c('Not Late\nin April','Late in\nApril')) +
scale_y_continuous(breaks = c(0,1), labels = c('Not Late in\nSeptember','Late in\nSeptember')) +
scale_x_continuous(labels = function(x) paste0('NT$',scales::number(x, accuracy = 1, big.mark = ',')))+
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position.inside = c(.7, .8),
legend.background = element_rect(color = 'black')
)
ggsave('Matching/Matching/matched_values.pdf', width = 6, height = 5,device=cairo_pdf)uci %>%
mutate(id = row_number()) %>%
mutate(BillApril = BillApril*1000) %>%
# obtain the control group that was not overdue in April
dplyr::filter(!LateApril) %>%
# based on distance matching
mutate(dif = abs(BillApril - 91630)) %>%
arrange(dif) %>%
# 27281/27719
slice(1:2)## # A tibble: 2 × 6
## LateSept LateApril BillApril AGE id dif
## <lgl> <lgl> <dbl> <dbl> <int> <dbl>
## 1 FALSE FALSE 91609 37 27281 21
## 2 TRUE FALSE 91596 24 27719 34另一种主流匹配方法是倾向得分匹配法\(（Propensity \space Score \space Matching）\)。其核心思想是：若观测对象接受处理的概率相同（即具有同等处理倾向性），则视为相似样本13。
Marco Caliendo and Sabine Kopeinig. Some practical guidance for the implementation of propensity score matching. Journal of Economic Surveys, 22(1): 31–72, 2008.与距离匹配类似，倾向得分匹配同样具有理论合理性。其核心不在于消除匹配变量的所有变异，而在于准确识别处理变量对结果变量的因果效应——关键在于阻断匹配变量构成的后门路径1415。倾向得分匹配严格遵循这一原则，认为只要确保处理倾向性匹配，即可满足因果识别条件16。
如\(图14.3\)所示，匹配变量\(A、B、C\)均位于后门路径上，但这些路径都经过 处理倾向性\(（Treatment Propensity）\)——即接受处理的概率。虽然处理概率本身不可直接观测，但通过匹配变量对处理变量的回归分析，我们可以直接估计该倾向得分17。
回到信用卡案例：当使用 \(BillApril\)（单位：千新台币）作为自变量进行\(logit\)回归时，得到截距项为\(0.091\)，\(BillApril\)系数为\(0.0003\)。若为第\(10305\)行处理组样本\(（BillApril=91630新台币）\)寻找匹配对象：
变量标准化：\(91630\)新台币 → \(91.630\)（千新台币）
处理概率预测：\(logit(0.091+0.0003*91.630)=0.116\)18。
# a logit regression of treatment on in thousands of New Taiwan dollars
m <- glm(LateApril ~ BillApril, data = uci, family = binomial(link = 'logit'))
msummary(m, stars = TRUE, out='markdown', gof_omit = 'AIC|BIC|F|Lik')| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| (Intercept) | -2.290*** |
| (0.023) | |
| BillApril | 0.003*** |
| (0.000) | |
| Num.Obs. | 30000 |
| RMSE | 0.30 |
ind <- coef(m)[1] + coef(m)[2]*91.630
pv <- exp(ind)/(1+exp(ind))
pv## (Intercept)
## 0.1158829接下来计算所有对照组样本的预测处理概率，筛选与处理组概率值\(（0.116）\)最接近的匹配对象。结果显示，第\(27281\)行对照样本的预测概率同样为\(0.116\)——这堪称完美匹配19！
uci <- uci %>%
mutate(ind = coef(m)[1] + coef(m)[2]*BillApril) %>%
mutate(pred = exp(ind)/(1+exp(ind)),pred2 = predict(m, type='response'))
all.equal(unname(uci$pred), unname(uci$pred2))## [1] TRUEuci %>%
mutate(id = row_number()) %>%
# obtain the control group that was not overdue in April
dplyr::filter(!LateApril) %>%
mutate(dif = abs(pred - pv)) %>%
arrange(dif) %>%
slice(1)## # A tibble: 1 × 9
## LateSept LateApril BillApril AGE ind pred pred2 id dif
## <lgl> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 FALSE FALSE 91.6 37 -2.03 0.116 0.116 27281 0.00000605\(（2）\)我们是在挑选匹配个体，还是在构建一个匹配加权样本？我有几次在”寻找匹配”的情境中提及匹配过程，还有几次提到了权重的构建。这两种方式都代表了匹配处理组和对照组的不同方法。
两者都不是默认方法——它们各有优缺点，进行匹配必然意味着要在不同选项间做一些选择，没有一种方法能完全优于另一种。
样本选择式匹配的过程就是决定将哪些对照样本纳入最终匹配样本：匹配度达标则入选，不达标则淘汰。所有入选样本均被赋予相同权重20。
在上一节中，我们挑选了一个匹配。我们查看了第\(10305\)行的处理组观测值，然后注意到第\(27281\)行的对照组观测值，在匹配变量的距离以及倾向得分值方面，都是最接近的。要是我们只选一个匹配来 “纳入”，就会选第\(27281\)行，该行会得到权重\(1\)。其他所有行都 “不纳入”，权重为\(0\)（至少在它们不是其他处理组观测值的最佳匹配时是这样）。
但这只是一种方法。我们还能怎么确定哪些对照组纳入、哪些不纳入呢？要了解这个，你得看步骤\(（3）\) 。
另一种构建匹配加权样本的方法怎么样呢？和挑选匹配个体的过程一样，这种过程需要查看所有对照组观测值，看它们与处理组观测值的接近程度。不过，每个对照组观测值不会简单地被纳入或排除，而是会根据其与处理组观测值的接近程度，或者在使匹配对照组看起来与处理组相似方面起到多大作用，获得不同的权重。
cv <- cov(uci %>% dplyr::select(BillApril, AGE) %>% mutate(BillApril = BillApril*1000))
cv## BillApril AGE
## BillApril 3.546692e+09 26137.64855
## AGE 2.613765e+04 84.96976通常（虽非绝对），构建匹配加权样本的方法，会让观测值离某一给定处理组观测值越近，权重就越大；离得越远，权重就越小。比如，我们发现第\(10305\)行和第\(27281\)行在匹配变量上的差值是\(\vert91,630 - 91,609\vert = 21\) 。而第二好的对照组匹配在第\(27719\)行，差值为\(\vert91,630 - 91,596\vert = 34\) 。我们可能想把这两个样本点都纳入，但给差值小的观测值\(（27281行）\)的权重，比差值大的观测值\(（27719行）\)的权重要高。或许我们给第\(27281行\)权重\(\frac{1}{21}\) ，给第\(27719行\)权重\(\frac{1}{34}\)21。
uci[27281,]## # A tibble: 1 × 7
## LateSept LateApril BillApril AGE ind pred pred2
## <lgl> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FALSE FALSE 91.6 37 -2.03 0.116 0.116uci[27719,]## # A tibble: 1 × 7
## LateSept LateApril BillApril AGE ind pred pred2
## <lgl> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 TRUE FALSE 91.6 24 -2.03 0.116 0.116我们究竟要如何构建这些权重呢？这是步骤\(（4）\)要解决的问题。
我们可以在 \(图14.4\) 中对比这两种方法。左边，我们寻找与处理组观测值非常接近的匹配。这种情况下，我们进行一对一匹配，所以只选第\(27281\)行的观测值，也就是最接近的匹配，将其视为 “纳入”。其他所有观测值都被置灰，完全不计入。右边，我们用所有未处理的数据来构建针对第\(10305\)行的比较估计（其他处理组观测值在此处被排除，但之后我们得回过头来为它们也做匹配 ）。但这些数据并非都被平等计入。在\(X\)轴上靠近我们处理组观测值的观测值，会被大量计入（大的气泡 ）。随着距离变远，它们被计入的程度越来越小。远到一定程度，权重就降到零。
# 定义需要高亮的ID
highlight_ids <- c(10305, 27281)
# 提前筛选需要高亮的数据点
highlight_data <- sub %>% dplyr::filter(id %in% highlight_ids)
# 创建基础绘图
p1 <- ggplot(sub, aes(x = BillApril, y = yjit)) +
# 添加散点，按LateApril区分形状
geom_point(aes(shape=LateApril)) + # 添加透明度区分高亮/非高亮点
# 直接使用filter后的数据添加高亮点
geom_point(aes(shape=LateApril),data=highlight_data, size=3) +
# 为高亮点添加标签
geom_text(
aes(label = scales::number(id, big.mark = ",")),
data = highlight_data,
family = "sans",
size = 10 / .pt,
hjust = -0.2,
) +
# 添加参考线
geom_hline(yintercept = 0.5, linetype = "dashed") +
# 设置坐标轴和标题
labs(x = "Bill in April",shape = "Late in April", title = "(a) One-to-one Matching") +
# 自定义图形元素
scale_shape_manual(values = c(1, 19),labels = c("Not Late\nin April", "Late in\nApril")) +
scale_y_continuous(breaks = c(0, 1),labels = c("Not Late in\nSeptember", "Late in\nSeptember")) +
scale_x_continuous(labels = function(x) paste0("NT$", scales::number(x, accuracy = 1, big.mark = ",")),breaks = c(91600, 91800)) +
# 应用主题
theme_pubr() +
theme(
text = element_text(size = 10, family = "sans"),
axis.title.x = element_text(size = 10, family = "sans"),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = c(0.7, 0.8),
legend.background = element_rect(color = "black")
)## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# 基础绘图：使用过滤后的数据 (排除 LateApril=TRUE 的记录，但保留 id=10305 的记录)
p2 <- ggplot(sub %>% dplyr::filter(!LateApril | id == 10305), aes(x = BillApril, y = yjit)) +
# 散点层：点的大小映射到 wt 变量，形状映射到 LateApril 变量
geom_point(aes(size = wt, shape = LateApril)) +
# 高亮显示：突出显示 wt > 0 的数据点，用 bank 变量作为标签
gghighlight::gghighlight(wt > 0, label_key = blank, use_group_by = FALSE) +
# 水平参考线：在 y=0.5 处添加虚线
geom_hline(aes(yintercept = .5), linetype = 'dashed') +
# 文本标注：为特定 id (10305) 添加数值标签
geom_text(aes(x = BillApril, y = yjit, label = scales::number(id, big.mark = ',')),
data = sub %>% dplyr::filter(id %in% c(10305)),
family = 'sans',
size = 10/.pt,
hjust = -.4) +
# 标签设置：x轴、图例和标题
labs(x = 'Bill in April', shape = 'Late in April', title = '(b) Kernel Matching') +
# 图例控制：隐藏 size 和 shape 图例
guides(size='none', shape='none') +
# 形状映射：设置形状值及其标签
scale_shape_manual(values = c(1,19), labels = c('Not Late\nin April','Late in\nApril')) +
# y轴设置：刻度位置和标签
scale_y_continuous(breaks = c(0,1), labels = c('Not Late in\nSeptember','Late in\nSeptember')) +
# x轴设置：格式化标签为货币格式，设置刻度位置
scale_x_continuous(labels = function(x) paste0('NT$',scales::number(x, accuracy = 1, big.mark = ',')), breaks = c(91600,91800)) +
# 主题应用：使用 pubr 主题
theme_pubr() +
# 主题调整：字体、轴标题等细节
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()
)## Warning in ggrepel::geom_label_repel(mapping, layer$data, fill = "white", :
## Ignoring unknown aesthetics: shapeplot_grid(p1,p2)
ggsave('Matching/Matching/select_vs_kernel.pdf', width = 7.5, height = 3.5,device=cairo_pdf)那么该如何选择呢？采用样本选择法还是差异化权重法22？这个问题没有标准答案——若有，我也无需同时介绍两种方法了。
样本选择法具有显著的直观优势：通过精准筛选个体构建可比对照组，不仅实施更简便，还能清晰追踪匹配过程。该方法还能避免极端权重问题——设想\(10305\)行与\(27281\)行差异仅为\(0.000001\)（而非\(21\)），该对照将获得\(1/0.000001\)的巨额权重，导致对照组均值几乎完全由该样本决定。虽然存在解决方案，但此类风险需时刻警惕。
另一方面，加权样本法具有某些优良的统计特性。它的敏感性较低——“入选”或”淘汰”是非此即彼的刚性区分，因此像卡钳值（后续会讨论）或测量误差等微小变化都可能显著改变结果，这使得样本选择法的结果波动更大。差异化的权重还能很好地处理以下事实：不同观测值在数据中能找到的匹配质量本身就有优劣之分。
当使用倾向得分时，必须采用加权样本法。尽管”入选/淘汰”式匹配颇为流行，但该方法实际上可能扩大处理组与对照组在匹配变量上的差异——因为不同变量对倾向得分的复合影响，可能导致高值变量相互匹配的失真情况。
Gary King and Richard Nielsen. Why propensity scores should not be used for matching. Political Analysis, 27(4):435–454, 2019.\(（3）\)若我们在挑选匹配个体，选多少个？要是我们决定采用挑选匹配个体的方法，就需要确定为每个处理组观测值挑选多少个对照组匹配个体。
这里有三种主要方法：挑选最佳匹配（一对一匹配 ）、挑选排名前\(K\)的最佳匹配（\(K\)近邻匹配 ），或者挑选所有可接受的匹配，也叫半径匹配，因为会接受在可接受半径范围内的所有匹配。
这些流程的运作方式基本和你预想的一样。一对一匹配时，为每个处理组观测值挑选单个最佳对照组匹配个体。\(K\)近邻匹配时，挑选最佳对照组匹配个体…… 以及第二佳、第三佳…… 直到选够\(K\)个匹配个体（或用完潜在匹配个体 ）。而挑选所有可接受匹配，就是确定什么样的匹配是可接受的\(（见步骤5）\)，然后用所有可接受匹配进行匹配。
\(1\)对\(1\)匹配：一种匹配流程，每个处理组观测值仅匹配一个对照组观测值，不同于”一对多”匹配（每个处理组观测值可能匹配多个对照组）。
\(K\)近邻匹配：一种匹配流程，为每个处理组观测值选取\(K\)个最佳对照组匹配个体作为匹配。
对于这些流程中的每一种，我们还需要决定是有放回匹配还是无放回匹配。要是某个对照组观测值恰好是两个不同处理组观测值的最佳匹配（或相应地，是\(K\)个最佳匹配之一 ），该怎么办呢？如果是无放回匹配，那个对照组只能匹配其中一个处理组观测值，另一个处理组观测值就得找其他对照组。如果是有放回匹配，我们可以多次使用同一个对照组，给它一个权重，数值等于它被匹配的次数。
我们该如何在这些不同选项间做选择呢？这很大程度上要在偏差和方差间做权衡。
每个处理组观测值对应的匹配数越少，这些匹配能越好。显然，最佳匹配比第二佳匹配更优。比如，对比\(1对1匹配\)和\(2近邻匹配\)，\(2近邻匹配\)会纳入对照组中一些与处理组并非那么相近可比的观测值。这会引入偏差，因为你没法那么笃定地说已经真正阻断了后门路径。
另一方面，每个处理组观测值对应的匹配数越多，对对照组均值的估计就越不易受随机波动影响，处理效应估计也就越精确。要是有\(100\)个处理组观测值，\(100\)个匹配对照组观测值的均值，比\(200\)个匹配对照组观测值的均值，抽样分布范围更广。就这么简单！匹配数越多，抽样差异越小，标准误也就越低。
那么，选择进行多少个匹配，要基于偏差和精度在你的估计中有多重要，以及要是尝试进行更多匹配，匹配质量会变得多差。要是你有无数个对照组观测值，而且能为每个处理组观测值挑选一大堆超级好的匹配，那么允许一大堆匹配可能不会引入太多偏差，还能通过这样做降低方差。那就这么做！但要是你的对照组很小，你的第三好匹配可能是个相当差的匹配，会引入很多偏差。不值得这么做。
关于有放回 / 无放回匹配的选择呢？这里也涉及偏差 / 方差的权衡。有放回匹配确保每个处理组观测值都能使用其最佳（或前\(K\)佳，或所有可接受的 ）匹配。这会减少偏差，因为同样，这种方法让我们能挑选最佳匹配。然而，这种方法意味着我们会反复使用相同的对照组观测值 —— 每个对照组观测值对均值的影响更大，所以抽样差异会更大（要是一个非常好的匹配在一个样本中与\(30\)个处理组匹配，但在另一个样本中不存在，会怎样呢？ ）。极端情况下，想象只有一个对照组观测值，把它与所有处理组匹配 —— 对照组的样本均值就只是那个观测值的结果值，标准误就等于结果的标准差。要是\(N = 1\) ，除以\(\sqrt{N}\) 来计算标准误没什么用。
有替换匹配除了偏差较低外还有另一优势：顺序无关性。假设进行无替换匹配时，处理观测\(1\)和\(2\)的最佳匹配都是对照组\(A\)，但二者的次优匹配分别是\(B\)和\(C\)。此时谁将获得\(A\)？若观测\(1\)匹配\(A\)，则\(C\)进入对照组；若观测\(2\)匹配\(A\)，则\(B\)进入对照组。对照组构成（进而估计结果）取决于我们如何”分配”优质对照组\(A\)。虽然存在某些分配原则（如将\(A\)分配给次优匹配较差的处理观测），但本质上仍具有一定任意性。
\(（4）\)若我们在构建匹配加权样本，权重随距离如何衰减？匹配的两种主要方法 —— 挑选样本或构建权重，在执行方式上都要做一些重要选择。在匹配加权样本方法中，我们要采用距离度量或倾向得分，并用其为每个观测值单独加权。但我们如何依据距离或倾向得分来转化为权重呢？
我们又有几种选择。使用权重的两种主要方法，一方面是核匹配，或更宽泛地说，基于核的匹配估计量；另一方面是逆概率加权。
基于核的匹配估计量使用核函数生成权重。在匹配背景下，核函数接收差异值并返回权重——最大值位于\(0\)差异处，随后向两侧平滑递减，最终归零并保持为\(0\)。该方法对优质匹配赋予高权重，次优匹配权重较低，劣质匹配权重为\(0\)23。
理论上存在无限多种核函数可供选择，其中\(Epanechnikov核\)尤为常用24。该核函数的优势在于计算简便25：
其中\(K(x)\)表示”核函数”，该函数仅在\(x \in [-1,1]\)区间内有效，超出区间则取\(0\)。\(图14.5\)直观展示了\(Epanechnikov核函数\)的图像。
# epan <- function(x, r=1) (3/4)*(1-r*(x^2)) %>% purrr::map(function(x) max(x,0)) %>% unlist()
epan <- function(x, r = 1) {
kernel <- (3/4) * (1 - r * x^2)
pmax(kernel, 0)
}
x <- seq(-2, 2, by = 0.1)
plot(x, epan(x), type = "l", main = "Epanechnikov Kernel")
需注意其硬编码范围——该核函数仅在差异值\(x \in [-1,1]\)区间有效。因此标准做法是：先将距离标准化，再输入核函数计算26。这可确保不会因变量尺度差异而导致匹配结果不同。
sd(uci$BillApril)## [1] 59.55411sd(uci$AGE)## [1] 9.217904x/sd(uci$BillApril)## [1] -0.033582906 -0.031903761 -0.030224615 -0.028545470 -0.026866325
## [6] -0.025187180 -0.023508034 -0.021828889 -0.020149744 -0.018470598
## [11] -0.016791453 -0.015112308 -0.013433162 -0.011754017 -0.010074872
## [16] -0.008395727 -0.006716581 -0.005037436 -0.003358291 -0.001679145
## [21] 0.000000000 0.001679145 0.003358291 0.005037436 0.006716581
## [26] 0.008395727 0.010074872 0.011754017 0.013433162 0.015112308
## [31] 0.016791453 0.018470598 0.020149744 0.021828889 0.023508034
## [36] 0.025187180 0.026866325 0.028545470 0.030224615 0.031903761
## [41] 0.033582906x/sd(uci$AGE)## [1] -0.21696906 -0.20612061 -0.19527216 -0.18442370 -0.17357525 -0.16272680
## [7] -0.15187834 -0.14102989 -0.13018144 -0.11933298 -0.10848453 -0.09763608
## [13] -0.08678762 -0.07593917 -0.06509072 -0.05424227 -0.04339381 -0.03254536
## [19] -0.02169691 -0.01084845 0.00000000 0.01084845 0.02169691 0.03254536
## [25] 0.04339381 0.05424227 0.06509072 0.07593917 0.08678762 0.09763608
## [31] 0.10848453 0.11933298 0.13018144 0.14102989 0.15187834 0.16272680
## [37] 0.17357525 0.18442370 0.19527216 0.20612061 0.21696906当核函数生成权重后，您可据此计算加权均值27。不过前文提及的是”基于核的匹配估计量”——还存在更复杂的方法：先通过核函数构建基础，再采用更精细的方式估计均值。这些方法我将在后续”匹配数据估计”章节详述。
在我们的信用卡债务示例中，这是如何运作的呢28？我们已经计算了两个差值。第 \(10305\) 行和第 \(27281\) 行的 \(BillApril\) 匹配变量差值为\(\vert91,630 - 91,609\vert = 21\) ，第 \(10305\) 行和第 \(27719\) 行的差值为\(\vert91,630 - 91,596\vert = 34\) 。接下来我们进行标准化 ——\(BillApril\)的标准差是\(59994\)，这样这两个差值就分别转化为\(\frac{21}{59,994} = 0.00035\) 和\(\frac{34}{59,994} = 0.00056\) 个标准差。接着我们用 \(Epanechnikov 核函数\) 处理这些值，得到权重分别为\(\frac{3}{4}(1 - 0.00035^2) = 0.999999\) 和\(\frac{3}{4}(1 - 0.00056^2) = 0.999997\) （显然，核函数认为这两个对照组几乎一样好 ）。
现在，我们要对所有对照组重复这个过程，但如果假设只有这两个对照组足够接近、有影响，那么我们将第 \(10305\) 行的处理后平均结果\(LateSpet\)值 \(0（false）\)，与第 \(27281\) 行的 \(0\) 值和第 \(27719\) 行的 \(1\) 值进行比较，用我们的权重将它们平均，得到\(\frac{0.999999×0 + 0.999997×1}{0.999999 + 0.999997} = 0.4999999\) 。最后我们得到处理效应为\(0 - 0.4999999 = -0.4999999\) 。
set.seed(1000)
sub <- uci %>%
mutate(id = row_number()) %>%
mutate(BillApril = BillApril*1000) %>%
# the distance is within 200
dplyr::filter(abs(BillApril - uci$BillApril[10305]*1000) < 200) %>%
# convert the LateSept variable to numerical form and add random perturbations
mutate(yjit = jitter(as.numeric(LateSept), amount = .5)) %>%
# standardize the BillApril variable
mutate(BA_std = (BillApril - mean(BillApril))/sd(BillApril)) %>%
mutate(blank = '')
refba <- sub %>% dplyr::filter(id == 10305) %>% pull(BA_std)
# obtain weights using the Epanechnikov function
sub <- sub %>% mutate(wt = epan(BA_std - refba, r = 1))
head(sub)## # A tibble: 6 × 12
## LateSept LateApril BillApril AGE ind pred pred2 id yjit BA_std
## <lgl> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 TRUE FALSE 91829 37 -2.03 0.116 0.116 2097 0.828 2.03
## 2 FALSE FALSE 91459 31 -2.03 0.116 0.116 2111 0.259 -1.20
## 3 FALSE FALSE 91753 26 -2.03 0.116 0.116 2947 -0.386 1.37
## 4 FALSE FALSE 91825 43 -2.03 0.116 0.116 3877 0.191 2.00
## 5 FALSE TRUE 91783 39 -2.03 0.116 0.116 4049 0.0164 1.63
## 6 FALSE TRUE 91534 33 -2.03 0.116 0.116 10098 -0.432 -0.544
## # ℹ 2 more variables: blank <chr>, wt <dbl>x <- uci %>% dplyr::select(BillApril, AGE) %>% mutate(BillApril = BillApril*1000) %>% slice(10305,27281) %>% as.matrix()
x <- x[2,] - x[1,]
x## BillApril AGE
## -21 14那么逆概率加权是怎样的呢？逆概率加权可追溯至\(Horvitz 和 Thompson（1952）\)的研究，经 \(Hirano、Imbers 和 Ridder（2003）\)发展而来，它是专门为倾向得分设计的。它会用观测值实际接受处理情况的概率的倒数来对每个观测值加权。所以，要是你实际接受了处理，且估计的倾向得分为\(0.6\) ，那你的权重就是\(\frac{1}{0.6}\) 。要是你实际没接受处理，且估计的倾向得分为\(0.6\) ，那你的权重就是1 除以你未接受处理的概率，也就是\(\frac{1}{1 - 0.6} = \frac{1}{0.4}\)
Horvitz, Daniel G., and Donovan J. Thompson. 1952. “A Generalization of Sampling Without Replacement from a Finite Universe.” Journal of the American Statistical Association 47 (260): 663–85. Hirano, Keisuke, Guido W. Imbens, and Geert Ridder. 2003. “Efficient Estimation of Average Treatment Effects Using the Estimated Propensity Score.” Econometrica 71 (4): 1161–89.我们用 “实际所处状态的概率的倒数” 来加权，让处理组和对照组更相似。处理组中权重最大的观测值，是那些最像未处理组的 —— 倾向得分小（所以倾向得分的倒数大 ）、本最不可能接受处理却接受了处理的观测值。同理，对照组中权重最大的观测值，是那些最像处理组的、本最可能接受处理却因某些原因没接受处理的观测值。
逆概率加权有几个不错的优势。除了指定倾向得分估计回归之外，你不需要做太多选择，这一点很好。而且，你完全可以不进行任何形式的匹配就进行逆概率加权 —— 无需将每个处理组观测值与每个对照组观测值进行比对，只需估计倾向得分，你就已经知道权重了。此外，正如\(Hirano、Imbers 和 Ridder（2003）\)所表明的，只要你有足够大的样本量，并且有灵活的方法来估计倾向得分，加权就是估计因果效应的最精确方式29。
这也有缺点，不足为奇。尤其是，有时会出现非常意外的处理或未处理情况，进而权重会变得极大。想象一下，有人倾向得分是\(0.999\) ，结果却未被处理 —— 千分之一的概率（但在大样本中，有时可能发生 ）。那个人的权重会是\(\frac{1}{1 - 0.999} = 1000\) 。这真的很大！只要倾向得分可能接近\(0\)或\(1\)，这种情况就会让逆概率加权不稳定。
有解决这些问题的办法。最常见的就是直接 “修剪” 数据，去掉任何倾向得分过于接近\(0\)或\(1\)的观测值。但这可能会以一种难以修正的方式导致标准误出错。其他办法包括，先把倾向得分转化为优势比（\(\frac{p}{1 - p}\) ，而非仅仅\(p\)），然后在处理组和对照组内分别对权重进行缩放，使它们的总和都为\(1\)\(（Busso、DiNardo 和 McCrary, 2014）\)。
Matias Busso, John DiNardo, and Justin McCrary. New evidence on the finite sample properties of propensity score reweighting and matching estimators. Review of Economics and Statistics, 96(5):885–897, 2014.让我们把逆概率加权应用到信用卡债务示例中。之前我们用逻辑回归模型计算出，第 \(10305\) 行的处理组观测值接受处理的概率是\(0.116\) ，第 \(27281\) 行的对照组观测值也是如此。我们还可以加入第 \(27719\) 行，它的倾向得分也是\(0.116\)（毕竟我们就是因为它们相似才选的。由于\(BillApril\)值如此相似，预测的处理概率差异要到小数点后第五位才会显现 ）。
那么每个观测值会得到怎样的权重呢？第 \(10305\) 行实际接受了处理，所以我们给该观测值的权重是 \(1\) 除以接受处理的概率，即\(\frac{1}{0.116} = 8.621\) 。它被赋予很大权重，因为它原本不太可能被处理，所以它看起来像未处理组，却接受了处理！对照组的行都会得到 \(1\) 除以未处理概率的权重，即\(\frac{1}{1 - 0.116} = 1.131\) 。它们原本就可能不被处理，所以不像处理组，权重不会太大。由此，我们可以得到处理组和对照组的加权均值。
\(（5）\)可接受的最差匹配是怎样的？匹配的目的是通过挑选与处理组观测值最为相似的对照组观测值，尝试让处理组和对照组具有可比性。但要是存在某个处理组观测值，没有任何对照组观测值与它有丝毫相似之处，那该怎么办呢？那样的话，我们就没法很好地用匹配为那个处理组观测值找到匹配了，不是吗？
大多数匹配方法会运用某种 “卡钳\(（caliper）\)” 或 “带宽\(（bandwidht）\)”，来判定一个匹配差异多大时就该舍弃。
基本思路很简单明了 —— 选一个数。这就是你的卡钳 / 带宽。要是距离（或者倾向得分的差异）的绝对值比这个数大，那就不把它算作匹配。要是这样一来，某个处理组观测值没有任何匹配了，那也没办法。那个观测值会被舍弃。
通常，卡钳 / 带宽是依据你要匹配的变量值的标准差（要是依据倾向得分匹配，就是倾向得分的标准差）来定义的，而非变量值本身，以此避免量纲问题。如前一节所述，在我们的信用卡债务示例中，处理组观测值 \(10305\) 与对照组观测值 \(27281\) 在 “四月账单\(（BillApril）\)”方面的距离为\(0.00035\)个标准差，与观测值 \(27719\) 的差异为\(0.00056\)个标准差。要是我们把卡钳定义为\(0.0005\)个标准差（要是你问我的话，这可能有点保守），那么那些匹配中只有一个会被算入，另一个会被舍弃。
有些匹配方法会自然地用到卡钳 / 带宽。任何基于核的匹配方法，都会对距离足够远的匹配赋予 \(0\) 权重，远到你所选的核函数会将其设为 \(0\)。对于 \(Epanechnikov核\)（埃潘尼科夫核 ），是距离为 \(1\) 或更大时。
另一种自然会这样做的匹配方法是精确匹配。在精确匹配中，卡钳 / 带宽设为 \(0\)！没有任何例外。要是你的匹配变量有丝毫差异，那就不算匹配。通常，精确匹配会以 “\(coarsened \space exact \space matching（粗化精确匹配,lacus、King和Porro \space 2012）\)” 的形式出现，它是精确匹配的一种，只是当你用连续变量作为匹配变量时，会先把它 “分组简化”。否则，你没法用连续变量作为匹配变量，因为别人不太可能，比如，和你有分毫不差的 “四月账单\(（BillApril）\)”金额，精确到分的百分位。\(coarsened\) 精确匹配实际上是一种很流行的方法，原因我会在 “多变量匹配” 部分讨论。
Stefano M Iacus, Gary King, and Giuseppe Porro. Causal inference without balance checking: Coarsened exact matching. Political Analysis, 20(1):1–24, 2012要是我们的方法不会贴心地把带宽设为\(0\)来帮我们选带宽，那我们该怎么选带宽呢？嗯，好问题。可惜，这里存在一种权衡。而且，你猜对了，是偏差和方差之间的权衡。
带宽设得越宽，我们能得到的潜在匹配就越多。要是我们允许多个匹配，就能基于更多观测值来计算处理组和对照组的均值，让均值更精确。要是我们只允许一个匹配，那么因没有匹配而被舍弃的观测值会更少，同样能让我们纳入更多观测值。
然而，当我们把带宽设得更宽时，会纳入更多差的匹配，这会让匹配质量变差，而且 “我们在阻断后门路径” 这一说法也更站不住脚。这就会把偏差重新引入到估计中。
关于”最优带宽选择”已有大量文献研究30，提出了诸多基于特定准则的带宽选取方法。常见准则如：“使样本对照组均值与实际对照组均值的平方误差最小化”。统计软件通常内置最优带宽选择程序，但建议查阅文档以明确其具体算法原理。
通常在进行任何匹配决策（包括带宽选择）时，选择往往介于以下两者之间：
匹配数量较少但质量较高，能产生偏差较小但精度较低的估计；
匹配数量较多但质量较差，能产生偏差较大但精度较高的估计。
以上就是核心要点。虽然还存在一些细节待探讨，但这些已是匹配估计中必须解决的关键概念性问题——或许还需考虑若干技术细节。至少这些内容能让您理解匹配分析的目标与原理，当然后续仍有更多技术细节需要展开。
14.4 多个匹配变量
呼——上一节内容真是信息量巨大！距离匹配 vs 倾向得分、匹配数量、偏差与方差的权衡…没完没了。而这还仅仅是单变量匹配的情况。现在我们要扩展到多变量匹配（这才是常规操作），复杂度想必会直线上升。
其实未必。实际上，从单变量匹配扩展到多变量匹配，核心问题只有一个：如何将这些多变量”压缩”成单一变量，以便直接套用前文的所有方法？
当然，这个问题还涉及诸多细节，但核心思路就是如此。我们面临多个变量，而匹配依赖于挑选”尽可能接近”的样本或根据接近程度加权——当”接近”只有一个定义时，操作会简单得多。因此，我们需要将多个变量的差异综合为单一度量。实现方式包括：距离匹配、倾向得分匹配或精确匹配，亦可混合使用——对某些变量要求精确匹配，对其他变量采用距离或倾向得分匹配31。每种选择都需审慎考量，这也正是本节要探讨的内容。
在本节及本章后续部分，我会使用来自布鲁克曼\(（2013）\)的数据添加示例代码。在该研究中，布鲁克曼想要探究美国政客的内在动机，即即便无明显回报，他们会做何事。这种动机的一个例子是，美国的黑人政客或许对支持美国黑人社区尤为感兴趣。这是一系列研究的一部分，这些研究旨在探究一般情况下，政客是否会向与他们类似的人提供额外支持。
David E Broockman.Black politicians are more intrinsically motivated to advance blacks’ interests: A field experiment manipulating political incentives.American Journal of Political Science, 57(3):521–536, 2013.为开展此项研究，布鲁克曼进行了一项实验。\(2010\)年，他向州立法议员（政客）发送了大量电子邮件32，仅为询问有关失业救济的信息。每封邮件均由虚构的“蒂龙·华盛顿”发出，在美国，该名字强烈暗示其归属于一名黑人男性 。
问题的核心在于议员是否会回复邮件。这如何揭示内在动机？因为布鲁克曼通过控制变量进行了巧妙设计：部分邮件声称发件人居住在该议员选区，部分则声明来自外地。回复选区外邮件对议员几乎没有直接利益——对方无法为你投票！但出于责任感或善意仍可能回复。
布鲁克曼进而提出三组问题：
黑人议员对选区外黑人邮件的回复率是否低于选区内的？是！
非黑人议员对选区外黑人邮件的回复率是否低于选区内的？也是！
关键发现：黑人议员的选区内外回复率差异是否小于非黑人议员？若成立，则表明黑人议员存在帮助黑人发件人的额外内在动机，同时验证了”内在动机假说”与”议员倾向帮助同类群体假说”33。
匹配是在哪个环节发挥作用的呢？是在最后一步，也就是我们对黑人和非黑人员工的选区内外回复差异进行比较的时候。这些群体（黑人和非黑人员工）往往在差异极大的地区当选。后门（混杂因素）大量存在。匹配可用于构建一个更优的对照组。在最初的研究中，布鲁克曼将选区的家庭收入中位数、选区人口中黑人的占比以及议员是否为民主党人用作匹配变量 34。
14.4.1 多个匹配变量下的距离匹配
让我们从距离度量开始：如何将多个匹配变量综合为单一距离指标，从而判定哪些观测”最接近”？这需要构建统一的距离计算公式。
终于可以松口气了：学界对马氏距离作为统一度量方法已形成基本共识35。
马氏距离的计算相当直观。我们先从稍简化的版本开始：首先将每个匹配变量的值除以其标准差，此时所有匹配变量的标准差都变为\(1\)。这能确保不会因变量量纲较大而导致其权重被放大。
在除以标准差后，我们可以计算距离。对于给定的处理组观测值\(A\)和给定的对照组观测值\(B\)，马氏距离是\(A\)和\(B\)之间所有差异的平方和。然后，在计算出总和后，再取平方根。换句话说，如果尝试用\(B\)的值来预测\(A\)的匹配变量，且每个匹配变量的标准差为\(1\)，那么它就是我们得到的平方残差和。这就是马氏距离。
或者说，无论如何，这是它的简化版本。我遗漏的一点是，我们不只是将每个变量除以其标准差。相反，我们要从变量的平方值中去除整个协方差矩阵的影响。如果所有匹配变量彼此不相关，这就简化为仅将每个变量除以其标准差。但如果它们彼此相关，这会消除匹配变量之间的关系，使它们在基于这些变量进行匹配之前相互独立。
这需要用一些矩阵代数来表示，但即便你不懂矩阵代数，你真正需要知道的是，\(x^Tx\)可以理解为 “将所有\(x\)平方，然后相加” 。然后，如果你眯着眼并记住取逆（\(^{-1}\)）意味着 “除以” ，你就能弄清楚这里发生了什么。两个观测值\(x_1\)和\(x_2\)之间的马氏距离为：
其中，\(x_1\) 是观测\(1\)所有匹配变量构成的向量，\(x_2\) 同理，\(S\)是所有匹配变量的协方差矩阵。
为何消除不同变量间的关系是有益的呢？因为这样能避免我们的匹配过度依赖某种潜在特征 —— 这种特征可能恰好在多个变量中体现。比如，假设你依据胡须长度、臂毛浓密程度和 “男子气概指数” 得分进行匹配。这三者在不同程度上都与 “身为男性” 相关。若不消除协方差，你实际上相当于三次依据 “身为男性” 进行匹配。这会造成冗余，而且匹配过程可能会拒绝将男性与非男性进行匹配，即便他们在其他特征上很匹配。若消除协方差，我们仍会依据 “男性特质” 匹配，但不会重复计算。
再回到我们的信用卡债务示例，我们将匹配变量扩展为两个：四月份账单\(（BillApril）\)以及个人的年龄（以年为单位，\(Age\)）。对比处理组观测值\(10305\) 和未处理组观测值\(27281\)，我们已知它们在\(BillApril\)上的距离为 \(|91630 - 91609| = 21\) 。年龄的差异为 \(|23 - 37| = 14\) 。
首先，我们用简化方法计算马氏距离，即忽略匹配变量间的关系。\(BillApril\)的标准差是新台币\(59554\)，将\(BillApril\)的所有值除以该标准差后，新的差异为 \(|1.5386 - 1.5382| = 0.0004\) 个标准差，变化不大！然后，\(Age\)的标准差是 \(9.22\) ，最终得到的距离为 \(|2.49 - 4.01| = 1.52\) 个标准差。
我们得到了以标准差为单位的差异。然后将它们平方后相加，即 \(0.0004^2 + 1.52^2 = 2.310\) ，再取平方根 ——\(\sqrt{2.310} = 1.520\) —— 即得到马氏距离为\(1.520\) 。随后，我们将该马氏距离与观测值\(10305\)和其他所有未处理组观测值之间的距离进行比较，最终与能使距离最小的未处理组观测值完成匹配。
t(x) %*% solve(cv) * x## BillApril AGE
## [1,] 2.568173e-05 2.31197# calculate the Mahalanobis distance
sqrt(t(x) %*% solve(cv) %*% x)## [,1]
## [1,] 1.520525让我们重新计算，但这次规范地加入协方差矩阵的运算步骤：已知\(BillApril\)差异值为\(21\)，年龄差异值为\(-14\)，将其与协方差矩阵结合运算后，得到马氏距离为
好的，在这个特定案例中，与简化版本的结果差异不大（\(1.521\) 对比 \(1.520\) ）。但这仍然绝对值得一做。毕竟，在我们用两种方法计算之前，我们并不知道这样做会不会有影响。
马氏距离也有一种便捷的图形化阐释。我们对数值进行平方、相加，然后取平方根。这有没有让你想起什么？嗯，对于两个匹配变量而言，这就如同我们从毕达哥拉斯定理（勾股定理）中所知的，是两点之间的距离。马氏距离是一个三角形斜边的长度，其中一条直角边为第一个匹配变量的差值，另一条直角边为第二个匹配变量的差值 。
您可以在\(图14.6\)中看到距离计算——它就是标准化坐标轴后两点之间的直线距离。这个概念同样适用于多维情况：增加第\(3\)个匹配变量时，您仍然是在两点之间画线，只不过这次是在三维空间中。观察周围任意两个物体并想象连接它们的线条，若匹配变量是宽度、长度和高度这三个空间维度，则该连线长度即为它们的马氏距离。
虽然四维及以上空间难以可视化，但计算原理相同——可自由扩展至任意维度。
sub2 <- uci %>%
mutate(id = row_number()) %>%
mutate(BA_std = (BillApril-mean(BillApril))/sd(BillApril),
Age_std = (AGE - mean(AGE))/sd(AGE)) %>%
mutate(blank = '')
# 筛选最相似样本
sub2 <- sub2 %>%
# 按BA_std与id=10305样本的绝对距离排序
arrange(abs(BA_std - sub2$BA_std[10305])) %>%
# 取前200个最相似的样本
slice(1:200) %>%
# 在200个样本中，再按Age_std与id=10305样本的绝对距离排序
arrange(abs(Age_std - sub2$Age_std[10305])) %>%
# 取前100个最相似的样本
slice(1:100) %>%
# 特殊处理：将id=11874和21733的样本优先排列（TRUE=1会排在前面）
arrange(-(id %in% c(11874, 21733))) %>%
# 最终选取前20个样本
slice(1:20)
# 提取两个特定样本(id=21733和11874)用于对比
pts <- sub2 %>% dplyr::filter(id %in% c(21733, 11874))
# 计算两个样本的标准化值与第一个样本的差值
pts %>%
dplyr::select(BA_std, Age_std) %>%
reframe(
BA = BA_std - first(BA_std), # BillApril标准化值的差值
AS = Age_std - first(Age_std) # AGE标准化值的差值
)## # A tibble: 2 × 2
## BA AS
## <dbl> <dbl>
## 1 0 0
## 2 -0.0178 0.325# 使用ggplot创建散点图，x轴为标准化后的4月账单(BA_std)，y轴为标准化后的年龄(Age_std)
ggplot(sub2) +
# 添加散点图层
geom_point(
aes(x = BA_std,
y = Age_std,
shape = LateApril, # 形状映射LateApril变量
size = id %in% c(21733, 11874)
)
) +
# 高亮显示特定ID(21733和11874)的点
gghighlight::gghighlight(
id %in% c(21733, 11874),
label_key = blank,
use_group_by = FALSE
) +
# 自定义形状样式和标签
scale_shape_manual(values = c(1,19), labels = c('Not Late\nin April','Late in\nApril')) +
# 自定义点的大小
scale_size_manual(values = c(1, 3)) +
# 添加三条线段注释(构成一个直角三角形)
annotate(geom = 'segment', x = pts$BA_std[2], xend = pts$BA_std[2], y = pts$Age_std[1], yend = pts$Age_std[2]) + # 垂直边
annotate(geom = 'segment', x = pts$BA_std[1], xend = pts$BA_std[2], y = pts$Age_std[1], yend = pts$Age_std[1]) + # 水平边
annotate(geom = 'segment', x = pts$BA_std[1], xend = pts$BA_std[2], y = pts$Age_std[1], yend = pts$Age_std[2]) + # 斜边
# 添加三个文本注释
annotate(geom = 'text', family = 'sans', size = 10/.pt, angle = 90,x = pts$BA_std[2], y = mean(pts$Age_std), label = 'Age Difference', vjust = -.4) + # 垂直边标签
annotate(geom = 'text', family = 'sans', size = 10/.pt, x = mean(pts$BA_std), y = pts$Age_std[1], label = 'April Bill Difference', vjust = 1.4) + # 水平边标签
annotate(geom = 'text', family = 'sans', size = 10/.pt, angle = -45, x = pts$BA_std[2], y = pts$Age_std[2], label = 'Mahalanobis', vjust = -.35, hjust = -1) + # 斜边标签
expand_limits(x = .91, y = -.6) + # 扩展坐标轴范围
# 设置坐标轴标签
labs(x = 'April Bill, Standardized',y = 'Age, Standardized') +
# 使用pubr主题并自定义主题元素
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
legend.position = c(.9, .2), # 图例位置
legend.background = element_rect(color = 'black')) # 图例背景## Warning in ggrepel::geom_label_repel(mapping, layer$data, fill = "white", :
## Ignoring unknown aesthetics: shape
# 保存图形为PDF
ggsave('Matching/Matching/mahalanobis.pdf', width = 6, height = 5, device=cairo_pdf)这确实揭示了马氏匹配以及一般意义上的距离匹配的一个缺点 —— 维度灾难。你添加的匹配变量越多，事物之间似乎就 “越远”。可以这样理解：想象一张办公楼的俯视图，这个俯视图只有两个维度：南北和东西。位于大楼西南角的旺达和同样在西南角的泰勒离得很近，是绝佳的匹配！现在将镜头移到侧面，这样我们就能在三维空间中看到同一栋楼：南北、东西和上下。结果发现，旺达在一楼的西南角，而泰勒在\(22\)楼的西南角，不再是理想的匹配了。
维度灾难意味着，你添加的匹配变量越多，为任何给定的处理组观测值找到可接受匹配的可能性就越小。有几种方法可以解决这个问题。一种方法是尝试限制匹配变量的数量，尽管这意味着我们无法关闭那些后门路径 —— 这可不太好。另一种方法是获取大量观测值。如果你有极多的观测值可供筛选来找到一个匹配，那么匹配难以找到也没关系。这听起来不错，但问题是，我们能获得的观测值数量是有限的。
第三种方法是，随着维度数量增加，对匹配质量、至少对卡尺 / 带宽\(（caliper/bandwidth）\)的值，适度放宽要求。有时，人们会不断尝试不同带宽，直到匹配数量看起来合适。一种更规范的方法是，使用常规带宽，但将马氏距离得分除以匹配变量数量的平方根\(（Mahony \space 2014）\)。
Colin Mahony. Effects of dimensionality on distance and probability density in climate space. The Seasons Alter, 2014.让我们将布鲁克曼数据与马氏距离方法结合！我们将使用\(1\space vs \space 1\)最近邻匹配，基于\(3\)个匹配变量\(（medianhhincme, blackpercent, lag\_democrat）\)计算马氏距离，将黑人议员选区与非黑人议员选区进行匹配。代码包含两部分：匹配权重计算和结果变量的处理效应估计。在布鲁克曼研究中，结果变量\(responded\)表示议员是否回复邮件。实际研究设计比单纯分析\(leg\_black\)对\(responded\)的影响更为复杂，但后续会详述。当前仅展示在\(medianhhincme, blackpercent\)和\(lag\_democrat\)上匹配后\(leg\_black\)与\(responded\)的关系36。
library(Matching)
library(tidyverse)
br <- causaldata::black_politicians
Y <- br %>% pull(responded)
D <- br %>% pull(leg_black)
X <- br %>% dplyr::select(medianhhincom , blackpercent , leg_democrat) %>% as.matrix()
# Weight = 2, oddly , denotes Mahalanobis distance
M <- Match(Y, D, X, Weight = 2, caliper = 1)
# See treatment effect estimate
summary(M)##
## Estimate... -0.0073462
## AI SE...... 0.072683
## T-stat..... -0.10107
## p.val...... 0.91949
##
## Original number of observations.............. 5593
## Original number of treated obs............... 364
## Matched number of observations............... 363
## Matched number of observations (unweighted). 405
##
## Caliper (SDs)........................................ 1 1 1
## Number of obs dropped by 'exact' or 'caliper' 1# obtain the matched data for subsequent analysis (Note: This method will duplicate records of observations for each match)
matched_treated <- tibble(id = M$index.treated , weight = M$weights)
matched_control <- tibble(id = M$index.control , weight = M$weights)
matched_sets <- bind_rows(matched_treated , matched_control)
# sum the weights for each ID
matched_sets <- matched_sets %>% group_by(id) %>% summarize(weight = sum(weight))
matched_br <- br %>% mutate(id = row_number()) %>% left_join(matched_sets, by = 'id')
# Note: The calculation of standard errors here is incorrect.
lm(responded ~ leg_black , data = matched_br, weights = weight)##
## Call:
## lm(formula = responded ~ leg_black, data = matched_br, weights = weight)
##
## Coefficients:
## (Intercept) leg_black
## 0.398531 -0.007346import pandas as pd
import numpy as np
# The more-popular matching tools in sklearn
# are more geared towards machine learning than statistical inference
from causalinference.causal import CausalModel
from causaldata import black_politicians
br = black_politicians.load_pandas().data
# Get our outcome, treatment, and matching variables
# We need these as numpy arrays
Y = br['responded'].to_numpy()
D = br['leg_black'].to_numpy()
X = br[['medianhhincom', 'blackpercent', 'leg_democrat']].to_numpy()
# Set up our model
M = CausalModel(Y, D, X)
# Fit, using Mahalanobis distance
M.est_via_matching(weights = 'maha', matches = 1)
print(M.estimates)##
## Treatment Effect Estimates: Matching
##
## Est. S.e. z P>|z| [95% Conf. int.]
## --------------------------------------------------------------------------------
## ATE -0.050 0.175 -0.285 0.776 -0.392 0.293
## ATC -0.052 0.186 -0.282 0.778 -0.418 0.313
## ATT -0.012 0.099 -0.120 0.904 -0.206 0.182# Note it automatically calcultes average treatments on average, on treated, and on untreated/control (ATC)causaldata black_politicians.dta, use clear download
* Create an id variable based on row number,
* which will be used to locate match IDs
g id = _n
* teffects nnmatch does nearest-neighbor matching using first the outcome
* variable and the matching variables then the treatment variable.
* The generate option creates variables with the match identities
teffects nnmatch (responded medianhhincom blackpercent leg_democrat) (leg_black), nneighbor(1) metric(mahalanobis) gen(match)
* Note this will produce a treatment effect by itself
* Doing anything else with the match in Stata is pretty tricky
* (or using a caliper). If you want either of those you may want to
* look into the mahapick or psmatch2 packages
* Although these won't give the most accurate standard errors或许马氏距离还不够彻底。我们可以摒弃将一堆距离归为一个的做法，也不用再纠结某个变量是否比其他变量更重要，而是假定所有变量都极其重要，不存在任何不匹配情况。我们又回到了精确匹配，具体来说是粗化精确匹配\(（coarsened \space exact \space matching）\)我在匹配概念那部分首次提到过\(lacus\)、\(King\)和 \(Porro\) \(2012\)年的研究。
Stefano M Iacus, Gary King, and Giuseppe Porro. Causal inference without balance checking: Coarsened exact matching. Political Analysis, 20(1):1–24, 2012.在粗化精确匹配中，只有在每个匹配变量上都完全匹配，才算作匹配成功。“粗化”这一环节的存在，是因为如果有连续变量需要匹配，你得先把它们 “粗化”，也就是分组，而不是依据精确数值匹配，否则根本无法找到匹配。这种粗化是怎样操作的呢？以信用卡债务示例里的 “四月账单” 为例，我们可以把这个变量分成\(10\)个十分位数组。对于第 \(10305\) 行那个常用的处理组观测值，不再匹配该行 “新台币 \(91630\) 元” 这个精确数值（ 没有其他观测值是这个数值 ），而是确认新台币 \(91630\) 元处于 “\(BillApril\)” 变量的第 \(80\) 到 \(90\) 百分位数区间，然后把它归到 “介于新台币 \(63153\)元（ 第 \(80\) 百分位数 ）和新台币\(112110\)元（ 第 \(90\) 百分位数 ）之间” 这一类别（ 所有观测值中约 \(10\%\) 属于该类别 ）37。
# check the q distribute
quantile(uci$BillApril, (0:10)/10)## 0% 10% 20% 30% 40% 50% 60% 70%
## -339.6030 0.0000 0.4760 2.7017 8.7702 17.0710 25.5084 39.2522
## 80% 90% 100%
## 63.1506 112.1104 961.6640一旦你确定了如何对每个连续变量进行分组，就去寻找精确匹配。每个处理组观测值，只有在找到至少一个精确匹配时才会被保留，否则就舍弃。每个对照组观测值，同样只有在找到至少一个精确匹配时才会被保留，并且会被赋予一个权重，该权重等于它与处理组观测值的匹配数量，除以与该处理组观测值匹配的对照组观测值数量。然后，再将这个结果乘以匹配的对照组观测值总数除以匹配的处理组观测值总数。换句话说，就是\((MyTreatedMatches / MyControlMatches) \times (TotalControlMatches / TotalTreatedMatches)\) 。这个过程无疑确保了，处理组和对照组在匹配变量（分组后 ）上绝对没有差异38。
粗化精确匹配确实需要强大的计算能力才能实施，尤其是在维度灾难出现的情况下。\(Z.Zhao（2004年）\)研究了一个使用相当标准数据的例子，尝试对\(12\)个变量进行匹配。在将这\(12\)个变量全部交互后，产生了大约\(600\)万个单元。而这个数据集的观测值数量远少于\(600\)万。不是每个观测值都能找到匹配。
Zhong Zhao. Using matching to estimate treatment effects: Data requirements, matching metrics, and Monte Carlo evidence. Review of Economics and Statistics, 86(1):91–107, 2004.对非常大的样本量的需求也不是儿戏，不能敷衍了事。粗化精确匹配如果应用于中等规模的样本（或者匹配变量过多的任何规模样本 ），会导致大量处理组观测值被舍弃。\(Black、lalkiya和Lerner（2020）\)重复了五项粗化精确匹配研究，发现很多处理组观测值最终都被舍弃了，这使得处理效应估计值的噪声大得多，而且如果某些类型的处理组观测值比其他类型更易找到匹配，还可能导致结果无法很好地反映平均处理效应。其中一项研究是布鲁克曼\(（2013）\)的 —— 哎呀！如果你一直好奇为什么布鲁克曼代码示例中的匹配变量列表这么短，是因为即便这个列表，在粗化精确匹配时也已经导致不少观测值被舍弃了，再增加变量会让问题更严重。\(Black、lalkiya和Lerner（2020）\)建议完全摒弃这种方法，不过我不会走那么极端。我认为你大概应该把它的使用限制在超大数据集上，而且你必须检查有多少处理组观测值因为找不到匹配而丢失。如果丢失很多，就该换用其他方法。
对大量观测值的需求，以及计算相对简便、计算要求较低39，使得粗化精确匹配在大数据数据科学领域颇受欢迎。在该领域，你可能会听到精确匹配被称作 “\(Doppelgangers\)（可理解为’\(Doppelgangers\)‘，实际是德语外来词，常指’极为相似的人’ ）”40。它在超大规模行政数据集中也很常用，比如研究人员设法获取到一国全体人口的政府记录这类数据时41。
咱们来看 \(Klemick、Mason和Sullivan（2020）\)的研究，这是一个粗化精确匹配在既非大数据也非超大规模行政数据集场景中应用的例子。他们想探究美国环境保护署\(（EPA）\)的超级基金清理项目对儿童健康的影响。该项目中，\(EPA\)先确定遭受严重环境污染的区域，然后对其进行清理。清理过程是否改善了儿童健康呢？
Heather Klemick, Henry Mason, and Karen Sullivan. Superfund cleanups and children’s lead exposure. Journal of Environmental Economics and Management,100:1022–1089, 2020他们的基本研究设计是对比一个社区在清理前后儿童血液中的铅含量，预期不仅清理后铅含量会下降，而且在距离清理场地非常近（不到\(2\)公里）区域的铅含量改善程度，会比距离稍远（\(2-5\)公里）区域的更大，因为后者在污染存在时受影响可能没那么大42。
一如往常，距离超级基金场地的远近和包括血铅水平在内的总体健康指标之间存在后门路径。作者在社区层面使用粗化精确匹配来关闭部分后门，依据社区靠近的超级基金场地、区域内住房的年代（\(1940\)年前建成住房的占比）、接受福利救助人口的占比、非裔美国人的占比，以及接受用于测量血铅水平（作为研究因变量）的血液筛查人口的占比，对社区进行匹配。
他们最初的数据涵盖约\(38\)万名居住在超级基金场地\(2\)公里内的人口，以及\(90\)万名居住在场地\(2-5\)公里外的人口。经过匹配，大量观测值被舍弃。匹配后的样本中，居住在超级基金场地\(2\)公里内的约有\(20.1\)万人，居住在\(2-5\)公里外的约有\(35.3\)万人 。
样本数量的大幅减少或许令人担忧，但我们确实看到，居住在场地附近和稍远区域人群之间的差异显著降低。匹配变量看起来当然更相似了，而且那些可能处于后门路径上的其他变量也是如此。西班牙裔人口占比并非匹配变量，但匹配前差异的一半在匹配后消失了。交通密度和教育水平方面的差异也有类似的降低。
他们有什么发现呢？该项目似乎有效降低了儿童的血铅水平。根据不同的设定，在不进行任何匹配时，他们发现铅水平降低了\(5.1\%-5.5\%\)。而进行匹配后，结果数值更大，为\(7.1\%-8.3\%\) 。
我们自己该如何进行粗化精确匹配呢？咱们回到布鲁克曼\(（2013）\)的研究，他在最初的论文中运用了粗化精确匹配 。
library(cem)## Loading required package: tcltk## Loading required package: lattice##
## How to use CEM? Type vignette("cem")library(tidyverse)
br <- causaldata::black_politicians
brcem <- br %>%
dplyr::select(responded, leg_black, medianhhincom, blackpercent, leg_democrat) %>%
na.omit() %>%
# Must be a data.frame, not a tibble
as.data.frame()
inc_bins <- quantile(brcem$medianhhincom, (0:6)/6)
create_even_breaks <- function(x, n) {
minx <- min(x)
maxx <- max(x)
return(minx + ((0:n)/n)*(maxx-minx))
}
bp_bins <- create_even_breaks(brcem$blackpercent, 6)
# A binary variable only needs to be divided into 2 bins.
ld_bins <- create_even_breaks(brcem$leg_democrat,2)
allbreaks <- list('medianhhincom' = inc_bins,'blackpercent' = bp_bins,'leg_democrat' = ld_bins)
# ensure that matching is not based on the outcome. Note that the baseline.group is the treated group.
c <- cem(treatment='leg_black', data=brcem, baseline.group='1', drop='responded', cutpoints=allbreaks, keep.all=TRUE)##
## Using 'leg_black'='1' as baseline group# obtain the matching weights for other analytical purposes.
brcem <- brcem %>% mutate(cem_weight = c$w)
lm(responded ~ leg_black, data = brcem, weights = cem_weight)##
## Call:
## lm(formula = responded ~ leg_black, data = brcem, weights = cem_weight)
##
## Coefficients:
## (Intercept) leg_black
## 0.34680 0.02302# Or use their (matching packages’) built - in effect estimation function att(). Note: These functions offer various options, such as using logit regression or machine learning methods to estimate treatment effects.
att(c, responded ~ leg_black, data = brcem)##
## G0 G1
## All 5229 364
## Matched 4491 338
## Unmatched 738 26
##
## Linear regression model on CEM matched data:
##
## SATT point estimate: 0.023020 (p.value=0.391783)
## 95% conf. interval: [-0.029659, 0.075699]14.4.2 熵平衡
距离匹配领域出现了一个新方法，虽然在实际应用中，相较于其他方法，它没那么流行，但似乎具备一些非常不错的特性。这个新方法就是熵平衡\(（entropy \space balancing）\)，源自 \(Hainmueller（2012）\)的研究。
Jens Hainmueller. Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies. Political Analysis, 20(1):25–46, 2012.在其他匹配方法中，你会以某种方式把匹配变量整合起来，然后基于这个整合结果进行匹配。接着，你期望所做的匹配能消除处理组和对照组之间的所有差异。而且，正如我会在 “检查匹配质量” 部分讨论的，你要查看是否还有差异残留，然后满心期待这些差异都消失了。
一旦你确定了一组限制条件，熵平衡就会开始运作，寻找一组满足这些限制条件的权重。只要这些限制条件有可能满足，你也会得到这样的权重44这差不多就是全部内容了！熵平衡能带来和粗化精确匹配一样的 “确保无差异” 结果，但无需限制匹配变量的数量，也不用舍弃大量处理组观测值。那我们之前怎么没这么做呢？
这部分内容简短有两个原因：第一，熵平衡的直观解释相当简单，但要深入探讨任何细节，都会变得复杂得多；第二，在撰写本文时，熵平衡不像其他方法那样广为人知。或许是因为相较于其他方法，它需要更多的技术细节45。
但它仍值得学习。也许未来它会变得更流行、更受关注。也许你会使用它，而这会让它在未来更流行、更受关注。或者，至少它可以作为众多未被纳入本书的其他匹配方法的替代。这类方法还有一些呢。怎么，我们现有的匹配程序的八百万种选项组合还不够吗？不！还有遗传匹配、子类化、最优匹配等等。
咱们来编写一个熵平衡估计的代码。示例会使用 \(R\) 和 \(Stata\)。有一个可用的 \(Python\) 库 “\(empirica\_calibration\)”，其中包含熵平衡功能，但你得费些周折，因为就连安装它都不像其他软件包那么简单。我目前就略过\(python\)的例子了，但你可以自己去研究46。
library(ebal)## ##
## ## ebal Package: Implements Entropy Balancing.## ## See http://www.stanford.edu/~jhain/ for additional information.library(tidyverse)
library(modelsummary)
br <- causaldata::black_politicians
Y <- br %>% pull(responded)
D <- br %>% pull(leg_black)
X <- br %>%
dplyr::select(medianhhincom, blackpercent, leg_democrat) %>%
mutate(incsq = medianhhincom^2, bpsq = blackpercent^2) %>%
as.matrix()
eb <- ebalance(D, X)## Converged within tolerancebr_treat <- br %>% dplyr::filter(leg_black == 1) %>% mutate(weights = 1)
br_con <- br %>% dplyr::filter(leg_black == 0) %>% mutate(weights = eb$w)
br <- bind_rows(br_treat, br_con)
m <- lm(responded ~ leg_black, data = br, weights = weights)
msummary(m, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 0.302*** |
| (0.009) | |
| leg_black | 0.091*** |
| (0.013) | |
| Num.Obs. | 5593 |
| R2 | 0.009 |
| R2 Adj. | 0.009 |
| AIC | 24483.2 |
| BIC | 24503.1 |
| Log.Lik. | -12238.623 |
| F | 51.253 |
| RMSE | 0.51 |
14.5 基于多匹配变量的倾向得分加权法
倾向得分\(（propensity \space score）\)可能是将多个匹配变量整合为一个可用于匹配的单一数值的最常用方法。记住，倾向得分是对给定观测值接受处理的概率的估计值。倾向得分匹配通常意味着选择一组倾向得分值相近的匹配对照观测值。不过，如前所述，现代建议倾向于使用倾向得分进行逆概率加权\(（King和Nielsen,2019）\)。
King, Gary, and Richard Nielsen. 2019. “Why Propensity Scores Should Not Be Used for Matching.” Political Analysis 27 (4): 435–54.在本节中，我将重点关注倾向得分本身的估计。然后，正如本节标题所示，你们将继续使用这些倾向得分来构建匹配权重，而非选择一组匹配观测值。关于如何从倾向得分构建这些权重的具体细节将在代码示例中展示，也会在本章末尾的 “匹配与处理效应” 部分进一步详述，因为权重的计算取决于你想要的处理效应类型 。
逆概率权重\(（inverse \space probability \space weights）\)：基于倾向得分的匹配权重，通常为\(1\)除以实际接受的处理水平的接受倾向。
倾向得分通常通过回归来估计。这对我来说非常方便，因为我刚写了一整章关于回归的内容，你们可以直接去读那一章47。
具体而言，倾向得分通常通过\(Logit\)回归或\(Probit\)回归进行估计48。基于回归的倾向得分估计，其核心工作在于构建恰当的回归模型。
在构建回归模型时，有哪些选择要做呢？你会遇到标准的模型构建问题 —— 任何影响处理的决定因素肯定都应该包含在模型里。但也有重要的函数形式方面的考虑。哪些匹配变量应该以多项式形式纳入？或者以交互项形式纳入？
你在这里是要预测某个东西，也就是倾向得分，但你得把思考放在因果图的范畴内。毕竟，你是在尝试关闭后门路径。所以在构建模型时，若变量处于后门路径上，就将其作为预测变量纳入；若处于前门路径上，就将其排除，这仍然是关键。
不过，在选择模型时，你也不是全然盲目（或者完全依赖理论 ）。毕竟，一个好的倾向得分应该能起到关闭后门路径的作用，对吧？我们可以检验这一点。\(图14.3\)中的后门路径形式是\(Treatment \leftarrow TreatmentPropensity \leftarrow A \rightarrow Outcome\)。控制\(TreatmentPropensity\)应该能关闭\(Treatment\)和匹配变量\(A\)之间的所有路径。所以，正如 \(Dehejia和Wahba（2002）\)所建议的，你可以把倾向得分划分成区间。然后，在每个区间内（也就是控制 \(treatmentPropensity\)后 ），你可以检验每个匹配变量是否与处理相关。如果相关，你可能得尝试为有问题的匹配变量添加更多多项式项或交互项。这被称为 “分层检验”，是一种平衡性检验，我会在下面的 “平衡性检验” 部分详细介绍。
Dehejia, Rajeev H., and Sadek Wahba. 2002. “Propensity Score-Matching Methods for Nonexperimental Causal Studies.” Review of Economics and Statistics 84 (1): 151–61.一旦完成\(Logit\)或\(Probit\)模型估计后，您将获得预测概率49——这些预测概率即为所需的倾向得分。
library(causalweight)## Loading required package: rangerlibrary(tidyverse)
br <- causaldata::black_politicians
# ps est
m <- glm(leg_black ~ medianhhincom + blackpercent + leg_democrat, data = br, family = binomial(link = 'logit'))
# ps value
br <- br %>% mutate(ps = predict(m, type = 'response'))
# Trim control group samples that fall outside the range of the treatment group's propensity scores (see the discussion on "common support" later in the text).
minps <- br %>% dplyr::filter(leg_black == 1) %>% pull(ps) %>% min(na.rm = TRUE)
maxps <- br %>% dplyr::filter(leg_black == 1) %>% pull(ps) %>% max(na.rm = TRUE)
br <- br %>% dplyr::filter(ps >= minps & ps <= maxps)
# create inverse probability weighting weights
br <- br %>% mutate(ipw = case_when(leg_black == 1 ~ 1/ps, leg_black == 0 ~ 1/(1-ps)))
# For weighted regression (the standard errors here will be wrong unless we bootstrap the entire process).
# See the section on doubly robust estimation or the code examples in the simulation chapter.
lm(responded ~ leg_black, data = br, weights = ipw)##
## Call:
## lm(formula = responded ~ leg_black, data = br, weights = ipw)
##
## Coefficients:
## (Intercept) leg_black
## 0.4035 0.1564# Or we can use thecausalweight package! First, extract the variables: outcome variable (Outcome)/treatment variable (Treatment)/covariates (Covariates)
Y <- br %>% pull(responded)
# Treatment
D <- br %>% pull(leg_black)
X <- br %>% dplyr::select(medianhhincom, blackpercent, leg_democrat) %>% as.matrix()
# Note: The default output is the average treatment effect (ATE), not the average treatment effect on the treated (ATT)
# and its trimming of propensity scores is based on extreme values rather than the range of scores matching the treatment group
IPW <- treatweight(Y, D, X, trim = .001, logit = TRUE)
IPW$effect## [1] 0.1426612IPW$se## [1] 0.1394989import pandas as pd
import numpy as np
import statsmodels.formula.api as sm
# The more-popular matching tools in sklearn are more geared towards machine learning than statistical inference
from causalinference.causal import CausalModel
from causaldata import black_politicians
br = black_politicians.load_pandas().data
# Get our outcome, treatment, and matching variables as numpy arrays
Y = br['responded'].to_numpy()
D = br['leg_black'].to_numpy()
X = br[['medianhhincom', 'blackpercent', 'leg_democrat']].to_numpy()
# Set up our model
M = CausalModel(Y, D, X)
# Estimate the propensity score using logit
M.est_propensity()
# Trim the score with improved algorithm trim_s to improve balance
M.trim_s()
# If we want to use the scores elsewhere, export them (we could have also done this with sm.Logit)
br['ps'] = M.propensity['fitted']
# We can estimate the effect directly (note this uses "doubly robust" methods
# as will be later described, which is why it doesn't match the sm.wls result)
M.est_via_weighting()
print(M.estimates)##
## Treatment Effect Estimates: Weighting
##
## Est. S.e. z P>|z| [95% Conf. int.]
## --------------------------------------------------------------------------------
## ATE 0.047 0.080 0.588 0.557 -0.109 0.203# Or we can do our own weighting
br = br.assign(ipw = lambda x: x.leg_black*(1/x.ps) + (1-x.leg_black)*(1/(1-x.ps)))
# Now, use the weights to estimate the effect (this will produce incorrect standard errors unless we bootstrap the whole process, as in the doubly robust section later, or the Simulation chapter)
m = sm.wls(formula = 'responded ~ leg_black', weights = br['ipw'], data = br).fit()
m.summary()| Dep. Variable: | responded | R-squared: | 0.022 |
|---|---|---|---|
| Model: | WLS | Adj. R-squared: | 0.022 |
| Method: | Least Squares | F-statistic: | 128.6 |
| Date: | 周日, 17 8月 2025 | Prob (F-statistic): | 1.75e-29 |
| Time: | 22:07:52 | Log-Likelihood: | -5768.6 |
| No. Observations: | 5593 | AIC: | 1.154e+04 |
| Df Residuals: | 5591 | BIC: | 1.155e+04 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 0.4100 | 0.009 | 43.596 | 0.000 | 0.392 | 0.428 |
| leg_black | 0.1499 | 0.013 | 11.340 | 0.000 | 0.124 | 0.176 |
| Omnibus: | 2184.371 | Durbin-Watson: | 1.936 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 2062976.257 |
| Skew: | 0.078 | Prob(JB): | 0.00 |
| Kurtosis: | 97.087 | Cond. No. | 2.63 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
causaldata black_politicians.dta, use clear download
* We can estimate our own propensity score with probit or logit
logit leg_black medianhhincom blackpercent leg_democrat
* Get predicted values for our propensity score
predict ps, pr
* "Trim" control observations outside of
* treated propensity score range
* (we'll discuss this later in Common Support)
summ ps if leg_black
replace ps = . if ps < r(min) | ps > r(max)
* Get inverse probability weights
g ipw = 1/ps if leg_black
replace ipw = 1/(1-ps) if !leg_black
* Use inverse weights in regression
reg responded leg_black [pw = ipw]
* (or see the bootstrap standard error example in the section on doubly
* robust estimation, or the simulation chapter) Note that simple vce(bootstrap)
* won't work since it won't bootstrap the matching process
* Or do everything in teffects
* (This will use improved standard errors,
* unless you've done bootstrap,
* but won't let us do the trim)
teffects ipw (responded) (leg_black medianhhincom blackpercent leg_democrat, logit)还有其他估算倾向得分的方法。实际上，任何能将二元变量对多个预测变量进行回归，且得出 \(0\) 到 \(1\) 之间预测值的方法，理论上都可行。但这并不意味着所有这些方法都是好主意。
在众多替代方法中，有些其他方法被证明特别有用。总体而言，除了\(logit\)和\(probit\)之外，那些有用的方法，做到了\(logit\)和\(probit\)做不到的事50。\(logit\)和\(probit\)难以应对的两类情况是高度非线性和高维度。
高度非线性是个问题？\(logit\)和\(probit\)难道不是非线性回归模型吗？嗯，确实是，但如果你回顾第\(13\)章的 “非线性回归” 部分，就会发现它们要求指数函数具有线性关系。当然，你可以添加多项式项，但该加多少呢？加够了吗？加太多的话，情况就会开始变得怪异。
当存在大量匹配变量时，就会出现高维度问题。或者更准确地说，当回归模型中所有匹配变量之间存在大量交互作用时，这种情况尤为显著。理论上，回归分析确实可以处理这类问题，但随着变量项数不断增加，回归方法终将不再是解决问题的最佳工具。
那么，究竟什么才是解决这一问题的理想工具？现有多种备选方案，其中在医疗领域较为知名的是高维倾向得分估计法\(（high-dimmensional \space propensity \space score \space estimation）\)。该方法的核心在于：在进行\(logit\)回归前，先从大量候选变量中筛选出兼具理论重要性和高预测效力的”最优”协变量组合。
Sebastian Schneeweiss,Jeremy A Rassen, Robert J Glynn, Jerry Avorn, Helen Mogun, and M Alan Brookhart. Highdimensional propensity score adjustment in studies of treatment effects using health care claims data. Epidemiology, 20(4):512,2009稍远但成效显著的是机器学习方法。机器学习特别擅长处理高维数据的预测问题——而这恰恰符合我们利用多维变量预测处理效应的需求。
运用机器学习方法估计倾向得分时，有两种主流方法表现突出：其一是正则化回归（如第\(13\)章所述），该方法通过纳入所有变量并利用惩罚参数自动决定变量的降权或剔除；其二是提升回归\(（boosted \space regression）\)。
提升回归是一种迭代算法，其工作原理如下：首先构建一个二元回归模型（如\(logit\)模型），识别预测误差较大的观测样本；随后在后续迭代中加大这些预测不佳样本的权重，迫使模型重点关注并降低其预测误差。该过程循环执行多次，最终通过整合所有迭代生成的子模型，输出一个综合评分。当应用于倾向得分估计时，这种方法能获得相当精准的估计结果\(（Griffn \space et \space al., 2017）\)。
Beth Ann Griffn, Daniel F McCaffrey, Daniel Almirall, Lane F Burgette, and Claude Messan Setodji. Chasing balance and other recommendations for improving nonparametric propensity score models. Journal of Causal Inference,5(2), 2017.14.6 匹配的假设
14.6.1 条件独立性假设
与其他统计方法一样，无论我们的方法设计得多么周密，它们始终依赖于某些假设。采用匹配方法可以放宽部分假设（例如线性假设——除非我们使用回归来估计倾向得分），但这绝不意味着匹配方法无需任何假设。
某些假设延续了我们先前的研究框架，其中最关键的是条件独立假设\(（conditional \space independence \space assumption）\)。还记得回归章节中关于”遗漏变量偏差”的讨论吗？——它实质上是回归术语中对”后门未闭合”情况的表述。而条件独立假设，正是匹配方法中”所有后门路径均已闭合”这一核心要求的专业表述。
换言之，条件独立假设的核心要义在于：你所选取的匹配变量集合足以阻断所有后门路径。处理变量与结果变量之间的关联必须满足以下二者之一：\(（1）\) 确实存在你希望识别的前门因果路径；\(（2）\)源于你已测量并纳入匹配变量集的混杂因素。
这不过是用新术语重新表述了我们已讨论过的核心假设。那么除此之外，还需要考虑哪些关键因素呢？
14.6.2 共同支持
无论采用何种匹配方法，其核心前提是必须存在可供匹配的合适对照组样本。若该条件不满足，将导致匹配失效——这种现象在因果推断中称为共同支持域失败\(（failure \space of \space common \space support）\)。
为何这会成为问题？我们以一个极端案例说明：假设研究者希望评估”提前释放计划”对累犯率的影响——该计划允许服刑人员在刑满前提前出狱（通常附带特定条款）51。现有数据包括两组人员：实验组（参与提前释放的刑满人员）与对照组（服满刑期的刑满人员），以及他们出狱后十年内的再犯罪记录。
研究者推断，囚犯在服刑期间的行为表现可能同时影响其获得提前释放的几率及再犯罪概率。因此，研究采用监狱审查委员会制定的”行为评分”作为匹配变量——该委员会正是决定是否批准提前释放的决策机构。
然而问题出现了！监狱评审委员会的\(1-10\)分行为评分与提前释放决定存在完全对应关系：所有\(8-10\)分者均获准提前释放，\(1-7\)分者则全部遭拒。这导致研究者试图为实验组（\(8-10\)分获释者）寻找匹配对照组时，需要找出那些评分相同（\(8-10\)分）却未获释的样本——而此类样本根本不存在。此时分析就面临共同支持域缺失问题，即缺乏可比对照组样本这一根本缺陷。
再设想另一种情况：若采用加权分析法而非样本匹配，我们应如何为\(1-7\)分的对照组样本分配权重，才能使其与处理组（\(8-10\)分）具有可比性？实际上，这些样本的合理权重只能是零——因为低分段群体中完全不存在符合处理组特征的个体。最终结果与匹配法殊途同归：仍因共同支持域缺失而无法进行有效比较。
即使情况稍缓，问题依然存在。假设\(9-10\)分者全部获准提前释放，而\(8\)分者部分获释、部分未获释——我们仍无法为\(9-10\)分的处理组找到理想匹配对象。即便存在极少数\(9-10\)分未获释的样本，其数量也过于稀少，难以支撑可靠的因果推断。
共同支持域假设的核心要求是：处理组与对照组在匹配变量的分布上必须存在实质性重叠。若采用倾向得分匹配，则体现为两组的倾向得分分布必须具备显著重合区间。
检验共同支持域的方法相对简单，且存在多种实现方式：
检验共同支持域的第一种方法是直接对比处理组与对照组在关键变量上的分布差异。如\(图14.7（a）\)所示，该图呈现了我们基于布鲁克曼\(（2013）\)研究中的匹配变量——选区非裔人口占比\(（blackpercent）\)——在两组的分布情况。
## BROOCKMAN STUFF
br <- read_dta('Matching/broockman2013.dta')
m <- glm(leg_black ~ medianhhincom + blackpercent + leg_democrat, data = br, family = binomial(link = 'logit'))
br <- br %>%
mutate(ps = predict(m, type = 'response')) %>%
mutate(ipw = case_when(leg_black == 1 ~ 1/ps,leg_black == 0 ~ 1/(1-ps)))
# Common support
p1 <- ggplot(br, aes(x = blackpercent, linetype = factor(leg_black))) +
geom_density(size = 1) +
guides(linetype = FALSE) +
scale_x_continuous(labels = function(x) scales::percent(x, accuracy = 1)) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .3, y = 15, label = 'Non-Black Legislator') +
annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .6, y = 4, label = 'Black Legislator') +
labs(x = 'Percentage Black in District',y = 'Density',title = '(a) Percentage Black Common Support') +
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = 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.## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p2 <- ggplot(br, aes(x = ps, linetype = factor(leg_black))) +
geom_density(size = 1) +
guides(linetype = FALSE) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .3, y = 40, label = 'Non-Black Legislator') +
annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .6, y = 6, label = 'Black Legislator') +
labs(x = 'Propensity Score',y = 'Density',title = '(b) Propensity Score Common Support') +
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"))
plot_grid(p1,p2)
ggsave('Matching/Matching/common_support.pdf', width = 8, height = 3.5,device=cairo_pdf)这个变量是否具有良好的共同支持域？这有点不确定。分布不同并非问题 —— 我们可以使用匹配或重新加权的方法，将对照组分布中相应部分放大，直到它们与处理组分布相同，但是前提是对照组分布中存在这些可被放大的部分。如果分布的某部分没有数据，我们就无能为力了。显然，这里有些存在一些缺乏足够对照的区域。黑人人口占比超过\(50\%\)但立法会议员不是黑人的选区并不多。这对共同支持域来说是个问题，因为我们所研究的许多黑人议员代表的是黑人人口占比为\(50\%-80\%\)的地区。
如果存在黑人人口占比为 \(50\%-80\%\) 但议员不是黑人的地区，这些地区将被赋予极大权重，从而增加估计中的噪声。或者，如果不存在这样的地区，那么那些黑人人口占比为 \(50\%-80\%\) 且议员是黑人的地区就不得不被剔除。而这部分数据占了处理组数据的大部分！这样一来，对照结果就无法很好地代表那些有黑人议员的实际群体了。
这种分布对比方法同样适用于倾向得分分析。如\(图14.7（b）\)所示，处理组与对照组的倾向得分分布呈现预期差异：处理组得分整体右偏，这本身并不构成问题52。真正的方法论问题在于——当处理组某些区间的得分在对照组中完全或几乎不存在对应分布时（尤见于高倾向得分区间）。本研究确实存在此类问题：具有高倾向得分的非非裔议员选区样本严重不足。事实上，对照组的倾向得分密集堆积在\(0\)值附近，导致该区域密度异常膨胀，不仅压缩了图表其他区域的显示空间，更造成关键分布特征的辨识困难。
在倾向得分分析中，我们不妨借此机会检验另一个关键假设：真实的倾向得分必须满足绝对非确定性条件（即得分值永远不等于\(0\)或\(1\)）53。若发现得分分布大量堆积在\(0\)或\(1\)附近，则需警惕模型设定可能存在根本缺陷。
另一种检验共同支持域的直观方法是评估匹配成功率：若仅有\(1\%\)的处理组样本未能找到合格匹配对象，表明数据质量良好；反之若\(90\%\)的处理组样本匹配失败，则确证共同支持域缺失。
如何处理支持问题？我们在处理支持问题时可以采取的第一步是避免在没有支持的情况下尝试匹配。也就是说，我们希望舍弃那些匹配效果不佳的观测数据。
某些匹配方法会为我们自动完成这一步骤。如果我们采用卡钳值进行匹配筛选，那么不良匹配项就不会被选为匹配对象。缺乏共同支持的观测数据将无法进入匹配样本。搞定！
其他匹配方法则需要在事后进行一些额外处理。倾向得分方法（尤其是逆概率加权法）不会自动剔除共同支持范围之外的观测值。因此，通常需要对倾向得分进行”修剪”。具体操作是：分别对处理组和对照组，找出倾向得分密度分布为\(0\)（或接近\(0\)，低于某个截断值）的数值范围，然后删除另一组中所有倾向得分落在本组不存在的数值范围内的观测值。换句话说，就是将数据限制在共同支持的范围内。
Marco Caliendo and Sabine Kopeinig. Some practical guidance for the implementation of propensity score matching. Journal of Economic Surveys, 22(1):31–72, 2008.\(图14.8\)复现了\(图14.7（b）\)的结果，但进行了共同支持区域的修剪处理，要求每个宽度为\(0.02\)的倾向得分区间内至少包含\(10\)个对照组观测值54。在我们保留的数据中，每个观测值都至少存在若干个倾向得分相近的匹配对象。现在我们可以确信：当我们通过加权使两组样本更加相似时，我们不会强行纳入那些本质上不可比较的观测值——只保留那些能用于构建可比组的样本。
br <- br %>% mutate(cutps = cut(ps, breaks = (0:50)/50))
table(br %>% dplyr::filter(leg_black==0) %>% pull(cutps))##
## (0,0.02] (0.02,0.04] (0.04,0.06] (0.06,0.08] (0.08,0.1] (0.1,0.12]
## 4466 312 95 56 45 38
## (0.12,0.14] (0.14,0.16] (0.16,0.18] (0.18,0.2] (0.2,0.22] (0.22,0.24]
## 24 23 15 15 14 15
## (0.24,0.26] (0.26,0.28] (0.28,0.3] (0.3,0.32] (0.32,0.34] (0.34,0.36]
## 8 4 10 5 5 5
## (0.36,0.38] (0.38,0.4] (0.4,0.42] (0.42,0.44] (0.44,0.46] (0.46,0.48]
## 3 0 5 3 2 0
## (0.48,0.5] (0.5,0.52] (0.52,0.54] (0.54,0.56] (0.56,0.58] (0.58,0.6]
## 5 3 2 2 0 3
## (0.6,0.62] (0.62,0.64] (0.64,0.66] (0.66,0.68] (0.68,0.7] (0.7,0.72]
## 2 1 4 0 2 2
## (0.72,0.74] (0.74,0.76] (0.76,0.78] (0.78,0.8] (0.8,0.82] (0.82,0.84]
## 3 1 3 3 4 1
## (0.84,0.86] (0.86,0.88] (0.88,0.9] (0.9,0.92] (0.92,0.94] (0.94,0.96]
## 2 3 2 4 3 4
## (0.96,0.98] (0.98,1]
## 2 0ggplot(br %>% dplyr::filter(ps <= .24 ), aes(x = ps, linetype = factor(leg_black))) +
geom_density(size = 1) +
guides(linetype = FALSE) +
expand_limits(x = 1) +
annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .18, y = 50, label = 'Non-Black Legislator') +
annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .2, y = 10, label = 'Black Legislator') +
annotate(geom = 'rect', xmin = .24, xmax = 1, ymin = -1, ymax = 1, color = 'white', fill = 'white') +
labs(x = 'Propensity Score',y = 'Density') +
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"))
ggsave('Matching/Matching/trimmed_common_support.pdf', width = 6, height = 5,device=cairo_pdf)这并非修剪共同支持区域的唯一方法。常见的做法是修剪异常值——任何处于倾向得分最高或最低\(X\%\)分位数的观测值都会被剔除。也可以不采用百分位数，直接剔除接近\(0\)或\(1\)的数值（比如\(0-0.01\)和\(0.99-1\)之间的观测值），毕竟我们假设真实数据中本就不应存在绝对为\(0\)或\(1\)的倾向得分。另一种方法是仅要求区域内存在任意匹配项（而不要求足够数量的匹配项），例如：先确定处理组观测值的倾向得分最小值和最大值，然后删除对照组中所有超出该范围的观测值。
所有这些确保共同支持的方法都不可避免地需要舍弃部分数据。这让我们触及整个匹配过程中最棘手的难题——如果这一过程导致大量观测值（尤其是处理组的观测值）被剔除，那么往好了说，我们只能宣称自己获得了”匹配良好”个体的处理效应估计（尽管”匹配良好”的定义本身可能存疑）；往坏了说，这甚至意味着匹配方法彻底失败。匹配分析的核心前提在于存在可比较的观测值，而现实中这样的观测值可能根本不存在。或许换用其他匹配方法能解决问题，但更可能的情况是：匹配方法根本无法为你有效关闭那些混杂因素的后门。
14.6.3 平衡性检验
匹配分析的根本目标是通过匹配变量关闭后门路径。那么问题来了：我们成功了吗？平衡性\(（Balance）\)这一概念的核心假设就是：我们的匹配方法已经针对感兴趣的变量关闭了后门路径。理论上（在理想情况下），我们选择的权重应该能使处理组和对照组在所有后门变量上具有完全相同的取值分布。当某个后门变量在处理组和对照组之间不再存在差异时，就意味着该后门路径已被彻底关闭。
“完全相同的取值分布”在现实中其实很难实现。那么，我们该如何判断平衡性是否达到可接受的水平55？如果未能达标又会产生什么后果？与其他允许低质量匹配的情况类似，平衡性不佳会导致估计偏差——当平衡性不理想时，后门路径实际上未被关闭，处理效应的估计结果就会受到所有这些棘手后门变量的污染。
值得庆幸的是，我们无需过多理论推演来判断是否存在问题——通过已匹配的数据，我们可以直接检验处理组与对照组在理论上本应平衡的变量上是否仍存在显著差异。
若通过上述任一方法检测到平衡性问题，我们需重新审视匹配流程——基于匹配的估计方法本质上是一个迭代过程：先执行匹配，检验平衡性；若结果不理想，则调整匹配方案（如增加变量、改变倾向得分模型函数形式、缩小卡钳值等），并循环往复直至达标56。
检验平衡性的一种常用方法是制作平衡性检验表\(（balance \space table）\)。这类表格的构建逻辑相当直观：首先选取理论上应保持平衡的若干变量，然后分别展示处理组与对照组的描述性统计量（通常包括均值、标准差，有时也会补充其他指标），最后对每个变量进行均值差异检验，判断组间差异是否显著。更完善的表格可能还会添加一些进阶指标，例如标准化均值差\(（standardized \space difference \space in \space means）\)——即用组间均值差值除以合并标准差（听起来很高大上吧！）57。此外，您可能在某些平衡性检验表中看到名为“标准化偏差”\(（standardized \space bias）\)的变体指标，其计算原理与标准化均值差类似，仅在某些标准化细节上存在差异58。
平衡性检验表通常需制作两次：首次使用原始未匹配数据，用以展示匹配方法所需解决的平衡性问题严重程度；第二次则使用匹配后数据，理想情况下应显示平衡性问题已消除。良好的平衡性应呈现何种状态？其核心特征在于：所有变量的组间均值差异均无实际意义\(（meaningfully \space differences）\)。若进行统计检验，则不应发现任何具有统计学显著性的均值差异——或至少，显著差异的数量不应超过随机波动导致的预期值59。
平衡表通常很容易编写代码。在\(R\)中，\(Matching\)包中的\(MatchBalance\)和\(vtable\)包中的\(sumtable\)都有很多不错的选择，这两个函数既可以在不使用匹配权重和使用匹配权重的情况下运行，也可以结合\(Match\)的输出来比较匹配前后的平衡性。在\(Stata\)中，你可以在\(teffects\)分析之后使用 \(tebalance \space summarize\)（尽管这不会进行显著性检验）。要在\(Stata\)中获得更标准的平衡表，你可以使用\(balancetable\)包，添加权重来检查匹配后的平衡性。在\(Python\)中，如果你使用\(causalinference\)包创建了一个名为\(M\)的\(CausalModel()\)，可以使用\(print(M.summary_stats)\)来获取平衡表。
我们可以通过布鲁克曼\(（2013）\)研究的案例来观察平衡性检验表示例——该表格截取自\(R\)语言\(MatchBalance\)函数输出的部分结果（基于本章前文所述的马氏距离匹配法）。需要特别说明的是，平衡性检验表并无统一标准格式，您所用软件生成的表格可能与此示例存在差异。
我们在这里能看到什么呢？首先，在匹配进行之前，存在严重的不平衡情况。处理组和对照组之间的均值在实际意义和统计上都有显著差异。匹配之后情况如何呢？平衡性看起来相当不错。对于我们关注的三个变量，处理组和对照组之间的均值没有大的、有实际意义的差异。就 “立法者是民主党人”（代码中的\(leg\_democrat\)）而言，均值差异实际上为\(0\)。
我们确实看到，匹配后，选区中黑人占比在统计上有显著差异（在\(95\%\)的置信水平下 ），\(t\) 统计量的 \(p\) 值低于 \(0.05\)。不过，其均值非常接近\(（0.515 \space vs \space 0.513）\)，实际意义不大。此外，要判断我们是否应该担忧，我们可能得查看完整的分布情况，看看分布是否存在差异60。
平衡性检验表通常以均值平衡为核心检测目标，主要关注处理组与对照组的变量均值是否一致。虽然理论上也可用于检验其他统计量（如方差、分位数等）的平衡性，但这类扩展应用在实际研究中较为少见。
这确实是个严重问题！试想：若处理组与对照组的匹配变量分布如\(图14.9\)所示——即便均值相同，但分布形态存在明显差异——您还能确信自己真的关闭了后门路径吗？
set.seed(1000)
# 生成处理组数据：1000个来自标准正态分布的随机数
tib <- tibble(x = rnorm(1000))
# 构建处理组的密度图
p1 <- tib %>%
ggplot(aes(x = x)) +
geom_density() + # 绘制密度曲线
geom_vline(aes(xintercept = mean(x)), linetype = 'dashed') + # 添加均值垂直线
# 在指定位置添加均值标签（使用scales包格式化数字）
annotate(x = mean(tib$x), y = .2, geom = 'label', family = 'sans',size = 10/.pt,
label = paste0('Mean ', scales::number(mean(tib$x), accuracy = .01))) +
labs(x = 'Matching Variable X',y = 'Density',title = '(a) Distribution of X in Treated Group') +
theme_pubr() + # 应用pubr主题
# 自定义文本样式
theme(text = element_text(size = 10, family = "sans"),axis.title.x = element_text(size = 10, family = "sans"))
# 生成对照组数据：500个来自N(-2,1)和500个来自N(2,1)的混合数据
tib2 <- tibble(x = c(rnorm(500, -2), rnorm(500, 2)))
# 构建对照组的密度图（结构与处理组类似）
p2 <- tib2 %>%
ggplot(aes(x = x)) +
geom_density() +
geom_vline(aes(xintercept = mean(x)), linetype = 'dashed') +
# 注意：这里仍使用tib的均值（可能需要改为tib2的均值？）
annotate(x = mean(tib$x), # <<<< 潜在问题：应改为mean(tib2$x)
y = .05,
geom = 'label',
family = 'sans',
size = 10/.pt,
label = paste0('Mean ', scales::number(mean(tib$x), accuracy = .01))
) +
labs(x = 'Matching Variable X',y = 'Density',title = '(b) Distribution of X in Control Group') +
theme_pubr() +
theme(text = element_text(size = 10, family = "sans"),axis.title.x = element_text(size = 10, family = "sans"))
# 使用cowplot包将两个图并排排列
cowplot::plot_grid(p1, p2)
# 保存为PDF（使用Cairo设备确保字体兼容性）
ggsave('Matching/bad_balance_distribution.pdf', width = 8,height = 3.5,device = cairo_pdf)这暴露出匹配方法的重大缺陷！当处理组与对照组的匹配变量呈现\(图14.9\)所示的分布差异时——即便两组均值完全相同——我们怎能确信真正消除了混杂偏误？均值平衡仅是必要条件，而分布形态的一致性才是充分保障。
检验处理组与对照组变量分布差异的方法，其实在\(图14.9\)中已有所体现：只需绘制两组（匹配后）变量的分布图并进行直观比对。与\(图14.9\)的唯一区别在于，通常我们会将两组分布叠加显示（如\(图14.10\)所示），这样能更清晰地识别潜在差异——通过叠加呈现，分布不匹配的问题便一目了然。
虽然分布存在差异并不必然导致匹配方法失效，但如\(图14.10\)所示的严重不匹配情况，可能需要重新审视匹配方案的设计。这种分布对比图适用于所有连续型匹配变量，而对于类别型或二元匹配变量，则可借鉴第\(3\)章介绍的离散分布分析工具进行类似检验。
tib3 <- tib %>%
mutate(Group = 'Treated') %>%
bind_rows(tib2 %>% mutate(Group = 'Control'))
ggplot(tib3, aes(x=x, linetype = Group, color = Group)) + geom_density(linewidth = 1) +
guides(linetype = FALSE, color = FALSE) +
scale_linetype_manual(values = c('solid','dashed')) +
scale_color_manual(values = c('darkgray','black')) +
annotate(x = 1.4, y = .3, geom = 'text', family = 'sans',size = 10/.pt, label = 'Treated') +
annotate(x = 2.5, y = .2, geom = 'text', family = 'sans',size = 10/.pt, label = 'Control', color = 'darkgray') +
labs(x = 'Matching Variable X',y = 'Density') +
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"))
ggsave('Matching/Matching/overlaid_balance.pdf', width = 6, height = 5,device=cairo_pdf)有一种情况几乎必须要用到这种叠加密度图，那就是在处理倾向得分的时候。你实际上只有一个真正重要的匹配变量，也就是倾向得分，所以你最好希望（处理组和对照组倾向得分的）密度分布能很好地重合。
你可以在\(图14.11\)中看到实际效果，该图展示了布鲁克曼研究中匹配变量 “选区黑人占比” 以及倾向得分的叠加密度图，二者都是在应用逆概率权重（并修剪掉倾向得分在 \(0\) 到 \(0.02\) 以及 \(0.98\) 到 \(1\) 之间的数据 ）之后呈现的。
p1 <- ggplot(br %>% dplyr::filter(ps >= .02 & ps < .98), aes(x = blackpercent, linetype = factor(leg_black), weight = ipw)) +
geom_density(size = 1) +
guides(linetype = 'none') +
ggplot2::annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .2, y = 1.3, label = 'Non-Black') +
ggplot2::annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .38, y = 3, label = 'Black Legislator') +
labs(x = 'Propensity Score',y = 'Density',size = 10/.pt, title = '(a) Percentage Black after Weighting') +
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"))
p2 <- ggplot(br %>% dplyr::filter(ps >= .02 & ps < .98), aes(x = ps, linetype = factor(leg_black), weight = ipw)) +
geom_density(size = 1) +
guides(linetype = 'none') +
ggplot2::annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .3, y = 5, label = 'Non-Black Legislator') +
ggplot2::annotate(geom = 'text', family = 'sans', size = 10/.pt,x = .35, y = 2, label = 'Black Legislator') +
labs(x = 'Propensity Score',y = 'Density',title = '(b) Propensity Score After Trimming\nand Inverse Probability Weighting') +
theme_pubr() +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"))
plot_grid(p1,p2)
ggsave('Matching/Matching/ps_after_weight.pdf', width = 8, height = 3.5,device=cairo_pdf)现在这些（分布情况）看起来怎样呢？它们并不完美。要是有一个非常好的倾向得分模型，处理组和对照组观测值的倾向得分分布真的应该几乎完全相同。这些分布比匹配前好多了，但仍然差别挺大。选区黑人占比的分布也是如此。这可能是个很好的信号，提示我们回到倾向得分模型，尝试添加一些变量来改善我们看到的平衡性。
叠加密度图自己动手制作相当容易，但在某些情况下，也有现成的打包好的版本可供使用。
library(tidyverse)
br <- causaldata::black_politicians
# We can estimate our own propensity score
m <- glm(leg_black ~ medianhhincom + blackpercent + leg_democrat ,data = br, family = binomial(link = 'logit'))
# Get predicted values
br <- br %>% mutate(ps = predict(m, type = 'response'))
# Create IPW weights
br <- br %>% mutate(ipw = case_when(leg_black == 1 ~ 1/ps,leg_black == 0 ~ 1/(1-ps)))
# Density plots for raw data
ggplot(br, aes(x = medianhhincom , color = factor(leg_black))) + geom_density()
# And using our matching weights
ggplot(br, aes(x = medianhhincom , color = factor(leg_black), weight = ipw)) + geom_density()
import pandas as pd
import numpy as np
import seaborn as sns
from causaldata import black_politicians
br = black_politicians.load_pandas().data
# Let's check the raw pre-matching balance of a given matching variable
# fig1 = sns.kdeplot(data = br, x = 'medianhhincom',hue = 'leg_black',common_norm = False)
# fig1.plot()
# Start the new plot
# fig1.get_figure().clf()
# Add weights from any weighting method to check post-matching density
# Here we have the ipw variable from our propensity score matching in previous code blocks (make sure you run those first!)
# fig2 = sns.kdeplot(data = br, x = 'medianhhincom',hue = 'leg_black', common_norm = False, weights = 'ipw')
# fig2.plot()倾向得分的应用为平衡性检验提供了另一重要维度：当样本实现充分平衡且匹配有效关闭后门路径时，在控制倾向得分后，匹配变量与处理状态应不存在显著关联（即达到平衡）。这一逻辑关系清晰体现在\(图14.12\)的因果图中——若阻断”\(TreatmentPropensity \rightarrow Treatment\)“路径，则协变量\(A\)、\(B\)、\(C\)与处理变量之间应不存在任何连通路径。
\(Dehejia和Wahba（2002）\)提出了一种使用分层检验\(（stratification \space test）\)来检验我们图表中这一含义的方法。在分层检验中，我们反复查看倾向得分的有限范围。通过将自己限制在倾向得分的一个狭窄区间内，我们实际上是在对其进行控制。所以，把倾向得分的整个范围划分成区间，然后检查在每个区间内，匹配变量和处理是否相关61。
Rajeev H Dehejia and Sadek Wahba. Propensity score-matching methods for nonexperimental causal studies. Review of Economics and Statistics, 84(1):151–161, 2002.14.7 基于匹配数据的效应估计
14.7.1 均值差异估计
确实，虽然本章已接近尾声，我们才刚进入处理效应估计的环节——这恰恰反映了匹配方法的核心特征：超过\(80\%\)的工作量都集中在前期准备阶段。具体而言，研究人员需要：\(（1）\)选择合适的匹配方法（如\(K\)近邻匹配或卡钳匹配）；\(（2）\)执行匹配并计算权重；\(（3）\)检验共同支持域和平衡性；\(（4）\)必要时调整匹配方法直至达到满意的平衡性。只有当这些前期工作全部完成后，处理效应估计本身反而变得直接而简单，通常只需计算加权后的均值差异即可。这种工作重心的分布充分说明，在匹配分析中，研究设计的严谨性远比最终的估计技术更为重要62。
手中有了一组匹配权重后，如何估计处理效应呢？只需计算处理组和对照组结果变量的加权均值，然后进行比较即可。好了，就这样，结束了。
好吧，其实并非真的这么简单。不过，这确实占了很大一部分。如果你唯一想做的就是得出处理对结果变量均值的效应，那么比较加权均值就能实现。
需要牢记的一个重要细节是，虽然处理效应的估计很简单，但该处理效应的标准误就没那么简单了。计算加权均值的标准误很容易。但你计算出的标准误不会考虑到权重是经过估计的，而且可能还会忽略匹配过程中有些观测值被剔除了这一情况。对数据进行预处理以匹配观测值并创建权重的步骤会引入一些不确定性，从而增大标准误。将这种不确定性纳入你的标准误计算中是很重要的。
有几种不同的方法可用于计算合适的标准误调整，其中大多数方法都与确定结果变量在给定处理和匹配变量（或来自匹配变量的权重 ）条件下的方差有关。然而，这些方法都相当复杂，而且并非本书的重点。在你准备深入研究那些描述你想要进行的标准误调整的计量经济学论文63之前，你可能会使用能自动完成这项工作的软件。所以，这里的关键要点是：“除非你愿意花时间把它做对，否则不要在匹配情形下自行计算处理效应的标准误，让基于匹配的估计软件为你代劳吧。”
有一种例外情况是使用自助法标准误，正如第\(13\)章关于回归的内容中所讨论的。自助法标准误可应用于匹配估计。在将自助法应用于匹配时，你首先有放回地随机对数据重新抽样（这样你得到的样本量与初始样本量相同，但有些观测值会被其他观测值的多个副本替代 ），然后从头开始重新执行匹配，最后在自助法样本中估计处理效应。多次重复这个过程，所有自助法样本中处理效应的标准差就是你处理效应的标准误。你可以看出这个过程是如何自然地纳入匹配过程带来的任何不确定性的 —— 它让这种不确定性直接融入我们的处理效应分布中。
从直观上看（因为你在主动纳入匹配过程对抽样分布的影响 ），以及从操作层面看（因为自助法标准误自己手动计算相当简单，如果所用命令中没有该功能，只需编写少量代码即可 ），自助法标准误对于匹配估计而言可能是个很棒的选择。话虽如此，自助法标准误只能用于”构建加权匹配样本”的匹配方法。对于”选择匹配对象”的方法，它们根本无法正常发挥作用，因为该过程所做的明确的纳入 / 排除决策无法转化为具有恰当抽样分布的自助法（阿巴迪和因本斯，2008 ）。所以，如果你使用的是加权方法，那就用自助法标准误吧。如果不是，要么在你的软件中使用适当的标准误调整方法，要么如果你胆子够大，你可以去找一种为配合选择匹配对象方法而设计的自助法变体，使其能正常发挥作用。
Abadie, Alberto, and Guido W. Imbens. 2008. “On the Failure of the Bootstrap for Matching Estimators.” Econometrica 76 (6): 1537–57.如果使用的命令本身不提供匹配过程的自助法功能，那么编写代码实现匹配过程的自助法会有点棘手，因为你不是在对单个命令进行自助法操作（或者说不是在回归命令中选择 “自助法标准误” 选项这么简单 ）。你需要对整个匹配过程从头到尾进行自助法操作。在本章后面，我会在双重稳健估计的应用中提供一些代码示例，第\(15\)章也会有。对这些代码稍作调整，应该就能帮助你为各类匹配过程计算自助法标准误了。
14.7.2 匹配与回归的融合应用
在很多情况下，我们想做的事比对比处理组和对照组的均值要复杂一点。比如，或许我们的研究设计取决于处理和结果之间关系的特定函数形式，这种函数形式很难用匹配法捕捉，但用回归法很容易。但我们仍想通过匹配法关闭后门路径64。那该怎么办呢？我们可以把匹配权重和回归模型结合起来，进行回归调整。
这里的方法相当直接。从匹配过程中，我们已经有了一组选定的匹配观测值，或者一组匹配权重。我们可以把这个子样本或者那些权重，应用到我们的回归模型中，就像我们在第\(13\)章处理样本权重那样。
实际上，我们已经做过这个流程 —— 之前提到的加权均值的基本差异，就可以通过把匹配权重应用到仅以处理为预测变量的回归中实现，即：\(\text{Outcome} = \beta_0 + \beta_1\text{Treatment} + \varepsilon\) 。效应的估计值就是 \(\hat{\beta_1}\) 。这也是我在代码示例中一直用的方法。在此基础上给回归添加其他内容，跨度并不大。
使用权重并非唯一的办法（而且，要是你已对那些权重有其他安排，比如用于调查样本权重，那使用匹配权重可能就不理想了 ）。将匹配与回归结合的另一种方法是，仔细思考倾向得分是什么 —— 它是一个变量，若对其加以控制，理应阻断从处理到结果的所有后门路径。所以对它进行控制！倾向得分（或者基于倾向得分的逆概率权重 ）可以作为控制变量加入回归。
实际上，既然我们讲到这一部分了，终于可以回过头来，进行布鲁克曼\(（2013）\)实际开展过的分析。他不只是做了匹配，然后据此估计处理效应。不是的！他采用了那些匹配权重，将其用于回归，这样一来，他不仅能看出黑人议员和其他议员在回复率上的差异，还能看出，对于黑人发出的信件，选区内外回复率的差异，在黑人议员和其他议员之间是否存在不同。这就需要在回归中加入交互项。更不用说，我们还可以加入大量其他控制变量。
library(cem); library(tidyverse); library(modelsummary)
br <- causaldata::black_politicians
# This copies the CEM code from the CEM section See that section's code for comments and notes
# Limit to just the relevant variables and omit missings (of which there are none in this data)
brcem <- br %>%
dplyr::select(responded, leg_black, medianhhincom, blackpercent, leg_democrat) %>%
na.omit() %>%
as.data.frame()
# Create breaks
inc_bins <- quantile(brcem$medianhhincom, (0:6)/6)
create_even_breaks <- function(x, n) {
minx <- min(x)
maxx <- max(x)
return(minx + ((0:n)/n)*(maxx-minx))
}
bp_bins <- create_even_breaks(brcem$blackpercent, 6)
ld_bins <- create_even_breaks(brcem$leg_democrat,2)
allbreaks <- list('medianhhincom' = inc_bins,'blackpercent' = bp_bins,'leg_democrat' = ld_bins)
c <- cem(treatment = 'leg_black', data = brcem, baseline.group = '1', drop = 'responded',
cutpoints = allbreaks, keep.all = TRUE)##
## Using 'leg_black'='1' as baseline group# Get weights for other purposes. Note this exact code only
# works because we didn't have to drop any NAs. If we did, lining things up would be trickier
br <- br %>% mutate(cem_weight = c$w)
m1 <- lm(responded ~ leg_black*treat_out + nonblacknonwhite + black_medianhh + white_medianhh +
statessquireindex + totalpop + urbanpercent,
data = br,
weights = cem_weight)
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | 0.561*** |
| (0.027) | |
| leg_black | -0.142*** |
| (0.035) | |
| treat_out | -0.351*** |
| (0.013) | |
| nonblacknonwhite | 0.089** |
| (0.037) | |
| black_medianhh | 0.191*** |
| (0.017) | |
| white_medianhh | -0.169*** |
| (0.015) | |
| statessquireindex | -0.516*** |
| (0.087) | |
| totalpop | 0.013*** |
| (0.001) | |
| urbanpercent | 0.123*** |
| (0.024) | |
| leg_black × treat_out | 0.229*** |
| (0.048) | |
| Num.Obs. | 4829 |
| R2 | 0.192 |
| R2 Adj. | 0.191 |
| AIC | 15973.7 |
| BIC | 16045.0 |
| Log.Lik. | -7975.827 |
| F | 127.431 |
| RMSE | 0.53 |
import pandas as pd
import statsmodels.formula.api as smf
from causaldata import black_politicians
br = black_politicians.load_pandas().data
# This copies the CEM code from the CEM section
# See that section's code for comments and notes
br = (br.assign(inc_bins = pd.qcut(br['medianhhincom'], 6),bp_bins = pd.qcut(br['blackpercent'], 6)))
treated = (br.query('leg_black == 1').groupby(['inc_bins','bp_bins','leg_democrat'])
.size().to_frame('treated').query('treated > 0'))## <string>:1: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.control = (br.query('leg_black == 0').groupby(['inc_bins','bp_bins','leg_democrat'])
.size().to_frame('control').query('control > 0'))## <string>:1: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.# Merge the counts back in
br = (br.join(treated, on = ['inc_bins','bp_bins','leg_democrat']).join(control, on = ['inc_bins','bp_bins','leg_democrat']))
# Create weights
br['weight'] = 0
br.loc[(br['leg_black'] == 1) & (br['control'] > 0), 'weight'] = 1
totalcontrols = sum(br.loc[br['leg_black']==0]['treated'] > 0)
totaltreated = sum(br.loc[br['leg_black']==1]['control'] > 0)
br['controlweights'] = (br['treated']/br['control'])*(totalcontrols/totaltreated)
br.loc[(br['leg_black'] == 0), 'weight'] = br['controlweights']## <string>:1: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise an error in a future version of pandas. Value '[nan nan nan ... nan nan nan]' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.# Now, use the weights to estimate the effect
m = smf.wls(formula = '''responded ~ leg_black*treat_out + nonblacknonwhite + black_medianhh + white_medianhh +
statessquireindex + totalpop + urbanpercent''',weights = br['weight'],data = br).fit()
m.summary()| Dep. Variable: | responded | R-squared: | 0.140 |
|---|---|---|---|
| Model: | WLS | Adj. R-squared: | 0.136 |
| Method: | Least Squares | F-statistic: | 39.83 |
| Date: | 周日, 17 8月 2025 | Prob (F-statistic): | 2.91e-66 |
| Time: | 22:07:55 | Log-Likelihood: | -2947.6 |
| No. Observations: | 2221 | AIC: | 5915. |
| Df Residuals: | 2211 | BIC: | 5972. |
| Df Model: | 9 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 0.3741 | 0.038 | 9.743 | 0.000 | 0.299 | 0.449 |
| leg_black | -0.1793 | 0.038 | -4.732 | 0.000 | -0.254 | -0.105 |
| treat_out | -0.3420 | 0.021 | -15.920 | 0.000 | -0.384 | -0.300 |
| leg_black:treat_out | 0.1943 | 0.053 | 3.682 | 0.000 | 0.091 | 0.298 |
| nonblacknonwhite | -0.1127 | 0.046 | -2.425 | 0.015 | -0.204 | -0.022 |
| black_medianhh | 0.1060 | 0.034 | 3.075 | 0.002 | 0.038 | 0.174 |
| white_medianhh | -0.0070 | 0.023 | -0.307 | 0.759 | -0.052 | 0.038 |
| statessquireindex | 0.3432 | 0.132 | 2.597 | 0.009 | 0.084 | 0.602 |
| totalpop | 0.0031 | 0.002 | 2.024 | 0.043 | 9.62e-05 | 0.006 |
| urbanpercent | 0.0428 | 0.040 | 1.059 | 0.290 | -0.036 | 0.122 |
| Omnibus: | 266.053 | Durbin-Watson: | 1.963 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 1259.837 |
| Skew: | 0.477 | Prob(JB): | 2.69e-274 |
| Kurtosis: | 6.564 | Cond. No. | 165. |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
causaldata black_politicians.dta, use clear download
* Use any matching method to get weights
* We'll use coarsened exact matching like in the original
cem medianhhincom blackpercent leg_democrat(#2), tr(leg_black)
reg responded i.treat_out##i.leg_black nonblacknonwhite leg_democrat ///
leg_senator south blackpercent black_medianhh white_medianhh ///
statessquire totalpop urbanpercent [iweight=cem_weights]从这些回归中，我们得到\(leg\_black\)的系数为\(-0.142\)，\(treat\_out\)的系数为\(-0.351\)，\(leg\_black \times treat\_out\)的系数为\(0.229\)。回顾一下我们在第\(13\)章学到的对含交互项模型的解释方法，我们知道，这意味着与白人议员相比，黑人议员回复选区內邮件的可能性低 \(14.2\) 个百分点；对于白人议员来说，选区外邮件得到回复的可能性比选区內邮件低 \(35.1\) 个百分点。但就选区內与选区外的回复差距而言，黑人议员比白人议员小 \(22.9\) 个百分点，这一结果支持了”黑人议员有帮助黑人的内在动机”这一观点。
如果我们操作得当，回归和匹配结合起来的效果会大于它们各自效果之和。双重稳健估计是一种将回归和匹配相结合的方法，即便匹配存在问题或者回归存在问题（只要不是两者同时有问题 ），这种方法依然有效。
双重稳健估计\(（Doubly \space robust \space estimation）\)： 一种即便回归模型或匹配模型中的一个存在设定错误，仍能识别效应的估计方法。
我说的 “存在问题” 是什么意思呢？它并非涵盖所有情况 —— 要知道，你仍然需要实际测量并纳入那些能让你关闭后门路径的变量65。双重稳健估计可没法让你绕过这一点。但回归和匹配都依赖于假设。比如，你的回归模型可能设定有误 —— 模型中纳入了所有恰当的变量，但函数形式或许不太对。在设定用于估计倾向得分的模型时，也可能出现同样的情况。
双重稳健估计是这样一种方法：只要回归模型和倾向得分模型这两个模型中有一个设定恰当，该方法就仍然有效。另一个模型可以是错误的。
双重稳健性是某些估计量所具备的性质，而非一种特定的估计量。实际上，有很多方法都具有双重稳健性66。我不用费太多脑筋，现在就能列出十来种。但咱们没那么多时间。所以，我们来看三种方法，第一种会详细讲解。
我们要考察的第一个估计量是伍德里奇\(（Wooldridge,2010,第21.3.4章）\)中描述的那个。这个估计量很有吸引力，因为要是你刚听说双重稳健估计，然后得自己去构思它的实现过程，那实际操作起来，几乎就和你猜想的一样，而且它确实有效！
步骤相当简单明了：
利用匹配变量，为每个观测值估计倾向得分\(p\)。
用倾向得分生成逆概率权重：对于处理组观测值，权重为 \(1/p\) ；对于对照组观测值，权重为 \(1/(1 - p)\) 。
仅使用处理组观测值，用匹配变量作为预测因子（根据具体研究，也可能加入其他预测因子 ），结合逆概率权重，估计回归模型。
重复步骤\(3\)，但仅使用对照组观测值。
用步骤\(3\)和\(4\)得到的模型，为整个样本生成 “处理后” 和 “未处理” 的预测观测值。
对比 “处理后” 均值与 “未处理” 均值，得到因果效应的估计值。
如需标准误，可使用自助法标准误，或伍德里奇\(（Wooldridge,2010）\)中介绍的异方差稳健标准误（并非普通的异方差稳健标准误 ）。
虽然步骤不少，但每个步骤都不复杂，手工操作也很容易。实际上，咱们来看看这种方法的代码实现过程。
library(boot); library(tidyverse)
br <- causaldata::black_politicians
# Function to do IPW estimation with regression adjustment
ipwra <- function(br, index = 1:nrow(br)) {
# Apply bootstrap index
br <- br %>% slice(index)
# estimate and predict propensity score
m <- glm(leg_black ~ medianhhincom + blackpercent + leg_democrat, data = br, family = binomial(link = 'logit'))
br <- br %>% mutate(ps = predict(m, type = 'response'))
# Trim control observations outside of treated PS range
minps <- br %>%
dplyr::filter(leg_black == 1) %>%
pull(ps) %>%
min(na.rm = TRUE)
maxps <- br %>%
dplyr::filter(leg_black == 1) %>%
pull(ps) %>%
max(na.rm = TRUE)
br <- br %>% dplyr::filter(ps >= minps & ps <= maxps)
# Create IPW weights
br <- br %>% mutate(ipw = case_when(leg_black == 1 ~ 1/ps,leg_black == 0 ~ 1/(1-ps)))
# Estimate difference with IPW weights
w_means <- br %>%
group_by(leg_black) %>%
summarize(m = weighted.mean(responded, w = ipw)) %>%
arrange(leg_black)
return(w_means$m[2] - w_means$m[1])
}
b <- boot(br, ipwra, R = 200)
# See estimate and standard error
b##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = br, statistic = ipwra, R = 200)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.1563906 -0.01709766 0.1345327import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from causaldata import black_politicians
br = black_politicians.load_pandas().data
# As mentioned, est_via_weighting from causalinference already does doubly robust estimation but it's a different kind!
# Let's do Wooldridge here.
def ipwra(br):
# Estimate propensity
m = smf.logit('leg_black ~ medianhhincom + blackpercent + leg_democrat', data = br)
# Get fitted values and turn them into probabilities
m = m.fit().fittedvalues
br['ps'] = np.exp(m)/(1+np.exp(m))
# Trim observations outside of treated range
minrange = np.min(br.loc[br['leg_black']==1, 'ps'])
maxrange = np.max(br.loc[br['leg_black']==1, 'ps'])
br = br.loc[(br['ps'] >= minrange) & (br['ps'] <= maxrange)]
# Get inverse probability score
br = br.assign(ipw = lambda x: x.leg_black*(1/x.ps) + (1-x.leg_black)*(1/(1-x.ps)))
# Regress treated and nontreated separately, then predict for whole sample
mtreat = smf.wls('''responded ~ medianhhincom + blackpercent + leg_democrat''',
weights = br.loc[br['leg_black'] == 1, 'ipw'],
data = br.loc[br['leg_black'] == 1])
mcontrol = smf.wls('''responded ~ medianhhincom + blackpercent + leg_democrat''',
weights = br.loc[br['leg_black'] == 0, 'ipw'],
data = br.loc[br['leg_black'] == 0])
treat_predict = mtreat.fit().predict(exog = br)
con_predict = mcontrol.fit().predict(exog = br)
# Compare means
diff = np.mean(treat_predict) - np.mean(con_predict)
return diff
# And a wrapper function to bootstrap
def ipwra_b(br):
n = br.shape[0]
br = br.iloc[np.random.randint(n, size=n)]
diff = ipwra(br)
return diff
# Run once on the original data to get our estimate
est = ipwra(br)## Optimization terminated successfully.
## Current function value: 0.089701
## Iterations 10# And then a bunch of times to get the sampling distribution
dist = [ipwra_b(br) for i in range(0,50)]## Optimization terminated successfully.
## Current function value: 0.092291
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.086002
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.094522
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.092213
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.087896
## Iterations 9
## Optimization terminated successfully.
## Current function value: 0.094587
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.092912
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.084807
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.067998
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.085098
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089612
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.085776
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089443
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089732
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.088962
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.081135
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.091747
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089611
## Iterations 9
## Optimization terminated successfully.
## Current function value: 0.080298
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.092124
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089043
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.087070
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.094060
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.094840
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.085571
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089137
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.084418
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.086574
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.087447
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.086830
## Iterations 9
## Optimization terminated successfully.
## Current function value: 0.086121
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.095475
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.087215
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.082761
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.095160
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.086958
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.098115
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.077682
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.093067
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.093612
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.092436
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.099829
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.077854
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.088307
## Iterations 11
## Optimization terminated successfully.
## Current function value: 0.092550
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.094784
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.082754
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.097212
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.095991
## Iterations 10
## Optimization terminated successfully.
## Current function value: 0.089562
## Iterations 10# Our estimate
est## np.float64(0.04992616176452569)# and its standard error
np.std(dist)## np.float64(0.0942115105679956)causaldata black_politicians.dta, use clear download
* Wrap this in a program so we can bootstrap
capture program drop ipwra
program def ipwra, rclass
quietly{
* Estimate propensity
logit leg_black medianhhincom blackpercent leg_democrat
* Get predicted values for our propensity score
predict ps, pr
* "Trim" control observations outside of
* treated propensity score range
summ ps if leg_black
replace ps = . if ps < r(min) | ps > r(max)
* Get inverse probability weights
g ipw = 1/ps if leg_black
replace ipw = 1/(1-ps) if !leg_black
* Regress treated and nontreated separately,
* then predict for whole sample
regress responded medianhhincom blackpercent ///
leg_democrat [pw = ipw] if !leg_black
predict untreated_outcome
regress responded medianhhincom blackpercent ///
leg_democrat [pw = ipw] if leg_black
predict treated_outcome
* Compare means
summarize untreated_outcome
local unt_mean = r(mean)
summarize treated_outcome
local tr_mean = r(mean)
local diff = `tr_mean' - `unt_mean'
}
drop ps ipw untreated_outcome treated_outcome
return scalar diff = `diff'
end
* And a version for bootstrapping
capture program drop ipwra_b
program def ipwra_b, rclass
preserve
* Create bootstrap sample
bsample
* Estimate
ipwra
* Bring back original data
restore
return scalar diff = r(diff)
end
* Run once with original data
ipwra
local effect = r(diff)
* Then boostrap
simulate diff = r(diff), reps(2000): ipwra_b
local N = _N
bstat, stat(`effect') n(`N')
* Or skip all that and just use teffects ipwra (no trim here though)
causaldata black_politicians.dta, use clear download
teffects ipwra (responded medianhhincom blackpercent leg_democrat) ///
(leg_black medianhhincom blackpercent leg_democrat, logit)我提到的另外两种方法怎么样呢？我简单讲讲这些，算是大致了解一下。
一种是增强逆概率加权估计量\(（Augmented \space Inverse \space Probability \space Weighted \space Estimator,AIPWE）\)（谭，\(2010\) ）。和伍德里奇的方法类似，\(AIPWE\)既要估计一个以处理为因变量的匹配模型，也要估计一个以结果为因变量的回归模型。\(AIPWE\)的基本思路是，在回归模型中纳入一个增强项（想不到吧 ），该增强项会根据匹配模型的设定误差程度对结果进行调整。
Tan, Zhiqiang. 2010. “Bounded, Efficient and Doubly Robust Estimation with Inverse Weighting.” Biometrika 97 (3): 661–82.\(AIPWE\)有一些不错的性质。即便两个模型都存在轻微设定误差，在某些情况下它也能表现不错（当然，还是更希望至少有一个模型是正确设定的 ）。此外，如果结果是处理模型不存在设定误差，\(AIPWE\)会使结果模型的标准误更小。厉害吧！
第三种也是最后一种方法可以简单讲讲，因为我之前已经讨论过。熵平衡在之前的 “距离匹配” 部分出现过，它本身就具有双重稳健性（赵启、珀西瓦尔，\(2017\) ）。尽管没有直接纳入回归，但熵平衡专注于消除处理组和对照组之间的均值差异，结果其作用方式和回归专注于消除均值差异很相似。双重稳健性就这么体现出来了！
Zhao, Qingyuan, and Daniel Percival. 2017. “Entropy Balancing Is Doubly Robust.” Journal of Causal Inference 5 (1).14.8 匹配与处理效应
和任何设计及估计方法一样，我们必须弄清楚，我们得到的是何种处理效应均值（正如第\(10\)章所讲 ）。
谈到匹配时，这成了一个非常有意思的问题。匹配最妙的地方之一在于，和回归不同，要让匹配估计量给出你想要的几乎任何一种处理效应均值，都超级容易67。
这怎么可能呢？嗯，首先，咱们先想想本章一直讨论的匹配程序 —— 这些程序会给我们带来哪种处理效应均值呢？
我们可以从逻辑上思考这个问题68。我们用匹配来生成一组与处理组相似可比的对照组观测值。这些对照组观测值应该和 “若处理组未接受处理时的情况” 一样。
所以，将我们对”处理组若未接受处理会得到何种结果”的估计，与”处理组实际得到的结果”进行对比？用第\(10\)章里的经验法则逻辑来看，在我看来，这听起来就像是”对处理组的平均处理效应”。事实上，本章到目前为止描述的大多数方法，给出的就是这个69。
或者说，至少非常接近它们会给出的结果。匹配不会给你完全精确的 “对处理组的平均处理效应”。和回归一样，它给出的均值会受其关闭后门路径方式的影响。就像回归给出的是方差加权的处理效应，匹配给出的是分布加权的处理效应，即观测值的匹配变量组合越常见，其处理效应的权重就越大。
为了实际举例说明这一点，咱们像第\(10\)章那样，设想一些不同类型的人群。\(表14.2\)展示了我们可能会抽样的两种不同类型的人群。
现在，假设我们最终得到一个样本，其中有\(500\)名未接受处理的阿尔弗雷德\(（Alfreds）\)和\(500\)名未接受处理的布里安娜\(（Briannas）\)，还有\(1000\)名接受处理的阿尔弗雷德\(（Alfreds）\)和\(200\)名接受处理的布里安娜\(（Briannas）\)。我们要基于性别进行匹配。
在处理组中，我们会看到平均结果为 \((1000 \times 4 + 200 \times 5)/(1200) \approx 4.17\) 。那么，匹配后的未处理组情况如何呢？比如，如果我们进行有放回的一对一匹配，最终会得到和处理组一样的情况，即\(1000\)名阿尔弗雷德和\(200\)名布里安娜，但这次他们是未接受处理的。所以，未处理组的平均值是 \((1000 \times 2 + 200 \times 1)/(1200) \approx 1.83\) 。这样，我们得到的处理效应就是 \(4.17 - 1.83 = 2.34\) 。这个结果更接近男性的处理效应\(（2）\)，而非女性的处理效应\(（4）\)，因为处理组样本中男性数量远多于女性 —— 它的权重取决于特定观测的匹配变量值的常见程度70。
14.8.1 从匹配中获得其他处理效应
啊，可我之前没把关键信息点出来71。 “对处理组的平均处理效应”并非匹配能给出的唯一结果。我们还能得到平均处理效应，甚至”对未处理组的平均处理效应”72。那要怎么操作呢？
其实相当简单。要是你想得到”对处理组的平均处理效应”，就像我们之前做的那样，构建一个尽可能与处理组相似的对照组。
想得到”对未处理组的平均处理效应”？简单！把操作反过来就行。不是让对照观测值匹配处理观测值，而是让处理观测值匹配对照观测值。流程一样，只是把两组的角色互换。这样得到的处理组会尽可能接近未处理组。将其与未处理组对比，就是把”未处理时实际得到的结果”和”我们认为该组若接受处理会得到的结果”进行对比。这就得到了”对未处理组的平均处理效应” ！73
从这里，你可以相当轻松地得到平均处理效应。如果处理组的均值是”对处理组的平均处理效应\(（ATT）\)“，未处理组的均值是”对未处理组的平均处理效应\(（ATUT）\)“，那么总体均值（也称为平均处理效应）是多少呢？嗯，如果有比例为\(X\)的人接受了处理，那么平均处理效应\(（ATE）\)的计算方式就是 \(ATE = (x \times ATT) + ((1 - x) \times ATUT)\) 。你其实就是对来自”两边”（处理组和未处理组 ）的均值做平均，按照两边人群规模进行加权，从而得到总体均值。
利用匹配来获取平均处理效应的做法相当常见，有时候在那些专门用于通过匹配分析处理效应的软件命令里，这甚至是默认选项。\(Stata\)软件中的 \(teffects\)命令就是例子 —— 本章里\(teffects\)的示例默认生成平均处理效应（而且这些示例都有\(atet\)选项，用于获取”对处理组的平均处理效应”）。
这和逆概率加权法有什么关联呢？毕竟，在进行逆概率加权时，你其实并没有真正去匹配什么，只是依据逆倾向得分权重来加权而已。
事实证明，你将倾向得分转化为逆倾向权重的具体方式，决定了你得到的处理效应均值的类型。加权是为了让对照组看起来像处理组？得到的就是”对处理组的平均处理效应”！反过来操作？得到的就是”对未处理组的平均处理效应”！两边都操作？得到的就是”平均处理效应”！
实际上，我一直讲的逆概率加权方法 —— 处理组权重为 \(1/p\) ，未处理组权重为 \(1/(1 - p)\) ，其中 \(p\) 是倾向得分 —— 会同时增加”与处理组最相似的未处理组成员”以及”与未处理组最相似的处理组成员”的权重。实际上，这给我们的是平均处理效应，而非”对处理组的平均处理效应”。
要是我们确实想要其他某种均值（指\(ATT\)或\(ATUT\)），该怎么做呢？我们只需使用不同的权重就行。
要是我们想要”对处理组的平均处理效应”，直观来看，我们要让未处理组变得像处理组。所以，首先给处理组所有成员相同的权重\(1\)。我们不需要调整他们去像其他组的人 —— 我们要调整其他所有人来像他们，这样其他人就能作为参照组。在我看来，这就对应”对处理组的平均处理效应”。
处理组权重设为\(1\)。那未处理组呢？不用 \(1/(1 - p)\) 这个权重，给他们 \(p/(1 - p)\) 。之前用的 \(1/(1 - p)\) ，有点像把未处理组的代表性拉回 “中性”，也就是所有人被平等计数。但我们不想要平等计数。我们想让未处理组像处理组，这意味着要依据他们与处理组的相似程度，给他们更高的权重。所以 \(p/(1 - p)\) 比 \(1/(1 - p)\) 对 \(p\) 的反应更强烈。
同理，“对未处理组的平均处理效应”也很容易推导出来。这种情况下，未处理组权重设为\(1\)，处理组不用能得到\(ATE\)的 \(1/p\) ，而是用 \((1 - p)/p\) ，这样就能让情况与未处理组合理匹配，得到我们想要的\(ATUT\)（对未处理组的平均处理效应 ）。
与回归分析的“找到所有与\(W\)的变化相关的变化，并通过删除它来关闭\(W\)所在的任何开放后门”相对。↩︎
尽管实际上大多数老板（和观众）要么对统计足够精通，可以理解回归，要么对统计一窍不通，以至于你可以简单地说“调整\(X\)中的差异”，他们不会真正关心你是怎么做的。↩︎
匹配方法同样适用于非二元处理变量，现有若干技术方案可实现这一点——例如选取处理变量的特定取值或取值区间进行匹配。但匹配方法的应用仍主要集中于二元处理情形，本章后续讨论也将聚焦于此种情况。↩︎
匹配方法中还存在另一种常见操作：通过选取与控制组样本最相似的处理组样本进行匹配。为简化说明，本章前半部分将暂时搁置这种方法，但后续会展开讨论。↩︎
在假定不存在其他未阻断的后门路径的前提下↩︎
关于”哪种方法更值得学习”这个问题，有趣的是，答案取决于你所处的领域。在某些学科（如经济学），回归分析是标准方法，匹配法则非常罕见；而在另一些领域（如社会学），两种方法都很常见；还有些学科则主要使用匹配法，回归分析反而少见。你可以选择只学习本领域的主流方法，忽略那些冷门技术。或者做个酷炫的反叛者——掌握所有方法，然后根据具体问题选择最佳工具。太帅了！让那些老古板瞧瞧——你们可以打压我们，但永远杀不死摇滚精神！呃，我是说杀不死匹配法。↩︎
你可以通过其他方式利用权重来估计效应，我将在本章后续部分进行阐述。↩︎
某些加权平均方法会省略“除以权重总和”的步骤，而是要求权重本身之和等于\(1\)↩︎
这些数字从哪儿来的呢？0.16×500 是这样来的：每个男性，“男性”取值为1，权重是0.16 。所以每个男性就是0.16×1 。把500个男性加起来，得到0.16×500 。女性“男性”取值为0，权重是0.04 ，把500个女性加起来，得到0.04×0 。 然后看分母，不管你是男性还是女性，权重都要算进去。所以我们把男性的500个0.16权重加起来，得到0.16×500 ，再加上女性的500个0.04权重，即0.04×500 。↩︎
这就是\(350\)名就业男性，每名权重为\(0.16\) ，加上\(275\)名就业女性，每名权重为\(0.04\) ，除以所有权重之和 ——\(500\)名男性，每名权重为\(0.16\) ，加上\(500\)名女性，每名权重为\(0.04\) 。↩︎
更常见的情况是：我们需要明确计算机执行匹配时的具体筛选标准。↩︎
当然，只挑选一个匹配对象并非唯一的匹配方式——我们之后会谈到这一点。↩︎
笔者推荐Caliendo和Kopeinig（2008）的权威指南，该文献对倾向得分匹配方法进行了深入浅出的阐述。本章写作过程中亦多次参考其方法论框架。↩︎
需特别注意的是：若匹配变量未能阻断后门路径，即使其能预测处理变量，反而可能导致估计偏差——这一观点在Bhattacharya和Vogt（2007）的研究中已获实证支持。↩︎
如”匹配质量检验”章节将详述，倾向得分匹配的核心特征在于：其不追求匹配变量的直接平衡。例如：
当变量\(A\)每增加\(1\)单位使处理概率提升\(0.1\)
变量\(B\)每增加\(1\)单位使处理概率提升\(0.05\)时，倾向得分匹配会判定\(A=2/B=1\)的样本与\(A=1/B=3\)的样本高度匹配——这与距离匹配的判断标准存在本质差异。
倾向得分匹配方法存在一个重要限制：当仅能阻断部分后门路径时，若未控制的开放路径变量与已匹配变量存在关联，反而可能通过倾向得分加剧剩余路径的匹配偏差。Brooks和Ohsfeldt（2012）对此现象进行了系统阐释。↩︎
这一看似简单的步骤，实则揭示了倾向得分匹配相较于距离匹配的核心难点：必须预先确定回归模型规格。鉴于我们刚用完整章篇幅讨论回归分析，您应该能预见这将引发哪些建模争议。↩︎
虽然\(logit\)或\(probit\)模型是估计倾向得分的标准方法，但该方法存在参数限制。在大样本情况下，采用非参数估计方法反而可能提高最终估计的精确度——详见Hahn（1998）的研究。↩︎
当仅使用单一匹配变量时，距离匹配与倾向得分匹配的差异微乎其微——因此两次匹配得到相同最优结果（\(第27281\)）实属必然。若重绘\(图14.2\)，仅需将\(X\)轴标签改为”倾向得分”即可，图示结构完全不变。这种等价性将在引入多变量匹配后发生本质变化。↩︎
通常匹配样本会获得相同权重，但存在例外情况：
单匹配场景：若处理样本\(1\)匹配\(1\)个对照，处理样本\(2\)匹配\(3\)个对照
常规方法：\(4\)个对照均分权重
替代方法：处理样本\(1\)的对照权重\(3\)倍于其他
多重匹配场景：同一对照样本匹配多个处理样本时，该对照将被重复计数
这种精确的加权方案虽非标准方法，但能帮助理解核心概念。↩︎
从技术层面看，样本选择法本质上是加权匹配的一种特殊形式——当所有权重非\(0\)即\(1\)时，即构成一种特殊的加权方案。但鉴于二者实操差异显著，本文仍将作为独立方法论述。↩︎
此处的例外是均匀核函数：对特定范围内的所有样本赋予相同权重，超出范围则骤降至零。↩︎
Epanechnikov核不仅广泛应用于基于核的匹配，在密度函数估计等其它核方法中也同样常用。↩︎
除了计算简便（许多核函数都具备此特点）外，Epanechnikov核的流行更因其能最小化渐近均方积分误差——该指标衡量变量估计密度与实际密度的差异。不过您可能无需深究这点。↩︎
标准化即减去均值后除以标准差。不过在本例中，减去均值的步骤实际可省略——因为后续计算涉及差值运算，均值最终会被抵消。↩︎
若某对照观测匹配多个处理样本，则将其所有权重累加。↩︎
我们终于完成了处理效应的估计！↩︎
需要说明的是，我并非宣称发现了唯一正确的方法——但当我自己进行匹配分析，或被问及匹配方法选择时，通常首推倾向得分逆概率加权法或熵平衡法（本章”熵平衡”章节将详述）。正如全章案例所示，虽存在其他方法更适用的场景，但这两种始终是我的优先选择。↩︎
我本想在此引用几篇代表性论文，但实话说这个领域文献浩如烟海。建议直接登录您常用的期刊检索系统，搜索”匹配最优带宽选择”关键词——相信您能找到大量相关研究。↩︎
例如研究某项政策变革对企业的影响时，您可能需要在多维度匹配相似企业，但严格要求行业属性必须精确匹配——如此可避免将电影院与摩托车修理店进行比较，即便它们在其它指标上极为相似。↩︎
这个邮件主题选择相当巧妙（恕我直言）：回复邮件能切实改善对方生活境况，而不仅关乎政策议题推进。↩︎
更强有力的证据需通过对照实验获得：当发件人改为白人姓名时：
预期结果：效应方向相反
若未出现：则可能表明黑人议员本身就更乐于助人
具体而言，他采用了粗化精确匹配结合回归分析的两步法——本章将详细探讨这两个步骤。↩︎
我必须声明这并非唯一方法——数据科学家们此刻可能倍感困惑，因为他们更常用欧氏距离（该指标忽略协方差调整步骤，后文将详述）。但无论如何，请允许教科书作者保留这难得的共识性结论。↩︎
关于本章代码示例的重要说明：我将展示这些命令的使用方法，但在匹配分析中，不同软件包的关键默认设置差异很大——使用何种标准误？如何处理最近邻匹配中的平局情况（当多个观测同等匹配时）？应该计算平均处理效应\(ATE\)还是处理组的平均处理效应\(ATT\)？请勿期望不同语言会输出相同结果。若需了解导致差异的具体参数及如何使其一致，必须深入研究各命令的文档。↩︎
如此宽的区间（64,000至112,000新台币）是否真能阻断\(BillApril\)的后门路径值得商榷。或许需要更精细的分箱方案——构建区间时务必审慎！↩︎
这并不代表真正消除了所有差异。如前所述，粗化处理的连续变量仅在分箱层面完全匹配，而箱内实际值仍存在差异。此外，未用于匹配的变量也可能存在差异。更严重的是，若因维度灾难而在粗化精确匹配中减少匹配变量——糟糕！↩︎
无需计算平方根或\(logit\)模型——非匹配即剔除。更重要的是：无需逐一样本匹配，仅需统计各分箱中处理组与对照组的观测数即可生成权重。↩︎
粗化精确匹配的另一应用场景：当数据集存在缺失值时，可在非缺失变量上寻找精确匹配，参照匹配样本的取值进行填补。更详细的缺失值处理方法参见第22章。↩︎
是的，这类数据集确实存在。北欧国家常见从出生到死亡的全生命周期追踪数据，其他国家也有类似但更有限的版本，例如使用美国IRS全部税收数据库的研究。无需多问，每位研究者都对能获取这些数据的人深感嫉妒——用一根手指换数据访问权？反正我还有九根呢。↩︎
该设计方法正是第\(18\)章将要讨论的双重差分法的雏形。↩︎
就个人理解，我将熵平衡法视为”某种用于匹配的矩估计方法”。这句话对多数读者而言可能完全不知所云——但若您恰好具备相关背景，那么这短短九个字其实相当精准地概括了熵平衡法的本质。事实上，学界确实存在与熵平衡法高度相似的实践：通过矩估计方法生成匹配权重。↩︎
若您还保留着多元微积分的知识记忆，就会发现熵平衡法本质上只是拉格朗日乘数法（Lagrange multiplier）的一种应用。↩︎
此外，该方法对软件实现可靠性的依赖度更高——难道您打算自行编写权重搜索算法？好吧，或许您确实有此能力。但对多数研究者而言，熵平衡法的技术门槛明显高于本章介绍的其他方法。当向统计学基础较弱的对象解释时，“我施加了矩条件约束”这样的表述，远不如”我找到了匹配变量值相似的样本”来得直观易懂，前者听起来简直像在施展魔法。↩︎
See https://github.com/google/empirical_calibration.
↩︎若您因回归章节篇幅过长而试图仅阅读匹配章节——恐怕要令您失望了。（因果推断中的匹配方法本质上仍依赖于回归技术框架）↩︎
在此场景下，线性概率模型（Linear Probability Model, LPM）的缺陷尤为明显：由于直接使用预测值作为倾向得分，任何超出[0,1]区间的预测结果均无法使用。因此，线性概率模型不适用于倾向得分估计。↩︎
注意！在某些编程语言中，对\(Logit\)或\(Probit\)模型使用标准的预测函数，不会默认输出概率预测值，反而会给出线性预测指数。查看你所用预测函数的文档，同时检查输出值，确保其在 \(0\) 到 \(1\) 之间。↩︎
有道理↩︎
假释\(（parole）\)即是提前释放计划的典型范例。↩︎
若实际情况并非如此（即对照组倾向得分反而更高），则将产生逻辑悖论——这意味着本应作为对照的群体反而比处理组更有可能接受干预。↩︎
这一限制的根本原因在加权计算中显而易见：若某处理组样本的倾向得分恰好为\(1\)（即\(100\%\)确定接受处理），其权重计算公式将出现\(1/0\)的分母为零错误；同理，对照组中得分为\(0\)的样本也会导致相同的数学悖论。↩︎
我采用的这种”最低数量即视为共同支持”的做法其实是主观设定的——完全是我自己制定的标准。实际上，当支持度不为零时，究竟达到多少才算足够，并不存在公认的经验法则。↩︎
本节所述内容同样适用于其他需要确保两组变量无差异的研究场景。例如：
随机实验——我们预期随机化分组自然会产生平衡性
第20章将讨论的断点回归——要求特定数据子集必须满足平衡性条件
这正是粗化精确匹配（coarsened exact matching）与熵平衡法（entropy balancing）等方法的优势所在——它们通过设计原理规避了平衡性问题，其效果之显著，甚至使得那篇开创性的粗化精确匹配论文直接以《无需平衡性检验的因果推断》为题。这些方法彻底摆脱了”匹配→检验→调整”的迭代循环！↩︎
若您未系统学习过基础统计课程却能坚持读到本书此处——说实话，我不知该作何感想。不过无论如何，对于未接触过均值差异检验的读者，其实质仅仅是检验某变量在两组的均值是否存在差异。您甚至可以通过回归实现：将目标变量对处理指示变量进行回归即可完成检验。↩︎
具体而言，该标准化差值等于组间均值差除以合并标准差（其计算方式为：先求处理组与对照组方差之和，乘以0.5后取平方根）。
（文献引用：Paul R. Rosenbaum与Donald B. Rubin在1985年发表于《美国统计学家》（The American Statistician）第39卷第1期33-38页的论文《运用包含倾向得分的多元匹配抽样方法构建对照组》中首次系统阐述了该方法。）↩︎
若平衡性检验表中包含大量变量，即使其中少数变量存在显著差异（特别是当均值差异本身实际幅度很小时），通常并不构成严重问题。↩︎
该平衡性检验表除t检验外，还包含Kolmogorov-Smirnov检验（简称K-S检验）结果。该非参数检验的核心功能是评估两组数据的整体分布一致性。值得注意的是，“选区非裔人口比例”变量的K-S检验p值显著高于临界水平（具体表现为p≈0.82），这从分布形态角度进一步佐证了前文所述”均值差异虽显著但实际影响甚微”的结论。↩︎
划分成多少个区间呢？这些区间应该窄到足以让倾向得分本身在每个区间内达到平衡。如果没达到，就划分得更窄些。但也不能窄到每个分层（strata ，这里指划分出的区间 ）内没有足够的观测值来进行检验。↩︎
紧急警示↩︎
再次强调，Caliendo 和 Kopeinig（2008）是个不错的起点。↩︎
这一点将在\(第18章\)关于双重差分法的内容中再次涉及，该方法依赖于对处理组和对照组（可能已进行匹配）进行比较，并精准识别结果与时间之间关系中存在的特定断点。↩︎
不这样做会导致两种方法都无法奏效 —— 记住，至少需要有一种方法能起作用。↩︎
我们如何知道某个特定的估计量具有双重稳健性呢？有人会撰写一篇包含大量方程式的论文来证明这一点。这类论文有很多。↩︎
至少在你通过关闭后门来识别效应的情况下是这样。当把匹配与前门隔离方法结合时，情况会变得更复杂一点。不过，（匹配）还是比回归更灵活。↩︎
这种解释并不适用于那些同样对处理组进行加权的方法，比如我所描述的逆概率加权法——这一点我稍后会讲到。↩︎
至少就我所描述的大多数方法而言是这样。有些代码示例给出的是平均处理效应，因为那些命令的默认设置就是生成平均处理效应。请查阅文档！我在下一节也会谈到这一点。↩︎
这是一种按处理组内的分布进行加权的方法，因为它是 “对处理组的平均处理效应（average - treatment - on - the - treated ）” 方法。不过，在计算平均处理效应、对未处理组的平均处理效应等情况时，分布加权会依据其他相关组中匹配变量的分布来进行。↩︎
对，就是这么拼写的。↩︎
当然，这要受我刚才提到的分布加权的影响。↩︎
这一切都假设，那些能预测处理效应变异的变量已被纳入匹配变量中。而且，如果你使用倾向得分进行匹配，这进一步假设，那些能预测处理效应变异的变量与那些能预测处理本身的变量是相同的。↩︎