有缺失值怎么办——多重填补技术的思路与软件实现

多重填补是针对单一填补而言，所谓单一填补，就是对每一缺失值填一个值，这类方法不能反映缺失数据的不确定性，因而容易导致标准误低估。比如上一篇文章说的回归方法填补，填补后，由于填补值是根据直线回归估计出来的，他们都在一条直线上，没有误差，标准误肯定低估了，因为正常情况下，所有的数据距离回归线总是有点误差的。

比如下图，可以看出，估计的值，都在直线，误差肯定小。

而多重填补是模拟生成一个缺失数据的随机分布，然后从这一分布中随机抽取数据作为缺失值的填补，因而可以反映不确定性。正因为是随机抽取，导致多重填补技术存在一个副作用：用该法填补数据会同时产生多组填补值，而且你的填补值跟我的填补值永远都不一样（除非指定同样的种子数），换句话说，如果你发表一篇填补数据的文章，即使我用你的程序重新跑一遍，结果也跟你的不一样（你只能祈祷结果差别不要太大）。

下面介绍一下MI法的思想。

设有缺失值的变量为y，用于预测填补的变量为x（简单起见，设只有一个变量），利用回归填补法，可以得到：

y=a+bx

根据x值可以填补y值。但这样填补后的值，如果再分析y与x的关系，会使二者关联人为增大，标准误变低（原因前面刚说了）。

为了解决这一问题，很自然的一个想法就是，把y加上一个随机波动项，使y和x之间不是严格的直线关系。即把二者关系变成：

y=a+bx+e

e是一个随机数字（通常可以从残差中随机抽取）。这样的话，y和x之间就不是严格的线性关系，填补的值就不在y=a+bx的直线上了，而是有一定的偏离（偏离大小为e）。通过这种思路，可以使得标准误变大一些，更接近实际。（事实上，大家仔细想想线性回归，也就是这么回事，回归线与数据点总是有点差距的，这就是残差）

上述方式尽管会使标准误略有增大，但这仍不够，因为它们都是在a+bx的基础上添加一个随机波动。也就是说，是把a和b这两个系数当做真实参数，而实际上它们只是样本估计值。

我们可以想象一下实际情况：总体的截距和斜率是无法知道的，只能根据当前的一次抽样数据获得样本截距a和斜率b。理论上，如果重复10次抽样，应该会得到10个a和b，每次计算的a和b肯定会有所不同，它们都只是围绕总体截距和斜率波动的随机值。

因此，为了更加符合实际情况，还需要把a和b也作为随机波动。那怎么把a和b作为随机抽取的值呢？这就需要有一个关于a和b的分布，然后从这个分布中随机抽取a和b。这一分布通常称为贝叶斯后验分布，也就是说，我们要从贝叶斯后验分布中随机抽取a和b，每次抽取的都会有所不同，这样建立的模型也会不同，从而估计的缺失值也是随机的。这相当于用统计的方法来模拟实际，让得到的结果最大可能地接近实际情形。

现在就面临这样的问题：如何产生贝叶斯后验分布并从该分布中随机抽取？目前大多数统计软件都是用马尔科夫链蒙特卡洛（Markov chain Monte Carlo，MCMC）模拟的方法。所谓MCMC，就是两种技术的组合，蒙特卡洛主要是用来模拟抽样，马尔科夫链主要是用来产生平稳的分布链。

比如下图就是一个MCMC模拟产生的一个分布，缺失值的填补就是从产生的分布中随机抽取，而且通常是随机抽取多个，所以是多重填补。

可能有人会说，我得到多个填补值，怎么用到最终的分析呢？一般是这样：比如取5个，你可以计算其均值，也可以用这5次填补的结果重新做统计分析（如回归），最后得到一个综合的结果。如下图：

MCMC的算法极其复杂（幸亏我们有统计软件），但思路还是比较容易理解的，主要步骤如下：

第一，先利用无缺失的数据作为初始的均数和协方差矩阵，建立一个y对x的回归方程，得到回归系数。

第二，将y中的缺失值利用回归方程填补上，并对每一填补值加上一个从残差中随机抽取的数值。

第三，利用填补的“完整”数据，重新计算一个新的均数和协方差矩阵。

第四，根据第三步中得到的新的均数和协方差矩阵的后验分布，从中随机抽取一个均数和协方差。

第五，利用第四步中抽取的均数和协方差，再回到第一步，把刚才步骤再次执行一遍。如此多次循环，直到最终的收敛（通常收敛的意思是指，上一次的结果与本次结果差异非常小甚至无差异）。用最后收敛所得到的填补值作为最终填补值。

所以，其实多重填补技术也并不是难以捉摸，只要仔细想想，其实思路还是挺清楚的。当然，可能有的人更关心的是到底如何在软件上实现。下面就给出SAS和R的实现语句，感兴趣的自己练习一下。

SAS可通过proc mi和procmianalyze实现EM算法填补、多重填补等：

/* 以下程序产生EM算法的填补结果*/

proc mi data=aa nimpute=0;

em itprint out=emimp;

varycourse ht wt; /* 要填补的变量*/

run;

proc print data=emimp;

run;

/* 以下程序产生多重填补结果*/

proc mi data=aa nimpute=5 out=outmi; /*指定产生5次填补结果*/

mcmc displayinit initial=em(itprint) plots=trace;

/* initial=em(itprint)表示以EM算法作为初始估计值，并输出迭代结果*/

vary course ht wt;

run;

/*下面reg过程分别根据5次填补值输出回归结果，为后面的mianalyze过程准备*/

proc reg data=outmi outest=outreg covout;

mode ly=course ht wt;

by _Imputation_;

run;

/*下面mianalyze过程对5次填补结果进行综合，给出最终的综合估计结果*/

proc mianalyze data=outreg;

model effectsIntercept course ht wt;

run;

R软件可通过mice包执行多重填补，设含缺失值的数据集为f1，则主要语句包括：

library(mice)

im=mice(f1,seed=123,m=10) #指定输出10套填补值#

fit=with(im,lm(y~course+ht+wt))

ps=pool(fit) #10套填补结果的合并#

im$im$y #输出y的填补结果，其它以此类推#

complete(im,action=2) #输出第2套填补值，其它以此类推#

如果觉得文章有用，不用赞赏，花一秒钟帮忙点击下面的广告就可以了（点开再关闭就行）。谢谢！