|
多重插补我们已经有过介绍：《缺失值的处理：多重插补》、《缺失值分析基础》。在进行缺失值的多重插补时，软件中常默认插补5次（如SPSS和R的mice包），有时我们有会根据缺失值的缺失比例确定缺失值的插补次数，比如某个数据集平均缺失10%，则MI时可以进行10次插补。但插补从来都不是我们的目的，而是为了更好的进行后面的分析。既然多重插补会产生一系列的插补数据集，自然会有一个问题：我们应该用哪个数据集进行后续的分析呢？ 从逻辑上讲，我们需要对每一个填充好的数据集进行分析，然后对这些结果进行合并。在SPSS中直接提供了一些常用的直接分析多重插后补数据集的分析过程：这些过程会单独分析每一个插补集并给出了汇总结果。这些我们在《缺失值的处理：多重插补》有过介绍，但目前不是所有的SPSS分析过程都能实现这种分析。R中有很多程序包都可以实现多重插补，如mice、Hmisc、mi等。我们以常用的mice包为例介绍多重插补及其汇总分析。
多重插补的分析过程的核心步骤有三：首先通过用合理的数据值替换缺失的值，创建多个完整的数据集，然后估计每个数据集的模型参数，最后将多个数据集的参数值合并为一个。mice包通过几个重要的函数[mice、with和pool]来完成这些步骤。mice函数可以生成多重插补的数据集，结果储存在mids的类中；将mids作为一个参数，可用complete函数获取每次插补后的数据集，结果储存在mils的类中；with函数用于多重插补重复分析（对每次插补数据进行分析），结果储存在mira的类中；pool函数可以合并多次重复分析的结果。示例：来自R语言VIM包的diabetes数据集（Indian Prime Diabetes Data）。该数据集最初来自美国国家糖尿病、消化和肾脏疾病研究所，用来根据数据集中的某些诊断测量值来诊断性的预测患者是否患有糖尿病，所有患者都是至少21岁的皮马印度血统的女性。 使用mice进行多重插补并不难，包括插补后的分析主要有三个核心命令： data(diabetes,package="VIM") #载入数据 impt<-mice(diabetes,m = 5,print = FALSE) #数据插补，默认插补5次fit<-with(impt, glm(Outcome~Age+BMI+BloodPressure+Glucose+DiabetesPedigreeFunction+Insulin+SkinThickness+Pregnancies, family = binomial)) #对插补后的每个数据集分别进行分析summary(est) 但一个完整的插补分析过程不止于此，它应该包含更多的分析内容。【1】缺失值分析 data(diabetes,package="VIM") 数据基本信息包括最小值、第一分位数、中位值、均值、第三分位数、最大值及缺失值数量。各变量的缺失值数量也用以下简单命令查看：miss.nu<-mice(diabetes, maxit=0)md.pattern(diabetes,rotate.names=TRUE) #检查数据缺失模式（missing data pattern） 缺失类型分析显示共有763个缺失值，其中胰岛素（Insulin）缺失374例，皮肤厚度（SkinThickness）缺失 227例，怀孕次数（Pregnancies）缺失111例，血压（BloodPressure）缺失35例，体质指数（BMI）缺失11例，葡萄糖（Glucose）缺失5例。另外有336条记录是无缺失的，在缺失病例中，胰岛素和皮肤厚度共同缺失的病例最多，有170条，其次是单变量胰岛素缺失有119例。 为了明确数据的缺失类型，我们可能还会想知道变量A的数据缺失是否取决于变量B，可以通过制作A变量无缺失和缺失对应的B变量柱状图来验证一下。比如看下缺失比较多的皮肤厚度（SkinThickness）和胰岛素（Insulin）。 library(lattice) histogram(~Insulin|is.na(SkinThickness), data=diabetes) #绘制胰岛素的柱状图命令：histogram(diabetes$Insulin)或histogram(~ Insulin, data=diabetes)很显然，如果Insulin是完全随机缺失，则不论SkinThickness缺失与否，其分布是相同的，本例Insulin缺失取决于SkinThickness，因此Insulin缺失不可能是完全随机缺失。【2】多重插补 mice(data,m=5,method=NULL,predictorMatrix,ignore=NULL,where=NULL,blocks,visitSequence=NULL,formulas,blots=NULL,post=NULL,defaultMethod=c("pmm","logreg","polyreg","polr"),maxit=5,printFlag=TRUE,seed=NA,data.init=NULL ,...)m：指定插补次数，默认是5次；method：指定数据的插补方法，默认方法由参数defaultMethod指定，defaultMethod指定的方法有4个，分别对应连续型变量、二分类变量、无序多分类和有序多分类变量，可用方法可通过命令methods(mice)来查看，具体包括：impt<-mice(diabetes,m=5,print=FALSE) 结果包括两部分内容：插补方法和预测矩阵。 1）插补方法Imputation methods 我们可以看到本例中各变量缺失值插补方法（Imputation methods）都是PMM，对无缺失的列结果显示的是" "。命令impt$meth可以单独查看每一列的缺失值插补方法。当然你可能从专业或者经验上认为，某个变量更合适另外一种插补方法，只需要单独对该变量的插补方法进行修改就可以了，如将Insulin的缺失插补方法改为Bayesian线性回归插补，则命令如下：impt2<-mice(diabetes,seed=20221114,print = FALSE) impt2$meth["Insulin"]<-"norm"然后重新运行插补命令即可： impt2<-mice(diabetes, method=impt2$meth, seed=20221114,print=F)2）预测矩阵PredictorMatrix 命令impt$pred可以单独查看预测矩阵。通过预测矩阵（PredictorMatrix）可以了解缺失变量插补时用到了哪些变量，每个缺失变量在预测矩阵中都有一行和一列，值1表示该列变量用于插补行变量，变量不允许对自身进行插补所以矩阵对角线为零，没有在行变量中显示的变量表示该变量无缺失值，本例DiabetesPedigreeFunction、Age和Outcomg无缺失。我们可能不想用所用的变量进行插补，而是只用/不用某一些，或者只用相关性比较高的一些变量进行缺失值插补，那就需要对这个预测矩阵进行修改。比如我不想用怀孕次数（Pregnancies）来预测其他缺失变量：impt3<-mice(diabetes,seed=20221114,print = FALSE) impt3$pred[ ,"Pregnancies"]<- 0然后用新的预测矩阵重新进行插补就可以啦： impt3<-mice(diabetes, predictorMatrix=impt3$pred,seed=20221114,print=F)还有人想能不能只用某些相关变量进行插补呢？当然！这里就用到了mice包里的quickpred函数： quickpred(data,mincor = 0.1,minpuc = 0,include="",exclude = "",method = "pearson")data：含缺失值的矩阵或数据框；mincor：向量或矩阵，指定最小的阈值，利用被预测和预测变量的相关系数绝对值大于等于该值的的预测变量进行插补，默认为0.1；minpuc：指定可用个案数的最小比例；include：指定必须包含的变量；exclude：指定排除的变量；method指定相关的类型，有pearson、kendall和spearman。impt4<-mice(diabetes, pred=quickpred(diabetes, mincor=0.3), seed=20221114,print=F) #利用绝对相关系数≥0.3的变量进行预测；predictorMatrix可简写为predmice包中还提供了流入和流出（Influx and outflux）统计量，可以用来发现缺失值插补时重要的预测变量，用这两个统计量可以绘制通量图（flux plot）。 Influx and outflux的理论可参见Flexible Imputation of Missing Data. Second Edition.Chapter4.比较复杂，这里只简单说一下这两个统计量的意义。这两个变量取决于变量的缺失值比例，当变量无缺失时Influx取值为0，变量完全缺失其值为1，而outflux则是无缺失其值为1，完全缺失为0。高Influx值的因素更容易被其他变量预测，而高outflux值的因素常作为潜在的预测因子。通量图可用于发现潜在的预测因子、更容易被预测的变量以及干扰插补模型的变量。一般来说，离对角线较近的变量通常比距离较远的变量关联性更好（关联性好的预测性更好），位于较低区域（尤其是左下角附近）且对后续分析不感兴趣的变量最好在插补之前从数据中删除。本例SkinThickness和Insulin的outflux值较低，丢失数据比较严重，预测能力有限，Influx值也较低，被预测的效果可能也一般。另外，通过改变插补方法，结合插补次数还可以实现均值填充、回归填充、随机回归填充、随机森林填充等。 impt5<-mice(diabetes, method ="mean", m=1, maxit=1, seed=20221114) impt6<-mice(diabetes, method="norm.predict", m=1, maxit=1, seed=20221114) impt7<-mice(diabetes, method ="norm.nob", m=1, maxit=1, seed=20221114) impt8<-mice(diabetes, method ="rf", m=1, maxit=10, seed=20221114) 数据插补完成，你可能想看一下插补的数据啥样： impt$imp #查看所有缺失变量的插补数据（注意这个命令仅显示缺失值的插补结果）impt$imp$BMI #查看BMI的插补数据（仅显示缺失值的插补结果）impt$data #查看原始数据 
【3】插补结果检查 impt<-mice(diabetes,maxit=30, seed=20221113,print = FALSE)结果显示的是插补值得均值（左）和标准差（右），如果线迹充分混合且没有趋势，表明算法收敛良好，一般默认的5次迭代就可以达到很好的效果。进一步的诊断检查： stripplot(impt) #单独绘制血压的分类散点图：stripplot(impt, BloodPressure~.imp, pch=8, cex=2)结果中蓝色绘制的是观察数据，红色绘制插补数据。如果数据是完全随机缺失（MCAR），则插补数据应具有与观测数据具有相同的分布。如果分布不同，缺失值可能是是MAR（甚至是MNAR）。或者使用箱线图： 密度图也常用来查看插补数据集与观察数据的分布情况： densityplot(impt) #单独绘制血压的分类散点图可用densityplot(impt,~Insulin)相对而言，插补结果与原数据的分布类型比较相似，Glucose和BMI偏离相对大一些。【4】提取插补数据集 提取某一次的插补数据集，这样可以单独对某一次的插补结果进行分析。complete(data, action = 1L, include = FALSE, mild = FALSE, ...)data：函数mice()产生的mids类对象；action：小于等于插补次数，如1表示提取第1次的插补数据集，0表示原始数据（未进行插补的含缺失值的数据），默认为首次插补数据集。也可以是字符"all"、"long"、"stacked"、"broad" 或 "repeated"，long和stacked都是堆叠型数据格式，两者的区别在于stacked结果不显示插补次数和个案的所在行数，而broad和repeated都是每次插补结果占一列的格式，两者的区别在于列顺序不一样。include：是否包含原始数据集。complete(impt,3) #提取第3次插补结果complete(impt,"long") #提取插补结果，格式为“long”【5】插补数据集分析 fit<-with(impt, glm(Outcome~Age+BMI+BloodPressure+Glucose+DiabetesPedigreeFunction+Insulin+SkinThickness+Pregnancies, family = binomial)) #对插补后的每个数据集都进行模型拟合fit #显示每个变量的缺失值数量，以及每个插补数据集的logistic回归参数估计值summary(fit$analyses[[3]]) #显示第3个插补数据集的logistic回归参数估计值及统计学检验结果我们会发现每个插补数据集的分析结果都不太一样，那应该以哪次结果为准呢？这时就需要用到mice包中pool()或pool.syn()函数来进行结果汇总了。【6】汇总分析 m：插补次数；ubar：插补组内方差的均值；b：插补间方差；t：合并估计的总方差（笔者觉得应该是标准误的平方）；dfcom：插补数据集的自由度；df：假设检验的残余自由度；riv：因缺失导致的方差相对增加；lamba：由于缺失而导致的总方差的比例；fmi：因缺失导致的信息（系数）丢失比例。我们可以看到Insulin和Glucose的lambda和fmi值偏大，Insulin可能是因为缺失值太多导致，而Glucose偏大原因未知，不过从前面箱线图的插补结果来看，该变量第二次插补结果值有些偏高的。 summary(est) 理论上，只要函数符合R中的coeff方法的准则规范，用其进行多次插补后的结果都可以用pool函数进行汇总，而且pool.scalar函数可以“手动”汇总m个重复完整数据的单变量分析估计值。 pool.scalar(Q,U,n=Inf, k=1,rule=c("rubin1987","reiter2003"))Q:单次估计值向量；U：单次估计值的方差向量：n：样本量；k：需要估计的参数数量；rule：汇总规则。summary(fit$analyses[[1]]) #获取第1次插补结果的参数估计值 summary(fit$analyses[[2]]) #获取第2次插补结果的参数估计值summary(fit$analyses[[3]]) #获取第3次插补结果的参数估计值summary(fit$analyses[[4]]) #获取第4次插补结果的参数估计值summary(fit$analyses[[5]]) #获取第5次插补结果的参数估计值以合并BMI估计值为例： bmi.pool<-pool.scalar(Q=c(0.0906097,0.0819811,0.083386,0.0925563,0.0847294), U=c(0.0199806^2,0.0194731^2,0.019599^2,0.0197666^2,0.0199036^2),n=768)结果跟前面直接汇总的结果参数一致，除了m、ubar、b、t、dfcom、df、riv、lamba、fmi外，多了qhat和u两个指标，其中qhat即输入的Q参数，是每次插补的参数估计结果，u即输入的U参数，说明文件上指的是每次插补的参数估计的方差，但从后面的P值计算时推断这个应该是均数之差的标准误。 我们感兴趣的值汇总估计值及标准误： 结果显示，手动合并结果同与直接进行summary(pool(fit))的结果完全一致。这样我们就可以对任意的单变量分析结果进行汇总分析啦。比较遗憾的是手动合并的结果中并未显示合并统计量的P值，但可以计算： statistic=bmi.pool$qbar/sqrt(bmi.pool$t) P<-2*pnorm(statistic,lower.tail=FALSE)还有系数的95%的置信区间，采用正态渐进法如下： bmi.pool$qbar-qnorm(0.05/2,lower.tail=FALSE)*sqrt(bmi.pool$t) #95%CI下限bmi.pool$qbar+qnorm(0.05/2,lower.tail=FALSE)*sqrt(bmi.pool$t) #95%CI上限不过一般软件中用计算置信区间的默认方法好像不是该法，所以具体数字会略有差异。 在汇总分析时还有最后一个重要的问题：如何解决变量选择的问题，比如逐步回归、lasso回归等。 多重插补包括三个重要的步骤：插补生成多个数据集，分析每个数据集，汇总每个数据集的估计结果。如果在分析每个数据集时用到了变量选择，那么就可能碰到这样一种情况：不同的数据集选择的变量不一致。如是，接下来的汇总分析就没法进行了。这个问题至少有三种解决策略可以考虑：（1）多数服从少数（Majority）法。单个数据集分析时，利用一半以上的模型都选择的变量进行建立；（2）堆叠加权（Stack）法。将m个插补的数据集合并成一个堆叠数据集，每条记录赋予一个固定的权重1/m;（3）Wald统计量法。这一部分的具体实操以后另起一文吧。 — END —
|
