较小的x,m 表示“观察到”，而较大的值表示随机量。大家正在讨论的缺失机制是（X*，M）:的关系或联合分布的假设

随机完全缺失(MCAR) : 值缺失的概率就像抛硬币，与数据集中的任何变量无关。缺失值只是一个麻烦。您可以忽略它们，只关注数据集中完全观察到的部分，这样就不会有偏差。在数学中，对于所有m 和x:

随机缺失(MAR) : 缺失的概率现在可以取决于数据集中观察到的变量。一个典型的例子是收入和年龄等两个变量，其中总是观察到年龄，但对于某些年龄值可能会缺少收入。这听起来可能很合理，但事情可能会变得复杂。在数学中，对于所有m 和x:

非随机缺失(MNAR) : 这里一切皆有可能，我们不能一概而论。但最终我们需要学习在一种模式m’ 中给出的观测值的缺失值的条件分布，以便在另一种模式m’ 中进行插补。

实现此目的的一种众所周知的方法称为“链式多重插补”(MICE)：首先使用简单的插补方法（例如均值插补）填充值。然后，对于每次迭代t，对于每个变量j，对所有其他估算变量（已估算）执行回归分析。然后将这些变量的值填充到所有未观察到的X_j 的学习插值器中。在R语言中，可以使用mice包轻松实现。我的方法在实际应用中效果非常好，并且MICE 重现某些实例的底层分布的能力令人惊叹。下面我们将使用一个非常简单的示例（仅缺少一个变量，因此我们可以手动编码）来模拟MICE 通常迭代执行的过程，以更好地理解其工作原理。

数据

我们将看一个非常简单但说明性的示例：考虑具有两个联合正态变量X_1、X_2 的数据集。我们假设两个变量的方差为1，正相关系数为0.5。为了便于说明，我们假设X_1是收入的对数，X_2是年龄。另外，我们在设定收入X_1时存在一个缺失机制，即当年龄“高”时，X_1往往会缺失。这是我们设置的：

当X_2（年龄）较大（即大于0）时，X_1（收入）缺失的概率为0.8。既然我们假设.虽然这个例子很简单，但如果我们假设老年人收入较高，那么从一种模型切换到另一种模型时，收入和年龄的分布会发生显着变化。在模型m2 中，缺少收入，观察到的年龄和（未观察到的）收入的值往往较高。

我们先看一下生成数据的代码：

库(MASS)库(小鼠)set.seed(10)n-3000Xstar – mvrnorm(n=n, mu=c(0,0), Sigma=matrix( c(1,0.7,0.7,1), nrow=2, byrow=T ))colnames(Xstar) – Paste0(‘X’,1:2)## 引入缺失机制M-matrix(0, ncol=ncol(Xstar), nrow=nrow(Xstar))M[Xstar[,2 ] 0, 1]- sample(c(0,1), size=sum(Xstar[,2] 0), Replace=T, prob=c(1-0.8,0.8) )## 这会产生观察到的通过使用M:X-XstarX[M==1] 屏蔽X^* 的数据集- NA## 绘制分布shiftpar(mfrow=c(2,1))plot(Xstar[!is.na(X[,1]), 1:2], xlab=”, main=”, ylab=”, cex=0.8, col=’深蓝’, xlim=c(-4,4), ylim=c(-3,3))plot( Xstar[is.na(X[,1]),1:2], xlab=”, main=”, ylab=”, cex=0.8, col=’深蓝’, xlim=c(-4,4) , ylim=c(-3,3))

插补是一个分布预测问题

如何处理缺失值引起的分布变化？发现填充缺失值的最佳插补算法

研究该图立即表明均值和回归插值效果不好，因为它们完全无法再现原始数据分布。相比之下，高斯插值看起来相当不错，事实上我认为很难将它与真实的东西区分开来。这可能看起来只是一个技术概念，但它会产生后果。想象一下，如果给您任何这些插补数据集，现在您想要在X_1 上回归X_2 时找到回归系数（与我们用于插补的相反）。在这种情况下，真实值由beta=cov(X_1, X_2)/var(X_1)=0.7 给出。

## 将X_2 回归到=data.frame(impnormpredict))$coefficients[‘X1’],2)## beta=0.90## 高斯插补估计(lm(X2~X1, data=data.frame(impnorm))$系数[‘X1’],2)## beta=0.71## 真值插补估计round(lm(X2~X1, data=data.frame(Xstar))$coefficients[‘X1’],2)## beta=0.71高斯插值的结果非常接近0.7 (0.71)，更重要的是，非常接近使用完整（未观察到的）数据获得的估计！均值插值法低估了beta 值，而回归插值法则高估了beta 值。回归插补因为条件均值插补人为地强化了变量之间的关系，这会导致科学和（数据科学）实践中的估计效果被高估！

回归插补可能看起来过于简单，但这正是机器学习和其他领域常用的插补方法的工作原理。例如，knn插值和随机森林插值（即missForest）。特别是随机森林插值在多篇基准论文中得到了赞扬和推荐，并被广泛使用。 missForest 是在观察到的数据上拟合随机森林，然后通过条件均值进行简单插值。使用它的结果将与回归插值非常相似，导致估计中变量和偏差之间的关系被人为强化！

如何评估插补方法？

上面我们说了插值应该看成是一个分布预测问题，那么这个分布预测问题应该如何评价呢？

假设我们开发了一种新的插值方法，现在想要将其与现有方法（例如missForest、MICE 或GAIN）进行基准测试。在这种情况下，我们人为地引入缺失值。然后将该真实数据集与我们的插补结果进行比较。我们假设上面的回归插值是一种新方法，我们想将其与均值插值和高斯插值进行比较。

即使在最负盛名的会议中，也是通过计算均方根误差(RMSE) 来完成的：

但是使用RMSE 来评估我们的插补会偏向那些插补条件均值的方法，例如回归插值、knn 插值和missForest。

而就像我们要解释的问题一样，由于插值本质上是一个分布预测问题，除了使用RMSE等统计指标外，还应该评估插值方法是否能够恢复数据的整体分布。这可以通过比较不同插值方法生成的数据分布的统计特性（例如均值、方差、偏度等）或使用更复杂的分布相似性度量（例如Earth Mover 距离或Kullback-Leibler 散度）来实现。这里我们使用能量距离来衡量完全观测数据的分布与估算的“分布”之间的差异

能量距离：能量距离是衡量两个分布之间差异的统计量，基于从两个分布中随机抽取的样本对之间的距离。在插补的背景下，它有效地评估了插补数据分布与原始数据分布之间的相似性。

通过使用能量距离等分布敏感的评估指标，我们可以更准确地评估插值方法是否能够恢复数据的原始统计特性，而不仅仅是平均水平上的接近程度。这种方法特别适合评估旨在捕获数据的复杂依赖性和分布特征的先进插值技术。

如何处理缺失值引起的分布变化？发现填充缺失值的最佳插补算法

library(energy)## 计算能量距离的函数：# impX 是插补数据集# Xstar 是完全观测的数据集## 使用energy 包的eqdist.e 函数计算能量距离energycalc – function(impX, Xstar) {# Note: eqdist.e 计算测试的能量统计量，实际上是#=n^2/(2n)*energydistance(impX,Xstar)，但我们只对相对值感兴趣round(eqdist.e( rbind( Xstar,impX), c(nrow(Xstar), nrow(impX)) ),2)} 现在我们将这两个分数应用到这个研究项目中，并测试我们的回归插补是否优于其他两个：

par(mfrow=c(2,2))## 与之前相同的图，但现在带有RMSE 和能量距离## 添加了plot(meanimp[!is.na(X[,1]),c(‘X2′,’ X1′)], main=paste(‘平均插补’, ‘\nRMSE’, RMSEcalc(meanimp, Xstar), ‘\nEnergy’, energycalc(meanimp, Xstar)), cex=0.8, col=’darkblue’, cex .main=1.5)points(meanimp[is.na(X[,1]),c(‘X2′,’X1′)], col=’darkred’, cex=0.8 )plot(impnormpredict[!is.na (X[,1]),c(‘X2′,’X1’)], main=paste(‘回归插补’,’\nRMSE’, RMSEcalc(impnormpredict, Xstar), ‘\nEnergy’, energycalc(impnormpredict, Xstar)), cex=0.8, col=’深蓝’, cex.main=1.5)points(impnormpredict[is.na(X[,1]),c(‘X2′,’X1′)], col=’ darkred’, cex=0.8 )plot(impnorm[!is.na(X[,1]),c(‘X2′,’X1’)], main=paste(‘高斯插补’,’\nRMSE’, RMSEcalc (impnorm, Xstar), ‘\nEnergy’, energycalc(impnorm, Xstar)), col=’darkblue’, cex.main=1.5)points(impnorm[is.na(X[,1]),c(‘X2 ‘,’X1′)], col=’深红色’, cex=0.8 )plot(Xstar[!is.na(X[,1]),c(‘X2′,’X1′)], main=’真相’, col=’深蓝’, cex.main=1.5)points(Xstar[is.na(X[,1]),c(‘X2′,’X1′)], col=’深绿’, cex=0.8 ）

使用均方根误差(RMSE) 进行评估时，回归插值看起来表现非常好！它超过了平均插值和高斯插值。但这与我们上面的分析相冲突，选择回归插补可能会导致结果存在很大偏差。而（缩放的）能量距离正确地将高斯插值识别为最佳方法，与视觉和更好的参数估计一致。

因此，在评估插补方法时（当有真实数据时），应避免使用RMSE 和MAE 等度量。相反，应该将其视为分布预测问题并使用能量距离等分布度量进行评估。过度使用RMSE 作为评估工具对该领域的研究产生了一些严重影响。

当基础观测数据不可用时，评估问题变得更加困难。在本文[1] 中，开发了一个分数，即使在这种情况下也可以对插补方法进行排名！

library(mice)source(‘Iscore.R’)methods-c(‘mean’, #mice-mean’norm.predict’, #mice-sample’norm.nob’) # 高斯插补## 我们首先定义函数允许三种方法的插补：imputationfuncs-list()imputationfuncs[[‘mean’]] – function(X,m){ # 1. 估计平均值meanX-mean(X[!is.na(X[,1]), 1])## 2. Imputemeanimp-Xmeanimp[is.na(X[,1]),1] -meanXres-list()for (l in 1:m){res[[l]] -meanimp}return(res) }imputationfuncs[[‘norm.predict’]] – function(X,m){ # 1. 估计回归lmodelX1X2-lm(X1~. data=as.data.frame(X[!is.na(X[,1] ]),]) )## 2. Imputeimpnormpredict-Ximpnormpredict[is.na(X[,1]),1] -predict(lmodelX1X2, newdata=as.data.frame(X[is.na(X[,1]) ]),]) )res-list()for (l in 1:m){res[[l]] – impnormpredict}return(res)}imputationfuncs[[‘norm.nob’]] – function(X,m){ # 1. 估计回归lmodelX1X2-lm(X1~. data=as.data.frame(X[!is.na(X[,1]),]) )## 2. Imputeimpnorm-Xmeanx-predict(lmodelX1X2, newdata=as.data.frame(X[is.na(X[,1]),]) )var – var(lmodelX1X2$residuals)res-list()for (l in 1:m){impnorm[is.na(X) [,1]),1] -rnorm(n=length(meanx),mean=meanx,sd=sqrt(var))res[[l]] – impnorm}return(res)}scoreslist – Iscores_new(X,imputations=NULL, imputationfuncs=imputationfuncs, N=30) Scores-do.call(cbind,lapply(scoreslist, function(x) x$score ))names(scores)-methodsscores[order(scores)]#meannorm.predictnorm .nob # -0.7455304 -0.5702136 -0.4220387Scores 可以在不看到缺失数据的值的情况下识别分布，特别是当数据具有两个以上维度时。

随机缺失比你想象的更奇怪

在阅读缺失值插补的文献时，人们往往会认为当缺失数据机制为MAR（Missing At Random，随机缺失）时问题就已经解决了，而所有缺失问题都来自于是否可以假设是三月。虽然在使用标准程序（例如最大似然）时这可能是正确的，但如果想要找到一种好的（非参数）插值方法，则情况并非如此。

论文[1]讨论了MAR下可能发生的复杂的分布变化。当从完全观测的模型变为需要插值的模型时，观测变量的分布是否可能发生变化？回到上面的例子，其中X_1是收入，X_2是年龄。正如您在第一张图中所看到的，分布看起来完全不同。以及条件分布X_1 | X_2保持不变！原则上这允许识别正确的估算分布。

这一发现对于理解和处理MAR 下的缺失数据具有重要意义。尽管数据在完全观察和部分缺失时可能会出现不同的分布，但通过关注条件分布的稳定性，可以更准确地估算缺失值。这意味着非参数插值技术可以适应数据的这种固有复杂性，而不是简单地应用参数方法或假设数据具有均匀分布。在实际应用中，这种深入的分析揭示了为什么简单的插补方法（例如均值插补或条件平均插补）可能不足以捕获数据的真实结构，特别是当涉及复杂的相关性和动态分布时。

库(drf)## mouse-DRF ##par(mfrow=c(2,2))#Fit DRFDRF – drf(X=X[!is.na(X[,1]),2, drop=F] , Y=X[!is.na(X[,1]),1, drop=F], num.trees=100)impDRF-X# 预测未观测点的权重wx-predict(DRF, newdata=X[is.na (X[,1]),2, drop=F] )$weightsimpDRF[is.na(X[,1]),1] -apply(wx,1,函数(wxi) 样本(X[!is.na (X[,1]),1, drop=F], size=1, 替换=T, prob=wxi))plot(impDRF[!is.na(X[,1]),c(‘X2’, ‘X1’)], main=paste(‘DRF 插补’, ‘\nRMSE’, RMSEcalc(impDRF, Xstar), ‘\nEnergy’, energycalc(impDRF, Xstar)), cex=0.8, col=’darkblue’, cex.main=1.5)points(impDRF[is.na(X[,1]),c(‘X2′,’X1′)], col=’darkred’, cex=0.8 )## mouse-cart## impcart-Ximpcart[is.na(X[,1]),1] -mice.impute.cart(X[,1], ry=!is.na(X[,1]), X[,2, drop=F], wy=NULL)plot(impDRF[!is.na(X[,1]),c(‘X2′,’X1’)], main=paste(‘购物车插补’, ‘\nRMSE’, RMSEcalc(impcart, Xstar), ‘\nEnergy’, energycalc(impcart, Xstar)), cex=0.8, col=’darkblue’, cex.main=1.5)points(impDRF[is.na(X[,1]) ,c(‘X2′,’X1′)], col=’darkred’, cex=0.8 )plot(impnorm[!is.na(X[,1]),c(‘X2′,’X1’)] , main=paste(‘高斯插补’,’\nRMSE’, RMSEcalc(impnorm, Xstar), ‘\nEnergy’, energycalc(impnorm, Xstar)), col=’darkblue’, cex.main=1.5)points(impnorm [is.na(X[,1]),c(‘X2′,’X1′)], col=’darkred’, cex=0.8 )

如何处理缺失值引起的分布变化？发现填充缺失值的最佳插补算法

虽然mouse-cart 和mouse-DRF 都做得很好，但它们仍然不如高斯插值。这本身并不奇怪，因为高斯插值是这种情况下的理想插值方法（因为(X_1, X_2) 确实是高斯插值）。 X_2的分布变化可能会导致mice-cart和mice-DRF难以恢复3000个观测值的分布（这些方法通常非常有效）。

本文还讨论了一个类似但更极端的涉及两个变量（X_1，X_2）的示例。在这个例子中，分布变化更加显着，基于森林的方法相应地面临挑战：

分布变化的影响：当底层数据的分布发生显着变化时，基于模型的插值方法（例如基于决策树或随机森林的插值方法）可能难以准确恢复数据的真实分布。这是因为这些模型依赖于数据中存在的模式和结构，当这些模式和结构由于丢失数据而发生变化时，模型的有效性可能会降低。

高斯插值的适用性：高斯插值假设数据遵循高斯分布，这使得它在处理本质上是高斯分布的数据时非常有效。对于此类数据，高斯插值可以更好地捕获数据的整体统计特性，包括均值和协方差结构。

方法选择的重要性：选择正确的插补方法对于保证插补结果的质量至关重要。当面对复杂的数据和分布时，了解每种方法的优点和局限性可以帮助研究人员和从业者选择最适合其数据特征的插补方法。

总结

缺失值确实是一个棘手的问题。处理缺失值的最好方法就是尽可能避免它们的出现，但这几乎是不可能的，所以即使只考虑随机缺失（MAR），寻找插补方法的工作还远远没有结束。我们仍然缺乏一种方法可以（1）进行非参数分布预测和（2）适应MAR 下可能发生的分布变化。有时感觉人们把问题过于复杂化了，因为一些MICE 方法表现得非常好，可能足以解决许多缺失值问题。

有一些非常先进的机器学习方法，例如GAIN 及其变体，它们尝试使用神经网络来插值数据。这些方法可能是正确的，因为它们遵循正确的想法：估算观察到的情况中缺失的条件分布。但使用一段时间后，其表现并不尽如人意，尤其是与MICE相比。

所以如果我遇到缺失值问题，我首先可以尝试的是mice-cart或者论文中开发的新方法mice-DRF。他们重现数据的能力是惊人的。但这些方法至少需要适度的样本量，例如超过200 或300。如果您的观测值少于200 个，则可以选择更简单的方法，例如高斯插值。如果您想找到这些方法中最好的方法，请尝试使用本文中的评分标准来找到最佳方法。

最后，不要忘记这些方法都不能有效地处理插补不确定性！因为我们在本文中只讨论单一插补。正确的多重插补需要考虑插补方法本身的不确定性，这通常是使用贝叶斯方法完成的，这可能是另一篇论文。