邊際結構模型基本原理及其應用實例介紹*

田丹平 張 敏 黃淵秀 高 林 董 晶 李 黎 鄧 欣 楊 芳 胡國清

由于觀察性研究無法隨機分配不同處理組的暴露水平，故在暴露變量和結局變量間的因果關聯推斷方面無法獲得與隨機對照試驗同等質量的證據。針對觀察性研究的這一不足，國外學者在過去的幾十年中先后提出了一些新的方法。邊際結構模型(marginal structural models,MSMs)是近二十年中出現的一類新的因果推斷方法，它由哈佛大學的Robins于1999年提出。邊際結構模型主要適用于存在時間依賴協變量與時間依賴暴露變量的觀察性研究。鑒于目前國內絕大多數醫學工作者對邊際結構模型尚不熟悉，本文擬簡要介紹其基本原理，并結合實例說明其應用。

1.基本概念

(1)時間依賴混雜因素(time-dependent confounders)：變量取值隨時間變化而變化的混雜因素。它既影響結局變量，也影響后期的暴露變量。

(2)反事實變量或潛結局變量(counterfactual or potential outcomes)：反事實變量是理論上可能出現但現實卻未被觀測到的變量。假設欲估計某暴露變量A與結局變量Y的因果關聯，A所有可能的取值有J+1種，用a=(0,1,…,J)表示。現實中，任一觀察對象某時刻的暴露變量A只能取J+1中的一種。但可以想象，假定第i個觀察對象某時刻的暴露值a為j(j=0,…,J)，將有一個與之相對應的結局變量，表示為Yi,a=j。變量Yi,a=j是一個虛擬變量，因為在現實中，觀察對象i的暴露變量A只會取0,1,…,J中的一個值，此時，現實中觀測不到暴露變量a分別取剩余J-1個值時結果變量Yi,a=j的取值。所有這些暴露值的結局變量Yi,a=j都被稱作反事實變量，研究對象i的反事實變量可以用向量(Yi,a=0,Yi,a=1,…,Yi,a=j)表示。通常，我們可以合理假定，當對象i接受暴露Xi=j，實際的結局事件Yi就等于Yi,a=j。對于每一個觀察對象，現實中都只能觀察到其反事實結局變量向量中一個元素，而向量中其他值在現實中則無法觀察到[2,3]。

在時點研究中，令a0=1和a0=0分別表示研究對象接受治療和未接受治療，用Ya0=1和Ya0=0分別表示研究對象接受治療和未接受治療的結局，即反事實變量。在臨床治療中，無法同時觀察到同一研究對象既接受治療又不接受治療的結局。如果某研究對象接受了治療，則此研究對象能觀察到的結局Y就為Ya0=1，而無法同時觀察到反事實變量向量中的另一取值Ya0=0[1]。

當暴露變量為時間依賴變量時，也可以定義相對應的反事實變量：將時間依賴暴露變量表示為At，t=0,1,…,K，并且在時間點t時At的取值at為0,1,…,J中的一個。對應于任何一個可能的暴露歷史，at，t=0,1,…,K，有一個相應的反事實變量，表示為Ya0,a1,…,ak。同樣，現實中，因為一個研究對象只有一種暴露歷史，因此無法觀察到其相應于其他可能的暴露歷史的反事實變量。如果t=1，對于試驗對象i，當Ai,0=1,Ai,1=0時，那么從時刻0至時刻1期間觀察到的結局變量為Yi,a0=1,a1=0，而反事實變量向量的其余部分，如Yi,a0=1,a1=1、Yi,a0=0,a1=1、Yi,a0=0,a1=0則觀察不到。下文主要基于時點研究，即t=0時闡述邊際結構模型的基本原理及其應用。

2.觀察性研究的因果推斷圖

為幫助理解邊際結構模型的基本原理，先簡要介紹觀察性研究中的因果關聯推斷圖。令：A表示暴露變量，L表示已觀測協變量，U表示未觀測協變量，Y表示結局變量。圖1展示了基線期上述各變量間的3種不同關系。

與圖1a相比，圖1b中基線期0中缺少由U指向A的箭頭，說明A不受U的影響，但受L的影響。在圖1c中，既缺少由U指向A的箭頭，也缺少由L指向A的箭頭，說明A不受混雜因素(即U和L)的影響。僅在圖1c的情況下，混雜因素才不影響A與Y之間的因果關聯推斷，但此類情況在現實中，除非是隨機試驗，否則基本不存在。

事實上，在觀察性研究中，由于無法得到未觀測協變量U的數據，當存在U0的影響時(圖1a)，如果沒有合理的前提假設，任何方法都不能得到A0與Y之間因果關聯的無偏估計。故下文介紹的邊際結構模型僅針對于只存在已觀測協變量L影響的情況(圖1b)。

圖1 觀察性研究中因果關聯推斷簡單示意圖(下標0表示基線期，此圖即時點研究的因果關聯示意圖。資料來源：文獻[1])

3.邊際結構模型的基本數學模型

在時點研究中，可采用線性模型、對數線性模型和logistic回歸模型估計暴露變量與結局變量間的因果關聯強度，并可獲得各研究對象的兩個反事實變量P[Ya0=1=1]和P[Ya0=0=1]的估計值，具體模型為：

P[Ya0=1]=ψ0+ψ1a0

(1)

log(P[Ya0=1])=θ0+θ1a0

(2)

logit(P[Ya0=1])=β0+β1a0

(3)

在上述模型中，若研究對象接受治療(a0=1)時，Ya0=Ya0=1；反之，若研究對象未接受治療(a0=0)時，Ya0=Ya0=0。模型(1)、(2)和(3)中因果關聯強度指標分別為：率差RD=ψ1、相對危險度RR=eθ1、優勢比OR=eβ1。

上述模型之所以被稱為邊際結構模型是因為：模型(1)、(2)和(3)均是建立在Ya0 = 1與Ya0=0的邊際分布基礎之上，而非二者的聯合分布基礎之上，故將此類模型稱為邊際模型。在經濟學和社會科學研究領域中常將基于反事實變量的概率模型稱為結構模型，故將此類模型稱為邊際結構模型。

在模型(1)、(2)和(3)中，每個模型都包含兩個未知參數且只有兩個未知概率P[Ya0=1=1]和P[Ya0=0=1]，因此是飽和模型。

對于存在已觀測協變量L影響的情況，可類似構建飽和模型(4)、(5)和(6)獲得未調整混雜因素的因果關聯強度：

(4)

(5)

(6)

通常，模型(1)、(2)和(3)的參數估計值與模型(4)、(5)和(6)的估計值不相同。僅在暴露變量不受混雜因素影響的情況下，模型(4)、(5)和(6)中的粗估計值才等價于真實的因果關聯估計值。

4.邊際結構模型的基本原理

在觀察性研究中，暴露變量可能受已觀測協變量(混雜因素)的影響。邊際結構模型正是通過構造已觀測協變量與暴露變量的模型，提出校正已觀測混雜偏倚的暴露變量逆概率權重，消除已觀測協變量的影響，從而獲得暴露變量與結局變量間真實的因果關聯。

下文以時點研究為例，重點介紹邊際結構模型逆概率權重的計算原理。

(1)時點研究中逆概率權重的估算

在時點研究中，研究對象i暴露的逆概率權重的計算公式為：

wi=1/P[A0=a0i|L0=l0i]

其中：loi表示研究對象i的已觀測協變量L0的觀測值；a0i為研究對象i的暴露變量A0的觀測值；P[A0=a0i|L0=l0i]表示研究對象i在L0=l0i的情況下，出現A0=a0i的概率。

可通過擬合A0與L0的logistic回歸模型獲得觀察對象接受暴露的概率P[A0=1|L0=l0]。估計研究對象接受暴露概率的模型為：

logit(P[A0=1|L0=l0])=α0+α1l0

(7)

可通過常用統計軟件得到logistic回歸模型中參數α0和α1的估計值。若某研究對象實際接受的暴露水平為A0=0，則其接受自身暴露水平的概率為1-P[A0=1|L0=l0]，進而得到各研究對象的逆概率權重wi。

(2)暴露逆概率權重的意義

在無未觀測混雜因素影響的假定下，可通過對各研究對象賦予逆概率權重w，來消除已觀測協變量L0的混雜影響。例如：對于某研究對象i，若其逆概率權重wi=4，則表示在虛擬人群中該研究對象將會被重復4次。

采用暴露逆概率權重，可在實際觀測人群的基礎上構造一個虛擬人群，此虛擬人群具有兩個重要性質：第一、與在實際人群中不同，在虛擬人群中A0不受L0的影響；第二、在虛擬人群中，P(Ya0=1=1)和P(Ya0=0=1)的概率與在實際人群中相同。因此，可通過逆概率權重控制L0造成的混雜偏倚，得到因果關聯的無偏估計。這正是邊際結構模型的精髓之所在。

5.邊際結構模型與傾向評分法的異同

傾向評分法(propensity score method)由Rosebaum和Rubin于1983年首次提出，是指在一定協變量條件下，一個觀察對象可能接受某種暴露的可能性。在非隨機研究中，處理組與對照組某些背景特征分布不同，每個個體是否接受“暴露”的概率受其他因素的影響[4]。在傾向評分法中，當一組觀測協變量L影響研究對象接受某感興趣的暴露時(A=1)，在該組協變量的影響下研究對象接受感興趣的暴露的概率為P(A=1|L=li)，此概率即為研究對象i的傾向評分pi。

傾向評分是反映所觀測協變量在兩組間均衡性的一個近似函數，它最大限度地概括了協變量的作用。若分別來自暴露組和對照組的兩個研究對象具有相同的傾向評分，可認為他們是被隨機分配到的兩組，接受暴露或者對照的概率相同[4-6]。傾向評分法的假定條件為：計算傾向評分的協變量包括所有影響分組的混雜因素，即不存在未觀測的混雜因素[4-5]。

邊際結構模型與傾向評分法均是控制混雜偏倚的有力工具，其應用都建立在不存在未觀測混雜因素影響的基礎之上[1,4]。

與邊際結構模型相比，盡管對于接受感興趣的暴露(A=1)的研究對象而言，1/pi就是邊際結構模型中的暴露逆概率權重；但對于未接受感興趣的暴露(A=0)的研究對象而言，其逆概率權重卻不是1/pi，而是1/(1-pi)[1,4]。Robins等人認為，傾向評分法可能受到匹配不充分或層內混雜控制欠佳的影響，特別是在樣本量較小時，在控制觀測混雜偏倚方面不如邊際結構模型[1,6]。

6.應用示例

美國密歇根大學心血管聯盟開展的一個多中心合作的觀察研究，其主要目的是研究因ST段抬高心肌梗死的住院患者在接受皮冠狀動脈介入治療前使用β受體阻滯劑對治療結局的影響。其中，術前使用β受體阻滯劑為暴露變量，術中死亡率為結果變量，吸煙狀態、高血壓、心臟衰竭等因素為協變量，變量賦值見表1。

考慮到本研究為觀察研究，且數據初步分析顯示協變量明顯影響到了暴露變量的取值，故考慮選用邊際結構模型來控制協變量的影響。

本例采用邊際結構模型，具體步驟如下：

(1)估算研究對象術前使用β受體阻滯劑的概率。參考公式(7)，采用logistics 回歸模型擬合邊際結構模型，以表1所列的協變量作為自變量，以術前使用β受體阻滯劑與否為因變量，模型估計參數見表2。依據此logistic回歸模型估計各研究對象術前使用β受體阻滯劑的概率。

(2)求逆概率權重。按公式wi=1/P[A0=a0i|L0=l0i]，可求出研究對象的逆概率權重。

(3)加權調整。采用第(2)步算得的權重進行加權，采用廣義線性模型中的重復測量logit模型分析數據。

結果顯示，采用邊際結構模型得到的暴露變量的效應值為：OR=0.3798，95%CI：0.2112-0.6828，P=0.0012。而不控制協變量，直接擬合暴露變量和結果變量的logistic 模型得到的暴露變量粗效應值為：cOR=0.4320，95%CI：0.1916-0.9739，P=0.0430。此結果表明，β受體阻滯劑能降低病人術中死亡率，但是，未調整協變量之前，β受體阻滯劑降低術中死亡率的效果被低估了。

6.復雜數據的邊際結構模型

當暴露變量為多水平變量(連續或等級)時，當觀察性研究中的暴露變量為時間依賴變量時，同樣可以基于保守的劑量反應關系構建線性logistic邊際結構模型，通過計算平穩權重來控制已觀測協變量的混雜影響。當存在刪失數據時，把刪失也看作是隨時間變化的處理，同樣可以采用邊際結構模型來調整失訪所致的偏倚[1-3]。由于本文篇幅有限，在此不作介紹。

討 論

邊際結構模型是近些年新出現的一類因果推斷方法。該法與傳統因果推斷方法有較大區別，它通過構造已觀測協變量與暴露變量的模型，通過逆概率權重反映已觀測協變量對暴露變量和結局變量的影響，然后根據逆概率權重消除已觀測協變量的影響，從而獲得暴露變量與結局變量間真實的因果關聯[1-2]。

上述實例顯示，邊際結構模型能很好地控制已測協變量所致的混雜偏倚。Robins等人指出，邊際結構模型尤其適合于暴露變量和協變量均為時間依賴變量的觀察性研究。由于缺乏此類實例，本文未能介紹此類應用。

參 考 文 獻

1. Robins J,Hernán M,Brumback B.Marginal structural models and causal inference in epidemiology.Epidemiology，2000,11(5):550-560.

2.章玨.基于邊際均值的多項分布數據的因果效應估計.南京:東南大學,2006.

3.Robins J.Marginal Structural Models versus Structural Nested Models as Tools for Causal Inference.Statistical Models in Epidemiology:The Environment and Clinical Trials.New York:Springer-Verlag,1999:95-134.

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

5.李智文,任愛國.傾向評分法在調整混雜中的應用.詹思延 主編.流行病學進展.第十二卷 .北京:人民衛生出版社,2010:358-376.

6.王永吉,蔡宏偉,夏結來,等.傾向指數第一講傾向指數的基本概念和研究步驟.中華流行病學雜志,2010,31 (3):347-348.

7.Valle J,Zhang M,Dixon S,et al.Impact of pre-procedural beta blockade on inpatient mortality in patients undergoing primary percutaneous coronary intervention for ST elevation myocardial infarction.Am J Cardiol，2013，111(12):1714-1720.