[R-sig-Geo] [Local.moran.exact] Error in integrate : evaluation of function gave a result of wrong type

[Roger Bivand]

On Mon, 17 Jul 2017, Lisa Menez wrote:

> Good morning to all 
>
> I am currently trying to use local.moran.exact function from spdep package.
> More precisely, I am willing to use the alternative version (local.moran.exact.alt) 
> but I obtain an error message I can’t figure out :
>
> "Error in integrate(integrand, lower = 0, upper = upper) : evaluation of function gave a result of wrong type"
> 
>
> #======================================
> #My issue is somewhere is those lines : 
> #======
> setwd(".")
> library(spdep)
> EU_NUTS <- read.table("./mydata.txt")
>
> ###distance based neighbors - - - 
> Coord <- cbind(EU_NUTS$"x", EU_NUTS$"y")
> EU_NUTS_kn0 <- knn2nb(knearneigh(Coord, k = 1), row.names =EU_NUTS$NUTS_ID)
> dist <- unlist(nbdists(EU_NUTS_kn0, Coord))
> summary(dist)
> max_k0 <- max(dist)
> EU_NUTS_kd1 <- dnearneigh(Coord, d1 = 0, d2 = 1 * max_k0, row.names =EU_NUTS$NUTS_ID)
> Kd1list<-nb2listw(EU_NUTS_kd1, glist=NULL, style="W", zero.policy=FALSE)
>
>  ###OLS model
>  OLS<-lm(EU_NUTS$"GrowthR" ~ log(EU_NUTS$"X2000"))
>
>  ###SAR model
>  ERROR<-errorsarlm(OLS,, Kd1list , na.omit, etype="emixed", method="eigen", quiet=NULL,tol.solve=1.0e-10, zero.policy= TRUE)
>  lm.error <- lm(ERROR$tary ~ ERROR$tarX - 1)
>
>  Omega <- invIrW(Kd1list, rho=ERROR$lambda)
>  Omega <- tcrossprod(Omega)
>
>
> localmoran.exact.alt(lm.error, nb=EU_NUTS_kd1, Omega=Omega, style = "S",  zero.policy = TRUE,  alternative = "greater", spChk = NULL,  resfun = residuals,
> 				 save.Vi = FALSE, useTP=TRUE,  truncErr=1e-6,  zeroTreat=0.1)
>
>
> #======
> #The error message is get is  : 
> #======
>
> Error in integrate(integrand, lower = 0, upper = upper) : evaluation of function gave a result of wrong type
>
> #=======================
>
>
> Does someone have any idea of what I am doing wrong ?
>

Thanks for providing a reprex.

First, the points you are using are in geographical coordinates, which you 
ignore when calculating distances, so the distance criterion and the 
first nearest neighbours are wrong.

Second, you include all NUTS regions, including Reunion and other very 
overseas territories, with the consequence that most regions are defined 
as each other's neighbours (there are ~ 50% non-zero weights). This is far 
from sparse! Please motivate your choice of regions - it is very hard to 
argue that Reunion influences its nearest neighbours (Crete, Cyprus??). 
With so many neighbours on average, the spatial weights are all small.

Third, you do not motivate using the error Durbin model. If you do not use 
it, localmoran.exact.alt() does complete. However, few of the exact 
p-values for alternative="greater" are "significant" (t0 with Omega=Omega, 
t1 with Omega=diag(276) - the identity matrix):

> table(findInterval(as.data.frame(t0)[,5], c(0, 0.001, 0.01, 0.1, 0.5, 
0.9, 0.99, 0.999, 1)))

   1   2   3   4   5   6   7
   1   6  34 108 104  21   2
> table(findInterval(as.data.frame(t1)[,5], c(0, 0.001, 0.01, 0.1, 0.5, 
0.9, 0.99, 0.999, 1)))

  1  2  3  4  5  6  7
  7 13 48 94 90 22  2

So using the alternative model taking account of global autocorrelation 
does have an effect (7 not 20 regions with unadjusted p < 0.01), but for 
an odd selection of units and weights).

Note that p-values must be adjusted for multiple comparisons, for example:

> table(findInterval(p.adjust(as.data.frame(t0)[,5], "fdr"), c(0, 0.001, 
0.01, 0.1, 0.5, 0.9, 0.99, 0.999, 1)))

   3   4   5   6   7
   1  13 133 125   4
>
> table(findInterval(p.adjust(as.data.frame(t1)[,5], "fdr"), c(0, 0.001, 
0.01, 0.1, 0.5, 0.9, 0.99, 0.999, 1)))

  1  2  3  4  5  6  7
  1  2 11 91 71 96  4


So unless we know that you actually wanted to do what you did, 1-3 above 
suggest that your choices need to be motivated.

Hope this clarifies,

Roger

>
>
> The full dataset is available here : 
> http://www.i3s.unice.fr/~menez/RCode/mydata.txt <http://www.i3s.unice.fr/~menez/RCode/mydata.txt>
>
> The Rcode is here : http://www.i3s.unice.fr/~menez/RCode/moran.R <http://www.i3s.unice.fr/~menez/RCode/moran.R>
>
>
>
> Thank you for your attention and your help.
> Best Regards
>
> Lisa
>
>
>
>
>
> 	[[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; 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