Error when using simulateResiduals() on a glmmTMB object

[emzepeda]

I fit this model with no warnings from glmmTMB: fitSubMale <- glmmTMB(caseN ~ medIncST + propWST + distST + natST + popDensST + (1|step_id_), data = dataSub[dataSub$Sex == "m",], family = nbinom2, start=list(theta=c(log(1e3))), map = list(theta=factor(c(NA))))

Model summary:

Next, I used simulateResiduals(fitSubMale) to run model diagnostics which results in the following error:

"Simulations from your fitted model produce infinite values. Consider if this is sensible Error in securityAssertion("Simulations from your fitted model produce NaN values. DHARMa cannot calculated residuals for this. This is nearly certainly an error of the regression package you are using", : Message from DHARMa: During the execution of a DHARMa function, some unexpected conditions occurred. Even if you didn't get an error, your results may not be reliable. Please check with the help if you use the functions as intended. If you think that the error is not on your side, I would be grateful if you could report the problem at https://github.com/florianhartig/DHARMa/issues

Context: Simulations from your fitted model produce NaN values. DHARMa cannot calculated residuals for this. This is nearly certainly an error of the regression package you are using"

I've searched quite a bit and have been unable to figure out what the issue with glmmTMB might be. I believe I've used all the functions as intended. Do you have any ideas about what might cause this error?

[florianhartig]

Hello emzepeda,

sorry, I think I forgot to respond to your issue.

a) would it be possible that you share your data with me, either here or privately (via email)

b) why are you using map = list(theta=factor(c(NA))). It seems you fix the theta of the neg binom to log(1000)? I'm also a bit confused why the dispersion parameter is 267 then.

[emzepeda]

Hi, thanks for getting back to me.

a) I can share the data over email. Is your email with the domain "biologie.uni-regensburg.de" the best one to use?

b) I was following the recommendations of Muff et al. 2020 when I fixed the intercept variance for step_id_. They suggest doing this to avoid bias in random effect intercepts resulting from the shrinkage of estimates towards the overall mean.

[florianhartig]

Hi,

yes, any of my many Uni Regensburg email aliases will work, the easiest is ur.de

However, I think I don't need the data, the reason of the problem is quite clear - if you fix the RE variance to such a large value and try simulating from this model, you will get numerical over / underflows because you have exp(rnorm(sd = very high)). I have just tried this out with the code from the Muff paper. Even if the simulations worked, with a fixed RE variance, residuals would look horrible (glmmTMB doesn't allow to condition on the fitted REs, so it will always re-simulate the REs, which obviously doesn't fit to your data at all).

I have skimmed the paper but I don't really get the details of the motivation to do this - if you set (1|stepID) with a super high variance, estimates will be effectively identical to a fixed effect on stepID. There will probably be a very slight regularisation, but I wonder if there is any difference. Wouldn't just fitting a fixed effect on stepID do the same (which would also work with DHARMa). @jmsigner, if you have time, could you may comment?

Regardless of the exact motivation, if you used REs as a trick for shrinkage estimation, models cannot be checked with DHARMa, because DHARMa assumes that the fitted model is the data-generating model (which in this case it isn't because we don't think the data comes from a RE model with variance, we use the RE as a "trick" to impose a light regularising prior).

[florianhartig]

p.s.: for future reference: if you want to set shrinkage regularisation on parameters, doing this via a Bayesian prior in brms would be an option that wouldn't confuse DHARMa. In brms, you would just fit a fixed effect on step_id_ and set a normal prior with variance e^4 on the parameters. In this case, DHARMa would simulate correctly, because it's clear that the model doesn't assume an RE with an e^4 variance.

[zepedae]

Thanks for clarifying. I'll give the Bayesian approach a try.

[jmsigner]

Hi Florian, thanks for bringing me in. The motivation for setting the variance to a large fixed value was, that we wanted to avoid estimating numerous fixed effects, since in a typical study there are many thousands or ten thousand of steps.

[florianhartig]

Ah, OK, I understand. But is this a numeric problem, or just to get the output more clean? If it's just the latter, @emzepeda could simply fit the fixed effect model and ignore all outputs on the steps?

If you want to keep the model structure like it is, with the current version of glmmTMB, which always re-simulates the fitted REs, this will not work out of the box with DHARMa. The reason is that the glmmTMB crew has not yet implemented the option to simulate conditional on the fitted REs, see https://github.com/glmmTMB/glmmTMB/issues/322

The only option would be to do the simulations by hand (i.e. predict, simulate the neg.binom) and then use createDHARMa() to read them back in DHARMa.