轉(zhuǎn)發(fā)：有缺失值怎么辦——多重填補技術(shù)的思路與軟件實現(xiàn) (sohu.com)
多重填補是針對單一填補而言怖现，所謂單一填補耳高，就是對每一缺失值填一個值提佣，這類方法不能反映缺失數(shù)據(jù)的不確定性楚午，因而容易導(dǎo)致標準誤低估。比如上一篇文章說的回歸方法填補矾策，填補后甲喝，由于填補值是根據(jù)直線回歸估計出來的，他們都在一條直線上雳旅，沒有誤差，標準誤肯定低估了间聊，因為正常情況下攒盈，所有的數(shù)據(jù)距離回歸線總是有點誤差的。
而多重填補是模擬生成一個缺失數(shù)據(jù)的隨機分布哎榴，然后從這一分布中隨機抽取數(shù)據(jù)作為缺失值的填補型豁，因而可以反映不確定性。正因為是隨機抽取尚蝌，導(dǎo)致多重填補技術(shù)存在一個副作用：用該法填補數(shù)據(jù)會同時產(chǎn)生多組填補值迎变，而且你的填補值跟我的填補值永遠都不一樣（除非指定同樣的種子數(shù)），換句話說飘言，如果你發(fā)表一篇填補數(shù)據(jù)的文章衣形，即使我用你的程序重新跑一遍，結(jié)果也跟你的不一樣（你只能祈禱結(jié)果差別不要太大）姿鸿。
下面介紹一下MI法的思想谆吴。
設(shè)有缺失值的變量為y，用于預(yù)測填補的變量為x（簡單起見苛预，設(shè)只有一個變量）句狼，利用回歸填補法，可以得到：
y=a+bx
根據(jù)x值可以填補y值热某。但這樣填補后的值腻菇，如果再分析y與x的關(guān)系，會使二者關(guān)聯(lián)人為增大昔馋，標準誤變低（原因前面剛說了）芜繁。
為了解決這一問題，很自然的一個想法就是绒极，把y加上一個隨機波動項骏令，使y和x之間不是嚴格的直線關(guān)系。即把二者關(guān)系變成：
y=a+bx+e
e是一個隨機數(shù)字（通陈⑻幔可以從殘差中隨機抽壤拼）。這樣的話铡俐，y和x之間就不是嚴格的線性關(guān)系凰兑，填補的值就不在y=a+bx的直線上了，而是有一定的偏離（偏離大小為e）审丘。通過這種思路吏够，可以使得標準誤變大一些，更接近實際。（事實上锅知，大家仔細想想線性回歸播急，也就是這么回事，回歸線與數(shù)據(jù)點總是有點差距的售睹，這就是殘差）
上述方式盡管會使標準誤略有增大桩警，但這仍不夠，因為它們都是在a+bx的基礎(chǔ)上添加一個隨機波動昌妹。也就是說捶枢，是把a和b這兩個系數(shù)當做真實參數(shù)，而實際上它們只是樣本估計值飞崖。
我們可以想象一下實際情況：總體的截距和斜率是無法知道的烂叔，只能根據(jù)當前的一次抽樣數(shù)據(jù)獲得樣本截距a和斜率b。理論上固歪，如果重復(fù)10次抽樣长已，應(yīng)該會得到10個a和b，每次計算的a和b肯定會有所不同昼牛，它們都只是圍繞總體截距和斜率波動的隨機值术瓮。
因此，為了更加符合實際情況贰健，還需要把a和b也作為隨機波動胞四。那怎么把a和b作為隨機抽取的值呢？這就需要有一個關(guān)于a和b的分布伶椿，然后從這個分布中隨機抽取a和b辜伟。這一分布通常稱為貝葉斯后驗分布，也就是說脊另，我們要從貝葉斯后驗分布中隨機抽取a和b导狡，每次抽取的都會有所不同，這樣建立的模型也會不同偎痛，從而估計的缺失值也是隨機的旱捧。這相當于用統(tǒng)計的方法來模擬實際，讓得到的結(jié)果最大可能地接近實際情形踩麦。
現(xiàn)在就面臨這樣的問題：如何產(chǎn)生貝葉斯后驗分布并從該分布中隨機抽让渡摹？目前大多數(shù)統(tǒng)計軟件都是用馬爾科夫鏈蒙特卡洛（Markov chain Monte Carlo谓谦，MCMC）模擬的方法贫橙。所謂MCMC，就是兩種技術(shù)的組合反粥，蒙特卡洛主要是用來模擬抽樣卢肃，馬爾科夫鏈主要是用來產(chǎn)生平穩(wěn)的分布鏈疲迂。
比如下圖就是一個MCMC模擬產(chǎn)生的一個分布，缺失值的填補就是從產(chǎn)生的分布中隨機抽取莫湘，而且通常是隨機抽取多個尤蒿，所以是多重填補。
可能有人會說逊脯，我得到多個填補值优质，怎么用到最終的分析呢竣贪？一般是這樣：比如取5個军洼，你可以計算其均值，也可以用這5次填補的結(jié)果重新做統(tǒng)計分析（如回歸）演怎，最后得到一個綜合的結(jié)果匕争。如下圖：
MCMC的算法極其復(fù)雜（幸虧我們有統(tǒng)計軟件），但思路還是比較容易理解的爷耀，主要步驟如下：
第一甘桑，先利用無缺失的數(shù)據(jù)作為初始的均數(shù)和協(xié)方差矩陣，建立一個y對x的回歸方程歹叮，得到回歸系數(shù)跑杭。
第二，將y中的缺失值利用回歸方程填補上咆耿，并對每一填補值加上一個從殘差中隨機抽取的數(shù)值德谅。
第三，利用填補的“完整”數(shù)據(jù)萨螺，重新計算一個新的均數(shù)和協(xié)方差矩陣窄做。
第四，根據(jù)第三步中得到的新的均數(shù)和協(xié)方差矩陣的后驗分布慰技，從中隨機抽取一個均數(shù)和協(xié)方差椭盏。
第五，利用第四步中抽取的均數(shù)和協(xié)方差吻商，再回到第一步掏颊，把剛才步驟再次執(zhí)行一遍。如此多次循環(huán)艾帐，直到最終的收斂（通常收斂的意思是指蚯舱，上一次的結(jié)果與本次結(jié)果差異非常小甚至無差異）。用最后收斂所得到的填補值作為最終填補值掩蛤。
所以枉昏，其實多重填補技術(shù)也并不是難以捉摸，只要仔細想想揍鸟，其實思路還是挺清楚的兄裂。當然句旱，可能有的人更關(guān)心的是到底如何在軟件上實現(xiàn)。下面就給出SAS和R的實現(xiàn)語句晰奖，感興趣的自己練習(xí)一下谈撒。
SAS可通過proc mi和procmianalyze實現(xiàn)EM算法填補、多重填補等：
/*?以下程序產(chǎn)生EM算法的填補結(jié)果*/
proc?mi?data=aa?nimpute=0;
em?itprint?out=emimp;
varycourse ht wt;?/*?要填補的變量*/
run;
proc?print?data=emimp;
run;
/*?以下程序產(chǎn)生多重填補結(jié)果*/
proc?mi?data=aa?nimpute=5?out=outmi;?/*指定產(chǎn)生5次填補結(jié)果*/
mcmc?displayinit?initial=em(itprint)?plots=trace;
/* initial=em(itprint)表示以EM算法作為初始估計值匾南，并輸出迭代結(jié)果*/
vary course ht wt;
run;
/*下面reg過程分別根據(jù)5次填補值輸出回歸結(jié)果啃匿，為后面的mianalyze過程準備*/
proc?reg?data=outmi?outest=outreg?covout;
mode ly=course ht wt;
by?_Imputation_;
run;
/*下面mianalyze過程對5次填補結(jié)果進行綜合，給出最終的綜合估計結(jié)果*/
proc?mianalyze?data=outreg;
model effectsIntercept course ht wt;
run;
