多重插补

一、直观理解：与其编一个“假真相”，不如生成多本合理的完整账本

多重插补（Multiple Imputation, MI）承认一件事：缺失值不是有一个已知的正确答案，而是有一个条件分布。与其像均值插补那样填死一个数，不如：

1. 根据已观测数据，对缺失部分进行多次随机抽取，得到 $m$ 本彼此不同、但都合理的“完整数据集”；
2. 对每一本完整数据用同一套分析模型（回归、检验、评分）得到估计；
3. 再用 Rubin 规则把 $m$ 份估计合并，使最终的标准误同时包含“普通抽样波动”和“因为缺失而多出来的不确定”。

好比档案室丢了几页账：你不会只请一个人凭印象补一页就当真，而是请几位懂业务的会计各自补一版，业务结论若在多版之间都稳定，才更可信；若有的版本赚、有的版本亏，你必须把这种分歧写进不确定性里。

多重插补的产出重点往往不是“那一张被填好的表”，而是合并后的参数估计与区间。若目标只是给机器学习模型喂一个矩阵，MI 仍然有用（可对 $m$ 个模型预测做平均），但论文写法应与经典推断场景有所区分。

二、适用场景与优缺点

类型	说明
适用场景	缺失达到不可忽视的比例，且可假设 MCAR 或更常见的 MAR；需要参数估计、置信区间或假设检验，而不是单纯排行榜预测；有足够辅助变量帮助解释缺失；希望相对直接删除减少信息浪费，相对单值插补恢复不确定度。
核心优点	在合理模型下可得到近似无偏的点估计与更诚实的方差；能利用辅助变量；比“填完当完整数据”科学；结论可做插补间变异诊断。
局限与缺点	步骤与表述更复杂；插补模型定错会系统性偏差；MNAR 不能靠默认 MI 自动解决；实现细节多（$m$、迭代、变量类型、兼容性）；竞赛时限紧时工程成本高于均值 / KNN。

使用多重插补前通常需要完成以下检查：

1. 缺失比例、模式，以及与关键变量 / 结果的关系（支持 MAR 的证据或理由）；
2. 插补模型中是否纳入了分析模型里的变量及重要辅助变量（包容性原则）；
3. 变量类型（连续、二分类、有序、多项）是否有对应的条件模型；
4. 计划的 $m$、迭代次数、随机种子与收敛诊断；
5. 最终分析是合并参数，还是合并预测——二者流程不同。

三、核心原理：从条件分布抽样，再按 Rubin 规则合并

1. 插补阶段

在 MAR 假设下，缺失数据的条件分布由观测数据决定。实践中常用：

• 联合模型法：假设多元正态等联合分布，对缺失块做条件抽样；
• 全条件设定 / 链式方程（MICE）：对每个含缺失变量，用其余变量回归（或分类），轮转迭代，收敛后抽一版完整数据；重复 $m$ 次得 $m$ 版。

对连续变量的一个示意（概念上）：

$$x_{\mathrm{mis}}\mid x_{\mathrm{obs}} \sim \mathcal{N}\big(\hat\mu(x_{\mathrm{obs}}),\,\hat\sigma^2(x_{\mathrm{obs}})\big)$$

每次插补先考虑参数的不确定性（不是只用点估计 $\hat\mu$），再抽缺失值——这正是“多重”能够补上方差的关键。

2. 分析阶段

对第 $\ell$ 个完整数据集，拟合同一分析模型，得到目标量估计 $Q^{(\ell)}$ 及其（在该完整数据下的）方差 $U^{(\ell)}$。例如 $Q$ 可以是回归系数、组间均值差、相关系数等。

3. 合并阶段（Rubin 规则）

点估计：

$$\bar Q=\frac{1}{m}\sum_{\ell=1}^{m}Q^{(\ell)}$$

完整数据内平均方差：

$$\bar U=\frac{1}{m}\sum_{\ell=1}^{m}U^{(\ell)}$$

插补间方差：

$$B=\frac{1}{m-1}\sum_{\ell=1}^{m}\big(Q^{(\ell)}-\bar Q\big)^2$$

总方差：

$$T=\bar U+\Big(1+\frac{1}{m}\Big)B$$

$\sqrt{T}$ 即合并后的标准误。$B$ 相对 $\bar U$ 越大，说明缺失造成的额外不确定越大。

缺失信息分数的一个常用近似：

$$\lambda \approx \frac{(1+1/m)B}{T}$$

$\lambda$ 接近 0 表示缺失几乎不增加不确定；接近 1 表示结论高度受插补支配，需谨慎。

4. 与单值插补的本质差别

项目	均值 / KNN 等单值插补	多重插补
填补次数	1	$m$
是否抽随机	通常否	是
标准误	往往偏小	含插补不确定性
主要交付物	一张完整表	合并推断 + 诊断

四、关键决策 / 如何选参数

1. $m$ 取多少？

经典文献常用 $m=5$。更现代的经验是：

• 缺失较少、分析简单：$m=5\sim10$；
• 缺失较多或 $\lambda$ 较大：$m=20$ 甚至更高；
• 经验规则之一：$m$ 大致不低于“平均缺失率的百分数”（如 15% 缺失 → $m\approx 15+$），再结合计算成本。

应报告：增大 $m$ 后 $\bar Q$ 与 $T$ 是否稳定。

2. 插补模型里放哪些变量？

遵循包容性：分析模型中的因变量、自变量、交互所涉变量，以及与缺失强相关的辅助变量，都应进入插补模型。漏掉因变量再去插自变量，可能把真正关联“洗白”。

3. 迭代与收敛（MICE）

每个变量一轮条件更新，通常 10–20 次迭代后再导出一版；可通过轨迹图查看均值 / 方差是否稳定。不收敛时检查高度共线、完美预测、稀疏类别。

4. 变换与约束

偏态正值可在对数尺度插补再还原；比例数据注意边界；需保持的逻辑约束（如部分之和为总量）要在插补后检查或改用更结构化模型。

5. 预测任务怎么用 MI？

对每个测试样本可得到 $m$ 个填补特征，从而 $m$ 个预测，再平均（分类可平均概率）。不要只随便留一版完整矩阵当“真理”。

6. 随机种子

必须固定种子以保证复现，并在附录给出。

五、MATLAB 示例

MATLAB 没有像 R mice 那样开箱即用的完整 MI 生态，但可以用自举 + 条件回归残差抽样演示 Rubin 合并的完整流程（教学与小规模数值问题足够；正式研究可对接专门实现或自定义 MICE）。下列示例目标量 $Q$ 取“$x_2$ 对 $x_1$ 的回归斜率”。

rng(11);

% ----- 生成 MAR 缺失：x2 缺失概率随 x1 增大 -----
n = 120;
x1 = randn(n, 1);
x2 = 0.7 * x1 + 0.6 * randn(n, 1);
y  = 1.5 + 0.8 * x1 + 0.5 * x2 + 0.5 * randn(n, 1);

pMiss = 1 ./ (1 + exp(-(0.3 + 0.9 * x1)));  % MAR
x2_obs = x2;
x2_obs(rand(n, 1) < pMiss) = NaN;

fprintf("x2 缺失率 = %.1f%%\n", 100 * mean(isnan(x2_obs)));

% ----- 多重插补：m 次 -----
m = 20;
Q = zeros(m, 1);   % 斜率：y ~ x1 + x2 中 x2 的系数
U = zeros(m, 1);   % 对应估计方差

for ell = 1:m
    % 1) 在观测到 x2 的样本上拟合 x2 ~ x1
    obs = ~isnan(x2_obs);
    Xdes = [ones(sum(obs),1), x1(obs)];
    beta = Xdes \ x2_obs(obs);
    resid = x2_obs(obs) - Xdes * beta;
    sigma = std(resid, 0);

    % 2) 对参数做简单贝叶斯/自助式扰动，体现参数不确定
    %    这里用残差自助 + 重新回归
    bootIdx = randi(sum(obs), sum(obs), 1);
    beta_b = Xdes(bootIdx, :) \ x2_obs(obs);
    resid_b = x2_obs(obs) - Xdes * beta_b;
    sigma_b = std(resid_b, 0);
    if sigma_b < 1e-8, sigma_b = sigma; end

    % 3) 对缺失位条件抽样
    x2_fill = x2_obs;
    mis = isnan(x2_obs);
    muMis = [ones(sum(mis),1), x1(mis)] * beta_b;
    x2_fill(mis) = muMis + sigma_b * randn(sum(mis), 1);

    % 4) 在第 ell 本完整数据上做目标分析
    Xanal = [ones(n,1), x1, x2_fill];
    [b, ~, ~, ~, stats] = regress(y, Xanal);
    Q(ell) = b(3);                 % x2 的系数
    % regress 的 stats(4) 不是系数协方差；用手动估计方差
    yhat = Xanal * b;
    s2 = sum((y - yhat).^2) / (n - 3);
    covB = s2 * inv(Xanal' * Xanal);
    U(ell) = covB(3, 3);
end

% ----- Rubin 合并 -----
Qbar = mean(Q);
Ubar = mean(U);
B = var(Q, 0);                    % 样本方差，分母 m-1
T = Ubar + (1 + 1/m) * B;
se = sqrt(T);
lambda = (1 + 1/m) * B / T;       % 缺失信息分数近似

fprintf("合并后的 x2 系数 Qbar = %.4f\n", Qbar);
fprintf("SE = %.4f\n", se);
fprintf("95%% 近似区间 = [%.4f, %.4f]\n", Qbar - 1.96*se, Qbar + 1.96*se);
fprintf("插补间方差 B = %.4f, 平均期内方差 Ubar = %.4f\n", B, Ubar);
fprintf("缺失信息分数 lambda ≈ %.3f\n", lambda);

% 对照：均值单值插补后的“天真”回归 SE（通常偏小）
x2_mean = x2_obs;
x2_mean(isnan(x2_mean)) = mean(x2_obs, "omitnan");
Xsingle = [ones(n,1), x1, x2_mean];
bS = Xsingle \ y;
yhatS = Xsingle * bS;
s2S = sum((y - yhatS).^2) / (n - 3);
seS = sqrt(s2S * inv(Xsingle' * Xsingle));
fprintf("单值均值插补: 系数=%.4f, SE=%.4f（往往小于 MI 的 SE）\n", bS(3), seS(3,3)^0.5);

% 插补间估计的分散情况
figure;
histogram(Q, 10);
xline(Qbar, "r", "LineWidth", 1.5);
xlabel("Q^{(\ell)}");
title("m 次插补下的目标估计分布");
grid on;

说明：上例的插补模型有意写得可读，用于掌握 MI 思想；变量多、类型杂时，应把链式方程、类别变量逻辑回归等写完整，并做收敛检查。

六、如何解释结果

重点解释合并后的推断，而不是某一版填补表上的某一个格子。

建议报告：

1. 缺失机制的判断与依据（至少说明为何 MAR 可接受或局限）；
2. 插补方法（联合模型 / MICE / 自定义条件回归）、进入插补的变量列表；
3. $m$、迭代次数、随机种子；
4. 每个目标参数的 $\bar Q$、$\bar U$、$B$、$T$、置信区间；
5. 缺失信息分数 $\lambda$ 或等价诊断；
6. 与直接删除、单值插补在点估计与区间宽度上的对照；
7. 敏感性分析：去掉某辅助变量、改变 $m$、或简单 MNAR 情景下结论是否翻盘。

若 $m$ 次得到的 $Q^{(\ell)}$ 彼此差很远，应原话告诉评委：证据被缺失严重稀释，而不是只给一个漂亮的点估计。

七、常见误区

1. 做了 $m$ 次插补却只留用第 1 份数据

这与单值插补无异，还浪费了随机性；必须分析全部 $m$ 份并合并。

2. 插补模型漏掉分析模型中的变量

最典型：插补自变量时不用 $y$，导致 $y$ 与 $x$ 的关系衰减。

3. 把插补值当成测量真值写进原始数据附录

应保留缺失指示，并说明完整数据是分析用中间产物。

4. 默认 MI 能解决 MNAR

若“越高越不填”，需要更强假设、外部信息或专门模型，不能靠增加 $m$ 硬扛。

5. 只比较点估计，不看区间与 $\lambda$

MI 的核心贡献之一是更诚实的不确定度；忽略 $T$ 就丢掉了方法灵魂。

八、小结

多重插补通过从缺失数据的条件分布中抽取 $m$ 个合理完整版本，分别分析后再用 Rubin 规则合并，使结论同时反映抽样误差与缺失带来的额外不确定。

它的核心思想可以概括为：

不要伪造一个唯一的填补真相；生成多本说得通的账，若结论仍稳定，才值得相信，并把版本之间的分歧写进误差里。

在 MAR 且插补模型足够包容时，MI 通常优于直接删除和单值插补的天真推断。使用时把机制假设、$m$ 与合并统计量写清楚，并保留简单方法作对照。对竞赛预测，可以把 $m$ 次预测平均；对科研式数模推断，务必报告 $\bar Q$ 与 $T$，而不是只贴一张“填好的表”。