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.modifed.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.
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.modifed.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.modifed.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.
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.