如何用Python（优先statsmodels）获取混合线性模型预测置信区间？

核心解决方案：利用statsmodels的get_prediction()方法 #

statsmodels的线性混合模型（MixedLM）内置了get_prediction()方法，可直接计算多自变量场景下的预测值、标准误及置信区间，无需手动推导公式。以下是完整实现步骤：

步骤1：拟合线性混合模型 #

先加载数据并拟合模型，以经典的sleepstudy数据集为例（类似R中的同名数据）：

import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 加载示例数据
data = sm.datasets.get_rdataset("sleepstudy", "lme4").data

# 拟合线性混合模型（反应时~天数，按受试者分组）
model = smf.mixedlm("Reaction ~ Days", data=data, groups=data["Subject"])
result = model.fit()

步骤2：生成预测数据网格 #

针对要可视化的自变量生成连续取值的网格，其他自变量固定为均值或参考水平（多变量场景下需注意）：

# 生成预测网格：仅变化Days，固定Subject为任意值（边际预测无需随机效应）
pred_grid = pd.DataFrame({
    "Days": np.linspace(data["Days"].min(), data["Days"].max(), 50),
    "Subject": data["Subject"].iloc[0]
})

步骤3：计算预测值与置信区间 #

使用get_prediction()获取预测结果，re_formula="0"指定计算边际预测（population-level，即不考虑随机效应的群体平均预测）；若需条件预测（包含随机效应的个体水平预测），可去掉该参数：

# 计算边际预测及置信区间
pred_marginal = result.get_prediction(pred_grid, re_formula="0")
pred_df = pred_marginal.summary_frame(alpha=0.05)  # alpha=0.05对应95%置信区间

summary_frame()直接返回包含mean（预测值）、mean_ci_lower/mean_ci_upper（置信区间上下限）的DataFrame，无需手动计算。

步骤4：可视化结果 #

结合matplotlib/seaborn绘制原始数据、预测线及置信区间：

plt.figure(figsize=(10,6))
# 绘制原始数据散点
sns.scatterplot(x="Days", y="Reaction", data=data, alpha=0.3, label="Raw Data")
# 绘制预测线
plt.plot(pred_grid["Days"], pred_df["mean"], color="#e74c3c", label="Population-Level Prediction")
# 绘制置信区间填充
plt.fill_between(pred_grid["Days"], 
                 pred_df["mean_ci_lower"], 
                 pred_df["mean_ci_upper"], 
                 color="#e74c3c", alpha=0.2, label="95% Confidence Interval")

plt.xlabel("Days of Sleep Deprivation")
plt.ylabel("Reaction Time (ms)")
plt.legend()
plt.title("Linear Mixed Model Prediction")
plt.show()

多自变量场景示例 #

若模型包含多个自变量（如添加性别变量），只需在生成预测网格时固定其他变量的取值：

# 添加模拟性别变量
data["Gender"] = np.random.choice(["Male", "Female"], size=len(data))
# 拟合多自变量混合模型
model_multi = smf.mixedlm("Reaction ~ Days + Gender", data=data, groups=data["Subject"])
result_multi = model_multi.fit()

# 生成预测网格：固定Gender为Male，变化Days
pred_grid_male = pd.DataFrame({
    "Days": np.linspace(data["Days"].min(), data["Days"].max(), 50),
    "Gender": "Male",
    "Subject": data["Subject"].iloc[0]
})
pred_male = result_multi.get_prediction(pred_grid_male, re_formula="0").summary_frame()

# 同理生成Female的预测网格与结果
pred_grid_female = pred_grid_male.copy()
pred_grid_female["Gender"] = "Female"
pred_female = result_multi.get_prediction(pred_grid_female, re_formula="0").summary_frame()

# 绘图
plt.figure(figsize=(10,6))
sns.scatterplot(x="Days", y="Reaction", hue="Gender", data=data, alpha=0.3)
plt.plot(pred_grid_male["Days"], pred_male["mean"], color="#3498db", label="Prediction (Male)")
plt.fill_between(pred_grid_male["Days"], pred_male["mean_ci_lower"], pred_male["mean_ci_upper"], color="#3498db", alpha=0.2)
plt.plot(pred_grid_female["Days"], pred_female["mean"], color="#f39c12", label="Prediction (Female)")
plt.fill_between(pred_grid_female["Days"], pred_female["mean_ci_lower"], pred_female["mean_ci_upper"], color="#f39c12", alpha=0.2)

plt.xlabel("Days of Sleep Deprivation")
plt.ylabel("Reaction Time (ms)")
plt.legend()
plt.title("Prediction with Multiple Predictors")
plt.show()

关键说明 #

• 无需手动计算标准误：get_prediction()已自动处理多自变量的协方差矩阵计算，覆盖复杂模型场景。
• 边际/条件预测切换：re_formula="0"生成群体平均预测，去掉该参数则生成包含随机效应的个体预测。
• 可视化灵活性：可根据需求调整绘图样式，完全实现R中plot_model()的类似效果。

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