Predictions from hurdle models in R

Question

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

Answer 1

I figured it out, so for those who are interested: the model.matrix command does not include the offset, so it needed to be added manually. Once the offset is taken into account, the predictions of predict, type = "response" and the manual approach provided by Zuur coincide.