基于有向無環圖和線性回歸理論定量分析常見混雜的影響*

荊 典 段重陽 譚 銘，2△ 陳平雁△

【提 要】 目的 結合有向無環圖(directed acyclic graphs，DAGs)與線性回歸模型，提出常見混雜的定量分析方法。方法 針對典型的兩種DAGs(情形1：X←C→Y；情形2：X←C1→M←C2→Y)，基于線性回歸理論給出偏倚大小的定量表達式，并探討各參數對因果效應估計的影響。結果 對于情形1的DAG，暴露X與結局Y的因果效應估計需通過關閉二者間后門路來控制混雜C，否則X與Y間因果估計會受到混雜偏倚的影響。理論推導結果顯示，若只改變C對X的效應則同時會增加X的方差，此效應的強弱對混雜偏倚的影響是非單調的，除非在改變C與X的效應的同時控制X的方差不變，而C與Y間的效應強弱對混雜偏倚的影響則是單調的。對于情形2，又稱M-DAGs無需做變量控制即可做X與Y的因果估計，但當錯誤地控制碰撞點M后會導致后門路徑打開而出現混雜，此時需進一步控制C1和(或)C2來關閉后門路徑。我們用回歸理論解釋了該結果，并且得到當C1與X相關性較高時，同時控制C1和M的方法結果會不穩定。結論 本研究基于有向無環圖，根據線性回歸理論推導出常見混雜的定量分析方法，該方法也適用于無法觀測的混雜，為多因素因果推斷提供了一種實用工具。

在多變量的因果推斷中，如何選擇納入回歸分析的協變量始終是一個研究的熱點問題[1-4]。有向無環圖(directed acyclic graphs，DAGs)應用甚廣，能夠基于預先假定的因果關系找到應校正的混雜因素集合[5-8]。然而，DAGs只能定性表達混雜的影響[9-11]。如何結合DAGs的定性方式，又能定量表達混雜的影響是一個值得研究的課題。因此，本研究將通過DAGs說明傳統的混雜偏倚和由于校正碰撞點所引起的選擇偏倚，并經由線性回歸理論推導討論校正偏倚及改變變量間關系對于偏倚大小的影響，以期對DAGs的結論用回歸理論做解釋說明，并同時給出DAGs無法給出的混雜定量影響的相關結論。

方 法

DAGs利用圖示方法描述待研究的處理或暴露(X)、結局(Y)以及協變量之間的因果關系[10，12]。“后門路徑(backdoor path)”是DAGs的一個重要概念，其定義為即使移除所有從暴露指向其子代變量的路徑，在暴露和結局間仍存在一條非因果效應的路徑[12-13]。如果X和Y之間除因果路徑還存在開放的后門路徑，在估計X和Y之間的因果效應時就會產生偏倚。由X和Y的共同原因導致開放的后門路徑所引起的偏倚稱為混雜偏倚(confounding bias)，而由兩個變量(其中之一為X或/與X相關的變量，另一為Y或/與Y相關的變量)的共同效應導致開放的后門路徑所引起的偏倚稱為選擇偏倚(selection bias)或碰撞點偏倚(collider bias)。需要注意的是，對于碰撞點施加干預會打開額外的后門路徑并引起偏倚，如碰撞點DAG(M-DAG)[8，12-14]。

圖1為探討混雜常用的DAGs，其中圖(a)、(b)、(c)表示暴露X與結局Y之間除了因果路徑X→Y外，還存在開放的后門路徑，即(a)：X←C→Y；(b)：X←C←U→Y；(c)：X←U→C→Y，因此X與Y之間存在混雜。若要控制混雜則需要關閉后門路徑，如采用協變量校正、分層、匹配等方法。圖(d)、(e)常用于說明并非所有變量都可以校正，有時校正變量反而會引起混雜。由于碰撞點的存在，如(d)：C1→M←C2；(e)：X→C←U，除了因果路徑X→Y外，并沒有開放的后門路徑，因此不存在混雜。但是，如果過度校正了碰撞點，則反而打開了后門路徑，引起選擇偏倚。上述5個DAGs中我們可以發現(b)和(c)為(a)的變形，(e)為(d)的變形，所以我們選擇具代表性的DAG(a)和DAG(d)，結合線性回歸模型建立控制混雜的定量分析方法[12，15-16]。

圖1 常見混雜的DAGs

1.DAG(a)

DAG(a)為最常見的一種混雜偏倚，即由于共同原因C的存在導致X和Y之間存在開放的后門路徑(X←C→Y)。例如，在研究體重(X)是否會影響收縮壓(Y)時，發現年齡(C)均會對二者造成影響，即使體重和收縮壓不存在真實的因果效應，仍能通過年齡(C)發現關聯[12]。如下我們用回歸理論解釋該關系。

表示X受C的影響，令X=γ0+γ1C+εx；表示Y受X和C的共同影響，令Y=λ0+λ1C+λ2X+εy，其中γ1、λ1和λ2分別為C對X、C對Y和X對Y的效應；εx和εy分別為X和Y的噪聲。

由線性回歸模型Y=α+βX+ε擬合X和Y的效應，可得到β的估計值如下：

(1)

其中，Cov(X，Y)是X和Y的協方差，Var(X)是X的方差，SD(X)和SD(Y)是X和Y的標準差，rxy是X和Y的相關系數，并且rxy=Cov(X，Y)/[SD(Y)SD(X)]。由上述方程可進一步得到C同X和Y的相關關系如下：

Cov(C，X)=γ1Var(C)

Cov(C，Y)=λ1Var(C)+λ2Cov(C，X)=

(λ1+λ2γ1)Var(C)

Cov(X，Y)=λ1Cov(C，X)+λ2Var(X)=

2λ1λ2Cov(C，X)+Var(εy)=

(2)

由式(1)和式(2)可得β的期望為：

(3)

rxy的期望為：

(4)

(5)

(6)

更特殊的，當X對Y的效應為0時，即λ2=0時，式(2)可簡化為：

Cov(C，X)=γ1Var(C)

Cov(C，Y)=λ1Var(C)

Cov(X，Y)=λ1γ1Var(C)

(7)

以及

(8)

2.DAG(d)

DAG(d)被稱之為M-DAG，其中存在一個同時受到暴露X和結局Y的父代(ancestor)C1和C2影響的碰撞點M，此時校正M將會打開原本閉合的后門路徑(X←C1→M←C2→Y)。例如，研究總膽固醇(X)對收縮壓(Y)的影響時，發現年齡(C1)會影響總膽固醇(X)，也可影響空腹血糖(M)。與此同時，心率(C2)也會影響收縮壓(Y)和空腹血糖(M)。當研究對空腹血糖進行校正時，就會打開原本閉合的后門路徑從而引起偏倚[12]。如下我們同樣用回歸理論解釋該關系。

令M=α0+α1C1+α2C2+εm，X=γ0+γ1C1+εx，Y=λ0+λ1C2+λ2X+εy，其中α1、α2、γ1、λ1、λ2分別為C1對M、C2對M、C1對X、C2對Y、X對Y的效應。εm、εx和εy分別為M、X和Y的噪聲。由線性回歸模型Y=β0+β1X+ε擬合未校正M時X和Y的效應，用模型Y=β0+β1X+β2M+ε擬合校正M的情況，以及用模型Y=β0+β1X+β2M+β3C+ε(C為C1或C2)擬合同時校正M和C1或C2時的情況。基于回歸理論可知估計值β為：

(9)

由公式(9)可得由模型Y=β0+β1X+β2M+ε得出校正M后β1的估計值：

(10)

對于同時校正M和C1或C2，可得：

(11)

式中C表示C1或C2。

令H=XTX，可得β1的估計值為：

(12)

及：

(13)

更特殊的，當X對Y的效應為0時，即λ2=0有：

Cov(C1，X)=γ1Var(C1)

Cov(C1，M)=α1Var(C1)

Cov(C1，Y)=0

Cov(C2，X)=0

Cov(C2，M)=α2Var(C2)

Cov(C2，Y)=λ1Var(C2)

Cov(X，Y)=0

Cov(X，M)=α1γ1Var(C1)

Cov(M，Y)=α2λ1Var(C2)

(14)

若校正C1，可知Cov(X，Y)=0且Cov(C1，Y)=0，則公式(13)可簡化為：

(15)

若校正C2，可知Cov(X，Y)=0且Cov(C2，X)=0，則公式(13)可簡化為：

(16)

結 果

基于前面理論推導結果，我們可討論各參數對偏倚大小的影響。為形象展示各參數間關系，在具體討論中我們同時基于一假定的參數繪制了關系曲線，結果見圖2。

1.DAG(a)

(1)不對C校正時X與Y效應估計偏倚大小

由公式(3)可知，不對混雜因素C進行校正將會恒定錯誤估計真實因果效應為

(17)

具體可見圖2(A)。

針對之前提到的實例，當研究體重(X)對收縮壓(Y)的因果效應時，若不校正年齡(C)這一混雜因素，則會錯誤估計真實因果效應。并且當體重(X)確定后，其方差Var(X)相對確定，此時偏倚大小與年齡(C)的方差Var(C)、年齡(C)對體重(X)的因果效應γ1，年齡(C)對收縮壓(Y)的因果效應λ1在其余二者固定不變的情況下成線性正相關關系。

進一步研究改變C對X和C對Y的相關關系所帶來的可能影響。本研究考慮兩種方式改變其相關，其一是直接增加C對X或C對Y的因果效應，如增大年齡(C)對體重(X)的影響或是其對收縮壓(Y)的影響；其二是通過對X或者Y施加噪聲(給定SDs)。

(2)C對X和C對Y效應大小對X與Y效應估計偏倚大小的影響

當僅增加C對Y效應λ1時，Cov(C，Y)、Cov(X，Y)、Var(Y)、E(rxy)和E(rcy)隨之增加，從公式(17)知，此時除了λ1變化其余各項均不變，因此偏倚(17)隨之線性增加。具體可見圖2(C)。

(3)X、Y噪聲大小對X與Y效應估計偏倚大小的影響

以特殊情形，即不考慮X和Y存在因果效應為例，如假設體重(X)不會影響收縮壓(Y)，則可由公式(7)、(17)知，當僅增加X噪聲標準差SD(εx)時，Var(εx)和Var(X)隨之增加，故偏倚(17)隨之下降。當僅增加Y噪聲標準差SD(εy)時，偏倚(17)不變。當X和Y因果效應不為零時，即假設體重(X)會影響收縮壓(Y)時，則基于公式(2)上述結論依然成立。具體可見圖2(D)、(E)。

(4)對變量間相關系數估計的影響

對于DAG(a)，首先探討暴露X對結局Y真實因果效應改變所帶來的可能影響。由公式(2)～(5)可知固定其他變量不變時，增大X對Y的因果效應λ2，如體重(X)對收縮壓(Y)影響λ2增大時，Var(X)、Var(C)、Var(εy)、Cov(C，X)和λ1不變，Cov(C，Y)、Cov(X，Y)和Var(Y)逐漸增加。年齡(C)和體重(X)的相關系數rcx的期望不隨λ2的增加而改變，而體重(X)和收縮壓(Y)相關系數rxy的期望及年齡(C)和收縮壓(Y)相關系數rcy的期望隨λ2的增加而增大。

當X和Y的因果效應為0時，如假設體重(X)不會影響收縮壓(Y)，當僅增加體重(X)噪聲標準差SD(εx)時，可由公式(7)、(8)知E(rxy)、E(rcx)會隨著SD(εx)的增加而降低，而E(rcy)不變。當僅增加收縮壓(Y)噪聲標準差SD(εy)時，Var(εy)和Var(Y)增加，故僅E(rxy)和E(rcy)隨著SD(εy)的增加而降低。當X和Y的因果效應不為0時，即體重(X)會影響收縮壓(Y)，僅增加體重(X)噪聲標準差SD(εx)時Var(Y)也會增加，此時E(rxy)隨之不再單調下降，E(rcy)隨之下降，其余結論同上述保持一致。具體可見圖2(A)～(E)。

表1 DAG(a)結果匯總

2.DAG(d)

對于DAG(d)，主要研究混雜因素的校正對X和Y因果效應估計的影響。此時若不考慮X和Y存在因果效應，則由公式(3)可知通過模型Y=β0+β1X+ε所得的β1估計值是無偏的，即不對任何變量進行校正時不存在偏倚。由公式(10)可知對M進行校正時，Cov(X，M)和Cov(M，Y)會導致偏倚，且偏倚的期望為

(18)

上述關系示例圖見圖2(F)。

圖2 改變各參數對變量間關系的影響

例如研究總膽固醇(X)對收縮壓(Y)的影響時納入空腹血糖(M)會存在偏倚。當總膽固醇(X)對收縮壓(Y)不存在影響時，由公式(15)可知納入所在年齡(C1)后，X和C1以及M和C1的協方差會抵消X和M的協方差，即同時校正M和C1所得的估計值無偏；由公式(16)可知納入心率(C2)后，Y和C2以及M和C2的協方差會抵消Y和M的協方差，即同時校正M和C2所得的估計值無偏。當X對Y的效應不為0時，即總膽固醇(X)對收縮壓(Y)存在影響時，由公式(13)知該結論同樣成立，即在存在選擇偏倚時要想得到無偏的效應估計可將所在年齡(C1)或心率(C2)納入模型分析。

由于X和C2的相關系數恒定為0，而X和C1的高度相關使公式(15)中H的行列式接近于0，同時校正M和C1相比于校正M和C2會更加不穩定。

實例分析

由于實際應用中我們無法獲得真實的效應，因此無法對偏倚進行準確評價，本實例分析我們基于真實數據背景產生模擬數據，對本文中理論推導結果進行驗證。模擬數據的參數設置基于東莞寮步社區高血壓人群調查體檢數據庫得到的匯總數據。針對DAG(a)和DAG(d)模擬研究體重(X)或總膽固醇(X)對血壓收縮壓(Y)的影響[17-18]。本研究不關注抽樣誤差影響因此設置較大樣本量，基于檢驗效能大于95%設置樣本量為800。

對于DAG(a)，年齡(C)服從N(60，13.7)的正態分布，體重(X)(N(64.5，12.2))和收縮壓(Y)(N(136.1，16.3))的噪聲均服從N(0，1)的正態分布，并分別由回歸方程X=86.4-0.36C+εx和Y=117+0.22C+0.1X+εy產生。

對于DAG(d)，年齡(C1)服從N(60，13.7)的正態分布，心率(C2)服從N(72，11)的正態分布，總膽固醇(X)(N(199.5，45.0))、空腹血糖(M)(N(5.35，1.65))和收縮壓(Y)(N(136.1，16.3))的噪聲均服從N(0，0.1)的正態分布，并分別由回歸方程X=187+0.22C1+εx、Y=126+0.083C2+0.03X+εy和M=4.1+0.013C1+0.03C2+εm產生。

由表2可知，DAG(a)中校正混雜因素年齡(C)后得到的體重(X)系數估計值結果無偏，而不校正則會低估0.587個單位(公式(17)計算結果為0.591)。而DAG(d)中，僅納入總膽固醇(X)以及同時校正空腹血糖(M)和年齡(C1)或心率(C2)所得的系數估計值均無偏，但同時校正C1所得的結果不穩定。若直接納入總膽固醇(X)和空腹血糖(M)則會低估0.151個單位(公式(18)計算結果為0.140)。實例分析情況均同本研究理論推導結果一致。

表2 回歸分析結果

討論與結論

本研究基于兩種常見的DAGs，通過線性回歸理論對混雜偏倚和因選擇偏倚導致的混雜進行了分析，結果與DAGs的圖示規則保持一致。此外，本研究給出了偏倚大小的定量表達式(見公式(17)，(18))，探討了不同方式改變混雜因素同處理或結局的效應時對偏倚大小的影響。

對于經典的混雜DAG，當我們考慮混雜與暴露的效應對偏倚影響的大小時，需同時考慮暴露的方差的控制，使其更符合研究者的設定。比如現實數據分析中，若要考慮未觀測的混雜導致的偏倚影響大小時，對于觀測到的暴露，我們應該假設其方差固定，然后改變假定的混雜與暴露的效應，以評價其可能的影響。如基于數據分析得出X與Y的回歸系數為2，此時可以討論若真實回歸系數為0時，存在影響多大的潛在混雜C，可使得估計出的回歸系數為2，即當存在影響多大的潛在混雜C時可扭轉當前回歸系數2到0或者負值?；诠?17)，此時X的方差已知，則假定真實的回歸系數為0時，λ1γ1Var(C)作為一個整體取值為λ1γ1Var(C)=2Var(X)可使得估計出的有偏的回歸系數為2，此時可將該值與觀測到的專業上具有意義的各混雜λ1γ1Var(C)相應值進行比較，來判斷該λ1γ1Var(C)=2Var(X)存在的可能性，如果該可能性較小，則認為當前估計結果受潛在混雜影響而扭轉結局的可能性較小[19]。

對于M-DAG，不對任何協變量做處理時估計值無偏。若對碰撞點進行校正，還需校正位于碰撞點和處理或結局間的變量以得到真實因果效應的無偏估計。在變量間相關性較強的情況下同時校正碰撞點和位于碰撞點和結局間的變量所得的估計結果會更穩定。

由于篇幅限制，本研究僅考慮變量服從正態分布的情況。雖然結局指標的類型不會影響DAGs的混雜理論，但在考慮實際應用的回歸分析策略時，結局的類型不能忽略，例如二分類結局利用logistic回歸對協變量進行校正可能同標準線性回歸大不相同。同時，變量間成非線性關系時混雜因素的影響也有待探討。

本研究建立了將DAGs與線性回歸模型結合估計混雜偏倚的定量分析方法，該方法還適用于無法觀測的混雜，為多因素因果推斷提供了一種實用工具。