我有半连续数据(许多精确的零和连续的正结果),我正在尝试建模。我从 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"
我的障碍模型的伽马分量和 zigamma 模型的固定效应分量的系数是相同的,但 SE 不同,这在我的实际数据中对我感兴趣的预测变量的显着性具有重大影响。零膨胀模型的系数不同,我还注意到二项式分量中的 z 值是二项式模型中 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
drop1
)。
glmmTMB
的较新版本还可以使用
Tweedie distribution作为响应的条件分布;该分布允许零和连续正值的组合。