问题描述
我有我正在尝试建模的半连续数据（许多精确的零和连续的正结果）。我从 Zuur 和 Ieno 的 R 零膨胀模型初学者指南 中学到了大量关于零质量的建模数据，它区分了零膨胀伽马模型和他们所谓的“零改变模型”。 " 伽马模型，他们将其描述为障碍模型，它结合了零点的二项式分量和正连续结果的伽马分量。我一直在探索使用 ziGamma 包中的 glmmTMB 选项，并将得到的系数与我按照 Zuur 的书（第 128-129 页）中的说明构建的障碍模型进行比较，他们确实做到了不重合。我无法理解为什么不这样做，因为我知道伽马分布不能取零值，所以我认为每个零膨胀伽马模型在技术上都是一个障碍模型。谁能为我照亮这个？在代码下方查看有关模型的更多评论。
library(tidyverse)
library(boot)
library(glmmTMB)
library(parameters)
### DATA
id <- rep(1:75000)
age <- sample(18:88,75000,replace = TRUE)
gender <- sample(0:1,replace = TRUE)
cost <- c(rep(0,30000),rgamma(n = 37500,shape = 5000,rate = 1),sample(1:1000000,7500,replace = TRUE))
disease <- sample(0:1,replace = TRUE)
time <- sample(30:3287,replace = TRUE)
df <- data.frame(cbind(id,disease,age,gender,cost,time))
# create binary variable for non-zero costs
df <- df %>% mutate(cost_binary = ifelse(cost > 0,1,0))
### HURDLE MODEL (MY VERSION)
# gamma component
hurdle_gamma <- glm(cost ~ disease + gender + age + offset(log(time)),data = subset(df,cost > 0),family = Gamma(link = "log"))
model_parameters(hurdle_gamma,exponentiate = T)
# binomial component
hurdle_binomial <- glm(cost_binary ~ disease + gender + age + time,data = df,family = "binomial")
model_parameters(hurdle_binomial,exponentiate = T)
# predicted probability of use
df$prob_use <- predict(hurdle_binomial,type = "response")
# predicted mean cost for people with any cost
df_bin <- subset(df,cost_binary == 1)
df_bin$cost_gamma <- predict(hurdle_gamma,type = "response")
# combine data frames
df2 <- left_join(df,select(df_bin,c(id,cost_gamma)),by = "id")
# replace NA with 0
df2$cost_gamma <- ifelse(is.na(df2$cost_gamma),df2$cost_gamma)
# calculate predicted cost for everyone
df2 <- df2 %>% mutate(cost_pred = prob_use * cost_gamma)
# mean predicted cost
mean(df2$cost_pred)
### glmmTMB with ziGamma
zigamma_model <- glmmTMB(cost ~ disease + gender + age + offset(log(time)),family = ziGamma(link = "log"),ziformula = ~ disease + gender + age + time,data = df)
model_parameters(zigamma_model,exponentiate = T)
df <- df %>% predict(zigamma_model,new data = df,type = "response") # doesn't work
# "no applicable method for "predict" applied to an object of class "data.frame"
我的障碍模型的 gamma 分量和 zigamma 模型的固定效应分量的系数相同，但 SE 不同，在我的实际数据中，这对我感兴趣的预测变量的显着性有重大影响。零膨胀模型的系数不同，我还注意到二项式分量中的 z 值是我的二项式模型中的负倒数。我认为这与我的二项式模型对存在概率建模（1 表示成功）和 glmmTMB 可能对缺席概率建模（0 表示成功）有关？
总而言之，谁能指出我在 glmmTMB ziGamma 模型上做错了什么？
解决方法
glmmTMB 包可以做到这一点：
glmmTMB(formula,family=ziGamma(link="log"),ziformula=~1,data= ...)
应该这样做。也许 VGAM 中也有一些东西？
回答有关系数和标准误的问题：
- 二项式系数的符号变化正是您所怀疑的（估计 0 [glmmTMB] 的概率与不为零的概率 [您/Zuur 的代码] 之间的差异）
- 模型二项式部分的标准误差接近但不相同：使用
broom.mixed::tidy，
round(1-abs(tidy(hurdle_g,component="zi")$statistic)/
abs(tidy(hurdle_binomial)$statistic),3)
## [1] 0.057 0.001 0.000 0.000 0.295
截距为 6%，年龄影响高达 30% ...
- 条件 (
cost>0) 组件的标准误差的近两倍差异绝对让我感到困惑；如果我们简单地在 glmmTMB 与 glm 中实现 Gamma/log-link，它就会成立。在这种情况下，很难知道如何检查哪个是正确的/黄金标准应该是什么。在这种情况下，我可能不信任 Wald p 值，而是尝试使用似然比检验来获取 p 值（通过drop1）。
在这种情况下，模型被严重错误指定（即成本是均匀分布的，与 Gamma 完全不同）；我想知道这是否会使事情变得更难/更糟？