多重插补的mice包学习
Congyuan He (hcy)
2026-04-03
1 MCAR MAR MNAR
完全随机缺失（MCAR，Missing Completely At Random）：数据的缺失与任何变量（包括已观测和未观测变量）无关，即缺失是完全随机的。在这种情况下，缺失数据不会引入系统性偏差，因此分析结果仍具有代表性。然而，现实中很少出现完全随机缺失的情况。
随机缺失（MAR，Missing At Random）：数据的缺失仅与已观测到的变量有关，而与未观测到的变量无关。这意味着，缺失的发生可以通过其他已观测变量来解释。在这种情况下，使用适当的统计方法（如多重插补）可以有效地处理缺失数据。最一般的假设。
非随机缺失（MNAR，Missing Not At Random）：数据的缺失与未观测到的变量相关，即缺失本身与数据的值有关。例如，收入较高的人可能不愿透露他们的收入水平。这种情况下，处理缺失数据变得更加复杂，可能需要对缺失机制进行建模或收集更多数据来解决。
2 常规处理方式
Listwise deletion ： 删除整行的观测，在MCAR假设（强假设）下可以得到无偏估计，优点简单易行，缺点是数据损失大。如果不满足MCAR将产生估计偏差。
Pairwise deletion ： 不直接删除包含缺失值的整个观测记录，而是根据每个分析所需变量的实际可用数据进行计算，例如算均值加NA.rm=TRUE。优点在于最大限度的利用数据 缺点在于样本量不一致，同样的若不满足MCAR估计是有偏的
Mean imputation：用均值替代缺失值。优点是简便易行，缺点是低估方差扰乱变量之间的关系，使得除了平均值意外的几乎所有估计值产生偏差，若不满足MCAR，均值的估计也会有偏，这种策略应该在存在少量的缺失时快速修补并且应该尽量避免使用
说明的例子：
## Wind Temp Month Day Solar.R Ozone
## 111 1 1 1 1 1 1 0
## 35 1 1 1 1 1 0 1
## 5 1 1 1 1 0 1 1
## 2 1 1 1 1 0 0 2
## 0 0 0 0 7 37 44
- Regression imputation: 利用回归模型的预测值来替代缺失值,优点：均值和回归系数在特定条件下无偏，预测接近完美时能较好地重构缺失数据；缺点：相关性容易被高估，可能产生虚假的关系，数据变异性被低估
fit <- lm(Ozone ~ Solar.R, data = airquality, na.action = na.omit)
pred <- predict(fit, newdata = ic(airquality)) # 预测填补了ozone
data <- airquality[, c("Ozone", "Solar.R")]
imp <- mice(data, method = "norm.predict",print = FALSE)# 正态预测法（基于线性回归模型）来预测缺失值
xyplot(imp, Ozone ~ Solar.R)
- Stochastic regression imputation：随机回归插补是回归插补尝试通过向预测添加噪声来解决相关偏差的改进版本 优点：保持回归系数和变量相关性，添加噪声扩展了预测分布，保留数据变异性；缺点：可能产生不合理的负值
data <- airquality[, c("Ozone", "Solar.R")]
imp <- mice(data, method = "norm.nob",print = FALSE)
xyplot(imp, Ozone ~ Solar.R)
- LOCF and BOCF：LOCF,Last observation carried forward （最后观测值前移）和BOCF,baseline observation carried forward （基线观测值前移）是处理纵向数据缺失的一种临时插补方法。其基本思想是用之前的观测值来替代缺失数据。当连续多个值缺失时，该方法会寻找上一次的观测值作为替代。
# tidyr包中的fill函数就是应用LOCF
airquality2 <- tidyr::fill(airquality, Ozone)
plot(c(1:nrow(airquality2)),airquality2$Ozone,type = 'l')
- Indicator method 假设我们想拟合一个回归模型，但某个解释变量存在缺失值。指示器法将每个缺失值替换为零或均值，并在回归模型中添加一个响应指标来扩展模型。该方法适用于每个不完整的变量，最终用户分析的是扩展后的模型，而不是原始模型。
指标法的优点在于它能够保留完整的数据集，不会丢失样本信息；同时，通过加入响应指示符，该方法能捕捉观察数据与缺失数据之间的系统性差异，在特定条件下（如随机试验、模型为线性且无交互作用时）还可以无偏估计治疗效果。然而，对于多个不完整的基线协变量，其置信区间覆盖率可能不理想，而且该方法依赖于较严格的假设条件，若条件不满足则可能引入偏倚。
library(data.table)
imp <- mice(airquality, method = "mean", m = 1,
maxit = 1, print = FALSE)
airquality2 <- cbind(complete(imp),
r.Ozone = is.na(airquality[, "Ozone"]))
setDT(airquality2)
airquality2[r.Ozone==TRUE,Ozone:=0]
fit <- lm(Wind ~ Ozone + r.Ozone, data = airquality2)
summary(fit)##
## Call:
## lm(formula = Wind ~ Ozone + r.Ozone, data = airquality2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.557 -2.257 -0.261 1.760 10.504
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.608430 0.453088 27.828 < 2e-16 ***
## Ozone -0.065189 0.008482 -7.686 1.84e-12 ***
## r.OzoneTRUE -2.351673 0.669776 -3.511 0.00059 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3 on 150 degrees of freedom
## Multiple R-squared: 0.2842, Adjusted R-squared: 0.2747
## F-statistic: 29.78 on 2 and 150 DF, p-value: 1.285e-11##
## Call:
## lm(formula = Wind ~ Ozone, data = airquality)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2521 -2.2792 -0.2376 1.7474 10.5036
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.608430 0.433126 29.11 < 2e-16 ***
## Ozone -0.065189 0.008108 -8.04 9.27e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.868 on 114 degrees of freedom
## (37 observations deleted due to missingness)
## Multiple R-squared: 0.3619, Adjusted R-squared: 0.3563
## F-statistic: 64.64 on 1 and 114 DF, p-value: 9.272e-13- 总结
3 Multiple imputation
- 为什么用多重插补？
因为其解决了之前常规策略中标准误too small的问题（前表），其提供一种机制来处理插补本身的不确定性
它将缺失数据问题的解与完整数据问题的解分开。首先解决缺失数据问题，然后解决完整数据问题。将这两个阶段分开的能力简化了统计建模，因此有助于更好地了解科学研究现象。
- 一个简单的案例
imp <- mice(airquality, seed = 1, m = 20, print = FALSE) # imputation
# 默认的方法根据缺失数据来判定 可见帮助文档 defaultMethod
fit <- with(imp, lm(Ozone ~ Wind + Temp + Solar.R)) # analysis
summary(pool(fit)) # pooling## term estimate std.error statistic df p.value
## 1 (Intercept) -65.87829658 23.09377412 -2.852643 69.97033 5.696702e-03
## 2 Wind -3.01897171 0.66252377 -4.556775 70.51194 2.125022e-05
## 3 Temp 1.63483547 0.25107557 6.511328 75.99913 7.203792e-09
## 4 Solar.R 0.05861581 0.02267832 2.584662 90.10797 1.135441e-02##
## Call:
## lm(formula = Ozone ~ Wind + Temp + Solar.R, data = airquality)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.485 -14.219 -3.551 10.097 95.619
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -64.34208 23.05472 -2.791 0.00623 **
## Wind -3.33359 0.65441 -5.094 1.52e-06 ***
## Temp 1.65209 0.25353 6.516 2.42e-09 ***
## Solar.R 0.05982 0.02319 2.580 0.01124 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.18 on 107 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.6059, Adjusted R-squared: 0.5948
## F-statistic: 54.83 on 3 and 107 DF, p-value: < 2.2e-16- 不推荐用IM的场景
- 若缺失仅出现在回归模型的Y变量，且分析目标为回归系数，直接删除缺失样本与多重插补效果等同，此时优先选择更简单高效的CC方法
2.当数据缺失概率仅依赖X变量（与Y无关），且模型设定正确时，直接删除缺失样本分析回归系数无偏。此时不推荐多重插补。
- 若缺失仅出现在Y或X单一方（不同时缺失），且缺失概率仅依赖Y（如病例对照研究设计），直接删除缺失样本分析的回归系数仍无偏，可替代多重插补。
- 插补次数的选择
传统建议：当缺失信息比例较低时，如≤30%，3-5次即可计算成本低
如果需要高计算精度，那么有一个规则 γ=0.3 → m=24；γ=0.5 → m=59；γ=0.9 → m=258
实用好记的 插补次数可以设置为缺失百分比
实际操作：模型测试调整阶段 m=5, 高精度需求最终计算阶段m=20~100
4 Imputation in practice
- 默认的根据数据类型可够选择的变量插补方法
- 查看插补变量的predictors
对于中小型数据集（最多包含 20-30 个变量，没有派生变量、交互效应和其他复杂性），选择尽可能多的预测因子通常是有益的。包括尽可能多的预测因子往往会使 MAR 假设更合理，从而减少对 MNAR 机制进行特殊调整的需要
## age bmi hyp chl
## age 0 1 1 1
## bmi 1 0 1 1
## hyp 1 1 0 1
## chl 1 1 1 0对于包含众多变量的大型数据库，不必将所有变量都纳入插补模型。一般来说，在用于线性回归的完整数据模型中，最佳预测变量数量通常在15个左右。以下是选择预测变量的策略总结：
必须包含完整数据模型中的变量
加入与缺失相关的变量
选择解释力强的变量
剔除在不完整子集中缺失过多的变量
- 插补的诊断
##
## iter imp variable
## 1 1 bmi hyp chl
## 1 2 bmi hyp chl
## 1 3 bmi hyp chl
## 1 4 bmi hyp chl
## 1 5 bmi hyp chl
## 2 1 bmi hyp chl
## 2 2 bmi hyp chl
## 2 3 bmi hyp chl
## 2 4 bmi hyp chl
## 2 5 bmi hyp chl
## 3 1 bmi hyp chl
## 3 2 bmi hyp chl
## 3 3 bmi hyp chl
## 3 4 bmi hyp chl
## 3 5 bmi hyp chl
## 4 1 bmi hyp chl
## 4 2 bmi hyp chl
## 4 3 bmi hyp chl
## 4 4 bmi hyp chl
## 4 5 bmi hyp chl
## 5 1 bmi hyp chl
## 5 2 bmi hyp chl
## 5 3 bmi hyp chl
## 5 4 bmi hyp chl
## 5 5 bmi hyp chl
- 再来一个例子
# 处理数据 一个儿童发育的数据
select <- with(boys, age >= 8 & age <= 21.0)
djm <- boys[select, -4]
djm$gen <- as.integer(djm$gen)
djm$phb <- as.integer(djm$phb)
djm$reg <- as.integer(djm$reg)
str(djm)## 'data.frame': 424 obs. of 8 variables:
## $ age: num 8.83 8.86 8.87 8.91 8.93 ...
## $ hgt: num 134 125 145 138 140 ...
## $ wgt: num 31.2 31 38.2 30 37.2 24.9 26.9 48.2 29.4 30 ...
## $ hc : num 53.4 51.6 54.8 54.7 55.9 53.8 53.4 50.5 53.9 53.6 ...
## $ gen: int NA 1 2 1 NA NA 1 2 1 NA ...
## $ phb: int NA 1 1 1 NA NA 1 1 1 NA ...
## $ tv : int NA 1 2 2 NA NA 4 2 2 NA ...
## $ reg: int 2 1 4 4 4 4 2 3 2 4 ...# age 表示十进制年龄（0-21 岁）；
# hgt 为身高（单位：厘米）；
# wgt 为体重（单位：公斤）；
# bmi 代表身体质量指数；
# hc 表示头围（单位：厘米）；
# gen 为生殖器 Tanner 分期（G1-G5）；
# phb 指阴毛 Tanner 分期（P1-P6）；
# tv 为睾丸体积（单位：毫升）；
# reg 则记录地区信息，包括北部、东部、西部、南部和城市。
dfcs <- boys[select, -4]
str(dfcs)## 'data.frame': 424 obs. of 8 variables:
## $ age: num 8.83 8.86 8.87 8.91 8.93 ...
## $ hgt: num 134 125 145 138 140 ...
## $ wgt: num 31.2 31 38.2 30 37.2 24.9 26.9 48.2 29.4 30 ...
## $ hc : num 53.4 51.6 54.8 54.7 55.9 53.8 53.4 50.5 53.9 53.6 ...
## $ gen: Ord.factor w/ 5 levels "G1"<"G2"<"G3"<..: NA 1 2 1 NA NA 1 2 1 NA ...
## $ phb: Ord.factor w/ 6 levels "P1"<"P2"<"P3"<..: NA 1 1 1 NA NA 1 1 1 NA ...
## $ tv : int NA 1 2 2 NA NA 4 2 2 NA ...
## $ reg: Factor w/ 5 levels "north","east",..: 2 1 4 4 4 4 2 3 2 4 ...
## age reg hgt wgt hc gen phb tv
## 223 1 1 1 1 1 1 1 1 0
## 18 1 1 1 1 1 1 1 0 1
## 1 1 1 1 1 1 1 0 1 1
## 1 1 1 1 1 1 0 1 0 2
## 146 1 1 1 1 1 0 0 0 3
## 33 1 1 1 1 0 0 0 0 4
## 1 1 1 0 0 1 1 1 1 2
## 1 1 1 0 0 0 1 1 1 3
## 0 0 2 2 34 180 180 198 596# 用以下三种不同方式进行插补
# normal model
jm <- mice(djm, method = "norm", seed = 93005, m = 10,
print = FALSE)
# predictive mean matching
pmm <- mice(djm, method = "pmm", seed = 71332, m = 10,
print = FALSE)
# the proportional odds model
fcs <- mice(dfcs, seed = 81420, m = 10, print = FALSE) # 看表默认的是plor
# 看一下插补结果
xyplot(jm, gen ~ age | as.factor(.imp), subset = .imp < 6,# 取前5个数据集展示
xlab = "Age", ylab = "Genital stage", col = mdc(1:2),
ylim = c(0, 6))
# 上图你会发现对于gen的插补是连续的你还需要进一步的将其转化为分类的变量比如4舍5入
xyplot(pmm, gen ~ age | as.factor(.imp), subset = .imp < 6,# 取前5个数据集展示
xlab = "Age", ylab = "Genital stage", col = mdc(1:2),
ylim = c(0, 6))
# 直接就是分类的变量了
xyplot(fcs, gen ~ age | as.factor(.imp), subset = .imp < 6,# 取前5个数据集展示
xlab = "Age", ylab = "Genital stage", col = mdc(1:2),
ylim = c(0, 6))
