貝葉斯傾向評分方法及其與傳統傾向評分方法的比較*

第二軍醫大學衛生統計學教研室(200433)

張 筱 葉小飛 張新佶 郭曉晶 吳美京 張天一 李 慧 賀 佳△

近年來，由于隨機對照試驗(randomized controlled trials,RCT)通常存在所選人群有限、樣本量較少、隨訪時間較短、價格昂貴等局限性，它的大規模應用受到一定的限制。而觀察性研究(observational study,OBS)通常不對研究人群進行嚴格的限制、樣本量較大、觀察時間較長、研究成本較低、可以收集到足夠的信息觀察到特殊人群(如孕婦、兒童、有合并其他疾病的患者)的暴露情況，也可以對罕見事件的發生做出評價[1]，因而近年來國內外越來越多的研究者采用觀察性研究方法在醫學領域進行大規模人群研究。但觀察性研究不能像RCT那樣采用隨機化設計對研究對象隨機分配，可能出現混雜因素在組間分布不均衡現象，產生混雜偏倚。傳統的控制混雜因素的方法如分層法(當分層數目較多時會產生過度分層的問題)及logistic回歸法(受模型線性假設條件的限制)存在一定的局限性[2]，因此，迫切需要更加有效的方法來均衡不同特征數據之間的差異，更為準確地控制混雜因素的影響，使不同特征的數據間具有可比性，從而更好地挖掘出數據中隱藏的信息，獲得更準確的信號檢測結果。目前，傾向評分法及貝葉斯傾向評分法是解決這一問題的較有力工具。

1.傾向評分法

傾向評分法(propensity score analysis,PSA)作為均衡組間混雜因素的新方法由Rosenbaum和Rubin在1983年首次提出[3]，其基本原理是將多個協變量的影響用一個傾向評分值來表示(相當于降低了協變量的維度)，然后根據傾向評分值進行不同對比組間的分層、匹配或加權，即均衡對比組間協變量的分布，最后在協變量分布均衡的層內或者匹配組中估計處理效應。在大樣本的情況下，經過傾向評分值調整的組間個體，除了處理因素和結局變量分布不同外，其他協變量應該均衡可比，相當于進行了“事后隨機化”，使觀察性數據達到“接近隨機分配數據”的效果[4-5]。

估計傾向性評分值是傾向評分法的第一步，也是核心步驟。其估計的準確與否，直接關系到組間均衡的效果，進而影響到對處理效應的正確推斷。目前，常用的估計傾向評分值的模型有logistic回歸、probit回歸、判別分析以及數據挖掘中的神經網絡、支持向量機、分類與回歸樹、Boosting算法等機器學習方法[6]。判別分析要求協變量服從多元正態分布，而流行病學資料中存在著較多的分類變量，因此該方法在流行病學中較少使用[7]；probit回歸函數表示累積標準正態分布函數的逆函數或反函數，其結果不易解釋，限制了該方法的廣泛使用；數據挖掘的方法因其穩定性差、難以理解、結構復雜等缺陷，較少有研究將其應用于實際數據中[8]。

傳統的傾向評分法即借助logistic回歸模型估計傾向評分值，其應用最為廣泛[9]。采用傳統的logistic回歸模型估計傾向評分值具有模型簡單、容易實現、穩健性好、結果易于解釋等優勢。然而，其在應用中存在的問題也不容忽視：(1)通過logistic回歸模型估計的傾向評分分值與其真實值的偏倚較大，而研究者未考慮到傾向評分值的不準確性對混雜因素組間均衡性的影響，進而影響到處理效應估計的準確性，特別是在傾向評分調整法和匹配法中，這個問題更為突出[10]；(2)傳統的logistic回歸模型估計傾向評分值時，連續型協變量與logit(y)需要滿足線性關系的限制條件，當此條件不被滿足時，結果的準確性將受到影響，而在實際應用中，研究者往往忽略了對此假設條件進行檢驗[11]；(3) 當事件發生數與協變量數之比小于10時，傳統的logistic回歸估計得到的傾向評分值也存在較大的偏倚[12]；(4)估計處理效應時無法利用先驗信息。如果可以利用已有的信息，可有效提高處理效應估計值的精確性和可靠性；(5) 無法很好地處理缺失數據、高維數據等問題。

因此，國外越來越多的研究者開始將貝葉斯統計的思想引入到傾向評分法中，構建貝葉斯傾向評分模型，估計傾向評分值及處理效應，以彌補傳統的傾向評分法無法解決的問題。

2.貝葉斯傾向評分法(Bayesian propensity score analysis，BPSA)

貝葉斯統計是將關于未知參數的先驗信息與樣本信息綜合，并不斷通過樣本數據更新先前認知的統計方法。它采用馬爾可夫鏈蒙特卡羅算法(MCMC)對估計參數的后驗分布進行抽樣，并結合所構建的似然函數對所估計的先前認知(即MCMC法上一次迭代的先驗信息與樣本信息的綜合)不斷地進行修正，最終得到估計參數的一個穩定的后驗分布，根據后驗信息去推斷總體參數。貝葉斯統計是基于總體信息、樣本信息及先驗信息進行的統計推斷。1985年，Rubin等人論證了傾向評分法與貝葉斯思想結合的合理性，首次在專著中提出將貝葉斯統計的思想引入到傾向性評分法中，但并未構建相應的模型[13]。目前，國外學者普遍認為傾向評分模型中的參數是未知的，具有不確定性及隨機性，因此可以構建貝葉斯傾向評分模型，且貝葉斯傾向評分模型可以與貝葉斯因果推斷模型或傳統的概率統計因果推斷模型相結合。圍繞此種觀點，現已提出了多種貝葉斯傾向評分的模型。

Hoshino在2008年提出一種半參數貝葉斯傾向評分模型，并與復雜模型相結合(如結構方程模型)，用于處理潛在變量的影響或解決多組間比較的問題[14]，但因其復雜的數理推理過程限制了模型的廣泛應用。2009年，McCandless、Gustafson、Austin等人構建了另一種貝葉斯傾向評分模型[10]，可以同時估計傾向分值和處理效應。目前，此模型的應用較多。如公式(1)、公式(2)所示，其中公式(1)為結局變量的估計，X代表處理因素，β為處理因素的回歸系數(即處理效應的估計)，Z(C,γ)代表公式(1)中所估計的傾向分值。公式(2)為傾向分值的估計，C代表混雜因素。

logit[Pr(Y=1/X,C)]=βX+ζTg(z(C,γ))

(1)

logit[Pr(Y=1/X,C)]=γTC

(2)

(3)

(4)

(5)

McCandless先應用貝葉斯logistic回歸模型納入候選協變量，利用先驗信息，根據貝葉斯統計的思想，得出每個個體的傾向性評分值，將PS值分為五段，作為潛變量引入貝葉斯回歸模型，并構建混雜因素條件下結局變量與處理因素聯合分布的似然函數。他提出的BPSE算法有以下三個特點：(1)采用無信息先驗分布，假定先驗信息β、γ、ξ方差很大且服從正態分布，見公式(3)、(4)、(5)；(2)對傾向分值進行分層(分層節點的選擇可直接選用五分位值或采用立方樣條法擬合節點)，然后將傾向分值作為潛變量納入處理效應估計的模型中；(3)借助MCMC法，利用所構建聯合分布的似然函數，同時估計出處理效應的回歸系數(β)、協變量的回歸系數(γ)以及傾向分值的回歸系數(ξ)。

如上文所述，McCandless等人構建的貝葉斯傾向評分模型將傾向分值進行分層后再作為協變量納入結局變量的似然函數，并同時估計出傾向分值及處理效應。在其結果解釋時，他將傾向分值的回歸系數(ξ)當作冗余參數。而事實上傾向分值的回歸系數(ξ)可以反映結局變量的估計對傾向評分估計的影響，而上述研究中忽略了此種影響關系。因此，2010年，McCandless等人又探討了兩步進行的貝葉斯回歸調整法，以控制結局變量的估計對傾向分值估計的影響，更加精確地估計傾向評分值[15]。并將其應用于時依性結局變量的數據分析中(如生存分析資料)混雜因素的控制。如公式(6)表示為同時估計模型的概率密度函數，公式(7)表示分兩步進行時傾向分值的概率密度函數。

(6)

(7)

(8)

2012年，David Kaplan、Jianshen Chen等學者認為先前McCandless、An等學者所構建的貝葉斯聯合分布忽略了結局變量對傾向評分分布的影響，進一步假設聯合分布中若不存在結局變量，那么會產生不同的傾向分值的分布，于是提出了分兩步完成的貝葉斯傾向評分法[17]。簡單地說，在David Kaplan的研究中，兩步的BPSA就是兩個獨立的過程，即估計傾向評分值時采用貝葉斯logistic模型(與McCandless、An的估計傾向分值的似然函數相同)，而用最小二乘法(如公式(9)所示，X代表處理因素，β為處理因素的回歸系數即所要估計的處理效應)或貝葉斯回歸估計處理效應。David Kaplan還將他提出的兩步BPSA法與McCandless等人的一步BPSA法進行了比較，研究發現兩步完成的貝葉斯傾向評分分層法所估計的處理效應的方差小于一步完成的貝葉斯傾向評分的分層法，但最優匹配法的結果較接近。此外，貝葉斯傾向評分法與傳統的方法相比較，發現貝葉斯傾向評分法估計的處理效應的方差較大，與McCandless等人的研究結果一致。他認為傳統的傾向評分法估計的處理效應的方差較小，而BPSA提供的處理效應的方差雖然略有增大然而更加準確。

Y=a+βX+ε

(9)

貝葉斯傾向評分法與傳統的logistic傾向評分法對處理效應的估計哪個更加精確呢？這個問題目前尚未有統一的結論。不同方法特點不同：(1) 處理效應的估計受不同傾向評分值的利用方式的影響。如直接選用五分位值將傾向分值分為五層[10]或采用立方樣條法擬合節點分為四層[18]，選用最鄰近匹配法(nearest neighborhood matching)[16]或最優匹配法(optimal full matching)[17]，對結果都有影響。(2)不同貝葉斯傾向評分法選擇的估計結局變量的似然函數不同[10,16,18]。如McCandless 的研究中處理效應的似然函數采用貝葉斯logistic回歸，而An的研究中處理效應的似然函數采用一般線性模型。(3)先驗信息分布的選擇(如有研究選擇有信息先驗分布，也有研究選擇無信息)以及超參數的設定不同。David Kaplan等人通過模擬研究證明先驗信息設定的越準確，處理效應的估計越接近真值。因此，需要基于自己研究的數據特征，對不同學者提出的貝葉斯傾向評分法進行橫向比較研究，篩選最優的貝葉斯傾向評分法。

3.貝葉斯傾向評分模型的實際應用

4.貝葉斯傾向評分法與logistic傾向評分法的比較

現從以下九個方面對傳統的傾向評分法及貝葉斯傾向評分法進行簡要的比較(如表1所示)：(1)PS值估計的準確性 可通過協變量的回歸系數(γ)來反映，貝葉斯傾向評分法可以更加精確地估計γ的可信區間[10,15-16,18]。(2)模型假設 傳統的logistic回歸模型需要滿足對數線性假設的條件，結果才更可靠，而貝葉斯傾向評分模型對此假設條件并不敏感[10]。(3)均衡性 根據國外多項研究結果，與結局相關的變量均為重要的混雜因素，需要對其進行組間均衡性檢驗。而與處理因素強相關、與結局變量弱相關的混雜，對結果的影響不大，可以忽略。傳統的logistic回歸模型主要均衡與處理因素相關的混雜，可能會遺漏掉某些與結局相關的重要的混雜；而貝葉斯傾向評分主要均衡與結局變量相關的混雜[18]。(4)缺失數據 若某個或某幾個協變量存在缺失，logistic回歸等傳統統計方法便無法得到傾向評分值。而貝葉斯傾向評分法允許缺失數據的存在。(5)高維數據 對于高維數據，變量之間可能存在各種各樣的線性及非線性關系或交互作用，貝葉斯傾向評分法在處理這些問題方面有著明顯的優勢。(6)潛變量(latent variable) 也稱為不可測量的變量(unmeasured variable)。貝葉斯傾向評分法可以與復雜模型相結合(如結構方程模型等)，處理潛在變量的問題[14,20]。(7)樣本大小或事件發生數目 當樣本數較小或事件發生數較小時，使用logistic傾向評分法，結果不夠穩定，而貝葉斯傾向評分法在此種情況下有著明顯的優勢[16-17]。(8)軟件應用 傳統logistic回歸模型可在SAS、R、Stata等多種軟件中實現，貝葉斯傾向評分則只能通過R軟件計算。(9)難易程度 logistic回歸模型具有簡單、容易實現等特點。貝葉斯傾向評分法的實現比較復雜，需要預先確定先驗分布(無信息先驗分布、有信息先驗分布)、設定參數的初始值、選取恰當的抽樣方式(Metropolis或Gibbs抽樣)，還需要借助MCMC法才能得以實現。

表1 貝葉斯傾向評分法及傳統的傾向評分法的比較

5.展望

隨著信息化技術的推進，在日常業務中可以通過信息系統收集大量的觀察性數據，如不良反應自發呈報系統(spontaneous reporting system,SRS)[21]、醫院信息系統(HIS)、電子病歷(EMR)等。如何將這些數據有效利用，提供有價值的關于干預因素與結局之間因果關系的“證據”或“線索”，為醫學與政策問題的研究及解決提供巨大的數據支持及循證支持，已成為統計方法學研究中面臨的巨大挑戰。貝葉斯傾向評分法是近年來新提出的一種處理觀察性研究中混雜偏倚的有力工具，它可以有效地利用先前研究或系統累積數據中的大量信息、充分考慮傾向分值的隨機性而更加精確地估計傾向評分值、與復雜模型相結合處理復雜結構的數據以及借助MCMC法快速的估計出各項參數的后驗分布等，因而較傳統的傾向評分法具有廣泛的應用前景。目前，貝葉斯傾向評分法在應用過程中還存在著較多的問題：不同學者提出的聯合分布及似然函數也不盡相同；不同研究中先驗信息的設定方式不同；傾向分值的利用方式也不盡相同。以上因素都會影響到處理效應的估計結果。因此，在對觀察性數據進行研究時，需要充分考慮數據特征及研究目的，選擇最佳的模型進行分析。此外，目前的貝葉斯傾向評分法僅限于兩分類的處理因素及結局變量，對于多分類的處理因素、隨時間變化的處理因素、連續型結局變量等觀察性數據中常見的問題，還需要進一步的研究。

參 考 文 獻

1.Perkins SM,Tu W,Underhill MG,et al.The use of propensity scores in pharmacoepidemiologic research.Pharmacoepidemiol Drug Saf,2000,9(2):93-101.

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

3.Rosenbaum PR,Rubin DB.The central role of the propensity score in observational studies for causal effects.Biometrika,1983,70:41-55.

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

5.李智文,李宏田,張樂.用SPSS宏程序實現觀察對象的傾向評分配比.中國衛生統計,2011,28(1):89-92.

6.Westreich D,Lessler J,Funk MJ.Propensity score estimation:neural networks,support vector machines,decision trees (CART),and meta-classifiers as alternatives to logistic regression.J Clin Epidemiol,2010,63(8):826-33.

7.D′Agostino RB Jr.Propensity score methods for bias reduction in the comparison of a treatment to a non-randomized control group.Stat Med,1998,17(19):2265-81.

8.Setoguchi S,Schneeweiss S,Brookhart MA,et al.Evaluating uses of data mining techniques in propensity score estimation:a simulation study.Pharmacoepidemiol Drug Saf,2008,17(6):546-55.

9.Felix J.A Systematic Review of Propensity Score Methods in the Social Sciences.Multivariate Behavioral Research,2011,46:1,90-118.

10.McCandless LC,Gustafson P,Austin PC.Bayesian propensity score analysis for observational data.Stat Med,2009,28(1):94-112.

11.Bagley SC,White H,Golomb BA.Logistic regression in the medical literature:standards for use and reporting,with particular attention to one medical domain.J Clin Epidemiol,2001,54(10):979-85.

12.Cepeda MS,Boston R,Farrar JT,et al.Comparison of logistic regression versus propensity score when the number of events is low and there are multiple confounders.Am J Epidemiol,2003,158(3):280-7.

13.Rubin DB.The use of propensity scores in applied Bayesian inference.In Bayesian Statistics 2,Bernardo JM,De Groot MH,Lindley DV,Smith AFM (eds).Valencia University Press,North-Holland:Amsterdam,1985,63-72.

14.Hoshino T.A Bayesian propensity score adjustment for latent variable modeling and MCMC algorithm.Computational Statistics & Data Analysis,2008,52,1413-1429.

15.McCandless LC,Douglas IJ,Evans SJ,et al.Cutting feedback in Bayesian regression adjustment for the propensity score.Int J Biostat,2012,6(2):Article.

16.An W.Bayesian propensity score estimators:incorporating uncertainties in propensity scores into causal inference.Sociological Methodology,2010,40,151-189.

17.David K,Jianshen C.A two-step Bayesian approach for propensity score analysis:simulations and case study.Psychometrika,2012,77(3):581-609.

18.McCandless LC,Gustafson P,Austin PC,et al.Covariate balance in a Bayesian propensity score analysis of beta blocker therapy in heart failure patients.Epidemiol Perspect Innov,2009,6(5).

19.McCandless LC,Gustafson P,Austin PC.Code for fitting Bayesian propensity analysis to a toy synthetic dataset [CP/OL].http://www.biomedcentral.com/content/supplementary/1742-5573-6-5-S1.R.

20.McCandless LC,Richardsonand S,Nicky GB.Propensity Score Adjustment for Unmeasured Confounding in Observational Studies.ESRC National Center for Research Methods NCRM Working Paper Series,2,2008.

21.錢維,王超,吳騁,等.運用隨機森林分析藥品不良反應發生的影響因素.中國衛生統計,2013,30(2):209-213.