I'm struggling to evaluate the spatial autocorrelation of several databases. These databases consist of coordinates and different environmental information associated with those coordinates. Strangely, I can get very low Morans's I values (almost zero), suggesting no clustering, while ps are also very low, suggesting that clustering cannot be rejected... Trying to understand the situation, I compared (following different tutorials published online) the performance of three different libraries in the R language. Firstly, I created an ad hoc database that, I know beforehand, is not clustered:

rm(list = ls())
# Let's create some data
dta <- data.frame(LON = runif(50, 30, 60), LAT = runif(50, 30, 60), 
                            X = round(runif(50, 1, 10),1))
plot(dta$LON, dta$LAT, type = "n")
text(dta$LON, dta$LAT, labels = dta$X)

Plot of the random data

# Let's calculate the distances among points
dta.dists <- pointDistance(dta[, c("LON", "LAT")], dta[, c("LON", "LAT")], lonlat=TRUE, allpairs = T)
diag(dta.dists) <- NA
d1 <- 0
d2 <- max(dta.dists, na.rm = T)
dta.dists.inv <- 1/dta.dists
diag(dta.dists.inv) <- 0

Then, I calculated the Moran's I index using the package ape

Moran.I(dta$X, dta.dists.inv)

Which yields:

[1] 0.02912351

[1] -0.02040816

[1] 0.03655231

[1] 0.1753889

Then, I continued using the package spdep

coo <- coordinates(cbind(dta$LON, dta$LAT))
nb  <-  dnearneigh(coo, d1, d2)
moran.test(dta$X, nb2listw(nb, style="C"))

Obtaining the following results:

        Moran I test under randomisation

data:  dta$X  
weights: nb2listw(nb, style = "C")    

Moran I statistic standard deviate = -1.4754e-09, p-value = 0.5
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
    -2.040816e-02     -2.040816e-02      2.211772e-17 

Finally, I applied the elsa package:

coordinates(dta) <- ~LON +LAT
elsa::moran(dta[,1], d1, d2)


[1] -0.02040816

That is: three packages and two different values.

Which one is the correct one?

  • 1
    I would start by specifying for the sdped and ape Morans whether you want a one or two sided alternative hypotheses. As described here, rdrr.io/rforge/spdep/man/moran.test.html spdep defaults to one sided. The docs here rdrr.io/cran/ape/man/MoranI.html show ape defaults to a two sided alternative hypothesis.
    – jgm_GIS
    Sep 22 at 14:14
  • 3
    The discrepancy in the observed Moran's I values stems from the different weights used in the moranI() function and both the moran.test() and moran() functions. In the former, you are adopting inverse weights for each neighbor. In the latter two examples you are assigning equal weights to all neighbors (1/n). In the spdep package, you might want to look at the nb2listwdist() function using idw style and alpha =1 to try and mimic the inverse weighting strategy used with ape. You might also need to transform the lat/lon values if spdep does not have a geodesic option in its weighting functions.
    – MannyG
    Sep 22 at 14:46

2 Answers 2


To answer your question in brief, I would use spdep. It is tested against geoda and pysal and gives the same answer as both of those tools.

The errors you were running into were likely caused by using different weights matrixes. Your inverse distance weights are not the same as the C style weights created by spdep. Lets compare the first row of dta.dists.inv to the weights created by nb2listw(nb, style = "C")

# all of the above is the same 
coo <- coordinates(cbind(dta$LON, dta$LAT))
nb  <-  dnearneigh(coo, d1, d2)
listw <- nb2listw(nb, style="C")
moran.test(dta$X, nb2listw(nb, style="C"))

#> [1] 0.02040816 0.02040816 0.02040816 0.02040816 0.02040816 0.02040816 0.02040816
 [8] 0.02040816 0.02040816 0.02040816

#> [1] 0.000000e+00 6.310735e-07 4.059586e-07 5.128319e-07 4.767474e-07 5.054279e-07
 [7] 7.016061e-07 4.394086e-07 6.369202e-07 3.639265e-07

So that alone is enough to ensure different results.

Let's take a better example using the famous guerry dataset. We can compare {ape} and {spdep}. elsa is a library I've never heard of and does not support weight matrixes or anything other than distances—so I wouldn't use it.


# spdep
df <- Guerry::gfrance85
x <- df$Crime_pers
nb <- poly2nb(df)
listw <- nb2listw(nb)

moran.test(x, listw)
#>  Moran I test under randomisation
#> data:  x  
#> weights: listw    
#> Moran I statistic standard deviate = 6.0484, p-value = 7.316e-10
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic       Expectation          Variance 
#>       0.411459718      -0.011904762       0.004899501

# ape
ape::Moran.I(x, listw2mat(listw))
#> Registered S3 method overwritten by 'ape':
#>   method   from 
#>   plot.mst spdep
#> $observed
#> [1] 0.4114597
#> $expected
#> [1] -0.01190476
#> $sd
#> [1] 0.06999644
#> $p.value
#> [1] 1.463168e-09

Those results are identical. There are differences in calculating variance but the I is the same.


The p value is telling you if the uncertainty around the estimate includes zero. If you've got a very small estimate and a very very small uncertainty then you'll get a small p-value. It says there is clustering, but not much of it.

You usually need a lot of data to get very very small uncertainties when estimates are very small though. Normally the uncertainty swamps out the estimate if the estimate is very small so that zero can't be excluded.

The different packages use slightly different methods for defining and computing the index and these are detailed in the documentation (and ultimately, the code).

  • Thanks, @Spacedman. My datasets have very variable numbers of data, from eight to 465... I understand that different packages can yield different results but in the case of the example, while the expectations are the same, the observed values are very different. These differences are huge when running my dataset..
    – perep1972
    Sep 22 at 12:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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