胡純嚴(yán) ，胡良平 ，2*

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

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

1 因果圖過(guò)程的局限性

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

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

1.2 不能評(píng)估動(dòng)態(tài)處理方案

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

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

1.3 因果效應(yīng)識(shí)別是一個(gè)總體概念

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

1.4 因果效應(yīng)識(shí)別是一個(gè)非參數(shù)概念

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

1.5 不能識(shí)別某些因果圖模型中的因果效應(yīng)

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

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

2.1 問(wèn)題與背景信息

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

圖1 血清尿酸鹽對(duì)心血管疾病風(fēng)險(xiǎn)影響的因果圖模型Figure 1 Causal graph model of the effect of serum urate on risk of cardiovascular disease

2.2 分析因果圖模型的思路

2.2.1 基本方法

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

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

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

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

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

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

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

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

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

proc means data=CVDdata；

var Urate；

ods output Summary=SampleMeansOutput；

run；

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

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

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

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

2.3.1 基本情況

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

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

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

式（1）中，μ是Urate的總體平均值。在因果關(guān)系定義中，對(duì)照條件定義為低于總體平均Urate的半個(gè)單位，處理?xiàng)l件定義為高于總體平均Urate的半個(gè)單位。其次，考慮Urate中一個(gè)標(biāo)準(zhǔn)差變化的因果效應(yīng)，也以總體平均值為中心。因果效應(yīng)現(xiàn)在定義為預(yù)期風(fēng)險(xiǎn)差異，見(jiàn)式（2）：

式（2）中，σ是Urate的總體標(biāo)準(zhǔn)差。

根據(jù)前面提到的真實(shí)數(shù)據(jù)生成過(guò)程，通過(guò)生成大量潛在結(jié)果（100 000 000次重復(fù)）來(lái)計(jì)算兩個(gè)群體因果效應(yīng)。通過(guò)該方法，總體效應(yīng)UnitEff為0.007 6，標(biāo)準(zhǔn)化總體效應(yīng)StdEff為0.006 8。這些值是研究者根據(jù)隨機(jī)樣本估計(jì)的目標(biāo)因果效應(yīng)。

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

給定數(shù)據(jù)的因果圖模型，研究者可以使用因果圖過(guò)程分析變量Urate對(duì)CVD因果效應(yīng)的可識(shí)別性。以下程序使用該過(guò)程列出可用于識(shí)別此因果效應(yīng)的有效調(diào)整集。為簡(jiǎn)潔起見(jiàn)，使用MAXSIZE=2選項(xiàng)僅構(gòu)造不超過(guò)兩個(gè)元素的調(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)容的說(shuō)明】第2列的“大小”指協(xié)變量的個(gè)數(shù)（各行均有2個(gè)）；第3列的“最小”指所找到的調(diào)整集是否為最小的調(diào)整集（各行上的調(diào)整集包含2個(gè)協(xié)變量，均為最小的調(diào)整集）。

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

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

2.4 平均處理效應(yīng)或預(yù)期風(fēng)險(xiǎn)差異

2.4.1 基本情況

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

2.4.2 計(jì)算處理水平與對(duì)照水平的樣本值

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

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

【SAS輸出解釋】

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

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

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

proc logistic data=CVDdata noprint；

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

score data=ScoreData out=ProbNaive；

run；

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

預(yù)期CVD值的相應(yīng)估計(jì)值如表4所示。

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

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

兩個(gè)估計(jì)器計(jì)算的因果效應(yīng)UnitEff和StdEff的估計(jì)值如表5所示。

表5 因果效應(yīng)估計(jì)的匯總Table 5 Causa effect estimation summary

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

3 討論與小結(jié)

3.1 討論

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

3.2 小結(jié)

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