R/esv.R

In spmodel: Spatial Statistical Modeling and Prediction

#' Compute the empirical semivariogram
#'
#' @description Compute the empirical semivariogram for varying bin sizes and
#'   cutoff values.
#'
#' @param formula A formula describing the fixed effect structure. \code{.} on
#'   the right-hand side represents every variable in \code{data} except the
#'   response and the x-coordinate/y-coordinate columns (\code{xcoord}/\code{ycoord},
#'   or, for an \code{sf} object, the geometry column), which are never
#'   included via \code{.} (though they may still be given explicitly).
#' @param data A data frame or \code{sf} object containing the variables in \code{formula}
#'   and geographic information.
#' @param xcoord Name of the variable in \code{data} representing the x-coordinate.
#'   Can be quoted or unquoted. Not required if \code{data} is an \code{sf} object.
#' @param ycoord Name of the variable in \code{data} representing the y-coordinate.
#'   Can be quoted or unquoted. Not required if \code{data} is an \code{sf} object.
#' @param cloud A logical indicating whether the empirical semivariogram should
#'   be summarized by distance class or not. When \code{cloud = FALSE} (the default), pairwise semivariances
#'   are binned and averaged within distance classes. When \code{cloud} = TRUE,
#'   all pairwise semivariances and distances are returned (this is known as
#'   the "cloud" semivariogram).
#' @param robust A logical indicating whether the robust semivariogram
#' (Cressie and Hawkins, 1980) is used. The default is \code{FALSE}.
#' @param bins The number of equally spaced bins. The default is 15. Ignored if
#'   \code{cloud = TRUE}.
#' @param cutoff The maximum distance considered.
#'   The default is half the diagonal of the bounding box from the coordinates.
#' @param dist_matrix A distance matrix to be used instead of providing coordinate names.
#' @param partition_factor An optional formula specifying the partition factor.
#'   If specified, semivariances are only computed for observations sharing the
#'   same level of the partition factor.
#'
#' @details The empirical semivariogram is a tool used to visualize and model
#'   spatial dependence by estimating the semivariance of a process at varying distances.
#'   For a constant-mean process, the
#'   semivariance at distance \eqn{h} is denoted \eqn{\gamma(h)} and defined as
#'   \eqn{0.5 * Var(z1  - z2)}. Under second-order stationarity,
#'   \eqn{\gamma(h) = Cov(0) - Cov(h)}, where \eqn{Cov(h)} is the covariance function at distance \code{h}. Typically the residuals from an ordinary
#'   least squares fit defined by \code{formula} are second-order stationary with
#'   mean zero. These residuals are used to compute the empirical semivariogram.
#'   At a distance \code{h}, the empirical semivariance is
#'   \eqn{1/N(h) \sum (r1 - r2)^2}, where \eqn{N(h)} is the number of (unique)
#'   pairs in the set of observations whose distance separation is \code{h} and
#'   \code{r1} and \code{r2} are residuals corresponding to observations whose
#'   distance separation is \code{h}. The robust version is described by
#'   Cressie and Hawkins (1980). In spmodel, these distance bins actually
#'   contain observations whose distance separation is \code{h +- c},
#'   where \code{c} is a constant determined implicitly by \code{bins}. Typically,
#'   only observations whose distance separation is below some cutoff are used
#'   to compute the empirical semivariogram (this cutoff is determined by \code{cutoff}).
#'
#'   When using [splm()] with \code{estmethod} as \code{"sv-wls"}, the empirical
#'   semivariogram is calculated internally and used to estimate spatial
#'   covariance parameters.
#'
#' @name esv
#'
#' @return If \code{cloud = FALSE}, a tibble (data.frame) with distance bins
#'   (\code{bins}), the average distance (\code{dist}), the average semivariance (\code{gamma}), and the
#'   number of (unique) pairs (\code{np}). If \code{cloud = TRUE}, a tibble
#'   (data.frame) with distance (\code{dist}) and semivariance (\code{gamma})
#'   for each unique pair.
#'
#' @export
#'
#' @examples
#' esv(sulfate ~ 1, sulfate)
#' plot(esv(sulfate ~ 1, sulfate))
#' @references Cressie, N & Hawkins, D.M. 1980. Robust estimation of the variogram.
#' \emph{Journal of the International Association for Mathematical Geology},
#' \strong{12}, 115-125.
esv <- function(formula, data, xcoord, ycoord, cloud = FALSE, robust = FALSE, bins = 15, cutoff, dist_matrix, partition_factor) {
  # covert sp to sf
  attr_sp <- attr(class(data), "package")
  if (!is.null(attr_sp) && length(attr_sp) == 1 && attr_sp == "sp") {
    stop("sf objects must be used instead of sp objects. To convert your sp object into an sf object, run sf::st_as_sf().", call. = FALSE)
  }

  # non standard evaluation for the x and y coordinates -- missing() is
  # checked here, on the original (not yet substituted) argument
  # as.character(substitute()) normalizes both quoted ("x") and
  # unquoted (x) column-name references into a plain string. Downstream
  # code checks is.null(xcoord) (not missing()) to see whether the
  # argument was supplied.
  xcoord <- if (missing(xcoord)) NULL else as.character(substitute(xcoord))
  ycoord <- if (missing(ycoord)) NULL else as.character(substitute(ycoord))

  ## convert sf to data frame (point geometry) (1d objects obsolete)
  ### see if data has sf class
  if (inherits(data, "sf")) {
    # collapse any polygon/line geometries to their centroid so distances can
    # be computed between simple x/y points
    data <- suppressWarnings(sf::st_centroid(data))
    data <- sf_to_df(data)
    ### name xcoord ".xcoord" to be used later
    xcoord <- ".xcoord"
    ### name ycoord ".ycoord" to be used later
    ycoord <- ".ycoord"
  }

  # a "." in formula must be expanded before any model.frame()/model.matrix()/
  # lm() call below, excluding the coordinate columns -- see expand_formula_dot()
  formula <- expand_formula_dot(formula, data, c(xcoord, ycoord))

  # filter out missing response values
  na_index <- is.na(data[[all.vars(formula)[1]]])
  data <- data[!na_index, , drop = FALSE]
  # finding model frame
  data_model_frame <- model.frame(formula, data, drop.unused.levels = TRUE, na.action = na.pass)
  # model matrix with potential NA
  # na.action = na.pass above keeps NA rows so we can detect and reject them
  # explicitly here (with a clearer error) rather than silently dropping them
  ob_predictors <- complete.cases(model.matrix(formula, data_model_frame))
  if (any(!ob_predictors)) {
    stop("Cannot have NA values in predictors.", call. = FALSE)
  }

  # compute spatial distances
  if (missing(dist_matrix)) {
    if (is.null(xcoord)) {
      stop("The xcoord argument must be specified.", call. = FALSE)
    }

    if (!is.null(xcoord)) {
      if (!xcoord %in% colnames(data)) {
        stop("The xcoord argument must match the name of a variable in data.", call. = FALSE)
      }
    }

    if (is.null(ycoord)) {
      dist_matrix <- spdist(data, xcoord)
    } else {
      if (!ycoord %in% colnames(data)) {
        stop("The ycoord argument must match the name of a variable in data.", call. = FALSE)
      }
      dist_matrix <- spdist(data, xcoord, ycoord)
    }
  }


  dist_matrix <- as.matrix(dist_matrix)
  # distance matrix is symmetric with a zero diagonal, so only the upper
  # triangle (unique pairs) is needed going forward
  dist_matrix <- dist_matrix[upper.tri(dist_matrix)]

  if (any(dist_matrix == 0)) {
    warning("Zero distances observed between at least one pair. Ignoring pairs. If using splm(), consider a different estimation method.", call. = FALSE)
  }

  if (missing(partition_factor)) partition_factor <- NULL

  if (!is.null(partition_factor)) {
    # partition_matrix_val is a 0/1 indicator of whether a pair shares the same
    # partition level; multiplying zeroes out distances for pairs in different
    # levels so they get filtered out below alongside genuine zero distances
    partition_matrix_val <- as.matrix(partition_matrix(partition_factor, data = data))
    partition_matrix_val <- partition_matrix_val[upper.tri(partition_matrix_val)]
    dist_matrix <- dist_matrix * partition_matrix_val
  }

  if (missing(cutoff)) cutoff <- NULL
  if (is.null(cutoff)) {
    # default cutoff of half the maximum observed distance keeps enough pairs
    # per bin while avoiding the sparsely populated tail of far-apart pairs
    cutoff <- max(dist_matrix) / 2
  }

  dist_vector <- dist_matrix
  if (any(dist_vector == 0)) {
    # zero distances are either coincident locations or cross-partition pairs
    # zeroed out above; either way they should be excluded, not treated as a
    # real distance-zero bin
    dist_index <- dist_vector > 0 & dist_vector <= cutoff
  } else {
    dist_index <- dist_vector <= cutoff
  }

  dist_vector <- dist_vector[dist_index]

  # compute squared differences in the residuals
  lmod <- lm(formula = formula, data = data)
  residuals <- residuals(lmod)
  # reuse spdist() (normally for spatial coordinates) on the 1-d residual
  # vector to get pairwise |r_i - r_j| differences in matrix form for free
  residual_matrix <- as.matrix(spdist(xcoord_val = residuals))
  residual_matrix <- residual_matrix[upper.tri(residual_matrix)]
  if (!is.null(partition_factor)) {
    residual_matrix <- residual_matrix * partition_matrix_val
  }

  residual_vector <- residual_matrix
  residual_vector <- residual_vector[dist_index]
  # semivariance uses squared residual differences, (r_i - r_j)^2
  residual_vector2 <- residual_vector^2

  # cloud returns every pairwise value; otherwise pairs are binned by distance
  # and averaged within each bin, optionally with the robust estimator
  if (cloud) {
    esv_out <- get_esv_cloud(residual_vector2, dist_vector)
  } else {
    if (robust) {
      # Cressie-Hawkins robust estimator works on sqrt(|difference|) rather
      # than the squared difference, which downweights the influence of outliers
      residual_vector12 <- sqrt(residual_vector)
      esv_out <- get_esv_robust(residual_vector12, dist_vector, bins, cutoff)
    } else {
      esv_out <- get_esv(residual_vector2, dist_vector, bins, cutoff)
    }
  }


  esv_out <- structure(esv_out, class = c("esv", class(esv_out)), call = match.call(), cloud = cloud)
  esv_out
}

#' @rdname esv
#' @method plot esv
#' @param x An object from \code{esv()}.
#' @param ... Other arguments passed to other methods.
#' @export
plot.esv <- function(x, ...) {
  cal <- attr(x, "call")
  if (!is.na(m.f <- match("formula", names(cal)))) {
    cal <- cal[c(1, m.f)]
    names(cal)[2L] <- ""
  }
  cc <- deparse(cal, 80)
  nc <- nchar(cc[1L], "c")
  abbr <- length(cc) > 1 || nc > 75
  sub.caption <- if (abbr) {
    paste(substr(cc[1L], 1L, min(75L, nc)), "...")
  } else {
    cc[1L]
  }

  dotlist <- list(...)
  dotlist <- get_esv_dotlist_defaults(x, dotlist, cloud = attr(x, "cloud"))

  do.call("plot", c(list(x = x$dist, y = x$gamma), dotlist))
  title(sub = sub.caption)
}