在R中使用pscl包绘制Hurdle模型计数部分及结果的方法

用predict.hurdle提取计数模型预测值绘图 #

完全可以通过predict()函数（pscl包的predict.hurdle方法）单独提取计数模型的预测结果，步骤如下：

1. 构造预测数据集
   生成包含自变量所有组合的新数据集，确保覆盖depthm的取值范围和所有Season水平：

# 假设你的原始数据框名为df
new_data <- expand.grid(
  depthm = seq(min(df$depthm, na.rm = TRUE), max(df$depthm, na.rm = TRUE), length.out = 100),
  Season = unique(df$Season)
)
# 若模型有其他控制变量，需设置为典型值（如分类变量取参考水平，连续变量取均值）
# 示例：new_data$covariate <- mean(df$covariate, na.rm = TRUE)

2. 提取计数模型的预测值
   使用type = "count"参数指定提取截断负二项模型的预测结果（默认是log尺度，需用exp()转换为原始计数尺度）：

# 假设你的hurdle模型名为hurdle_model
pred_results <- predict(hurdle_model, newdata = new_data, type = "count", se.fit = TRUE)
# 转换为原始计数尺度，并计算95%置信区间
new_data$pred_count <- exp(pred_results$fit)
new_data$lower_ci <- exp(pred_results$fit - 1.96 * pred_results$se.fit)
new_data$upper_ci <- exp(pred_results$fit + 1.96 * pred_results$se.fit)

3. 用ggplot2绘制计数模型结果
   可视化预测曲线及置信区间：

library(ggplot2)
ggplot(new_data, aes(x = depthm, y = pred_count, color = Season)) +
  # 添加置信区间带
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, fill = Season), alpha = 0.2, color = NA) +
  # 添加预测曲线
  geom_line(linewidth = 1) +
  # 自定义标签和主题
  labs(x = "水深 (m)", y = "预测Byths计数（截断负二项模型）", color = "季节", fill = "季节") +
  theme_bw()

更优的Hurdle模型可视化方案 #

除了单独绘制计数模型，还可以结合零模型的结果，展示总期望预测值（即零模型非零概率 × 计数模型条件期望），完整呈现Hurdle模型的预测效果：

# 提取零模型的非零概率（type = "zero"得到的是logit尺度的预测，需用plogis转换为概率）
new_data$non_zero_prob <- plogis(predict(hurdle_model, newdata = new_data, type = "zero"))
# 计算总期望预测值
new_data$total_exp <- new_data$non_zero_prob * new_data$pred_count

# 绘制总期望曲线
ggplot(new_data, aes(x = depthm, y = total_exp, color = Season)) +
  geom_ribbon(aes(ymin = non_zero_prob * lower_ci, ymax = non_zero_prob * upper_ci, fill = Season), alpha = 0.2, color = NA) +
  geom_line(linewidth = 1) +
  labs(x = "水深 (m)", y = "Byths总预测期望", color = "季节", fill = "季节") +
  theme_bw()

关键注意事项 #

• 确保新数据集的变量名与模型公式中的变量完全一致，交互项会自动由expand.grid生成的组合计算
• 若模型包含随机效应，需添加re.form = NA参数排除随机效应的影响（仅预测固定效应部分）
• 当depthm和Season的交互效应显著时，曲线会呈现不同季节下水深效应的差异，这正是模型要展示的核心结果

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