Addition of moran.plot.drop and moran.plot.seismogram

NAMESPACE

@@ -45,8 +45,9 @@ export(listw2sn, sn2listw, read.gwt2nb,	write.sn2gwt,
         read.swmdbf2listw, read_swm_dbf, write.swmdbf, write_swm_dbf, write.sn2DBF,
         lm.LMtests, lm.RStests,
 	lm.morantest, localG, localG_perm, localmoran, localmoran_perm, moran,
-	moran.test, moran.mc, moran.plot, localmoran.sad, lm.morantest.sad,
-	nb2listw, nb2listwdist, nb2mat, listw2mat, mat2listw, nbdists, nblag,
+	moran.test, moran.mc, moran.plot, moran.plot.drop, moran.plot.seismogram,
+	localmoran.sad, lm.morantest.sad, nb2listw, nb2listwdist, nb2mat,
+	listw2mat, mat2listw, nbdists, nblag,
 	nblag_cumul, poly2nb, read.gal, write.nb.gal, read.geoda, relativeneigh,
 	soi.graph, sp.correlogram, sp.mantel.mc, set.spChkOption, chkIDs,
 	get.spChkOption, spNamedVec, tri2nb,

R/moran.plot.drop.R

@@ -0,0 +1,133 @@
+moran.plot.drop <- function(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL, significant = TRUE, xlab = NULL, ylab = NULL, return_df = TRUE, spChk = NULL, labels = NULL, zero.policy=attr(listw, "zero.policy")) {
+  if (!inherits(listw, "listw")) 
+    stop(paste(deparse(substitute(listw)), "is not a listw object"))
+  if (!inherits(locmoran, "localmoran")) 
+    stop(paste(deparse(substitute(locmoran)), "is not a localmoran object"))
+  stopifnot(is.vector(x))
+  stopifnot(is.logical(significant))
+  stopifnot(is.logical(return_df))
+  stopifnot(is.numeric(alpha))
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  xname <- deparse(substitute(x))
+  if (!is.numeric(x)) 
+    stop(paste(xname, "is not a numeric vector"))
+  if (anyNA(x)) 
+    stop("NA in X")
+  n <- length(listw$neighbours)
+  if (n != length(x)) 
+    stop("objects of different length")
+  usePadj <- FALSE
+  if (!is.null(adjusted_p)) {
+    if(n != length(adjusted_p))
+      stop("objects of different length")
+    if (anyNA(adjusted_p)) 
+      stop("NA in vector of adjusted p values")
+    usePadj <- TRUE
+  }
+  if (is.null(spChk)) 
+    spChk <- get.spChkOption()
+  if (spChk && !chkIDs(locmoran, listw)) 
+    stop("Check of data and weights ID integrity failed")
+  labs <- TRUE
+  if (is.logical(labels)) {
+    if(!labels)
+      labs <- FALSE
+    labels <- as.character(attr(listw, "region.id"))
+  } else if (!is.logical(labels) && !is.null(labels)) {
+    if(length(labels) != n) {
+      warning("Length of the labels vector does not match number of regions. region.id from the listw object is used instead.")
+      labels <- as.character(attr(listw, "region.id"))
+    }
+  } else if (is.null(labels)) {
+    labs <- FALSE
+    labels <- as.character(attr(listw, "region.id"))
+  }
+  if (is.null(xlab))
+    xlab <- xname
+  if (is.null(ylab)) 
+    ylab <- paste("spatially lagged centred", xname)
+
+  WX <- lag.listw(listw, scale(x, scale = F), zero.policy = zero.policy)
+  if (anyNA(WX)) warning("no-neighbour observation(s) found - use zero.policy=TRUE")
+  if(usePadj)
+    cv <- min(abs(locmoran[(which(adjusted_p <= alpha/2)), 4]))
+  else
+    cv <- abs(qnorm(1 - alpha, lower.tail = FALSE))
+
+  b <- ((cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b2 <- ((-cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))] <- b2[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))]
+
+  if(significant) {
+    x_q1 <- x[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    y_q1 <- WX[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    b_q1 <- b[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    labels_q1 <- labels[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    x_q2 <- x[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    y_q2 <- WX[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    b_q2 <- b[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    labels_q2 <- labels[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    x_q3 <- x[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    y_q3 <- WX[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    b_q3 <- b[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    labels_q3 <- labels[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    x_q4 <- x[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    y_q4 <- WX[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    b_q4 <- b[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    labels_q4 <- labels[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+  } else {
+    x_q1 <- x[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] < cv)]
+    y_q1 <- WX[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] < cv)]
+    b_q1 <- b[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] < cv)]
+    x_q2 <- x[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] > (-1) * cv)]
+    y_q2 <- WX[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] > (-1) * cv)]
+    b_q2 <- b[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] > (-1) * cv)]
+    x_q3 <- x[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] < cv)]
+    y_q3 <- WX[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] < cv)]
+    b_q3 <- b[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] < cv)]
+    x_q4 <- x[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] > (-1) * cv)]
+    y_q4 <- WX[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] > (-1) * cv)]
+    b_q4 <- b[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] > (-1) * cv)]
+  }
+
+  lw.lm <- lm(WX ~ x)
+  plot(x, WX, xlab=xlab, ylab=ylab, pch = 20, cex = 0.5, col = "gray70", xlim = c(min(x),max(x)), ylim = c(min(WX, b),max(WX, b)))
+  abline(h = mean(WX), lty = "dashed", col = "grey30")
+  abline(v = mean(x), lty = "dashed", col = "grey30")
+  abline(lw.lm, lty = "dotted", col = "grey40")
+
+  for(i in 1:length(x_q1)) {
+    lines(c(x_q1[i], x_q1[i]), c(y_q1[i], b_q1[i]), type = "l", col = "lightpink", lwd = 1)
+    points(c(x_q1[i]), c(y_q1[i]), pch = 20, cex = 0.5, col = "firebrick")
+    if (labs) 
+      text(x_q1[i], y_q1[i], labels = labels_q1[i], pos = 2, cex = 0.5, col = "firebrick")
+  }
+
+  for(i in 1:length(x_q2)) {
+    lines(c(x_q2[i], x_q2[i]), c(y_q2[i], b_q2[i]), type = "l", col = "lightblue1", lwd = 1)
+    points(c(x_q2[i]), c(y_q2[i]), pch = 20, cex = 0.5, col = "royalblue")
+    if (labs) 
+      text(x_q2[i], y_q2[i], labels = labels_q2[i], pos = 2, cex = 0.5, col = "royalblue")
+  }
+
+  for(i in 1:length(x_q3)) {
+    lines(c(x_q3[i], x_q3[i]), c(y_q3[i], b_q3[i]), type = "l", col = "lightpink", lwd = 1)
+    points(c(x_q3[i]), c(y_q3[i]), pch = 20, cex = 0.5, col = "firebrick")
+    if (labs) 
+      text(x_q3[i], y_q3[i], labels = labels_q3[i], pos = 2, cex = 0.5, col = "firebrick")
+  }
+
+  for(i in 1:length(x_q4)) {
+    lines(c(x_q4[i], x_q4[i]), c(y_q4[i], b_q4[i]), type = "l", col = "lightblue1", lwd = 1)
+    points(c(x_q4[i]), c(y_q4[i]), pch = 20, cex = 0.5, col = "royalblue")
+    if (labs) 
+      text(x_q4[i], y_q4[i], labels = labels_q4[i], pos = 2, cex = 0.5, col = "royalblue")
+  }
+
+  if(return_df) {
+    res <- data.frame(labels = labels, x = x, WX = WX, b = b, line_lengths = abs(WX - b))
+    invisible(res)
+  }
+}

R/moran.plot.seismogram.R

@@ -0,0 +1,91 @@
+moran.plot.seismogram <- function(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL, xlab = NULL, ylab = NULL, return_df = TRUE, spChk = NULL, zero.policy=attr(listw, "zero.policy")) {
+  if (!inherits(listw, "listw")) 
+    stop(paste(deparse(substitute(listw)), "is not a listw object"))
+  if (!inherits(locmoran, "localmoran")) 
+    stop(paste(deparse(substitute(locmoran)), "is not a localmoran object"))
+  stopifnot(is.vector(x))
+  stopifnot(is.logical(return_df))
+  stopifnot(is.numeric(alpha))
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  xname <- deparse(substitute(x))
+  if (!is.numeric(x)) 
+    stop(paste(xname, "is not a numeric vector"))
+  if (any(is.na(x))) 
+    stop("NA in X")
+  n <- length(listw$neighbours)
+  if (n != length(x)) 
+    stop("objects of different length")
+  usePadj <- FALSE
+  if (!is.null(adjusted_p)) {
+    if(n != length(adjusted_p))
+      stop("objects of different length")
+    if (any(is.na(adjusted_p))) 
+      stop("NA in vector of adjusted p values")
+    usePadj <- TRUE
+  }
+  if (is.null(spChk)) 
+    spChk <- get.spChkOption()
+  if (spChk && !chkIDs(locmoran, listw)) 
+    stop("Check of locmoran and weights ID integrity failed")
+  if (is.null(xlab)) 
+    xlab <- xname
+  if (is.null(ylab)) 
+    ylab <- paste("spatially lagged centred", xname)
+
+  WX <- lag.listw(listw, scale(x, scale = F), zero.policy = zero.policy)
+  if (anyNA(WX)) warning("no-neighbour observation(s) found - use zero.policy=TRUE")
+  if(usePadj)
+    cv <- min(abs(locmoran[(which(adjusted_p <= alpha/2)),4]))
+  else
+    cv <- abs(qnorm(1 - alpha, lower.tail = FALSE))
+  b <- ((cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b2 <- ((-cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))] <- b2[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))]
+
+  x_q1 <- x[which(x > mean(x) & WX > mean(WX))]
+  y_q1 <- WX[which(x > mean(x) & WX > mean(WX))]
+  b_q1 <- b[which(x > mean(x) & WX > mean(WX))]
+  x_q2 <- x[which(x > mean(x) & WX < mean(WX))]
+  y_q2 <- WX[which(x > mean(x) & WX < mean(WX))]
+  b_q2 <- b[which(x > mean(x) & WX < mean(WX))]
+  x_q3 <- x[which(x < mean(x) & WX < mean(WX))]
+  y_q3 <- WX[which(x < mean(x) & WX < mean(WX))]
+  b_q3 <- b[which(x < mean(x) & WX < mean(WX))]
+  x_q4 <- x[which(x < mean(x) & WX > mean(WX))]
+  y_q4 <- WX[which(x < mean(x) & WX > mean(WX))]
+  b_q4 <- b[which(x < mean(x) & WX > mean(WX))]
+
+  lw.lm <- lm(WX ~ x)
+  plot(x, WX, xlab=xlab, ylab=ylab, pch = 20, cex = 0.5, col = "gray70", xlim = c(min(x),max(x)), ylim = c(min(WX, b),max(WX, b)))
+  abline(h = mean(WX), lty = "dashed", col = "grey30")
+  abline(v = mean(x), lty = "dashed", col = "grey30")
+  abline(lw.lm, lty = "dotted", col = "grey40")
+
+  df_q1 <- data.frame(x_q1, b_q1)
+  df_q1 <- df_q1[order(x_q1),]
+  df_q2 <- data.frame(x_q2, b_q2)
+  df_q2 <- df_q2[order(x_q2),]
+  df_q3 <- data.frame(x_q3, b_q3)
+  df_q3 <- df_q3[order(x_q3),]
+  df_q4 <- data.frame(x_q4, b_q4)
+  df_q4 <- df_q4[order(x_q4),]
+
+  for(i in 1:(length(df_q1$x_q1) - 1))
+    lines(c(df_q1$x_q1[i], df_q1$x_q1[i+1]), c(df_q1$b_q1[i], df_q1$b_q1[i+1]), type = "l", col = "firebrick")
+
+  for(i in 1:length(df_q2$x_q2)) 
+    lines(c(df_q2$x_q2[i], df_q2$x_q2[i+1]), c(df_q2$b_q2[i], df_q2$b_q2[i+1]), type = "l", col = "royalblue")
+
+  for(i in 1:length(df_q3$x_q3))
+    lines(c(df_q3$x_q3[i], df_q3$x_q3[i+1]), c(df_q3$b_q3[i], df_q3$b_q3[i+1]), type = "l", col = "firebrick")
+
+  for(i in 1:length(df_q4$x_q4))
+    lines(c(df_q4$x_q4[i], df_q4$x_q4[i+1]), c(df_q4$b_q4[i], df_q4$b_q4[i+1]), type = "l", col = "royalblue")
+
+  if(return_df) {
+    res <- data.frame(labels=as.character(attr(listw, "region.id")), x=x, wx=WX, b=b)
+    invisible(res)
+  }
+}

man/moran.plot.Rd

@@ -19,7 +19,7 @@ moran.plot(x, listw, y=NULL, zero.policy=attr(listw, "zero.policy"), spChk=NULL,
   \item{spChk}{should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use \code{get.spChkOption()}}
   \item{labels}{character labels for points with high influence measures, if set to FALSE, no labels are plotted for points with large influence}
   \item{xlab}{label for x axis}
-  \item{ylab}{label for x axis}
+  \item{ylab}{label for y axis}
   \item{quiet}{default NULL, use !verbose global option value; if TRUE, output of summary of influence object suppressed}
   \item{plot}{default TRUE, if false, plotting is suppressed}
   \item{return_df}{default TRUE, invisibly return a data.frame object; if FALSE invisibly return an influence measures object}

man/moran.plot.drop.Rd

@@ -0,0 +1,57 @@
+\name{moran.plot.drop}
+\alias{moran.plot.drop}
+
+\title{Moran drop plot}
+
+\description{
+A version of the Moran scatterplot (see \code{\link{moran.plot}}) supplemented by lines indicating \emph{p} values for visual inspection of statistical significance.
+}
+\usage{
+moran.plot.drop(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL,
+  significant = TRUE, xlab = NULL, ylab = NULL, return_df = TRUE,
+  spChk = NULL, labels = NULL, zero.policy=attr(listw, "zero.policy"))
+}
+
+\arguments{
+  \item{x}{a numerical vector holding the attribute of interest}
+  \item{listw}{a \code{listw} spatial weights object}
+  \item{locmoran}{a fitted object of type \code{localmoran}}
+  \item{alpha}{default 0.05; the desired significance level regarding local Moran's \emph{I}}
+  \item{adjusted_p}{default NULL; an optional vector of \emph{p} values adjusted to account for multiple testing as is returned by \code{\link{p.adjustSP}}; if NULL, standard normal approximation is used to determine critical values}
+  \item{significant}{default TRUE; a parameter indicating whether to display critical value distances of significant (default) or non-significant observations}
+  \item{xlab}{default NULL; an optional label for the x-axis}
+  \item{ylab}{default NULL; an optional label for the y-axis}
+  \item{return_df}{default TRUE; invisibly return a data.frame object, if FALSE do not return anything}
+  \item{spChk}{default NULL to use \code{get.spChkOption()}; should the locmoran names be checked against the listw spatial objects for identity integrity, TRUE, or FALSE}
+  \item{labels}{default NULL; no labels are plotted by default; region IDs from the \code{listw} object are used as labels of significant observations if set to TRUE; custom labels are used if a character vector is provided}
+  \item{zero.policy}{default option stored in the listw object; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors}
+}
+
+\details{
+The Moran drop plot is a version of the Moran scatterplot supplemented by visual indications of \emph{p} values. The standard Moran scatterplot provides information about the effect size but not about the level of confidence to determine whether the effects shown might be more than just random outcomes. The Moran drop plot marks significant points in red (positive) and blue (negative spatial autocorrelation) and adds so-called drop lines connecting the significant observations to the scatterplot positions of their associated critical values. The coordinates of the latter are determined either under the assumption of approximate standard normality of the z-scores of the local Moran's \emph{I} values, or based on \emph{p} values provided through \code{adjusted_p}. In the latter case, the critical value is approximated based on the observation with the highest adjusted \emph{p} value that still satisfies the selected significance level. The longer the lines, the lower the associated \emph{p} values. The visualisation thus enables visual inspection of statistical significance while maintaining the relationship to both attribute value ranges and scatterplot quadrants. It is also possible to invert the visualised relationship and display the distances of non-significant observations to their corresponding critical values (if significant is set to FALSE).
+}
+
+\value{
+When return_df is TRUE, a data frame object with the following members is returned:
+  \item{labels}{either the labels provided or the region ids, if not specified}
+  \item{x}{the attribute values}
+  \item{wx}{the spatial lags of the centred attribute values}
+  \item{b}{the y-coordinates (i.e. hypothetical spatial lags) of the critical values given x}
+  \item{line_lengths}{the absolute distances between b and x}
+}
+
+\references{Westerholt, R. (2024): Extending the Moran scatterplot by indications of critical values and \emph{p}-values: introducing the Moran seismogram and the drop plot. Proceedings of the 32nd Annual GIS Research UK Conference (GISRUK), Leeds, UK. \url{https://doi.org/10.5281/zenodo.10897792}}
+
+\author{Rene Westerholt \email{rene.westerholt@tu-dortmund.de}}
+
+\seealso{\code{\link{moran.plot}}}
+
+\examples{
+# Boston example (CMEDV; owner-occupied housing in USD)
+data(boston)
+boston.tr <- sf::st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1])
+boston.nb <- poly2nb(boston.tr)
+boston.listw <- nb2listw(boston.nb)
+moran.plot.drop(boston.c$CMEDV, boston.listw, localmoran(boston.c$CMEDV, boston.listw), 0.01,
+  significant = TRUE, labels = NULL)
+}