Test for distance-based (spatial, phylogenetic or similar) autocorrelation
This function performs a Moran's I test for distance-based (spatial, phylogenetic or similar) autocorrelation on the calculated quantile residuals
testSpatialAutocorrelation(simulationOutput, x = NULL, y = NULL,
distMat = NULL, alternative = c("two.sided", "greater", "less"),
plot = T)simulationOutput |
an object of class DHARMa, either created via |
x |
the x coordinate, in the same order as the data points. Must be specified unless distMat is provided. |
y |
the y coordinate, in the same order as the data points. Must be specified unless distMat is provided. |
distMat |
optional distance matrix. If not provided, euclidean distances based on x and y will be calculated. See details for explanation |
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 Moran.I test from the package ape, based on the provided distance matrix of the data points.
There are several ways to specify this distance. If a distance matrix (distMat) is provided, calculations will be based on this distance matrix, and x,y coordinates will only used for the plotting (if provided) If distMat is not provided, the function will calculate the euclidean distances between x,y coordinates, and test Moran.I based on these distances.
Testing for spatial autocorrelation requires unique x,y values - if you have several observations per location, either use the recalculateResiduals function to aggregate residuals per location, or extract the residuals from the fitted object, and plot / test each of them independently for spatially repeated subgroups (a typical scenario would repeated spatial observation, in which case one could plot / test each time step separately for temporal autocorrelation). Note that the latter must be done by hand, outside testSpatialAutocorrelation.
Important to note for all autocorrelation tests (spatial / temporal): the autocorrelation tests are valid to check for residual autocorrelation in models that don't assume such a correlation (in this case, you can use conditional or unconditional simulations), or if there is remaining residual autocorrelation after accounting for it in a spatial/temporal model (in that case, you have to use conditional simulations), but if checking unconditional simulations from a model with an autocorrelation structure on data that corresponds to this model, they will be significant, even if the model fully accounts for this structure.
This behavior is not really a bug, but rather originates from the definition of the quantile residuals: quantile residuals are calculated independently per data point, i.e. without consideration of any correlation structure between data points that may exist in the simulations. As a result, the simulated distributions from a unconditional simulation will typically not reflect the correlation structure that is present in each single simulation, and the same is true for the subsequently calculated quantile residuals.
The bottomline here is that spatial / temporal / other autoregressive models should either be tested based on conditional simulations, or (ideally) custom tests should be used that are not based on quantile residuals, but rather compare the correlation structure in the simulated data with the correlation structure in the observed data.
Florian Hartig
testData = createData(sampleSize = 40, family = gaussian())
fittedModel <- lm(observedResponse ~ Environment1, data = testData)
res = simulateResiduals(fittedModel)
# Standard use
testSpatialAutocorrelation(res, x = testData$x, y = testData$y)
# Alternatively, one can provide a distance matrix
dM = as.matrix(dist(cbind(testData$x, testData$y)))
testSpatialAutocorrelation(res, distMat = dM)
# You could add a spatial variogram via
# library(gstat)
# dat = data.frame(res = residuals(res), x = testData$x, y = testData$y)
# coordinates(dat) = ~x+y
# vario = variogram(res~1, data = dat, alpha=c(0,45,90,135))
# plot(vario, ylim = c(-1,1))
# if there are multiple observations with the same x values,
# create first ar group with unique values for each location
# then aggregate the residuals per location, and calculate
# spatial autocorrelation on the new group
# modifying x, y, so that we have the same location per group
# just for completeness
testData$x = as.numeric(testData$group)
testData$y = as.numeric(testData$group)
# calculating x, y positions per group
groupLocations = aggregate(testData[, 6:7], list(testData$group), mean)
# calculating residuals per group
res2 = recalculateResiduals(res, group = testData$group)
# running the spatial test on grouped residuals
testSpatialAutocorrelation(res2, groupLocations$x, groupLocations$y)
# careful when using REs to account for spatially clustered (but not grouped)
# data. this originates from https://github.com/florianhartig/DHARMa/issues/81
# Assume our data is divided into clusters, where observations are close together
# but not at the same point, and we suspect that observations in clusters are
# autocorrelated
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)
# It's a good idea to use a RE to take out the cluster effects. This accounts
# for the autocorrelation within clusters
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 simulated, 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, it's also possible to simulated residuals conditional on fitted
# REs (re.form). Conditional on the fitted REs (i.e. accounting for the clusters)
# the residuals should now be indepdendent. The remaining RSA we see here is
# probably due to the RE shrinkage
res = simulateResiduals(fittedModel, re.form = NULL)
testSpatialAutocorrelation(res, x = testData$x, y = testData$y)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.