I am using hurdle models to predict a continuous cost variable that has many exact zeros. I have fitted a hurdle model with a binomial component and a gamma component, but when I am trying to combine the two components of the model to predict average costs, I seem to be calculating the predicted probabilities incorrectly. Specifically, I used predict, type = "response"
in R and compared those predictions to the manual approach outlined in Zuur and Ieno's Beginner's Guide to Zero-Inflated Models in R on pages 128-129, and the two do not match. What am I doing wrong? Example data and code below.
library(tidyverse)
library(boot)
### DATA
disease <- sample(0:1, 75000, replace=TRUE)
age <- sample(18:88, 75000, replace=TRUE)
score <- sample(0:23, 75000, replace=TRUE)
gender <- sample(0:1, 75000, replace=TRUE)
costs <- c(sample(0:100000, (75000/2), replace=TRUE), rep(0, (75000/2)))
time <- sample(30:3287, 75000, replace=TRUE)
df <- data.frame(cbind(disease, age, score, gender, costs, time))
# create binary variable for non-zero costs
df <- df %>% mutate(costs_binary = ifelse(costs > 0, 1, 0))
### HURDLE MODEL
# gamma component
hurdle_gamma <- glm(costs ~ disease + gender + (score * age)^2 +
offset(log(time)),
data = subset(df, costs > 0),
family = Gamma(link = "log"))
# binomial component
hurdle_binomial <- glm(costs_binary ~ disease + gender + (score * age)^2 +
offset(log(time)),
data = df,
family = binomial)
# my estimate of predicted probability of use
df$prob_use <- predict(hurdle_binomial, type = "response")
# Zuur's method
gamma <- coef(hurdle_binomial)
Xb <- model.matrix(~disease + gender + (score * age)^2 +
offset(log(time)), data = df)
eta.binary <- Xb %*% gamma
pi <- exp(eta.binary) / (1 + exp(eta.binary))
# compare predicted means
mean(pi) #0.5, which is (I presume not coincidentally) the proportion of zero costs in the data
mean(df$prob_use) #much smaller than 0.5