陳亨貴 潘雄飛 黃媛媛 何保昌 陳 法 潘 安△

·計(jì)算機(jī)應(yīng)用·

傾向評(píng)分配比法及其在Stata軟件實(shí)現(xiàn)

陳亨貴1潘雄飛1黃媛媛2何保昌2陳 法2潘 安1△

傾向評(píng)分法是由Rosenbaum和Rubin在1983年首次提出[1]，主要用于觀察性研究組中混雜因素的事后均衡，對(duì)非研究混雜因素進(jìn)行類似隨機(jī)化的均衡處理。近年來，傾向評(píng)分法的應(yīng)用越來越廣泛，尤其在臨床[2]，許多學(xué)者利用其來降低混雜偏倚，達(dá)到協(xié)變量在處理組和對(duì)照組的分布均衡可比[3]，以獲得類似前瞻性隨機(jī)對(duì)照試驗(yàn)的效果[3-4]。目前已有學(xué)者撰文如何在SPSS、SAS中實(shí)現(xiàn)傾向評(píng)分配比法[5-6]，部分Stata書籍中也介紹過此方法，但最多僅講述到如何算出傾向評(píng)分，并未提及算出傾向評(píng)分后如何將配對(duì)的個(gè)體篩選出來，且目前國內(nèi)尚無關(guān)于在Stata實(shí)現(xiàn)傾向評(píng)分配比的報(bào)道。因此，本文旨在通過實(shí)例演示如何在Stata上實(shí)現(xiàn)傾向評(píng)分配比法，從而為更好地運(yùn)用傾向評(píng)分配比法提供參考。

原理和方法

傾向評(píng)分是指在給定條件下，研究對(duì)象進(jìn)入到處理組或者對(duì)照組的條件概率[7]。定義如下：

e(xi)=P(Zi=1|X=xi)

(1)

假定Zi和Xi相互獨(dú)立，則可利用logit或probit模型計(jì)算出每個(gè)對(duì)象的傾向評(píng)分[8]，公式如下：

(2)

其中i為研究對(duì)象(1～N)，Z為分組變量，X為特征變量，P為傾向評(píng)分(0～1)。

傾向評(píng)分配比法則是在此基礎(chǔ)上，將處理組與對(duì)照組中傾向評(píng)分接近或相同的研究對(duì)象進(jìn)行配對(duì)，從而使兩組間的協(xié)變量分布均衡。其配比比例可視處理組和對(duì)照組的樣本量而定，如果兩者相差較多可進(jìn)行1∶n配比，但建議n最高不得超過4。關(guān)于配比的算法主要有全局最優(yōu)匹配法(global optimal algorithms)和局部最優(yōu)匹配法(local optimal algorithms)兩種[7]，前者是以保證整個(gè)數(shù)據(jù)集傾向評(píng)分的總體差值最小為出發(fā)點(diǎn)，但卻不能保證處理組和對(duì)照組中的兩個(gè)研究對(duì)象的差值是最小的，因而在實(shí)際應(yīng)用中較少；后者是從對(duì)照組中選擇與處理組的某個(gè)研究對(duì)象傾向評(píng)分最為接近或相同的個(gè)體進(jìn)行配比，直至所有的處理組對(duì)象配對(duì)完成。局部最優(yōu)匹配主要通過設(shè)定卡鉗值(δ)來完成，其中卡鉗值是指處理組和對(duì)照組中兩個(gè)研究對(duì)象的傾向評(píng)分差值，表示形式如下：

δ>|Pi-Pj|=min{|Pi-Pj|}

(3)

其中Pi為處理組對(duì)象的傾向評(píng)分，Pj為對(duì)照組對(duì)象的傾向評(píng)分。

該值設(shè)的越大，處理組可以匹配到的對(duì)照組的研究對(duì)象越多，但同時(shí)也會(huì)加大估計(jì)處理效應(yīng)的偏倚；相反，如果卡鉗值設(shè)的越小，處理組匹配到的對(duì)象就會(huì)越少，雖然匹配后的均衡性會(huì)更好，但有可能損失較大的樣本量從而降低估計(jì)處理效應(yīng)的準(zhǔn)確性，Austin等通過大量的模擬實(shí)驗(yàn)研究結(jié)果顯示最優(yōu)卡鉗值是0.02或0.03[9]，而王永吉等利用蒙特卡洛模擬研究發(fā)現(xiàn)使用logit模型算出的傾向評(píng)分合并標(biāo)準(zhǔn)偏差的20%是最優(yōu)的卡鉗值[10]。此外，匹配過程中對(duì)照是否允許被放回，即被配對(duì)的對(duì)照組的研究對(duì)象是否準(zhǔn)予被重復(fù)利用目前尚存有爭議[8]，但一般情況下不允許放回，即不能重復(fù)利用已配對(duì)的對(duì)象。

實(shí)例分析

數(shù)據(jù)來源于方積乾主編的衛(wèi)生統(tǒng)計(jì)學(xué)(第7版)[11]，原數(shù)據(jù)為探討鼻咽癌發(fā)病的危險(xiǎn)因素，研究人員對(duì)某醫(yī)院腫瘤防治中心放療科的105例鼻咽癌新發(fā)病例與130名健康人進(jìn)行病例對(duì)照研究，結(jié)果發(fā)現(xiàn)吸煙是鼻咽癌的危險(xiǎn)因素之一。本文為更好地介紹如何在Stata中實(shí)現(xiàn)傾向評(píng)分配比法，用該數(shù)據(jù)作為例子，以有無吸煙作為 “處理因素”，探討吸煙對(duì)鼻咽癌發(fā)病的影響。具體變量命名及賦值見表1。

表1 各變量命名及賦值

*：原始數(shù)據(jù)中飲茶變量中有個(gè)觀測值為“3”，本例將其替換為中位數(shù)“0”。

操作步驟

首先在Stata軟件中安裝實(shí)現(xiàn)傾向評(píng)分配比法的程序包psmatch2，其次以傾向評(píng)分1:1配比法為例，運(yùn)用logit模型計(jì)算出每個(gè)研究對(duì)象的傾向評(píng)分，并利用卡鉗值匹配為每個(gè)處理組找到適合的配對(duì)個(gè)體，接著把處理組和對(duì)照組匹配成功的對(duì)象篩選出來，最后運(yùn)用條件logistic模型估計(jì)配比后的處理效應(yīng)并與配比前的處理效應(yīng)進(jìn)行比較。具體如下：

ssc install psmatch2//聯(lián)網(wǎng)條件下安裝psmatch2程序包，版本要求Stata11.0及以上

gen temp=runiform()

sort temp//以上兩步是對(duì)觀測值進(jìn)行隨機(jī)排序

psmatch2 x6 x1 x2 x3 x4 x5 x7 x8 x9,out(y)logit neighbor(1)common caliper(0.03)ties noreplacement//計(jì)算傾向評(píng)分并為處理組找到相應(yīng)的對(duì)照，neighbor是指與處理組匹配的對(duì)照的個(gè)數(shù)，本文以最常用的1∶1配比法為例，因而是neighbor(1)，若要做1∶2配比則可相應(yīng)地變?yōu)閚eighbor(2)，以此類推；caliper是本例所用的卡鉗值,noreplacement是指不重復(fù)利用對(duì)照

sort _n1//對(duì)_n1變量進(jìn)行排序

pstest,both//配比前后各協(xié)變量的均衡性比較

psgraph//以圖片形式展示匹配結(jié)果

save “C:UsersadminDesktopPSMpscore.dta”//至此，傾向評(píng)分的計(jì)算及為每個(gè)處理組尋找合適的對(duì)照已完成，保存一份算出傾向評(píng)分的數(shù)據(jù)庫pscore.dta

keep if x6==1//保留吸煙組，刪除不吸煙組

drop if _support==0//刪除未配比成功的

keep _n1//保留不吸煙組的id編號(hào)

duplicates drop _n1,force//刪除重復(fù)的對(duì)照

ren _n1 _id//重命名_n1為_id以便后續(xù)合并作為關(guān)鍵變量

save “C:UsersadminDesktopPSM\_idduizhao.dta”//選出對(duì)照組的id號(hào)并保存至相應(yīng)文件夾

use “C:UsersadminDesktopPSMpscore.dta”,clear

//打開之前保存的pscore.dta數(shù)據(jù)庫

joinby _id using “C:UsersadminDesktopPSM\_idduizhao.dta”,unmatched(none)

//與對(duì)照組的_idduizhao.dta數(shù)據(jù)庫合并，篩選出對(duì)照的數(shù)據(jù)庫

sort _id//對(duì)_id進(jìn)行排序

egen match=fill(1 2)//根據(jù)id大小給對(duì)照組的每個(gè)觀測值編號(hào)，并生成一個(gè)新變量match

save “C:UsersadminDesktopPSMduizhao.dta”//保存篩選出的對(duì)照數(shù)據(jù)庫

use “C:UsersadminDesktopPSMpscore.dta”,clear

//打開之前保存的pscore.dta數(shù)據(jù)庫

keep if x6==1//保留吸煙組，刪除不吸煙組

drop if _support==0//刪除未配比成功的

duplicates drop _n1,force//刪除重復(fù)利用對(duì)照的研究對(duì)象

sort _id//對(duì)_id進(jìn)行排序

egen match=fill(1 2)//根據(jù)id大小給處理組的每個(gè)觀測值編號(hào)，并生成一個(gè)新變量match

append using “C:UsersadminDesktopPSMduizhao.dta”

//將處理組與對(duì)照組合并，該數(shù)據(jù)庫即為傾向評(píng)分1∶1配比后的數(shù)據(jù)庫

結(jié) 果

1．安裝實(shí)現(xiàn)傾向評(píng)分配比法的程序包psmatch2：出現(xiàn)all files already exist and are up to date提示該程序包已經(jīng)安裝成功。

2．計(jì)算傾向評(píng)分并為處理組找到相應(yīng)的對(duì)照。新生成的幾個(gè)關(guān)鍵變量解釋如下：temp起到對(duì)所有觀測對(duì)象隨機(jī)排序的作用；_pscore是每個(gè)觀測對(duì)象對(duì)應(yīng)的傾向評(píng)分；_id是針對(duì)每個(gè)觀測對(duì)象的唯一編號(hào)(以傾向評(píng)分的大小來排序，本例共有235個(gè)樣本，故_id也是從1編到235)；_treated表示某個(gè)研究對(duì)象是否屬于處理組(treated代表處理組，untreated代表對(duì)照組)；_support表示某個(gè)研究對(duì)象是否被成功匹配(on support 表示匹配成功，off support表示未匹配成功)；_n1表示處理組被匹配到的對(duì)照對(duì)象的_id(如圖1第一行表示的是_id號(hào)為157的研究對(duì)象與_id號(hào)為4的研究對(duì)象配對(duì)，若命令中neighbor選項(xiàng)內(nèi)填的是2，還會(huì)生成_n2，以此類推)；_pdif表示配比的一組觀察對(duì)象概率值的差(圖1)。

圖1 傾向評(píng)分結(jié)果

3．傾向評(píng)分1∶1配比結(jié)果圖示化：Treated代表吸煙者，Untreated代表非吸煙者；On support表示匹配成功，Off support表示未匹配成功,本例有12個(gè)研究對(duì)象未匹配成功(圖2)。

4．傾向評(píng)分配比前后處理組與對(duì)照組在年齡、性別、鼻咽癌家族史、慢性鼻炎史、職業(yè)接觸有害物質(zhì)、飲茶、長期鍛煉、生活工作壓力等方面的分布差異比較，配比后，各協(xié)變量全部均衡可比(表2)。

圖2 傾向評(píng)分1∶1配比結(jié)果

變量配比前處理組(n=79)對(duì)照組(n=156)P值配比后處理組(n=45)對(duì)照組(n=45)P值年齡(x±s)49.16±11.8047.04±12.390.2148.72±11.0850.70±10.800.30性別*<0.0011.00 否6(7.59)65(41.67)5(8.96)6(8.96) 是73(92.41)91(58.33)61(91.04)61(91.04)鼻咽癌家族史*0.070.83 否58(73.42)130(83.33)53(79.10)52(77.61) 是21(26.58)26(16.67)14(20.90)15(22.39)慢性鼻炎史*0.080.40 否60(75.95)133(85.26)35(82.09)51(76.12) 是19(24.05)23(14.74)12(17.91)16(23.88)職業(yè)接觸有害物質(zhì)*0.550.44 否58(73.42)120(76.92)46(68.66)50(74.63) 是21(26.58)36(23.08)21(31.34)17(25.37)飲茶*0.060.73 否36(45.57)91(58.33)30(44.78)32(47.76) 是43(54.43)65(41.67)37(55.22)35(52.24)長期鍛煉*0.581.00 否52(65.82)97(62.18)42(62.69)42(62.69) 是27(34.18)59(37.82)25(37.31)25(37.31)生活工作壓力*0.311.00 否50(63.29)109(69.87)42(62.69)42(62.69) 是29(36.71)47(30.13)25(37.31)25(37.31)

*：數(shù)據(jù)為例數(shù)及構(gòu)成比

討 論

近年來，Stata軟件以其簡單易懂和操作靈活、可以編寫do文件或者直接通過一條命令實(shí)現(xiàn)數(shù)據(jù)處理等功能得到了許多用戶的青睞，尤其在中國，有越來越多的人開始學(xué)習(xí)Stata軟件。對(duì)于傾向評(píng)分配比法，雖然在1983年就被提出，但直到21世紀(jì)初它才開始在國外盛行，在中國被廣泛運(yùn)用更是到了2010年之后[12]。由于Stata最初開發(fā)時(shí)并沒有設(shè)計(jì)傾向評(píng)分配比的程序，因而在Stata中沒有現(xiàn)成模塊，這也導(dǎo)致了很多用戶不知道該如何通過Stata來實(shí)現(xiàn)傾向評(píng)分配比法。因此，本文首先簡要介紹了傾向評(píng)分配比法的基本原理，然后結(jié)合實(shí)際例子演示如何在Stata中安裝傾向評(píng)分配比程序包及如何計(jì)算每個(gè)觀察對(duì)象的傾向評(píng)分，再者演示了如何將匹配到的處理和非處理對(duì)象進(jìn)行1:1配對(duì)并刪除了未匹配和重復(fù)出現(xiàn)的對(duì)象，最后對(duì)配比前后各協(xié)變量的均衡性做了比較，從中發(fā)現(xiàn)Stata能夠較為直觀快速地實(shí)現(xiàn)傾向評(píng)分配比法，并能一定程度上控制非研究混雜因素的干擾。

本次研究結(jié)果發(fā)現(xiàn)，經(jīng)過傾向評(píng)分配比之后，納入模型的協(xié)變量在處理組與對(duì)照組的均衡性均有一定程度的提高，其中性別在配比前差異有統(tǒng)計(jì)學(xué)意義(P<0.05)，經(jīng)過配比后差異無統(tǒng)計(jì)學(xué)意義(P>0.05)，且經(jīng)過配比后處理組與對(duì)照組在年齡、鼻咽癌家族史、慢性鼻炎史、飲茶、長期鍛煉、生活工作壓力等方面的均衡性均有所提高，表明傾向評(píng)分配比能夠起到控制混雜的作用[13-14]。但是，由于傾向評(píng)分配比能夠平衡的是可觀察的混雜因素，對(duì)可能存在的潛在混雜因素卻無能為力，因此傾向評(píng)分配比法只能作為現(xiàn)有條件下控制混雜影響的較優(yōu)途徑[15-16]。

此外，由于本次研究的目的是為了介紹如何在Stata中實(shí)現(xiàn)傾向評(píng)分配比法，所以選取的數(shù)據(jù)庫樣本量較小，其中處理組和對(duì)照組的數(shù)量較為接近，導(dǎo)致最后匹配的成功率較低且出現(xiàn)了重復(fù)利用對(duì)照的現(xiàn)象，經(jīng)過處理后損失了較大一部分樣本量，這有可能會(huì)對(duì)實(shí)際的研究結(jié)果造成影響。因此，在未來的研究中，研究者可以進(jìn)行敏感性分析[17]、加大對(duì)照組的樣本數(shù)或適當(dāng)放寬匹配的精度來估計(jì)、減少其對(duì)結(jié)果的影響。

[1] Rosenbaum PR,Rubin DB.The central role of the propoensity score in observational studies for causal effects.Biometrika,1983,70(1):41-55.

[2] 張華,李楠,趙一鳴.傾向評(píng)分法在臨床研究中的應(yīng)用.中華兒科雜志,2015,53(6):480-480.

[3] 陸俊,黃昌明,鄭朝輝,等.腹腔鏡根治性全胃切除術(shù)治療老年原發(fā)性胃癌患者的傾向評(píng)分配比預(yù)后分析.中華消化外科雜志,2016,15(3):221-227.

[4] 張亮,李嬋娟,夏結(jié)來,等.傾向得分區(qū)間匹配法用于非隨機(jī)對(duì)照試驗(yàn)的探索與研究.中國衛(wèi)生統(tǒng)計(jì),2012,29(1):53-57.

[5] 李智文,李宏田,張樂,等.用SPSS宏程序?qū)崿F(xiàn)觀察對(duì)象的傾向評(píng)分配比.中國衛(wèi)生統(tǒng)計(jì),2011,28(1):89-90.

[6] 李智文,任愛國.傾向評(píng)分法在SAS軟件中的實(shí)現(xiàn).中國生育健康雜志,2010,21(5):320-320.

[7] Rosenbaum PR.Optimal matching for observational studies.Journal of the American Statistical Association.1989,84(408):1024-1032.

[8] Austin PC.A critical appraisal of propensity-score matching in the medical literature between 1996 and 2003.Statistics in Medicine,2008,27(12):2037-2049.

[9] Austin PC,Mamdani MM.A comparison of propensity score methods:a case-study estimating the effectiveness of post-AMI statin use.Statistics in Medicine,2006,25(12):2084-2106.

[10]Wang Y,Cai H,Li C,et al.Optimal caliper width for propensity score matching of three treatment groups:a Monte Carlo study.Plos One，2013,8(12):e81045.

[11]方積乾.衛(wèi)生統(tǒng)計(jì)學(xué).人民衛(wèi)生出版社,2012.

[12]王永吉,蔡宏偉,夏結(jié)來,等.傾向指數(shù)第一講傾向指數(shù)的基本概念和研究步驟.中華流行病學(xué)雜志,2010,31(3):347-348.

[13]王超,吳騁,許金芳,等.傾向性評(píng)分匹配法在不良反應(yīng)信號(hào)檢測中的應(yīng)用.中國衛(wèi)生統(tǒng)計(jì),2012,29(6):855-858.

[14]李智文,張樂,劉建蒙,等.傾向評(píng)分配比在流行病學(xué)設(shè)計(jì)中的應(yīng)用.中華流行病學(xué)雜志,2009,30(5):514-517.

[15]王永吉,蔡宏偉,夏結(jié)來,等.傾向指數(shù)第三講應(yīng)用中的關(guān)鍵問題.中華流行病學(xué)雜志,2010,31(7):823-825.

[16]Korhonen P,Heintjes EM,Williams R,et al.Pioglitazone use and risk of bladder cancer in patients with type 2 diabetes:retrospective cohort study using datasets from four European countries.BMJ,2016,354:i3903.

[17]Chiu M,Rezai MR,Maclagan LC,et al.Moving to a Highly Walkable Neighborhood and Incidence of Hypertension:A Propensity-Score Matched Cohort Study.Environmental Health Perspectives,2016,124(6):754-760.

1.華中科技大學(xué)同濟(jì)醫(yī)學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(430030) 2.福建醫(yī)科大學(xué)公共衛(wèi)生學(xué)院

△ 通信作者：潘安，E-mail:panan@hust.edu.cn

鄧 妍)