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