[R-sig-Geo] different Moran's I results using R and matlab

[Roger Bivand]

On Thu, 3 Mar 2016, Roger Bivand wrote:

> On Thu, 3 Mar 2016, Qiuhua Ma wrote:
>
>>  Hi,
>>
>>  I run the exactly same regression using R and matlab and get the same
>>  regression results.
>
> Well, you need to give the exact reference to the matlab code you have used - 
> are they the functions in Spatial Econometrics Toolbox under spatial/stats? I 
> do not see the robust variants there.
>
> Most likely the data or the weights are different. Make your data set 
> available on a link, and I'll take a look.

Having not heard back, I ran moran() and lmerror() in Matlab using the 
Spatial Econometrics Toolbox on the data presented in stats/moran_d.m:

load anselin.dat;

y = anselin(:,1);
n = length(y);

x = [ones(n,1) anselin(:,2:3)];

xc = anselin(:,4);
yc = anselin(:,5);
[j W j] = xy2cont(xc,yc);

result = moran(y,x,W);

> tests$result
, , 1

        [,1]
meth   "moran"
nobs   49
nvar   3
morani 0.2861962
istat  4.019423
imean  -0.03180435
ivar   0.006259337
prob   5.834084e-05

(R result after moving the input data from Matlab to R with 
R.matlab::readMat())

> lm.morantest(lm_obj, lw, alternative="two.sided")

 	Global Moran I for regression residuals

data:
model: lm(formula = y ~ x - 1)
weights: lw

Moran I statistic standard deviate = 4.0194, p-value = 5.834e-05
alternative hypothesis: two.sided
sample estimates:
Observed Moran I      Expectation         Variance
      0.286196209     -0.031804345      0.006259337

[SAME RESULT]

and similarly with lmerror()

> tests$result2
, , 1

      [,1]
meth "lmerror"
lm   10.7656
prob 0.001034042
chi1 17.611
nobs 49
nvar 3


> lm.LMtests(lm_obj, lw)

 	Lagrange multiplier diagnostics for spatial dependence

data:
model: lm(formula = y ~ x - 1)
weights: lw

LMErr = 10.766, df = 1, p-value = 0.001034

[SAME RESULT]

The result for lmlag() differs, probably because the Matlab code uses two 
different values for sigma and epe:

sigma  = (e'*e)/(n-k);

epe = (e'*e)/n;
lm1 = (e'*W*y)/epe;

t1 = trace((W+W')*W);
D1=W*x*b;
M=eye(n)-x*inv((x'*x))*x';
D=(D1'*M*D1)*(1/sigma)+t1;

lmlag = (lm1*lm1)*(1/D);

so does not match the R code which follows Eq. 13 in Anselin et al. (1996, 
p. 84) in dividing the sum of squared errors by n, not n-k.

I have not found lmlag_robust() or lmerror_robust() anywhere.

Please provide the code of these functions if you need further 
clarification - the R code for the Moran's I test on OLS residuals is OK, 
as is the LM error test. The LM lag test chooses a version using ML sigma 
in harmony with the source article.

Roger

>
> Roger
>
>>
>>  However I got different results for Moran't I test and LM test results.
>>
>>  *R command:*
>>  nb4 <- knn2nb(knearneigh(sale.sp, k = 4))
>>  knn4listw <- nb2listw(nb4, style="W")
>>
>>  lm.morantest(near7_comrisk_semilog.out, knn4listw)
>>  lm.LMtests(near7_comrisk_semilog.out, knn4listw, test=c("LMerr", "LMlag",
>>  "RLMerr", "RLMlag"))
>>
>>  *R results:*
>>  Moran I statistic standard deviate = 1.4135, p-value = 0.07875
>>  Observed Moran's I        Expectation           Variance
>>       8.780168e-03      -2.762542e-04       4.104967e-05
>>
>>  LMerr = 1.8752, df = 1, p-value = 0.1709
>>  LMlag = 0.14335, df = 1, p-value = 0.705
>>  RLMerr = 2.1193, df = 1, p-value = 0.1455
>>  RLMlag = 0.38741, df = 1, p-value = 0.5337
>>
>>  *Matlab command:*
>>  W_sale_4 = make_nnw(xc,yc,4);
>>
>>  moran4= moran(y,x, W_sale_4)
>>  error_4 = lmerror(y,x,W_sale_4)
>>  lag_4 = lmlag(y,x,W_sale_4)
>>  error_4r = lmerror_robust(y,x,W_sale_4)
>>  lag_4r = lmlag_robust(y,x,W_sale_4)
>>
>>  *Matlab results:*
>>  Moran I-test for spatial correlation in residuals
>>  Moran I                    0.13979528
>>  Moran I-statistic         22.05018073
>>  Marginal Probability       0.00000000
>>  mean                      -0.00126742
>>  standard deviation         0.00639735
>>
>>  error_4 =
>>     meth: 'lmerror'
>>       lm: 475.1839
>>     prob: 0
>>     chi1: 17.6110
>> 
>>
>>  lag_4 =
>>     meth: 'lmlag'
>>       lm: 465.4518
>>     prob: 0
>>     chi1: 17.6110
>>
>>  error_4r =
>>     meth: 'lmerror_robust'
>>       lm: 76.5845
>>     prob: 0
>>     chi1: 6.6400
>>
>>  lag_4r =
>>
>>     meth: 'lmlag_robust'
>>       lm: 66.1547
>>     prob: 4.4409e-16
>>     chi1: 6.6400
>>
>>  Did I do anything wrong? Any thought on this problem?
>>
>>  thanks,
>>
>>  qiuhua
>>
>>   [[alternative HTML version deleted]]
>>
>>  _______________________________________________
>>  R-sig-Geo mailing list
>>  R-sig-Geo at r-project.org
>>  https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412