3

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.

1
  • A reproducible example would still be useful, but we can use the one in help(lm.morantest).
    – Spacedman
    Jun 9, 2020 at 12:01

1 Answer 1

5

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.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.