帶有不可觀測混雜因素的二值結果的處理效應的估計

長春中醫藥大學藥學院(130117) 韓曦英 孫 健 李秀昌

醫藥學研究中常常要研究對于某種處理(treatment)對二值結局變量的效應(effect)的大小。如研究某種藥物的副作用時，是否服用該藥物為兩值處理變量(服用該藥物，T=1；未服用該藥物，T=0)，其結局變量為兩值變量(有副作用Y=1；沒有副作用Y=0)，等等[1-2]。但是在觀察研究中往往存在著混雜因素，最常見的如年齡、性別、文化程度、認知水平等等。混雜因素既與處理因素相關又與結局變量之間存在相關關系,導致處理因素與結局變量之間的真實關系受到干擾，即處理效應的估計變得不準確。如何控制混雜因素是研究者在觀察性研究中主要關注的問題[3]。

當混雜因素是可以完全觀測時，可以通過加權、匹配、分層等方法將可觀測的混雜進行均衡，消除混雜偏倚。目前廣泛使用的是傾向得分(propensity score)法[4-7]。如別玉坤[8]等利用傾向得分匹配法對兩種方法治療膽總管結石的療效進行評價；KD SONG[9]等應用傾向得分匹配法比較反復肝切除術與射頻(RF)消融術治療復發性肝細胞癌(HCC)的長期療效；Takuma[10]等采用傾向得分匹配法比較符合米蘭標準肝癌病人的肝動脈化療栓塞聯合射頻消融術與手術切除等等。但是有些混雜因素是不可觀測的，如認知能力、文化背景、家庭環境影響、家族基因等等。如何通過統計學方法將不可觀測的混雜的效應從處理效應中剝離掉是目前處理效應估計問題中的熱點。

Gao[11]提出一種動態離散模型中處理效應的估計方法，無需假定不可觀測混雜的分布，只要其方差較大，即個體之間的差異較大，就能基于極大似然估計方法給出處理效應較好的估計[12-15]。本文將這種方法使用到帶有不可觀測混雜的兩值結果的處理效應模型中(相當于T=2的兩期動態離散模型)，給出具體的估計量形式。

方 法

1.模型

帶有不可觀測混雜的二值結局變量的處理效應模型由兩個方程構成：

模型中第一個方程為選擇方程，表示第i個個體是否接受了處理(如是否服用了研究藥物)；第二個方程為結局方程，表示第i個個體是否出現了待評價的二值結局，如是否出現了副作用等。

2.估計方法與主要結論

參數估計的框架工作來源于Gao[11]的文章中T=2的情況。其主要引理和定理如下：

由引理1和引理2，可直接獲得以下定理：

由定理2有：

其中

由上述似然函數并不能得到參數的顯示表達解，可以通過R程序利用牛頓算法進行迭代同時得到參數β和λ的極大似然估計量。具體過程可以通過R程序實現。

3.估計性質

估計量與Gao[11]中的估計量具有相同的近似性質。

結 果

1.模擬結果

使用R程序對本文的工作進行模擬。其中x1～N(0,1)，x2=x1+N(0,1)，ε1、ε2～N(0,1)，β=0.5，進行模擬，結果如表1。

表1 模型中不同的αi分布下處理效應參數的模擬結果(n=1000)

2.實例

某中醫院采用針刺配合雷火灸、常規針刺結合電針兩種不同的針灸方法分別作為對照和處理方法治療陽虛寒凝型膝骨關節炎臨床癥狀，對療效進行考察評價。

根據中西醫診斷和入組標準共納入雷火灸組32例、電針組33例，患者共65例，兩組均為每日治療1次，7天為一療程，每療程休息1天，連續4個療程。4療程結束后評價兩組膝骨性關節炎的療效。療效評價指標確定為二值結果變量，即

雖然在進行臨床觀察入組前經檢驗兩組協變量(性別、患病部位、年齡、病程、VAS積分、WOMAC疼痛評分、WOMAC僵硬評分、WOMAC功能活動積分、WOMAC總分)均達到平衡，結果見表2。但臨床研究中可能有不可觀測的混雜因素影響分組結果或者療效結果，如患者以往經驗、學習背景等影響的對于某種療法的認同感等等，這些因素都在本模型中作為不可觀測的混雜因素αi，故可適用于本文所提之模型。

表2 兩組治療前協變量平衡統計分析結果

表3 兩組治療后臨床療效

結 論

如何處理不可觀測的混雜因素在處理效應的估計問題中一直是人們比較關心的問題。特別在觀察性研究中，進行非隨機化的對比研究，不能人為地設置處理因素，受試對象是否接受處理也不是由隨機化確定的。

在這個過程中不可避免地會有可觀測的和不可觀測的混雜因素同處理因素共同對結局變量產生影響。可觀測的混雜因素產生的偏倚已經有較多的手段去處理。而不可觀測的混雜因素，往往是個體異質性，人們通常對其分布狀況是沒有先驗信息，實驗設計和觀察期間也沒有辦法對其進行有效的控制。在統計過程中直接估計出來的難度非常大，因此在處理效應的估計中目前主流的做法是將不可觀測的混雜在估計的過程中去掉。本文采用的處理效應方法在假設個體異質性的離散程度較大的前提下，對帶有不可觀測混雜的兩值結局變量的處理效應進行估計具有較好的估計性質，比較簡便可行。

[1] Rubin DB.Inference and missing data (with discussion).Biometrika,1976(63):581-592.

[2] Rubin DB.Bayesian inference for causal effects:The role of randomization.The Annals of Statistics,1978(6):34-58.

[3] Wooldridge JM.Econometric analysis of cross section and panel data.The MIT press，2002.

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

[5] Luo Z,Gardiner JC,Bradley CJ.Applying Propensity Score Methods in Medical Research:Pitfalls and Prospects.Medical Care Research and Review:formerly Medical Care Review,2010,67(5):528.

[6] Vansteelandt S，Daniel RM.On regression adjustment for the propensity score.Statistics in Medicine,2014(33):4053-4072,

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

[8] 別玉坤,楊成林,曹衛，等.傾向得分匹配法對兩種方法治療膽總管結石的療效再評價.中國普通外科雜志,2014,23(8):1059-1062.

[9] KD Song,HK Lim,,H Rhim,等.肝切除術后復發性肝細胞癌反復肝切除術與射頻消融術的比較:傾向得分匹配研究.國際醫學放射學雜志,2015,38(4):378.

[10]Takuma Y,Takabatake H,Morimoto Y,et al.采用傾向得分匹配法比較符合米蘭標準肝癌病人的肝動脈化療栓塞聯合射頻消融術與手術切除.國際醫學放射學雜志,2014,(1):82.

[11]Gao W,Bergsma W,Yao Q.Estimation for Dynamic and Static Panel Probit Models with Large Individual Effects.Journal of Time,2016.

[12]Arellano M.Discrete choices with panel data .Investigaciones Economicas,2001,27:423-458

[13]Bartolucci F,Farcomeni A.A multivariate extension of the dynamic logit model for longitudinal data based on a latent Markov heterogeneity structure.Journal of the American Statistical Association,2009,104(486):816-831.

[14]Honoré BE,Kyriazidou E.Panel data discrete choice models with lagged dependent variables.Econometrica,2000,68(4):839-874.

[15]Lancaster T.The incidental parameter problem since 1948.Journal of econometrics,2000,95(2):391-413.