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

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

隨著PS方法的廣泛應(yīng)用，逐漸暴露出傳統(tǒng)理論的一些缺陷，如：默認(rèn)PS值固定，忽略了其不確定性，影響處理效應(yīng)估計的準(zhǔn)確性[4]；采用logistic等回歸模型估計PS時，模型假設(shè)有時無法滿足[5]；在估計模型參數(shù)時，無法利用先驗信息；沒有處理缺失數(shù)據(jù)的有效手段等[6]。基于以上問題，國外學(xué)者提出了貝葉斯傾向性評分(bayesian propensity score,BPS)，可有效彌補(bǔ)傳統(tǒng)方法的不足[4,7-10]。國內(nèi)近年來已有部分學(xué)者開始研究BPS相關(guān)理論[11-13]，但如何實施BPS分析尚缺乏專門介紹。本文將結(jié)合BPS理論介紹如何應(yīng)用R軟件實現(xiàn)BPS分析。

方法原理

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

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

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

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

(1)

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

(2)

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

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

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

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

(3)

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

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

(4)

(5)

(6)

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

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

處理效應(yīng)則通過γ的后驗均值估計如下：

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法的處理效應(yīng)方差估計值如下：

(13)

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

模 擬

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

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

(14)

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

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

(15)

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

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

(1)整合實際數(shù)據(jù)。本文數(shù)據(jù)通過模擬產(chǎn)生，真實處理效應(yīng)預(yù)設(shè)為0.5。

(2)加載函數(shù)包，實現(xiàn)BPS分析。其中，IUPS包提供了“聯(lián)合”BPS最鄰近匹配和回歸的函數(shù)，BayesPSA包可實現(xiàn)“分步”BPS最優(yōu)匹配和分層法。“分半”BPS法以加權(quán)法為例進(jìn)行了展示，讀者可參考程序?qū)崿F(xiàn)其他幾種PS利用方式。

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

表1 BPS分析結(jié)果匯總表

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

討 論

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

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

此外，BPS對于協(xié)變量的均衡效果以及BPS在多分類處理因素中的應(yīng)用均是當(dāng)前的研究熱點[12,17]，讀者可在實現(xiàn)簡單BPS分析的基礎(chǔ)上，進(jìn)一步了解和掌握BPS的相關(guān)理論，推廣其在實際數(shù)據(jù)分析工作中的應(yīng)用。

附錄：

###模擬數(shù)據(jù)

gendata<-function(size){

#產(chǎn)生自變量

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))

#產(chǎn)生處理因素

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

for(i in 1:size){

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

}

#產(chǎn)生結(jié)局變量

e<-rnorm(size,0,0.1)

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

#產(chǎn)生數(shù)據(jù)

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

return(simudata)

}

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

###聯(lián)合BPS法

library(IUPS)

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

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

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

#BPS回歸

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

###分半BPS加權(quán)法

library(MCMCpack)

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

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))

}

#構(gòu)造加權(quán)函數(shù)(采用逆處理概率加權(quán))

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]))

}

#估計處理效應(yīng)及其標(biāo)準(zhǔn)誤

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分別用于指定結(jié)局變量名，處理變量名，treatment.success.level定義接受處理時所取變量值，vars為協(xié)變量名，prior.b0、prior.B0分別為先驗均值和先驗精度，mcmc.burnin為“消火”次數(shù)，mcmc.iter為迭代次數(shù)，mcmc.thin為間隔長度，mcmc.tune為Metropolis調(diào)整參數(shù)，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最優(yōu)匹配

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)