# Getting different Moran indexes?

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())
library(raster)
# Let's create some data
set.seed(1234)
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)
`````` ``````# 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`

``````library(ape)
Moran.I(dta\$X, dta.dists.inv)
``````

Which yields:

``````\$observed
 0.02912351

\$expected
 -0.02040816

\$sd
 0.03655231

\$p.value
 0.1753889
``````

Then, I continued using the package `spdep`

``````library(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:

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

Obtaining:

`````` -0.02040816
``````

That is: three packages and two different values.

Which one is the correct one?

• 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. Sep 22, 2022 at 14:14
• 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. Sep 22, 2022 at 14:46

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
library(spdep)
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"))

listw\$weights[]
#>  0.02040816 0.02040816 0.02040816 0.02040816 0.02040816 0.02040816 0.02040816
 0.02040816 0.02040816 0.02040816

dta.dists.inv[1,]
#>  0.000000e+00 6.310735e-07 4.059586e-07 5.128319e-07 4.767474e-07 5.054279e-07
 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.

``````library(spdep)

# 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
#>  0.4114597
#>
#> \$expected
#>  -0.01190476
#>
#> \$sd
#>  0.06999644
#>
#> \$p.value
#>  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.. Sep 22, 2022 at 12:10