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

Sat Mar 5 20:39:06 CET 2016

On Fri, 4 Mar 2016, Roger Bivand wrote:

> 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:

This is definitely the case - setting

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

sigma = epe;

gives exactly the same results as test="LMlag" in spdep::lm.LMtests()

So the reported differences are due to 1) the user not using the same 
data in R and Matlab, and 2) the Matlab code being undecided about 
dividing the sum of squared errors by n or n-k, and consequently making 
choices that might seem sensible but aren't fully supported in the 
underlying article.

Roger

> 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