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

楊金英 崔朝杰

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

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

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

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

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

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

表2 Group by smoke analyses under different missing data assumptions

這里

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

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

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

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

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

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

在圖(1)中，個(gè)體響應(yīng)與否依賴于其是否吸煙，并且給定吸煙狀態(tài)后，個(gè)體響應(yīng)與否與分組無關(guān)，則變量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)。

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

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

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

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

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

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

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

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

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

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

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

用EM算法進(jìn)行補(bǔ)值和參數(shù)估計(jì)

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

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

結(jié)果分析

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

首先針對(duì)時(shí)刻2響應(yīng)與否與分組的二維聯(lián)列聯(lián)表:

=0處理組X觀測(cè)R2=1 缺失R2=1 156 34對(duì)照組X =0 216 83

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

如果個(gè)體在初次干預(yù)后是不吸煙的:

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

通過列聯(lián)表的獨(dú)立性檢驗(yàn)，χ2=1.2389，P=0.5213，檢驗(yàn)結(jié)果無統(tǒng)計(jì)學(xué)意義。

如果個(gè)體在初次干預(yù)后仍然吸煙:

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

通過列聯(lián)表的獨(dú)立性檢驗(yàn)，χ2=1.8115，P=0.025，檢驗(yàn)結(jié)果有統(tǒng)計(jì)學(xué)意義。

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

綜合上述分析我們得到的結(jié)論是:

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

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

第三，如果初次干預(yù)后個(gè)體不吸煙了，則在24個(gè)月后，個(gè)體是否吸煙也與所在組無關(guān);

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

回到試驗(yàn)背景，我們看到初次干預(yù)結(jié)果是非常重要的，通過自我控制來戒煙的方式從長(zhǎng)期來看是更加有效的。這也為今后的戒煙訓(xùn)練工作提供了一個(gè)理論指導(dǎo)。

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

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í)截尾下具有部分缺失數(shù)據(jù)兩個(gè)指數(shù)總體參數(shù)估計(jì)與檢驗(yàn)．吉林大學(xué)學(xué)報(bào)(理學(xué)版)，2009，47(1):26-30．