[R-sig-Geo] [localmoran.exact]

Wed Aug 15 12:50:27 CEST 2018

> Le 14 août 2018 à 20:09, Roger Bivand <Roger.Bivand using nhh.no> a écrit :
> 
> If you post, you then reply in-thread on the list. Otherwise, nobody else knows what the resolution was. : 
Sorry. My bad. 
> 
> Your explanation below contains lots of inaccuracies - all Spatial*DataFrames have standard data.frame access methods as I showed.
No doubt about it ! I will try do better now that I know about it.  

Thank you Sir.
Lisa 

> 
> Roger Bivand
> 
> On Tue, 14 Aug 2018, Lisa Menez wrote:
> 
>> Sir Bivand,
>> 
>> The complete example was here :
>> 
>> My database is available here : http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt> and my code : http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt>
>> 
>> The localmoran.sad() actually works. Is it because « If we assume global spatial independence, we can give for any projection matrix M the spectrum of eigenvalues as well as the saddle point approximation of local Moran’s Ii in analytical terms. This increases the computational efficiency as it avoids reverting to numerical calculations of the spectrum of eigenvalues and finding iteratively the root of the saddle point equation » ?
>> 
>> To the best of my knowledge, I am constrained by the settings of my data (extracted from Eurostat) on calling EU_NUTS using data rather than EU_NUTS. Indeed EU_NUTS are spatial polygons and EU_NUTS using data their relative attributes. I failed to change this name.
>> 
>> for each i, EU_NUTS using data[,year[i]] is a numerical vector. year is a character vector containing the sample years, allowing to go through the corresponding EU_NUTS using data' columns.
>> 
>> On your advice, I tried to go through debug(localmoran.exact). As I could not find any good reason why (2*t2<t1^2), I turned to your last remark : Changing NUTS_L1 into a factor before using it into the loop. It actually was the issue. I miscoded the factor, keeping character chains and not transforming them into numbers.
>> 
>> 
>> Thank you very much for your help and recommandations.
>> Lisa
>> 
>> 
>>> Le 13 août 2018 à 16:12, Roger Bivand <Roger.Bivand using nhh.no> a écrit :
>>> 
>>> On Mon, 13 Aug 2018, Lisa Menez wrote:
>>> 
>>>> 
>>>> Dear all,
>>>> 
>>>> I have been trying to obtain local Moran’s I with the exact approach (localmoran.exact) but I have some issues and would need your help.
>>> 
>>> A complete reproducible example is required (script and data, data may be downloaded, do not attach). Does the same thing happen with localmoran.sad()?
>>> 
>>>> 
>>>> I am trying to bring to light how spatial patterns of wealth have evolved from 2000 to 2016 in the European Union.
>>>> 
>>>> Looking at the Regional Gross Domestic Product (GDP) of European regions, my lm model looks like this :
>>>> 
>>>> lm(log(EU_NUTS using data[,year[i]])~1+factor(EU_NUTS using data$NUTS_L1))
>>>> 
>>>> where log(EU_NUTS using data[,year[i]]) are log of regional GDPs
>>>> 
>>> 
>>> Never, ever use @ to access data. EU_NUTS using data$NUTS_L1 should be EU_NUTS$NUTS_L1 and should be (made a) factor first. It isn't clear what log(EU_NUTS using data[,year[i]]) will end up being - is "year" a character vector? If so, EU_NUTS[[year[i]]] is the correct syntax.
>>> 
>>>> and factor(EU_NUTS using data$NUTS_L1) is the vector of country dummies.
>>>> 
>>>> Sadly, I obtain an error message that says :
>>>> 
>>>> Error in integrate(integrand, lower = 0, upper = upper) :
>>>> a limit is missing
>>>> Moreover : Warning messages:
>>>> 1: In sqrt(2 * t2 - t1^2) : production of NaN
>>>> 2: In sqrt(2 * t2 - t1^2) : production of NaN
>>>> 
>>>> According to the GitHub source code ,
>>> 
>>> No, as you can see from https://cran.r-project.org/package=spdep, the development site of spdep is https://github.com/r-spatial/spdep/. Never use non-authorised copies. To see the source, simply type the function name. Use debug(localmoran.exact) to step through the function while it runs, to see the values of the objects you are in doubt about. Do not do this in RStudio, its debugging interface tries to be far too clever.
>>> 
>>> The following is illegible in plain text; post in plain text only.
>>> 
>>> Roger
>>> 
>>>> 
>>>>      Vi <- listw2star <https://rdrr.io/cran/spdep/man/localmoran.sad.html>(B, select[i], style=style, n, D <https://rdrr.io/r/stats/deriv.html>, a,
>>>> 	    zero.policy=zero.policy)
>>>>      Viu <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, u, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
>>>> 	Ii <- c <https://rdrr.io/r/base/c.html>(crossprod <https://rdrr.io/r/base/crossprod.html>(u, Viu) / utu)
>>>>      ViX <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, X, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
>>>>      MViM <- t <https://rdrr.io/r/base/t.html>(X) %*% ViX %*% XtXinv
>>>>      t1 <- -sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM))
>>>>      sumsq.Vi <- function <https://rdrr.io/r/base/function.html>(x) {
>>>>          if <https://rdrr.io/r/base/Control.html> (is.null <https://rdrr.io/r/base/NULL.html>(x)) NA <https://rdrr.io/r/base/NA.html>
>>>> 	    else <https://rdrr.io/r/base/Control.html> sum <https://rdrr.io/r/base/sum.html>(x^2)
>>>> 	}
>>>> 	trVi2 <- sum <https://rdrr.io/r/base/sum.html>(sapply <https://rdrr.io/r/base/lapply.html>(Vi$weights <https://rdrr.io/r/stats/weights.html>, sumsq.Vi), na.rm=TRUE <https://rdrr.io/r/base/logical.html>)
>>>> 	t2a <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(t <https://rdrr.io/r/base/t.html>(ViX) %*% ViX %*% XtXinv))
>>>> 	t2b <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM %*% MViM))
>>>> 	t2 <- trVi2 - 2*t2a + t2b
>>>> 	e1 <- 0.5 * (t1 + sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
>>>> 	en <- 0.5 * (t1 - sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
>>>> 
>>>> My guess is that my setting causes (2*t2<t1^2) : Is it ? Why would it ?
>>>> 
>>>> My database is available here : http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt>
>>>> and my code : http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt>
>>>> 
>>>> I thank you very much for your attention and your help.
>>>> 
>>>> Best Regards
>>>> Lisa
>>>> 
>>>> PhD Student
>>>> GREDEG, Université Côté d’Azur
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 	[[alternative HTML version deleted]]
>>>> 
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo using 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 using nhh.no
>>> http://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>> 
>> 
> 
> -- 
> 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 using nhh.no
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en