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

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

·計算機應用·

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

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

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

原理和方法

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

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

(1)

假定Zi和Xi相互獨立，則可利用logit或probit模型計算出每個對象的傾向評分[8]，公式如下：

(2)

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

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

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

(3)

其中Pi為處理組對象的傾向評分，Pj為對照組對象的傾向評分。

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

實例分析

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

表1 各變量命名及賦值

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

操作步驟

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

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

gen temp=runiform()

sort temp//以上兩步是對觀測值進行隨機排序

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

sort _n1//對_n1變量進行排序

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

psgraph//以圖片形式展示匹配結果

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

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

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

keep _n1//保留不吸煙組的id編號

duplicates drop _n1,force//刪除重復的對照

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

save “C:UsersadminDesktopPSM\_idduizhao.dta”//選出對照組的id號并保存至相應文件夾

use “C:UsersadminDesktopPSMpscore.dta”,clear

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

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

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

sort _id//對_id進行排序

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

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

use “C:UsersadminDesktopPSMpscore.dta”,clear

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

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

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

duplicates drop _n1,force//刪除重復利用對照的研究對象

sort _id//對_id進行排序

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

append using “C:UsersadminDesktopPSMduizhao.dta”

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

結 果

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

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

圖1 傾向評分結果

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

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

圖2 傾向評分1∶1配比結果

變量配比前處理組(n=79)對照組(n=156)P值配比后處理組(n=45)對照組(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è)接觸有害物質*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ù)及構成比

討 論

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

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

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

[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] 張華,李楠,趙一鳴.傾向評分法在臨床研究中的應用.中華兒科雜志,2015,53(6):480-480.

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

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

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

[6] 李智文,任愛國.傾向評分法在SAS軟件中的實現(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)計學.人民衛(wèi)生出版社,2012.

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

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

[14]李智文,張樂,劉建蒙,等.傾向評分配比在流行病學設計中的應用.中華流行病學雜志,2009,30(5):514-517.

[15]王永吉,蔡宏偉,夏結來,等.傾向指數(shù)第三講應用中的關鍵問題.中華流行病學雜志,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.華中科技大學同濟醫(yī)學院公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學系(430030) 2.福建醫(yī)科大學公共衛(wèi)生學院

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

鄧 妍)