Differences between moran.test and lm.morantest

Mon, Jul 31, 2017 1:37 PM

Thanks a lot for such a detailed response, Roger.

Best
Javi

-----Mensaje original-----
De: Roger Bivand [mailto:Roger.Bivand at nhh.no] 
Enviado el: domingo, 30 de julio de 2017 18:42
Para: Javier Garc?a
CC: r-sig-geo at r-project.org
Asunto: Re: [R-sig-Geo] Differences between moran.test and lm.morantest

On Sat, 29 Jul 2017, Javier Garc?a wrote:

variance differ dramatically, getting opposite conclusions.

moran.test() is for "primary" variables only - read the reference in
?moran.test. The mean model applied to the this variable is the intercept,
that is the mean only:

+ randomisation=FALSE, alternative="two.sided")

 	Moran I test under normality

data:  COL.OLD$CRIME
weights: nb2listw(COL.nb, style = "W")

Moran I statistic standard deviate = 5.6754, p-value = 1.384e-08
alternative hypothesis: two.sided
sample estimates:
Moran I statistic       Expectation          Variance
       0.510951264      -0.020833333       0.008779831

under the Normal assumption.

lm.morantest() can reproduce the same results for the same model:

+ alternative="two.sided")

 	Global Moran I for regression residuals

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

Moran I statistic standard deviate = 5.6754, p-value = 1.384e-08
alternative hypothesis: two.sided
sample estimates:
Observed Moran I      Expectation         Variance
      0.510951264     -0.020833333      0.008779831

Note that the Normal VI for moran.test() is:

         VI <- (wc$nn * wc$S1 - wc$n * wc$S2 + 3 * S02)/(S02 *
             (wc$nn - 1))

and for lm.morantest():

     XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
     X <- model.matrix(terms(model), model.frame(model))
...
     Z <- lag.listw(listw.U, X, zero.policy = zero.policy)
     C1 <- t(X) %*% Z
     trA <- (sum(diag(XtXinv %*% C1)))
     EI <- -((N * trA)/((N - p) * S0))
     C2 <- t(Z) %*% Z
     C3 <- XtXinv %*% C1
     trA2 <- sum(diag(C3 %*% C3))
     trB <- sum(diag(4 * (XtXinv %*% C2)))
     VI <- (((N * N)/((S0 * S0) * (N - p) * (N - p + 2))) * (S1 +
         2 * trA2 - trB - ((2 * (trA^2))/(N - p))))

where you see easily that the issue of dependent residuals (only n-p are 
independent even in the aspatial case, so you see n-p instead of n) is 
handled, as are products of X, WX, and (X'X)^{-1}.

Reading the references is essential and should not be neglected. Reading 
the code may help, but doesn't explain the reasoning behind the code and 
the output. If you need to test model residuals, use the appropriate test. 
Using moran.test() on extracted residuals if covariates are included in 
the model is never justified.

Different outcomes from Moran's I and LM suggest severe mis-specification 
in your model, see for example:

https://groups.google.com/forum/#!msg/openspace-list/k4F4jI9cU1I/s5bj8r4nwn4
J

for a simplified flowchart for choosing models.

Hope this clarifies,

Roger

data

http://www.unibertsitate-hedakuntza.ehu.es/p268-content/es/contenidos/inform

acion/manual_id_corp/es_manual/images/firma_email_upv_euskampus_bilingue.gif

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

Differences between moran.test and lm.morantest

Thread (3 messages)