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

```