Moran's I z-value in spdep

Question

I wish to incrementally find the neighbour distance which maximises Moran's I test statistic. I have read that the I values are not comparable across varying spatial weight matrices but that the z-value should be maximised instead. Where can I find/how can I compute the z-value in lm.morantest from spdep package? I don't see a z-value either in the value or in the structure of lm.morantest printed to the console (same applies to moran.test for a single variable):

library(spdep)

lm.morantest(lm(pat_ct~pub_ct, sca_ct), w_dist, alternative="two.sided")
#Global Moran I for regression residuals

#data:
#model: lm(formula = pat_ct ~ pub_ct, data = sca_ct)
#weights: w_dist

#Moran I statistic standard deviate = 12.155, p-value < 2.2e-16
#alternative hypothesis: two.sided
#sample estimates:
#Observed Moran I      Expectation         Variance
#    3.556003e-02    -8.665431e-05     8.600474e-06

str(lm.morantest(lm(pat_ct~pub_ct, sca_ct), w_dist, alternative="two.sided"))
#List of 6
# $ statistic  : num [1, 1] 12.2
#  ..- attr(*, "names")= chr "Moran I statistic standard deviate"
# $ p.value    : num [1, 1] 5.39e-34
# $ estimate   : Named num [1:3] 3.56e-02 -8.67e-05 8.60e-06
#  ..- attr(*, "names")= chr [1:3] "Observed Moran I" "Expectation" "Variance"
# $ method     : chr "Global Moran I for regression residuals"
# $ alternative: chr "two.sided"
# $ data.name  : chr "\nmodel: lm(formula = pat_ct ~ pub_ct, data = sca_ct)\nweights: w_dist\n"
# - attr(*, "class")= chr "htest"

I assume a minimal reproducible example is not required to answer this generic question.

Answer 1

A "Z-score" (I wouldn't call it a Z-value) is the number of standard deviations that a statistic is away from its expected mean.

The ESRI documentation https://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-how-spatial-autocorrelation-moran-s-i-spatial-st.htm shows the computation of Z_I.

For the first example in help(lm.morantest), we see:

> m

    Global Moran I for regression residuals

data:  
model: lm(formula = CRIME ~ HOVAL + INC, data = COL.OLD)
weights: nb2listw(COL.nb, style = "W")

Moran I statistic standard deviate = 2.9539, p-value = 0.003138
alternative hypothesis: two.sided
sample estimates:
Observed Moran I      Expectation         Variance 
     0.235638354     -0.033302866      0.008289408 

And so we can compute Z_I from the ESRI formula and those three last printed values (ignore the "Observed Moran I" label - this is hung over from the name of the first element in the calculation):

> (m$estimate["Observed Moran I"] - m$estimate["Expectation"])/sqrt(m$estimate["Variance"])
Observed Moran I 
        2.953899 

Hmmm that number looks familiar...

    Global Moran I for regression residuals

data:  
model: lm(formula = CRIME ~ HOVAL + INC, data = COL.OLD)
weights: nb2listw(COL.nb, style = "W")

Moran I statistic standard deviate = ***--> 2.9539 <--***, p-value = 0.003138
alternative hypothesis: two.sided

If you look at the code you'll even see:

ZI <- (I - EI)/sqrt(VI)
...
statistic <- ZI
attr(statistic, "names") <- "Moran I statistic standard deviate"

It looks like Roger has preferred to call it the statistic standard deviate, even though its stored in a variable ZI, but I think its the same. There's probably some reason to not call it a Z-score related to the distribution of the statistic.