使用ggpmisc的stat_poly_eq提取logit模型预测方程时遇报错

问题描述 #

在基于真实数据绘制logit模型预测曲线并提取方程时，使用ggpmisc::stat_poly_eq函数，不添加use_label参数时能正常运行（仅显示R²=1），但添加use_label(labels=c("eq"))以显示方程时触发报错：

Error in `check_subclass()`:
! `x` must be either a string or a <Geom> object, not a <uneval> object.

报错来自ggplot2的check_subclass函数。

报错原因 #

1. 语法错误：stat_poly_eq的use_label参数接受字符串或命名向量，而非函数调用。错误写法use_label(labels=c("eq"))会传入未求值的表达式（uneval对象），触发类型检查报错。
2. 逻辑偏差：stat_poly_line和stat_poly_eq是用于拟合多项式回归的统计层，但代码中newdat已经是logit模型的非线性预测结果（S型曲线），用这两个函数会重新拟合多项式，既不符合logit模型的预测逻辑，也无法得到目标方程。

解决方案 #

步骤1：修正语法错误（仅解决报错，不匹配logit方程） #

将错误的函数调用改为参数赋值：

stat_poly_eq(data=newdat, mapping=aes(x=newdat[[1]], y=newdat[[2]]),
             use_label = "eq")

注意：此方法仅解决报错，得到的是多项式拟合方程，不是logit模型的方程

步骤2：正确绘制logit预测曲线并显示logit方程 #

方法一：手动提取logit方程并添加

直接从模型中提取系数构建方程文本，用annotate添加到图中：

plot_logit_lab = function(log_mod){
  
  mod_frame = model.frame(log_mod)
  var_names = names(mod_frame)
  # 提取logit模型系数并保留三位小数
  coefs = round(coef(log_mod), 3)
  intercept = coefs[1]
  slope = coefs[2]
  # 构建概率形式的logit方程文本
  eq_text = paste0("P(y=1) = 1/(1+e^(-(", intercept, " + ", slope, "x)))")
  
  newdat = setNames(data.frame(seq(min(mod_frame[[2]]), max(mod_frame[[2]]), 
                                   len=100)),
                    var_names[2])

  newdat[var_names[1]] = predict(log_mod, newdata = newdat, type="response")
  
  the_plot <- ggplot() +
    geom_point(data=mod_frame, aes(x=.data[[var_names[2]]], y=.data[[var_names[1]]])) +   
    # 用geom_line绘制logit预测曲线，替代stat_poly_line
    geom_line(data=newdat, aes(x=.data[[var_names[2]]], y=.data[[var_names[1]]]), color="red") +
    # 在图中添加logit方程
    annotate("text", x=max(mod_frame[[2]])*0.8, y=max(mod_frame[[1]])*0.8, 
             label=eq_text, hjust=0, size=4)
  
  return(the_plot)
}

方法二：用ggpmisc::stat_fit_glm直接拟合并显示方程

无需手动计算预测值，直接用stat_fit_glm拟合logit模型并输出方程：

plot_logit_lab = function(log_mod){
  
  mod_frame = model.frame(log_mod)
  var_names = names(mod_frame)
  
  the_plot <- ggplot(mod_frame, aes(x=.data[[var_names[2]]], y=.data[[var_names[1]]])) +
    geom_point() +
    # 拟合logit模型并绘制曲线
    stat_fit_glm(method = "glm", method.args = list(family = binomial(link = "logit")),
                 geom = "line", aes(color="Logit Fit")) +
    # 解析并显示logit方程
    stat_fit_glm(method = "glm", method.args = list(family = binomial(link = "logit")),
                 geom = "text", aes(label = after_stat(eq.label)),
                 parse = TRUE, hjust=0, vjust=1) +
    scale_color_manual(values="red") +
    theme(legend.position="none")
  
  return(the_plot)
}

测试验证 #

运行原有测试代码，调用修正后的函数即可得到正确的logit预测曲线和方程：

car <- cars %>%
  mutate(d = case_when(dist < 50 ~ "less",
                        dist >= 50 ~ "more"),
         s = speed)

car$d <- factor(car$d, levels=c("less","more"))

car_sum <- car %>%
  count(s, d) %>%
  group_by(s) %>%
  mutate(s_n = sum(n),
         d_prob = n/sum(n))

car_logit <- glm(formula = d_prob ~ s,
                weight = s_n,
                family = binomial(link = "logit"),
                data = car_sum)

car_plot <- plot_logit_lab(car_logit)

car_plot

内容的提问来源于stack exchange，提问作者atlas-wants