基于工具變量識別因果效應(yīng)以及用數(shù)據(jù)區(qū)分不同模型

胡純嚴 ，胡良平 ，2*

（1.軍事科學(xué)院研究生院，北京 100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029*通信作者：胡良平，E-mail：lphu927@163.com）

因果圖理論的4個主要內(nèi)容包括關(guān)聯(lián)的來源、因果圖模型的統(tǒng)計性質(zhì)、識別和調(diào)整以及工具變量，本文在介紹這些理論的基礎(chǔ)上，針對兩個實例，并借助SAS軟件完成用工具變量識別因果效應(yīng)以及用數(shù)據(jù)區(qū)分不同模型的任務(wù)。

1 因果圖理論的4個主要內(nèi)容

1.1 關(guān)聯(lián)的來源

在因果圖過程中，每個因果圖模型都必須是有向無環(huán)圖（Directed acyclic graph，DAG）。由DAG表示的因果圖模型對信息在底層數(shù)據(jù)生成過程中的流動方式具有明確的定義。該信息流由三個圖形結(jié)構(gòu)封裝，可用于在DAG中組裝每條路徑［1］。這三個結(jié)構(gòu)對應(yīng)因果圖模型中關(guān)聯(lián)的三個基本來源（即因果關(guān)系、混淆和內(nèi)生選擇），這些結(jié)構(gòu)如圖1所示。

圖1 關(guān)聯(lián)的3個基本來源Figure 1 Three fundamental sources of association

在因果關(guān)系U→V→W中，變量U和W是關(guān)聯(lián)的，這種關(guān)聯(lián)是因果鏈的結(jié)果。如果要對中介變量V進行調(diào)節(jié)，那么這將阻塞信息流，從而使變量U和W不再關(guān)聯(lián)。

在混淆的關(guān)系U←V→W中，沒有與變量U和W相關(guān)的因果路徑，然而，U和W仍然是相關(guān)的，這種關(guān)聯(lián)是由混雜變量V引起的，它是變量U和W的共同父項，如果研究者要以共同原因V為條件，那么這將阻塞信息流，從而使變量U和W不再關(guān)聯(lián)。

在內(nèi)生選擇結(jié)構(gòu)U→V←W中，變量U和W共同決定其共同子變量V的值，但變量U和W不相關(guān)，然而，如果研究者要以共同的結(jié)果變量V為條件，就可以創(chuàng)建一個信息流，變量U和W就會關(guān)聯(lián)起來［2］。

粗略地說，如果在因果效應(yīng)分析中有一個處理變量（如U）和一個結(jié)果變量（如W），那么目標是消除變量U和W之間的非因果關(guān)聯(lián)，并保持因果關(guān)聯(lián)不變。因此，這三個基本圖形結(jié)構(gòu)不僅對應(yīng)關(guān)聯(lián)的三個基本來源，也對應(yīng)偏差的三個基本來源。一般來說，當研究者有一組處理變量和結(jié)果變量時，如果控制因果路徑上的一個變量，就會阻塞流經(jīng)該因果路徑的信息流，這被稱為過度控制偏差；同樣，如果研究者無法控制一個令人困惑的共同原因，那么處理變量和結(jié)果變量之間的一些關(guān)聯(lián)就是混淆的結(jié)果，被稱為混淆偏差；而如果研究者控制變量的共同結(jié)果，就會在處理變量和結(jié)果變量之間建立非因果關(guān)系的關(guān)聯(lián)，被稱為內(nèi)生選擇偏差。

1.2 因果圖模型的統(tǒng)計性質(zhì)

1.2.1 局部馬爾科夫性

有兩種方法可以解釋DAG中的假設(shè)：①DAG是“組織有關(guān)外部干預(yù)及其互動的一種結(jié)構(gòu)”［3］；②DAG定義了一組變量之間的信息流。這兩種解釋在另外兩種假設(shè)下是等效的［1］：①DAG中的變量滿足局部馬爾可夫性；②DAG滿足弱忠實性屬性。局部馬爾可夫性質(zhì)表明，DAG中的每個變量在統(tǒng)計上是獨立的，取決于它的父變量和非退化變量集［4］。

1.2.2 d分離

如果以下任一條件成立，則DAG中的路徑稱為由一組變量Z形成的d分離：①路徑包含鏈U→V→W或分叉U←V→W，使得V∈Z；②路徑包含一個碰撞器U→V←W，使得V?Z以及V的后代不在Z中。

1.2.3 阻塞/非阻塞

一條d分離的路徑被稱為阻塞。如果X中的一個節(jié)點和Y中的一個節(jié)點之間的每條路徑都被阻塞，則變量集X通過變量集Z與變量集Y進行d分離。阻塞/非阻塞術(shù)語反映了因果圖模型中的信息流，如果路徑被阻塞，則信息不會流經(jīng)該路徑；如果路徑未被阻塞，則信息可能會流經(jīng)該路徑。d分離和信息流之間的聯(lián)系體現(xiàn)在弱忠實性假設(shè)中。弱信度表示，如果一個DAG中的兩個變量X和Y不是d分離的，那么這兩個變量至少依賴于一個在DAG上分解的分布［5-6］。

1.2.4 全局馬爾科夫性

通過將因果圖模型解釋為代表變量之間關(guān)聯(lián)流的DAG，可以將DAG背后的因果假設(shè)轉(zhuǎn)化為條件獨立性。具體來說，如果兩個變量在DAG中由集合Z進行d分離，那么這兩個變量必須在統(tǒng)計上獨立于Z。換言之，d分離是一個全局馬爾可夫性質(zhì)。如果條件獨立性只包含觀察到的變量，則可以使用觀察到的數(shù)據(jù)執(zhí)行統(tǒng)計檢驗，以查看獨立性是否成立。因此，d分離標準確定了因果圖模型具有可觀測和可檢驗的含義［1］。

事實上，DAG的全局馬爾可夫性質(zhì)和局部馬爾可夫性質(zhì)在邏輯上是等價的［4］。如果研究者有一個局部或全局馬爾可夫性質(zhì)的完整列表，就可以使用semigraphoid公理［7-8］推導(dǎo)出另一個列表。在CAUSALGRAPH過程中，可以使用PROC CAUSALGRAPH語句中的IMAP選項來請求這些屬性的列表。

1.3 識別和調(diào)整

一對變量之間的統(tǒng)計關(guān)聯(lián)可以分為兩個部分：因果部分和非因果部分（虛假部分）。如果所有虛假關(guān)聯(lián)都能被消除，那么因果關(guān)系就會被識別出來。因此，一種可能的識別方法是調(diào)整識別，這是回歸和匹配因果效應(yīng)識別的基礎(chǔ)［2］。

當研究者使用調(diào)整識別時，會尋找一個調(diào)整集，當在分析中進行控制時，它會阻塞DAG中的所有非因果路徑，而不會阻塞任何因果路徑。路徑的因果屬性是從模型中邊的方向繼承的。也就是說，因果屬性是因果圖模型的屬性，在分析過程中不會改變。然而，路徑是否被阻塞不僅取決于代表因果圖模型的DAG結(jié)構(gòu)，還取決于調(diào)整集中包含的變量集。因此，研究者必須謹慎地選擇一個調(diào)整集，以便在不引入任何過度控制或內(nèi)生選擇偏差的情況下消除所有混淆偏差。

1.4 工具變量

在實際的醫(yī)學(xué)研究中，一般來說，變量可分為自變量、中介變量和因變量（結(jié)果變量）。然而，統(tǒng)計學(xué)上還提出了一種“工具變量”［9］，其定義如下：某個變量Z與模型中某個自變量X高度相關(guān)，但卻不與隨機誤差項相關(guān)，那么就可以用變量Z與模型中相應(yīng)回歸系數(shù)得到一個一致估計量，在模型的參數(shù)估計過程中，變量Z被作為一個工具使用，故稱為工具變量。為了加深對工具變量的理解，舉例如下［10］：研究者將抽煙的孕婦隨機分成兩組，試驗組接受“減少或停止抽煙”的勸告或鼓勵，而對照組未接受勸告或鼓勵。研究者記錄每位受試者兩個結(jié)果變量的取值，即孕婦在8個月的孕期內(nèi)每天抽煙支數(shù)（記為S）及其嬰兒的出生體重（記為B）。研究者關(guān)心的是S對B的因果效應(yīng)（雖然S與B之間的關(guān)聯(lián)可能存在混淆，如對S可能存在測量誤差），為了解決S對B的因果效應(yīng)的估計問題，研究者利用隨機化方法（R），假設(shè)R可能與S高度相關(guān)，但R僅通過對S的影響進而對B產(chǎn)生影響（即以S為條件，隨機化對B沒有影響）。在這種情況下，變量R被稱為工具變量。

2 因果圖過程的應(yīng)用

2.1 使用工具變量識別因果效應(yīng)

2.1.1 實例與背景信息

【例1】沿用文獻［9］中的“Example 34.2”，此例的結(jié)論為：“不可能使用調(diào)整集來確定持久性全氟烷基物質(zhì)（PFAS）對母乳喂養(yǎng)持續(xù)時間（Duration）的因果效應(yīng)”。試問：原因是什么？如何解決所面臨的問題？

【分析與解答】此例采用的是Timmermann等［11］的因果圖模型，研究了法羅群島居民中母親接觸PFAS與Duration之間的關(guān)系。該例表明，研究者無法構(gòu)建調(diào)整集來估計處理變量PFAS對結(jié)果變量Duration的因果效應(yīng)。這是因為處理變量和結(jié)果變量之間存在混淆偏差，這些混淆來自未觀察到的變量“行為被認為是健康的程度”（HealthBehavior）和“先前是否有母乳喂養(yǎng)經(jīng)驗”（PrevBF）。

在許多存在未測量混淆的情況下，如果愿意假設(shè)因果圖模型中的某些邊具有特定的參數(shù)形式，通過使用工具變量［12-13］，仍然可以估計因果效應(yīng)。

2.1.2 用SAS實現(xiàn)因果效應(yīng)分析

2.1.2.1 初步分析

以下語句調(diào)用PROC CAUSALGRAPH過程列出可用于估計因果效應(yīng)的工具變量。設(shè)所需要的SAS程序如下；

【SAS輸出結(jié)果及解釋】

輸出結(jié)果為兩類變量：第一類為“工具變量”，即孩子出生時母親的年齡（Age）；第二類為“條件變量”，包括飲酒（Alcohol）、母親孕前體重指數(shù)（BMI）、母親接受初等教育的時間（Education）、胎次（Parity）和吸煙（Smoking）。

輸出結(jié)果表明：變量Age可用于確定處理變量PFAS對結(jié)果變量Duration的因果效應(yīng)。

2.1.2.2 改變條件變量再分析

構(gòu)造工具變量時產(chǎn)生的條件集可能不是最小的。例如，以下檢驗表明，如果僅調(diào)整Education和Parity兩個變量，也可以使用變量Age作為工具變量。

所需要的SAS程序與“2.1.2.1節(jié)”的SAS程序相同，只需要在UNMEASURED語句之后增加以下語句：

testid "Minimal CIV" Age/conditional=（Education Parity）；

【SAS程序說明】這里的新功能是使用TESTID語句。研究者想調(diào)查當條件變量為Education和Parity時，Age是否可以作為因果效應(yīng)分析的工具變量。這組條件變量是之前分析中提出的條件變量的適當子集。

【SAS輸出結(jié)果及解釋】

當條件變量為Education和Parity時，以Age作為工具變量，研究處理變量PFAS對結(jié)果變量Duration的因果效應(yīng)是有效的。

2.2 用數(shù)據(jù)區(qū)分模型

2.2.1 實例與背景信息

【例2】沿用文獻［9］中的“Example 34.3”，研究者設(shè)定了兩個模型，見圖2、圖3。

圖2 血清尿酸鹽對心血管疾病風險影響的第一個可能因果圖模型Figure 2 The first possible causal models of the effect of serum urate on risk of cardiovascular disease

圖3 血清尿酸鹽對心血管疾病風險影響的第二個可能因果圖模型Figure 3 The second possible causal models of the effect of serum urate on risk of cardiovascular disease

在圖2中，假設(shè)血壓（PreviousBP）和使用抗高血壓藥物（AntiHypertensiveUse）介導(dǎo)了變量尿酸鹽（Urate）對心血管疾?。–VD）的影響。在圖3中，因果方向相反，并且假設(shè)使用AntiHypertensiveUse對Urate直接產(chǎn)生因果效應(yīng)（對CVD而言，Urate就成了中介變量）。

兩個模型之間的差異以灰色（特指變量Anti-HypertensiveUse到Urate之間的灰色箭頭）突出顯示。兩個模型中的處理變量Urate和結(jié)果變量CVD均位于陰影框內(nèi)。注意，變量營養(yǎng)（Nutrition）對應(yīng)于潛在結(jié)構(gòu)（在圖中以虛線框表示），故不進行測量或觀察。還假設(shè)未測量變量先前血壓（PreviousBP）。當研究者有多個可能的因果圖模型時，若能找到一個對所有模型都有效的公共調(diào)整集，便可使用調(diào)整技術(shù)來估計數(shù)據(jù)的因果效應(yīng)。

試基于圖2和圖3中設(shè)定的因果圖模型，構(gòu)建它們的公共調(diào)整集。

2.2.2 用SAS實現(xiàn)因果效應(yīng)分析

以下語句調(diào)用PROC CAUSALGRAPH來構(gòu)造公共調(diào)整集。設(shè)所需要的SAS程序如下：

【SAS輸出結(jié)果及解釋】

NOTE：沒有滿足指定準則的適用于所有模型的調(diào)整設(shè)置。

輸出結(jié)果中的注釋表明，在本例中，未能找到共用調(diào)整集。因此，研究者必須找到單獨的調(diào)整集，然后使用每個模型分別估計因果效應(yīng)，或者必須確定最能代表數(shù)據(jù)生成過程的模型。PROC CAUSALGRAPH語句中的IMAP（條件獨立性假設(shè)）選項可對模型屬性進行分析。

2.3 枚舉和檢驗隱含的統(tǒng)計特性

2.3.1 隱含的統(tǒng)計特性

研究者可以根據(jù)可用數(shù)據(jù)檢驗因果圖模型中隱含的統(tǒng)計特性。如果隱含的統(tǒng)計特性在數(shù)據(jù)中不存在，則應(yīng)考慮修改或放棄模型。如果有多個模型，則可以比較這些模型的統(tǒng)計含義，以找到一個模型中具有的隱含統(tǒng)計特性，而不是其他模型中隱含的統(tǒng)計特性。然后可以在數(shù)據(jù)中檢驗此特性，并使用相應(yīng)的檢驗結(jié)果來確定哪個模型最能代表真實的數(shù)據(jù)生成過程。

下面的語句調(diào)用PROC CAUSALGRAPH過程來枚舉本例中兩個模型具有的統(tǒng)計特性。對于每個模型，過程生成的條件獨立屬性表由ODS OUTPUT語句保存到一個數(shù)據(jù)集。與“用數(shù)據(jù)區(qū)分模型”的SAS程序基本相同，下面僅列出不同之處：

proc causalgraph imap=global；

ods output imap=SimpleBPIndep；

proc causalgraph imap=global；

ods output imap=AltBPIndep；

【SAS程序說明】PROC CAUSALGRAPH過程中的IMAP=GLOBAL選項都會為使用MODEL語句指定的每個模型生成一個全局馬爾可夫?qū)傩员?。每個全局馬爾可夫?qū)傩杂蓛蓚€變量組成，這兩個變量在統(tǒng)計上獨立于另一個變量集（可能為空，在這種情況下，獨立性是無條件的）。如果觀察到馬爾可夫?qū)傩灾械拿總€變量，則可以使用數(shù)據(jù)執(zhí)行統(tǒng)計檢驗（例如可以檢驗零偏相關(guān)），以查看該屬性是否可以被修改。而涉及一個或多個未測量變量的獨立性屬性無法檢驗。

【SAS輸出結(jié)果及解釋】

與圖2對應(yīng)的輸出結(jié)果顯示，可以找到僅包含變量肥胖（Obesity）的調(diào)整集，而且，它是最小的調(diào)整集?；谧兞縊besity，可有效地進行Urate對CVD的因果效應(yīng)分析。

與圖3對應(yīng)的輸出結(jié)果顯示，可以找到11個調(diào)整集，因篇幅所限，具體輸出結(jié)果從略。

2.3.2 輸出數(shù)據(jù)集中的觀測

以下程序輸出每個模型的前10個觀察到的條件獨立性屬性。為了簡潔起見，本例重點介紹前10個結(jié)果。設(shè)所需要的SAS程序如下：

【SAS輸出結(jié)果及解釋】

因篇幅所限，與模型1（見圖2）和模型2（見圖3）對應(yīng)的“條件獨立性屬性”的輸出結(jié)果從略，現(xiàn)概要解釋如下。

對這兩個模型進行比較，兩個模型有4個條件集，使得變量AntiHypertensiveUse的使用在條件上獨立于變量CVD，并且這4個條件集對于兩個模型都是相同的。因此，研究者無法通過檢驗變量Anti-HypertensiveUse和CVD之間的條件獨立屬性來區(qū)分這兩個模型。接下來，比較變量肌酐（Creatinine）和CVD的條件獨立性屬性。Thor12SimpleBP模型（模型1）有4個這樣的屬性，但Thor12AltBP模型（模型2）有5個這樣的屬性。Thor12AltBP模型給出了Creatinine和CVD在集合上獨立的統(tǒng)計含義（Anti-HypertensiveUse、CurrentBP 及 Obesity），但 Thor12 SimpleBP模型中未給出此含義。

如果研究者要找到變量Creatinine和CVD之間的非零偏相關(guān)（在排除AntiHypertensiveUse、CurrentBP及Obesity后），研究者將有證據(jù)拒絕Thor12AltBP模型。研究者可以繼續(xù)對兩個模型中唯一的獨立屬性進行類似分析。最好的模型是其條件獨立性屬性與可用數(shù)據(jù)中的零偏相關(guān)最為匹配的模型。

3 討論與小結(jié)

3.1 討論

工具變量在因果效應(yīng)分析中起著重要作用，如果在某個實際問題中確實存在一個或多個工具變量，它們必將對其他變量起混雜效應(yīng)。因此，找出全部工具變量，并在其后的統(tǒng)計分析中充分發(fā)揮其作用（例如工具變量回歸分析［9-10］，它有別于通常的回歸分析［14-15］），將有助于獲得正確的統(tǒng)計結(jié)果和結(jié)論。

3.2 小結(jié)

本文介紹了因果圖理論的4個主要內(nèi)容，包括關(guān)聯(lián)的來源、因果圖模型的統(tǒng)計性質(zhì)、識別和調(diào)整以及工具變量；針對兩個實例并基于SAS軟件［9］，完成了使用工具變量識別因果效應(yīng)以及用數(shù)據(jù)區(qū)分不同模型的因果效應(yīng)分析任務(wù)。