draft bivariate moran hotspots #155

R/local-moran-bv.R

@@ -6,10 +6,16 @@ local_moran_bv_calc <- function(x, y, listw) {
 
 
 localmoran_bv <- function(x, y, listw, nsim = 199, scale = TRUE,
-  alternative="two.sided", iseed=1L, no_repeat_in_row=FALSE) {
+  alternative="two.sided", iseed=1L, no_repeat_in_row=FALSE,
+  zero.policy=attr(listw, "zero.policy")) {
   stopifnot(length(x) == length(y))
   if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
     "is not a listw object"))
+  if (is.null(zero.policy))
+    zero.policy <- get("zeroPolicy", envir = .spdepOptions)
+  stopifnot(is.logical(zero.policy))
+  n <- length(listw$neighbours)
+  if (n != length(x)) stop("Different numbers of observations")
 # FIXME is listw assumed to be row-standardized?
   n <- length(listw$neighbours)
   stopifnot(n == length(x))
@@ -21,12 +27,28 @@ localmoran_bv <- function(x, y, listw, nsim = 199, scale = TRUE,
   if (missing(nsim)) stop("nsim must be given")
   stopifnot(all(!is.na(x)))
   stopifnot(all(!is.na(y)))
+
+  xx <- mean(x)
+  ly <- lag.listw(listw, y, zero.policy=zero.policy)
+  lyy <- mean(ly)
+  lbs <- c("Low", "High")
+  quadr <- interaction(cut(x, c(-Inf, xx, Inf), labels=lbs), 
+    cut(ly, c(-Inf, lyy, Inf), labels=lbs), sep="-")
+  xmed <- median(x)
+  lymed <- median(ly)
+  quadr_med <- interaction(cut(x, c(-Inf, xmed, Inf), labels=lbs),
+    cut(ly, c(-Inf, lymed, Inf), labels=lbs), sep="-")
 
   # the variables should be scaled and are by default
   if (scale) {
     x <- as.numeric(scale(x))
     y <- as.numeric(scale(y))
   }
+
+  ly <- lag.listw(listw, y, zero.policy=zero.policy)
+  quadr_ps <- interaction(cut(x, c(-Inf, 0, Inf), labels=lbs), 
+    cut(ly, c(-Inf, 0, Inf), labels=lbs), sep="-")
+
   cards <- card(listw$neighbours)
   stopifnot(all(cards > 0L))
 # FIXME no zero.policy handling
@@ -103,6 +125,8 @@ localmoran_bv <- function(x, y, listw, nsim = 199, scale = TRUE,
   Prname_sim <- "Pr(folded) Sim"
   colnames(res) <- c("Ibvi", "E.Ibvi", "Var.Ibvi", "Z.Ibvi", Prname,
     Prname_rank, Prname_sim)
+  attr(res, "quadr") <- data.frame(mean=quadr, median=quadr_med, 
+    pysal=quadr_ps)
   class(res) <- c("localmoran", class(res))
   res
 }

man/localmoran_bv.Rd

@@ -3,7 +3,7 @@
 \title{Compute the Local Bivariate Moran's I Statistic}
 \usage{
 localmoran_bv(x, y, listw, nsim = 199, scale = TRUE, alternative="two.sided",
- iseed=1L, no_repeat_in_row=FALSE)
+ iseed=1L, no_repeat_in_row=FALSE, zero.policy=attr(listw, "zero.policy"))
 }
 \arguments{
 \item{x}{a numeric vector of same length as \code{y}.}
@@ -19,7 +19,7 @@ localmoran_bv(x, y, listw, nsim = 199, scale = TRUE, alternative="two.sided",
 \item{alternative}{a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less".}
 \item{iseed}{default NULL, used to set the seed for possible parallel RNGs.}
 \item{no_repeat_in_row}{default \code{FALSE}, if \code{TRUE}, sample conditionally in each row without replacements to avoid duplicate values, \url{https://github.com/r-spatial/spdep/issues/124}}
-
+  \item{zero.policy}{default default \code{attr(listw, "zero.policy")} as set when \code{listw} was created, if attribute not set, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA}
 }
 \value{
 a \code{data.frame} containing two columns \code{Ib} and \code{p_sim} containing the local bivariate Moran's I and simulated p-values respectively.
@@ -42,6 +42,11 @@ nb <- poly2nb(columbus)
 listw <- nb2listw(nb)
 set.seed(1)
 (res <- localmoran_bv(columbus$CRIME, columbus$INC, listw, nsim = 499))
+columbus$hs <- hotspot(res, Prname="Pr(folded) Sim", cutoff=0.05,
+ quadrant.type="pysal", p.adjust="none")
+if (require("tmap", quietly=TRUE)) {
+  tm_shape(columbus) + tm_fill("hs")
+}
 }
 \references{
 Anselin, Luc, Ibnu Syabri, and Oleg Smirnov. 2002. “Visualizing Multivariate Spatial Correlation with Dynamically Linked Windows.” In New Tools for Spatial Data Analysis: Proceedings of the Specialist Meeting, edited by Luc Anselin and Sergio Rey. University of California, Santa Barbara: Center for Spatially Integrated Social Science (CSISS).