如何在Python中实现R ggplot2带glm误差估计的分组平滑线图

核心逻辑说明 #

原R端绘图的两个核心特征是seaborn默认regplot无法覆盖的：

• 拟合过程为带多协变量调整的GLM模型，而非简单一元OLS回归
• 置信区间由GLM预测逻辑计算生成，而非模型残差的简单区间估计
  seaborn regplot的内置逻辑不支持自定义模型拟合方式，因此无法直接匹配需求
  实现时可直接复刻R端的语法逻辑，避免逐行换算的适配误差。

分步实现代码 #

1. 依赖引入 #

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from plotnine import (
    ggplot, aes, geom_smooth, theme_bw,
    scale_y_continuous, labs
)
from matplotlib.ticker import PercentFormatter

2. 对齐R端的GLM拟合与置信区间计算 #

这一步完全对应R里glm(group~age,height,weight, data=df)的计算逻辑，按SCD分组生成拟合结果：

# 对单个SCD分组计算GLM拟合结果与置信区间
def calc_group_glm(group_data, formula):
    # 拟合与R端完全一致的GLM模型
    glm_model = smf.glm(formula=formula, data=group_data).fit()
    # 生成覆盖该组年龄范围的预测网格
    age_range = np.linspace(group_data["Age"].min(), group_data["Age"].max(), 100)
    predict_grid = pd.DataFrame({"Age": age_range})
    # 补全模型中其余协变量的取值（取组内均值，和R端predict默认处理逻辑一致）
    for cov in ["height", "weight"]:
        predict_grid[cov] = group_data[cov].mean()
    # 计算预测均值与95%置信区间
    pred_result = glm_model.get_prediction(predict_grid).summary_frame(alpha=0.05)
    return pd.concat([predict_grid, pred_result], axis=1)

# 遍历所有SCD分组，汇总所有拟合结果
all_fit_result = df.groupby("SCD", group_keys=False).apply(
    lambda g: calc_group_glm(g, "group ~ Age + height + weight")
).reset_index(drop=True)

3. 复刻R端绘图参数 #

plotnine的绘图语法和R ggplot2完全同源，参数可一一对应：

final_plot = (
    ggplot(data=all_fit_result, mapping=aes(x="Age", y="mean", color="SCD"))
    # 绘制拟合线+置信区间带，对应原R端geom_smooth的误差展示
    + geom_smooth(
        mapping=aes(ymin="ci_low", ymax="ci_high"),
        stat="identity"
    )
    # 对齐R端的theme_bw黑白主题
    + theme_bw()
    # 对齐y轴百分比格式
    + scale_y_continuous(labels=PercentFormatter(xmax=1))
    # 对齐坐标轴标签
    + labs(x="Age", y="Weight")
)

# 渲染展示图表
print(final_plot)

常见适配说明 #

• 若不想使用plotnine，可直接用matplotlib基于all_fit_result绘制：按SCD分组调用plt.plot()画拟合线，用plt.fill_between()填充置信区间范围，手动调整y轴为百分比格式即可
• 若需要调整GLM的分布族（比如二项分布对应logit连接），只需在smf.glm里传入family=smf.families.Binomial()，和R端glm(family=binomial())的计算结果完全对齐
• 计算出的GLM模型结果可直接通过glm_model.summary()输出，和R端summary(glm(...))的统计量输出逻辑一致

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