testTemporalAutocorrelation | R Documentation |

This function performs a standard test for temporal autocorrelation on the simulated residuals

testTemporalAutocorrelation(simulationOutput, time, alternative = c("two.sided", "greater", "less"), plot = T)

`simulationOutput` |
an object of class DHARMa, either created via |

`time` |
the time, in the same order as the data points. |

`alternative` |
a character string specifying whether the test should test if observations are "greater", "less" or "two.sided" compared to the simulated null hypothesis |

`plot` |
whether to plot output |

The function performs a Durbin-Watson test on the uniformly scaled residuals, and plots the residuals against time. The DB test was originally be designed for normal residuals. In simulations, I didn't see a problem with this setting though. The alternative is to transform the uniform residuals to normal residuals and perform the DB test on those.

Testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use recalculateResiduals function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testTemporalAutocorrelation.

Standard DHARMa simulations from models with (temporal / spatial / phylogenetic) conditional autoregressive terms will still have the respective temporal / spatial / phylogenetic correlation in the DHARMa residuals, unless the package you are using is modelling the autoregressive terms as explicit REs and is able to simulate conditional on the fitted REs. This has two consequences

If you check the residuals for such a model, they will still show significant autocorrelation, even if the model fully accounts for this structure.

Because the DHARMa residuals for such a model are not statistically independent any more, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.

There are three (non-exclusive) routes to address these issues when working with spatial / temporal / other autoregressive models:

Simulate conditional on the fitted CAR structures (see conditional simulations in the help of simulateResiduals)

Rotate simulations prior to residual calculations (see parameter rotation in simulateResiduals)

Use custom tests / plots that explicitly compare the correlation structure in the simulated data to the correlation structure in the observed data.

Florian Hartig

`testResiduals`

, `testUniformity`

, `testOutliers`

, `testDispersion`

, `testZeroInflation`

, `testGeneric`

, `testTemporalAutocorrelation`

, `testSpatialAutocorrelation`

, `testQuantiles`

, `testCategorical`

testData = createData(sampleSize = 40, family = gaussian(), randomEffectVariance = 0) fittedModel <- lm(observedResponse ~ Environment1, data = testData) res = simulateResiduals(fittedModel) # Standard use testTemporalAutocorrelation(res, time = testData$time) # If you have several observations per time step, e.g. # because you have several locations, you will have to # aggregate timeSeries1 = createData(sampleSize = 40, family = gaussian(), randomEffectVariance = 0) timeSeries1$location = 1 timeSeries2 = createData(sampleSize = 40, family = gaussian(), randomEffectVariance = 0) timeSeries2$location = 2 testData = rbind(timeSeries1, timeSeries2) fittedModel <- lm(observedResponse ~ Environment1, data = testData) res = simulateResiduals(fittedModel) # Will not work because several residuals per time # testTemporalAutocorrelation(res, time = testData$time) # aggregating residuals by time res = recalculateResiduals(res, group = testData$time) testTemporalAutocorrelation(res, time = unique(testData$time)) # testing only subgroup location 1, could do same with loc 2 res = recalculateResiduals(res, sel = testData$location == 1) testTemporalAutocorrelation(res, time = unique(testData$time)) # example to demonstrate problems with strong temporal correlations and # how to possibly remove them by rotating residuals ## Not run: set.seed(1) C <- exp(-as.matrix(dist(seq(0,50,by=.5)))) obs <- as.numeric(mvtnorm::rmvnorm(1,sigma=C)) opar <- par(mfrow = c(1,2)) image(C, main = "Specified autocorrelation (covariance)") plot(obs, type = "l", main = "Time series") par(opar) # calculate standard DHARMa residuals ## simulations from the model: x = replicate(1000, as.numeric(mvtnorm::rmvnorm(1,sigma=C))) res <- createDHARMa(x, obs, integerResponse = F) plot(res) testTemporalAutocorrelation(res, time = 1:length(res$scaledResiduals)) # calculated rotated DHARMa residuals to remove temporal correlation # this only works if the autocorrelation is homogeneous / stationary res <- createDHARMa(x, obs, integerResponse = F, rotation = C) testUniformity(res) testTemporalAutocorrelation(res, time = 1:length(res$scaledResiduals)) # the same, but with a covariance based on simulations res <- createDHARMa(x, obs, integerResponse = F, rotation = "estimated") testUniformity(res) testTemporalAutocorrelation(res, time = 1:length(res$scaledResiduals)) ## End(Not run)

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.