3

我有要建模的半连续数据（许多精确的零和连续的积极结果）。我从 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, 75000, replace = TRUE)
cost <- c(rep(0, 30000), rgamma(n = 37500, shape = 5000, rate = 1), 
          sample(1:1000000, 7500, replace = TRUE))
disease <- sample(0:1, 75000, replace = TRUE)
time <- sample(30:3287, 75000, 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), 0, 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 模型上做错了什么？

4

1 回答 1

2

该glmmTMB软件包可以这样做：

glmmTMB(formula, family=ziGamma(link="log"), ziformula=~1, data= ...)

应该这样做。也许VGAM也有什么？

要回答有关系数和标准误的问题：

  • 二项式系数的符号变化正是您所怀疑的（估计 0 [glmmTMB] 的概率与非零概率 [your/Zuur's code] 之间的差异）
  • 模型二项式部分的标准误差接近但不相同：使用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 完全不同）；我想知道这是否会使事情变得更难/更糟？

于 2021-01-16T00:18:44.850 回答