如何为pscl包构建的hurdle模型预测结果生成置信区间

pscl包的hurdle模型为两部分混合模型：零模型（二项分布，拟合物种出现概率）、计数模型（截断泊松/负二项分布，拟合出现后的条件丰度），可信区间需要分别对两部分参数做不确定性传递后合并结果，最常用的实现方法为参数Bootstrap法，具体步骤如下：

实现步骤 #

• 拟合初始hurdle模型，明确两部分的协变量结构：hurdle()函数中公式竖线左侧为计数部分协变量，右侧为零部分（出现概率）协变量
• 从拟合模型的系数多元正态分布中抽样参数，抽样次数建议不低于1000次，保证分位点估计稳定性
• 基于每次抽样的系数，分别计算新数据的出现概率、条件丰度，相乘得到该次抽样的总期望丰度
• 对所有抽样结果的分位点取2.5%、97.5%，即为95%可信区间

示例代码 #

# 加载依赖包
library(pscl)
library(MASS)

# 1. 拟合hurdle模型，可根据需求替换为poisson分布
mod <- hurdle(abundance ~ temp + elevation | temp + precipitation, 
              data = 野外采样数据集, dist = "negbin")

# 2. 准备非采样区域的环境变量数据，列名需和模型协变量完全匹配
new_dat <- data.frame(temp = 非采样区域温度值, 
                      elevation = 非采样区域海拔值, 
                      precipitation = 非采样区域降水值)

# 3. 生成模型系数抽样样本
抽样次数 <- 1000
coef_samps <- MASS::mvrnorm(n = 抽样次数, mu = coef(mod), Sigma = vcov(mod))

# 4. 拆分两部分的系数抽样结果
coef_zero <- coef_samps[, grep("zero_", colnames(coef_samps))]
coef_count <- coef_samps[, grep("count_", colnames(coef_samps))]

# 5. 生成新数据的两部分设计矩阵
X_zero <- model.matrix(~ temp + precipitation, data = new_dat)
X_count <- model.matrix(~ temp + elevation, data = new_dat)

# 6. 逐次计算预测值
pred_mat <- matrix(NA, nrow = nrow(new_dat), ncol = 抽样次数)
for(i in 1:抽样次数){
  # 计算物种出现概率
  p_occur <- plogis(X_zero %*% coef_zero[i,])
  # 计算出现后的条件丰度
  cond_abund <- exp(X_count %*% coef_count[i,])
  # 计算总期望丰度
  pred_mat[,i] <- p_occur * cond_abund
}

# 7. 提取点估计和95%可信区间
new_dat$预测丰度均值 <- apply(pred_mat, 1, mean)
new_dat$预测丰度下限 <- apply(pred_mat, 1, quantile, 0.025, na.rm = T)
new_dat$预测丰度上限 <- apply(pred_mat, 1, quantile, 0.975, na.rm = T)

# 若仅需出现概率的可信区间，提取plogis(X_zero %*% coef_zero[i,])的结果计算分位点即可
# 若需全区域种群总估计的可信区间，将每次抽样的所有新数据预测值求和后再计算分位点即可

注意事项 #

• 如果你的模型零部分和计数部分协变量结构完全一致，可共用同一个设计矩阵，无需单独生成
• 若需要更高精度的区间估计，可将抽样次数调整为5000次，仅需修改抽样次数参数即可
• 如果要对预测结果做空间绘图，直接将new_dat对应的值匹配到空间栅格的像元位置即可

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