glm()预测值与ggplot2 geom_smooth()绘图结果不匹配问题求解

问题原因说明 #

• 问题1报错原因：glm()的predict()方法是输入自变量预测因变量，你传入的newdata仅包含因变量C3_proportion，缺少模型必须的自变量elevation，因此触发变量未找到的报错。且预测函数本身不支持直接传入因变量反推自变量，需要额外编写逻辑实现。
• 问题2预测值和绘图结果不一致原因：二项分布GLM默认predict()返回的是链接函数（logit）尺度的对数几率值，而绘图中的拟合线是*响应尺度（0-1比例）*的结果，指定type = "response"参数即可得到和图一致的比例值。

解决步骤 #

1. 获得与绘图一致的预测结果 #

在predict()调用中添加type = "response"参数：

predict(fit, newdata = data.frame(elevation=c(1000,2000,3000)), type = "response")
#        1         2         3 
#0.2919801 0.5184085 0.7375279

输出结果与图中拟合线的数值完全匹配。

2. 反推指定C3属占比对应的海拔 #

首先需要注意：y=1经过logit转换后结果为无穷大，你的样本最大C3属占比仅为0.875，模型无法外推得到y=1对应的海拔，我们以y=0.25、0.5为例实现反推，可通过uniroot()函数求解方程获得结果，核心逻辑是找到使得模型预测值等于目标y的海拔值：

# 定义反推函数
inv_predict <- function(target_y, model = fit, x_range = range(df$elevation)) {
  # 目标函数：计算给定海拔下预测值与目标y的差值
  diff_fun <- function(x) {
    pred <- predict(model, newdata = data.frame(elevation = x), type = "response")
    return(pred - target_y)
  }
  # 求解区间内的根
  res <- uniroot(diff_fun, interval = x_range)
  return(res$root)
}

# 执行反推
target_y <- c(0.25, 0.5)
sapply(target_y, inv_predict)
# [1] 1180.749 1968.083

得到的结果即为目标C3属占比对应的海拔值，其中0.25对应约1181m，0.5对应约1968m，与图中拟合线趋势完全一致。

内容的提问来源于stack exchange，提问作者J. D. Vidal