testTemporalAutocorrelation {DHARMa}

Test for temporal autocorrelation

Description

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

Usage

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

Arguments

simulationOutput

an object of class DHARMa, either created via simulateResiduals for supported models or by createDHARMa for simulations created outside DHARMa, or a supported model. Providing a supported model directly is discouraged, because simulation settings cannot be changed in this case.

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

Details

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.

Note

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

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

2. 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:

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

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

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

Author(s)

Florian Hartig

See Also

testResiduals, testUniformity, testOutliers, testDispersion, testZeroInflation, testGeneric, testTemporalAutocorrelation, testSpatialAutocorrelation, testQuantiles, testCategorical

Examples

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)