荊 典 段重陽(yáng) 譚 銘，2△ 陳平雁△

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

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

方 法

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

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

圖1 常見(jiàn)混雜的DAGs

1.DAG(a)

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

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

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

(1)

其中，Cov(X，Y)是X和Y的協(xié)方差，Var(X)是X的方差，SD(X)和SD(Y)是X和Y的標(biāo)準(zhǔn)差，rxy是X和Y的相關(guān)系數(shù)，并且rxy=Cov(X，Y)/[SD(Y)SD(X)]。由上述方程可進(jìn)一步得到C同X和Y的相關(guān)關(guān)系如下：

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)

更特殊的，當(dāng)X對(duì)Y的效應(yīng)為0時(shí)，即λ2=0時(shí)，式(2)可簡(jiǎn)化為：

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，其中存在一個(gè)同時(shí)受到暴露X和結(jié)局Y的父代(ancestor)C1和C2影響的碰撞點(diǎn)M，此時(shí)校正M將會(huì)打開(kāi)原本閉合的后門路徑(X←C1→M←C2→Y)。例如，研究總膽固醇(X)對(duì)收縮壓(Y)的影響時(shí)，發(fā)現(xiàn)年齡(C1)會(huì)影響總膽固醇(X)，也可影響空腹血糖(M)。與此同時(shí)，心率(C2)也會(huì)影響收縮壓(Y)和空腹血糖(M)。當(dāng)研究對(duì)空腹血糖進(jìn)行校正時(shí)，就會(huì)打開(kāi)原本閉合的后門路徑從而引起偏倚[12]。如下我們同樣用回歸理論解釋該關(guān)系。

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

(9)

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

(10)

對(duì)于同時(shí)校正M和C1或C2，可得：

(11)

式中C表示C1或C2。

令H=XTX，可得β1的估計(jì)值為：

(12)

及：

(13)

更特殊的，當(dāng)X對(duì)Y的效應(yīng)為0時(shí)，即λ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)可簡(jiǎn)化為：

(15)

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

(16)

結(jié) 果

基于前面理論推導(dǎo)結(jié)果，我們可討論各參數(shù)對(duì)偏倚大小的影響。為形象展示各參數(shù)間關(guān)系，在具體討論中我們同時(shí)基于一假定的參數(shù)繪制了關(guān)系曲線，結(jié)果見(jiàn)圖2。

1.DAG(a)

(1)不對(duì)C校正時(shí)X與Y效應(yīng)估計(jì)偏倚大小

由公式(3)可知，不對(duì)混雜因素C進(jìn)行校正將會(huì)恒定錯(cuò)誤估計(jì)真實(shí)因果效應(yīng)為

(17)

具體可見(jiàn)圖2(A)。

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

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

(2)C對(duì)X和C對(duì)Y效應(yīng)大小對(duì)X與Y效應(yīng)估計(jì)偏倚大小的影響

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

(3)X、Y噪聲大小對(duì)X與Y效應(yīng)估計(jì)偏倚大小的影響

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

(4)對(duì)變量間相關(guān)系數(shù)估計(jì)的影響

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

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

表1 DAG(a)結(jié)果匯總

2.DAG(d)

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

(18)

上述關(guān)系示例圖見(jiàn)圖2(F)。

圖2 改變各參數(shù)對(duì)變量間關(guān)系的影響

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

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

實(shí)例分析

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

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

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

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

表2 回歸分析結(jié)果

討論與結(jié)論

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

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

對(duì)于M-DAG，不對(duì)任何協(xié)變量做處理時(shí)估計(jì)值無(wú)偏。若對(duì)碰撞點(diǎn)進(jìn)行校正，還需校正位于碰撞點(diǎn)和處理或結(jié)局間的變量以得到真實(shí)因果效應(yīng)的無(wú)偏估計(jì)。在變量間相關(guān)性較強(qiáng)的情況下同時(shí)校正碰撞點(diǎn)和位于碰撞點(diǎn)和結(jié)局間的變量所得的估計(jì)結(jié)果會(huì)更穩(wěn)定。

由于篇幅限制，本研究?jī)H考慮變量服從正態(tài)分布的情況。雖然結(jié)局指標(biāo)的類型不會(huì)影響DAGs的混雜理論，但在考慮實(shí)際應(yīng)用的回歸分析策略時(shí)，結(jié)局的類型不能忽略，例如二分類結(jié)局利用logistic回歸對(duì)協(xié)變量進(jìn)行校正可能同標(biāo)準(zhǔn)線性回歸大不相同。同時(shí)，變量間成非線性關(guān)系時(shí)混雜因素的影響也有待探討。

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