USEPA · michaeldumelle · Nov 7, 2024 · Sep 19, 2024 · Sep 20, 2024 · Sep 20, 2024
diff --git a/NAMESPACE b/NAMESPACE
@@ -38,6 +38,7 @@ S3method(cov_initial_search,cubic)
 S3method(cov_initial_search,exponential)
 S3method(cov_initial_search,gaussian)
 S3method(cov_initial_search,gravity)
+S3method(cov_initial_search,ie)
 S3method(cov_initial_search,jbessel)
 S3method(cov_initial_search,magnetic)
 S3method(cov_initial_search,matern)
@@ -57,6 +58,7 @@ S3method(cov_initial_search_glm,cubic)
 S3method(cov_initial_search_glm,exponential)
 S3method(cov_initial_search_glm,gaussian)
 S3method(cov_initial_search_glm,gravity)
+S3method(cov_initial_search_glm,ie)
 S3method(cov_initial_search_glm,jbessel)
 S3method(cov_initial_search_glm,magnetic)
 S3method(cov_initial_search_glm,matern)
@@ -95,6 +97,7 @@ S3method(get_initial_range,cubic)
 S3method(get_initial_range,exponential)
 S3method(get_initial_range,gaussian)
 S3method(get_initial_range,gravity)
+S3method(get_initial_range,ie)
 S3method(get_initial_range,jbessel)
 S3method(get_initial_range,magnetic)
 S3method(get_initial_range,matern)
@@ -125,6 +128,7 @@ S3method(gloglik_products,cubic)
 S3method(gloglik_products,exponential)
 S3method(gloglik_products,gaussian)
 S3method(gloglik_products,gravity)
+S3method(gloglik_products,ie)
 S3method(gloglik_products,jbessel)
 S3method(gloglik_products,magnetic)
 S3method(gloglik_products,matern)
@@ -156,6 +160,7 @@ S3method(laploglik_products,cubic)
 S3method(laploglik_products,exponential)
 S3method(laploglik_products,gaussian)
 S3method(laploglik_products,gravity)
+S3method(laploglik_products,ie)
 S3method(laploglik_products,jbessel)
 S3method(laploglik_products,magnetic)
 S3method(laploglik_products,matern)
@@ -236,6 +241,7 @@ S3method(spcov_matrix,cubic)
 S3method(spcov_matrix,exponential)
 S3method(spcov_matrix,gaussian)
 S3method(spcov_matrix,gravity)
+S3method(spcov_matrix,ie)
 S3method(spcov_matrix,jbessel)
 S3method(spcov_matrix,magnetic)
 S3method(spcov_matrix,matern)
@@ -257,6 +263,7 @@ S3method(spcov_optim2orig,cubic)
 S3method(spcov_optim2orig,exponential)
 S3method(spcov_optim2orig,gaussian)
 S3method(spcov_optim2orig,gravity)
+S3method(spcov_optim2orig,ie)
 S3method(spcov_optim2orig,jbessel)
 S3method(spcov_optim2orig,magnetic)
 S3method(spcov_optim2orig,matern)
@@ -276,6 +283,7 @@ S3method(spcov_orig2optim,cubic)
 S3method(spcov_orig2optim,exponential)
 S3method(spcov_orig2optim,gaussian)
 S3method(spcov_orig2optim,gravity)
+S3method(spcov_orig2optim,ie)
 S3method(spcov_orig2optim,jbessel)
 S3method(spcov_orig2optim,magnetic)
 S3method(spcov_orig2optim,matern)
@@ -294,6 +302,7 @@ S3method(spcov_vector,cubic)
 S3method(spcov_vector,exponential)
 S3method(spcov_vector,gaussian)
 S3method(spcov_vector,gravity)
+S3method(spcov_vector,ie)
 S3method(spcov_vector,jbessel)
 S3method(spcov_vector,magnetic)
 S3method(spcov_vector,matern)
@@ -312,6 +321,7 @@ S3method(sprnorm,cubic)
 S3method(sprnorm,exponential)
 S3method(sprnorm,gaussian)
 S3method(sprnorm,gravity)
+S3method(sprnorm,ie)
 S3method(sprnorm,jbessel)
 S3method(sprnorm,magnetic)
 S3method(sprnorm,matern)
@@ -423,6 +433,7 @@ importFrom(stats,model.frame)
 importFrom(stats,model.matrix)
 importFrom(stats,model.offset)
 importFrom(stats,model.response)
+importFrom(stats,na.fail)
 importFrom(stats,na.omit)
 importFrom(stats,na.pass)
 importFrom(stats,pchisq)

diff --git a/NEWS.md b/NEWS.md
@@ -2,7 +2,20 @@
 
 ## Major Updates
 
-* Added the `range_constrain` argument to `splm()` and `spglm()` to contrain the range parameter to enhance numerical stability.
+* Added the `range_constrain` argument to `splm()` and `spglm()` to constrain the range parameter to enhance numerical stability. The default for `range_constrain` is `FALSE`, implying the range is not constrained.
+* Updated the `seal` data with additional polygons and a factor variable, `stock`, with two levels (`8` and `10`) that indicates seal stock (i.e., seal type).
+
+## Minor Updates
+
+* Changed diagonal tolerance threshold for `spglm()` and `spgautor()` model objects. See [this link](https://usepa.github.io/spmodel/articles/technical.html#sec:computational) for details.
+* Added the `"ie"` spatial covariance type to `splm()` and `spglm()` models. For `splm()` models, `"ie"` is an alias for `"none"`. For `spglm()` models, `"none"` now fixes both the `de` and `ie` covariance parameters at zero, while `"ie"` fixes the `de` covariance parameter at zero but allows the `ie` covariance parameter to vary. Thus, `"none"` from `spmodel $\le$ v0.8.0` matches `"ie"` from `spmodel` v0.9.0 and but is different from `"none"` from `spmodel v0.9.0`.
+* Added the `na.action` argument to `predict.spmodel()` functions to clarify that missing values in `newdata` return an error.
+* Minor documentation updates.
+
+## Bug Fixes
+
+* Fixed a bug that caused incorrect degrees of freedom for the likelihood ratio test (`anova(model1, model2)`) when `estmethod` is `"ml"` for both models.
+* Fixed a bug that caused an error in `anova(object1, object2)` when the name of `object1` had special characters (e.g., `$`).
 
 # spmodel 0.8.0
 
@@ -149,7 +162,7 @@
 * Fixed a bug in `spautor()` that prevented an error from occurring when a partition factor was not categorical or not a factor
 * Fixed a bug in `covmatrix(object, newdata)` that returned a matrix with improper dimensions when `spcov_type` was `"none"`.
 * Fixed a bug in `predict()` that caused an error when at least one level of a fixed effect factor was not observed within a local neighborhood (when the `local` method was `"covariance"` or `"distance")`.
-* Fixed a bug in `cooks.distance()` that used the Pearson residuals instead of the standarized residuals.
+* Fixed a bug in `cooks.distance()` that used the Pearson residuals instead of the standardized residuals.
 
 # spmodel 0.3.0
 

diff --git a/R/anova.R b/R/anova.R
@@ -81,7 +81,6 @@
 #' tidy(anova(spmod, lmod))
 anova.splm <- function(object, ..., test = TRUE, Terms, L) {
 
-
   # see if one or two models
   object2_list <- list(...)
 
@@ -139,13 +138,19 @@ anova.splm <- function(object, ..., test = TRUE, Terms, L) {
       stop("The fixed effect coefficients must be the same when performing a likeihood ratio test using the reml estimation method. To perform the likelihood ratio tests for different fixed effect and covariance coefficients simultaneously, refit the models using the ml estimation method.", call. = FALSE)
     }
     Chi2_stat <- abs(-2 * (logLik(object2) - logLik(object)))
-    df_diff <- abs(object2$npar - object$npar)
+
+    # df for ml vs reml
+    df1 <- object$npar
+    df2 <- object2$npar
+    if (object$estmethod == "ml") df1 <- df1 + object$p
+    if (object2$estmethod == "ml") df2 <- df2 + object2$p
+    df_diff <- abs(df1 - df2)
     p_value <- pchisq(Chi2_stat, df_diff, lower.tail = FALSE)
     if (object2$npar < object$npar) {
-      full_name <- as.character(substitute(object))
+      full_name <- deparse(substitute(object)) # replace as.character with deparse
       reduced_name <- as.character(as.list(substitute(list(...)))[-1])
     } else {
-      reduced_name <- as.character(substitute(object))
+      reduced_name <- deparse(substitute(object)) # replace as.character with deparse
       full_name <- as.character(as.list(substitute(list(...)))[-1])
     }
     if (test) {

diff --git a/R/cov_betahat_adjust.R b/R/cov_betahat_adjust.R
@@ -18,7 +18,7 @@ cov_betahat_adjust <- function(invcov_betahat_list, betahat_list,
                                randcov_params, cov_betahat_noadjust, var_adjust) {
   P <- length(betahat_list)
   # reset var_adjust if only one partition
-  if (P == 1 || inherits(spcov_params, "none")) {
+  if (P == 1 || inherits(spcov_params, c("none", "ie"))) {
     var_adjust <- "none"
   }
   # reset var_adjust if partitioning used but no local option used

diff --git a/R/cov_estimate_gloglik.R b/R/cov_estimate_gloglik.R
@@ -19,7 +19,7 @@ cov_estimate_gloglik_splm <- function(data_object, formula, spcov_initial, estme
   if (data_object$anisotropy) {
     dist_matrix_list <- NULL
   } else {
-    if (inherits(spcov_initial, "none") && is.null(data_object$randcov_initial)) {
+    if (inherits(spcov_initial, c("none", "ie")) && is.null(data_object$randcov_initial)) {
       dist_matrix_list <- NULL
     } else {
       dist_matrix_list <- lapply(data_object$obdata_list, function(x) spdist(x, data_object$xcoord, data_object$ycoord))

diff --git a/R/cov_initial_search.R b/R/cov_initial_search.R
@@ -585,6 +585,9 @@ cov_initial_search.none <- function(spcov_initial_NA, estmethod, data_object,
   best_params
 }
 
+#' @export
+cov_initial_search.ie <- cov_initial_search.none
+
 #' @export
 cov_initial_search.matern <- function(spcov_initial_NA, estmethod, data_object,
                                       dist_matrix_list, weights,

diff --git a/R/cov_initial_search_glm.R b/R/cov_initial_search_glm.R
@@ -461,6 +461,9 @@ cov_initial_search_glm.none <- function(spcov_initial_NA, dispersion_initial_NA,
   best_params
 }
 
+#' @export
+cov_initial_search_glm.ie <- cov_initial_search_glm.none
+
 #' @export
 cov_initial_search_glm.matern <- function(spcov_initial_NA, dispersion_initial_NA, estmethod, data_object,
                                           dist_matrix_list, weights,

diff --git a/R/data.R b/R/data.R
@@ -62,9 +62,12 @@
 #'
 #' @description Estimated harbor-seal trends from abundance data in southeast Alaska, USA.
 #'
-#' @format A \code{sf} object with 62 rows and 2 columns:
+#' @format A \code{sf} object with 149 rows and 2 columns:
 #' \itemize{
 #'   \item log_trend: The log of the estimated harbor-seal trends from abundance data.
+#'   \item stock: A seal stock factor with two levels: 8 and 10. The factor levels indicate the
+#'     type of seal stock (i.e., type of seal). Stocks 8 and 10 are two distinct stocks
+#'     (out of 13 total stocks) in southeast Alaska.
 #'   \item geometry: \code{POLYGON} geometry representing polygons in an Alaska
 #'     Albers projection (EPSG: 3338).
 #' }
@@ -124,7 +127,7 @@
 #'   \item strat: A factor representing strata (used for sampling). Can take values \code{L} and \code{M}.
 #'   \item count: The count (number) of moose observed.
 #'   \item presence: A binary factor representing whether no moose were observed (value \code{0}) or at least one moose was observed
-#'     (va ue \code{1}).
+#'     (value \code{1}).
 #'   \item geometry: \code{POINT} geometry representing coordinates in an Alaska
 #'     Albers projection (EPSG: 3338). Distances between points are in meters.
 #' }

diff --git a/R/get_data_object.R b/R/get_data_object.R
@@ -31,7 +31,7 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e
     data_sf <- NULL
   }
 
-  if (!is_sf && missing(xcoord) && !inherits(spcov_initial, "none")) {
+  if (!is_sf && missing(xcoord) && !inherits(spcov_initial, c("none", "ie"))) {
     stop("The xcoord argument must be specified.", call. = FALSE)
   }
 
@@ -52,7 +52,7 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e
   ycoord_orig_name <- NULL
   ycoord_orig_val <- NULL
   # find coordinate dimension and set defaults
-  if (inherits(spcov_initial, "none") && estmethod %in% c("reml", "ml")) {
+  if (inherits(spcov_initial, c("none", "ie")) && estmethod %in% c("reml", "ml")) {
     dim_coords <- 0
     if (missing(xcoord)) {
       xcoord <- ".xcoord"
@@ -213,15 +213,15 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e
   # }
 
   # range constrain
-  max_range_scale <- 5
+  max_range_scale <- 4
   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")) {
+  if (inherits(spcov_initial, c("none", "ie"))) {
     range_constrain <- FALSE
   }
 

diff --git a/R/get_data_object_glm.R b/R/get_data_object_glm.R
@@ -32,7 +32,7 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord,
     data_sf <- NULL
   }
 
-  if (!is_sf && missing(xcoord) && !inherits(spcov_initial, "none")) {
+  if (!is_sf && missing(xcoord) && !inherits(spcov_initial, c("none", "ie"))) {
     stop("The xcoord argument must be specified.", call. = FALSE)
   }
 
@@ -53,7 +53,7 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord,
   ycoord_orig_name <- NULL
   ycoord_orig_val <- NULL
   # find coordinate dimension and set defaults
-  if (inherits(spcov_initial, "none") && estmethod %in% c("reml", "ml")) {
+  if (inherits(spcov_initial, c("none", "ie")) && estmethod %in% c("reml", "ml")) {
     dim_coords <- 0
     if (missing(xcoord)) {
       xcoord <- ".xcoord"
@@ -222,8 +222,8 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord,
   betahat <- backsolve(R_val, qr.qty(qr_val, y_trans))
   resid <- y_trans - X %*% betahat
   s2 <- sum(resid^2) / (n - p)
-  # diagtol <- 1e-4
-  diagtol <- min(1e-4, 1e-4 * s2)
+  diagtol <- 1e-4
+  # diagtol <- min(1e-4, 1e-4 * s2)
 
 
 
@@ -233,15 +233,15 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord,
   max_halfdist <- sqrt((max(x_range) - min(x_range))^2 + (max(y_range) - min(y_range))^2) / 2
 
   # range constrain
-  max_range_scale <- 5
+  max_range_scale <- 4
   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")) {
+  if (inherits(spcov_initial, c("none", "ie"))) {
     range_constrain <- FALSE
   }
 

diff --git a/R/get_initial_range.R b/R/get_initial_range.R
@@ -33,6 +33,10 @@ get_initial_range.circular <- function(spcov_type, max_halfdist, ...) {
 get_initial_range.none <- function(spcov_type, max_halfdist, ...) {
   Inf
 }
+
+#' @export
+get_initial_range.ie <- get_initial_range.none
+
 #' @export
 get_initial_range.cubic <- function(spcov_type, max_halfdist, ...) {
   max_halfdist

diff --git a/R/get_spcov_params.R b/R/get_spcov_params.R
@@ -7,7 +7,7 @@
 #'
 #' @noRd
 get_spcov_params <- function(spcov_type, spcov_orig_val) {
-  if (spcov_type %in% c("exponential", "spherical", "gaussian", "triangular", "circular", "none", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic")) {
+  if (spcov_type %in% c("exponential", "spherical", "gaussian", "triangular", "circular", "none", "ie", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic")) {
     spcov_params_val <- spcov_params(
       spcov_type = spcov_type,
       de = spcov_orig_val[["de"]],

diff --git a/R/gloglik_products.R b/R/gloglik_products.R
@@ -85,6 +85,8 @@ gloglik_products.circular <- gloglik_products.exponential
 #' @export
 gloglik_products.none <- gloglik_products.exponential
 #' @export
+gloglik_products.ie <- gloglik_products.none
+#' @export
 gloglik_products.cubic <- gloglik_products.exponential
 #' @export
 gloglik_products.pentaspherical <- gloglik_products.exponential

diff --git a/R/laploglik_products.R b/R/laploglik_products.R
@@ -121,6 +121,8 @@ laploglik_products.circular <- laploglik_products.exponential
 #' @export
 laploglik_products.none <- laploglik_products.exponential
 #' @export
+laploglik_products.ie <- laploglik_products.none
+#' @export
 laploglik_products.cubic <- laploglik_products.exponential
 #' @export
 laploglik_products.pentaspherical <- laploglik_products.exponential

diff --git a/R/loocv.R b/R/loocv.R
@@ -64,7 +64,7 @@ loocv.splm <- function(object, cv_predict = FALSE, se.fit = FALSE, local, ...) {
   }
 
   # iid if relevant otherwise pass
-  if (inherits(coef(object, type = "spcov"), "none") && is.null(object$random)) {
+  if (inherits(coef(object, type = "spcov"), c("none", "ie")) && is.null(object$random)) {
     return(loocv_iid(object, cv_predict, se.fit, local))
   }