圖模型方法用于二值變量相關性分析中缺失數(shù)據(jù)的估計*

楊金英 崔朝杰

在社會調(diào)查、經(jīng)濟研究和醫(yī)學試驗中，人們會經(jīng)常遇到缺失數(shù)據(jù)的情況，如何對缺失數(shù)據(jù)進行補值，長久以來備受統(tǒng)計界關注。本文將通過建立圖模型對含有缺失數(shù)據(jù)的兩值戒煙訓練數(shù)據(jù)進行分析。

文獻〔1〕分析了醫(yī)學上對吸煙者采取不同的戒煙措施的效果，對缺失數(shù)據(jù)提出了三種補值方法，并利用優(yōu)比(OR)值進行了敏感性分析。

這項研究中共有489名吸煙志愿者自愿參加，共進行了四次觀測(干預后，6個月后，12個月后，24個月后)。參加試驗的個體隨機分為三個組:對照組，社會支持Ⅰ組和社會支持Ⅱ組。然而被隨機分到社會支持Ⅰ組和社會支持Ⅱ組中的個體，試驗后期電話通知他們聚會，大約有一半從來沒有出現(xiàn)過，這就導致了數(shù)據(jù)缺失。把社會支持Ⅰ組和社會支持Ⅱ合并為一組，統(tǒng)稱為處理組。對缺失數(shù)據(jù)提出了三種補值方法(missing=smoking，last observation carried forward(LOCF)，a little multiple imputation)。主要分析了兩個時刻的觀測數(shù)據(jù)，即干預后的觀察時刻(時刻1)，24個月后的觀察時刻(時刻2)。為了使數(shù)據(jù)表述方便，我們引入下面的記號:

(t=1，2)，用以上記號可將在時刻2響應，時刻1、時刻2都吸煙且處于處理組的試驗個體頻數(shù)表示為nR2Y1Y2X=n1111，其他情況可類似表示。這樣試驗數(shù)據(jù)可匯總在表1中。

表1 數(shù)據(jù)匯總情況表

針對時刻2的缺失數(shù)據(jù)，用前述三種補值辦法，并就優(yōu)比(OR)的不同取值進行敏感性分析，結果如表2所示。

表2 Group by smoke analyses under different missing data assumptions

這里

表示吸煙與分組處理之間關系的強弱，OR越大相關關系越強。Marginal為所有個體時刻2的邊緣分布，Stratified1為考慮時刻1的分層分布(LOCF補值)，Stratified2為考慮時刻2的分層分布(a little multiple imputation)。

由于數(shù)據(jù)缺失比例較大，23.93%幾乎為個體數(shù)目的四分之一，從上述分析可以看出統(tǒng)計分析結果受到缺失數(shù)據(jù)補值方法的影響很大，所以補值方法的好壞直接影響著統(tǒng)計分析結果的可信度。

在文獻〔1〕提出的三種補值辦法中，方法1(Missing=Smoking)是比較冒進的、理想化的處理辦法，認定缺失的個體更有可能吸煙，而且吸煙會導致個體缺失，顯然這種方法對處理組明顯有利。而方法2(LOCF)認為丟失時刻的吸煙狀態(tài)與最后觀測時刻的吸煙狀態(tài)一致，這種假設也是不符合情理的。上述兩種方法都會導致數(shù)據(jù)分析結果是有偏的。文獻〔1〕通過設定不同的OR值進行了敏感性分析，但是該文作者沒有指出OR在這項研究中到底取何值最為合適。

針對文獻〔1〕提出的補值方法及研究結果，我們所關心的問題是(1)缺失的個體吸煙的可能性是否更大?(2)時刻1和時刻2的吸煙狀況是否相關?(3)戒煙效果是否和個體所在的組(處理組或?qū)φ战M)相關?針對以上問題，本文對含有缺失數(shù)據(jù)的兩值吸煙數(shù)據(jù)再次進行分析，具體方法是:首先對數(shù)據(jù)缺失機制建立了三個圖模型，討論其可識別性之后，又確定了三個圖模型完全數(shù)據(jù)的聯(lián)合密度函數(shù);其次用EM算法對缺失數(shù)據(jù)進行補值，對參數(shù)進行估計;最后利用補值前后的數(shù)據(jù)評估考察缺失與吸煙的關系、時刻1與時刻2吸煙狀況的關系、戒煙效果與個體所在組之間的關系。

對缺失數(shù)據(jù)機制建立圖模型

假設不考慮時刻1的吸煙狀況，只考慮分組、時刻2是否響應及其是否吸煙之間的關系，則可建立圖(1)、圖(2)所示的兩種圖模型:

在圖(1)中，個體響應與否依賴于其是否吸煙，并且給定吸煙狀態(tài)后，個體響應與否與分組無關，則變量X、Y2、R2的聯(lián)合分布為 P(R2，Y2，X)=P(R2|Y2)P(Y2|X)P(X)或 P(R2|Y2)P(X|Y2)P(Y2)。

假設現(xiàn)在我們考慮時刻1的吸煙狀態(tài)，因為在時刻1時所有個體都是可觀測的，個體在此時刻響應與否不依賴于Y1的取值，則Y1與響應示性變量R1中間沒有邊。而在時刻2，個體響應與否僅依賴于Y2的取值，可建立圖模型(3)，則針對圖模型(3)，變量Y1、Y2、R1、R2的聯(lián)合分布 P(Y1，Y2，R1，R2)=P(Y2|Y1)P(R2|Y2)P(Y1)P(R1)

由文獻〔2〕中的引理1易證圖模型(1)、(2)、(3)均是可識別的，此處就不詳細討論了。

下面針對三個圖模型來確定完全數(shù)據(jù)的聯(lián)合密度函數(shù)，三個圖模型中參數(shù)的設定見表4。

在圖模型(1)、(2)中，假設完全數(shù)據(jù)頻數(shù)nijk服從參數(shù)為(n;P(R2=i，Y2=j，X=k))的多項分布，觀測頻數(shù)記為:

則完全數(shù)據(jù)的聯(lián)合概率密度函數(shù)為

在圖模型(3)中，假設完全數(shù)據(jù)頻數(shù)nijk服從參數(shù)為(n;P(R2=i，Y1=j，Y2=k))的多項分布，觀測頻數(shù)記為:

為潛在的缺失數(shù)據(jù)。

為處理組中可觀測數(shù)據(jù)。

為對照組中可觀測數(shù)據(jù)。

則處理組、對照組完全數(shù)據(jù)的聯(lián)合概率密度函數(shù)分別為

表3 缺失數(shù)據(jù)的估計值

用EM算法進行補值和參數(shù)估計

EM算法是對缺失數(shù)據(jù)補值、參數(shù)估計的一種強有利的工具〔3－5〕。針對三個圖模型，我們通過SAS語言編程用EM算法迭代直至設定的精度后，對缺失數(shù)據(jù) x1、x2、x3、x4進行補值結果如表3 所示，參數(shù) αi(i=1，2，…，k)估計結果如表4 所示。

表4 參數(shù)的估計值

結果分析

導致數(shù)據(jù)缺失可能和處理有關，也可能和當時的吸煙狀態(tài)有關，下面我們就來探討數(shù)據(jù)缺失與哪個因素相關。

首先針對時刻2響應與否與分組的二維聯(lián)列聯(lián)表:

=0處理組X觀測R2=1 缺失R2=1 156 34對照組X =0 216 83

通過列聯(lián)表獨立性檢驗，計算得χ2=6.8218，P=0.0127，拒絕缺失和分組是獨立的，說明在不同的組中數(shù)據(jù)缺失的頻率是不同的，對照組的缺失頻率27.76%高于處理組的缺失頻率17.89%，顯然對照組的個體更易于缺失。

其次，我們根據(jù)補值后的數(shù)據(jù)來看吸煙與缺失的二維列聯(lián)表:

294 23.65不吸煙Y2=0=0吸煙Y 2=1觀測R2=1 缺失R2 78 83

通過列聯(lián)表獨立性檢驗，計算得χ2=122.3，P＜0.001，拒絕缺失和吸煙是獨立的，且吸煙個體缺失的頻率7.45%低于不吸煙個體缺失頻率51.55%，這說明多數(shù)的吸煙者還是愿意參加戒煙訓練，渴望戒掉煙癮的。

接下來，我們討論分組的影響。

在時刻1，所有的個體吸煙狀態(tài)都可以觀測得到，數(shù)據(jù)沒有丟失，此時我們考察分組與吸煙狀態(tài)的獨立性，得到二維列聯(lián)表如下:

110 229不吸煙Y2=0=0吸煙Y 2=1處理組X=1 對照組X 80 70

通過列聯(lián)表的獨立性檢驗，計算得χ2=18.288，P＜0.001，所以在初次干預后，拒絕分組和吸煙狀態(tài)是獨立的。在初次干預之后處于對照組的吸煙個體頻率76.59%要高于處理組中的吸煙個體頻率57.89%，說明初次干預效果顯著。

在試驗即將結束時(也就是在24個月后)，由于一部分個體沒有參加聚會，此時我們得不到那些個體是否吸煙的信息，從而導致數(shù)據(jù)的缺失。首先在可以觀測到的372個個體中，僅用不完全數(shù)據(jù)來分析，得到的分組與吸煙狀態(tài)的二維列聯(lián)表如下:

118 176不吸煙Y 2=0=0吸煙Y2=1處理組X=1 對照組X 38 40

通過列聯(lián)表的獨立性檢驗，χ2=1.9506，P=0.1721，檢驗結果無統(tǒng)計學意義，此時我們不能拒絕吸煙和分組是獨立這一假設。

把缺失個體通過EM算法補值的數(shù)據(jù)加上，得到的分組與吸煙狀態(tài)的二維列聯(lián)表如下:

131.95 197.202不吸煙Y 2=0=0吸煙Y2=1處理組X=1 對照組X 58.05 101.789

通過列聯(lián)表的獨立性檢驗，χ2=0.6456，P=0.415，同樣不能拒絕分組和吸煙是獨立的假設。

如果我們在t1時刻分層，加上個體在時刻1是否吸煙的信息，再次觀察個體在不同分組情況下戒煙率是否是有區(qū)別，根據(jù)模型(3)中參數(shù)估計和補值后的缺失數(shù)據(jù)，我們得到下列二維列聯(lián)表:

如果個體在初次干預后是不吸煙的:

Y1=1處理組X=0不吸煙Y2=0 吸煙Y2=1 35.102 44.898對照組X =0 27.094 42.906

通過列聯(lián)表的獨立性檢驗，χ2=1.2389，P=0.5213，檢驗結果無統(tǒng)計學意義。

如果個體在初次干預后仍然吸煙:

Y1=1處理組X=1不吸煙Y2=0 吸煙Y2=1 22.948 87.052對照組X =0 74.704 154.296

通過列聯(lián)表的獨立性檢驗，χ2=1.8115，P=0.025，檢驗結果有統(tǒng)計學意義。

由以上分析我們得到了一個非常有趣的結果，就是如果不考慮初次干預后是否吸煙，則分組效果是不顯著的;如果個體在初次干預后就不吸煙了，那么再進行干預也就失去了意義，當然分組效果是不顯著的。但是如果個體在初次干預后仍然吸煙，則時刻2是否吸煙受到分組的影響，而且處理組的吸煙率高于對照組，這說明初次干預沒有戒掉煙癮的吸煙者，對于后面的多次干預是沒有顯著效果的，從長期來看，反倒是對照組中依靠自我控制來戒煙的效果更明顯。

綜合上述分析我們得到的結論是:

第一，缺失與分組有關，位于對照組的個體缺失的比例更大，缺失個體吸煙的可能性低于可觀測到個體，這說明多數(shù)缺失的個體已戒除煙癮，不需要繼續(xù)參加戒煙訓練，同時說明文獻〔1〕中missing=smoking的補值方法是不合理的;

第二，在時刻1干預效果顯著，吸煙與分組有關，處理組的戒煙率高于對照組;如果不考慮初次干預后個體吸煙的狀態(tài)，則無論是對有缺失數(shù)據(jù)的情況，還是對于通過EM算法補值后得到完全數(shù)據(jù)的情況，戒煙效果與個體所在組無關;

第三，如果初次干預后個體不吸煙了，則在24個月后，個體是否吸煙也與所在組無關;

第四，如果初次干預后個體仍然吸煙，則在24個月后，個體是否吸煙與所在組有關，處理組的戒煙方法沒有明顯的效果。

回到試驗背景，我們看到初次干預結果是非常重要的，通過自我控制來戒煙的方式從長期來看是更加有效的。這也為今后的戒煙訓練工作提供了一個理論指導。

本文的結論更符合實際背景，這在一定程度上說明本文所建立的圖模型是合理的。與文獻〔1〕相比，本文的補值方法可信度更高，可用于分析其他含有缺失數(shù)據(jù)的兩值數(shù)據(jù)，本文得到的結論對戒煙訓練工作更有參考價值。

1．Hedeker D，Robin J，Demirtas H．Analysis of binary outcomes with missing data:missing=smoking，last observation carried forward，and a little multiple imputation．Society for the Study of Addiction，2007，102:1564-1573．

2．Ma WQ，Geng Z，Hu YH．Identification of graphical models for nonignorable nonresponse of binary outcome in longitudinal studies．Journal of Multivariate Analysis，2003，87:24-45．

3．Little RJA，Rubin DR．Statistical analysis with missing data．New York:Wiley，1987．

4．Dempster AP，Laird NM，Rubin DB．Maximum likelihood from incomplete data via the EM algorithm．Journal of the Royal Statistical Society，Series B，1977，39:1-38．

5．趙志文，王思洋，王瑞庭，等．定時截尾下具有部分缺失數(shù)據(jù)兩個指數(shù)總體參數(shù)估計與檢驗．吉林大學學報(理學版)，2009，47(1):26-30．