From 9d6228a92bead517a661096750ea63b191bd5128 Mon Sep 17 00:00:00 2001 From: Roger Bivand Date: Thu, 19 Dec 2024 09:56:50 +0100 Subject: [PATCH] re-instate rgeoda --- DESCRIPTION | 4 +- NEWS.md | 1 + man/localC.Rd | 467 +++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 469 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index addcabc..e1dc947 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: spdep Version: 1.3-9 -Date: 2024-12-02 +Date: 2024-12-20 Title: Spatial Dependence: Weighting Schemes, Statistics Encoding: UTF-8 Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), @@ -46,7 +46,7 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), person("Danlin", "Yu", role = "ctb")) Depends: R (>= 3.3.0), methods, spData (>= 2.3.1), sf Imports: stats, deldir, boot (>= 1.3-1), graphics, utils, grDevices, units, s2, e1071, sp (>= 1.0) -Suggests: spatialreg (>= 1.2-1), Matrix, parallel, dbscan, RColorBrewer, lattice, xtable, foreign, igraph, RSpectra, knitr, classInt, tmap, spam, ggplot2, rmarkdown, tinytest +Suggests: spatialreg (>= 1.2-1), Matrix, parallel, dbscan, RColorBrewer, lattice, xtable, foreign, igraph, RSpectra, knitr, classInt, tmap, spam, ggplot2, rmarkdown, tinytest, rgeoda (>= 0.0.11.1) URL: https://github.com/r-spatial/spdep/, https://r-spatial.github.io/spdep/ BugReports: https://github.com/r-spatial/spdep/issues/ Description: A collection of functions to create spatial weights matrix diff --git a/NEWS.md b/NEWS.md index 2ed14ee..a2ab357 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,6 @@ # Version 1.3-9 (development) +* re-instate **rgeoda** references # Version 1.3-8 (2024-12-02) diff --git a/man/localC.Rd b/man/localC.Rd index 58c5ba4..2365b87 100644 --- a/man/localC.Rd +++ b/man/localC.Rd @@ -83,7 +83,472 @@ x <- spData::afcon$totcon 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 +\donttest{ +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) { + tmap4 <- packageVersion("tmap") >= "3.99" + if (tmap4) { + tm_shape(g) + tm_polygons(fill="PC1", + fill.scale=tm_scale(values="brewer.rd_yl_gn", breaks=brks, + midpoint=0), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="PC2", + fill.scale=tm_scale(values="brewer.rd_yl_gn", breaks=brks, + midpoint=0), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lm_PC1", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lm_PC1_spdep", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lg_PC1", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lg_PC1a", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + tm_shape(g) + tm_fill("lg_PC1a", 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc_PC1", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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 <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L +# test with single core + if (ncpus > 1L) ncpus <- 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc_PC1_spdep", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="both_PC1", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc005_PC1", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc005_PC1_spdep", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc005_PC1", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc005_PC1_spdep", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lc_PC2", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lmc_PC", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lmc_PC_spdep", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lmc_vars", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + 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 <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L +# test with single core + if (ncpus > 1L) ncpus <- 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) { + if (tmap4) { + tm_shape(g) + tm_polygons(fill="lmc_vars_spdep", + fill.scale=tm_scale(values="brewer.set3", value.na="gray95", + label.na="Insignificant"), + fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), + frame=FALSE, item.r=0)) + } else { + tm_shape(g) + tm_fill("lmc_vars_spdep", textNA="Insignificant", + colorNA="gray95") + tm_borders() # rep. Fig. 11 + } + } else { + plot(g["lmc_vars_spdep"]) + } +} +} \dontrun{ library(reticulate)