0

I am computing Moran’s Index from a very large distance matrix between geographic points. The time to execute the code is too high. I used the dism function in the geosphere package.

Here a reproducible example in which the Moran's index is calculated for three covariates:

## Generate a dataset
N = 100000
dat <- data.frame(long = runif(N, min = -130, max = -60), lat = runif(N, min = 16, max = 60), 
                  var1 = runif(N, min = 0, max = 40),
                  var2 = runif(N, min = 0, max = 40),
                  var3 = runif(N, min = 0, max = 40))
                  
## Generate a distance matrix
mat_dist <- as.matrix(geosphere::distm(cbind(dat[,c("long")], dat[,c("lat")])))/1000 ## distance in km
## summary(mat_dist)

## Define the data frame containing the Moran's index for each covariate
moranI <- data.frame(matrix(NA, nrow = 1, ncol = 3))
names(moranI) <- paste0("var", 1:3)
## print(moranI)

## Compute the Moran's index for each covariate
ncol_moranI <- 1
while(ncol_moranI <= dim(moranI)[2]){
  
  ## Modify the distance matrix
  mat_dist_moranI <- mat_dist
  
  ## Generate a matrix of inverse distance weights
  mat_dist_moranI <- 1/mat_dist_moranI
  
  ## Replace -Inf with 0
  diag(mat_dist_moranI) <- 0
  
  ## Calculate the Moran's index
  moranI[ncol_moranI] <- abs(ape::Moran.I(dat[, c(names(moranI)[ncol_moranI])], mat_dist_moranI, scaled = TRUE)$observed)
  
  ## Update the incremental variable
  ncol_moranI <- ncol_moranI + 1
}
## print(moranI)

I tested the distance function in the terra package, but the execution time is also too high.

Is there any way to speed up the code ? Any help would be greatly appreciated. Thanks so much.

7
  • 1
    "The time to execute the code is too high." is a bit vague - how long does it currently take for a given N and what would be acceptable? Is it purely the distance calculation that takes all the time? Or is it the Moran computation?
    – Spacedman
    Apr 26, 2023 at 7:44
  • Thanks Spacedman. It's the distance calculation that takes all the time. For N = 1000, the time execution is 1.36s and N = 10000, it's 137.63 s (a bit too long)
    – Pierre
    Apr 26, 2023 at 15:46
  • Is 10,000 the most you need? Because your example has N = 100000 - ten times that...
    – Spacedman
    Apr 26, 2023 at 15:59
  • I need to run the code using a dataframe with 100,000 rows.
    – Pierre
    Apr 26, 2023 at 16:07
  • If this works for you I think we can close this Q as a duplicate! gis.stackexchange.com/questions/291304/… from github.com/mcooper/moranfast although it doesn't do geospatial distance...
    – Spacedman
    Apr 26, 2023 at 20:30

1 Answer 1

2

For 100,000 points your distance matrix is going to be 10,000,000,000 elements, and on typical hardware it will be size rather than time that will stop it. Any code that requires a pre-computed distance matrix for its weights will fail unless you have around 100Gb of RAM on your system. I don't know of any Moran computation code that takes a functional approach, where instead of passing a matrix you pass a function, and when the code needs a weight it calls the function. This is much lighter in memory usage but potentially slower with a function call overhead and loses any parallel advantage of doing loads of distance calculations. But if memory is a constraint and you really need Moran's I for 100,000 points with a full distance range consideration, then that's the way to do it.

But... do you really need to do that? You're doing a hypothesis test which results in either rejection of the null hypothesis or inability to reject the null hypothesis. The hypotheses in this case are spatially invariant - if you can reject the null in one corner of the data, then that's sufficient to reject it everywhere. Suppose my null hypothesis is "there's no gold on this planet" and I dig a small field and find some gold, I can reject that hypothesis. If there's spatial autocorrelation in a subset of your data, then that's the hypothesis rejected for all the data.

Which leans towards an approach - subset your data into small regions and run multiple Moran I tests, do the appropriate multiple-hypothesis testing correction and that should do. You could even just take a spatially at-random subset of your data - 5,000 points runs in pretty quick time on my desktop - and do the test on that, and repeat with resampling.

Essentially I think you're doing the wrong thing here, and there's probably some other approach to answer whatever the underlying question is about your data.

4
  • Thanks Spacedman. I need to select covariates based on Moran I. What is the best solution? Solution 1: subset the data and run the Moran I test for each covariate (e.g. divide into 1000 subsets of 1000 rows and take the mean of Moran I indexes across the 1000 subsets) or solution 2: take the mean across N (e.g. 5000) random subsets sampled with replacement) ? Thanks so much for your help.
    – Pierre
    Apr 26, 2023 at 16:04
  • Ah so you're not using Moran's I as a hypothesis test, you're planning to rank the covariates by Moran's I and use that for covariate selection in another stage? Moran's I should be invariant under random thinning and spatial subsetting but I can't get my brain around the relative uncertainties of those as estimators of the population Moran I... Another approach may be to do some (sampled) variograms to assess the range of correlation in your data then compute Moran I on a sparse matrix cutoff at some distance...
    – Spacedman
    Apr 26, 2023 at 16:12
  • Thanks Spacedman, Yes, exactly. But the variograms should be defined from subsets of my dataframe in my case. That’s right ?
    – Pierre
    Apr 26, 2023 at 16:33
  • 1
    Yes. Explore a bit. Compute Moran's I on spatially random subsets of varying sizes and see if the sample variance decreases as subset size increases, you might then extrapolate to N=100,000 to get an estimate with variance for the whole data. If the extrapolated variances are small enough then you can rank the covariates with confidence. There may even be literature on this...
    – Spacedman
    Apr 26, 2023 at 20:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.