【因果推断】一文读懂倾向匹配得分Stata及R操作应用
本文将介绍Stata及R软件进行倾向匹配得分操作应用，主要包括倾向匹配得分命令简介、语法格式、倾向匹配得分操作步骤 思路，涉及倾向匹配得分应用、平衡性检验、共同取值范围检验、核密度函数图等内容。
1
命令简介
Stata does not have a built-in command for propensity score matching, a non-experimental method of sampling that produces a control group whose distribution of covariates is similar to that of the treated group. However, there are several user-written modules for this method. The following modules are among the most popular:
Stata没有一个内置的倾向评分匹配的命令，一种非实验性的抽样方法，它产生一个控制组，它的协变量分布与被处理组的分布相似。但是，这个方法有几个用户编写的模块。以下是最受欢迎的模块（主要有如下几个外部命令）
psmatch2.ado
pscore.ado
nnmatch.ado
psmatch2.ado was developed by Leuven and Sianesi (2003) and pscore.ado by Becker and Ichino (2002). More recently, Abadie, Drukker, Herr, and Imbens (2004) introduced nnmatch.ado. All three modules support pair-matching as well as subclassification.
You can find these modules using the .net command as follows:
net search psmatch2
net search pscore
net search nnmatch
You can install these modules using the .ssc or .net command, for example:
ssc install psmatch2, replace
After installation, read the help files to find the correct usage, for example:
help psmatch2
上述主要介绍了如何获得PSM相关的命令，总结一下目前市面上用的较好的命令为psmatch2.
PSM 相关命令
help psmatch2
help nnmatch
help psmatch
help pscore
持续获取最新的 PSM 信息和程序
findit propensity score
findit matching
psmatch2 is being continuously improved and developed. Make sure to keep your version up-to-date as follows
ssc install psmatch2, replace
where you can check your version as follows:
which psmatch2
2
语法格式
语法格式为：
help psmatch2
psmatch2 depvar [indepvars] [if exp] [in range] [, outcome(varlist) pscore(varname) neighbor(integer) radius caliper(real) mahalanobis(varlist) ai(integer) population altvariance kernel llr kerneltype(type) bwidth(real) spline nknots(integer) common trim(real) noreplacement descending odds index logit ties quietly w(matrix) ate]选项含义为：
depvar因变量 ；
indepvars表示协变量；
outcome(varlist)表示结果变量；
logit指定使用logit模型进行拟合，默认的是probit模型；
neighbor(1)指定按照1:1进行匹配，如果要按照1:3进行匹配，则设定为neighbor(3)；
radius表示半径匹配
核匹配 (Kernel matching)
其他匹配方法
广义精确匹配(Coarsened Exact Matching) || help cem
局部线性回归匹配 (Local linear regression matching)
样条匹配 (Spline matching)
马氏匹配 (Mahalanobis matching)
pstest $X, both做匹配前后的均衡性检验，理论上说此处只能对连续变量做均衡性检验，对分类变量的均衡性检验应该重新整理数据后运用χ2检验或者秩和检验。但此处对于分类变量也有一定的参考价值。
psgraph对匹配的结果进行图示。
3
Stata操作与应用
政策背景：国家支持工作示范项目（ National Supported Work,NSW ）
研究目的：检验接受该项目（培训）与不接受该项目（培训）对工资的影响。基本思想：分析接受培训组（处理组， treatment group ）接受培训行为与不接受培训行为在工资表现上的差异。但是，现实可以观测到的是处理组接受培训的事实，而处理组没有接受培训会怎样是不可能观测到的，这种状态也成为反事实（ counterfactual ）。
匹配法就是为了解决这种不可观测事实的方法。在倾向得分匹配方法（ Propensity Score Matching ）中，根据处理指示变量将样本分为两个 组，一是处理组，在本例中就是在 NSW 实施后接受培训的组；二是对照组 （ comparison group ），在本例中就是在 NSW 实施后不接受培训的组。倾向得分 匹配方法的基本思想是，在处理组和对照组样本通过一定的方式匹配后，在其他 条件完全相同的情况下，通过接受培训的组（处理组）与不接受培训的组（对照组）在工资表现上的差异来判断接受培训的行为与工资之间的因果关系。
1、首先进行数据结构查看
use "ldw_exper.dta", clear eddesc结果为：
2、描述性分析
tabulate t, summarize(re78) means standard结果为：
3、倾向匹配得分
3.1 首先进行排序，生成随机数种子
set seed 20180105 //产生随机数种子 gen u=runiform() sort u //排序 或者 order u3.2 倾向匹配得分
local v1 "t"local v2 "age edu black hisp married re74 re75 u74 u75"global x "`v1' `v2' "
psmatch2 $x, out(re78) neighbor(1) ate ties logit common // 1:1 匹配$表示引用宏变量，
等价于 psmatch2 t age edu black hisp married re74 re75 u74 u75, out(re78) neighbor(1) ate ties logit common结果为：
3.3 查看匹配后数据
结果为：
打开数据编辑窗口，会发现软件自动生成了几个新变量：
其中_pscore是每个观测值对应的倾向值；
_id是自动生成的每一个观测对象唯一的ID（事实上这列变量即是对_pscore排序）；
_treated表示某个对象是否试验组；
_n1表示的是他被匹配到的对照对象的_id（如果是1:3匹配，还会生成_n2, _n3）；
_pdif表示一组匹配了的观察对象他们概率值的差。
3.4 均衡性检验
pstest $v2, both graph结果为：
3.5 共同取值范围
psgraph结果为：
3.6 核密度函数图
twoway(kdensity _ps if _treat==1,legend(label(1 "Treat")))(kdensity _ps if _treat==0, legend(label(2 "Control"))),xtitle(Pscore> ) title("Before Matching")
. twoway(kdensity _ps if _treat==1,legend(label(1 "Treat")))(kdensity _ps if (_weight!=1&_weight!=.), legend(label(2 "Control"))),> xtitle(Pscore) title("After Matching")结果为：
3
R操作与应用
描述：
这是国家支持工作示范(NSW)处理组数据的子样本和当前人口调查(CPS)的比较样本。Lalonde(1986)、Dehejia和Wahba(1999)对这些数据进行了广泛的分析。
R软件操作倾向匹配得分的函数安装包比较多，下面介绍一种，叫做Matching，其他的还有MatchIt等。
使用Matching前，需要先进行安装：install.packages(“Matching”)就可以安装了。
1、首先进行数据查看
> data("lalonde")> summary(lalonde) age educ black Min. :17.00 Min. : 3.0 Min. :0.0000 1st Qu.:20.00 1st Qu.: 9.0 1st Qu.:1.0000 Median :24.00 Median :10.0 Median :1.0000 Mean :25.37 Mean :10.2 Mean :0.8337 3rd Qu.:28.00 3rd Qu.:11.0 3rd Qu.:1.0000 Max. :55.00 Max. :16.0 Max. :1.0000 hisp married Min. :0.00000 Min. :0.0000 1st Qu.:0.00000 1st Qu.:0.0000 Median :0.00000 Median :0.0000 Mean :0.08764 Mean :0.1685 3rd Qu.:0.00000 3rd Qu.:0.0000 Max. :1.00000 Max. :1.0000 nodegr re74 Min. :0.000 Min. : 0.0 1st Qu.:1.000 1st Qu.: 0.0 Median :1.000 Median : 0.0 Mean :0.782 Mean : 2102.3 3rd Qu.:1.000 3rd Qu.: 824.4 Max. :1.000 Max. :39570.7 re75 re78 u74 Min. : 0 Min. : 0 Min. :0.0000 1st Qu.: 0 1st Qu.: 0 1st Qu.:0.0000 Median : 0 Median : 3702 Median :1.0000 Mean : 1377 Mean : 5301 Mean :0.7326 3rd Qu.: 1221 3rd Qu.: 8125 3rd Qu.:1.0000 Max. :25142 Max. :60308 Max. :1.0000 u75 treat Min. :0.0000 Min. :0.0000 1st Qu.:0.0000 1st Qu.:0.0000 Median :1.0000 Median :0.0000 Mean :0.6494 Mean :0.4157 3rd Qu.:1.0000 3rd Qu.:1.0000 Max. :1.0000 Max. :1.0000 > View(lalonde)结果为：
数据介绍：包含对以下12个变量的445个观察值。
age表示年龄
educ：受教育年限
black：是否为黑人
hispan：是否为西班牙裔
married：是否已婚
nodegree：是否没有毕业文凭。
re74，1974年的实际收入。
re75，1975年的实际收入。
re78，1978年的实际收入。
u74，1974年的收入指标变量为零。
u75，1975年收入的指标变量为零。
treat：是否接受职业培训。
2、估计倾向匹配模型
glm1 <- glm(treat~age + I(age^2) + educ + I(educ^2) + black + hisp + married + nodegr + re74 + I(re74^2) + re75 + I(re75^2) + u74 + u75, family=binomial, data=lalonde)3、保存数据对象
X <- glm1$fittedY <- lalonde$re78Tr <- lalonde$treat3、一对一匹配，评估处理效应的影响(“estimand”选项默认为ATT)。
rr <- Match(Y=Y, Tr=Tr, X=X, M=1);summary(rr)结果为：
4、我们来看看协变量的平衡性，为了提高速度，nboots被设置为较小的值。正式分析时提高到至少500。
mb <- MatchBalance(treat~age + I(age^2) + educ + I(educ^2) + black + hisp + married + nodegr + re74 + I(re74^2) + re75 + I(re75^2) + u74 + u75, data=lalonde, match.out=rr, nboots=10)结果为：