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

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
LM
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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

LM
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 ?



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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: [Local.moran.exact] Error in integrate : evaluation of function gave a result of wrong type

Roger Bivand
Administrator
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
> [hidden email]
> 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: [hidden email]
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
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Loading...