R语言倾向性评分
什么是倾向性评分
倾向性评分（Propensity Score）是统计学和因果推断中的一个概念，主要用于减少混杂偏差，从而更准确地估计某种处理或干预的效果。它指的是在给定一组协变量的条件下，个体接受某种处理（例如药物治疗、政策干预等）的概率。
举例说明：
假设你想评估某种新药对患者康复效果的影响。你有一组患者数据，其中包括患者是否服用了新药（处理变量），以及患者的年龄、性别、病情严重程度等信息（协变量）。
直接比较服药组和未服药组的康复率可能不公平，因为服药组可能包含更多病情严重的患者，或者年龄更大的患者，这些因素都会影响康复结果。你可以用患者的年龄、性别、病情等协变量，建立一个模型（如逻辑回归）来预测每个患者服药的概率，这个概率就是倾向性评分。
即是说，在非随机试验（观察性研究）中，处理组和对照组可能存在系统性差异（混杂因素）。通过计算每个个体接受处理的概率（倾向性评分），可以在分析时调整这些混杂因素，减少偏差，使处理组和对照组在这些背景变量上更加可比。换句话说，倾向性评分可理解为一个介于0和1之间的数值，反映了“在相似背景下，该个体接受处理的概率(倾向)”。
利用MatchIt 包进行分析
R语言中，MatchIt是进行倾性评分（Propensity Score Matching，PSM）的常用包。通过匹配处理组和对照组的样本，使得两组在协变量上的分布更加相似，从而减少混杂偏差，提高因果推断的准确性。
数据集
library("MatchIt")
## Warning: 程辑包'MatchIt'是用R版本4.3.3 来建造的
data("lalonde") # 使用 Lalonde 数据集
dplyr::glimpse(lalonde)
## Rows: 614
## Columns: 9
## $ treat <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ age <int> 37, 22, 30, 27, 33, 22, 23, 32, 22, 33, 19, 21, 18, 27, 17, 1…
## $ educ <int> 11, 9, 12, 11, 8, 9, 12, 11, 16, 12, 9, 13, 8, 10, 7, 10, 13,…
## $ race <fct> black, hispan, black, black, black, black, black, black, blac…
## $ married <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
## $ nodegree <int> 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1…
## $ re74 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ re75 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ re78 <dbl> 9930.0460, 3595.8940, 24909.4500, 7506.1460, 289.7899, 4056.4…
这是一份关于国家支持就业示范项目（NSW）和人口收入动态调查（PSID）数据的子样本。
数据包含614条观测记录，其中185人为实验组（接受干预），429人为对照组（未接受干预）。每条记录包含9个变量：
treat：是否接受干预（1=接受，0=未接受） age：年龄（岁） educ：受教育年限（年） race：种族/族裔（黑人、拉美裔、白人） married：婚姻状况（1=已婚，0=未婚） nodegree：是否无高中文凭（1=无文凭，0=有文凭） re74：1974年收入（美元） re75：1975年收入（美元） re78：1978年收入（美元） 其中，treat是处理变量，re78是最终结果变量，其余为处理前的协变量。
匹配前初始不平衡
# 创建预匹配对象
m.out0 <- matchit(treat ~ age + educ + race + married +
nodegree + re74 + re75,
data = lalonde,
method = NULL,
distance = "glm")
# 查看不平衡情况
summary(m.out0)
##
## Call:
## matchit(formula = treat ~ age + educ + race + married + nodegree +
## re74 + re75, data = lalonde, method = NULL, distance = "glm")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.1822 1.7941 0.9211 0.3774
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## raceblack 0.8432 0.2028 1.7615 . 0.6404
## racehispan 0.0595 0.1422 -0.3498 . 0.0827
## racewhite 0.0973 0.6550 -1.8819 . 0.5577
## married 0.1892 0.5128 -0.8263 . 0.3236
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## distance 0.6444
## age 0.1577
## educ 0.1114
## raceblack 0.6404
## racehispan 0.0827
## racewhite 0.5577
## married 0.3236
## nodegree 0.1114
## re74 0.4470
## re75 0.2876
##
## Sample Sizes:
## Control Treated
## All 429 185
## Matched 429 185
## Unmatched 0 0
## Discarded 0 0
通过summary(m.out0)得到的结果显示，处理组（treated）和对照组（control）在多项协变量上存在明显差异，比如年龄（age）、种族（race）、婚姻状况（married）和过去收入（re74、re75）等。 这些差异用标准化均值差（Std. Mean Diff.）等指标衡量，数值未接近0表明不平衡，此时未经调整，直接比较处理组和对照组的结果可能存在偏差。
执行匹配
# 最近邻倾向匹配
m.out1 <- matchit(treat ~ age + educ + race + married +
nodegree + re74 + re75,
data = lalonde,
method = "nearest",
distance = "glm")
# 完全匹配，使用 probit 链接函数
m.out2 <- matchit(treat ~ age + educ + race + married +
nodegree + re74 + re75,
data = lalonde,
method = "full",
distance = "glm",
link = "probit")
评估匹配质量
# 查看匹配后平衡情况
summary(m.out1, un = FALSE) # 只显示匹配后结果
##
## Call:
## matchit(formula = treat ~ age + educ + race + married + nodegree +
## re74 + re75, data = lalonde, method = "nearest", distance = "glm")
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.3629 0.9739 0.7566 0.1321
## age 25.8162 25.3027 0.0718 0.4568 0.0847
## educ 10.3459 10.6054 -0.1290 0.5721 0.0239
## raceblack 0.8432 0.4703 1.0259 . 0.3730
## racehispan 0.0595 0.2162 -0.6629 . 0.1568
## racewhite 0.0973 0.3135 -0.7296 . 0.2162
## married 0.1892 0.2108 -0.0552 . 0.0216
## nodegree 0.7081 0.6378 0.1546 . 0.0703
## re74 2095.5737 2342.1076 -0.0505 1.3289 0.0469
## re75 1532.0553 1614.7451 -0.0257 1.4956 0.0452
## eCDF Max Std. Pair Dist.
## distance 0.4216 0.9740
## age 0.2541 1.3938
## educ 0.0757 1.2474
## raceblack 0.3730 1.0259
## racehispan 0.1568 1.0743
## racewhite 0.2162 0.8390
## married 0.0216 0.8281
## nodegree 0.0703 1.0106
## re74 0.2757 0.7965
## re75 0.2054 0.7381
##
## Sample Sizes:
## Control Treated
## All 429 185
## Matched 185 185
## Unmatched 244 0
## Discarded 0 0
# 可视化倾向得分分布
plot(m.out1, type = "jitter", interactive = FALSE)
# 可视化协变量密度图
plot(m.out1, type = "density", interactive = FALSE,
which.xs = ~age + married + re75)
plot(summary(m.out1))
- 最近邻匹配（m.out1）:
匹配后，协变量的差异有所减小，但仍存在不平衡的情况（部分变量Std. Mean Diff.仍较大）。 匹配后样本量减少，部分对照组样本未被匹配（244个未匹配），说明匹配剔除了部分样本以提高组间可比性。 散点图和密度图显示处理组和对照组的倾向得分和协变量分布更接近，但个别变量如年龄的平衡效果一般。 本例中最近邻匹配改善了平衡但效果有限，可能导致估计偏差和不稳定。
# 查看匹配后平衡情况
summary(m.out2, un = FALSE)
##
## Call:
## matchit(formula = treat ~ age + educ + race + married + nodegree +
## re74 + re75, data = lalonde, method = "full", distance = "glm",
## link = "probit")
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5773 0.5764 0.0045 0.9949 0.0043
## age 25.8162 25.5347 0.0393 0.4790 0.0787
## educ 10.3459 10.5381 -0.0956 0.6192 0.0253
## raceblack 0.8432 0.8389 0.0119 . 0.0043
## racehispan 0.0595 0.0492 0.0435 . 0.0103
## racewhite 0.0973 0.1119 -0.0493 . 0.0146
## married 0.1892 0.1633 0.0660 . 0.0259
## nodegree 0.7081 0.6577 0.1110 . 0.0504
## re74 2095.5737 2100.2150 -0.0009 1.3467 0.0314
## re75 1532.0553 1561.4420 -0.0091 1.5906 0.0536
## eCDF Max Std. Pair Dist.
## distance 0.0486 0.0198
## age 0.2742 1.2843
## educ 0.0730 1.2179
## raceblack 0.0043 0.0162
## racehispan 0.0103 0.4412
## racewhite 0.0146 0.3454
## married 0.0259 0.4473
## nodegree 0.0504 0.9872
## re74 0.1881 0.8387
## re75 0.1984 0.8240
##
## Sample Sizes:
## Control Treated
## All 429. 185
## Matched (ESS) 50.76 185
## Matched 429. 185
## Unmatched 0. 0
## Discarded 0. 0
# 创建 Love plot
plot(summary(m.out2))
- 全匹配（m.out2）:
全匹配后，协变量平衡显著改善，标准化均值差均接近0（通常小于0.1被视为良好平衡）。 匹配后所有样本均被保留，无样本丢失。 Love plot（平衡图）显示匹配后各协变量的平衡指标均大幅改善。 全匹配结合probit倾向得分模型更适合本数据，实现了良好的平衡。
平均处理效应（ATT）
# 提取匹配后的数据
m.data <- match_data(m.out2)
#dplyr::glimpse(m.data)
# 拟合结果模型
fit <- lm(re78 ~ treat * (age + educ + race + married + nodegree + re74 + re75),
data = m.data,
weights = weights)
利用全匹配后的数据拟合加权线性回归模型，估计处理组相对于对照组的差异（本例中为1978年收入re78）。
# 效应估计
library("marginaleffects")
avg_comparisons(fit,
variables = "treat",
vcov = ~subclass,
newdata = subset(m.data, treat == 1))
##
## Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1977 704 2.81 0.00501 7.6 596 3357
##
## Term: treat
## Type: response
## Comparison: 1 - 0
估计平均处理效应（ATT）是1977美元，表示接受处理（NSW项目）的人群，平均1978年的收入比未接受处理的人高1977美元。
Std=704，z=2.81，对应p=0.005，统计显著，说明该效应不太可能由随机波动引起。置信区间（596, 3357）。
注意事项
- 匹配不能替代良好的研究设计，协变量选择至关重要
- 可尝试多种匹配方法并比较结果，例如倾向性评分加权（IPTW），分层等
- 实际分析中匹配后仍需检查平衡性，不平衡时需要调整匹配方法
更多信息请可参考 MatchIt包示例：
vignette("matching-methods")`
vignette("assessing-balance")`
vignette("estimating-effects")