CODEWITHSUNDEEP
rregression
Andrew Bade
Andrew Bade

Reputation: 323

Passing New Data to simulate.glmmTMB

In a previous question (Generate a predicted count distribution from a ZINB model of class glmmTMB) I asked how to generate a predicted count distribution for a zero-inflated negative binomial model of class "glmmTMB". One solution that I have since found to that question is the function simulate.glmmTMB (https://www.rdocumentation.org/packages/glmmTMB/versions/0.2.3/topics/simulate.glmmTMB). However, I want to do the simulation on test data to validate a model for predictive ability and I only see how to run simulations on the same data used to fit the model.

In the example below, how could I simulate outcomes for the newdata data frame?

library(glmmTMB)
data("bioChemists", package = "pscl")
zinb <- glmmTMB(art ~ fem + mar + kid5 + phd + ment, ziformula = ~ ., data = 
bioChemists, family = nbinom2(link = "log"))
sim_1 <- simulate(zinb) #works as expected

#make new dataframe
newdata = unique(bioChemists[,c("fem","mar","kid5","phd","ment")])
sim_2 <- simulate(zinb, newdata = newdata) #ignores newdata

Upvotes: 3

Views: 848

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226577

I think this works (could be encapsulated in a function etc.):

n <- nrow(newdata)
## construct model matrices for conditional and Z-I
##  components (.[-2] drops the LHS of the formula)
X <- model.matrix(formula(zinb)[-2],newdata)
X_zi <- model.matrix(formula(zinb,component="cond")[-2],newdata)
## extract coefficients
beta <- fixef(zinb)$cond
beta_zi <- fixef(zinb)$zi
## draw random values for each component
cond <- rnbinom(n, mu=exp(X %*% beta), size=sigma(zinb))
zi <- rbinom(1, prob=plogis(X_zi %*% beta_zi), size=1)
cond*zi

The last step is a little too clever: ifelse(zi==0,0,cond) might be clearer, or for the last three steps you could use the rzinbinom function in the emdbook package ...

In general I think simulate() methods should allow both newdata and newparams - it opens up a big range of possibilities for parametric bootstrapping, posterior simulation, etc. etc..

More compactly (and probably more robustly):

cond <- predict(zinb, newdata=newdata, type="conditional")
zi <- predict(zinb, newdata=newdata, type="zprob")
emdbook::rzinbinom(nrow(newdata),
                   mu=cond, size=sigma(zinb),
                   zprob=zi)

Upvotes: 4

Related Questions