使用R进行倾向得分匹配（PSM）

【译文】使用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)

P.S.2：非常感谢我的同事Katharina Kuba向我推荐了MatchIt包！

注：本文原载与datascienceplus网站

原文链接： http://datascienceplus.com/how-to-use-r-for-matching-samples-propensity-score