参考资料 https://stefvanbuuren.name/fimd/ch-practice.html

1 缺失值的情况

  • 如果缺失概率在所有情况下都相同,那么数据就被称为完全随机缺失(MCAR)

  • 如果缺失概率仅在观察数据所定义的组内相同,那么数据就是随机缺失(MAR)

  • 如果 MCAR 和 MAR 都不成立,那么我们就称之为非随机缺失(MNAR),处理 MNAR 的策略是找到更多有关缺失原因的数据,或进行假设分析,以了解结果在不同情况下的敏感度

  • 大多数简单的修复方法只有在限制性且往往不现实的 MCAR 假设下才能起作用。如果 MCAR 不可信,这些方法就会提供有偏差的估计值。

2 缺失值的解决方案

2.1 Listwise deletion

函数 na.omit()和complete.cases() 删除分析变量中存在一个或多个缺失值的所有案例。

fit <- lm(Ozone ~ Wind, data = airquality, na.action = na.omit)
naprint(na.action(fit))#37个观测在拟合中被删除
## [1] "因为不存在,37个观察量被删除了"
colSums(is.na(airquality))#这是因为有37个ozone缺失值
##   Ozone Solar.R    Wind    Temp   Month     Day 
##      37       7       0       0       0       0
fit2 <- lm(Ozone ~ Wind + Solar.R, data = airquality)
naprint(na.action(fit2))#在多变量回归时会带来更多的观测被删除,导致解释起来不一致
## [1] "因为不存在,42个观察量被删除了"

Complete-case analysis (listwise deletion)完整案例分析的一个重要优点是方便。如果数据是 MCAR 数据,列表式删除会产生无偏的均值、方差和回归权重估计值。

列表删除法的缺点是可能造成浪费。在实际应用中,丢失一半以上原始样本的情况并不少见,尤其是在变量数量较多的情况下。

如果数据不是 MCAR,列表删除会严重影响均值、回归系数和相关系数的估计值。

然而,完整案例分析并不总是坏事。根据缺失数据出现的位置(结果或预测因子)以及完整数据模型的参数和模型形式,缺失数据的影响是不同的。在回归分析中,列表删除具有一些独特的特性,使其在特定情况下具有吸引力。在有些情况下,列表删除甚至比最复杂的程序都能提供更好的估计值。(https://stefvanbuuren.name/fimd/sec-when.html

2.2 Pairwise deletion

成对删除法又称可用案例分析法(available-case analysis),试图弥补列表删除法的数据丢失问题。该方法计算所有观测数据的均值和(共)方差。因此,变量 X 的均值基于所有观测到 X 数据的案例,变量 Y 的均值基于所有观测到 Y 值的案例,以此类推。至于相关性和协方差,则取 X 和 Y 均不缺失的所有数据。然后,将汇总统计矩阵输入程序,进行回归分析、因子分析或其他建模程序。

SPSS、SAS 和 Stata 中的许多程序都有成对删除选项。在 R 中,我们可以计算空气质量数据在成对删除情况下的均值和相关性:

data <- airquality[, c("Ozone", "Solar.R", "Wind")]
mu <- colMeans(data, na.rm = TRUE)
cv <- cov(data, use = "pairwise")

标准 lm() 函数不将均值和协方差作为输入,但 lavaan 软件包(Rosseel,2012 年)提供了这一功能:

library(lavaan)
## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.
fit <- lavaan("Ozone ~ 1 + Wind + Solar.R
              Ozone ~~ Ozone",
             sample.mean = mu, sample.cov = cv,
             sample.nobs = sum(complete.cases(data)))
summary(fit)
## lavaan 0.6.17 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         4
## 
##   Number of observations                           111
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   Ozone ~                                             
##     Wind             -5.545    0.644   -8.605    0.000
##     Solar.R           0.118    0.025    4.681    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Ozone            75.402    8.463    8.910    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Ozone           565.035   75.845    7.450    0.000

该方法也有一些缺点。首先,如果数据不是 MCAR,估计值可能会有偏差。此外,协方差和/或相关矩阵可能不是正定的,而这是大多数多变量程序的要求。对于高度相关的变量,问题通常更为严重。计算标准误差时应使用多大的样本量尚不清楚。取平均样本量得出的标准误差太小。此外,成对删除法要求数值数据遵循近似正态分布,而在实际中,我们经常会遇到混合类型的变量。

如果数据近似于多元正态分布,变量之间的相关性较低,且 MCAR 假设可信,则配对删除法效果最佳。其他情况下不推荐使用。

2.3 Mean imputation

快速解决缺失数据的方法是用平均值代替。对于分类数据,我们可以使用模式。假设我们想在 Ozone 和 Solar.R 中估算空气质量数据的平均值。现在可以通过以下方法计算每个变量的均值

#install.packages("mice")
library("mice")
## 
## 载入程辑包:'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
imp <- mice(airquality, method = "mean", m = 1, maxit = 1)
## 
##  iter imp variable
##   1   1  Ozone  Solar.R

参数 method = mean 指定了均值估算,参数 m = 1 要求使用单个估算数据集,而 maxit = 1 则将迭代次数设为 1(无迭代)。后两个选项可以保持默认值,结果基本相同。

均值估算是一种快速、简单的缺失数据修复方法。但是,它会低估方差,扰乱变量之间的关系,使除均值以外的几乎所有估计值产生偏差,并在数据不存在 MCAR 时使均值估计值产生偏差。

平均值估算也许只能在少数数值缺失的情况下作为一种快速修复方法使用,一般情况下应避免使用。

2.4 Regression imputation

回归归因结合了其他变量的知识,目的是产生更智能的归因。第一步是根据观测数据建立模型。然后,根据拟合模型计算不完整案例的预测值,作为缺失数据的替代。假设我们通过 Solar.R 的线性回归来预测臭氧。

fit <- lm(Ozone ~ Solar.R, data = airquality)
pred <- predict(fit, newdata = ic(airquality))#ic提取缺失值

另一种回归估算方法是使用mice:

data <- airquality[, c("Ozone", "Solar.R")]
imp <- mice(data, method = "norm.predict", seed = 1,
           m = 1, print = FALSE)
xyplot(imp, Ozone ~ Solar.R)

回归估算与均值估算一样,在 MCAR 条件下对均值的估算是无偏的,如果解释变量是完整的,对估算模型的回归权重的估算也是无偏的。此外,如果影响缺失的因素是回归模型的一部分,那么回归权重在 MAR 下也是无偏的。在本例中,Solar.R 可以解释臭氧缺失概率的任何差异。另一方面,相关性会向上偏移,估算数据的可变性会被系统性低估。低估的程度取决于解释的方差和缺失案例的比例。

回归估算,以及其在机器学习中的现代化身,可能是本文介绍的所有方法中最危险的一种。我们可能会认为,我们应该通过保留变量之间的关系来做好工作。但实际上,回归估算人为地强化了数据中的关系。相关性向上偏移。变异性被低估。估算结果好得不真实。回归估算是产生假阳性和虚假关系的秘诀。

2.5 Stochastic regression imputation

随机回归估算是回归估算的一种改进,试图通过在预测中加入噪声来解决相关性偏差问题。以下代码通过随机回归估算法对 Solar.R 中的臭氧进行估算。

data <- airquality[, c("Ozone", "Solar.R")]
imp <- mice(data, method = "norm.nob", m = 1, maxit = 1,
            seed = 1, print = FALSE)
xyplot(imp, Ozone ~ Solar.R)

method = norm.nob 参数请求使用普通的非贝叶斯随机回归方法。该方法首先估计线性模型下的截距、斜率和残差方差,然后计算每个缺失值的预测值,并从残差中随机抽取一个值添加到预测值中。具体细节(https://stefvanbuuren.name/fimd/sec-linearnormal.html#sec:linearnormal)。

seed参数使解决方案具有可重复性。如图 显示,在预测值中加入噪声后,估算值的分布就会分散。

请注意,会出现一些新的复杂情况。有一个估算值为负值。这些负值不一定是由于从残差分布中抽取的结果,也可能是对非负数据使用线性模型的结果。事实上,图中 显示了观测数据中的几个负预测值。

然而,随机回归估算是概念上的一大进步。有些分析师可能会认为,通过添加随机噪声来 “破坏”最佳预测结果有违直觉,但这恰恰是随机回归适用于估算的原因。执行良好的随机回归估算不仅保留了回归权重,还保留了变量之间的相关性。从残差中提取的主要思路非常强大,是更高级估算技术的基础。

2.6 LOCF and BOCF

最后观测值前移(LOCF)和基线观测值前移(BOCF)是纵向数据的临时估算方法。其原理是用前一个观测值替代缺失数据。当多个值连续缺失时,该方法会搜索最后一个观察值。

tidyr 软件包中的函数fill()通过填充最后一个已知值来应用 LOCF。这在记录数值变化的情况下非常有用,比如时间到事件数据。例如,我们可以使用 LOCF 在臭氧中填入以下值

airquality2 <- tidyr::fill(airquality, Ozone)

LOCF 很方便,因为它能生成一个完整的数据集。在我们可以确定缺失值的情况下,例如纵向数据中的管理变量,可以放心地使用 LOCF。但对于结果而言,LOCF 方法则值得商榷。偏差可能是双向的,即使在 MCAR 条件下,LOCF 也可能产生有偏差的估计值。在使用 LOCF 后,需要使用适当的统计分析方法来区分真实数据和估算数据。然而,这通常是做不到的。

2.7 Indicator method

假设我们想要拟合一个回归结果,但其中一个解释变量存在缺失值。指标法将每个缺失值替换为零,并通过响应指标扩展回归模型。该过程适用于每个不完整变量。用户分析的是扩展模型而不是原始模型。

imp <- mice(airquality, method = "mean", m = 1,
            maxit = 1, print = FALSE)
airquality2 <- cbind(complete(imp),#补充完整数据
                     r.Ozone = is.na(airquality[, "Ozone"]))#添加缺失值指标
fit <- lm(Wind ~ Ozone + r.Ozone, data = airquality2)#将指标作为变量

请注意,由于缺失的数据是臭氧数据,我们需要将回归模型的方向反过来。

指标法在公共卫生和流行病学领域很受欢迎。指标法的优点是保留了完整的数据集。此外,通过纳入反应指标,它还允许观察到的数据和未观察到的数据之间存在系统性差异,因此可能更有效。

对于多个不完整基线协变量,治疗估计值周围置信区间的覆盖率是否令人满意尚不清楚。指标法发挥作用的条件在实践中可能无法满足。指标法在特定情况下可能有其用途,但作为处理缺失数据的通用方法,指标法是失败的。

2.8 Multiple imputation

多重估算可创建 m>1 个完整数据集。每个数据集都由标准分析软件进行分析。根据集合规则(“鲁宾规则”),将 m 个结果集合为最终的点估计加标准误差。

继续使用空气质量数据集,可以直接应用多重估算。以下代码对缺失数据进行了二十次估算,拟合了一个线性回归模型来预测每个估算数据集中的臭氧,汇集了二十组估计参数,并计算了用于检验权重显著性的 Wald 统计量。

imp <- mice(airquality, seed = 1, m = 20, print = FALSE)
fit <- with(imp, lm(Ozone ~ Wind + Temp + Solar.R))
summary(pool(fit))

更多理论知识参见:https://stefvanbuuren.name/fimd/ch-mi.html