Creating a DHARMa object to test model residuals for a ZINB model that isn't supported by DHARMa::simulateResiduals

ZINB_DHARMa_example.r

# Downlad data
download.file("https://github.com/cran/pscl/raw/master/vignettes/DebTrivedi.rda", destfile = "DebTrivedi.rda")
load(file="DebTrivedi.rda")
df = DebTrivedi[, c(1, 6:8, 13, 15, 18)]
# this is some data about the number of outpatient appointments attended by 
# people in the US (it's from the pscl vignette)
# it doesn't really matter what this data is, just that it has a zero 
# inflated count outcome.

library(MASS)     # for glm.nb
library(pscl)     # for zeroinfl
library(DHARMa)   # for the DHARMa functions
library(emdbook)  # for rzinbinom


# try a regular GLM first
model_glm = glm.nb(ofp ~ ., data=df)
summary(model_glm)

glm_simulation = simulateResiduals(fittedModel = model_glm)
testZeroInflation(glm_simulation)
# From this plot we can see that the outcome is zero inflated and the model 
# does not account for the zero inflation.


# now try a ZINB
model = zeroinfl(ofp ~ ., data=df, dist="negbin")
summary(model)

# construct the DHARMa object manually
# 1) predictions from the two model sectors
p = predict(model, type="zero")
mu = predict(model, type="count")

# 2) simulated outcomes
sim_outcome = replicate(
  250, # this is the same number of replicates as the default for DHARMa::simulateResiduals
  rzinbinom(
    n = length(model$fitted.values),
    mu=mu,
    size=model$theta,
    zprob=p
  )
)

# 3) Bringing eveything together
simulation_model = createDHARMa(
  simulatedResponse=sim_outcome,
  observedResponse=model$y,
  fittedPredictedResponse=predict(model, type="response"),
  integerResponse=TRUE
)

# 4) Plots
plot(simulation_model)
testZeroInflation(simulation_model)
# Zero inflation test now passes :)