使用R進(jìn)行傾向得分匹配（PSM）

根據(jù)維基百科，傾向得分匹配（PSM）是一種用來(lái)評(píng)估處置效應(yīng)的統(tǒng)計(jì)方法。廣義說(shuō)來(lái)，它將樣本根據(jù)其特性分類，而不同類樣本間的差異就可以看作處置效應(yīng)的無(wú)偏估計(jì)。因此，PSM不僅僅是隨機(jī)試驗(yàn)的一種替代方法，它也是流行病研究中進(jìn)行樣本比較的重要方法之一。讓我們舉個(gè)栗子：

與健康相關(guān)的生活質(zhì)量（HRQOL）被認(rèn)為是癌癥治療的重要結(jié)果之一。對(duì)癌癥患者而言，最常用的HRQOL測(cè)度是通過(guò)歐洲癌癥研究與治療中心的調(diào)查問(wèn)卷計(jì)算得出的。EORTC QLD-C30是一個(gè)由30個(gè)項(xiàng)目組成，包括5個(gè)功能量表，9個(gè)癥狀量表和一個(gè)全球生活質(zhì)量量表的的問(wèn)卷。所有量表都會(huì)給出一個(gè)0-100之間的得分。癥狀量表得分越高代表被調(diào)查人生活壓力越大，其余兩個(gè)量表得分越高代表生活質(zhì)量越高。

然而，如果沒(méi)有任何參照，直接對(duì)數(shù)據(jù)進(jìn)行解釋是很困難的。幸運(yùn)的是，EORTC QLQ-C30問(wèn)卷也在一些一般人群調(diào)查中使用，我們可以對(duì)比患者的得分和一般人群的得分差異，從而判斷患者的負(fù)擔(dān)癥狀和一些功能障礙是否能歸因于癌癥治療。PSM在這里可以以年齡和性別等特征，將相似的患者和一般人群進(jìn)行匹配。

生成兩個(gè)隨機(jī)數(shù)據(jù)框

由于我不希望在本文使用真實(shí)數(shù)據(jù)，我需要生成一些仿真數(shù)據(jù)。使用Wakefield包可以很容易地實(shí)現(xiàn)這個(gè)功能。

第一步，我們創(chuàng)建一個(gè)名為df.patients的數(shù)據(jù)框，我希望它包含250個(gè)病人的年齡和性別數(shù)據(jù)，所有病人的年齡都要在30-78歲之間，并且70%的病人被設(shè)定為男性。

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函數(shù)會(huì)返回創(chuàng)建的數(shù)據(jù)框的基本信息，如你所見(jiàn)，患者平均年齡為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

第二步，我們需要?jiǎng)?chuàng)建另一個(gè)名為df.population的數(shù)據(jù)框。我希望這個(gè)數(shù)據(jù)集的數(shù)據(jù)和患者的有些不同，因此正常人群的年齡區(qū)間被設(shè)定為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

合并數(shù)據(jù)框

在匹配樣本之前，我們需要把兩個(gè)數(shù)據(jù)框合并。先生成一個(gè)新變量Group來(lái)代表觀測(cè)來(lái)自哪個(gè)全體（邏輯型變量），再添加另一個(gè)變量Distress來(lái)反應(yīng)個(gè)體的痛苦程度。Distress變量是利用Wakefield包中的age函數(shù)創(chuàng)建的，可以發(fā)現(xiàn)，女性承受的痛苦級(jí)別更高。

mydata <- rbind(df.patients, df.population) mydata$Group <- as.logical(mydata$Sample == 'Patients') mydata$Distress <- ifelsmydata <- 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'))

當(dāng)我們比較兩類樣本的年齡和性別分布時(shí)，我們可以發(fā)現(xiàn)明顯的區(qū)別：

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')

更進(jìn)一步，我們還發(fā)現(xiàn)一般人群的痛苦程度顯著較高。

樣本匹配

現(xiàn)在，我們已經(jīng)完成了全部的準(zhǔn)備工作，可以開(kāi)始使用MatchIT包中的matchit函數(shù)來(lái)匹配兩類樣本了。函數(shù)中method=‘nearest’的設(shè)定指明了使用近鄰法進(jìn)行匹配。其他方法包括，次分類，優(yōu)化匹配等。ratio=1意味著這是一一配對(duì)。同時(shí)也請(qǐng)注意Group變量需要是邏輯型變量。

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

為了后續(xù)工作的便利，我們將summary函數(shù)的輸出賦值給名為a的變量。

在匹配萬(wàn)樣本后，一般人群樣本量所見(jiàn)到了和患者樣本一致（250個(gè)觀測(cè)）。

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

根據(jù)輸出結(jié)果，匹配后的年齡和性別分布基本一致了。

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

傾向得分的分布可以使用MatchIt包中的plot函數(shù)進(jìn)行繪制。

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

輸出如下：

保存匹配樣本

最后，讓我們把匹配好的樣本保存在df.match數(shù)據(jù)框里。

df.match <- match.data(match.it)[1:ncol(mydata)]
rm(df.patients, df.population)

現(xiàn)在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')，我們可以對(duì)比兩類人群間痛苦程度的差異是否依舊顯著。

由于p值為0.222，學(xué)生t檢驗(yàn)的結(jié)果不再顯著。因此，PSM幫助我們避免犯下第一類錯(cuò)誤。

P.S.1：本文只用的所有包可通過(guò)如下代碼加載：數(shù)據(jù)分析師培訓(xùn)

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