26-因果推断中的协变量匹配
这里将继续使用 cfps2010.dta 这个数据演示上述命令的使用方法。对该数据基本情况的介绍请参见第四章第五节。
与第四章演示线性回归时相同，本章分析的因变量仍为 lninc，干预变量为 college。
一元线性回归结果显示，在不控制任何变量的情况下，college 对 lninc 的回归系数为 0.824，这意味着，考上大学的样本和没有考上大学的样本在 lninc 上的均值差为 0.824。不过，考虑到这两个子样本在 hukou、age、gender、race、sib-ling 和 fmedu 这几个协变量上存在很大差异，上述结果很可能是有偏的。
下面，我们将通过匹配法来校正干预组和控制组在这些协变量上的差异。
. use "../../Data/causal-inference/cfps2010.dta", clear
. reg lninc college,vce(cluster provcd)
Linear regression Number of obs = 4,137
F(1, 24) = 271.17
Prob > F = 0.0000
R-squared = 0.1095
Root MSE = 1.1498
(Std. err. adjusted for 25 clusters in provcd)
------------------------------------------------------------------------------
| Robust
lninc | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
college | .823612 .0500155 16.47 0.000 .7203851 .926839
_cons | 9.353189 .1084703 86.23 0.000 9.129317 9.577061
------------------------------------------------------------------------------
. 我们使用的第一种方法是精确匹配。首先，我们仅对 hukou 这一个变量实施精确匹配，具体如下：
Treatment-effects estimation Number of obs = 4,137
Estimator : nearest-neighbor matching Matches: requested = 1
Outcome model : matching min = 786
Distance metric: Mahalanobis max = 1594
------------------------------------------------------------------------------
| AI robust
lninc | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATE |
college |
(是 vs 否) | .8021249 .0337823 23.74 0.000 .7359127 .868337
------------------------------------------------------------------------------可以发现，在对 hukou 实施了精确匹配之后，college 对 lninc 的平均干预效应下降到了 0.802。
如果我们同时对 hukou、age、gender、race、sibling 和 fmedu 实施精确匹配，其结果又会如何呢？
--------------------------------------------------------------------------- SystemError Traceback (most recent call last) Cell In[5], line 1 ----> 1 get_ipython().run_cell_magic('stata', '', 'teffects nnmatch (lninc) (college),ematch(hukou age gender race sibling fmedu)\n') File c:\Users\asus\AppData\Local\Programs\Python\Python312\Lib\site-packages\IPython\core\interactiveshell.py:2543, in InteractiveShell.run_cell_magic(self, magic_name, line, cell) 2541 with self.builtin_trap: 2542 args = (magic_arg_s, cell) -> 2543 result = fn(*args, **kwargs) 2545 # The code below prevents the output from being displayed 2546 # when using magics with decorator @output_can_be_silenced 2547 # when the last Python token in the expression is a ';'. 2548 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False): File C:\Program Files/Stata18\utilities\pystata\ipython\stpymagic.py:276, in PyStataMagic.stata(self, line, cell, local_ns) 274 _stata.run(cell, quietly=True, inline=_config.stconfig['grshow']) 275 else: --> 276 _stata.run(cell, quietly=False, inline=_config.stconfig['grshow']) 278 if '-gw' in args or '-gh' in args: 279 _config.set_graph_size(gwidth, gheight) File C:\Program Files/Stata18\utilities\pystata\stata.py:313, in run(cmd, quietly, echo, inline) 311 _stata_wrk1("qui " + cmds[0], echo) 312 else: --> 313 _stata_wrk1(cmds[0], echo) 314 else: 315 if inline: File C:\Program Files/Stata18\utilities\pystata\stata.py:71, in _stata_wrk1(cmd, echo) 69 err = callback[0] 70 callback.clear() ---> 71 raise SystemError(err) 72 except KeyboardInterrupt: 73 outputter.done() SystemError: no exact matches for observation 4; use option osample() to identify all observations with deficient matches r(459);
从结果的输出看，Stata 给出了一个错误提示，提示内容无法为所有个案找到在所有协变量上取值都相同的匹配对象。
于此同时，Stata 还给出了一个建议，即使用 option osample() 去标记那些匹配失败的个案。
参照这个建议，我们在增加选项 osample() 后重新执行了上述命令。输出结果显示，有 593 名个案无法按照我们设定的要求实施精确匹配。
对标识变量 overlap 进行统计描述也可以发现，有 593 名个案违反了共同取值范围假定。
--------------------------------------------------------------------------- SystemError Traceback (most recent call last) Cell In[6], line 1 ----> 1 get_ipython().run_cell_magic('stata', '', 'teffects nnmatch (lninc) (college),ematch(hukou age gender race sibling fmedu) osample(overlap)\n') File c:\Users\asus\AppData\Local\Programs\Python\Python312\Lib\site-packages\IPython\core\interactiveshell.py:2543, in InteractiveShell.run_cell_magic(self, magic_name, line, cell) 2541 with self.builtin_trap: 2542 args = (magic_arg_s, cell) -> 2543 result = fn(*args, **kwargs) 2545 # The code below prevents the output from being displayed 2546 # when using magics with decorator @output_can_be_silenced 2547 # when the last Python token in the expression is a ';'. 2548 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False): File C:\Program Files/Stata18\utilities\pystata\ipython\stpymagic.py:276, in PyStataMagic.stata(self, line, cell, local_ns) 274 _stata.run(cell, quietly=True, inline=_config.stconfig['grshow']) 275 else: --> 276 _stata.run(cell, quietly=False, inline=_config.stconfig['grshow']) 278 if '-gw' in args or '-gh' in args: 279 _config.set_graph_size(gwidth, gheight) File C:\Program Files/Stata18\utilities\pystata\stata.py:313, in run(cmd, quietly, echo, inline) 311 _stata_wrk1("qui " + cmds[0], echo) 312 else: --> 313 _stata_wrk1(cmds[0], echo) 314 else: 315 if inline: File C:\Program Files/Stata18\utilities\pystata\stata.py:71, in _stata_wrk1(cmd, echo) 69 err = callback[0] 70 callback.clear() ---> 71 raise SystemError(err) 72 except KeyboardInterrupt: 73 outputter.done() SystemError: 593 observations have no exact matches; they are identified in the osample() variable r(459);
overlap |
violation |
indicator | Freq. Percent Cum.
------------+-----------------------------------
0 | 3,544 85.67 85.67
1 | 593 14.33 100.00
------------+-----------------------------------
Total | 4,137 100.00在同时对多个变量实施精确匹配，且变量中包含像 age 这样的连续变量的时候，精确匹配必然会遭遇维度诅咒，其最直接的表现就是匹配后的样本量大为下降。
为了避免损失过多样本，我们尝试对 age、race 和 sibling 这三个变量实施马氏匹配，其余变量依然采用精确匹配。具体如下：
Treatment-effects estimation Number of obs = 4,137
Estimator : nearest-neighbor matching Matches: requested = 1
Outcome model : matching min = 1
Distance metric: Mahalanobis max = 50
------------------------------------------------------------------------------
| AI robust
lninc | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATE |
college |
(是 vs 否) | .7699847 .041446 18.58 0.000 .6887521 .8512174
------------------------------------------------------------------------------可以发现，在对 age、race 和 sibling 这三个变量实施马氏匹配以后，软件不再提示有匹配失败的问题，且给出了 ATE 的估计值为 0.770.
上面演示的马氏匹配采用的是默认的 1 对 1 匹配，若要采用 1 对多匹配，需要使用选项 nneighbor()。具体来说，我们可以参照阿巴迪和因本斯的建议采用 1 对 4 匹配，同时使用 1 对 4 匹配下的稳健标准误，具体如下：
Treatment-effects estimation Number of obs = 4,137
Estimator : nearest-neighbor matching Matches: requested = 4
Outcome model : matching min = 4
Distance metric: Mahalanobis max = 50
------------------------------------------------------------------------------
| AI robust
lninc | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATE |
college |
(是 vs 否) | .7936157 .0392848 20.20 0.000 .7166189 .8706125
------------------------------------------------------------------------------可以发现，1 对 4 匹配时的 ATE 为 0.794，与 1 对 1 匹配是很接近，因此分析结果总体上来说是稳健的。
不过，马氏匹配并不能完全消除干预组和对照组在协变量上的差异。为了解决这个问题，研究者最好使用阿巴迪和因本斯提出的方法对估计结果进行偏差校正。
对 age、race 和 sibling 这三个变量实施偏差校正后的 1 对 4 匹配估计结果如下：
Treatment-effects estimation Number of obs = 4,137
Estimator : nearest-neighbor matching Matches: requested = 4
Outcome model : matching min = 4
Distance metric: Mahalanobis max = 50
------------------------------------------------------------------------------
| AI robust
lninc | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATE |
college |
(是 vs 否) | .7949957 .0392914 20.23 0.000 .717986 .8720054
------------------------------------------------------------------------------分析结果显示，偏差校正后的估计值为 0.795，与偏差校正前的结果几乎没有差异，这说明在这个例子中，马氏匹配的偏差很小，可以忽略不计。
最后，对上述命令使用选项 atet，即可得到偏差校正后得 ATT。分析结果显示，其估计值为 0.699，小于 ATE：
Treatment-effects estimation Number of obs = 4,137
Estimator : nearest-neighbor matching Matches: requested = 4
Outcome model : matching min = 4
Distance metric: Mahalanobis max = 50
------------------------------------------------------------------------------
| AI robust
lninc | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATET |
college |
(是 vs 否) | .6993424 .0475923 14.69 0.000 .6060633 .7926215
------------------------------------------------------------------------------接下来，我们将演示粗化精确匹配。在使用这种方法之前，我们可以先用命令 imb 计算干预组和控制组在所有协变量上的不平衡性：
. # 使用 imb 命令需要下载 cem 程序包，
Unknown #command
. ssc install cem,replace
checking cem consistency and verifying not already installed...
all files already exist and are up to date.
. imb hukou age gender race sibling fmedu,treatment(college)
Multivariate L1 distance: .42932057
Univariate imbalance:
L1 mean min 25% 50% 75% max
hukou .11753 .11753 0 0 0 0 0
age .32978 -6.5436 0 -6 -9 -9 0
gender .04055 -.04055 0 0 0 0 0
race .01072 -.01072 0 0 0 0 0
sibling .13284 .13284 0 0 0 0 0
fmedu .22349 .01915 0 0 1 -1 0
. 分析结果显示，干预组和控制组在所有 6 个协变量上的总体不平衡指数 \(L_1\) 为0.429 。此外，Stata 还汇报了干预组和控制组在每个协变量上的不平衡指数以及两组人在每个协变量的均值、最小值、25%分位数、50%分位数、75%分位数和最大值上的差异。
一般来说，在比较数据不平衡程度的时候，总体不平衡指数最重要，因此，接下来我们主要根据这个指标来评估匹配效果。
1 粗化精确匹配
现在，我们尝试对数据进行粗化精确匹配。
首先，我们以30、35、40、45和50为分割点，将年龄粗化为一个 6 分类变量；而其他变量则使用软件默认的方式进行粗化。
需要注意的是，fmedu 有 3 个取值：初中及以下、高中及以上和数据缺失。考虑到这 3 个类别不宜合并，因此，我们在其变量名后的括号中使用 #0 表示不对该变量进行粗化。
采用上述设置的命令和结果如下：
. local cut = "#0"
. cem hukou age(30 35 40 45 50) gender race sibling fmedu (`cut'),treatment(col
> lege)
Matching Summary:
-----------------
Number of strata: 206
Number of matched strata: 133
0 1
All 2494 1643
Matched 2409 1583
Unmatched 85 60
Multivariate L1 distance: .20630897
Univariate imbalance:
L1 mean min 25% 50% 75% max
hukou 8.3e-16 1.1e-16 0 0 0 0 0
age .07527 -.1677 0 0 1 0 0
gender 1.4e-15 6.7e-16 0 0 0 0 0
race 9.0e-16 1.8e-15 0 0 0 0 0
sibling 2.4e-15 -1.9e-16 0 0 0 0 0
fmedu 3.7e-16 2.8e-15 0 0 0 0 0
. 输出结果显示，按照上述方法进行粗化精确匹配之后，大多数干预组个案和控制组个案找到了与之精确匹配的对象，匹配失败的个案之后 145 个。其中干预组个案60个，控制组个案85个。
匹配之后，干预组和控制组在所有6个协变量上的总体不平衡指数下降到了 0.206，且单个协变量的不平衡指数都非常接近 0。
这说明与匹配之前相比，匹配后的样本在平衡性上有很大的提升。
以上所示的命令将 age 粗化为一个等间距的 6 分类变量。此外，我们也可以根据分位数将 age 粗化为等规模的 6 分类变量。具体命令和结果如下：
(using the scott break method for imbalance)
Matching Summary:
-----------------
Number of strata: 179
Number of matched strata: 118
0 1
All 2494 1643
Matched 2421 1579
Unmatched 73 64
Multivariate L1 distance: .22223609
Univariate imbalance:
L1 mean min 25% 50% 75% max
hukou 5.3e-15 4.6e-15 0 0 0 0 0
age .08636 -.24076 0 0 0 0 0
gender 6.0e-15 7.2e-15 0 0 0 0 0
race 1.2e-15 2.1e-15 0 0 0 0 0
sibling 4.1e-15 1.7e-15 0 0 0 0 0
fmedu 7.4e-15 5.3e-15 0 0 0 0 0分析结果显示，通过上述方法保留的匹配样本略多于前一种方法，但平衡性指数却比之前高。所以综合来看，匹配效果与之前相当。
最后，我们演示根据 cem 的自动粗化算法得到的匹配结果：
(using the scott break method for imbalance)
Matching Summary:
-----------------
Number of strata: 369
Number of matched strata: 214
0 1
All 2494 1643
Matched 2330 1523
Unmatched 164 120
Multivariate L1 distance: .09076546
Univariate imbalance:
L1 mean min 25% 50% 75% max
hukou 1.1e-15 -6.7e-16 0 0 0 0 0
age .02518 -.02531 0 0 0 0 0
gender 1.8e-15 -2.3e-15 0 0 0 0 0
race 8.7e-18 0 0 0 0 0 0
sibling 1.3e-15 -3.1e-16 0 0 0 0 0
fmedu 1.6e-15 1.2e-15 0 0 0 0 0可以发现，使用这种方法损失的样本量很多，但不平衡指数却比之前两种方法有大幅下降。
这体现了在粗化精确匹配过程中经常出现的一个现象，即排除在外的极端样本越多，匹配样本的平衡性通常就越能够得到保证。考虑到第三章方法的平衡性指数最小，且样本损失依然在可接受的范围内，我们接下来将用该方法生成的权重 cem_weights 进行后续分析
2 使用权重进行回归分析
首先，我们可以使用该权重进行加权后的一元线性回归：
Source | SS df MS Number of obs = 3,852
-------------+---------------------------------- F(1, 3850) = 333.41
Model | 450.866614 1 450.866614 Prob > F = 0.0000
Residual | 5207.69062 3,850 1.35264691 R-squared = 0.0797
-------------+---------------------------------- Adj R-squared = 0.0797
Total | 5658.55723 3,851 1.46937347 Root MSE = 1.1629
------------------------------------------------------------------------------
lninc | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
college | .6996738 .0383184 18.26 0.000 .6245475 .7748001
_cons | 9.47577 .0240912 393.33 0.000 9.428537 9.523003
------------------------------------------------------------------------------分析结果显示，通过这种方法得到的 college 的回归系数为 0.700，或者说 ATT（干预组的平均干预效应）为 0.700.这与之前偏差校正以后得到的 1 对 4 马氏匹配的估计结果非常接近。
接下来，我们可以在回归方程中纳入原始的协变量，以进一步消除匹配样本在这些变量上的不平衡性。
分析结果显示，纳入更多控制变量以后，college 的回归系数只发生了非常细微的变化。这主要是因为，粗化精确匹配已经在很大程度上消除了干预组和控制组在原始协变量上的不平衡性（匹配后 \(L_1\) 从 0.429 下降到了 0.091），此时，在回归方程中纳入更多控制变量的意义已经不大。
不过在统计分析时控制这些协变量是一种更加稳健的做法，因此，这里建议研究者在回归分析时始终对原始协变量进行统计控制。
Source | SS df MS Number of obs = 3,853
-------------+---------------------------------- F(8, 3844) = 59.25
Model | 621.146148 8 77.6432685 Prob > F = 0.0000
Residual | 5037.41109 3,844 1.31046074 R-squared = 0.1098
-------------+---------------------------------- Adj R-squared = 0.1079
Total | 5658.55723 3,852 1.46899201 Root MSE = 1.1448
------------------------------------------------------------------------------
lninc | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
college | .699765 .0377211 18.55 0.000 .6258097 .7737203
hukou | .1384669 .0417004 3.32 0.001 .0567099 .2202238
age | .0036032 .002706 1.33 0.183 -.0017021 .0089086
gender | .3631285 .0378722 9.59 0.000 .2888771 .4373799
race | .1562387 .1279194 1.22 0.222 -.0945576 .4070351
sibling | .1464629 .0537444 2.73 0.006 .0410927 .2518332
|
fmedu |
是 | .0025106 .043347 0.06 0.954 -.0824748 .087496
缺失 | .0241171 .0565652 0.43 0.670 -.0867836 .1350179
|
_cons | 8.892506 .162912 54.58 0.000 8.573104 9.211908
------------------------------------------------------------------------------