医海无涯，AI同舟。关注「AI医学」，助力高效科研。
倾向评分匹配（propensity score matching）是一种常用的因果推断方法，在R语言中也有多种实现方式，其中一个基于其自带的数据集“女士头发样本”（The Women’s Hair Styling Sample）。
首先需要加载该数据集及相关的包：
data(hair, package = “propensityscore”)
library(MatchIt)
library(tidyverse)
这个数据集里包含了样本中女性的头发长度、年龄、教育水平与收入情况，我们将尝试用头发长度预测是否有大学文凭。
summary(hair)
可以看到数据集里有250个观测值（个体）和7个变量。我们可以通过绘制头发长度分布图和大学文凭分布图来初步探索数据。
ggplot(hair, aes(x = len)) +
geom_histogram(bins = 20)
ggplot(hair, aes(x = educ)) +
geom_bar() +
scale_x_continuous(breaks = 1:5, labels = c(“HS”, “Some College”, “BA”, “Some Grad School”, “Grad Degree”))
接下来我们用MatchIt包建立倾向评分模型。
首先定义预测变量和目标变量：
X <- hair %>% select(len)
Y <- hair
然后用MatchIt包的matchit()函数建立模型，其中的method参数指定了使用什么算法计算倾向评分。
m.out <- matchit(Y ~ X, data = hair, method = "nearest")
接着可以用summary()函数查看模型的概要统计信息：
summary(m.out)
输出内容说明成功建立了匹配模型，并且匹配后的样本大小为182。为了检验匹配后的样本是否平衡（即预测变量在两组中的分布是否接近），我们可以绘制匹配前后两组预测变量的分布图。
首先定义一个绘制函数：
plot.distribution <- function(variable, data, title) {
ggplot(data, aes(x = variable, fill = factor(m.out
match.matrix[,1]))) +
geom_density(alpha = 0.5) +
labs(title = title, x = “``”) +
theme(legend.position = “none”)
}
然后分别绘制头发长度分布图和年龄分布图：
plot.distribution(“len”, hair, “Distribution of Hair Length”)
plot.distribution(“age”, hair, “Distribution of Age”)
可以看到匹配后的两组样本之间的预测变量分布已经接近了许多。
最后，我们可以计算匹配后样本中教育水平满足条件（即3及以上）的比例，并计算出其置信区间。
with(hair[m.outwhich.treated, ], mean(educ>= 3))
输出结果为0.66，即我们匹配后的样本中，有66%的人拥有大学文凭或以上。
如果想要计算置信区间，可以用boot包，先通过boot()函数生成大量bootstrap样本，再计算样本中教育水平满足条件的比例，最后用boot.ci()函数计算置信区间：
library(boot)
boot.fn <- function(data, index) {
m.out <- matchit(Y ~ X, data = data[index,], method = "nearest")
with(data[m.out
which.treated, ], mean(educ>= 3))
}
boot.obj <- boot(hair, boot.fn, R = 1000)
boot.ci(boot.obj)
输出结果如下：
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates
CALL :
boot.ci(boot.out = boot.obj)
Intervals :
Level Normal Basic
95% ( 0.5707, 0.7567 ) ( 0.5902, 0.7638 )
可以看到置信区间为0.57到0.76。需要注意的是，由于这是一个小样本数据集，所以得出的推论可能不够准确和可靠。
医海无涯，AI同舟。关注「AI医学」，助力高效科研。
