Stata 手动：各类匹配方法大全 A——理论篇

发布时间：2020-11-09 阅读 11153

Stata手动：各类匹配方法大全 A——理论篇

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者：黄俊凯 (中国人民大学)
E-Mail: kopanswer@126.com

匹配是研究处理效应的常见工具，本文总结了常见的匹配方法，并在第三部分给出 Stata 模拟以比较不同的匹配方法的优劣。

1. 单变量匹配 uni-variate match

单变量匹配的方法有精确匹配、粗糙精确匹配，k-近邻匹配和半径 (卡尺) 匹配。

1.1 精确匹配 exact match

顾名思义，当且仅当两个观测值的匹配变量相等时匹配成功。

1.2 粗糙精确匹配 coarsened exact match

粗糙精确匹配用途广泛，通常要求匹配变量是分类变量。比如公司金融中同行业的公司，又比如教育经济学中在同一个班的同学。粗糙精确匹配可以轻松的推广到多变量匹配的情形，如同行业同年度的公司，同班同性别的同学。

1.3 k-近邻匹配 k-nearest neighbor match

k-近邻匹配要求匹配变量是距离，它选取距离最近的 k 个观测值作为对照组。

1.4 radius (caliper) match

k-近邻匹配要求匹配变量也是距离，它事先设定半径 (上下半径可以不同)，找出设定范围内的全部观测值作为对照组。显然，随着半径的降低，匹配要求也更趋严格。

2. 多变量匹配 multi-variate match

多变量匹配的核心思路是降维 (dimension reduction)，将多变量降维为距离或得分，然后再运用单变量匹配的方法。多变量匹配的方法有欧氏距离、百分等级、马氏距离和倾向得分匹配等。

为能直观的图示，下文全部示例仅考虑双变量 (平面) 的情形。

2.1 欧氏距离匹配 euclidean distance match

最简单的测度空间就是欧氏距离空间。你可以轻松的用尺子直接量出点到原点、或任意两点之间的距离。

平面上任意点 $X = (X_{1}, X_{2}, \dots, X_{n})$ 到原点的距离公式为：

d_{e} (X) = \sqrt{X_{1}^{2} + X_{2}^{2} + \dots + X_{n}^{2}}

平面上任意两点 $X = (X_{1}, X_{2}, \dots, X_{n})$ 和 $Y = (Y_{1}, Y_{2}, \dots, Y_{n})$ 之间的距离公式为：

d_{e} (X, Y) = \sqrt{(X_{1} - Y_{1})^{2} + (X_{2} - Y_{2})^{2} + \dots + (X_{n} - Y_{n})^{2}}

在欧式距离空间中，坐标轴之间是相互垂直的，也就意味着随机变量之间的相关系数为零。这是一个非常强的假设，而实际研究中该假设往往不能成立。比如研究一个人的健康时，作为影响因素的身高和体重之间往往是正相关的。

在欧式距离空间中，坐标轴上的量纲也是默认相同的，也就意味着随机变量必须有相同的量纲。这同样是一个过于严苛的假设。同样是研究一个人的健康，身高的单位是厘米，体重的单位是公斤，无论如何长度单位和重量单位是不可比的。

2.2 马氏距离匹配 mahalanobis distance match

为了克服匹配时欧氏距离空间的上述两个缺陷，印度数学家 Mahalanobis 提出了著名的马氏距离空间。马氏距离对多匹配变量的联合分布的位置、形状做了标准化调整，也对多匹配变量之间的相关性做了正交化调整，从而将多匹配变量的联合分布转变为可以用欧氏距离计算的情形。

马氏距离实现上述功能的关键在于用协方差矩阵的逆矩阵做了调整。我们将在第三部分用 Stata 代码逐步展示它。任意点 $\vec{x}$ 到质心 (centroid) $\vec{c}$ 的距离公式如下：

d_{m} (\vec{x}) = \sqrt{(\vec{x} - \vec{c})^{'} {C o v}^{- 1} (\vec{x} - \vec{c})}

任意两点 ${\vec{x}}_{1}$ 和 ${\vec{x}}_{2}$ 之间的距离公式如下：

d_{m} ({\vec{x}}_{1}, {\vec{x}}_{2}) = \sqrt{({\vec{x}}_{1} - {\vec{x}}_{2})^{'} {C o v}^{- 1} ({\vec{x}}_{1} - {\vec{x}}_{2})}

其中 $C o v$ 为多维随机变量的样本协方差矩阵， $c$ 是质心或样本均值。若样本协方差矩阵是单位矩阵 (各维度之间独立同部分)，马氏距离退化为欧氏距离。

马氏距离测量相对于质心的距离，质心是一个基准点或中心点，可以认为是多元数据的总体平均值。质心是多元空间中所有变量的均值相交的点。马氏距离越大，数据点离质心越远。

2.3 百分等级匹配 percentile rank match

在多随机变量的联合分布中，距离近的点在概率分布函数 (PDF) 上往往不一定靠近。我们以如下标准正态分布 (钟形曲线) 为例,说明距离相近并不必然等于概率分布上的靠近。

不妨设样本容量为 $10000$ 。处理组观测值点 A 位于距离原点 $- 2 σ$ 处。则 $- σ$ 的点 B 和 $- 3 σ$ 处的点 C 与点 A 的距离相等，都等于 $- σ$ 。但是点 A 与点 B 之间有 $1359$ ( $13.59 % \times 10000$ ) 个观测值，而点 A 与点 C 之间仅有 $214$ ( $2.14 % \times 10000$ ) 个观测值。显然，尽管点 A 到点 B 和点 C 的距离是相等的，但从概率分布的角度来看在点 A 前后抽样更容易抽到点 C。

百分等级是一个相对位置量数。它是指将距离按照从小到大的顺序排列，小于某观测值距离的观测值个数与样本数的百分比，即为该观测值的百分等级。

不同于百分数 (percentile)，有相同距离的观测值总是有相同的百分等级 (percentile rank)，以消除指定百分等级时的武断性 (arbitrariness)。因此，百分等级的公式如下：

P R = \frac{L + 0.5 E}{N}

其中，PR 是百分等级，L 是距离小于该观测值距离的观测值数，E 是距离等于该观测值距离的观测值数，N 是样本数。例如，样本有 $100$ 个同学参加考试，你的物理得分是 $85$ 分，其中有 $60$ 个人分数低于你， $5$ 个人的分数和你相同，则你的百分等级数就等于：

\frac{60 + 0.5 \times 5}{100} = 62.5 %

百分等级匹配就是选择百分等级最接近的观测值作为反事实观测值，显然它是一种非参数方法。基于距离的百分等级不仅克服了距离不能体现随机变量分布函数的特点 (比如距离很远的点也可能有相似的百分等级)，而且将匹配变量之间的差异限制在 $[0, 1]$ 区间。然而它仍有自己的局限，就是无法体现个体进入处理组的选择过程，或者说并没有利用到底哪些个体进入了处理组这个关键信息。

2.4 倾向得分匹配 propensity score match

倾向得分匹配不是基于距离而是基于得分。倾向得分 (propensity score) 是指个体 i 在给定可观测变量 $x_{i}$ 的情况下进入处理组的条件概率。即 $p (x_{i}) = P r (D_{i} = 1 | x = x_{i})$ ，或简记为 $p (x)$ 。

对倾向得分的估计可使用参数估计或非参数估计，最流行的方法是 logit。Rosenbaum and Rubin (1983) 证明，如果可忽略性假定成立，则在给定倾向得分的情况下，结果变量在事前和事后的取值独立于是否个体进入处理组。

使用倾向得分匹配需要满足两个假定：1.共同支撑假设；2.平衡性假定。

共同支撑集是指这样一个集合，不妨设备择组观测值的倾向得分的取值范围为 $[p_{m i n}^{u}, p_{m a x}^{u}]$ , 设处理组观测值得倾向得分取值范围为 $[p_{m i n}^{t}, p_{m a x}^{t}]$ ，则共同支撑集内的任意观测值的倾向得分必需大于备择组和处理组的最小倾向得分中较大的那个，也必需小于最大倾向得分中较小的那个。即进入共同支撑集的观测值必需位于如下区间内：

[max (p_{m i n}^{u}, p_{m i n}^{t}), min (p_{m a x}^{u}, p_{m a x}^{t})]

其中 $p_{m i n}^{u}$ 是对照组倾向得分的最小值， $p_{m a x}^{u}$ 是对照组倾向得分的最大值； $p_{m i n}^{t}$ 是处理组倾向得分的最小值， $p_{m a x}^{t}$ 是处理组倾向得分的最大值。

在匹配时，为提供匹配质量，可仅保留倾向得分重叠部分的个体。共同支撑假设要求共同支撑集的取值范围不能太小，则匹配成功的样本可能不具有代表性。

再说平衡性假定，如果倾向得分估计是准确的，那么协变量 $x_{i}$ 在匹配后的处理组与控制组之间分布较均匀。例如，匹配后处理组的均值 $\bar{x^{t}}$ 与对照组的均值 $\bar{x^{u}}$ 应该较接近，也就是所谓的“数据平衡”。

多接近是接近？均值之间的差异与计量单位有关，你当然可以做组建军之差异或中位数差异检验，也可以使用标准化偏差 (standarized bias) 进行比较，公式如下：

\frac{{\bar{x}}^{t} - {\bar{x}}^{u}}{\sqrt{(s_{x}^{t})^{2} + (s_{x}^{u})^{2} / 2}}

经验法则告诉我们，标准化差距不得超过 10%。如果超过，则应重新估计倾向得分，或改变具体的匹配方法。

尽管倾向得分匹配体现了选择过程，但它仍有一些缺陷，比如它只能依可测变量选择，对共同支撑集有要求等。

3. 代码展示

我们用简单的二维模拟数据来展示不同多变量匹配方法之间的差异。

3.1 数据生成过程

clear
set obs 200
gen id = _n
set seed 10000
gen r0 = rnormal()
set seed 12345
gen r1 = rnormal()
set seed 65432
gen r2 = rnormal()
set seed 10101
gen re = rnormal()

*-选择变量
gen x1 = 2 * (r0 + r1)
gen x2 = 3 + (r0 + r2)

*-选择机制
gen group = 2 * x1 - 3 * x2 + 10 * re > 0		
label def group 0 "非对照组" 1 "处理组"
label val group group

*-数据分布
keep id x1 x2 group
tab group
correlate x1 x2
twoway (scatter x1 x2 if group == 0, mcolor(black))	///
	   (scatter x1 x2 if group == 1, mcolor(red)), 	///
	   title("数据分布") xlabel(-10(5)10) ///
       ylabel(-10(5)10, nogrid) ///
       aspect(1) legend(off)

下图种，红色的点为处理组观测值，黑色的点为对照组观测值

3.2 欧氏距离匹配的结果

变量 ec 记录欧氏距离匹配得到的匹配值

gen ec = .
gen ed = .
forvalues i = 1(1)200 {
  if group[`i'] == 1 {
    forvalues j = 1(1)200 {
      if group[`j'] == 0 {
        local d = sqrt((x1[`i'] - x1[`j'])^2 + (x2[`i'] - x2[`j'])^2)
        if `d' < ed[`i'] {
          qui replace ec = `j' in `i'
          qui replace ed = `d' in `i'
        }
      }
    }
  }
}
tab group
misstable sum ec ed

欧氏距离匹配的结果：

3.3 马氏距离匹配的结果

变量 mc 记录马氏距离匹配得到的匹配值

corr x1 x2, cov
mat cov = r(C)
mat cov = inv(cov)
mat list cov
scalar a_11 = cov[1, 1]
scalar a_22 = cov[2, 2]
scalar a_21 = cov[2, 1]
gen mc = .
gen md = .
forvalues i = 1(1)200 {
  if group[`i'] == 1 {
    forvalues j = 1(1)200 {
      if group[`j'] == 0 {
        local d = sqrt(a_11*(x1[`i'] - x1[`j'])^2 + a_22*(x2[`i'] - x2[`j'])^2 + 2 * a_21 * (x1[`i'] - x1[`j']))
        if `d' < md[`i'] {
          qui replace mc = `j' in `i'
          qui replace md = `d' in `i'
        }
      }
    }
  }
}

马氏距离匹配的结果：

3.4 倾向得分匹配的结果

变量 pc 记录倾向得分匹配得到的匹配值

gen pc = .
gen pd = .	
probit group x1 x2
predict score, xb
replace score = normal(score)
forvalues i = 1(1)200 {
  if group[`i'] == 1 {
    forvalues j = 1(1)200 {
      if group[`j'] == 0 {
        local d = abs(score[`i'] - score[`j'])
        if `d' < pd[`i'] {
          qui replace pc = `j' in `i'
          qui replace pd = `d' in `i'
        }
      }
    }
  }
}

倾向得分匹配的结果：

4. 参考文献和资料

Mahalanobis Distance: Simple Definition, Examples link
Mahalanobis, P. C. (1936). On the generalized distance in statistics. National Institute of Science of India. pdf
Percentiles, Percentile Rank & Percentile Range: Definition & Examples link
Rosenbaum, P. R., & Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1), 41-55. pdf

附：文中使用的 dofiles

*-3.1 数据生成过程

clear
set obs 200
gen id = _n
set seed 10000
gen r0 = rnormal()
set seed 12345
gen r1 = rnormal()
set seed 65432
gen r2 = rnormal()
set seed 10101
gen re = rnormal()

*-选择变量
gen x1 = 2 * (r0 + r1)
gen x2 = 3 + (r0 + r2)

*-选择机制
gen group = 2 * x1 - 3 * x2 + 10 * re > 0		
label def group 0 "非对照组" 1 "处理组"
label val group group

*-数据分布
keep id x1 x2 group
tab group
correlate x1 x2
twoway (scatter x1 x2 if group == 0, mcolor(black))	///
	   (scatter x1 x2 if group == 1, mcolor(red)), 	///
	   title("数据分布") xlabel(-10(5)10) ///
       ylabel(-10(5)10, nogrid) ///
       aspect(1) legend(off)
	   
	   
*-3.2 欧氏距离匹配的结果

gen ec = .
gen ed = .
forvalues i = 1(1)200 {
  if group[`i'] == 1 {
    forvalues j = 1(1)200 {
      if group[`j'] == 0 {
        local d = sqrt((x1[`i'] - x1[`j'])^2 + (x2[`i'] - x2[`j'])^2)
        if `d' < ed[`i'] {
          qui replace ec = `j' in `i'
          qui replace ed = `d' in `i'
        }
      }
    }
  }
}
tab group
misstable sum ec ed	   


*-3.3 马氏距离匹配的结果

corr x1 x2, cov
mat cov = r(C)
mat cov = inv(cov)
mat list cov
scalar a_11 = cov[1, 1]
scalar a_22 = cov[2, 2]
scalar a_21 = cov[2, 1]
gen mc = .
gen md = .
forvalues i = 1(1)200 {
  if group[`i'] == 1 {
    forvalues j = 1(1)200 {
      if group[`j'] == 0 {
        local d = sqrt(a_11*(x1[`i'] - x1[`j'])^2 + a_22*(x2[`i'] - x2[`j'])^2 + 2 * a_21 * (x1[`i'] - x1[`j']))
        if `d' < md[`i'] {
          qui replace mc = `j' in `i'
          qui replace md = `d' in `i'
        }
      }
    }
  }
}  


*-3.4 倾向得分匹配的结果

gen pc = .
gen pd = .	
probit group x1 x2
predict score, xb
replace score = normal(score)
forvalues i = 1(1)200 {
  if group[`i'] == 1 {
    forvalues j = 1(1)200 {
      if group[`j'] == 0 {
        local d = abs(score[`i'] - score[`j'])
        if `d' < pd[`i'] {
          qui replace pc = `j' in `i'
          qui replace pd = `d' in `i'
        }
      }
    }
  }
}

专题	嘉宾	直播/回看视频
⭐ 最新专题 ⭐		DSGE, 因果推断, 空间计量等
⭕ Stata数据清洗	游万海	直播, 2 小时，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]