Adding localC()

[rsbivand]

@JosiahParry has helpfully contributed a prototype implementation of localC() in #66 #67 .

A) One problem being examined is how to deal with no-neighbour observations, this may be perhaps handed by using C code already in place for handing the same problem for global Geary's C, but doubly scaling the input variable.

B) A further open problem is that the prototype permutation uses all of the observations (total permutation), but we know that for local Moran's I and local G permutatons, the current observation i is held out of the srt (conditional permutation). This relates to https://doi.org/10.1111/gean.12304 (https://doi.org/10.31219/osf.io/ugkhp), already contributed to spdep as #58 (and also very welcome @jeffcsauer.)

[rsbivand]

Addressing A) above: 732a6e2 adds zero.policy= and branching to legacy code looping in C when no-neighbour observations are present.

[rsbivand]

Addressing B) rgeoda and PySAL esda use conditional permutation. Have commited and pushed branching code, default conditional=FALSE 750510a. Analysis in help page is that rgeoda and conditional=TRUE agree well, but that rgeoda does not fold back two-sided >= 0.5 correctly. Using reticulate, PySAL and esda$Geary_Local() turn out even odder, with output re-ordered in addition to use of ddof=1 in np$std() creating intial differences. That taken into account, the local G values agree (re-calibrating to match esda in scale(x)), but the simulations are strangely ordered even given the order of the local Gs:

orig <- spData::africa.rook.nb
listw <- nb2listw(orig)
x <- spData::afcon$totcon

(A <- localC(x, listw))
listw1 <- nb2listw(droplinks(sym.attr.nb(orig), 3, sym=TRUE), zero.policy=TRUE)
(A1 <- localC(x, listw1, zero.policy=FALSE))
(A2 <- localC(x, listw1, zero.policy=TRUE))
if (require(rgeoda, quietly=TRUE)) {
  W <- create_weights(as.numeric(length(x)))
  for (i in 1:length(listw$neighbours)) {
    set_neighbors_with_weights(W, i, listw$neighbours[[i]], listw$weights[[i]])
    update_weights(W)
  }
  set.seed(1)
  B <- local_geary(W, data.frame(x))
  all.equal(A, lisa_values(B))
}
set.seed(1)
C <- localC_perm(x, listw, nsim = 499, conditional=TRUE, alternative="two.sided")
cor(ifelse(lisa_pvalues(B) < 0.5, lisa_pvalues(B), 1-lisa_pvalues(B)), attr(C, "pseudo-p")[,6])
# pseudo-p values probably wrongly folded
\dontrun{
library(reticulate)
use_python("/usr/bin/python", required = TRUE)
gp <- import("geopandas")
ps <- import("libpysal")
tf <- tempfile(fileext=".gal")
write.nb.gal(orig, tf)
fl <- ps$io$open(tf)
w$transform <- "R"
esda <- import("esda")
lM <- esda$Moran_Local(x, w)
all.equal(unname(localmoran(x, listw, mlvar=FALSE)[,1]), c(lM$Is))
# confirm x and w the same
lG <- esda$Geary_Local(connectivity=w)$fit(scale(x))
# np$std missing ddof=1 
n <- length(x)
D01 <- spdep:::geary.intern((x - mean(x)) / sqrt(var(x)*(n-1)/n), listw, n=n)
# lG components probably wrongly ordered
o <- match(round(D0, 6), round(lG$localG, 6))
all.equal(c(lG$localG)[o], D0)
# simulation order not retained
lG$p_sim[o]
attr(C, "pseudo-p")[,6]
}

@JosiahParry comments please on current main state. @ljwolf, @sjsrey, what am I not getting about the esda implementation?

[ljwolf]

Hey, thanks for asking!

At least for the local Geary ordering, I believe we have an open ticket, I just haven't had the time to get the patch wrangled and tested!