Nothing
#' @title Autocorrelation Plot
#' @description Autocorrelation plot (Anselin 1996), following Chen (2015),
#' aka, Moran's-I plot (univariate or bivariate)
#'
#' @param x Vector of x response variables
#' @param y Vector of y response variables
#' @param coords A matrix of coordinates corresponding to [x,y]
#' @param type.ac Type of autocorrelation plot ("xy", "yx")
#' @param dist.function ("inv.power", "neg.exponent")
#' @param scale.xy (TRUE/FALSE) scale the x,y vectors
#' @param scale.morans (FALSE/TRUE) standardize the Moran's index to an
#' expected [-1 to 1]?
#' @param ... Additional arguments passed to plot
#'
#' @details
#' The argument "type" controls the plot for x influencing y (type="xy") or y
#' influencing x (type="yx"). If y is not defined then the statistic is univariate
#' and only the "xy" plot will be available. The linear relationship between x and
#' its spatial lag (Wx) is indicative of the spatial autoregressive process, underlying
#' the spatial dependence. The statistic can be autocorrelation (univariate or
#' cross-correlation (bivariate). The quadrants are the zero intercept for random
#' autocorrelation and the red line represents the trend in autocorrelation. The quadrants
#' in the plot indicate the type of spatial association/interaction (Anselin 1996).
#' For example the upper-left quadrant represents negative associations of low values
#' surrounded by high and the lower-right quadrant represents negative associations of
#' high values surrounded by low.
#'
#' If y is not specified the univariate statistic for x is returned. the coords argument
#' is only used if k = NULL. Can also be an sp object with relevant x,y coordinate slot
#' (ie., points or polygons). If w = NULL, the default method for deriving spatial weights
#' matrix, options are: inverse power or negative exponent. If scale.xy = FALSE it is
#' assumed that they are already scaled following Chen (2015).
#'
#' @return A plot of the scaled variable against its spatially lagged values.
#'
#' @author Jeffrey S. Evans <jeffrey_evans@@tnc.org>
#'
#' @references
#' Chen., Y. (2015) A New Methodology of Spatial Cross-Correlation Analysis.
#' PLoS One 10(5):e0126158. doi:10.1371/journal.pone.0126158
#'
#' Anselin, L. (1996) The Moran scatterplot as an ESDA tool to assess local instability
#' in spatial association. pp. 111-125 in M. M. Fischer, H. J. Scholten and D. Unwin (eds)
#' Spatial analytical perspectives on GIS, London, Taylor and Francis
#'
#' Anselin, L. (1995) Local indicators of spatial association, Geographical Analysis,
#' 27:93-115
#'
#' @examples
#' p = c("sf", "sp", "spdep")
#' if(any(!unlist(lapply(p, requireNamespace, quietly=TRUE)))) {
#' m = which(!unlist(lapply(p, requireNamespace, quietly=TRUE)))
#' message("Can't run examples, please install ", paste(p[m], collapse = " "))
#' } else {
#' invisible(lapply(p, require, character.only=TRUE))
#'
#' data(meuse, package = "sp")
#' meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, agr = "constant")
#'
#' # Autocorrelation (univariate)
#' morans.plot(meuse$zinc, coords = st_coordinates(meuse)[,1:2])
#'
#' # Cross-correlation of: x influencing y and y influencing x
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(1,2))
#' morans.plot(x=meuse$zinc, y=meuse$copper, coords = st_coordinates(meuse)[,1:2],
#' scale.morans = TRUE)
#' morans.plot(x=meuse$zinc, y=meuse$copper, coords = st_coordinates(meuse)[,1:2],
#' scale.morans = TRUE, type.ac="yx")
#' par(opar)
#' }
#' @export morans.plot
morans.plot <- function(x, y = NULL, coords = NULL,
type.ac = c("xy", "yx"),
dist.function = "inv.power",
scale.xy = TRUE, scale.morans = FALSE, ...) {
if(missing(x)) stop("x must be specified")
if(is.null(y)) {
type.lab = paste0("Autocorrelation of ", deparse(substitute(x)))
y = x
type.ac = "xy"
message("Calculating autocorrelation for x")
} else {
type.lab = paste0("Cross-correlation of ", deparse(substitute(x)), " and ", deparse(substitute(y)))
message("Calculating cross-correlation for [x,y]")
}
if(length(y) != length(x)) stop("[X,Y] are not equal")
if( length(which(is.na(x))) != 0 | length(which(is.na(y))) != 0)
stop("NA's not permitted in [X,Y]")
dots <- as.list(match.call(expand.dots = TRUE)[-1])
if ( !is.null(dots[["x"]]) ) dots[["x"]] <- NULL
if ( !is.null(dots[["y"]]) ) dots[["y"]] <- NULL
if ( !is.null(dots[["coords"]]) ) dots[["coords"]] <- NULL
if ( !is.null(dots[["type.ac"]]) ) dots[["type.ac"]] <- NULL
if ( !is.null(dots[["dist.function"]]) ) dots[["dist.function"]] <- NULL
if ( !is.null(dots[["scale.xy"]]) ) dots[["scale.xy"]] <- NULL
if ( !is.null(dots[["scale.morans"]]) ) dots[["scale.morans"]] <- NULL
if(scale.xy == FALSE) warning("It is assumed that [x,y] vectors are already scaled")
type.ac = type.ac[1]
if( scale.xy ){
x.scale <- as.numeric(scale(x))
y.scale <- as.numeric(scale(y))
}
if( is.null(coords) ) stop("A coordinates matrix is required")
w <- sp::spDists( coords )
if( dist.function == "inv.power" ) {
message("Calculating spatial weights matrix using inverse power function")
w <- 1 / w
diag(w) <- 0
w <- w / sum(w)
} else if (dist.function == "neg.exponent") {
message("Calculating spatial weights matrix using negative exponent")
diag(w) <- NA
mu <- mean(w, na.rm=TRUE)
for(i in 1:nrow(w)) {
for(j in 1:nrow(w)) {
w[i,j] <- round(exp( (-2 * w[i,j]) / mu ),6)
}
}
diag(w) <- 0
} else {
stop("Not a valid matrix option")
}
# standardizes variable -1 to 1
scale.partial <- function(p) {
d <- (max(p) + min(p)) / 2
return((p - d) / (max(p) - d))
}
plot.ac <- function(x.plot,y.plot,dots) {
if (is.null(dots[["pch"]]) & "pch" %in% names(dots) == FALSE) dots[["pch"]] <- 20
if (is.null(dots[["cex"]]) & "cex" %in% names(dots) == FALSE) dots[["cex"]] <- 1
dots[["x"]] <- x.plot
dots[["y"]] <- y.plot
do.call("plot", dots)
}
if(type.ac[1] == "xy") {
if (is.null(dots[["xlab"]]) & "xlab" %in% names(dots) == FALSE) dots[["xlab"]] <- deparse(substitute(x))
if (is.null(dots[["ylab"]]) & "ylab" %in% names(dots) == FALSE) dots[["ylab"]] <- "nWy,xx^TWy"
if (is.null(dots[["main"]]) & "main" %in% names(dots) == FALSE) dots[["main"]] <- type.lab
# gsci.xy <- as.numeric(x*w%*%y)
gsci.xy <- as.numeric(x.scale*y.scale%*%w)
if(scale.morans) {
gsci.xy <- scale.partial(gsci.xy)
slag <- as.numeric(x.scale*y.scale%*%w)
slag <- scale.partial(slag)
} else {
slag <- as.numeric(x.scale*y.scale%*%w)
}
plot.ac(x.scale,gsci.xy,dots)
graphics::abline(v=0, lty=3)
graphics::abline(h=0, lty=3)
graphics::abline(stats::lm( slag ~ x.scale ), col="red")
} else if (type.ac[1] == "yx") {
type.lab <- paste0("Cross-correlation of ", deparse(substitute(y)), " and ", deparse(substitute(x)))
if (is.null(dots[["xlab"]]) & "xlab" %in% names(dots) == FALSE) dots[["xlab"]] <- deparse(substitute(y))
if (is.null(dots[["ylab"]]) & "ylab" %in% names(dots) == FALSE) dots[["ylab"]] <- "nWx,yy^TWx"
if (is.null(dots[["main"]]) & "main" %in% names(dots) == FALSE) dots[["main"]] <- type.lab
#gsci.xy <- as.numeric(y*w%*%x)
gsci.yx <- as.numeric(y.scale*x.scale%*%w)
if(scale.morans) {
gsci.yx <- scale.partial(gsci.yx)
slag <- as.numeric(x.scale*y.scale%*%w)
slag <- scale.partial(slag)
} else {
slag <- as.numeric(y.scale*x.scale%*%w)
}
plot.ac(y.scale,gsci.yx,dots)
graphics::abline(v=0, lty=3)
graphics::abline(h=0, lty=3)
graphics::abline(stats::lm( slag ~ y.scale ), col="red")
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.