Applying Moran's I across columns

Question

So I find myself needing to apply global Moran's I across ~20 variables (as in ~20 instances of univariate autocorrelation for each variable, not attempted multivariate spatial autocorrelation). I'm using R and the sf + spdep packages.

Where data_lisa is a sf object structured:

id | var_a | var_b | ... | var_n | geometry

...and lw the spatial weights list created with:

lw <- nb2listw(neighbours = poly2nb(data_lisa, 
                                    queen = TRUE), 
               style = "W",
               zero.policy = TRUE)

Using spdep, I can apply global Moran's I for a single variable as:

moran.mc(data_lisa$var_a,
         listw = lw, 
         nsim = 999, 
         zero.policy = TRUE)

...and receive all the expected results.

So I'm looking for help in how to programmatically apply this function across all of my variables.

The result of moran.mc is a list object which I suspect is where I'm encountering the greatest issues, as I don't have much experience interacting with lists.

Ideally the output would look something like this.

variable	moran_stat	pval
var_a	0.064	0.042
var_b	0.322	0.001
var_c	0.183	0.001

How can I do this?

Answer 1

Make a vector of your variable names either by subsetting the column names or programmatically, eg:

vars = names(data_lisa)[2:5]
vars = paste0("var_",letters[1:4])

For an example, I'm using the COL.OLD data you get from ?moran.mc and this set of variables:

vars = c("AREA_PL","PERIMETER","CRIME","POLYID")

then repeat using lapply over the variable names, and give the returned list the names of the vars:

> mcs = lapply(vars, function(v){moran.mc(COL.OLD[[v]], listw=colw, nsim=10)})
> names(mcs) = vars

Then extract the statistic and p-value, making a data frame with the correct names:

> stats = data.frame(
     lapply(names(mcs), 
         function(nm){
            m=mcs[[nm]]
            setNames(
               list(
                  stat=m$statistic,
                  p=m$p.value
                  ),
                  paste0(nm,c("_s","_p"))
             )
          }
        )
   )

Which produces:

> stats
          AREA_PL_s  AREA_PL_p PERIMETER_s PERIMETER_p   CRIME_s    CRIME_p
statistic 0.1404648 0.09090909   0.1741305  0.09090909 0.5109513 0.09090909
           POLYID_s   POLYID_p
statistic 0.8701314 0.09090909

I don't get how you want the output to be a table, since you only get one statistic per column, and not by id as implied in your sample table output. This is a global statistic.

Note this only uses base R packages (plus spdep) so should work in any R installation.

Answer 2

As a possible alternative to the approach suggested by Spacedman I propose one based on {dplyr} and iterating via a for cycle over a vector of names of your variables of interest.

As I don't have access to your data_lisa object I am using the well known & much loved North Carolina shapefile that ships with the {sf} package.

What this piece of code does is that it calculates the Moran's statistic for each variable in the variables vector, and appends the value to the results tibble. I find it easier to work with denormalized (long) data, but if you prefer wide ones there is {tidyr} to the rescue.

library(sf)
library(dplyr)
library(spdep)

# a shapefile
shape <- st_read(system.file("shape/nc.shp", package="sf")) 

# list weight object
lw <- nb2listw(neighbours = poly2nb(shape, 
                                    queen = TRUE), 
               style = "W",
               zero.policy = TRUE)

# variable names to be investigated
variables <- c("BIR74", "SID74", "NWBIR74")

# initate empty resultset
result <- tibble(NULL)

# now let's iterate!! :)
for (i in variables) {
  
# calculate the moran's object (as list)
res <- moran.mc(pull(shape, !!i),
         listw = lw, 
         nsim = 999, 
         zero.policy = TRUE)

# use the res object to create a new row in results dataset
result <- result %>% 
  bind_rows(tibble(variable = i,
                   statistic = res$statistic,
                   pvalue = res$p.value))

}

# check result
result
# A tibble: 3 × 3
  variable statistic pvalue
  <chr>        <dbl>  <dbl>
1 BIR74        0.139  0.015
2 SID74        0.148  0.018
3 NWBIR74      0.180  0.004