Merge pull request #33 from USEPA/develop

CRAN v0.10.0

Author: michaeldumelle

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: spmodel
 Title: Spatial Statistical Modeling and Prediction
-Version: 0.9.0
+Version: 0.10.0
 Authors@R: c(
     person(given = "Michael",
            family = "Dumelle",

NEWS.md

@@ -1,3 +1,12 @@
+# spmodel 0.10.0
+
+## Minor updates
+
+* Added a robust semivariogram option to `esv()`; see the `robust`argument to `esv()` ([#28](https://github.com/USEPA/spmodel/issues/28)).
+* Added `"none"` and `"ie"` spatial covariance types (via `spcov_type` or `spcov_initial`) to `spautor()`, `spgautor()`, and `spautorRF()` ([#27](https://github.com/USEPA/spmodel/issues/27)).
+* Fixed the aspect ratio at one for level curve plots of fitted model objects (anisotropic and isotropic).
+* Changed the title for level curve plots of fitted model objects from "Anisotropic level curve" to "Isotropic level curve" when the model is isotropic.
+
 # spmodel 0.9.0
 
 ## Major Updates

R/esv.R

@@ -15,11 +15,13 @@
 #'   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 dist_matrix A distance matrix to be used instead of providing coordinate names.
+#' @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.
@@ -36,7 +38,8 @@
 #'   \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}. In spmodel, these distance bins actually
+#'   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
@@ -59,7 +62,10 @@
 #' @examples
 #' esv(sulfate ~ 1, sulfate)
 #' plot(esv(sulfate ~ 1, sulfate))
-esv <- function(formula, data, xcoord, ycoord, cloud = FALSE, bins = 15, cutoff, dist_matrix, partition_factor) {
+#' @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) {
 
   # filter out missing response values
   na_index <- is.na(data[[all.vars(formula)[1]]])
@@ -166,9 +172,16 @@ esv <- function(formula, data, xcoord, ycoord, cloud = FALSE, bins = 15, cutoff,
   if (cloud) {
     esv_out <- get_esv_cloud(residual_vector2, dist_vector, formula)
   } else {
-    esv_out <- get_esv(residual_vector2, dist_vector, bins, cutoff)
+    if (robust) {
+      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)
+    }
   }
 
+
+
   # remove NA
   # esv_out <- na.omit(esv_out)
   esv_out <- structure(esv_out, class = c("esv", class(esv_out)), call = match.call(), cloud = cloud)

R/get_esv.R

@@ -29,6 +29,39 @@ get_esv <- function(residual_vector2, dist_vector, bins, cutoff, formula) {
   esv_out
 }
 
+get_esv_robust <- function(residual_vector12, dist_vector, bins, cutoff, formula) {
+
+  # compute semivariogram classes
+  dist_classes <- cut(dist_vector, breaks = seq(0, cutoff, length.out = bins + 1))
+
+  # compute squared differences within each class
+  gamma <- tapply(residual_vector12, dist_classes, function(x) {
+    1 / (0.914 + (0.988 / length(x))) * (mean(x)^4)
+  })
+
+  # compute pairs within each class
+  np <- tapply(residual_vector12, dist_classes, length)
+
+  # set as zero if necessary
+  np <- ifelse(is.na(np), 0, np)
+
+  # compute average distance within each class
+  dist <- tapply(dist_vector, dist_classes, mean)
+
+  # return output
+  esv_out <- tibble::tibble(
+    bins = factor(levels(dist_classes), levels = levels(dist_classes)),
+    dist = as.numeric(dist),
+    gamma = as.numeric(gamma),
+    np = as.numeric(np)
+  )
+
+  # set row names to NULL
+  # row.names(esv_out) <- NULL
+
+  esv_out
+}
+
 get_esv_cloud <- function(residual_vector2, dist_vector, formula) {
 
   esv_out <- tibble::tibble(dist = dist_vector, gamma = residual_vector2)

R/get_esv_dotlist.R

@@ -17,8 +17,13 @@ get_esv_dotlist <- function(..., max_halfdist) {
   if (!("cutoff" %in% names(dotlist))) {
     dotlist$cutoff <- max_halfdist
   }
+
+  if (!("robust" %in% names(dotlist))) {
+    dotlist$robust <- FALSE
+  }
+
   # make dotlist esv
-  dotlist_esv <- list(bins = dotlist$bins, cutoff = dotlist$cutoff)
+  dotlist_esv <- list(robust = dotlist$robust, bins = dotlist$bins, cutoff = dotlist$cutoff)
 }
 
 

R/plot.R

@@ -27,7 +27,7 @@
 #'   For [splm()] and [spglm()] fitted model objects, there are two additional values of \code{which}:
 #'   \itemize{
 #'     \item 7: Fitted spatial covariance function vs distance
-#'     \item 8: Fitted anisotropic level curve of equal correlation
+#'     \item 8: Fitted anisotropic (or isotropic) level curve of equal correlation
 #'   }
 #'
 #' @return No return value. Function called for plotting side effects.
@@ -195,17 +195,23 @@ plot.splm <- function(x, which, ...) {
       x_new <- x_orig
       y_new <- y_orig
     }
+    if (x$anisotropy) {
+      main_new <- "Anisotropic level curve"
+    } else {
+      main_new <- "Isotropic level curve"
+    }
     plot(
       x = x_new,
       y = y_new,
       xlab = "x-distance",
       ylab = "y-distance",
-      main = "Anisotropic level curve", # of equal correlation
+      main =  main_new,
       type = "l",
       xlim = c(-1, 1),
       ylim = c(-1, 1),
       xaxt = "n", # remove axis information
       yaxt = "n", # remove axis information
+      asp = 1, # set equal aspect ratio
       ...
     )
     title(sub = sub.caption)

R/spautor.R

@@ -11,14 +11,15 @@
 #'   the variables in \code{fixed}, \code{random}, and \code{partition_factor},
 #'   as well as potentially geographical information.
 #' @param spcov_type The spatial covariance type. Available options include
-#'   \code{"car"} and \code{"sar"}. Parameterizations of each spatial covariance type are
+#'   \code{"car"}, \code{"sar"}, and \code{"none"}. Parameterizations of each spatial covariance type are
 #'   available in Details. When \code{spcov_type} is specified, relevant spatial
 #'   covariance parameters are assumed unknown, requiring estimation.
 #'   \code{spcov_type} is not required (and is
 #'   ignored) if \code{spcov_initial} is provided.  Multiple values can be
 #'   provided in a character vector. Then \code{spautor()} is called iteratively
 #'   for each element and a list is returned for each model fit.
-#'   The default for \code{spcov_type} is \code{"car"}.
+#'   The default for \code{spcov_type} is \code{"car"}. When \code{spcov_type}
+#'   is \code{"none"}, [splm()] is called.
 #' @param spcov_initial An object from [spcov_initial()] specifying initial and/or
 #'   known values for the spatial covariance parameters.
 #'   Not required if \code{spcov_type} is provided. Multiple [spcov_initial()]
@@ -91,6 +92,7 @@
 #'      symmetry condition matrix \eqn{M}
 #'     \item sar: \eqn{[(I - range * W)(I - range * W)^T]^{-1}},
 #'      weights matrix \eqn{W}, \eqn{^T} indicates matrix transpose
+#'     \item none: \eqn{0}
 #'   }
 #'   If there are observations with no neighbors, they are given a unique variance
 #'   parameter called \code{extra}, which must be non-negative.
@@ -187,6 +189,19 @@ spautor <- function(formula, data, spcov_type, spcov_initial, estmethod = "reml"
     return(new_spautor_out)
   }
 
+
+
+  # call splm if spcov_type is none (works on an individual element of spcov_type or initial)
+  if ((!missing(spcov_type) && spcov_type %in% c("none", "ie")) || (!missing(spcov_initial) && inherits(spcov_initial, c("none", "ie")))) {
+    call_list <- as.list(match.call())[-1]
+    # remove spautor specific arguments
+    args_remove <- c("W", "row_st", "range_positive", "cutoff")
+    call_list <- call_list[!names(call_list) %in% args_remove]
+    penv <- parent.frame()
+    splm_out <- do.call("splm", call_list, envir = penv)
+    return(splm_out)
+  }
+
   # replace initial values with appropriate NA's
   if (missing(spcov_initial)) {
     spcov_initial <- spmodel::spcov_initial(spcov_type)

R/spautorRF.R

@@ -116,10 +116,14 @@ spautorRF <- function(formula, data, ...) {
     data <- data[order(c(which(!na_index), which(na_index))), , drop = FALSE]
     # perform splm
     spautor_out <- do.call(spmodel::spautor, c(list(formula = .ranger_resid ~ 1, data = data), spautor_args), envir = penv)
-    if (inherits(spautor_out, "spautor")) {
+    if (inherits(spautor_out, c("spautor"))) {
       spautor_out$call <- NA
       # output list with names and class
       sprf_out <- structure(list(call = match.call(), ranger = ranger_out, spautor = spautor_out, newdata = newdata), class = "spautorRF")
+    } else if (inherits(spautor_out, c("splm"))) { # splm for none and ie covariance
+      spautor_out$call <- NA
+      # output list with names and class
+      sprf_out <- structure(list(call = match.call(), ranger = ranger_out, spautor = spautor_out, newdata = newdata), class = "splmRF")
     } else {
       spautor_out <- lapply(spautor_out, function(x) {
         x$call <- NA

R/spgautor.R

@@ -19,14 +19,15 @@
 #'   the variables in \code{fixed}, \code{random}, and \code{partition_factor},
 #'   as well as potentially geographical information.
 #' @param spcov_type The spatial covariance type. Available options include
-#'   \code{"car"} and \code{"sar"}. Parameterizations of each spatial covariance type are
+#'   \code{"car"}, \code{"sar"}, \code{"none"}, and \code{"ie"}. Parameterizations of each spatial covariance type are
 #'   available in Details. When \code{spcov_type} is specified, relevant spatial
 #'   covariance parameters are assumed unknown, requiring estimation.
 #'   \code{spcov_type} is not required (and is
 #'   ignored) if \code{spcov_initial} is provided.  Multiple values can be
 #'   provided in a character vector. Then \code{spgautor()} is called iteratively
 #'   for each element and a list is returned for each model fit.
-#'   The default for \code{spcov_type} is \code{"car"}.
+#'   The default for \code{spcov_type} is \code{"car"}. When \code{spcov_type}
+#'   is \code{"none"} or \code{"ie"}, [splm()] is called.
 #' @param spcov_initial An object from [spcov_initial()] specifying initial and/or
 #'   known values for the spatial covariance parameters.
 #'   Not required if \code{spcov_type} is provided. Multiple [spcov_initial()]
@@ -165,6 +166,8 @@
 #'      symmetry condition matrix \eqn{M}
 #'     \item sar: \eqn{[(I - range * W)(I - range * W)^T]^{-1}},
 #'      weights matrix \eqn{W}, \eqn{^T} indicates matrix transpose
+#'     \item none: \eqn{0}
+#'     \item ie: \eqn{0} (see [spglm()] for differences compared to none)
 #'   }
 #'   If there are observations with no neighbors, they are given a unique variance
 #'   parameter called \code{extra}, which must be non-negative.
@@ -267,6 +270,17 @@ spgautor <- function(formula, family, data, spcov_type, spcov_initial, dispersio
     return(new_spgautor_out)
   }
 
+  # call spglm if spcov_type is none (works on an individual element of spcov_type or initial)
+  if ((!missing(spcov_type) && spcov_type %in% c("none", "ie")) || (!missing(spcov_initial) && inherits(spcov_initial, c("none", "ie")))) {
+    call_list <- as.list(match.call())[-1]
+    # remove spautor specific arguments
+    args_remove <- c("W", "row_st", "range_positive", "cutoff")
+    call_list <- call_list[!names(call_list) %in% args_remove]
+    penv <- parent.frame()
+    spglm_out <- do.call("spglm", call_list, envir = penv)
+    return(spglm_out)
+  }
+
   # set dispersion initial
   if (missing(dispersion_initial)) dispersion_initial <- NULL else family <- class(dispersion_initial)