Addition of two new functions: moran.plot.drop and moran.plot.seismogram + typo fix for moran.plot.Rd

R/moran.plot.drop.R

@@ -0,0 +1,115 @@
+moran.plot.drop <- function(x, listw, nsim = 999, cv = 2.58, significant = TRUE, plain = FALSE, zero.policy = FALSE, xlab = NULL, ylab = NULL, plot = TRUE, return_df = TRUE, spChk = NULL, labels = NULL) {
+  require(spdep)
+  if (!inherits(listw, "listw")) 
+    stop(paste(deparse(substitute(listw)), "is not a listw object"))
+  stopifnot(is.vector(x))
+  stopifnot(is.logical(significant))
+  stopifnot(is.logical(plain))
+  stopifnot(is.logical(return_df))
+  if (is.null(zero.policy)) 
+    zero.policy <- get("zeroPolicy", envir = .spdepOptions)
+  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")
+  if (is.null(spChk)) 
+    spChk <- get.spChkOption()
+  if (spChk && !chkIDs(x, listw)) 
+    stop("Check of data and weights ID integrity failed")
+  labs <- TRUE
+  if (is.logical(labels) && !labels) 
+    labs <- FALSE
+  if (is.null(labels) || length(labels) != n) 
+    labels <- as.character(attr(listw, "region.id"))
+  if (is.null(xlab)) 
+    xlab <- xname
+  if (is.null(ylab)) 
+    ylab <- paste("spatially lagged", xname)
+  Z <- as.vector(scale(x))
+  ZLXi <- lag.listw(listw, Z, zero.policy = zero.policy)
+  locmoran <- localmoran_perm(Z, listw, nsim = nsim)
+  ZIi <- locmoran[, 4]
+  b <- ((cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(Z) / (Z - mean(Z))
+  b2 <- ((-cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(Z) / (Z - mean(Z))
+  b[which((Z < 0 & ZLXi > 0) | (Z > 0 & ZLXi < 0))] <- b2[which((Z < 0 & ZLXi > 0) | (Z > 0 & ZLXi < 0))]
+  if(plot) {
+    if(!plain) {
+      if(significant) {
+        x_q1 <- Z[which(Z > 0 & ZLXi > 0 & ZIi >= cv)]
+        y_q1 <- ZLXi[which(Z > 0 & ZLXi > 0 & ZIi >= cv)]
+        b_q1 <- b[which(Z > 0 & ZLXi > 0 & ZIi >= cv)]
+        x_q2 <- Z[which(Z > 0 & ZLXi < 0 & ZIi <= (-1) * cv)]
+        y_q2 <- ZLXi[which(Z > 0 & ZLXi < 0 & ZIi <= (-1) * cv)]
+        b_q2 <- b[which(Z > 0 & ZLXi < 0 & ZIi <= (-1) * cv)]
+        x_q3 <- Z[which(Z < 0 & ZLXi < 0 & ZIi >= cv)]
+        y_q3 <- ZLXi[which(Z < 0 & ZLXi < 0 & ZIi >= cv)]
+        b_q3 <- b[which(Z < 0 & ZLXi < 0 & ZIi >= cv)]
+        x_q4 <- Z[which(Z < 0 & ZLXi > 0 & ZIi <= (-1) * cv)]
+        y_q4 <- ZLXi[which(Z < 0 & ZLXi > 0 & ZIi <= (-1) * cv)]
+        b_q4 <- b[which(Z < 0 & ZLXi > 0 & ZIi <= (-1) * cv)]
+      } else {
+        x_q1 <- Z[which(Z > 0 & ZLXi > 0 & ZIi < cv)]
+        y_q1 <- ZLXi[which(Z > 0 & ZLXi > 0 & ZIi < cv)]
+        b_q1 <- b[which(Z > 0 & ZLXi > 0 & ZIi < cv)]
+        x_q2 <- Z[which(Z > 0 & ZLXi < 0 & ZIi > (-1) * cv)]
+        y_q2 <- ZLXi[which(Z > 0 & ZLXi < 0 & ZIi > (-1) * cv)]
+        b_q2 <- b[which(Z > 0 & ZLXi < 0 & ZIi > (-1) * cv)]
+        x_q3 <- Z[which(Z < 0 & ZLXi < 0 & ZIi < cv)]
+        y_q3 <- ZLXi[which(Z < 0 & ZLXi < 0 & ZIi < cv)]
+        b_q3 <- b[which(Z < 0 & ZLXi < 0 & ZIi < cv)]
+        x_q4 <- Z[which(Z < 0 & ZLXi > 0 & ZIi > (-1) * cv)]
+        y_q4 <- ZLXi[which(Z < 0 & ZLXi > 0 & ZIi > (-1) * cv)]
+        b_q4 <- b[which(Z < 0 & ZLXi > 0 & ZIi > (-1) * cv)]
+      }
+    }
+    lw.lm <- lm(ZLXi ~ Z)
+    plot(Z, ZLXi, xlab="Z", ylab="WZ", pch = 20, cex = 0.33, col = "lightgrey", xlim = c(min(Z),max(Z)), ylim = c(min(ZLXi, b),max(ZLXi, b)))
+    abline(h = 0, lty = "dashed", col = "grey30")
+    abline(v = 0, lty = "dashed", col = "grey30")
+    abline(lw.lm, lty = "dotted", col = "grey40")
+    if(!plain) {
+      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 && length(x_q1) > 0) 
+          text(x_q1[i], y_q1[i], labels = labels[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  && length(x_q2) > 0) 
+          text(x_q2[i], y_q2[i], labels = labels[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  && length(x_q3) > 0) 
+          text(x_q3[i], y_q3[i], labels = labels[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  && length(x_q4) > 0) 
+          text(x_q4[i], y_q4[i], labels = labels[i], pos = 2, cex = 0.5, col = "royalblue")
+      }
+      if (zero.policy) {
+        n0 <- ZLXi == 0
+        if (any(n0)) {
+          symbols(x[n0], wx[n0], inches = FALSE, circles = rep(diff(range(x))/50, length(which(n0))), bg = "grey", add = TRUE)
+        }
+      }
+    }
+  }
+  if(return_df) {
+    res <- data.frame(z = Z, wz = ZLXi, b = b, line_lengths = abs(ZLXi - b))
+    invisible(res)
+  }
+}

R/moran.plot.seismogram.R

@@ -0,0 +1,85 @@
+moran.plot.seismogram <- function(x, listw, nsim = 999, cv = 2.58, plain = FALSE, zero.policy = FALSE, xlab = NULL, ylab = NULL, plot = TRUE, return_df = TRUE, spChk = NULL) {
+  require(spdep)
+  if (!inherits(listw, "listw")) 
+    stop(paste(deparse(substitute(listw)), "is not a listw object"))
+  stopifnot(is.vector(x))
+  stopifnot(is.logical(plain))
+  stopifnot(is.logical(return_df))
+  if (is.null(zero.policy)) 
+    zero.policy <- get("zeroPolicy", envir = .spdepOptions)
+  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")
+  if (is.null(spChk)) 
+    spChk <- get.spChkOption()
+  if (spChk && !chkIDs(x, listw)) 
+    stop("Check of data and weights ID integrity failed")
+  if (is.null(xlab)) 
+    xlab <- xname
+  if (is.null(ylab)) 
+    ylab <- paste("spatially lagged", xname)
+  Z <- as.vector(scale(x))
+  ZLXi <- lag.listw(listw, Z, zero.policy = zero.policy)
+  locmoran <- localmoran_perm(Z, listw, nsim = nsim)
+  ZIi <- locmoran[, 4]
+  b <- ((cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(Z) / (Z - mean(Z))
+  b2 <- ((-cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(Z) / (Z - mean(Z))
+  b[which((Z < 0 & ZLXi > 0) | (Z > 0 & ZLXi < 0))] <- b2[which((Z < 0 & ZLXi > 0) | (Z > 0 & ZLXi < 0))]
+  if(plot) {
+    if(!plain) {
+      x_q1 <- Z[which(Z > 0 & ZLXi > 0)]
+      y_q1 <- ZLXi[which(Z > 0 & ZLXi > 0)]
+      b_q1 <- b[which(Z > 0 & ZLXi > 0)]
+      x_q2 <- Z[which(Z > 0 & ZLXi < 0)]
+      y_q2 <- ZLXi[which(Z > 0 & ZLXi < 0)]
+      b_q2 <- b[which(Z > 0 & ZLXi < 0)]
+      x_q3 <- Z[which(Z < 0 & ZLXi < 0)]
+      y_q3 <- ZLXi[which(Z < 0 & ZLXi < 0)]
+      b_q3 <- b[which(Z < 0 & ZLXi < 0)]
+      x_q4 <- Z[which(Z < 0 & ZLXi > 0)]
+      y_q4 <- ZLXi[which(Z < 0 & ZLXi > 0)]
+      b_q4 <- b[which(Z < 0 & ZLXi > 0)]
+    }
+    lw.lm <- lm(ZLXi ~ Z)
+    plot(Z, ZLXi, xlab="Z", ylab="WZ", pch = 20, cex = 0.33, col = "lightgrey", xlim = c(min(Z),max(Z)), ylim = c(min(ZLXi, b),max(ZLXi, b)))
+    abline(h = 0, lty = "dashed", col = "grey30")
+    abline(v = 0, lty = "dashed", col = "grey30")
+    abline(lw.lm, lty = "dotted", col = "grey40")
+    if(!plain) {
+      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(z = Z, wz = ZLXi, b = b)
+    invisible(res)
+  }
+}

man/moran.plot.Rd

@@ -18,7 +18,7 @@ moran.plot(x, listw, zero.policy=attr(listw, "zero.policy"), spChk=NULL, labels=
   \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, nsim = 999, cv = 2.58,
+  significant = TRUE, plain = FALSE, zero.policy = FALSE,
+  xlab = NULL, ylab = NULL, plot = TRUE, return_df = TRUE,
+  spChk = NULL, labels = FALSE)
+}
+
+\arguments{
+  \item{x}{a numerical vector holding the attribute of interest}
+  \item{listw}{a \code{listw} spatial weights object}
+  \item{nsim}{default 999; number of conditonal permutation simulations}
+  \item{cv}{default 2.58; the desired critical value assuming approximate normality of standardised local Moran's \emph{I} values}
+  \item{significant}{default TRUE; a parameter indicating whether to display plot distances of significant (default) or non-significant observations}
+  \item{plain}{default FALSE; a plain Moran scatterplot is displayed if TRUE}
+  \item{zero.policy}{default NULL; use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors}
+  \item{xlab}{label for x axis}
+  \item{ylab}{label for y axis}
+  \item{plot}{default TRUE; if FALSE, plotting is suppressed}
+  \item{return_df}{default TRUE; invisibly return a data.frame object, if FALSE do not return anything}
+  \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}{default FALSE; no labels are plotted by default; character labels for points are assigned to significant observations if provided, region IDs are used if set to TRUE}
+}
+
+\details{
+The Moran drop plot is a version of the Moran scatterplot supplemented by indications of \emph{p} values. The standard Moran scatterplot provides indirect information about the effect size (distance from the trend line), 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 under the assumption of an approximate standard normality of the z-scores of the local Moran's \emph{I} values. The longer the lines, the higher 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{z}{the standardised attribute values}
+  \item{wz}{the standardised spatially lagged attribute values}
+  \item{b}{the y-coordinates of the critical values}
+  \item{line_lengths}{the absolute distances between b and z}
+}
+
+\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.shp", package="spData")[1])
+boston.nb <- poly2nb(boston.tr)
+boston.listw <- nb2listw(boston.nb)
+moran.plot.drop(boston.c$CMEDV, boston.listw, 999, 2.58, zero.policy = TRUE, significant = TRUE, plain = FALSE, labels = FALSE, plot = TRUE)
+}

man/moran.plot.seismogram.Rd

@@ -0,0 +1,53 @@
+\name{moran.plot.seismogram}
+\alias{moran.plot.seismogram}
+
+\title{Moran seismogram}
+
+\description{
+A variant of the Moran scatterplot (see \code{\link{moran.plot}}) supplemented by lines connecting location-wise critical values. The plot allows for visual inspection of potential spatial weights misspecifiation.
+}
+\usage{
+moran.plot.seismogram(x, listw, nsim = 999, cv = 2.58,
+  plain = FALSE, zero.policy = FALSE, xlab = NULL,
+  ylab = NULL, plot = TRUE, return_df = TRUE, spChk = NULL)
+}
+
+\arguments{
+  \item{x}{a numerical vector holding the attribute of interest}
+  \item{listw}{a \code{listw} spatial weights object}
+  \item{nsim}{default 999; number of conditonal permutation simulations}
+  \item{cv}{default 2.58; the desired critical value assuming approximate normality of standardised local Moran's \emph{I} values}
+  \item{plain}{default FALSE; a plain Moran scatterplot is displayed if TRUE}
+  \item{zero.policy}{default NULL; use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors}
+  \item{xlab}{label for x axis}
+  \item{ylab}{label for y axis}
+  \item{plot}{default TRUE; if FALSE, plotting is suppressed}
+  \item{return_df}{default TRUE; invisibly return a data.frame object, if FALSE do not return anything}
+  \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()}}
+}
+
+\details{
+The Moran seismogram is a version of the Moran scatterplot that is complemented by lines connecting the location-specific critical values of local Moran's \emph{I}. Lines in quadrants with positive spatial autocorrelation are shown in red and lines in quadrants with negative spatial autocorrelation are shown in blue. The representation is similar to a seismogram for detecting earthquakes and thus reveals potentially suspicious local spatial weights configurations by visualising spikes. The latter are displayed in an integrated manner with their positions in the attribute value range and in connection with the types of the associated spatial patterns (by the quadrants of the scatterplot).
+}
+
+\value{
+When return_df is TRUE, a data frame object with the following members is returned:
+  \item{z}{the standardised attribute values}
+  \item{wz}{the standardised spatially lagged attribute values}
+  \item{b}{the y-coordinates of the critical values}
+}
+
+\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.shp", package="spData")[1])
+boston.nb <- poly2nb(boston.tr)
+boston.listw <- nb2listw(boston.nb)
+moran.plot.seismogram(boston.c$CMEDV, boston.listw, 999, 2.58, zero.policy = TRUE, plain = FALSE, plot = TRUE)
+}