Spatial autocorrelation with nbinom1 and nbinom2

[florianhartig]

Question from a user:

I am using DHARMa 0.2.0 to look at the residuals and the spatial autocorrelation of various glmmTMB models. I know that DHARMa was only recently modified to accommodate glmmTMB, so I wanted to email you to make sure that there weren't any issues with the diagnostics I am looking at.

I am using zero-inflated, mixed effect, negative binomial models (using glmmTMB's nbinom1 and nbinom2 family). I have reason to believe that the residuals might be spatially autocorrelated, because the sampling tended to happen in spatial clusters. A random effect can be used to account for this clustered sampling pattern. I tested the spatial autocorrelation of models with and without this random effect. In nbinom1 models, spatial autocorrelation was non-significant after incorporating the random effect, but in nbinom2 models, spatial autocorrelation still appeared significant.

Is there any reason why glmmTMB models using the nbinom2 family might not work well with DHARMa?

[florianhartig]

No, I don't think that there is any reason to think that nbinom1 works better than nbinom2 - limitations of glmmTMB are described in #16 . The most serious issue is the current inability of glmmTMB to predict without REs, which could create problems such as #78, but this doesn't seem to be the issue here.

nbinom2 or nbinom1

To solve your problem, I'd say first that the decision of choosing nbinom2 or nbinom1 should not be based on the question of autocorrelation, but rather on how residual variance scales with the mean value (as these two distribution differ in how variance scales with the predicted value).

Usually, you should see which is better in the normal plot(res), but because of the glmmTMB issues, that's a bit dangerous if you have a model with random effects.

In your case, ideal would be if you would plot simulated residuals of your model against hand-calculated fixed effect predictions only for nbinom2 or nbinom1 and chose the option where you see less heteroscedasticity. Alternative (less preferred) is to do the same with the model without the REs.

spatial autocorrelation

About the spatial autocorrelation - here it becomes a bit more tricky. In many cases with spatial / temporal autocorrelation, fitting a spatial / temporal model doesn't remove the pattern from the residuals, it just accounts for the pattern in the calculation of the regression coefficients. You might have seen this in autoregressive models, for example.

Now, in your case you add a random effect to account for the spatial autocorrelation. The default in DHARMa, and the only option for glmmTMB (because of glmmTMB limitations) is to re-simulated all random effects. That means that (unlike when you simulate conditional on the REs, which would be possible for lme4) that if there was a spatial pattern in the REs, there should remain a (possibly weaker) spatial pattern in the residuals, also if you have the RE in there. You can see this if you simulate

clusters = 100
subsamples = 10
size = clusters * subsamples

testData = createData(sampleSize = size, family = gaussian(), numGroups = clusters )
testData$x  = rnorm(clusters)[testData$group] + rnorm(size, sd = 0.01)
testData$y  = rnorm(clusters)[testData$group] + rnorm(size, sd = 0.01)

library(lme4)
fittedModel <- lmer(observedResponse ~ Environment1 + (1|group), data = testData)

# DHARMa default is to re-simulted REs - this means spatial pattern remains
# because residuals are still clustered

res = simulateResiduals(fittedModel)
testSpatialAutocorrelation(res, x =  testData$x, y = testData$y)

# However, it should disappear if you just calculate an aggregate residuals per cluster
# Because at least how the data are simualted, cluster are spatially independent

res2 = recalculateResiduals(res, group = testData$group)
testSpatialAutocorrelation(res2, x =  aggregate(testData$x, list(testData$group), mean)$x, y = aggregate(testData$y, list(testData$group), mean)$x)

# For lme4, possible to simulated residuals conditional on fitted REs (re.form)
# This takes out most of the RSA - a remainder is probably due the shrinkage
# of the REs

res = simulateResiduals(fittedModel, re.form = NULL)
testSpatialAutocorrelation(res, x =  testData$x, y = testData$y)

I'll add this example to the help of the testSpatialAutocorrelation function for the next release.

Possibly this is all a bit confusing to a user, but I wouldn't view this as a bug, because simulations perform as expected, it's just about the question that you are asking.