Merge branch 'main' of https://github.com/USEPA/spmodel

Author: michaeldumelle

DESCRIPTION

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

NEWS.md

@@ -1,3 +1,9 @@
+# spmodel 0.9.0
+
+## Major Updates
+
+* Added the `range_constrain` argument to `splm()` and `spglm()` to contrain the range parameter to enhance numerical stability.
+
 # spmodel 0.8.0
 
 ## Major Updates

R/cov_estimate_sv.R

@@ -43,7 +43,8 @@ cov_estimate_sv <- function(data_object, formula, spcov_initial, estmethod,
   if (all(spcov_initial_val$is_known)) {
     spcov_estimate_val <- use_svloss_known(spcov_initial_val, dist_matrix_list, cov_initial_val$esv, weights)
   } else {
-    spcov_estimate_val <- use_svloss(spcov_initial_val, dist_matrix_list, cov_initial_val$esv, weights, optim_dotlist)
+    spcov_estimate_val <- use_svloss(spcov_initial_val, dist_matrix_list, cov_initial_val$esv, weights, optim_dotlist,
+                                     data_object = data_object)
   }
   spcov_estimate_val
 }

R/cov_initial_search.R

@@ -146,6 +146,7 @@ cov_initial_search.exponential <- function(spcov_initial_NA, estmethod, data_obj
       dist_vector = dist_vector,
       residual_vector2 = residual_vector2
     )
+
     # find minimum value and record parameters
     min_params <- unlist(cov_grid_splits[[which.min(objvals)]])
     spcov_params <- min_params[c("de", "ie", "range", "rotate", "scale")]

R/get_data_object.R

@@ -1,5 +1,6 @@
 get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, estmethod,
-                                 anisotropy, random, randcov_initial, partition_factor, local, ...) {
+                                 anisotropy, random, randcov_initial, partition_factor, local,
+                                 range_constrain, ...) {
 
 
   # covert sp to sf
@@ -198,6 +199,40 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e
   y_range <- range(obdata[[ycoord]])
   max_halfdist <- sqrt((max(x_range) - min(x_range))^2 + (max(y_range) - min(y_range))^2) / 2
 
+  # if (is.logical(range_constrain)) {
+  #   if (range_constrain) {
+  #     range_constrain_value <- 2 * max_halfdist * 5
+  #   } else {
+  #     range_constrain_value <- Inf
+  #   }
+  # } else if (is.numeric(range_constrain)) {
+  #   range_constrain <- TRUE
+  #   range_constrain_value <- 0
+  # } else {
+  #   stop("range_constrain must be logical or numeric.", call. = FALSE)
+  # }
+
+  # range constrain
+  max_range_scale <- 5
+  range_constrain_value <- 2 * max_halfdist * max_range_scale
+  if ("range" %in% names(spcov_initial$is_known)) {
+    if (spcov_initial$is_known[["range"]] || (spcov_initial$initial[["range"]] > range_constrain_value)) {
+      range_constrain <- FALSE
+    }
+  }
+
+  if (inherits(spcov_initial, "none")) {
+    range_constrain <- FALSE
+  }
+
+  if (is.logical(range_constrain)) {
+    if (!range_constrain) {
+      range_constrain_value <- NULL
+    }
+  } else {
+    stop("range_constrain must be logical.", call. = FALSE)
+  }
+
   # override anisotropy argument if needed
   anisotropy <- get_anisotropy_corrected(anisotropy, spcov_initial)
 
@@ -295,7 +330,8 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e
     randcov_list = randcov_list, randcov_names = randcov_names,
     sf_column_name = sf_column_name, terms = terms_val, var_adjust = local$var_adjust,
     X_list = X_list, xcoord = xcoord, xlevels = xlevels, y_list = y_list, ycoord = ycoord,
-    ycoord_orig_name = ycoord_orig_name, ycoord_orig_val = ycoord_orig_val, s2 = s2, diagtol = diagtol
+    ycoord_orig_name = ycoord_orig_name, ycoord_orig_val = ycoord_orig_val, s2 = s2, diagtol = diagtol,
+    range_constrain = range_constrain, range_constrain_value = range_constrain_value
   )
 }
 

R/get_data_object_glm.R

@@ -1,5 +1,6 @@
 get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord, ycoord, estmethod,
-                                  anisotropy, random, randcov_initial, partition_factor, local, ...) {
+                                  anisotropy, random, randcov_initial, partition_factor, local,
+                                  range_constrain, ...) {
 
 
 
@@ -231,6 +232,27 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord,
   y_range <- range(obdata[[ycoord]])
   max_halfdist <- sqrt((max(x_range) - min(x_range))^2 + (max(y_range) - min(y_range))^2) / 2
 
+  # range constrain
+  max_range_scale <- 5
+  range_constrain_value <- 2 * max_halfdist * max_range_scale
+  if ("range" %in% names(spcov_initial$is_known)) {
+    if (spcov_initial$is_known[["range"]] || (spcov_initial$initial[["range"]] > range_constrain_value)) {
+      range_constrain <- FALSE
+    }
+  }
+
+  if (inherits(spcov_initial, "none")) {
+    range_constrain <- FALSE
+  }
+
+  if (is.logical(range_constrain)) {
+    if (!range_constrain) {
+      range_constrain_value <- NULL
+    }
+  } else {
+    stop("range_constrain must be logical.", call. = FALSE)
+  }
+
   # override anisotropy argument if needed
   anisotropy <- get_anisotropy_corrected(anisotropy, spcov_initial)
 
@@ -334,7 +356,8 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord,
     randcov_list = randcov_list, randcov_names = randcov_names,
     sf_column_name = sf_column_name, size = size, terms = terms_val, var_adjust = local$var_adjust,
     X_list = X_list, xcoord = xcoord, xlevels = xlevels, y_list = y_list, ycoord = ycoord,
-    ycoord_orig_name = ycoord_orig_name, ycoord_orig_val = ycoord_orig_val, s2 = s2, diagtol = diagtol
+    ycoord_orig_name = ycoord_orig_name, ycoord_orig_val = ycoord_orig_val, s2 = s2, diagtol = diagtol,
+    range_constrain = range_constrain, range_constrain_value = range_constrain_value
   )
 }
 

R/glogclik.R

@@ -8,9 +8,10 @@
 #' @return (Gaussian log) composite likelihood while optimzing
 #'
 #' @noRd
-glogclik <- function(par, spcov_orig2optim, residual_vector2, dist_vector) {
+glogclik <- function(par, spcov_orig2optim, residual_vector2, dist_vector, data_object) {
   # transforming to original scale
-  spcov_orig_val <- spcov_optim2orig(spcov_orig2optim, par, spcov_profiled = FALSE)
+  spcov_orig_val <- spcov_optim2orig(spcov_orig2optim, par, spcov_profiled = FALSE,
+                                     data_object = data_object)
   # making a covariance parameter vector
   spcov_params_val <- get_spcov_params(spcov_type = class(spcov_orig2optim), spcov_orig_val = spcov_orig_val)
   glogclikloss_val <- get_glogclikloss(spcov_params_val, residual_vector2, dist_vector)

R/gloglik_anis.R

@@ -23,7 +23,8 @@ gloglik_anis <- function(par, spcov_orig2optim, data_object, estmethod,
                          randcov_profiled = NULL) {
 
   # transforming to original scale
-  spcov_orig_val <- spcov_optim2orig(spcov_orig2optim, par, spcov_profiled = spcov_profiled)
+  spcov_orig_val <- spcov_optim2orig(spcov_orig2optim, par, spcov_profiled = spcov_profiled,
+                                     data_object = data_object)
 
   # making a covariance parameter vector
   spcov_params_val <- get_spcov_params(spcov_type = class(spcov_orig2optim), spcov_orig_val = spcov_orig_val)

R/laploglik_anis.R

@@ -9,7 +9,8 @@ laploglik_anis <- function(par, spcov_orig2optim, dispersion_orig2optim, data_ob
   par <- dispersion_orig_val$new_par
 
   # transforming to original scale
-  spcov_orig_val <- spcov_optim2orig(spcov_orig2optim, par, spcov_profiled = spcov_profiled)
+  spcov_orig_val <- spcov_optim2orig(spcov_orig2optim, par, spcov_profiled = spcov_profiled,
+                                     data_object = data_object)
 
   # making a covariance parameter vector
   spcov_params_val <- get_spcov_params(spcov_type = class(spcov_orig2optim), spcov_orig_val = spcov_orig_val)

R/spcov_optim2orig.R

@@ -24,7 +24,12 @@ spcov_optim2orig.exponential <- function(spcov_orig2optim, par, spcov_profiled,
     ie <- exp(fill_optim_par_val[["ie_log"]])
   }
 
-  range <- exp(fill_optim_par_val[["range_log"]])
+  if (data_object$range_constrain) {
+    range <- expit(fill_optim_par_val[["range_logodds"]]) * data_object$range_constrain_value
+  } else {
+    range <- exp(fill_optim_par_val[["range_log"]])
+  }
+
   rotate <- pi * expit(fill_optim_par_val[["rotate_logodds"]])
   scale <- expit(fill_optim_par_val[["scale_logodds"]])
 
@@ -76,7 +81,12 @@ spcov_optim2orig.matern <- function(spcov_orig2optim, par, spcov_profiled, data_
   extra_t <- expit(fill_optim_par_val[["extra_logodds"]])
   # fix to be in [1/5, 5]
   extra <- extra_t * (5 - 1 / 5) + 1 / 5
-  range <- exp(fill_optim_par_val[["range_log"]])
+  # range <- exp(fill_optim_par_val[["range_log"]])
+  if (data_object$range_constrain) {
+    range <- expit(fill_optim_par_val[["range_logodds"]]) * data_object$range_constrain_value
+  } else {
+    range <- exp(fill_optim_par_val[["range_log"]])
+  }
   rotate <- pi * expit(fill_optim_par_val[["rotate_logodds"]])
   scale <- expit(fill_optim_par_val[["scale_logodds"]])
 
@@ -99,7 +109,12 @@ spcov_optim2orig.cauchy <- function(spcov_orig2optim, par, spcov_profiled, data_
 
 
   extra <- exp(fill_optim_par_val[["extra_log"]])
-  range <- exp(fill_optim_par_val[["range_log"]])
+  # range <- exp(fill_optim_par_val[["range_log"]])
+  if (data_object$range_constrain) {
+    range <- expit(fill_optim_par_val[["range_logodds"]]) * data_object$range_constrain_value
+  } else {
+    range <- exp(fill_optim_par_val[["range_log"]])
+  }
   rotate <- pi * expit(fill_optim_par_val[["rotate_logodds"]])
   scale <- expit(fill_optim_par_val[["scale_logodds"]])
 
@@ -122,7 +137,12 @@ spcov_optim2orig.pexponential <- function(spcov_orig2optim, par, spcov_profiled,
 
   extra_half <- expit(fill_optim_par_val[["extra_logodds"]])
   extra <- 2 * extra_half
-  range <- exp(fill_optim_par_val[["range_log"]])
+  # range <- exp(fill_optim_par_val[["range_log"]])
+  if (data_object$range_constrain) {
+    range <- expit(fill_optim_par_val[["range_logodds"]]) * data_object$range_constrain_value
+  } else {
+    range <- exp(fill_optim_par_val[["range_log"]])
+  }
   rotate <- pi * expit(fill_optim_par_val[["rotate_logodds"]])
   scale <- expit(fill_optim_par_val[["scale_logodds"]])
 

R/spcov_orig2optim.R

@@ -41,9 +41,16 @@ spcov_orig2optim.exponential <- function(spcov_initial, spcov_profiled, data_obj
 
   # range changes based on type
   range <- spcov_initial$initial[["range"]]
-  range_log <- log(range)
-  spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
-  spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  if (data_object$range_constrain) {
+    range_prop <- range / data_object$range_constrain_value
+    range_logodds <- logit(range_prop)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_logodds = range_logodds)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_logodds = spcov_initial$is_known[["range"]])
+  } else {
+    range_log <- log(range)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  }
 
   # anisotropy parameters
   ## rotate (between 0 and pi radians)
@@ -124,9 +131,20 @@ spcov_orig2optim.matern <- function(spcov_initial, spcov_profiled, data_object,
 
   # range changes based on type
   range <- spcov_initial$initial[["range"]]
-  range_log <- log(range)
-  spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
-  spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  if (data_object$range_constrain) {
+    range_prop <- range / data_object$range_constrain_value
+    range_logodds <- logit(range_prop)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_logodds = range_logodds)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_logodds = spcov_initial$is_known[["range"]])
+  } else {
+    range_log <- log(range)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  }
+  # range <- spcov_initial$initial[["range"]]
+  # range_log <- log(range)
+  # spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+  # spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
 
   # # extra p log (for now)
   # extra <- spcov_initial$initial[["extra"]]
@@ -191,9 +209,20 @@ spcov_orig2optim.cauchy <- function(spcov_initial, spcov_profiled, data_object,
 
   # range changes based on type
   range <- spcov_initial$initial[["range"]]
-  range_log <- log(range)
-  spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
-  spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  if (data_object$range_constrain) {
+    range_prop <- range / data_object$range_constrain_value
+    range_logodds <- logit(range_prop)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_logodds = range_logodds)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_logodds = spcov_initial$is_known[["range"]])
+  } else {
+    range_log <- log(range)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  }
+  # range <- spcov_initial$initial[["range"]]
+  # range_log <- log(range)
+  # spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+  # spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
 
   # extra p log
   extra <- spcov_initial$initial[["extra"]]
@@ -252,9 +281,20 @@ spcov_orig2optim.pexponential <- function(spcov_initial, spcov_profiled, data_ob
 
   # range changes based on type
   range <- spcov_initial$initial[["range"]]
-  range_log <- log(range)
-  spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
-  spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  if (data_object$range_constrain) {
+    range_prop <- range / data_object$range_constrain_value
+    range_logodds <- logit(range_prop)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_logodds = range_logodds)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_logodds = spcov_initial$is_known[["range"]])
+  } else {
+    range_log <- log(range)
+    spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+    spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
+  }
+  # range <- spcov_initial$initial[["range"]]
+  # range_log <- log(range)
+  # spcov_orig2optim_val <- c(spcov_orig2optim_val, range_log = range_log)
+  # spcov_orig2optim_is_known <- c(spcov_orig2optim_is_known, range_log = spcov_initial$is_known[["range"]])
 
   # extra p logodds (for now)
   extra <- spcov_initial$initial[["extra"]]

R/spglm.R

@@ -121,6 +121,14 @@
 #'   If \code{local} is \code{TRUE}, defaults for \code{local} are chosen such
 #'   that \code{local} is transformed into
 #'   \code{list(size = 100, method = "kmeans", var_adjust = "theoretical", parallel = FALSE)}.
+#' @param range_constrain An optional logical that indicates whether the range
+#'   should be constrained to enhance numerical stability. If \code{range_constrain = TRUE},
+#'   the maximum possible range value is 5 times the maximum distance in the domain.
+#'   If \code{range_constrain = FALSE}, then maximum possible range is unbounded.
+#'   The default is \code{FALSE}.
+#'   Note that if \code{range_constrain = TRUE} and the value of \code{range} in \code{spcov_initial}
+#'   is larger than \code{range_constrain}, then \code{range_constrain} is set to
+#'   \code{FALSE}.
 #' @param ... Other arguments to [esv()] or [stats::optim()].
 #'
 #' @details The spatial generalized linear model for point-referenced data
@@ -299,7 +307,7 @@
 #' McCullagh P. and Nelder, J. A. (1989) \emph{Generalized Linear Models}. London: Chapman and Hall.
 spglm <- function(formula, family, data, spcov_type, xcoord, ycoord, spcov_initial,
                   dispersion_initial, estmethod = "reml", anisotropy = FALSE,
-                  random, randcov_initial, partition_factor, local, ...) {
+                  random, randcov_initial, partition_factor, local, range_constrain, ...) {
 
   # set exponential as default if nothing specified
   if (missing(spcov_type) && missing(spcov_initial)) {
@@ -387,6 +395,10 @@ spglm <- function(formula, family, data, spcov_type, xcoord, ycoord, spcov_initi
     local <- NULL
   }
 
+  if (missing(range_constrain)) {
+    range_constrain <- TRUE
+  }
+
   # non standard evaluation for x and y coordinates
   xcoord <- substitute(xcoord)
   ycoord <- substitute(ycoord)
@@ -395,7 +407,7 @@ spglm <- function(formula, family, data, spcov_type, xcoord, ycoord, spcov_initi
   data_object <- get_data_object_spglm(
     formula, family, data, spcov_initial, xcoord, ycoord,
     estmethod, anisotropy, random, randcov_initial,
-    partition_factor, local, ...
+    partition_factor, local, range_constrain, ...
   )
 
   # parallel cluster if necessary