2
$\begingroup$

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
$\endgroup$

1 Answer 1

0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.