貝葉斯傾向性評分在R軟件中的實現*

四川大學華西公共衛生學院流行病與衛生統計學系(610041) 杜春霖 李宓兒 王 菊 李曉松 劉元元

常規的混雜因素校正方法包括分層、匹配或回歸分析，但當混雜因素較多時，均不再適用。傾向性評分(propensity score,PS)是解決此類問題的有效工具，它將個體所有協變量信息綜合為接受某種處理的條件概率，再利用此概率對樣本進行匹配、分層或加權等，可達到“事后隨機化”的效果[1]。2000年Imbens等人定義了廣義傾向性評分(generalized propensity score,GPS)，將傳統的PS理論擴展到處理因素為多分類及連續型變量的情形[2-3]。

隨著PS方法的廣泛應用，逐漸暴露出傳統理論的一些缺陷，如：默認PS值固定，忽略了其不確定性，影響處理效應估計的準確性[4]；采用logistic等回歸模型估計PS時，模型假設有時無法滿足[5]；在估計模型參數時，無法利用先驗信息；沒有處理缺失數據的有效手段等[6]?；谝陨蠁栴}，國外學者提出了貝葉斯傾向性評分(bayesian propensity score,BPS)，可有效彌補傳統方法的不足[4,7-10]。國內近年來已有部分學者開始研究BPS相關理論[11-13]，但如何實施BPS分析尚缺乏專門介紹。本文將結合BPS理論介紹如何應用R軟件實現BPS分析。

方法原理

Rubin早在1985年就指出PS實際上是一個隨機概率，在處理分配是強可忽略的前提下，PS是可測協變量的充分統計量，具有應用貝葉斯理論的條件。得益于計算科學的高速發展，貝葉斯統計得以廣泛應用，而BPS理論也逐漸趨于成熟，以下按照BPS理論的時序發展介紹3類BPS方法。

1.“聯合”BPS法(Bayesian Approach with Joint Likelihood of PS and Responses)

Hoshino在2008年首次提出了一種針對一般參數模型的擬貝葉斯估計方法，并與結構方程模型結合，用于處理潛變量模型等復雜問題[7]。McCandless等人和An先后在2009和2010年構造了BPS分層回歸法，BPS匹配法與回歸法[4,8]。以上方法均采用MCMC算法，利用PS條件下處理因素和結局變量的聯合分布似然函數，同時估計出處理效應、協變量和PS前的回歸系數。以McCandless的方法為例，建立兩個logistic回歸模型：

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

(1)

logit(Pr(Y=1|C))=γTC

(2)

“聯合”BPS法將PS估計和處理效應估計兩個獨立步驟對應的似然函數聯立為一組似然函數，此時PS作為一個未知量，其后驗樣本將受到結局模型的“反饋”，這類反饋可能導致處理效應的錯誤估計。Gelman和Corwin等人均指出PS應當只反映研究設計，而不應包含任何處理效應或結局的信息[14-15]。因此，當結局變量和PS關系被錯誤建立時，“聯合”BPS法的估計結果通常不準確。

2.“分半”BPS法(“Half-Bayesian” Approach for Propensity Score)

An在2010年的同一研究中提出了另一種“分半”BPS法可避免“模型反饋”的影響，該法采用貝葉斯logistic回歸模型估計PS，而在利用PS估計處理效應時，則采用傳統的OLS方法[8]。具體過程如下。

(1)構建協變量和處理因素的PS模型用于估計PS，即e(x)：

(3)

(3)構建PS、處理因素和結局變量的廣義線性模型用于估計處理效應，以回歸法為例：

Y=μ+γW+λe(X)+ε

(4)

(5)

(6)

3.“分步”BPS法(Two-step Bayesian Approach for Propensity Score)

為克服“模型反饋”的問題，McCandless等人采用Lun提出的一種“切斷反饋”技術對原方法進行了改進，其核心思想是將PS 模型參數的后驗分布作為協變量納入結局模型[9]。同McCandless等人想法類似，Kaplan和Chen進一步提出“分步”BPS法，該法第一步同“分半”BPS法，第二步則采用貝葉斯結局模型估計處理效應[10]。

處理效應則通過γ的后驗均值估計如下：

E(γ|W,Y,X)=E{E(γ|η,W,Y,X)|W,Y,X}

(7)

(8)

上式又可通過貝葉斯PS模型中η的后驗樣本均值估計得到，則：

(9)

γ的后驗方差可通過下式求得：

Var(γ|W,Y,X)=

E{Var(γ|η,W,Y,X)|W,Y,X}+

Var{E(γ|η,W,Y,X)|W,Y,X}

(10)

(11)

(12)

因此“分步”BPS法的處理效應方差估計值如下：

(13)

從(13)式可以看出處理效應的兩部分變異，分別對應后驗樣本方差的均值，以及后驗樣本均值的方差，同時綜合了PS和處理效應的不確定性。

模 擬

目前，“聯合”BPS法可通過R Studio下載An編寫的IUPS包實現?！胺职搿盉PS法可借助MCMCpack包編程實現，而“分步”BPS法的函數包BayesPSA需通過本地文件下載和安裝(鏈接：http://bise.wceruw.org/publications.html)。

本文數據仿照Chen和Kaplan研究中所用的模擬數據[10]。其中包含3個協變量(x1～N(1,1),x2～Poisson(2),x3～Bernoulli(0.5)) ，真實PS值計算如下：

(14)

通過比較e(x)i和Ui(U～Uniform(0,1))的大小，產生處理因素w(當Ui≤e(x)i時wi=1，Ui>e(x)i時wi=0)，結局變量y按下式生成：

y=0.4x1+0.3x2+0.2x3+0.5w+ε

(15)

其中，ε～N(0,0.1)，本例中樣本量設為100和250。

BPS分析的具體過程如下(R代碼可參考附錄)：

(1)整合實際數據。本文數據通過模擬產生，真實處理效應預設為0.5。

(2)加載函數包，實現BPS分析。其中，IUPS包提供了“聯合”BPS最鄰近匹配和回歸的函數，BayesPSA包可實現“分步”BPS最優匹配和分層法?！胺职搿盉PS法以加權法為例進行了展示，讀者可參考程序實現其他幾種PS利用方式。

(3)將第2步所得主要結果整理為表1。由于“聯合”BPS法中默認采用參數的無信息先驗，為使結果具有可比性，在“分半”BPS法和“分步”BPS法中同樣采用無信息先驗。

表1 BPS分析結果匯總表

從結果可看出“分半”BPS加權法偏倚最小，而“聯合”BPS回歸法的標準誤最小，3種BPS法的處理效應估計效果均隨樣本增大而更優。值得注意的是，此處假設正確建立了PS模型和結局模型，因此“聯合”BPS法的結果很少受到“模型反饋”的影響。對于本文BPS分析的結果，在實際應用中會受到先驗精度、PS利用方式以及模型不確定性等因素的影響，因此并不具備代表性，有待更為系統的模擬研究以探討這3類BPS方法的適用條件。

討 論

PS 方法作為一種處理觀察性研究中混雜偏倚的有力工具，已被廣泛應用于流行病學、心理學和社會學等多個領域。 PS理論從提出到目前已有30年的歷史，BPS作為其最新理論成果，具有廣闊的應用前景。BPS考慮了傾向分值的隨機性，可以更精確地估計傾向分值，并有效利用先驗知識，使結果更加真實可靠。此外，BPS可與復雜模型相結合處理相應的實際數據，尤其適用于小樣本情形[4]。

本文分別介紹了“聯合”BPS法、“分半”BPS法、“分步”BPS法，其中“聯合”法必須是在正確建立結局變量和PS模型的前提下，結果才具有參考價值，而不同學者提出的聯合分布似然函數不盡相同，尚無相關研究論證孰優孰劣。其次，有關先驗信息的設定僅Chen和 Kaplan在“分步”法中進行過討論，主要比較了無信息先驗和精確先驗在不同先驗精度下對結果的影響，其他形式先驗的估計效果有待進一步研究論證[10]。最后，PS模型正確與否和好壞與否極大程度影響到結果的準確性和可靠性，將貝葉斯模型平均(Bayesian Model Averaging,BMA)應用至BPS分析，可一定程度上綜合模型選擇的不確定性，得到更加合理的效應估計[16]。

此外，BPS對于協變量的均衡效果以及BPS在多分類處理因素中的應用均是當前的研究熱點[12,17]，讀者可在實現簡單BPS分析的基礎上，進一步了解和掌握BPS的相關理論，推廣其在實際數據分析工作中的應用。

附錄：

###模擬數據

gendata<-function(size){

#產生自變量

x1<-rnorm(size,1,1);x2<-rpois(size,2);x3<-rbinom(size,1,0.5)

#計算傾向性評分

ps<-exp(0.2*x1+0.3*x2-0.2*x3)/(1+exp(0.2*x1+0.3*x2-0.2*x3))

#產生處理因素

u<-runif(size,0,1);w<-NULL

for(i in 1:size){

w[i]<-ifelse(u[i]<=ps[i],1,0)

}

#產生結局變量

e<-rnorm(size,0,0.1)

y<-0.4*x1+0.3*x2+0.2*x3+0.5*w+e

#產生數據

simudata<-data.frame(y,x1,x2,x3,w)

return(simudata)

}

data1<-gendata(100);data2<-gendata(250)

###聯合BPS法

library(IUPS)

X<-as.matrix(data1$x1,data1$x2,data1$x3)

#BPS最鄰近匹配(“ATE”代表平均處理效應)

bpsm(data1$y,data1$w,X,method=“BPSM”,estimand=“ATE”)

#BPS回歸

bpsr(data1$y,data1$w,X)

###分半BPS加權法

library(MCMCpack)

#通過MCMC產生參數后驗樣本(burnin為“消火”次數，mcmc為迭代次數，thin為間隔長度，tune為Metropolis調整參數，，b0和B0分別對應服從多元正態分布參數的先驗均值和精度)

poster1<-MCMClogit(w～x1+x2+x3,data=data1,burnin=1000,mcmc=10000,thin=10,tune=0.25,b0=0,B0=0)

beta<-poster1

#計算PS

gmodel<-glm(w～x1+x2+x3,data=data1)

mat<-model.matrix(gmodel)

bps<-matrix(rep(0,nrow(beta)*nrow(data1)),nrow=nrow(data1))

for(i in 1:nrow(beta))

{

temp<-mat%*%beta[i,]

bps[,i]<-exp(temp)/(1+exp(temp))

}

#構造加權函數(采用逆處理概率加權)

psweight<-function(y,ps,trt)

{

wt<-rep(0,length(ps))

wt[trt==1]<-1/ps[trt==1]

wt[trt==0]<-1/(1-ps[trt==0])

mod<-lm(y～trt,weight=wt)

return(c(summary(mod)$coef[2,1],summary(mod)$coef[2,2]))

}

#估計處理效應及其標準誤

bwlm<-cbind(rep(0,nrow(beta)),rep(0,nrow(beta)))

for(i in 1:nrow(beta))

{

bwlm[i,]<-psweight(data1$y,bps[,i],data1$w)

}

mean(bwlm[,1])

sqrt(var(bwlm[,1])+mean(bwlm[,2]^2))

###分步BPS法

#計算PS (response和treatment分別用于指定結局變量名，處理變量名，treatment.success.level定義接受處理時所取變量值，vars為協變量名，prior.b0、prior.B0分別為先驗均值和先驗精度，mcmc.burnin為“消火”次數，mcmc.iter為迭代次數，mcmc.thin為間隔長度，mcmc.tune為Metropolis調整參數，method可選“BMA”，即貝葉斯模型平均。)

bps<-bpsa(data1,response='y',treatment='w',treatment.success.level=1,

vars=c('x1','x2','x3'),prior.b0=0,prior.B0=0,mcmc.burnin=1000,

mcmc.iter=10000,mcmc.thin=10,mcmc.tune=0.25,

method=c(“MCMC”))

#BPS最優匹配

bpsa.opt.match(bps,prior.b0=0,prior.B0=0,mcmc.burnin=1000,

mcmc.iter=10000,mcmc.thin=10,mcmc.tune=0.25)

#BPS分層

bpsa.strat(bps,prior.b0=0,prior.B0=0,mcmc.burnin=1000,

mcmc.iter=10000,mcmc.thin=10)