R – Applying Moran’s I Across Columns for Spatial Autocorrelation

autocorrelationrsfspdep

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?

Best Answer

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.

Related Question