倾向得分匹配(Propensity Score Matching, PSM)
在运用R进行生信分析时,我们经常会遇到需要我们确定某种蛋白或药物是否有效的课题,此时往往都需要用两组数据或更多去看其生效情况,但如果我们在数据清洗后直接查看结果,经常会得出其很有效的错觉。
其实这是因为很多时候,单从数据的定义本身,不一定能确定其效果以及要研究的事物真的作用到了其身上。有可能只是因为其他的潜在变量不一样,但由于其隐藏在看似不重要的数据后,不能直观地看出来。此时,PSM就显得尤为重要了,它的作用本质是起到控制变量的效果,也就是找出了与实验组相对的,我们耳熟能详的对照组。唯有尽量地减少其他因素对结果的干扰,我们才能更严谨地根据结果下判断。当然了,PSM不是真的给出两组有着接近相同的特征的数据,而是计算(倾向得分)概率值,其结果代表了这个数据基于它的特征会受到外界干预的可能性,其作用就是找出两组拥有相近倾向得分的数据。
我们用一个例子来说明:
library(dplyr)
set.seed(123)
# 生成数据:年龄、教育年限、性别是否影响参加培训
n <- 1000
data <- tibble(
age = rnorm(n, mean = 30, sd = 5), # 年龄
education = rnorm(n, mean = 12, sd = 2), # 教育年限
male = rbinom(n, 1, 0.5), # 性别(1=男,0=女)
# 倾向得分:年龄和教育越高,参加培训的概率越大
propensity_score = plogis(0.1 * age + 0.3 * education - 0.5 * male - 5),
treated = rbinom(n, 1, propensity_score), # 是否参加培训(干预)
# 工资:培训+2000元,但年龄和教育也影响工资
wage = 20000 + 500 * age + 1000 * education + 2000 * treated + rnorm(n, sd = 1000)
)
install.packages("Matching")
library(Matching)
# 使用 Matching 包进行倾向得分匹配
psm_model <- Match(
Y = data$wage, # 结果变量(工资)
Tr = data$treated, # 干预变量(是否培训)
X = data$propensity_score # 倾向得分
)
# 查看匹配效果
summary(psm_model) # 检查匹配后协变量平衡性
matched_indices <- c(psm_model$index.treated, psm_model$index.control)
matched_data <- data[matched_indices, ]
# 匹配前的比较(有偏差)
cat("匹配前效应估计:",
mean(data$wage[data$treated == 1]) - mean(data$wage[data$treated == 0]), "\n")
# 匹配后的比较(更准确)
cat("匹配后效应估计:",
mean(matched_data$wage[matched_data$treated == 1]) -
mean(matched_data$wage[matched_data$treated == 0]))
输出:
Estimate... 2226.8
AI SE...... 175.83
T-stat..... 12.665
p.val...... < 2.22e-16
Original number of observations.............. 1000
Original number of treated obs............... 761
Matched number of observations............... 761
Matched number of observations (unweighted). 878
匹配前效应估计: 4207.371
匹配后效应估计: 2216.294
首先,从直接比较和进行PSM比较后的效果比对可以发现,两组数据的工资差异不能直接用来说明是培训造成的,因为还有年龄和教育的差异被计算进去。而选择了相近得分后,可以发现差异只有两百多了;而从T-stat和p-val来看,前者的绝对值远大于2,且后者绝对值很接近于0,说明干预的效果很好,即培训的工资的提升很显著。当然,在日常对数据的分析中,不能单看这两个值,因为在样本很大的情况下,即便其本来的效果很小,也会体现出显著性。