零调整伽马（ZAG）模型连续与二分类部分联合效应的P值获取问题

可通过两种成熟方法获取两部分联合的估计值、标准误与p值，以下是具体实现思路和可运行代码：

方法1：Delta方法（轻量化实现） #

基于两个子模型参数的渐近协方差矩阵，通过泰勒展开推导联合统计量的方差，无需重拟合模型，适合大样本场景：

# 加载所需包
library(msm)

# 提取两个模型的所有参数，合并为一个参数向量
# 参数顺序：gamma(二分类截距), beta0(伽马截距), beta1(伽马X1系数), r(伽马形状参数)
pars <- c(
  gamma = coef(h2),
  beta0 = coef(h1)[1],
  beta1 = coef(h1)[2],
  r = 1/summary(h1)$dispersion
)

# 构建分块协方差矩阵（两个子模型独立，非对角块为0）
vcov_bin <- vcov(h2)
vcov_gamma <- vcov(h1)
# 离散参数的方差近似
vcov_disp <- (summary(h1)$dispersion^2)^2 / (length(datahurdle.pos$Y) - length(coef(h1)))
vcov_total <- diag(4)
vcov_total[1,1] <- vcov_bin[1,1]
vcov_total[2:3,2:3] <- vcov_gamma
vcov_total[4,4] <- vcov_disp

# 以X1取均值时的联合均值为例计算推断结果
x1_mean <- mean(datahurdle$X1)
# 联合均值表达式：exp(gamma)/(1+exp(gamma)) * exp(beta0 + beta1*x1_mean)
est_res <- deltamethod(
  ~ (exp(gamma)/(1+exp(gamma))) * exp(beta0 + beta1*x1_mean),
  mean = pars,
  cov = vcov_total
)

# 计算Wald检验p值
z_score <- est_res[1]/est_res[2]
p_value <- 2*pnorm(-abs(z_score))

# 输出结果
cat("联合均值估计值：", est_res[1], "\n",
    "联合标准误：", est_res[2], "\n",
    "联合p值：", p_value, "\n")

如果需要对每个样本的联合拟合值做推断，只需循环修改x1_mean为对应X1取值即可。

方法2：自助法（Bootstrap，稳健性更高） #

无需依赖参数渐近正态假设，样本量充足时推断结果更可靠：

set.seed(456)
n_boot <- 1000 # 可根据算力调整抽样次数
boot_est <- replicate(n_boot, {
  # 有放回抽样
  boot_idx <- sample(1:N, size = N, replace = TRUE)
  boot_data <- datahurdle[boot_idx,]
  
  # 拟合二分类子模型
  boot_data$Y0.1 <- as.numeric(boot_data$Y >0)
  h2_boot <- try(glm(Y0.1 ~1, family=binomial, boot_data), silent = TRUE)
  if(inherits(h2_boot, "try-error")) return(NA)
  
  # 拟合伽马子模型
  boot_pos <- subset(boot_data, Y>0)
  if(nrow(boot_pos) < 10) return(NA) # 避免正样本过少导致拟合失败
  h1_boot <- try(glm(Y ~X1, family = Gamma(link = "log"), data=boot_pos), silent = TRUE)
  if(inherits(h1_boot, "try-error")) return(NA)
  
  # 计算X1取均值时的联合均值
  pi_boot <- plogis(coef(h2_boot))
  mu_gamma_boot <- exp(coef(h1_boot)[1] + coef(h1_boot)[2]*mean(boot_data$X1))
  return(pi_boot * mu_gamma_boot)
})

# 去除拟合失败的样本
boot_est <- na.omit(boot_est)

# 计算推断结果
est_boot <- mean(boot_est)
se_boot <- sd(boot_est)
p_boot <- 2*pnorm(-abs(est_boot/se_boot))

# 输出结果
cat("自助法联合均值估计值：", est_boot, "\n",
    "自助法联合标准误：", se_boot, "\n",
    "自助法联合p值：", p_boot, "\n")

补充说明 #

• 若需检验某协变量在两部分的联合显著性，可提取该变量在两个子模型的系数，通过Wald卡方检验或自助法计算联合p值
• 也可直接使用glmmTMB包的zaga族直接拟合ZAG模型，包内置了联合边际效应的推断功能，无需手动拆分计算两部分参数

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