不同缺失機(jī)制并存時偏倚校正的模擬研究*

趙俊康 王 彤 榮惠英 孟繁龍

弱勢人群的醫(yī)療救助問題一直以來備受世界各國政府關(guān)注[1]。這部分弱勢人群的特點是收入偏低，極易陷入因貧致病和因病致貧的惡性循環(huán)中。根據(jù)第四次國家衛(wèi)生服務(wù)調(diào)查結(jié)果[2]，我國約有38%的居民生病不去看醫(yī)生，經(jīng)醫(yī)生診斷該住院治療而未住院的達(dá)21%，其中70.3%的人未住院的主要原因仍然是“經(jīng)濟(jì)困難”。這些潛在患者不選擇就醫(yī)使得從醫(yī)院收集數(shù)據(jù)仍然很難估計出這部分非從業(yè)人群的全部醫(yī)療費(fèi)用需求。低收入非從業(yè)人群更可能由于貧困等原因得到高于普通人群的致病機(jī)會，因此若用一個總的平均水平來估計弱勢人群醫(yī)療費(fèi)用的實際需求將明顯低估這種需求。

這種由于個體自我行為(因經(jīng)濟(jì)困難自主選擇不就醫(yī))所導(dǎo)致的樣本選擇偏倚，單靠好的抽樣設(shè)計是無法消除的。需要注意的是患病但自我選擇未就醫(yī)者應(yīng)答表現(xiàn)出的0消費(fèi)與真正未患病而不就醫(yī)者的消費(fèi)真值是不同的，即自我選擇未就醫(yī)者的醫(yī)療消費(fèi)真值未知，應(yīng)視為缺失數(shù)據(jù)[3]。將這類真值未知的0消費(fèi)數(shù)據(jù)刪除或者直接取因變量為0來應(yīng)用多元線性回歸等常規(guī)的統(tǒng)計學(xué)分析方法就忽視了這種無應(yīng)答偏倚；同時，像這類較大規(guī)模的社會學(xué)或流行病學(xué)調(diào)查中無應(yīng)答偏倚也是常態(tài)而不是偶然[4]，故而針對不同缺失機(jī)制下的無應(yīng)答偏倚探討其校正方法成為國內(nèi)外學(xué)者長期以來關(guān)注的問題。

Rubin等人于1976年提出的缺失機(jī)制主要包括完全隨機(jī)缺失MCAR(missing completely at random)、隨機(jī)缺失MAR (missing at random)和非隨機(jī)缺失NMAR(not missing at random)三類[5]。在MCAR假定下，對完全觀測個體使用的任何分析方法仍然有效；在MAR假定下,主流觀點是采用多重填補(bǔ)MI(multiple imputation)對隨機(jī)缺失進(jìn)行填補(bǔ)繼而得出無偏估計[6]；而對于由于自主選擇不就醫(yī)而導(dǎo)致的NMAR，本研究選用適合于該類型數(shù)據(jù)的受限因變量(limited dependent variable)統(tǒng)計模型來進(jìn)行校正[7-9]。

隨機(jī)缺失機(jī)制下的多重填補(bǔ)方法

1.多重填補(bǔ)的具體步驟

多重填補(bǔ)(MI)主要由三個獨立的步驟組成：填補(bǔ)階段、分析階段和合并階段。MI其實是包含了一組方法的一個廣義的術(shù)語，在其框架內(nèi)的所有的方法中都含有這三步過程。圖1描述了整個過程。

填補(bǔ)階段 分析階段 合并階段

圖1 MI的三個步驟

(1)填補(bǔ)階段為每個缺失值抽取m個估計值進(jìn)行填補(bǔ)，從而構(gòu)成m個完整數(shù)據(jù)集，這m個數(shù)據(jù)集中只有觀測數(shù)據(jù)是相同的，填補(bǔ)值一般不等。(2)分析階段：分析步的主要分析對象就是填補(bǔ)好的數(shù)據(jù)集，這一步將應(yīng)用數(shù)據(jù)原本完整時所用到的相同的方法來分析。唯一的區(qū)別在于要對每個完整數(shù)據(jù)集分別使用該方法處理，因此將分析m次。對于本研究中包含自主選擇性缺失的醫(yī)療費(fèi)用數(shù)據(jù)，這一步可使用填補(bǔ)后的多個數(shù)據(jù)集與選擇性偏倚導(dǎo)致的因變量為虛假0(視為缺失)的數(shù)據(jù)合并進(jìn)行樣本選擇模型分析以解決選擇性偏倚導(dǎo)致的那部分缺失，繼而得出m個樣本選擇模型擬合結(jié)果。(3)合并階段：綜合這m個擬合結(jié)果，根據(jù)Rubin(1987)提出的針對參數(shù)估計值與標(biāo)準(zhǔn)誤的合并準(zhǔn)則[10]，最終得到對目標(biāo)變量的統(tǒng)計推斷。

2.填補(bǔ)模型

(1)預(yù)測均數(shù)匹配法(PMM)

預(yù)測均數(shù)匹配法(PMM)是處理單調(diào)缺失模式中定量變量缺失的多重填補(bǔ)方法之一，PMM法的具體填補(bǔ)步驟如下：

令Yj為有缺失值的定量變量，用Yj及其協(xié)變量X1,X2,…,Xk均被觀測到的個體觀察值建立回歸模型：

Yj=β0+β1X1+β2X2+…+βkXk

(1)

每次填補(bǔ)產(chǎn)生填補(bǔ)值的步驟如下：

(2)

其中nj是變量Yj未缺失的觀測個體數(shù)。

然后，從標(biāo)準(zhǔn)正態(tài)分布N(0,1)中抽取k+1個獨立的變量，組成一個有k+1個元素的向量Z，得到新的回歸系數(shù)：

(3)

②對于每個缺失值，其預(yù)測值為

(4)

④最后從這k0個觀測值中隨機(jī)抽取一個值填補(bǔ)缺失值。

(2) 傾向性得分法

傾向性得分(PS)是指對給定的觀察到的協(xié)變量條件下，每個觀察值被分配到某特定處理組的條件概率。PS法的具體填補(bǔ)步驟如下：

①產(chǎn)生一個指示變量Rj，當(dāng)Rj為0時，表示變量Yj中有缺失值的個體；當(dāng)Rj為1時，則表示變量Yj中被觀測到的個體。

②擬合logistic回歸模型

logit(pj)=β0+β1X1+β2X2+…+βkXk

(5)

其中X1,X2,…,Xk是Yj的協(xié)變量，

pj=Pr(Rj=0|X1,X2,…,Xk)

logit(pi)=log(pi/(1-pi))

③根據(jù)模型計算變量Yj上每個個體數(shù)據(jù)缺失的傾向性得分logit(pj)，并根據(jù)該得分將所有的觀測分組，一般為5組，如果觀測數(shù)量較多，可分為更多的組。

⑤重復(fù)以上步驟，直到每一個缺失變量都得到填補(bǔ)。

(3)基于Bootstrap的EM算法

基于Bootstrap的EM算法(EMB)是通過bootstrap算法從參數(shù)的后驗密度中抽取新的參數(shù)，從而代替其他方法中復(fù)雜的抽取過程。該方法不必估計參數(shù)的方差矩陣，也不用像期望最大化重要性抽樣EMis(expectation maximization importance sampling)算法那樣進(jìn)行重要性抽樣，甚至不需要像數(shù)據(jù)增廣DA (data augmentation ) 算法一樣推導(dǎo)Markov鏈并檢查收斂性，而且還可以應(yīng)用于非常大型的數(shù)據(jù)[11]。

EM算法是一種迭代算法，廣泛運(yùn)用于尋找參數(shù)的最大似然估計值，尤其是在缺失數(shù)據(jù)的問題中非常有用[12]。它的每一次迭代都由兩步組成：E步(求期望)和M步(極大化)。其中E步是在給定已觀測到的數(shù)據(jù)和當(dāng)前參數(shù)下，求缺失數(shù)據(jù)的條件期望，然后用這些條件期望值去填補(bǔ)缺失值。M步是當(dāng)缺失數(shù)據(jù)被填補(bǔ)之后就像沒有缺失一樣進(jìn)行的極大似然估計。重復(fù)以上兩步,直至前后兩次計算結(jié)果達(dá)到規(guī)定的收斂標(biāo)準(zhǔn)。

將經(jīng)典的EM算法和Bootstrap算法相結(jié)合來從后驗分布中進(jìn)行抽取。具體填補(bǔ)步驟如下：

①從含量為n的完整數(shù)據(jù)集Yobs中有放回地抽取m個大小為n的樣本。

②對每一個樣本運(yùn)行穩(wěn)定而快速的EM算法，得到一組參數(shù)的點估計值μ和∑，共m組。

③用Ycom中的觀測數(shù)據(jù)分別結(jié)合每組參數(shù)估計值得出缺失數(shù)據(jù)的條件分布，并從中抽取填補(bǔ)值，繼而得到m個經(jīng)填補(bǔ)后的完整數(shù)據(jù)集。

(4)Markov Chain Monte Carlo方法

Markov Chain Monte Carlo(MCMC)方法在應(yīng)用于缺失值領(lǐng)域時稱為數(shù)據(jù)增廣DA算法[13]，同EM算法一樣，DA算法也是依次填補(bǔ)缺失值和推斷未知參數(shù)的一種迭代方法。區(qū)別就在于DA是以隨機(jī)的方式對缺失值和參數(shù)進(jìn)行抽取，而EM算法只是缺失值和參數(shù)的點估計。

該算法通過填充(imputation)及后驗(posterior)兩步迭代來實現(xiàn)：

①填補(bǔ)步(I-step)

填補(bǔ)的數(shù)據(jù)是從給定觀測數(shù)據(jù)、均數(shù)和協(xié)方差矩陣后的缺失數(shù)據(jù)的條件分布中隨機(jī)抽取得到，從貝葉斯角度來看，該分布又稱為后驗預(yù)測分布。

在每次迭代過程中，從事先給定的均值向量μ和協(xié)方差陣Σ的初始估計值開始，在給定Yobs下的條件分布P(Ymis|Yobs)中抽取Ymis。

在某個觀測個體具有類似的缺失模式時，令Yobs=y1，就得到均值向量及條件協(xié)方差矩陣分別為

(10)

(11)

的多元正態(tài)分布P(Ymis|Yobs=y1)，也就是Ymis的條件分布。整個I-step可以表述成：

(12)

②后驗步(P-step)

由于在多重填補(bǔ)過程中需要產(chǎn)生多個完整的數(shù)據(jù)集，因此在每個I-step需要不同的均數(shù)向量和協(xié)方差矩陣，因此P-step的目的就是輪流產(chǎn)生參數(shù)估計值。

(13)

(14)

(15)

這樣從后驗分布中抽取新的參數(shù)后，接下來的I-step使用這些新的參數(shù)值來產(chǎn)生新的填補(bǔ)值；然后新的填補(bǔ)數(shù)據(jù)繼續(xù)用于下一個P-step，繼而再抽取另一組新的參數(shù)估計值。如此循環(huán)往復(fù)重復(fù)這兩個步驟一定的次數(shù)，產(chǎn)生一個足夠長的隨機(jī)序列：

(16)

該隨機(jī)序列是就一條馬爾科夫鏈，并且在一定的正則條件下會收斂到一個穩(wěn)定分布[14]。當(dāng)該鏈?zhǔn)諗康揭粋€穩(wěn)定的分布P(Ymis,θ|Yobs)時，就可以近似獨立地從該分布中為缺失值抽取填補(bǔ)值。

兩種缺失機(jī)制并存時的兩階段校正模擬研究

1.模擬設(shè)計

首先需要構(gòu)建出含結(jié)果等式和選擇等式的樣本選擇模型，如下所示：

y0=x1+ε

(17)

d0=x2+v

(18)

d=1(d0≥c)，d=0(d0

(19)

y1=y0·d

(20)

取10000例觀測值(n=10000)，根據(jù)樣本選擇模型的結(jié)果等式(17)和選擇等式(18)模擬六個變量x1、x2、y0、d0、ε、v，x1和x2分別取自均值為0，標(biāo)準(zhǔn)差為1，相關(guān)系數(shù)為0的雙變量正態(tài)分布，而ε和v取自均值為0，標(biāo)準(zhǔn)差為1，相關(guān)系數(shù)為0.75的雙變量正態(tài)分布。y0和d0分別通過公式y(tǒng)0=e1+x1+ε和d0=1+x2+v求出。第一步首先對全部10000例觀測值的因變量分別進(jìn)行輕度、中度和重度截取，即以5%、30%和70%的比例向下截取產(chǎn)生對應(yīng)于調(diào)查中生病但自主選擇不就醫(yī)者發(fā)生的虛假0消費(fèi)因變量?？梢酝ㄟ^以上截取比例給y0定義一個相應(yīng)的界值c，當(dāng)d0>c時，令y1=y0且d=1；當(dāng)d0≤c時，令y1為缺失且d=0,最后對y1進(jìn)行對數(shù)轉(zhuǎn)換，對應(yīng)于使醫(yī)療費(fèi)用值近似服從正態(tài)分布。通過調(diào)整c的值就可以獲得針對樣本選擇模型不同程度缺失率的非隨機(jī)缺失數(shù)據(jù)。第二步分別在上述三個不同缺失率下，令d=1的個體以5%、30%和70%的比例隨機(jī)產(chǎn)生缺失(對應(yīng)于調(diào)查中可能發(fā)生的隨機(jī)缺失)。這樣就產(chǎn)生9種不同的組合數(shù)據(jù)。在上述不同組合下，首先對d=1的個體(僅存在隨機(jī)缺失)分別應(yīng)用PMM、PS、MCMC和EMB法進(jìn)行多重填補(bǔ)，然后把填補(bǔ)后的數(shù)據(jù)與d=0的數(shù)據(jù)合并，應(yīng)用樣本選擇模型的兩步似然估計來獲得各自的回歸系數(shù)估計量來校正虛假0消費(fèi)產(chǎn)生的選擇性偏倚。最后，重復(fù)抽樣100次,計算9種組合下兩階段校正方法所獲得結(jié)果等式中自變量x1的回歸系數(shù)和標(biāo)準(zhǔn)誤。本次模擬分析中，多重填補(bǔ)技術(shù)的EMB算法選用了R軟件，PMM、PS和MCMC法選用了SAS軟件，樣本選擇模型分析也選用了SAS軟件。

2．評價標(biāo)準(zhǔn)

在比較兩階段校正方法下樣本選擇模型結(jié)果等式回歸系數(shù)的優(yōu)劣時，選用以下三個評價標(biāo)準(zhǔn)[15-16]。

(1) 標(biāo)準(zhǔn)偏倚(Standardized bias)

當(dāng)標(biāo)準(zhǔn)偏倚落在±0.4區(qū)間之外時，偏倚就會對功效、可信區(qū)間覆蓋率和誤差率產(chǎn)生明顯的負(fù)面影響。因此，將±0.4作為評價標(biāo)準(zhǔn)偏倚的上下界值，即若某方法的標(biāo)準(zhǔn)偏倚絕對值超出0.4，此方法便無法接受。標(biāo)準(zhǔn)偏倚做為評價準(zhǔn)確度的指標(biāo)是方法評價指標(biāo)中的首要觀測指標(biāo)。

(2)可信區(qū)間平均長度(length)

如果一個方法與另一個方法相比，有相同的或更高的準(zhǔn)確度，但得出的可信區(qū)間平均卻更短，那么此方法的精確度就更高。

(3)均方誤差的平方根(RMSE)

3．模擬分析結(jié)果

表1 四種填補(bǔ)方法下樣本選擇模型結(jié)果等式的回歸系數(shù)估計值的各項評價標(biāo)準(zhǔn)比較(一)

從表1到表3可知，四種填補(bǔ)方法的標(biāo)準(zhǔn)偏倚絕對值均不等，其中PS法超過了所規(guī)定的界值，故該法效果相對不理想；其余三個方法中，均方誤差的平方根和可信區(qū)間平均長度均相差不大，因此根據(jù)標(biāo)準(zhǔn)偏倚絕對值大小便可判斷出不同缺失機(jī)制組合下的填補(bǔ)方法優(yōu)劣。

綜上，各種情況下不同方法的推薦結(jié)果如表4。

表2 四種填補(bǔ)方法下樣本選擇模型結(jié)果等式的回歸系數(shù)估計值的各項評價標(biāo)準(zhǔn)比較(二)

表3 四種填補(bǔ)方法下樣本選擇模型結(jié)果等式的回歸系數(shù)估計值的各項評價標(biāo)準(zhǔn)比較(三)

表4 不同缺失機(jī)制組合下的填補(bǔ)方法選擇

討 論

數(shù)據(jù)缺失現(xiàn)象在調(diào)查研究中非常普遍，它不僅會降低參數(shù)估計的效率，同時也給統(tǒng)計分析帶來很大偏倚。根據(jù)數(shù)據(jù)缺失機(jī)制，可將數(shù)據(jù)缺失分為三類：完全隨機(jī)缺失(MCAR)、隨機(jī)缺失(MAR)和非隨機(jī)缺失(NMAR)。針對MAR機(jī)制，統(tǒng)計學(xué)家們提出了多種方法來校正這種缺失帶來的偏倚，MI就是被廣為推崇的方法之一；針對NMAR機(jī)制，由于該機(jī)制的復(fù)雜性，當(dāng)前還沒有一種統(tǒng)一的方法來校正這種偏倚，不過當(dāng)回歸模型中的應(yīng)變量為非隨機(jī)缺失時，某些情況下可以應(yīng)用樣本選擇模型來糾正這種NMAR帶來的偏倚；但當(dāng)兩種缺失機(jī)制并存時的偏倚糾正方法尚未見有介紹。對此，本研究提出了一個兩階段策略以糾正不同缺失機(jī)制造成的偏倚。第一階段首先利用只包含隨機(jī)缺失數(shù)據(jù)的個體對單純無應(yīng)答缺失按照MAR機(jī)制進(jìn)行多重填補(bǔ)，在第二階段中使用填補(bǔ)后的多個數(shù)據(jù)集與選擇性偏倚導(dǎo)致的因變量為虛假0(視為缺失)數(shù)據(jù)合并進(jìn)行樣本選擇模型分析以校正由于非隨機(jī)缺失所造成的偏倚，最后對多個樣本選擇模型擬合結(jié)果進(jìn)行合并。

模擬研究結(jié)果表明：當(dāng)非隨機(jī)缺失為輕度時，PS法由于標(biāo)準(zhǔn)偏倚絕對值遠(yuǎn)遠(yuǎn)超過了規(guī)定的界值，所以該法的結(jié)果相對不理想；而MCMC、EMB和PMM法均得出較好的結(jié)果。不同程度隨機(jī)缺失情況下的填補(bǔ)方法選擇為：隨機(jī)缺失也為輕度時，MCMC法最好；隨機(jī)缺失為中度時，EMB法最好；在隨機(jī)缺失為重度時，PMM法最好。

當(dāng)非隨機(jī)缺失為中度時，PS法由于標(biāo)準(zhǔn)偏倚絕對值遠(yuǎn)遠(yuǎn)超過了規(guī)定的界值，所以仍不可取，而MCMC、EMB和PMM法均得出較好的結(jié)果。此時，無論隨機(jī)缺失程度如何，MCMC法都是最好的方法。

當(dāng)非隨機(jī)缺失為重度時，PS法由于標(biāo)準(zhǔn)偏倚絕對值遠(yuǎn)遠(yuǎn)超過了規(guī)定的界值，所以仍不可取，而MCMC、EMB和PMM法均得出較好的結(jié)果。此時，無論隨機(jī)缺失程度如何，PMM法都是最好的方法。

本文以醫(yī)療費(fèi)用調(diào)查研究中可能出現(xiàn)的兩種缺失為假設(shè)背景，探索性地提出兩階段策略糾正這兩種偏倚，希望能為以后在缺失值處理方面的應(yīng)用提供一些方法學(xué)依據(jù)。

參 考 文 獻(xiàn)

1.Fisher ES,Bynum JP,Skinner JS.Slowing the growth of health care costs—lessons from regional variation.New England Journal of Medicine,2009,360(9):849-852.

2.衛(wèi)生部統(tǒng)計信息中心.第四次國家衛(wèi)生服務(wù)調(diào)查主要結(jié)果.[cited 2010年3月16日];http://www.moh.gov.cn/publicfiles/business/htmlfiles/mohbgt/s3582/200902/39201.htm].

3.Baer OCJ.Bradley JC,et al.Testing and correcting for non-random selection bias due to censoring:an application to medical costs.Health Services and Outcomes Research Methodology,2003,4(2):93-107.

4.Peytchev A,Baxter RK,Carley-Baxter LR.Not All Survey Effort is Equal.Public Opinion Quarterly,2009,73(4):785-806.

5.Rubin DB.Inference and missing data.Biometrika,1976,63(3):581-592.

6.Little RJA.Rubin DB.Statistical analysis with missing data.2nd.Vol.2.2002:Wiley New York:2002.

7.薛小平,史東平，王彤.受限因變量模型及其半?yún)?shù)估計.中國衛(wèi)生統(tǒng)計,2007,24(2):211-213.

8.張磊.樣本選擇模型的似然估計與兩步估計.現(xiàn)代預(yù)防醫(yī)學(xué),2007,34(9):1607-1609.

9.張磊.樣本選擇模型及其在醫(yī)療費(fèi)用研究中的應(yīng)用.[碩士學(xué)位論文].山西：山西醫(yī)科大學(xué),2007.

10.Donald B.Multiple Imputation for Nonresponse in Surveys.American Journal of Sociology,1987,76:346.

11.Honaker J,King G.What to Do about Missing Values in Time‐Series Cross‐Section Data.American Journal of Political Science,2010,54(2):561-581.

12.Dempster AP,Laird NM,Rubin DB.Maximum likelihood from incomplete data via the EM algorithm.Journal of the Royal Statistical Society.Series B (Methodological),1977:1-38.

13.Tanner MA,Wong WH.The calculation of posterior distributions by data augmentation.Journal of the American Statistical Association,1987,82(398):528-540.

14.Schunk D.A Markov chain Monte Carlo multiple imputation procedure for dealing with item nonresponse in the German save survey.2007.

15.Co1lins LM，Schafer JL，Kam CM.A Comparison of Inclusive and Restrictive Strategies in Modem Missing Data Procedures.Psychological methods,2001,6(4):330.

16.Burton A,et al.The design of simulation studies in medical statistics.Statistics in Medicine,2006,25(24):4279-4292.