高山流水

但问耕耘，莫问收获

0%

倾向性评分匹配的R语言实现

发表于 更新于 分类于

说明
1. 译自Norbert Köhler的文章How to use R for matching samples (propensity score)

根据维基百科，倾向性评分匹配（propensity score matching,PSM）是一种用来评估处置效应的统计方法。

Propensity score matching (PSM) is a statistical matching technique that attempts to estimate the effect of a treatment, policy, or other intervention by accounting for the covariates that predict receiving the treatment

广义说来，倾向性评分分析假设只有当两组样本的研究对象具有相似特征时才能进行两组样本之间的无偏比较(unbiased comparison)。因此，PSM不仅仅是随机试验的一种替代方法，它也是流行病研究中进行样本比较的重要方法之一。让我们举个例子：

与健康相关的生活质量（Health-related quality of life,HRQOL）被认为是癌症治疗的重要结果之一。对癌症患者而言，最常用的HRQOL测度是通过欧洲癌症研究与治疗中心(European Organisation for Research and Treatment of Cancer,EORTC)的调查问卷计算得出的。

EORTC QLD-C30是一个由30个项目组成，包括5个功能量表，9个症状量表和一个整体生活质量量表的的问卷。所有量表都会给出一个0-100之间的得分。症状量表得分越高代表被调查人症状越重，其余两个量表得分越高代表生活质量越高。

然而，如果没有任何参照，直接对数据进行解释是很困难的。幸运的是，EORTC QLQ-C30问卷也在一些正常人群调查中使用，我们可以对比患者的得分和正常人群的得分差异，从而判断患者的一些症状和功能障碍是否能归因于癌症（或治疗）。PSM可以通过相似年龄和相同性别等特征，将患者和一般人群进行匹配。

生成两个随机数据框

本文用到的所有包可通过如下代码加载：

1
pacman::p_load(knitr, wakefield, MatchIt, tableone, captioner)

由于我们不希望在博客文章中使用真实数据，因此需要生成一些模拟数据。使用wakefield包可以很容易地实现这个功能。

  1. 创建一个名为df.patients的数据框，包含250个病人的年龄和性别数据，所有病人的年龄都要在30-78岁之间，并且70%的病人被设定为男性。

    1
    2
    3
    4
    5
    6
    7
    set.seed(180717)
    library(wakefield)
    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')

  2. summary函数会返回创建的数据框的基本信息。患者平均年龄为53.7岁，并且大约70%为男性(66.8%)。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    summary(df.patients)
         
          Age            Sex           Sample   
    Min.   :30.00   Male  :167   Patients:250  
    1st Qu.:44.00   Female: 83                 
    Median :57.00                              
    Mean   :55.24                              
    3rd Qu.:67.00                              
    Max.   :78.00

  3. 接下来，需要创建另一个名为df.population的数据框。我们希望这个数据集的数据和患者df.patients有些不同，因此正常人群的年龄区间被设定为18-80岁，并且男女各占一半。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    set.seed(180717)

    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')
     
    #下方表格显示样本平均年龄为50.1岁，男女比例也大致相等。

    summary(df.population)
         
    	Age           Sex             Sample    
    Min.   :18.0   Male  :489   Population:1000  
    1st Qu.:34.0   Female:511                    
    Median :51.0                                 
    Mean   :50.1                                 
    3rd Qu.:66.0                                 
    Max.   :80.0

合并数据框

在匹配样本之前，我们需要把两个数据框合并。首先，生成一个新变量Group来代表研究对象来自哪个全体（逻辑型变量），再添加另一个变量Distress来表示个体的痛苦程度。Distress变量是利用wakefield包中的age函数创建的，可以发现，女性承受的痛苦级别更高。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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'))


> mydata
# A tibble: 1,250 x 5
     Age Sex    Sample   Group Distress
   <int> <fct>  <fct>    <lgl>    <int>
 1    76 Female Patients TRUE        39
 2    49 Female Patients TRUE        25
 3    57 Male   Patients TRUE        39
 4    62 Male   Patients TRUE        32
 5    61 Male   Patients TRUE        21
 6    44 Male   Patients TRUE        27
 7    56 Female Patients TRUE        21
 8    54 Female Patients TRUE        28
 9    49 Male   Patients TRUE        17
10    64 Male   Patients TRUE        28
# ... with 1,240 more rows

if..else语句的基本语法是

1
2
3
4
5
if(boolean_expression) {
   // statement(s) will execute if the boolean expression is true.
} else {
   // statement(s) will execute if the boolean expression is false.
}

当我们比较两组样本的年龄和性别分布时，我们可以发现明显的区别：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(pacman)
pacman::p_load(tableone)
table1 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = mydata, 
                         factorVars = 'Sex', 
                         strata = 'Sample')
table1 <- print(table1, 
                printToggle = FALSE, 
                noSpaces = TRUE)

library(knitr)
kable(table1[,1:3],  
      align = 'c', 
      caption = 'Table 1: Comparison of unmatched samples')

Table 1: Comparison of unmatched samples

Patients Population p
n 250 1000
Age (mean (sd)) 55.24 (14.24) 50.10 (18.37) <0.001
Sex = Female (%) 83 (33.2) 511 (51.1) <0.001
Distress (mean (sd)) 25.12 (11.58) 24.35 (11.34) 0.344

样本匹配

现在，我们已经完成了全部的准备工作，可以开始使用MatchIT包中的matchit函数来匹配两组样本了。

1
2
3
4
5
6
7
8
9
library(MatchIt)
set.seed(180717)

match.it <- matchit(Group ~ Age + Sex, 
data = mydata, 
method="nearest", 
ratio=1)

a <- summary(match.it)
  1. 函数中method=‘nearest’的设定指明了使用近邻法(the nearest neighbors method)进行匹配。
  2. 其他方法包括exact matching,subclassification,optimal matching,genetic matching,ull matching等（method = c("exact", "subclass", "optimal", ""genetic", "full")）。
  3. ratio=1意味着这是一一配对。在我们的例子里，对于患者组中的每一个病例，将恰好匹配正常人群样本中的一个人。
  4. 同时也请注意Group变量需要是逻辑型变量（TRUE与FALSE）。

为了后续工作的便利，我们将summary函数的输出赋值给名为a的变量。在匹配完样本后，正常人群样本量减少到了和患者样本一致（n=250，详见Table 2）。

1
kable(a$nn, digits = 2, align = 'c', caption = 'Table 2: Sample sizes')

Table 2: Sample sizes

Control Treated
All 1000 250
Matched 250 250
Unmatched 750 0
Discarded 0 0

根据输出结果，匹配后的年龄和性别分布基本一致了。

1
kable(a$sum.matched[c(1,2,4)], digits = 2, align = 'c', caption = 'Table 3: Summary of balance for matched data')

Table 3: Summary of balance for matched data

Means Treated Means Control Mean Diff
distance 0.23 0.23 0.00
Age 55.24 55.22 0.02
SexMale 0.67 0.67 0.00
SexFemale 0.33 0.33 0.00

倾向性评分的分布可以使用MatchIt包中的plot函数进行绘制。

1
plot(match.it, type = 'jitter', interactive = FALSE)

保存匹配样本

最后，匹配的样本将保存到名为df.match的新数据框中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# 获得匹配数据集(含distance和weights列)的第1~ncol(mydata)列

df.match <- match.data(match.it)[1:ncol(mydata)]

# rm()函数：删除变量函数
rm(df.patients, df.population)

# 我们可以对比两类人群间痛苦程度的差异是否依旧显著
pacman::p_load(tableone)
table4 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = df.match, 
                         factorVars = 'Sex', 
                         strata = 'Sample')

# printToggle：是否打印输出。如果为FALSE，则不创建输出，并返回一个矩阵。
# noSpaces：是否删除为对齐添加的空格。如果您希望自己在其他软件中对齐数字，请使用此选项。

table4 <- print(table4, 
                printToggle = FALSE, 
                noSpaces = TRUE)

# align/列对齐: 'l' (left), 'c' (center) and/or 'r' (right). 默认则align = NULL, 数值列右对齐, 其他列左对齐.
# kable函数：Create Tables In LaTeX, HTML, Markdown And ReStructuredText

kable(table4[,1:3],  
      align = 'c', 
      caption = 'Table 4: Comparison of matched samples')
Patients Population p
n 250 250
Age (mean (sd)) 55.24 (14.24) 55.22 (14.21) 0.987
Sex = Female (%) 83 (33.2) 83 (33.2) 1.000
Distress (mean (sd)) 25.12 (11.58) 22.87 (11.72) 0.032