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:

         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?

2 Answers 2


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:


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(

Which produces:

> stats
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.

  • Yes my bad on the output formatting. Applying LISA would be the next step and I mixed them up. I've edited my question to remove that part as it's not relevant in this instance.
    – efor027
    May 24 at 22:50
  • I accepted Jindra's answer as it uses a tidyverse approach, but your solution also worked too. Thanks!
    – efor027
    May 24 at 22:52
  • 1
    @efor027, Spacedman's solution is more efficient and elegant than the tidy-based one. Saying that tidy is your primary criteria for a successful solution is a non-sequitur and illustrates a bit of a dogmatic attitude towards what constitutes a "correct" solution that is becoming prevalent in the tidyverse. May 24 at 23:03
  • Thanks for comments on the efficiency of the two approaches Jeffrey. I'll keep it in mind, especially as I may need to progress onto iterating across ~300 columns at a later date.
    – efor027
    May 24 at 23:53
  • 1
    Same re X vs Y. I was just nonplussed that my answer which returned the result in the precise form required by the Q (modulo the ID value in each row which didn't make any sense) wasn't accepted over an answer that returned a different output. FWIW I think the output format from @JindraLacko answer is a better way, but it wasn't what was asked for. Maybe the Q needs an edit to clarify the output format and then I'll tweak my code if needed? The latest edit doesn't specify an output format at all...
    – Spacedman
    May 25 at 9:06

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.


# 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
# 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

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.