DHARMa's temporal autocorrelation plot differs from acf plot

mlkirchner commented 1 year ago

Hi!

Thank you so much for this package! I use it all the time! Recently, I was analyzing a temporally autocorrelated dataset with a linear mixed model with glmmTMB and noticed that when I included the first order autocorrelation term, the DHARMA testTemporalAutocorrelation plot didn't seem to change much from the model without the term. Meanwhile if I plot the residuals with acf, the plots do differ. Why would this be happening? Are the residuals calculated by recalculateResiduals not accounting for the autocorrelation term somehow? Or am I missing an argument somewhere? See below for example:

First, the model without any correction for autocorrelation:

# Model w/o autocorrelation term
mean_lmm <- glmmTMB(mean_day ~ stratum + 
                      (1|tree_code) + 
                      (1|id) + (1|date),
                     data = tc_summ)

# DHARMA temporal autocorrelation results
simulationOutput <- simulateResiduals(fittedModel = mean_lmm) 
# Test for temporal autocorrelation
residrecal <- recalculateResiduals(simulationOutput, group = tc_summ$date)
testTemporalAutocorrelation(simulationOutput = residrecal, time = unique(tc_summ$date))

The results clearly show that the residuals are autocorrelated and agree with the results of acf.

# acf temporal autocorrelation results
acf(residuals(mean_lmm))

Now, the model that includes a first-order autocorrelation term:

# Model w/ AR1 term 
mean_lmm_ac <- glmmTMB(mean_day ~ stratum  + 
                         (1|tree_code)  + 
                         (1 | id) + (1 | date) + 
                         ar1(date + 0 | id), # First order autocorrelation 
                       data = tc_summ)

# DHARMA temporal autocorrelation results 
simulationOutput <- simulateResiduals(fittedModel = mean_lmm_ac) 
# Test for temporal autocorrelation
residrecal <- recalculateResiduals(simulationOutput, group = tc_summ$date)
testTemporalAutocorrelation(simulationOutput = residrecal, time = unique(tc_summ$date))

These results are nearly the same as the model without the AR1 term, and differ from the results of acf.

acf(residuals(mean_lmm_ac))

The models differ in their summaries, estimates, and other DHARMA diagnostic plots, so it's not that the models are the same, but I can't figure out why DHARMA results are differing from acf. I know that with lme models, if you don't specify the residuals(mod, type = "normalized"), the call to acf won't take an autocorrelation term into account. Is something similar happening here?

Thanks again! Michelle

florianhartig commented 1 year ago

Hello,

this is expected (however, admittedly, not ideal behaviour). The reason behind this is the following:

DHARMa residuals are not expected to change when including REs, in particular REs with covariance structures, because glmmTMB does not (yet) allow to condition simulations and predictions on the fitted REs. Basically DHARMa residuals are calculated around your fixed effect structure (though they re-simulate the REs). See https://github.com/florianhartig/DHARMa/issues/16 and just opened issue https://github.com/glmmTMB/glmmTMB/issues/888.
glmmTMB residuals are calculated conditional on fitted REs, so in this case you see the residual correlation minus the RE structure.

How to fix this:

First of all, the fact that your autocorrelation structure disappears with glmmTMB residuals (either unstructured or with a structured RE) means that the REs are temporally autocorrelated, so you should fit an AR1 structure.

Ifs you test a glmmTMB model with covariance in the REs (such as AR1) in DHARMa, you can currently either ignore the autocorrelation, or use the rotation argument, see https://github.com/florianhartig/DHARMa/issues/364. I hope that the glmmTMB developers will include an option to condition simulations on the fitted REs in the future, which should produce similar residuals as you produce with act(residuals(model)).

mlkirchner commented 1 year ago

Hi Florian,

Thanks for the quick and helpful response! I didn't realize that the DHARMa residuals for glmmTMB are not conditioned on the RE structure, but then it makes complete sense that the autocorrelation would still be present.

I'm fine with ignoring the autocorrelation in DHARMa until glmmTMB #888 is addressed.

But I did try the rotation argument and am still not getting the same results as acf. That's also to be expected, right? Because as you said in #301, "The problem, however, is that this will likely only work perfectly for symmetric linear homogenous cdfc, such as the multivariate normal." So, hypothetically, since I'm seeing a decrease in the autocorrelation structure by including the rotation argument, that's indicative that the AR1 structure is effective, even though the autocorrelation isn't completely solved by DHARMa's calculations?

This is mostly out of curiosity at this point, because as I said, I'm fine with continuing to use glmmTMB's residuals and acf to detect autocorrelation structure.

Thanks again! Michelle

florianhartig commented 1 year ago

Hello Michelle,

you tried this on your AR1 model, right? My comments referred to real GLMMs - in your case, you essentially fit an LMM, so the rotation should be exact, provided your nSim is large enough.

By exact, I mean that the autocorrelation part of the AR1 model is correctly removed, and if you see autocorrelation remaining in the DHARMa residuals with rotation, you should probably interpret this as some kind of signal remaining.

There are still differences between the residuals, because DHARMa does not condition on your other REs, so the residuals will differ in their exact numeric values. It could be that some of the other REs absorb autocorrelation. In particular, in your model, you have both

(1 | date) ar1(date + 0 | id)

The (1 | date) will likely absorb autocorrelation that does not fit the form of the AR1, which would show up in DHARMa but not in the glmmTMB plots. You could test this by calculating the ~~acdf~~ acf on the fitted REs of (1 | date)

Cheers, Florian

mlkirchner commented 1 year ago

Hi Florian,

Sorry for being uninitiated, but what do you mean by acdf? I'm not familiar with that term.

I suspect you're correct that the random intercept of date is absorbing some of the autocorrelation and causing the differences between the two.

I really appreciate your help on this!

Michelle

On Tue, Jan 31, 2023 at 11:15 AM Florian Hartig @.***> wrote:

Hello Michelle,

you tried this on your AR1 model, right? My comments referred to real GLMMs - in your case, you essentially fit an LMM, so the rotation should be exact, provided your nSim is large enough.

By exact, I mean that the autocorrelation part of the AR1 model is correctly removed, and if you see autocorrelation remaining in the DHARMa residuals with rotation, you should probably interpret this as some kind of signal remaining.

There are still differences between the residuals, because DHARMa does not condition on your other REs, so the residuals will differ in their exact numeric values. It could be that some of the other REs absorb autocorrelation. In particular, in your model, you have both

(1 | date) ar1(date + 0 | id)

The (1 | date) will likely absorb autocorrelation that does not fit the form of the AR1, which would show up in DHARMa but not in the glmmTMB plots. You could test this by calculating the acdf on the fitted REs of (1 | date)

Cheers, Florian

— Reply to this email directly, view it on GitHub https://github.com/florianhartig/DHARMa/issues/368#issuecomment-1410668743, or unsubscribe https://github.com/notifications/unsubscribe-auth/A5SADV4PHCTXBYI4MFL7X5DWVE3AZANCNFSM6AAAAAAULT3KEE . You are receiving this because you authored the thread.Message ID: @.***>

florianhartig commented 1 year ago

typo, I meant acf

florianhartig commented 1 year ago

OK, I will consider this closed for the moment, in case you have further questions feel free to re-open the issue!

florianhartig / DHARMa

DHARMa's temporal autocorrelation plot differs from acf plot #368