手把手教你多重填补法处理缺失数据! 其实R语言是你想不到的简单。。。
目前缺失数据的处理方法大致分为三类：删除法、填补法和不处理。多数统计软件进行缺失数据分析时，默认采用列表删除法。这种简单粗暴的方法在少数变量值缺失时影响不大，但涉及到大范围变量值缺失，实在是伤敌1000，自伤800。
填补的方法包括单一填补和多重填补。单一填补方法一般根据单个变量的特征进行各自填补，但是可能会导致目标变量的扭曲，并且无法解释其变量的不确定性。而多重填补（multiple imputation）则可以一定程度可以克服这些不足。
多重填补方法很多，近年来，MICE多重填补（Multivariate imputation by chained equations）成为流行方法。相比一般的MI方法，限定条件少，它可以处理多种数据类型，具有灵活的实用性。
没错，我今天介绍MICE方法，并利用R软件实现。R软件实现MICE的常用的包是mice包
MICE的基本过程
MICE主要过程包括：填补、统计分析、汇总
函数mice()首先从一个包含缺失数据的数据框开始，然后返回一个包含多个（默认为5个）完
整数据集的对象。每个完整数据集都是通过对原始数据框中的缺失数据进行插补而生成的。由于插补有随机的成分，因此每个完整数据集都略有不同。
然后，with()函数可依次对每个完整数据集应用统计模型（如线性模型或广义线性模型），
最 后 ，pool()函数将这些单独的分析结果整合为一组结果。最终模型的标准误和p值都将准确地反映出由于缺失值和多重插补而产生的不确定性。
第一步，数据准备
我的数据是主要包括11个变量，分别是性别和10条目的量表，量表每个条目是4等级变量，此外还有个10个条目计算出来的总分,由于原始变量存在着缺失，总分同样缺失严重。样本量是1290。
这里必须 要注意的是R语言关于变量的设置，到底是定量还是等级？
首先我产生一个数据集yiyu01，但是我把性别变成定量，看看一个二分类如果是定量形式进入模型，会是什么结果？
这样做的目的是了解不同变量的设定对R语言的影响。
第二步，数据导入
在R语言载入包，导入数据
library(lattice) #调入函数包
library(MASS)
library(nnet)
library(mice) #前三个包是mice的基础
library(foreign)
yiyu01 <- read.spss("E:/百度云同步盘/R/rhealth/yiyu01.sav",to.data.frame = T)
str(yiyu01)
第三步，多重填补
MICE方法中，先通过蒙特卡洛模拟产生多个数据框，然后通过Gibbs方法完成抽样和填补。每个包含缺失值的变量都默认可通过数据集的其它变量预测得来，于是这些预测方程便可以用来预测数据的 有效值，该过程不断迭代直到所有的缺失值都收敛为止。杜宇每个变量，研究者可以选择预测模式的形式和待待选入的变量。
默认地，预测的均值用来代替连续型变量中的缺失数据，而logistic或者多项logistic回归分别用来替换二分类变量和多分类变量的缺失值。
程序详解
可以通过md.pattern看缺失模式
md.pattern(yiyu01)
mice函数填补函数，mice函数
mp01<-mice(yiyu01, m=5，seed=6666)
mice函数详解
mice函数主要参数有mice(数据集，m=5，meth=，seed=)
m为默认插补数据集数量，5为默认值，meth为插补方法，seed为随机化种子数
methods插补方法很多，最常见的包括"pmm", "logreg", "polyreg", "polr"四种：
Predictive mean matching (continuous data)
Logistic regression imputation (binary data, factor with 2 levels)
Polytomous regression imputation for unordered categorical data (factor>= 2 levels)
Proportional odds model (ordered, >= 2 levels)
这四种方法是mice默认的方法，也就是当你method缺省的时候，系统根据每个变量的特征挑选某一种方法来进行填补。
但是默认法都是些常规的方法，可能会出现：
Error in nnet.default(X, Y, w, mask = mask, size = 0, skip = TRUE, softmax = TRUE, : too many (5148) weights
怎么办？一般加上nnet.MaxNWts = ，比如
mp03<-mice(yiyu03,m=2，nnet.MaxNWts = 10000)
其它包括，比如贝叶斯线性回归（norm）、基于bootstrap的线性回归（norm.boot）、线性回归预测值（norm.predict）、分类回归树（cart）、随机森林（rf）等等。（详情可help(mice)获取信息）。
使用这些插补方法对数据有严格的要求，比如贝叶斯线性回归等前三个模型都需要数据符合numeric格式，而PMM、cart、rf任意格式都行。
seed：随机抽样的种子。数字随意，数字一旦确定，每次抽样结果一样
好了，这样的话，R语言就开始填补了。。。。
数据库有点大，迭代时间真是漫长呀！
.
.
.
.
.
.
.
.
.
然后错了，继续改，继续运行
.
.
.
.
.
.
.
.
.
.
.
OK 终于跑完了。形成了填补库mp01，是list格式
现在可以用一些函数都研究下，填补后的基本情况
summary(mp01) #可以了解数据的情况
检查被填补的数据：mp01$imp$sex
1、2分别代表第一列的填补结果和第二列的填补结果
想看填补好的的原始数据？complete 函数
x1<-complete(mp02,action=01)
x2<-complete(mp02,action=02)
action=1 代表填补得到第一个数据库，action=2，代表填补得到的第二个数据库
我基于mp01产生了x1
第四步，分析
对于填补好的数据，可以利用complete函数单独拿出某一个数据库进行分析，也可以进行合并多个填补数据库分析。当然一般是推荐后者。
合并分析利用的是with函数、pool 函数
比如，研究性别sex与dc009_2011的线性关系
fit=with(mp01, lm(dc009_2011~sex))
xs<-pool(fit)
summary(xs)
当然，我们可以纳入更多的变量进行分析。
fit=with(mp01，lm(dc009_2011~sex+dc011_2011 +dc010_2011))
xs<-pool(fit)
summary(xs)
如果要计算置信区间，加个参数即可
summary(xs, conf.int = TRUE,conf.level = 0.95)
如果结局是分类数据，可以调动glm函数来进行
fit=with(mp03,glm(sex~total,family = binomial(link = logit)))
xs<-pool(fit)
xs
如果结局是生产数据，可以调动coxph函数来进行
fit=with(mp04,coxph(Surv(time,censor)~age+sex+bui+ch+p+stage+trt))
summary(fit)
xs<-pool(fit)
总的来说，R语言缺点是速度慢了点，下一次我利用stata做个尝试。