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))
run <- FALSE
if (require(rgeoda, quietly=TRUE)) run <- TRUE
if (run) {
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))
}
if (run) {
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 https://github.com/GeoDaCenter/rgeoda/issues/28
if (FALSE) {
tmap_ok <- FALSE
if (require(tmap, quietly=TRUE)) tmap_ok <- TRUE
if (run) {
# doi: 10.1111/gean.12164
guerry_path <- system.file("extdata", "Guerry.shp", package = "rgeoda")
g <- st_read(guerry_path)[, 7:12]
cor(st_drop_geometry(g)) #(Tab. 1)
lw <- nb2listw(poly2nb(g))
moran(g$Crm_prs, lw, n=nrow(g), S0=Szero(lw))$I
moran(g$Crm_prp, lw, n=nrow(g), S0=Szero(lw))$I
moran(g$Litercy, lw, n=nrow(g), S0=Szero(lw))$I
moran(g$Donatns, lw, n=nrow(g), S0=Szero(lw))$I
moran(g$Infants, lw, n=nrow(g), S0=Szero(lw))$I
moran(g$Suicids, lw, n=nrow(g), S0=Szero(lw))$I
}
if (run) {
o <- prcomp(st_drop_geometry(g), scale.=TRUE)
cor(st_drop_geometry(g), o$x[,1:2])^2 #(Tab. 2)
}
if (run) {
g$PC1 <- o$x[, "PC1"]
brks <- c(min(g$PC1), natural_breaks(k=6, g["PC1"]), max(g$PC1))
if (tmap_ok) tm_shape(g) + tm_fill("PC1", breaks=brks, midpoint=0) +
tm_borders() # Fig. 1
else pplot(g["PC1"], breaks=brks)
}
if (run) {
g$PC2 <- -1*o$x[, "PC2"] # eigenvalue sign arbitrary
brks <- c(min(g$PC2), natural_breaks(k=6, g["PC2"]), max(g$PC2))
if (tmap_ok) tm_shape(g) + tm_fill("PC2", breaks=brks, midpoint=0) +
tm_borders() # Fig. 2
else plot(g["PC2"], breaks=brks)
}
if (run) {
w <- queen_weights(g)
lm_PC1 <- local_moran(w, g["PC1"], significance_cutoff=0.01,
permutations=99999)
g$lm_PC1 <- factor(lisa_clusters(lm_PC1), levels=0:4,
labels=lisa_labels(lm_PC1)[1:5])
is.na(g$lm_PC1) <- g$lm_PC1 == "Not significant"
g$lm_PC1 <- droplevels(g$lm_PC1)
if (tmap_ok) tm_shape(g) + tm_fill("lm_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 3
else plot(g["lm_PC1"])
}
if (run) {
set.seed(1)
lm_PC1_spdep <- localmoran_perm(g$PC1, lw, nsim=9999)
q <- attr(lm_PC1_spdep, "quadr")$pysal
g$lm_PC1_spdep <- q
is.na(g$lm_PC1_spdep) <- lm_PC1_spdep[,6] > 0.02 # note folded p-values
g$lm_PC1_spdep <- droplevels(g$lm_PC1_spdep)
if (tmap_ok) tm_shape(g) + tm_fill("lm_PC1_spdep", textNA="Insignificant",
colorNA="gray95") + tm_borders() # rep. Fig. 3
else plot(g["lm_PC1_spdep"])
}
if (run) {
lg_PC1 <- local_g(w, g["PC1"], significance_cutoff=0.01,
permutations=99999)
g$lg_PC1 <- factor(lisa_clusters(lg_PC1), levels=0:2,
labels=lisa_labels(lg_PC1)[0:3])
is.na(g$lg_PC1) <- g$lg_PC1 == "Not significant"
g$lg_PC1 <- droplevels(g$lg_PC1)
if (tmap_ok) tm_shape(g) + tm_fill("lg_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 4 (wrong)
else plot(g["lg_PC1"])
g$lg_PC1a <- cut(g$PC1, c(-Inf, mean(g$PC1), Inf), labels=c("Low", "High"))
is.na(g$lg_PC1a) <- lisa_pvalues(lg_PC1) >= 0.01
g$lg_PC1a <- droplevels(g$lg_PC1a)
if (tmap_ok) tm_shape(g) + tm_fill("lg_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 4 (guess)
else plot(g["lg_PC1"])
}
if (run) {
lc_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.01,
permutations=99999)
g$lc_PC1 <- factor(lisa_clusters(lc_PC1), levels=0:4,
labels=lisa_labels(lc_PC1)[1:5])
is.na(g$lc_PC1) <- g$lc_PC1 == "Not significant"
g$lc_PC1 <- droplevels(g$lc_PC1)
if (tmap_ok) tm_shape(g) + tm_fill("lc_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 5
else plot(g["lc_PC1"])
}
if (run) {
set.seed(1)
system.time(lc_PC1_spdep <- localC_perm(g$PC1, lw, nsim=9999,
alternative="two.sided"))
}
if (run) {
if (require(parallel, quietly=TRUE)) {
ncpus <- detectCores()-1L
cores <- get.coresOption()
set.coresOption(ncpus)
system.time(lmc_PC1_spdep1 <- localC_perm(g$PC1, lw, nsim=9999,
alternative="two.sided", iseed=1))
set.coresOption(cores)
}
}
if (run) {
g$lc_PC1_spdep <- attr(lc_PC1_spdep, "cluster")
is.na(g$lc_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.01
g$lc_PC1_spdep <- droplevels(g$lc_PC1_spdep)
if (tmap_ok) tm_shape(g) + tm_fill("lc_PC1_spdep", textNA="Insignificant",
colorNA="gray95") + tm_borders() # rep. Fig. 5
else plot(g["lc_PC1_spdep"])
}
if (run) {
g$both_PC1 <- interaction(g$lc_PC1, g$lm_PC1)
g$both_PC1 <- droplevels(g$both_PC1)
if (tmap_ok) tm_shape(g) + tm_fill("both_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 6
else plot(g["both_PC1"])
}
if (run) {
lc005_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.005,
permutations=99999)
g$lc005_PC1 <- factor(lisa_clusters(lc005_PC1), levels=0:4,
labels=lisa_labels(lc005_PC1)[1:5])
is.na(g$lc005_PC1) <- g$lc005_PC1 == "Not significant"
g$lc005_PC1 <- droplevels(g$lc005_PC1)
if (tmap_ok) tm_shape(g) + tm_fill("lc005_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 7
else plot(g["lc005_PC1"])
}
if (run) {
g$lc005_PC1_spdep <- attr(lc_PC1_spdep, "cluster")
is.na(g$lc005_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.005
g$lc005_PC1_spdep <- droplevels(g$lc005_PC1_spdep)
if (tmap_ok) tm_shape(g) + tm_fill("lc005_PC1_spdep", textNA="Insignificant",
colorNA="gray95") + tm_borders() # rep. Fig. 7
else plot(g["lc005_PC1_spdep"])
}
if (run) {
lc001_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.001,
permutations=99999)
g$lc001_PC1 <- factor(lisa_clusters(lc001_PC1), levels=0:4,
labels=lisa_labels(lc001_PC1)[1:5])
is.na(g$lc001_PC1) <- g$lc001_PC1 == "Not significant"
g$lc001_PC1 <- droplevels(g$lc001_PC1)
if (tmap_ok) tm_shape(g) + tm_fill("lc001_PC1", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 8
else plot(g["lc001_PC1"])
if (run) {
g$lc001_PC1_spdep <- attr(lc_PC1_spdep, "cluster")
is.na(g$lc001_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.001
g$lc001_PC1_spdep <- droplevels(g$lc001_PC1_spdep)
if (tmap_ok) tm_shape(g) + tm_fill("lc001_PC1_spdep", textNA="Insignificant",
colorNA="gray95") + tm_borders() # rep. Fig. 8
else plot(g["lc001_PC1_spdep"])
}
}
if (run) {
lc_PC2 <- local_geary(w, g["PC2"], significance_cutoff=0.01,
permutations=99999)
g$lc_PC2 <- factor(lisa_clusters(lc_PC2), levels=0:4,
labels=lisa_labels(lc_PC2)[1:5])
is.na(g$lc_PC2) <- g$lc_PC2 == "Not significant"
g$lc_PC2 <- droplevels(g$lc_PC2)
if (tmap_ok) tm_shape(g) + tm_fill("lc_PC2", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 9
else plot(g["lc_PC2"])
}
if (run) {
lmc_PC <- local_multigeary(w, g[c("PC1","PC2")], significance_cutoff=0.00247,
permutations=99999)
g$lmc_PC <- factor(lisa_clusters(lmc_PC), levels=0:1,
labels=lisa_labels(lmc_PC)[1:2])
is.na(g$lmc_PC) <- g$lmc_PC == "Not significant"
g$lmc_PC <- droplevels(g$lmc_PC)
table(interaction((p.adjust(lisa_pvalues(lmc_PC), "fdr") < 0.01), g$lmc_PC))
}
if (run) {
if (tmap_ok) tm_shape(g) + tm_fill("lmc_PC", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 10
else plot(g["lmc_PC"])
}
if (run) {
set.seed(1)
lmc_PC_spdep <- localC_perm(g[c("PC1","PC2")], lw, nsim=9999, alternative="two.sided")
all.equal(lisa_values(lmc_PC), c(lmc_PC_spdep))
}
if (run) {
cor(attr(lmc_PC_spdep, "pseudo-p")[,6], lisa_pvalues(lmc_PC))
}
if (run) {
g$lmc_PC_spdep <- attr(lmc_PC_spdep, "cluster")
is.na(g$lmc_PC_spdep) <- p.adjust(attr(lmc_PC_spdep, "pseudo-p")[,6], "fdr") > 0.01
g$lmc_PC_spdep <- droplevels(g$lmc_PC_spdep)
if (tmap_ok) tm_shape(g) + tm_fill("lmc_PC_spdep", textNA="Insignificant",
colorNA="gray95") + tm_borders() # rep. Fig. 10
else plot(g["lmc_PC_spdep"])
}
if (run) {
lmc_vars <- local_multigeary(w, st_drop_geometry(g)[, 1:6],
significance_cutoff=0.00247, permutations=99999)
g$lmc_vars <- factor(lisa_clusters(lmc_vars), levels=0:1,
labels=lisa_labels(lmc_vars)[1:2])
is.na(g$lmc_vars) <- g$lmc_vars == "Not significant"
g$lmc_vars <- droplevels(g$lmc_vars)
table(interaction((p.adjust(lisa_pvalues(lmc_vars), "fdr") < 0.01),
g$lmc_vars))
}
if (run) {
if (tmap_ok) tm_shape(g) + tm_fill("lmc_vars", textNA="Insignificant",
colorNA="gray95") + tm_borders() # Fig. 11
else plot(g["lmc_vars"])
}
if (run) {
set.seed(1)
system.time(lmc_vars_spdep <- localC_perm(st_drop_geometry(g)[, 1:6], lw,
nsim=9999, alternative="two.sided"))
}
if (run) {
all.equal(lisa_values(lmc_vars), c(lmc_vars_spdep))
}
if (run) {
cor(attr(lmc_vars_spdep, "pseudo-p")[,6], lisa_pvalues(lmc_vars))
}
if (run) {
if (require(parallel, quietly=TRUE)) {
ncpus <- detectCores()-1L
cores <- get.coresOption()
set.coresOption(ncpus)
system.time(lmc_vars_spdep1 <- localC_perm(st_drop_geometry(g)[, 1:6], lw,
nsim=9999, alternative="two.sided", iseed=1))
set.coresOption(cores)
}
}
if (run) {
all.equal(lisa_values(lmc_vars), c(lmc_vars_spdep1))
}
if (run) {
cor(attr(lmc_vars_spdep1, "pseudo-p")[,6], lisa_pvalues(lmc_vars))
}
if (run) {
g$lmc_vars_spdep <- attr(lmc_vars_spdep1, "cluster")
is.na(g$lmc_vars_spdep) <- p.adjust(attr(lmc_vars_spdep1, "pseudo-p")[,6], "fdr") > 0.01
g$lmc_vars_spdep <- droplevels(g$lmc_vars_spdep)
if (tmap_ok) tm_shape(g) + tm_fill("lmc_vars_spdep", textNA="Insignificant",
colorNA="gray95") + tm_borders() # rep. Fig. 11
else plot(g["lmc_vars_spdep"])
}
}
if (FALSE) {
library(reticulate)
use_python("/usr/bin/python", required = TRUE)
gp <- import("geopandas")
ps <- import("libpysal")
W <- listw2mat(listw)
w <- ps$weights$full2W(W, rownames(W))
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
lC <- esda$Geary_Local(connectivity=w)$fit(scale(x))
# np$std missing ddof=1
n <- length(x)
D0 <- spdep:::geary.intern((x - mean(x)) / sqrt(var(x)*(n-1)/n), listw, n=n)
# lC components probably wrongly ordered https://github.com/pysal/esda/issues/192
o <- match(round(D0, 6), round(lC$localG, 6))
all.equal(c(lC$localG)[o], D0)
# simulation order not retained
lC$p_sim[o]
attr(C, "pseudo-p")[,6]
}
Run the code above in your browser using DataLab