vlambda博客
学习文章列表

一、倾向性评分匹配（Propensity Score Matching）

二、数据集介绍

同样选用lalonde数据集，如前一篇文章介绍，该数据集收集了614个观察对象的十个特征：

1. age：年龄，以岁为单位。
2. educ：学习时间，之前一共上过多少年学。
3. black：是否为黑人。
4. hispan：是否为西班牙裔。
5. married：是否已婚。
6. nodegree：是否没有毕业文凭。
7. re74：1974年的收入。
8. re75：1975年的收入。
9. re78：1978年的收入。
10. treat：1978年前是否接受职业培训。

研究假说是接受过职业培训会增加1978年的收入水平。

三、R语言实现

加载包与数据集

#MatchIt包做PSM，table1做表
library(MatchIt)
library(table1)
#载入数据集
data(lalonde)
#查看变量数据类型
#str(lalonde)

数据集处理

#转换变量数据类型
lalonde$black <- factor(lalonde$black)
lalonde$hispan <- factor(lalonde$hispan)
lalonde$married <- factor(lalonde$married)
lalonde$nodegree <- factor(lalonde$nodegree)
lalonde$black <- as.logical(lalonde$black == 1)
lalonde$hispan <- as.logical(lalonde$hispan == 1)
lalonde$married <- as.logical(lalonde$married == 1)
lalonde$nodegree <- as.logical(lalonde$nodegree == 1)
#更改变量标注
label(lalonde$black) <- "黑人"
label(lalonde$educ) <- "上学时间"
label(lalonde$hispan) <- "西班牙裔"
label(lalonde$married) <- "已婚"
label(lalonde$nodegree) <- "无毕业文凭"
label(lalonde$age) <- "年龄"
label(lalonde$re74) <- "1974年收入"
label(lalonde$re75) <- "1975年收入"
label(lalonde$re78) <- "1978年收入"
units(lalonde$age) <- "岁"
units(lalonde$educ) <- "年"

输出基线特征表格

#为了符合习惯用法，把table1函数修改为表示mean±sd和p值，并重新封装
func_table1<-function(df,ycol,xcol,xlabels=c("Control", "Treatment")){
 if( any(c(missing(df),missing(xcol),missing(ycol)))){
 return("parameters missing! arg1 data.frame, arg2 ycol, arg3 xcol")
 }
 else if(!all(c(is.data.frame(df),is.character(xcol),is.character(ycol)))==TRUE){
 return("parameter type error! arg1 data.frame, arg2,arg3 character vector")
 }
 df[[ycol]]<-as.numeric(df[[ycol]])
 df[[ycol]]<-factor(df[[ycol]],levels=c(0, 1, 2), labels=c(xlabels, "P-value"))
 
 rndr <- function(x, name, ...) {
 if (length(x) == 0) {
 y <- df[[name]]
 s <- rep("", length(render.default(x=y, name=name, ...)))
 if (is.numeric(y)) {
 p <- t.test(y ~ df[[ycol]])$p.value
 } else {
 p <- chisq.test(table(y, droplevels(df[[ycol]])))$p.value
 }
 s[2] <- sub("<", "&lt;", format.pval(p, digits=3, eps=0.001))
 s
 } else {
 render.default(x=x, name=name, ...)
 
 }
 }
 
 rndr.strat <- function(label, n, ...) {
 ifelse(n==0, label, render.strat.default(label, n, ...))
 }

 my.render.cont=c(.="Mean &plusmn; SD",.="Median [Min, Max]")
 
 my.render.cat <- function(x) {
 c("", sapply(stats.default(x), function(y) with(y,
 sprintf("%d (%0.0f %%)", FREQ, PCT))))
 }
 
 table1(as.formula(paste("~",paste(xcol,collapse = "+"),"|",ycol)),
 data=df, droplevels=F, render=rndr, render.strat=rndr.strat,
 render.continuous=my.render.cont,
 render.categorical=my.render.cat,
 overall=F)
}

#按是否接受职业培训分组，显示组间基本信息
func_table1(lalonde,"treat",names(lalonde)[-1],c("未职业培训组","职业培训组"))

可以看到组间并没有1978年收入没有差异，但是其他混杂因素均有较显著差异，现在就用R来行倾向性分析，匹配后可得到剔除混杂因素影响的样本。

set.seed(50)
#PSM
m.out <- matchit(data = lalonde,
 #以treat分组，匹配混杂因素(age,educ,black,hispan,married,nodegree,re74,re75)无差异的样本
 formula = treat ~ age + educ + black + hispan + married + nodegree + re74 + re75
 ,
 method = "nearest",
 distance = "logit",
 replace = FALSE,
 caliper = 0.05)
#匹配后样本数据
lalonde_matched<-match.data(m.out)

#更改变量标注
label(lalonde_matched$black) <- "黑人"
label(lalonde_matched$educ) <- "上学时间"
label(lalonde_matched$hispan) <- "西班牙裔"
label(lalonde_matched$married) <- "已婚"
label(lalonde_matched$nodegree) <- "无毕业文凭"
label(lalonde_matched$age) <- "年龄"
label(lalonde_matched$re74) <- "1974年收入"
label(lalonde_matched$re75) <- "1975年收入"
label(lalonde_matched$re78) <- "1978年收入"
units(lalonde_matched$age) <- "岁"
units(lalonde_matched$educ) <-"年"
#显示匹配后数据
func_table1(lalonde_matched,"treat",names(lalonde)[-1],c("未职业培训组","职业培训组"))

经过倾向性评分匹配获得的数据组间其他混杂因素均无统计学差异，而职业培训组1978年收入较对照组明显升高，提示接受职业培训可以增加收入水平。此处计算结果与上一篇文章SPSS实现PSM的计算结果有所不同，本文R语言中采用了 “nearest”算法计算倾向值，SPSS中算法与此不同，所以大家在选择方法及统计软件的时候需标注清楚，另外本文R语言的计算结果也是值得商榷的，并非采用了R语言得到了与回归分析一致的结果就万事大吉了，采用PSM的方法仍然要谨慎，上一篇文章中提到的PSM适用原则仍然很有必要（）。

四、倾向性评分匹配的原理

真正做倾向性评分的代码其实只有一行

m.out <- matchit(data = lalonde,
 formula = treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,
 method = "nearest",
 distance = "logit",
 replace = FALSE,
 caliper = 0.05)

我们在看看匹配的数据是什么样子

head(lalonde_matched)

## treat age educ black hispan married nodegree re74 re75 re78
## NSW1 1 37 11 TRUE FALSE TRUE TRUE 0 0 9930.0460
## NSW2 1 22 9 FALSE TRUE FALSE TRUE 0 0 3595.8940
## NSW5 1 33 8 TRUE FALSE FALSE TRUE 0 0 289.7899
## NSW6 1 22 9 TRUE FALSE FALSE TRUE 0 0 4056.4940
## NSW7 1 23 12 TRUE FALSE FALSE FALSE 0 0 0.0000
## NSW10 1 33 12 FALSE FALSE TRUE FALSE 0 0 12418.0700
## distance weights
## NSW1 0.63876993 1
## NSW2 0.22463424 1
## NSW5 0.70163874 1
## NSW6 0.69906990 1
## NSW7 0.65368426 1
## NSW10 0.04292461 1

其实倾向性评分就是给通过结局变量和混杂因素formula = treat ~ age + educ + black + hispan + married + nodegree + re74 + re75做了一个回归分析，得到了distance就是该数据的倾向指数，然后依据倾向指数来均衡组间非试验因素的分布。

五、讨论与疑问

我们知道从总体抽取样本时候会出现抽样误差，那么从所有样本里面抽取匹配后的子样本，也同样存在导致误差可能（本质上是一种筛选样本造成的选择性偏倚），如何评估抽取样本是否有问题？常用的做法是通过第二种方法比如敏感度分析再次验证下结果。

另外一个重要指标是stand mean different，计算公式为（实验组均值 - 对照组均值）/ 实验组标准差。各匹配变量的值最好都要＜0.1，我们的变量里面只有re75略大于0.1，匹配质量可以接受。

summary(m.out, standardize = T)$sum.matched

## Means Treated Means Control SD Control Std. Mean Diff.
## distance 0.5002592 0.4924892 0.2476051 0.035272557
## age 25.7064220 26.2201835 11.0272036 -0.071804344
## educ 10.2385321 10.2568807 2.6784589 -0.009125716
## blackFALSE 0.2660550 0.2844037 0.4532137 -0.050331164
## blackTRUE 0.7339450 0.7155963 0.4532137 0.050331164
## hispanTRUE 0.1009174 0.1009174 0.3026107 0.000000000
## marriedTRUE 0.2385321 0.2293578 0.4223617 0.023360849
## nodegreeTRUE 0.6880734 0.6697248 0.4724845 0.040249987
## re74 2732.5276193 2280.9101034 4410.6701851 0.092419194
## re75 2042.2571592 1667.3486967 2977.9886750 0.116458293
## eCDF Med eCDF Mean eCDF Max
## distance 0.018348624 0.022913183 0.082568807
## age 0.073394495 0.083355177 0.229357798
## educ 0.027522936 0.028134557 0.073394495
## blackFALSE 0.009174312 0.009174312 0.018348624
## blackTRUE 0.009174312 0.009174312 0.018348624
## hispanTRUE 0.000000000 0.000000000 0.000000000
## marriedTRUE 0.004587156 0.004587156 0.009174312
## nodegreeTRUE 0.009174312 0.009174312 0.018348624
## re74 0.036697248 0.057363592 0.220183486
## re75 0.027522936 0.031647168 0.091743119

参考文献

[1]. 周支瑞, 胡志德. 疯狂统计学. 长沙:中南大学出版社, 2018.

[2]. 周支瑞, 胡志德. 聪明统计学. 长沙:中南大学出版社, 2016.