含分类因子的二项glm交互图预测值计算及绘图异常问题

问题根源 #

你遇到的连线垂直串扰问题，来自预测集生成逻辑的匹配偏差：
之前用expand.grid()无差别纳入所有解释变量的全水平组合，把需要展示的核心交互效应（月相Phase × 调查时长Duration），和风向、海况等协变量水平带来的预测值波动混在了一起，最终导致同一调查时长取值对应多个不同预测值，ggplot连线时自然出现杂乱交叉。

注：分类协变量无法直接取均值做预测简化是正确的，但这不等于必须遍历所有协变量组合才能做效应图——固定协变量在代表性水平、或对协变量组合做加权平均，是学术领域制作效应图的通用标准逻辑。

调整步骤 #

• 先明确可视化的核心变量框架：X轴为连续变量Duration，分组（对应线条颜色）为分类变量Phase，这两类变量必须保留全维度网格：Duration序列覆盖15-45分钟的等距范围，每个Duration取值对应所有月相水平的预测结果，这部分是需要展示的交互效应本身。
• 对非核心的控制协变量（风向、海况）调整生成规则：放弃原有的全组合遍历，根据调查样本的实际分布固定在代表性水平：
  • 无序分类变量（如风向）固定为样本量占比最高的水平
  • 有序分类变量（如海况）固定为调查样本的中位水平
• 重新计算预测值并转换为概率尺度，此时每个Duration下每个Phase仅对应1个唯一预测值，不存在同横坐标多纵坐标的问题，连线串扰问题会直接消除。

可复用R实现代码 #

以固定风向为样本最常见的东南风、海况为样本中位2级为例，调整后的代码如下：

# 提取模型内置的月相分类水平
phase_levels <- levels(model$data$Phase)
# 生成15-45分钟的等距调查时长序列
duration_seq <- seq(from = 15, to = 45, by = 1)
# 构建预测数据集，非核心协变量固定在代表性水平
pred_data <- expand.grid(
  Duration = duration_seq,
  Phase = phase_levels,
  WindDir = "SE", # 替换为你样本中占比最高的风向水平
  Seastate = "2" # 替换为你样本中的中位海况水平
)
# 计算logit尺度预测值与标准误，转换为概率尺度
pred_out <- predict(model, newdata = pred_data, type = "link", se.fit = TRUE)
pred_data$occur_prob <- plogis(pred_out$fit)
pred_data$ci95_low <- plogis(pred_out$fit - 1.96 * pred_out$se.fit)
pred_data$ci95_high <- plogis(pred_out$fit + 1.96 * pred_out$se.fit)

调整完成后，直接用ggplot(pred_data, aes(x = Duration, y = occur_prob, color = Phase)) + geom_line() + geom_ribbon(aes(ymin = ci95_low, ymax = ci95_high, fill = Phase), alpha = 0.2)即可绘制出无串扰的交互效应图。

进阶可选方案 #

如果固定单一协变量水平的代表性不足，可采用加权平均的边际效应计算逻辑：

• 统计样本中每个风向、海况水平的出现占比，计算所有协变量组合的对应权重
• 遍历所有协变量组合，计算各组合下在不同Duration、Phase取值的预测概率
• 对同一Duration、Phase取值下的所有协变量组合的预测概率，按样本占比权重做加权平均，得到最终的边际预测值
  这种方式的结果不依赖某一个单一协变量水平，普适性更强，同样不会出现预测值串扰的问题。

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