【译文】使用R进行倾向得分匹配（PSM）
作者 Norbert Köhler
译者 钱亦欣
根据维基百科，倾向得分匹配（PSM）是一种用来评估处置效应的统计方法。广义说来，它将样本根据其特性分类，而不同类样本间的差异就可以看作处置效应的无偏估计。因此，PSM不仅仅是随机试验的一种替代方法，它也是流行病研究中进行样本比较的重要方法之一。让我们举个栗子：
与健康相关的生活质量（HRQOL）被认为是癌症治疗的重要结果之一。对癌症患者而言，最常用的HRQOL测度是通过欧洲癌症研究与治疗中心的调查问卷计算得出的。EORTC QLD-C30是一个由30个项目组成，包括5个功能量表，9个症状量表和一个全球生活质量量表的的问卷。所有量表都会给出一个0-100之间的得分。症状量表得分越高代表被调查人生活压力越大，其余两个量表得分越高代表生活质量越高。
然而，如果没有任何参照，直接对数据进行解释是很困难的。幸运的是，EORTC QLQ-C30问卷也在一些一般人群调查中使用，我们可以对比患者的得分和一般人群的得分差异，从而判断患者的负担症状和一些功能障碍是否能归因于癌症治疗。PSM在这里可以以年龄和性别等特征，将相似的患者和一般人群进行匹配。
本文我会演示如在在R中实现PSM。更为详尽的说明请参考： “A Step-by-Step Guide to Propensity Score Matching in R” 。
生成两个随机数据框
由于我不希望在本文使用真实数据，我需要生成一些仿真数据。使用Wakefield包可以很容易地实现这个功能。
第一步，我们创建一个名为df.patients的数据框，我希望它包含250个病人的年龄和性别数据，所有病人的年龄都要在30-78岁之间，并且70%的病人被设定为男性。
set.seed(1234)
df.patients <- r_data_frame(n = 250, age(x = 30:78, name = 'Age'), sex(x = c("Male", "Female"), prob = c(0.70, 0.30), name = "Sex"))
df.patients$Sample <- as.factor('Patients')summary函数会返回创建的数据框的基本信息，如你所见，患者平均年龄为53.7岁，并且大约70%为男性。
summary(df.patients)
## Age Sex Sample
## Min. :30.00 Male :173 Patients:250
## 1st Qu.:42.00 Female: 77
## Median :54.00
## Mean :53.71
## 3rd Qu.:66.00
## Max. :78.00第二步，我们需要创建另一个名为df.population的数据框。我希望这个数据集的数据和患者的有些不同，因此正常人群的年龄区间被设定为18-80岁，并且男女各占一半。
set.seed(1234)
df.population <- r_data_frame(n = 1000, age(x = 18:80, name = 'Age'), sex(x = c("Male", "Female"), prob = c(0.50, 0.50), name = "Sex"))
df.population$Sample <- as.factor('Population')下方表格显示样本平均年龄为49.5岁，男女比例也大致相等。
summary(df.population)
## Age Sex Sample
## Min. :18.00 Male :485 Population:1000
## 1st Qu.:34.00 Female:515
## Median :50.00
## Mean :49.46
## 3rd Qu.:65.00
## Max. :80.00合并数据框
在匹配样本之前，我们需要把两个数据框合并。先生成一个新变量Group来代表观测来自哪个全体（逻辑型变量），再添加另一个变量Distress来反应个体的痛苦程度。Distress变量是利用Wakefield包中的age函数创建的，可以发现，女性承受的痛苦级别更高。
mydata <- rbind(df.patients, df.population)
mydata$Group <- as.logical(mydata$Sample == 'Patients')
mydata$Distress <- ifelse(mydata$Sex == 'Male', age(nrow(mydata), x = 0:42, name = 'Distress'),age(nrow(mydata), x = 15:42, name = 'Distress'))当我们比较两类样本的年龄和性别分布时，我们可以发现明显的区别：
pacman::p_load(tableone)
table1 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), data = mydata, factorVars = 'Sex', strata = 'Sample')
table1 <- print(table1, printToggle = FALSE, noSpaces = TRUE)
kable(table1[,1:3], align = 'c', caption = 'Table 1: Comparison of unmatched samples')
更进一步，我们还发现一般人群的痛苦程度显著较高。
样本匹配
现在，我们已经完成了全部的准备工作，可以开始使用MatchIT包中的matchit函数来匹配两类样本了。函数中method=‘nearest’的设定指明了使用近邻法进行匹配。其他方法包括，次分类，优化匹配等。ratio=1意味着这是一一配对。同时也请注意Group变量需要是逻辑型变量。
set.seed(1234)
match.it <- matchit(Group ~ Age + Sex, data = mydata, method="nearest", ratio=1)
a <- summary(match.it)为了后续工作的便利，我们将summary函数的输出赋值给名为a的变量。
在匹配万样本后，一般人群样本量所见到了和患者样本一致（250个观测）。
kable(a$nn, digits = 2, align = 'c', caption = 'Table 2: Sample sizes')
根据输出结果，匹配后的年龄和性别分布基本一致了。
kable(a$sum.matched[c(1,2,4)], digits = 2, align = 'c', caption = 'Table 3: Summary of balance for matched data')
倾向得分的分布可以使用MatchIt包中的plot函数进行绘制。
plot(match.it, type = 'jitter', interactive = FALSE)输出如下：
保存匹配样本
最后，让我们把匹配好的样本保存在df.match数据框里。
df.match <- match.data(match.it)[1:ncol(mydata)]
rm(df.patients, df.population)现在pacman::p_load(tableone)
table4 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), data = df.match, factorVars = 'Sex', strata = 'Sample')
table4 <- print(table4, printToggle = FALSE, noSpaces = TRUE)
kable(table4[,1:3], align = 'c', caption = 'Table 4: Comparison of matched samples')，我们可以对比两类人群间痛苦程度的差异是否依旧显著。
由于p值为0.222，学生t检验的结果不再显著。因此，PSM帮助我们避免犯下第一类错误。
P.S.1：本文只用的所有包可通过如下代码加载：
pacman::p_load(knitr, wakefield, MatchIt, tableone, captioner)
