Skip to Main Content

How can I calculate Moran’s I in R? | R FAQ

Moran’s I is a measure of spatial autocorrelation–how related the values of a
variable are based on the locations where they were measured.  Using
functions in the ape library, we can calculate Moran’s I in R.  To download
and load this library, enter install.packages(“ape”) and then
library(ape)
.

Let’s look at an example. Our dataset,
ozone, contains ozone measurements from
thirty-two locations in the Los Angeles area aggregated over one month. The dataset
includes the station number (Station), the latitude and longitude of the station
(Lat and
Lon), and the average
of the highest eight hour daily averages (Av8top). This data, and other spatial datasets, can be
downloaded from the University of Illinois’ Spatial Analysis Lab. We can look
at a summary of our location variables to see the range of locations under
consideration.

ozone <- read.table("https://stats.idre.ucla.edu/stat/r/faq/ozone.csv", sep=",", header=T)
head(ozone, n=10)

   Station   Av8top      Lat       Lon
1       60 7.225806 34.13583 -117.9236
2       69 5.899194 34.17611 -118.3153
3       72 4.052885 33.82361 -118.1875
4       74 7.181452 34.19944 -118.5347
5       75 6.076613 34.06694 -117.7514
6       84 3.157258 33.92917 -118.2097
7       85 5.201613 34.01500 -118.0597
8       87 4.717742 34.06722 -118.2264
9       88 6.532258 34.08333 -118.1069
10      89 7.540323 34.38750 -118.5347

To calculate Moran’s I, we will need
to generate a matrix of inverse distance weights.  In the matrix, entries
for pairs of points that are close together are higher than for pairs of points
that are far apart.  For simplicity, we will treat the latitude and
longitude as values on a plane rather than on a sphere–our locations are close
together and far from the poles. When using latitude and longitude coordinates
from more distant locations, it’s wise to calculate distances based on spherical
coordinates (the geosphere package can be used).

We can first generate a distance matrix, then take
inverse of the matrix values and replace the diagonal entries with zero:

ozone.dists <- as.matrix(dist(cbind(ozone$Lon, ozone$Lat)))

ozone.dists.inv <- 1/ozone.dists
diag(ozone.dists.inv) <- 0
 
ozone.dists.inv[1:5, 1:5]
 
         1        2        3        4        5
1 0.000000 2.539795 2.446165 1.627570 5.391160
2 2.539795 0.000000 2.667061 4.531428 1.741071
3 2.446165 2.667061 0.000000 1.954357 2.002389
4 1.627570 4.531428 1.954357 0.000000 1.258716
5 5.391160 1.741071 2.002389 1.258716 0.000000

We have created a matrix where each off-diagonal entry [i,
j] in the
matrix is equal to 1/(distance between point
i
and point
j).
Note that this is just one of several ways in which we can calculate an inverse
distance matrix.  This is the formulation used by Stata.  In SAS,
inverse distance matrices have entries equal to 1/(1+ distance between point
i
and point
j) and there are numerous scaling options available.

We can now calculate Moran’s I using the command Moran.I.

Moran.I(ozone$Av8top, ozone.dists.inv)

$observed
[1] 0.2265501

$expected
[1] -0.03225806

$sd
[1] 0.03431138

$p.value
[1] 4.596323e-14

Based on these results, we can reject the null hypothesis that there is zero spatial autocorrelation present in the variable
Av8top at alpha = .05.

Alternatively, one could use a binary distance matrix.  This might be
used in looking at networking problems where two points are either connected or
not connected.  In this example, we can choose a distance d such
that pairs of points with distance less than d are considered connected
and pairs with distance greater than d are not and use this weights
matrix to calculate Moran’s I.  The code below does this for d =
.75.

ozone.dists.bin <- (ozone.dists > 0 & ozone.dists <= .75)
Moran.I(ozone$Av8top, ozone.dists.bin)

$observed
[1] 0.2045738

$expected
[1] -0.03225806

$sd
[1] 0.04771132

$p.value
[1] 6.910876e-07

In this example, we see similar results to those generated using the inverse Euclidean distance matrix
and we can reject the null hypothesis that there is zero spatial autocorrelation present in the variable
Av8top at alpha = .05.