Skip to content
Prev 25832 / 29559 Next
Original Raw

Differences between moran.test and lm.morantest

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

            

      
      
      
    

        
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/s5bj8r4nwn4J

for a simplified flowchart for choosing models.

Hope this clarifies,

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


  


  
  

    

  

    
      
    
    
Thread (3 messages)

  

  

    
            
        Javier García
        
        Differences between moran.test and lm.morantest
      
      Jul 28
    
    
      
            
        Roger Bivand
        
        Differences between moran.test and lm.morantest
      
      Jul 30
    
    
      
            
        Javier García
        
        Differences between moran.test and lm.morantest
      
      Jul 31