如何在Python中使用Rubin's Rules合并多重插补数据集的分析结果

Rubin规则本身是不绑定具体分析模型的通用合并逻辑，你只需要从每个插补数据集拿到对应指标的点估计值和标准误，就可以完成合并，完全可以覆盖线性回归、决策树之外的任意分析场景。

通用手动实现方案（适配所有模型类型） #

这个方案灵活性最高，不管是生存分析、时间序列模型、各类机器学习模型都可以使用，核心逻辑就是直接复现Rubin规则的计算步骤：

• 步骤1：遍历所有M个插补数据集，拟合你需要的分析模型，分别收集每个指标（如回归系数、特征重要性、风险比等）的点估计值Q_m，以及对应指标的方差估计值U_m（标准误的平方）
• 步骤2：计算合并后的点估计：所有插补集点估计的均值 $\bar{Q} = \frac{1}{M}\sum_{m=1}^M Q_m$
• 步骤3：计算插补内方差：所有插补集方差的均值 $\bar{U} = \frac{1}{M}\sum_{m=1}^M U_m$
• 步骤4：计算插补间方差：不同插补集点估计的样本方差 $B = \frac{1}{M-1}\sum_{m=1}^M (Q_m - \bar{Q})^2$
• 步骤5：计算合并后的总方差 $T = \bar{U} + (1 + \frac{1}{M})B$，对应的合并标准误为$\sqrt{T}$

以下是用逻辑回归做分析的示例代码，你可以替换为任意自定义模型：

import numpy as np
import statsmodels.api as sm
from miceforest import load_kernel

# 加载已经训练好的miceforest插补内核
kernel = load_kernel("your_imputation_kernel.dat")
M = kernel.dataset_count() # 插补数据集数量

coef_collect = []
se_collect = []

for m in range(M):
    # 取出第m个插补后的完整数据集
    imp_df = kernel.complete_data(dataset=m)
    X = sm.add_constant(imp_df.drop("target", axis=1))
    y = imp_df["target"]
    # 拟合分析模型，此处可替换为任意你需要的模型
    model = sm.Logit(y, X).fit(disp=False)
    coef_collect.append(model.params.values)
    se_collect.append(model.bse.values)

# 按Rubin规则合并结果
coef_arr = np.array(coef_collect)
se_arr = np.array(se_collect)

pooled_coef = np.mean(coef_arr, axis=0)
u_bar = np.mean(se_arr ** 2, axis=0)
b = np.var(coef_arr, axis=0, ddof=1)
total_var = u_bar + (1 + 1/M) * b
pooled_se = np.sqrt(total_var)

# 输出合并后的结果
for col, coef, se in zip(X.columns, pooled_coef, pooled_se):
    print(f"变量{col}：合并系数={coef:.4f}，合并标准误={se:.4f}")

如果你的模型没有解析标准误（比如随机森林、XGBoost这类非参模型），可以在每个插补集上用bootstrap抽样计算目标指标的标准误，再带入上述公式即可。

基于现有库扩展实现（省去规则复现步骤） #

你也可以基于Autoimpute提供的MiBaseRegressor基类扩展，自定义适配你的分析模型，不需要自己从头实现Rubin规则的计算逻辑：

from autoimpute.imputations import MultipleImputer
from autoimpute.analysis import MiBaseRegressor
from sklearn.ensemble import RandomForestRegressor
import numpy as np

# 自定义适配随机森林的多重插补分析类
class MiRandomForest(MiBaseRegressor):
    def _fit_model(self, X, y):
        # 拟合你自己的分析模型
        model = RandomForestRegressor(n_estimators=100, random_state=42)
        model.fit(X, y)
        # 收集要合并的指标：此处用特征重要性为例
        params = model.feature_importances_
        # Bootstrap计算标准误
        boot_imp = []
        for _ in range(100):
            idx = np.random.choice(len(X), len(X), replace=True)
            boot_m = RandomForestRegressor(n_estimators=100, random_state=42).fit(X.iloc[idx], y.iloc[idx])
            boot_imp.append(boot_m.feature_importances_)
        se = np.std(boot_imp, axis=0, ddof=1)
        return {"params": params, "se": se}

# 调用方式
mi = MultipleImputer(n=5, strategy="random forest")
imputed_data = mi.fit_transform(your_raw_data)
mi_analyzer = MiRandomForest(model_str="custom_rf")
pooled_result = mi_analyzer.fit(imputed_data, y="target")
print(pooled_result)

如果需要计算合并后的置信区间、p值，直接使用Barnard-Rubin自由度公式计算自由度后，基于t分布推导对应统计量即可。

内容的提问来源于stack exchange，提问作者Emily Robinson