【原】缺失值的插补：随机森林法

随机森林(Random Forest，RF)是机器学习（Machine Learning）领域中非常有名的集成学习（Ensemble Learning）算法，相当于套袋法（Bagging）和决策树的组合。简单地说，就是采用套袋法进行抽样来构建多棵决策树（分类树或回归树），然后采用投票或取均值的方式做出最终决策。套袋法的随机抽样采用的是自举法（Bootsrap），Bootsrap是从给定训练集中做有放回的均匀抽样，可以理解为每抽到一个样本后，都会将这个样本放回，之前采抽到的样本在放回后有可能继续被抽到。一般会随机抽取和训练集样本数相同的样本，这样得到的抽样集和训练集样本的个数相同，但是每个数据集的具体样本并不同。每个训练样本可以构建一棵决策树，多棵决策树就构成了森林，利用多棵树的预测结果投票选出最终分类或计算最终均值。随机森林就是套袋法的改进，与套袋法不同的是在构建决策树分枝的时候，RF只考虑部分预测变量，而Bagging是将所有变量都考虑进去。

因为是有回放的随机抽样，某个样本在Bagging的m次采样中都没有被采到的概率是36.8%（数学好的可以算一下），也就是说训练集中大约有36.8%的数据没有被采样，这大约36.8%的没有被采样到的教据就是袋外数据(Out Of Bag, 简称OOB)。这些数据没有参与训练集模型的拟合，可以取代测试集来估计误差，可用于模型的验证，检测模型的泛化能力。

missForest是一种利用随机森林来插补缺失值的非参数方法，适用于连续性变量，也适用于分类变量，其核心算法是利用已知的变量作为自变量，将含有缺失值的变量作为因变量建立随机森林来预测缺失值。

假定存在数据集X，对任意一个含有缺失值的变量Xs可以将数据集分成4个部分：变量Xs观测值（y_s.obs）、变量Xs的缺失值（y_s.mis）、变量Xs观测值对应的非Xs变量值（x_s.obs）、变量Xs缺失值对应的非Xs变量值（x_s.mis）。missForest是如何对Xs中的缺失值进行填充的呢？missForest会先对数据集中的缺失值用相应变量的中位数/均值、众数进行初始填充；然后按照缺失值的数量从低到高对变量进行排序。之后利用y_s.obs和x_s.obs作为训练集构建随机森林模型，训练好的模型就可以通过x_s.mis来预测缺失值y_s.mis了。当不止一个变量存在缺失值时，就需要不断重复迭代计算直到满足收敛标准。停止准则定义为新插补的数据矩阵与前一个数据矩阵之间的差异在两种类型的变量中出现首次增加。

缺失值插补后，还需要衡量插补值与真实值之间的偏差，这种插补效果的评估一般采用标准均方根误差（normalized root mean squared error，NRMSE）和缺失值错分比例（proportion of falsely classifified，PFC) 。NRMSE用于连续性缺失变量，PFC用于分类缺失变量。

示例：来自R语言VIM包的diabetes数据集（Indian Prime Diabetes Data）。该数据集最初来自美国国家糖尿病、消化和肾脏疾病研究所，用来根据数据集中的某些诊断测量值来诊断性的预测患者是否患有糖尿病，所有患者都是至少21岁的皮马印度血统的女性。

由于原数据集解释变量全部为连续变量，笔者按界值25、30、35、40将BMI转换成因子变量BMIc：正常(0)、超重(1)、一级肥胖(2)、二级肥胖(3)、三级肥胖(4)。

data(diabetes,package="VIM") #VIM是一个强大的缺失值处理程序包，是Visualization and Imputation of Missing Values的简称，里面含有各种缺失值处理方法及可视化处理函数，以后有机会再介绍这个函数

diabetes$BMIc<-ifelse (diabetes$BMI<24,0,diabetes$BMI)

diabetes$BMIc<-ifelse (diabetes$BMI>=24 & diabetes$BMI<30,1,diabetes$BMIc)

diabetes$BMIc<-ifelse (diabetes$BMI>=30 & diabetes$BMI<35,2,diabetes$BMIc)

diabetes$BMIc<-ifelse (diabetes$BMI>=35 & diabetes$BMI<40,3,diabetes$BMIc)

diabetes$BMIc<-ifelse (diabetes$BMI>=40,4,diabetes$BMIc)

diabetes$BMIc<-factor(diabetes$BMIc)

diabetes<-as.data.frame(diabetes)

missForest(xmis, maxiter=10, ntree=100, variablewise=FALSE, decreasing=FALSE, verbose=FALSE, mtry=floor(sqrt(ncol(xmis))), replace=TRUE, classwt=NULL, cutoff=NULL, strata=NULL, sampsize=NULL, nodesize=NULL,  maxnodes=NULL, xtrue=NA, parallelize=c('no', 'variables', 'forests'))

几个比较重要的参数：xmis指定需要填充的矩阵；maxiter在预先给定的停止准则未能满足时的最大迭代次数；ntree指定每个森林中成长树的数量；variablewise为每个变量分别返回OOB误差；decreasing设置填充变量的排列顺序；verbose提供迭代过程的额外输出，包括估计的插补误差、运行时间、真实的插补误差（需要跟参数xtrue指定完整的数据）；每次树分支成长时随机抽的变量数量；replace默认为bootstrap有回放抽样，FALSE为子抽样（无回放抽样）；classwt、cutoff、strata、sampsize用于非平衡数据、分层抽样及集中选择；nodesize和maxnodes用于设置终端节点的最大数量和终端节点的最小观测数；xtrue可以指定没有缺失的数据矩阵，用于估计真实的插补误差。

插补后可通过ximp、OOBerror和error 来查看插补后的数据、估算的OOB插补误差（连续性变量是NRMSE，分类变量是PFC）和真实的插补误差（xtrue参数提供时有效）。

在随机森林程序包randomForest中有一个内置的函数rfImpute也可以完成随机森林的缺失值插补，但具体算法上略微不同，rfImpute{randomForest}是通过对随机随机森林中的渐进值进行加权来完成的。

library(missForest)

set.seed(20220801)

dbt.impt<-missForest(diabetes) #数据插补过程其实就上面这一行代码，大部分参数按默认处理即可

dbt.impt

结果最后显示的NRMSE和PFC分别表示连续性变量和分类变量袋外数据误差。如果想查看每一次迭代的袋外误差，可修改命令参数如下：

missForest(diabetes, verbose=TRUE)

如果你想把插补后的完整数据保存下来，或者将数据导出进行后续的操作或者在其他软件中进行分析，下面的命令可以提供帮助：

diabetes_impt<-dbt.impt$ximp

write.csv(nmes,file = "D:/Temp/diabetes_impt.csv") #数据导出后可在excel中对数据进行清理

缺失值的插补从来都不是我们的最终分析目的，插补只是前期的数据处理，其后才是数据分析的开始。我们在《缺失值的处理：多重插补》一文中，曾构建了一个缺失值数据集，缺失类型囊括了完全随机缺失、随机缺失和非随机缺失。下面我们利用随机森林法对这个这个数据库的缺失值进下插补。

##导入构建的含缺失值数据集

library(haven)

mheart_1<- read_sav("D:/Temp/mheart_1.sav")

mheart_1<-mheart_1[,c(-1,-10)]

mheart_1$Attack<-factor(mheart_1$Attack)

mheart_1$Female<-factor(mheart_1$Female)

mheart_1$MarStatus<-factor(mheart_1$MarStatus)

mheart_1$HSGrad<-factor(mheart_1$HSGrad)

mheart_1$Alcohol<-factor(mheart_1$Alcohol)

mheart_1$Smoker<-factor(mheart_1$Smoker)

mheart_1<-as.data.frame(mheart_1)

##导入原始的完整数据作为参照

mheart_0<- read_sav("D:/Temp/mheart_0.sav")

mheart_0<-mheart_0[,c(-1,-10)]

mheart_0$Attack<-factor(mheart_0$Attack)

mheart_0$Female<-factor(mheart_0$Female)

mheart_0$MarStatus<-factor(mheart_0$MarStatus)

mheart_0$HSGrad<-factor(mheart_0$HSGrad)

mheart_0$Alcohol<-factor(mheart_0$Alcohol)

mheart_0$Smoker<-factor(mheart_0$Smoker)

mheart_0<-as.data.frame(mheart_0)

##对参照数据（未删除变量值的原始数据）进行logistic回归

lgrfit0<-glm(Attack~Female+Age+BMI+MarStatus+HSGrad+Alcohol+Smoker, family=binomial(link = logit), data=mheart_0)

summary(lgrfit0)

exp(cbind(OR=coef(lgrfit0),confint(lgrfit0)))

##直接删除缺失值个案，对数据完整的个案进行logistic回归

mheart_2<-na.omit(mheart_1)

lgrfit2<-glm(Attack~Female+Age+BMI+MarStatus+HSGrad+Alcohol+Smoker, family=binomial(link = logit), data=mheart_2)

summary(lgrfit2)

exp(cbind(OR=coef(lgrfit2),confint(lgrfit2)))

与参照数据的回归结果相比，暴力删除缺失会使很多变量的系数发生了较大变化。

##采用随机森林法对缺失值进行填插补，并对插补后的数据进行logistic回归

library(missForest)

set.seed(20220801)

mheart.imperror<- missForest(mheart_1, xtrue=mheart_0)

mheart_impt<-mheart.imperror$ximp #插补后的数据集

lgrfit1<-glm(Attack~Female+Age+BMI+MarStatus+HSGrad+Alcohol+Smoker, family=binomial(link = logit), data=mheart_impt)

summary(lgrfit1)

exp(cbind(OR=coef(lgrfit1),confint(lgrfit1)))

与原始的完整数据（参照）相比，有些变量的插补表现并不是太好，比如婚姻状况。在我们构建缺失值时，Age和Alcohol属于完全随机缺失（MCAR），BMI属于随机缺失（MAR），而MarStatus属于非随机缺失（MANR），对非随机缺失现有的插补方法都表现一般。当然真实世界中我们是没法知道数据缺失前本来的面貌的。不过从统计学意义上看，结果具有一致性。注意本示例为虚拟数据，不要太在意不合常理的专业结论。

感兴趣的也可以翻看《缺失值的处理：多重插补》，与多重插补结果进行一下比较。