How can a statistically significant difference between rasters be calculated?

Question

so after some gruelling use of MaxEnt I have four occurrence probability rasters for the four quarters of the year, each one is about 70MB+. I want to test whether they are significantly different from each other. I have been told that a Mantel test is the best option, however it seems to only be able to deal with two rasters at a time, and with such large rasters and numbers of permutations R has been failing me in terms of memory:

rsession(14001,0x7fff7be68000) malloc: * mach_vm_map(size=153281630035968) failed (error code=3) error: can't allocate region set a breakpoint in malloc_error_break to debug rsession(14001,0x7fff7be68000) malloc: mach_vm_map(size=153281630035968) failed (error code=3) error: can't allocate region * set a breakpoint in malloc_error_break to debug Error: cannot allocate vector of size 142754.6 Gb In addition: Warning message: In is.euclid(m1) : Zero distance(s)

I wonder if there is some way of doing a Mantel test or similar within QGIS, or whether anyone can help?

Answer 1

You might want to ask a more general question of the stats of this on Cross Validated.

If you look for significant differences between multiple pairs you need to correct for family wise error, to avoid this problem: https://www.xkcd.com/882/. Basically, at a 5% significance you should expect a false positive every twenty times.

For the specifics of your question, this answer is a good start: https://gis.stackexchange.com/a/132436/23160, but you'll notice it only compares two rasters. To get around this we can use Tukey's honest significant differences test (https://en.wikipedia.org/wiki/Tukey%27s_range_test), which will tell us how similar each pairing is based on an analysis of variance (ANOVA).

Brute force approach

To implement this in R, we can use the Pearson correlation answer (linked to above) as a base, but take a different statistical test:

install.packages(raster)
library(raster)

# list raster files
f = list.files("~/GIS/raster", pattern="tif", full.names=T)

# read and prep each raster
d = lapply(1:length(f), function(i){
   # read
   r = raster(f[i])
   # extract values
   data.frame(r=as.factor(i), val=getValues(r))
})

# Make single data frame
# - if you haven't already run out of ram, this might do it!
d = do.call("rbind.data.frame", d)

# ANOVA
d.aov = aov(val ~ r, d)

# Honest significant differences
result = TukeyHSD(d.aov)

# See output table
result

# Plot differences and confidence limits
plot(result)

Softly softly catchy monkey

Here we take the same approach as above, but aggregate our rasters so they are a coarser resolution and take up a manageable amount of space. I've only repeated the read and prep step, as all others are the same.

d = lapply(1:length(f), function(i){
   # read
   r = raster(f[i])
   # aggregate
   r = aggregate(r, fact=5)
   # extract values
   data.frame(r=as.factor(i), val=getValues(r))
})

Softly with a sledge hammer

You might also find that with 16 rasters the aggregation takes a while. If you're running on a Unix type system, you can easily parallelise the read and prep step with:

# note, mc.cores is the number of processor cores you want to use
library(parallel)
d = mclapply(1:length(f), mc.cores=4, function(i){
   # read
   r = raster(f[i])
   # aggregate
   r = aggregate(r, fact=5)
   # extract values
   data.frame(r=as.factor(i), val=getValues(r))
})

You could also try and consider how different each raster is to the other three as a group. However, you'd probably need to do this by cell and would only have a group of three to make a generalisation from. If you take the mean of a group of three the standard error is likely to be much bigger than the difference between the group and the subject raster.

Answer 2

This problem is inherently an NxN comparison and a Mantel or Partial Mantel are quite inappropriate here. A Mantel is a pairwise matrix correlation and is entirely dependent on distance (ie., ecological, geographic).

I have a function raster.modified.ttest in the spatialEco package that implements a moving window version of the Dutileul modified t-test, accounting for spatial autocorrelation (Dutilleul 1993; Clifford et al., 1989). This will return a raster (SpatialPixelsDataFrame) of the correlations, F-statistic, p-value and Moran's-I for x and y. The correlations can be interpreted as a Pearson product-moment correlation. Significance is derived via a Bootstrap. The statistics will be a function of the spatial scale defined in the kNN or neighbor distances so, avoid over localization with small neighborhoods/distances.

For computational tractability, I included a point sub-sampling option that allows for a defined sampling schema (random, hexagonal or systematic) and sample size. The statistics will be calculated for each point, based on the underlying rasters, and written to the associated SpatialPointsDataFrame object.

library(gstat)                                         
library(sp)
library(raster)
library(spatialEco)
data(meuse)

data(meuse.grid)

coordinates(meuse) <- ~x + y

coordinates(meuse.grid) <- ~x + y
GRID-1 log(copper):
v1 <- variogram(log(copper) ~ 1, meuse)

x1 <- fit.variogram(v1, vgm(1, "Sph", 800, 1))

G1 <- krige(zinc ~ 1, meuse, meuse.grid, x1, nmax = 30)
gridded(G1) <- TRUE

G1@data = as.data.frame(G1@data[,-2])
GRID-2 log(elev):
v2 <- variogram(log(elev) ~ 1, meuse)

x2 <- fit.variogram(v2, vgm(.1, "Sph", 1000, .6))

G2 <- krige(elev ~ 1, meuse, meuse.grid, x2, nmax = 30)
gridded(G2) <- TRUE

G2@data <- as.data.frame(G2@data[,-2])
G2@data[,1] <- G2@data[,1]
corr <- raster.modified.ttest(G1, G2)

  plot( raster(corr, 1) )

Example of using point sub-sampling with a sample size of 500 and a scale of 500m (radius).

corr.hex <- raster.modified.ttest(G1, G2, sub.sample = TRUE, d = 500, size = 100)   

You could also just grab a random sample of any two rasters and calculate the Root Mean Squared Error (RMSE).

rmse <- function(x, y) sqrt(mean((x - y)^2, na.rm=TRUE))
sr <- sample(1:nrow(G1), 100)
rmse( G1@data[sr,],  G2@data[sr,]) 

Please note that since the output is from multiple Maxent models there is some inherent stochasticity in the probabilities that cannot be accounted for without applying a simulation framework. Due to the nature of the algorithm and an arbitrary scaling of the initial probability distribution(s), you can get a different estimate, on the lower tail, every time you run the model.

Note that the function as of 2023 demands rasters in terra SpatRaster format

References

Clifford, P., S. Richardson, D. Hemon (1989), Assessing the significance of the correlation between two spatial processes. Biometrics 45:123-134.

Dutilleul, P. (1993), Modifying the t test for assessing the correlation between two spatial processes. Biometrics 49:305-314.