天津醫(yī)科大學(xué)流行病與衛(wèi)生統(tǒng)計學(xué)教研室(300070)

張洪璐 于澤洋 宋德勝 黃慧杰 李長平△

傾向評分(propensity score，PS)由Rosenbaum和Rubin于1983年提出，是一種對觀察性研究、非隨機(jī)化干預(yù)性研究或隨機(jī)化失敗的研究中混雜因素進(jìn)行校正的方法[1]，廣泛應(yīng)用于經(jīng)濟(jì)學(xué)、醫(yī)學(xué)等領(lǐng)域[2]。本文主要介紹三種傾向評分分析方法的SAS 9.4軟件實現(xiàn)，分別是傾向評分匹配、分層和加權(quán)，這三種方法可有效均衡混雜因素的組間分布差異，從而減少混雜偏倚，獲得干預(yù)效應(yīng)的有效估計[3]。既往使用SAS軟件進(jìn)行傾向評分分析時，需要較為復(fù)雜的程序或宏程序才可完成，代碼較為繁瑣[4-5]。本研究介紹的PSMATCH過程代碼簡潔，功能強大，可快速實現(xiàn)傾向評分匹配、分層和加權(quán)，為研究者在傾向評分分析中采用PSMATCH過程提供參考。

原理與方法

1.傾向評分分析的基本原理

傾向評分是指在給定一組協(xié)變量xi的情況下，每個個體i(i=1，…，N)被分配到干預(yù)組(Zi=1)的條件概率[1]：

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

xi是一組協(xié)變量，代表著許多維度，而e(xi)是xi的一個函數(shù)，將xi的所有維度簡化為一個單維的值，起到了“降維”的作用，e(xi)能反映xi的綜合作用。估計傾向評分的模型有很多種，如logistic回歸、probit模型、判別分析以及神經(jīng)網(wǎng)絡(luò)等。其中l(wèi)ogistic回歸是最常用且簡便的方式。

基于xi估計出e(xi)后，根據(jù)傾向評分選擇合適的傾向評分分析方法，如傾向評分匹配、分層、加權(quán)等。由于具有相同或相似傾向評分的研究對象有相同的概率接受干預(yù)，因此根據(jù)傾向評分匹配、分層、加權(quán)后，可以認(rèn)為進(jìn)入模型的研究對象的組間分布均衡，達(dá)到了“事后隨機(jī)化”的效果，可獲得干預(yù)效應(yīng)的無偏估計。

2.方法簡介

(1)傾向評分匹配

常見的兩類匹配方法是貪婪匹配和最佳匹配。貪婪匹配是一類“個體最優(yōu)”的算法，即每個干預(yù)組成員在對照組中尋找距離最近者進(jìn)行匹配，一旦配對成功就不再調(diào)整，因此匹配結(jié)果可能不是全局最優(yōu)的，包括馬氏距離匹配、最近鄰匹配和卡尺匹配等。最佳匹配是一種“整體最優(yōu)”的算法，即采用最小化總距離的方式對干預(yù)組成員和對照組成員進(jìn)行匹配，配對成功后可以調(diào)整，包括固定比例匹配、可變匹配和完全匹配[6]。

①馬氏距離匹配

對于干預(yù)組成員i和對照組成員j，其馬氏距離定義為：

d(i，j)=(xi-xj)TC-1(xi-xj)

其中xi和xj分別是i和j的匹配變量取值；C是對照組匹配變量的協(xié)方差矩陣，也可以采用其他定義方式。

②最近鄰匹配

其中ei和ej分別是干預(yù)組和對照組成員的傾向評分值。

③卡尺匹配

|ei-ej|<ε，j∈C

卡尺匹配是指干預(yù)組和對照組個體的傾向得分差值在預(yù)先設(shè)定的某個范圍內(nèi)(即卡尺ε)才會進(jìn)行下一步匹配，這樣的約束可以在一定程度上縮短總距離。通?？ǔ叽笮≡O(shè)置為ε≤0.25σ，σ是由樣本估計的傾向評分值的標(biāo)準(zhǔn)差。

④固定比例匹配

包括成對匹配(1∶1)和一對多匹配(1∶k，k≥2)。在不同的匹配組中，對照組成員的數(shù)量是固定的。

⑤可變匹配

一個干預(yù)組成員和一個或多個對照組成員進(jìn)行匹配。在不同的匹配組中，對照組成員的數(shù)量是變化的。

⑥完全匹配

每個干預(yù)組成員和一個或多個對照組成員進(jìn)行匹配，同時，每個對照組成員可以和一個或多個干預(yù)組成員進(jìn)行匹配。在不同的匹配組中，干預(yù)組和對照組成員的數(shù)量都是變化的。

(2)傾向評分分層

將所有研究對象分成具有相似傾向評分的層，從而均衡各層內(nèi)干預(yù)組和對照組個體的差異。研究表明，使用傾向評分將所有研究對象分為5層可以消除大約90%由測量的混雜因素引起的偏差[8]。

(3)傾向評分加權(quán)

傾向評分加權(quán)旨在對干預(yù)組和對照組成員進(jìn)行加權(quán)，使其能夠代表總體。PSMATCH過程提供三種計算權(quán)重的方法。

①逆概率加權(quán)(inverse probability of treatment weighting，IPTW)

采用第i個觀測的傾向評分e(xi)計算IPTW權(quán)重Wi：

②穩(wěn)定的逆概率加權(quán)(stabilized inverse probability of treatment weighting)

其中，Pt=Nt/(Nt+Nc)，Nt為干預(yù)組人數(shù)，Nc為對照組人數(shù)。

③概率比加權(quán)

采用第i個觀測的傾向評分e(xi)計算概率比權(quán)重Wi：

3.均衡效果的評價指標(biāo)

在完成上述傾向評分分析后，需判斷混雜因素在組間的分布是否均衡。PSMATCH過程能夠輸出兩種常用評價指標(biāo)(標(biāo)準(zhǔn)化均值差值和方差比)和圖形(條形圖、箱式圖等)的結(jié)果。

(1)標(biāo)準(zhǔn)化均值差值(standardized mean difference)

(2)方差比(variance ratio)

研究表明，d≤0.25，0.5≤VR≤2時，組間均衡效果較好[6]。而更嚴(yán)格的標(biāo)準(zhǔn)要求d≤0.1[10]。

4.其他概念

(1)共同支持域(common support region)

干預(yù)組和對照組傾向評分取值范圍的交集[1]。

(2)干預(yù)效應(yīng)

在一項評估某種干預(yù)措施效果的觀察性研究中，研究對象通常有兩種可能的結(jié)局[1，7]：

Y(1)：干預(yù)組個體的結(jié)局

Y(0)：與之情況相同的對照組個體的結(jié)局

將干預(yù)效應(yīng)定義為Y(1)-Y(0)，則平均干預(yù)效應(yīng)(average treatment effect，ATE)為：

ATE=E(Y(1)-Y(0))

干預(yù)組的平均干預(yù)效應(yīng)(average treatment effect for the treated，ATT)為：

ATT=E(Y(1)-Y(0)|Zi=1)

實例分析與SAS實現(xiàn)

SAS 9.4軟件是一種模塊化的大型應(yīng)用軟件系統(tǒng)，在最新的SAS/STAT15.1模塊中PSMATCH過程可以實現(xiàn)傾向得分分析。本研究通過SAS 9.4模擬生成一個數(shù)據(jù)集drugs，其中包括：患者編號(Patient-ID)、用藥情況[Drug：Drug_X為干預(yù)(n=113)、Drug_A為對照(n=373)]、性別(Gender，1=男，2=女)、年齡(Age，連續(xù)變量)、體質(zhì)量指數(shù)(BMI，連續(xù)變量)。在進(jìn)行不同用藥情況的組間比較時，首先應(yīng)均衡組間協(xié)變量的分布差異，本研究分別采用傾向評分匹配、分層和加權(quán)對不同用藥情況的各組間協(xié)變量的分布差異進(jìn)行了調(diào)整，并對三種方法進(jìn)行了比較。

1.傾向評分匹配

(1)SAS代碼

proc psmatch data=drugs region=cs；

class Drug Gender；

psmodel Drug(Treated=′Drug_X′)=Gender Age BMI；

match method=optimal(k=1) exact=Gender distance=lps caliper=0.25 weight=none；

assess lps allcov / plots=(barchart boxplot)；

run；

(2)主要程序及其他選項說明

CLASS語句指定分類變量。PSMODEL語句指定計算傾向評分的logistic回歸模型，TREATED選項指定干預(yù)組的值，可以用引號指定干預(yù)組的值，如本例中的′Drug_X′，也可以用FIRST或LAST分別指定第一個或最后一個類別為干預(yù)組，默認(rèn)情況下，TREATED=FIRST。

MATCH語句將干預(yù)組和對照組成員進(jìn)行匹配。METHOD選項指定匹配方法：①METHOD=GREEDY表示貪婪匹配，其中K選項指定匹配對照組成員的數(shù)量，order選項指定干預(yù)組成員的匹配順序；②METHOD=FULL表示完全匹配，其中KMAX和KMEAN分別指定每個干預(yù)組成員匹配對照組成員的最大數(shù)量和平均數(shù)量，NCONTROL和PCTCONTROL分別指定要匹配的對照組成員的數(shù)量和百分比；③METHOD=OPTIMAL表示最佳固定比例匹配，K選項同上。EXACT選項指定需要精確匹配的分類變量；④METHOD=VARRAITO表示可變匹配，其中KMIN指定每個干預(yù)組成員匹配對照組成員的最小數(shù)量，KMAX、KMEAN、NCONTROL和PCTCONTROL選項同上。DISTANCE選項指定匹配時要比較的距離類型，DISTANCE=LPS/PS/MAH分別表示傾向評分的logit變換/傾向評分/馬氏距離，可附加其他選項選擇計算馬氏距離的變量。CALIPER選項指定匹配的卡尺大小，可以附加MULT選項(指定卡尺乘數(shù)項)或MAHDISTANCE選項(上述DISTANCE=MAH時)。WEIGHT選項指定權(quán)重的類型：①WEIGHT=ATEWGT/MATCHATEWGT表示輸出ATE權(quán)重；②WEIGHT=ATTWGT/MATCHATTWGT/MATCHWGT表示輸出ATT權(quán)重；③WEIGHT=EQUAL/NONE表示兩組權(quán)重相同。

ACCESS語句用于評估組間差異：①ALLCOV選項對PSMODEL語句中指定的協(xié)變量差異進(jìn)行評估；②LPS選項對傾向評分的logit變換的差異進(jìn)行評估；③PS選項對傾向評分的差異進(jìn)行評估；④VAR選項對指定變量的差異進(jìn)行評估，這些變量可以是PSMODEL語句中沒有指定的變量。此外，可以在“/”后添加評估選項：①PLOTS選項指定評估變量均衡性的圖，如條形圖(BARCHART)、箱式圖(BOXPLOT)等；②STDBINVAR選項指定是否對標(biāo)準(zhǔn)化均值差值圖表中的二分類變量進(jìn)行標(biāo)化；③STDDEV選項指定標(biāo)準(zhǔn)化均值差值計算中使用的標(biāo)準(zhǔn)差類型；④VARINFO選項可顯示干預(yù)組和對照組的變量信息。

(3)主要結(jié)果解釋

圖1為匹配后各個變量的分布情況，最終的匹配集中有113對匹配對，匹配后各個變量的分布較均衡(圖1)。表1為匹配后的標(biāo)準(zhǔn)化均值差值。匹配后，干預(yù)組和對照組之間的標(biāo)準(zhǔn)化均值差值均有減少，傾向評分的logit變換值減少了98.35%，年齡減少了95.89%，BMI減少了83.19%，性別減少了100.00%。

圖1 匹配后各變量的分布

此時，各變量干預(yù)組和對照組之間的標(biāo)準(zhǔn)化均值差值分別為0.01375、0.02766、0.06461、0，均小于0.25，方差比分別為1.0155、1.1262、1.1967、1.0000，均處于0.5～2之間(表1)。

表1 匹配后標(biāo)準(zhǔn)化均值差值(干預(yù)-對照)

2.傾向評分分層

(1)SAS代碼

proc psmatch data=drugs region=allobs；

class Drug Gender；

psmodel Drug(Treated=′Drug_X′)=Gender Age BMI；

strata nstrata=5 key=treated stratumwgt=total；

assess ps var=(Gender Age BMI) / varinfo plots=(barchart cdfplot)；

run；

(2)程序說明

同樣調(diào)用PSMATCH過程。STRATA語句將所有觀測分層，NSTRATA選項指定層數(shù)。KEY 選項指定各層的構(gòu)成，KEY=TOTAL/TREATED分別表示各層包含大致相同數(shù)量的成員或干預(yù)組成員。STRATUMWGT指定層權(quán)重類型，STRATUMWGT=TOTAL表示以層中所有成員比例作為權(quán)重，用于估計ATE；STRATUMWGT=TREATED表示以層中干預(yù)組成員比例為權(quán)重，用于估計ATT。其他命令與匹配相同。

(3)主要結(jié)果解釋

表2為分層后SAS的輸出結(jié)果。分層后，干預(yù)組和對照組之間的標(biāo)準(zhǔn)化均值差值均有減少，傾向評分減少了89.47%，年齡減少了83.79%，BMI減少了74.88%，性別減少了87.45%。此時，各變量干預(yù)組和對照組之間的標(biāo)準(zhǔn)化均值差值分別為0.08121、-0.10922、0.09656、-0.00627，均小于0.25，方差比分別為0.8678、0.7981、0.8769、0.9999，均處于0.5～2之間(表2)。

表2 分層后標(biāo)準(zhǔn)化均值差值(干預(yù)-對照)

3.傾向評分加權(quán)

(1)SAS代碼

proc psmatch data=drugs region=allobs；

class Drug Gender；

psmodel Drug(Treated=′Drug_X′)=Gender Age BMI；

psweight weight=atewgt nlargestwgt=6；

assess lps var=(Gender Age BMI)/ varinfo plots=(barchart boxplot(display=(lps BMI)) wgtcloud)；

run；

(2)程序說明

同樣調(diào)用PSMATCH過程。PSWEIGHT語句基于傾向評分計算各個觀測的權(quán)重。WEIGHT選項指定權(quán)重的類型，WEIGHT=ATEWGT表示使用逆概率權(quán)重，可以在ATEWGT選項中附加STABILIZE選項，STABILIZE=YES表示使用穩(wěn)定的逆概率權(quán)重，STABILIZE=NO表示使用逆概率權(quán)重，用于估計ATE，WEIGHT=ATTWGT表示使用概率比權(quán)重，用于估計ATT。NLARGESTWGT選項可顯示權(quán)重最大的觀測。其他命令與匹配相同。

(3)主要結(jié)果解釋

表3為加權(quán)后SAS的輸出結(jié)果。加權(quán)后，干預(yù)組和對照組之間的標(biāo)準(zhǔn)化均值差值均有減少，傾向評分的logit變換值減少了84.85%，年齡減少了82.21%，BMI減少了89.14%，性別減少了25.21%。此時，各變量干預(yù)組和對照組之間的標(biāo)準(zhǔn)化均值差值分別為0.12636、-0.11985、0.04175、0.03735，均小于0.25，方差比分別為0.6874、0.7269、0.8957、1.0054，均處于0.5～2之間(表3)。

表3 加權(quán)后標(biāo)準(zhǔn)化均值差值(干預(yù)-對照)

討 論

本研究通過調(diào)用PSMATCH過程，使用相對簡單的SAS代碼，實現(xiàn)了傾向評分匹配、分層和加權(quán)。在PSMATCH過程中，通過添加不同語句和選項，可以指定共同支持域、分類變量、計算傾向評分的logistic回歸模型以及評估組間差異的方法。在MATCH語句中，可以指定匹配方法、需精確匹配的分類變量、距離類型、卡尺寬度及權(quán)重類型。在STRATA語句中，可以指定層數(shù)、各層構(gòu)成及權(quán)重類型。在PSWEIGHT語句中，可以指定權(quán)重類型。

從本數(shù)據(jù)的結(jié)果可知，在匹配、分層、加權(quán)后，各個變量的標(biāo)準(zhǔn)化均值差值都低于0.25，方差比都在0.5～2之間，可以認(rèn)為這三種方法在一定程度上均衡了組間差異。但以更嚴(yán)格的標(biāo)準(zhǔn)來說，僅在匹配后，各變量的標(biāo)準(zhǔn)化均值差值都低于0.1，而分層和加權(quán)后年齡的標(biāo)準(zhǔn)化均值差值均大于0.1，均衡的效果欠佳。這與Schroeder K[4]等人的研究結(jié)果一致，可能與匹配后總樣本量的改變有關(guān)。匹配過程中由于匹配條件導(dǎo)致的樣本剔除使得納入分析的匹配對間差異較小，而分層和加權(quán)過程中并未改變納入分析的總樣本量。

基于上述分析結(jié)果，本研究可以得出以下三個主要結(jié)論：

1.PSMATCH過程可以實現(xiàn)傾向評分匹配、分層和加權(quán)，使用的SAS代碼非常簡單，并可以根據(jù)研究需要添加多種語句和選項來估計ATE或ATT。

2.PSMATCH過程可以在匹配、分層和加權(quán)后，使用多種方式評估混雜因素在組間分布的差異。就本數(shù)據(jù)而言，這三種方法均可以均衡混雜因素在組間分布的差異。

3.將三種方法進(jìn)行比較時，匹配對組間差異的均衡效果最佳，分層和加權(quán)次之。但匹配造成了樣本量損失，所以要求較大的樣本例數(shù)和共同支持域。因此建議研究者在選擇方法時先測試多種傾向評分方法，再根據(jù)研究目的選擇合適的方法。