基于因果圖模型應用調(diào)整集估計數(shù)據(jù)的因果效應

胡純嚴 ，胡良平 ，2*

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

SAS/STAT中的PROC CAUSALGRAPH過程為因果圖過程［1］，該過程可以在不提供數(shù)據(jù)的前提條件下，基于設(shè)定的因果圖模型進行統(tǒng)計推斷。該過程有5個局限性：①因果圖過程不能處理有向循環(huán)的因果圖模型；②因果圖過程不能評估動態(tài)處理方案；③因果效應識別是一個總體概念；④因果效應識別是一個非參數(shù)概念；⑤因果圖過程不能識別某些因果圖模型中的因果效應。本文在介紹因果圖過程的局限性之后，針對一個實例并借助SAS軟件，實現(xiàn)基于因果圖模型應用調(diào)整集估計數(shù)據(jù)的因果效應。

1 因果圖過程的局限性

1.1 不能處理有向循環(huán)的因果圖模型

因果圖過程分析代表因果圖模型的有向無環(huán)圖（Directed acyclic graphs，DAG），這些DAG不能包含有向循環(huán)。在兩個變量（直接或間接）相互導致的情況下，基于DAG的因果圖過程分析可能存在困難。對于這種情況，一種常見的方法是引入額外的變量，以便在更精確的時間尺度上描述數(shù)據(jù)生成過程［2-3］。

1.2 不能評估動態(tài)處理方案

因果圖過程使研究者能夠在識別分析中指定多個處理變量和結(jié)果變量。當指定多個處理變量時，因果效應被解釋為聯(lián)合因果效應。也就是說，因果效應被解釋為同時對所有處理變量施加特定值的假設(shè)結(jié)果，研究者還可以將多個處理變量解釋為順序處理行動，前提是處理順序是預先確定的［2］。然而，研究者不能使用因果圖過程來評估動態(tài)處理方案的可識別性。

當研究者指定多個結(jié)果變量時，每個結(jié)果都被單獨解釋為一個獨特的因果效應。雖然解釋是獨立的，但因果圖過程只構(gòu)建對每個結(jié)果變量有效的調(diào)整集。在某些情況下，可能不存在此類的調(diào)整集，即使可以分別確定對每個結(jié)果的因果效應。例如，如果X對Y1的因果效應只能通過調(diào)整集Z1識別，而X對Y2的因果效應只能通過調(diào)整集Z2識別，其中，Z1和Z2是兩個不相交集，則不存在同時對兩個結(jié)果變量有效的調(diào)整集。

1.3 因果效應識別是一個總體概念

根據(jù)觀測數(shù)據(jù)估計的因果效應沒有有效的因果解釋，除非這些數(shù)據(jù)以因果圖模型的形式得到一組因果假設(shè)的補充［4］。然而，因果圖模型代表了在總體水平上變量之間的假設(shè)關(guān)系，而不是在個體水平上。因此，使用DAG描述因果效應識別的理論不考慮取樣變異性，識別條件在漸近極限下有效（隨著觀察次數(shù)的增加）［2］。成功的識別策略（使用調(diào)整集或條件工具變量）是使用非隨機試驗數(shù)據(jù)估計因果效應的第一步［5］。研究者應仔細考慮取樣變異在估計因果效應和檢驗模型的顯著性時的作用。

1.4 因果效應識別是一個非參數(shù)概念

因果效應的可識別性是一個完全非參數(shù)的概念，因為它不依賴于因果模型中變量和邊的分布或函數(shù)形式。然而，識別策略以及由該策略計算的任何估計值應被理解為以假設(shè)因果模型的有效性為條件［2］。此外，當因果效應被證明是確定的（例如使用調(diào)整集），這并不意味著研究者可以自由選擇一個參數(shù)估計器來量化效應，參數(shù)估計的適用性取決于參數(shù)假設(shè)，這些假設(shè)與因果圖模型的假設(shè)是分開的，必須針對每個具體情況進行證明［6］。

1.5 不能識別某些因果圖模型中的因果效應

在實踐中，常出現(xiàn)不能識別某些因果圖模型中的因果效應的情況。當在特定的因果圖模型中無法確定因果效應時，可采取一些補救措施：①研究者可以修改因果圖模型的假設(shè)，以查看數(shù)據(jù)生成過程是否可以由替代模型進行描述；②研究者可以考慮觀測其他變量，可能采取的形式是為以前未測量的變量添加觀測值，或為現(xiàn)有模型添加新變量和邊［4］，然而，在現(xiàn)有的一組變量中添加邊對識別因果關(guān)系不僅沒有幫助，甚至可能有害［4-5］。

2 應用調(diào)整集估計數(shù)據(jù)的因果效應

2.1 問題與背景信息

【例1】沿用文獻［1］中“Example 34.3”的問題和背景信息，模型中對處理變量尿酸鹽（Urate）和結(jié)果變量心血管疾病（CVD）進行了陰影處理。假設(shè)變量營養(yǎng)（Nutrition）對應于潛在結(jié)構(gòu)，故不進行測量。還假設(shè)變量先前高密度脂蛋白（PreviousHDL）為未測量變量。研究者設(shè)定變量之間的關(guān)系如圖1所示［6］。試使用因果圖過程來估計具備有效因果解釋的因果效應的大小。

圖1 血清尿酸鹽對心血管疾病風險影響的因果圖模型Figure 1 Causal graph model of the effect of serum urate on risk of cardiovascular disease

2.2 分析因果圖模型的思路

2.2.1 基本方法

要從數(shù)據(jù)集估計具備有效因果解釋的因果效應，可使用以下方法：①仔細考慮數(shù)據(jù)生成過程，并創(chuàng)建一個因果假設(shè)列表，以準確表示該過程，在因果圖模型中對這些假設(shè)進行編碼；②使用此圖形模型查找有效的識別策略；③利用識別結(jié)果構(gòu)造一個估計量，如分層估計量。

在大多數(shù)實際情況下，真正的數(shù)據(jù)生成過程并不明確。研究者必須定義一個假設(shè)，并用因果圖模型來呈現(xiàn)。要構(gòu)建這個因果圖模型，研究者可以依據(jù)專家意見、已建立的科學理論、先前的經(jīng)驗或其他可靠的知識來源。

2.2.2 產(chǎn)生模擬數(shù)據(jù)集

以下數(shù)據(jù)步創(chuàng)建了一個與圖1中的模型一致的模擬數(shù)據(jù)集，并定義了真正的數(shù)據(jù)生成過程。設(shè)所需要的SAS程序如下：

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

模擬數(shù)據(jù)集的前10行見表1。

表1 模擬數(shù)據(jù)集的前10行Table 1 The first 10 lines of the simulated data set

2.2.3 輸出Urate的匯總統(tǒng)計量

使用模擬數(shù)據(jù)集創(chuàng)建Urate的匯總統(tǒng)計量。設(shè)所需要的SAS程序如下：

proc means data=CVDdata；

var Urate；

ods output Summary=SampleMeansOutput；

run；

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

匯總統(tǒng)計量如表2所示。研究者可以使用ODS OUTPUT語句將變量Urate的匯總統(tǒng)計量存儲在輸出數(shù)據(jù)集中。在后面的分析中，研究者將使用此信息來定義感興趣的因果效應的處理和對照水平。

表2 Urate的匯總統(tǒng)計量Table 2 Summary statistics for urate

2.3 找出模型中可能的調(diào)整集

2.3.1 基本情況

在此例中，處理或暴露變量Urate是連續(xù)的。此外，該變量對中介變量痛風（Gout）和結(jié)果變量CVD的影響是非線性的。因為Urate沒有天然的處理和對照水平，所以研究者必須以某種方式定義感興趣的因果關(guān)系。常見的因果效應度量是平均處理效應或預期風險差異，即明確定義的處理和對照條件或水平之間的預期潛在結(jié)果值的差異。

在本例中，研究者認為感興趣的因果關(guān)系是CVD的預期風險差異，該風險差異與Urate從對照狀態(tài)變?yōu)樘幚頎顟B(tài)有關(guān)。這里考慮了定義對照和處理條件的兩種可能性。通過這種方式，研究者可以探索因果效應的大小如何取決于所考慮的處理變量的值。

首先，考慮Urate單位變化的因果效應，以總體平均值為中心。然后，在潛在結(jié)果表示法中，感興趣的因果效應是預期風險差異，見式（1）：

式（1）中，μ是Urate的總體平均值。在因果關(guān)系定義中，對照條件定義為低于總體平均Urate的半個單位，處理條件定義為高于總體平均Urate的半個單位。其次，考慮Urate中一個標準差變化的因果效應，也以總體平均值為中心。因果效應現(xiàn)在定義為預期風險差異，見式（2）：

式（2）中，σ是Urate的總體標準差。

根據(jù)前面提到的真實數(shù)據(jù)生成過程，通過生成大量潛在結(jié)果（100 000 000次重復）來計算兩個群體因果效應。通過該方法，總體效應UnitEff為0.007 6，標準化總體效應StdEff為0.006 8。這些值是研究者根據(jù)隨機樣本估計的目標因果效應。

2.3.2 列出可用于識別因果效應的調(diào)整集

給定數(shù)據(jù)的因果圖模型，研究者可以使用因果圖過程分析變量Urate對CVD因果效應的可識別性。以下程序使用該過程列出可用于識別此因果效應的有效調(diào)整集。為簡潔起見，使用MAXSIZE=2選項僅構(gòu)造不超過兩個元素的調(diào)整集。設(shè)所需要的SAS程序如下：

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

該例產(chǎn)生的調(diào)整集列表如表3所示。

表3 模型中可能的調(diào)整集Table 3 Possible adjustment sets for the model

【表3中有關(guān)內(nèi)容的說明】第2列的“大小”指協(xié)變量的個數(shù)（各行均有2個）；第3列的“最小”指所找到的調(diào)整集是否為最小的調(diào)整集（各行上的調(diào)整集包含2個協(xié)變量，均為最小的調(diào)整集）。

請注意，表3中不顯示空集。這意味著變量Urate和CVD之間的邊際關(guān)聯(lián)不能用來估計具備有效因果解釋的因果效應。相反，研究者必須使用另一種估算策略，例如，使用表3中的一個調(diào)整集的逐步調(diào)整估算。如本例后面的內(nèi)容所示，未能執(zhí)行此類調(diào)整會導致對因果效應的有偏估計。

研究者可以使用表3中的任何調(diào)整集來獲得變量Urate對CVD影響的估計，該估計具備有效的因果解釋。集合｛Smoking，StatinUse｝是一個有效的調(diào)整集，它還有一個特性，即集合中的兩個變量都是二值分類變量。因此，估計因果效應的一種可能方法是根據(jù)這兩個變量的水平進行分層分析。

2.4 平均處理效應或預期風險差異

2.4.1 基本情況

目前正在估計兩種因果效應。一個是Urate對CVD的未標準化單位效應，表示為UnitEff，另一個是Urate對CVD的標準化單位效應，表示為StdEff。這兩種因果效應都是根據(jù)預期CVD潛在結(jié)果值的差異來定義的，在某些Urate處理和對照水平上評估這些潛在結(jié)果，這些處理和對照水平是根據(jù)總體參數(shù)定義的。由于這些總體參數(shù)以及處理和對照水平未知，故需要從樣本中估計它們。

2.4.2 計算處理水平與對照水平的樣本值

下面的程序從本例前面創(chuàng)建的匯總統(tǒng)計表中計算Urate處理和對照水平的兩組樣本值。這些計算值存儲在數(shù)據(jù)集ScoreData中，研究者將使用該數(shù)據(jù)集來估計兩個因果效應。設(shè)所需要的SAS程序如下：

以下程序執(zhí)行Logistic回歸分析，該分析按因果圖過程結(jié)果建議的兩個調(diào)整變量的水平分層：

【SAS輸出解釋】

因篇幅所限，具體的輸出結(jié)果從略?，F(xiàn)對其主要內(nèi)容概要解說如下：

在上述兩個二值調(diào)整變量“｛Smoking，StatinUse｝”產(chǎn)生的4個層中進行了分析。在每個層中，可以通過UnitTreat和UnitControl之間的P_1差值計算未標準化單位效應，也可以通過StdTreat和StdControl之間的P_1差值計算標準化效應。然而，層內(nèi)的這些效應都不是因果效應估計值本身。必須使用層中UnitTreat和UnitControl之間P_1差值的加權(quán)平均值來計算因果效應UnitEff的估計值，其中，權(quán)重是層的樣本量。同樣，必須使用層中StdTreat和StdControl之間P_1差值的加權(quán)平均值來計算因果效應StdEff的估計值（注意：后文中表5中的分層估計列顯示了因果效應的這些估計）。

如前所述，如果研究者使用兩個變量之間的邊際關(guān)聯(lián)（即未調(diào)整）來估計變量Urate對CVD的效應，那么混雜的協(xié)變量會使估計結(jié)果產(chǎn)生偏差。嚴格地說，為了呈現(xiàn)這種有偏差的結(jié)果，以下PROC LOGISTIC過程步執(zhí)行不按任何協(xié)變量分層的Logistic回歸分析。設(shè)所需要的SAS程序如下：

proc logistic data=CVDdata noprint；

model CVD（event='1'）=Urate；

score data=ScoreData out=ProbNaive；

run；

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

預期CVD值的相應估計值如表4所示。

表4 未調(diào)整的后驗概率Table 4 Unadjusted posterior probabilities

由表4可知，研究者有兩組估計結(jié)果。一組結(jié)果是通過使用處理變量和結(jié)果變量之間的原始邊際關(guān)聯(lián)來計算的（見表4前兩行）；另一組結(jié)果是基于調(diào)整策略的分層估計器計算的（見表4后兩行）。

兩個估計器計算的因果效應UnitEff和StdEff的估計值如表5所示。

表5 因果效應估計的匯總Table 5 Causa effect estimation summary

由表5可知，使用分層估計計算的估計值（Stratified Estimation）非常接近真實值（True Effect）。因為集合｛Smoking，StatinUse｝是圖1所示數(shù)據(jù)生成過程的有效調(diào)整集（見前文表3最后一行）。然而，基于未經(jīng)調(diào)整的原始數(shù)據(jù)，使用Logistic回歸分析得到的估計值（Unadjusted Estimation）與True Effect不一致。因為基于因果圖模型分析的結(jié)果（見前文表3）表明，空集（指Logistic回歸模型中未列入任何協(xié)變量）不是有效的調(diào)整集。因此，這個例子表明，因果圖理論可以在混淆情況下識別因果效應；通過設(shè)計因果效應的分層估計，該例還展示了如何基于因果圖過程的識別結(jié)果實施良好的統(tǒng)計估計策略。

3 討論與小結(jié)

3.1 討論

采用多重回歸分析方法處理資料的前提是研究者已經(jīng)收集了各變量的具體數(shù)據(jù)［7-8］，而采用因果圖過程進行分析時，不需要提供各變量的具體數(shù)據(jù)，只需要研究者依據(jù)基本常識、專業(yè)知識和以往的研究經(jīng)驗對各變量之間的關(guān)系作出比較合理的設(shè)定，并將其呈現(xiàn)在因果圖上。由此可知，科學合理地運用因果圖過程，有助于探索出多因素研究課題中可能存在的協(xié)變量集合，從而為多因素多指標的研究課題的科研設(shè)計奠定良好基礎(chǔ)。

3.2 小結(jié)

本文介紹了因果圖過程的5個局限性，包括：①因果圖過程不能處理有向循環(huán)的因果圖模型；②因果圖過程不能評估動態(tài)處理方案；③因果效應識別是一個總體概念；④因果效應識別是一個非參數(shù)概念；⑤因果圖過程不能識別某些因果圖模型中的因果效應。同時，本文針對一個實例并基于SAS軟件，實現(xiàn)了用調(diào)整集估計數(shù)據(jù)的因果效應的目的。