Merge branch 'main' into project

.Rbuildignore

@@ -21,3 +21,4 @@
 ^CRAN-SUBMISSION$
 ^vignettes/web_only$
 revdep/
+^vignettes/articles$

.github/workflows/pkgdown.yaml

@@ -40,6 +40,7 @@ jobs:
           remotes::install_cran("cowplot")
           remotes::install_cran("rnaturalearth")
           remotes::install_cran("rnaturalearthdata")
+          remotes::install_cran("pROC")
           install.packages("Matrix", type = "source")
           install.packages("TMB", type = "source")
         shell: Rscript {0}

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.6.0.9001
+Version: 0.6.0.9010
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "[email protected]",
            role = c("aut", "cre"),
@@ -110,5 +110,5 @@ Config/testthat/parallel: true
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 SystemRequirements: GNU make

NAMESPACE

@@ -14,6 +14,7 @@ S3method(nobs,sdmTMB)
 S3method(plot,sdmTMBmesh)
 S3method(predict,sdmTMB)
 S3method(print,sdmTMB)
+S3method(print,sdmTMB_cv)
 S3method(ranef,sdmTMB)
 S3method(residuals,sdmTMB)
 S3method(simulate,sdmTMB)

NEWS.md

@@ -2,7 +2,31 @@
 
 * Add `project()` experimental function.
 
-* Add progress bar to `simulate.sdmTMB()`.
+* Add print method for `sdmTMB_cv()` output. #319
+
+* Add progress bar to `simulate.sdmTMB()`. #346
+
+* Add AUC and TSS examples to cross validation vignette. #268
+
+* Add `model` (linear predictor number) argument to coef() method. Also,
+  write documentation for `?coef.sdmTMB`. #351
+
+* Add helpful error message if some coordinates in make_mesh() are NA. #365
+
+* Add informative message if fitting with an offset but predicting with offset
+  argument left at NULL on newdata. #372
+
+* Fix passing of `offset` argument through in `sdmTMB_cv()`. Before it was being
+  omitted in the prediction (i.e., set to 0). #372
+
+* Fig bug in `exponentiate` argument for `tidy()`. Set `conf.int = TRUE` as
+  default. #353
+
+* Fix bug in prediction from `delta_truncated_nbinom1()` and 
+  `delta_truncated_nbinom2()` families. The positive component
+  needs to be transformed to represent the mean of the *un*truncated
+  distribution first before multiplying by the probability of a non-zero.
+  Thanks to @tom-peatman #350
 
 * Add `get_eao()` to calculate effective area occupied.
 
@@ -56,7 +80,7 @@
 # sdmTMB 0.5.0
 
 * Overhaul residuals vignette ('article') 
-  <https://pbs-assess.github.io/sdmTMB/articles/web_only/residual-checking.html>
+  <https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html>
   including brief intros to randomized quantile residuals, simulation-based
   residuals, 'one-sample' residuals, and uniform vs. Gaussian residuals.
 

R/cross-val.R

@@ -287,10 +287,10 @@ sdmTMB_cv <- function(
     cli_abort("`weights` cannot be specified within sdmTMB_cv().")
   }
   if ("offset" %in% names(dot_args)) {
-    .offset <- eval(dot_args$offset)
-    if (parallel && !is.character(.offset) && !is.null(.offset)) {
-      cli_abort("We recommend using a character value for 'offset' (indicating the column name) when applying parallel cross validation.")
+    if (!is.character(dot_args$offset)) {
+      cli_abort("Please use a character value for 'offset' (indicating the column name) for cross validation.")
     }
+    .offset <- eval(dot_args$offset)
   } else {
     .offset <- NULL
   }
@@ -369,7 +369,9 @@ sdmTMB_cv <- function(
 
     # FIXME: only use TMB report() below to be faster!
     # predict for withheld data:
-    predicted <- predict(object, newdata = cv_data, type = "response")
+    predicted <- predict(object, newdata = cv_data, type = "response",
+      offset = if (!is.null(.offset)) cv_data[[.offset]] else rep(0, nrow(cv_data)))
+
     cv_data$cv_predicted <- predicted$est
     response <- get_response(object$formula[[1]])
     withheld_y <- predicted[[response]]
@@ -451,7 +453,7 @@ sdmTMB_cv <- function(
   pdHess <- vapply(out, `[[`, "pdHess", FUN.VALUE = logical(1L))
   max_grad <- vapply(out, `[[`, "max_gradient", FUN.VALUE = numeric(1L))
   converged <- all(pdHess)
-  list(
+  out <- list(
     data = data,
     models = models,
     fold_loglik = fold_cv_ll,
@@ -460,9 +462,35 @@ sdmTMB_cv <- function(
     pdHess = pdHess,
     max_gradients = max_grad
   )
+  `class<-`(out, "sdmTMB_cv")
 }
 
 log_sum_exp <- function(x) {
   max_x <- max(x)
   max_x + log(sum(exp(x - max_x)))
 }
+
+#' @export
+#' @import methods
+print.sdmTMB_cv <- function(x, ...) {
+  nmods <- length(x$models)
+  nconverged <- sum(x$converged)
+  cat(paste0("Cross validation of sdmTMB models with ", nmods, " folds.\n"))
+  cat("\n")
+  cat("Summary of the first fold model fit:\n")
+  cat("\n")
+  print(x$models[[1]])
+  cat("\n")
+  cat("Access the rest of the models in a list element named `models`.\n")
+  cat("E.g. `object$models[[2]]` for the 2nd fold model fit.\n")
+  cat("\n")
+  cat(paste0(nconverged, " out of ", nmods, " models are consistent with convergence.\n"))
+  cat("Figure out which folds these are in the `converged` list element.\n")
+  cat("\n")
+  cat(paste0("Out-of-sample log likelihood for each fold: ", paste(round(x$fold_loglik, 2), collapse = ", "), ".\n"))
+  cat("Access these values in the `fold_loglik` list element.\n")
+  cat("\n")
+  cat("Sum of out-of-sample log likelihoods:", round(x$sum_loglik, 2), "\n")
+  cat("More positive values imply better out-of-sample prediction.\n")
+  cat("Access this value in the `sum_loglik` list element.\n")
+}

R/dharma.R

@@ -27,7 +27,7 @@
 #'
 #' @details
 #'
-#' See the [residuals vignette](https://pbs-assess.github.io/sdmTMB/articles/web_only/residual-checking.html).
+#' See the [residuals vignette](https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html).
 #'
 #' Advantages to these residuals over the ones from the [residuals.sdmTMB()]
 #' method are (1) they work with delta/hurdle models for the combined

R/families.R

@@ -222,6 +222,8 @@ nbinom1 <- function(link = "log") {
   add_to_family(x)
 }
 
+utils::globalVariables(".phi") ## avoid R CMD check NOTE
+
 #' @export
 #' @examples
 #' truncated_nbinom2(link = "log")
@@ -235,11 +237,21 @@ truncated_nbinom2 <- function(link = "log") {
     stats <- stats::make.link(linktemp)
   else if (is.character(link))
     stats <- stats::make.link(link)
-
+  linkinv <- function(eta, phi = NULL) {
+    s1 <- eta
+    if (is.null(phi)) phi <- .phi
+    s2 <- logspace_add(0, s1 - log(phi)) # log(1 + mu/phi)
+    log_nzprob <- logspace_sub(0, -phi * s2)
+    exp(eta) / exp(log_nzprob)
+  }
   structure(list(family = "truncated_nbinom2", link = linktemp, linkfun = stats$linkfun,
-    linkinv = stats$linkinv), class = "family")
+    linkinv = linkinv), class = "family")
 }
 
+logspace_sub <- function (lx, ly) lx + log1mexp(lx - ly)
+logspace_add <- function (lx, ly) pmax(lx, ly) + log1p(exp(-abs(lx - ly)))
+log1mexp <- function(x) ifelse(x <= log(2), log(-expm1(-x)), log1p(-exp(-x)))
+
 #' @export
 #' @examples
 #' truncated_nbinom1(link = "log")
@@ -253,9 +265,15 @@ truncated_nbinom1 <- function(link = "log") {
     stats <- stats::make.link(linktemp)
   else if (is.character(link))
     stats <- stats::make.link(link)
-
+  linkinv <- function(eta, phi = NULL) {
+    mu <- exp(eta)
+    if (is.null(phi)) phi <- .phi
+    s2 <- logspace_add(0, log(phi)) # log(1 + phi)
+    log_nzprob <- logspace_sub(0, -mu / phi * s2) # 1 - prob(0)
+    mu / exp(log_nzprob)
+  }
   structure(list(family = "truncated_nbinom1", link = linktemp, linkfun = stats$linkfun,
-    linkinv = stats$linkinv), class = "family")
+    linkinv = linkinv), class = "family")
 }
 
 #' @param df Student-t degrees of freedom fixed value parameter.

R/fit.R

@@ -1501,6 +1501,16 @@ sdmTMB <- function(
   sd_report <- TMB::sdreport(tmb_obj, getJointPrecision = get_joint_precision)
   conv <- get_convergence_diagnostics(sd_report)
 
+  ## save params that families need to grab from environments:
+  if (any(family$family %in% c("truncated_nbinom1", "truncated_nbinom2"))) {
+    phi <- exp(tmb_obj$par[["ln_phi"]])
+    if (delta) {
+      assign(".phi", phi, environment(out_structure[["family"]][[2]][["linkinv"]]))
+    } else {
+      assign(".phi", phi, environment(out_structure[["family"]][["linkinv"]]))
+    }
+  }
+
   out_structure$tmb_obj <- tmb_obj
   out <- c(out_structure, list(
     model      = tmb_opt,

R/mesh.R

@@ -88,6 +88,14 @@ make_mesh <- function(data, xy_cols,
     cli_abort(msg)
   }
 
+  all_x_non_na <- sum(is.na(data[[xy_cols[[1]]]])) == 0L
+  all_y_non_na <- sum(is.na(data[[xy_cols[[2]]]])) == 0L
+  if (!all_x_non_na || !all_y_non_na) {
+    msg <- c("Some coordinates in `xy_cols` were NA.", "
+      Remove or fix these rows before proceeding.")
+    cli_abort(msg)
+  }
+
   if (max(data[[xy_cols[1]]]) > 1e4 || max(data[[xy_cols[2]]] > 1e4)) {
     msg <- paste0(
       "The x or y column values are fairly large. ",
@@ -216,30 +224,11 @@ binary_search_knots <- function(loc_xy,
 #' @param ... Passed to [graphics::plot()].
 #'
 #' @importFrom graphics points
-#' @return `plot.sdmTMBmesh()`: A plot of the mesh and data points. If
-#'   \pkg{ggplot2} is installed, a \pkg{ggplot2} object is
-#'   returned, otherwise a base graphics R plot is returned. To make your own,
-#'   pass `your_mesh$mesh` to `inlabru::gg()`.
+#' @return `plot.sdmTMBmesh()`: A plot of the mesh and data points. To make your
+#'   own \pkg{ggplot2} version, pass `your_mesh$mesh` to `inlabru::gg()`.
 #' @rdname make_mesh
 #' @export
 plot.sdmTMBmesh <- function(x, ...) {
-  # r1 <- requireNamespace("inlabru", quietly = TRUE)
-  # r2 <- requireNamespace("ggplot2", quietly = TRUE)
-  # if (r1 && r2) {
-  #   dat <- data.frame(
-  #     x = x$loc_xy[,1,drop=TRUE],
-  #     y = x$loc_xy[,2,drop=TRUE]
-  #   )
-  #   ggplot2::ggplot() +
-  #     # inlabru::gg(x$mesh, ext.color = "grey20", ext.linewidth = 0.5, edge.color = "grey50") +
-  #     ggplot2::coord_sf() +
-  #     fmesher::geom_fm(data = x$mesh) +
-  #     ggplot2::geom_point(
-  #       data = dat,
-  #       mapping = ggplot2::aes(x = .data$x, y = .data$y), alpha = 0.4, pch = 20, colour = "#3182BD") +
-  #     # ggplot2::coord_fixed() +
-  #     ggplot2::labs(x = x$xy_cols[[1]], y = x$xy_cols[[2]])
-  # } else {
     plot(x$mesh, main = NA, edge.color = "grey60", asp = 1, ...)
     points(x$loc_xy, pch = 21, cex = 0.3, col = "#00000080")
     points(x$loc_centers, pch = 20, col = "red")

R/methods.R

@@ -45,12 +45,20 @@ fitted.sdmTMB <- function(object, ...) {
 #'
 #' @param object The fitted sdmTMB model object
 #' @param complete Currently ignored
+#' @param model Linear predictor for delta models. Defaults to the first
+#'   linear predictor.
 #' @param ... Currently ignored
 #' @importFrom stats coef
 #' @export
-#' @noRd
-coef.sdmTMB <- function(object, complete = FALSE, ...) {
-  x <- tidy(object)
+coef.sdmTMB <- function(object, complete = FALSE, model = 1, ...) {
+  if (is_delta(object)) {
+    assert_that(length(model) == 1L)
+    model <- as.integer(model)
+    assert_that(model %in% c(1L, 2L))
+    msg <- paste0("Returning coefficients from linear predictor ", model, " based on the `model` argument.")
+    cli_inform(msg)
+  }
+  x <- tidy(object, model = model)
   out <- x$estimate
   names(out) <- x$term
   out

R/predict.R

@@ -45,7 +45,7 @@
 #' @param mcmc_samples See `extract_mcmc()` in the
 #'   \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package for
 #'   more details and the
-#'   \href{https://pbs-assess.github.io/sdmTMB/articles/web_only/bayesian.html}{Bayesian vignette}.
+#'   \href{https://pbs-assess.github.io/sdmTMB/articles/bayesian.html}{Bayesian vignette}.
 #'   If specified, the predict function will return a matrix of a similar form
 #'   as if `nsim > 0` but representing Bayesian posterior samples from the Stan
 #'   model.
@@ -495,6 +495,15 @@ predict.sdmTMB <- function(object, newdata = NULL,
       cli_abort("`area` should be of the same length as `nrow(newdata)` or of length 1.")
     }
 
+    # newdata, null offset in predict, and non-null in fit #372
+    if (isFALSE(nd_arg_was_null) && is.null(offset) && !all(object$offset == 0)) {
+      msg <- c(
+        "Fitted object contains an offset but the offset is `NULL` in `predict.sdmTMB()`.",
+        "Prediction will proceed assuming the offset vector is 0 in the prediction.",
+        "Specify an offset vector in `predict.sdmTMB()` to override this.")
+      cli_inform(msg)
+    }
+
     if (!is.null(offset)) {
       if (nrow(proj_X_ij[[1]]) != length(offset))
         cli_abort("Prediction offset vector does not equal number of rows in prediction dataset.")

R/tidy.R

@@ -45,7 +45,7 @@
 #' tidy(fit, "ran_vals")
 
 tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model = 1,
-                 conf.int = FALSE, conf.level = 0.95, exponentiate = FALSE,
+                 conf.int = TRUE, conf.level = 0.95, exponentiate = FALSE,
                  silent = FALSE, ...) {
   effects <- match.arg(effects)
   assert_that(is.logical(exponentiate))
@@ -131,7 +131,6 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
   b_j_se <- se$b_j[!fe_names == "offset", drop = TRUE]
   fe_names <- fe_names[!fe_names == "offset"]
   out <- data.frame(term = fe_names, estimate = b_j, std.error = b_j_se, stringsAsFactors = FALSE)
-  if (exponentiate) out$estimate <- trans(out$estimate)
 
   if (x$tmb_data$threshold_func > 0) {
     if (x$threshold_function == 1L) {
@@ -152,6 +151,9 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
     out$conf.low <- as.numeric(trans(out$estimate - crit * out$std.error))
     out$conf.high <- as.numeric(trans(out$estimate + crit * out$std.error))
   }
+  # must wrap in as.numeric() otherwise I() leaves 'AsIs' class that affects emmeans package
+  out$estimate <- as.numeric(trans(out$estimate))
+  if (exponentiate) out$std.error <- NULL
 
   out_re <- list()
   log_name <- c("log_range")

R/tmb-sim.R

@@ -268,10 +268,30 @@ sdmTMB_simulate <- function(formula,
   if (!is.null(fit$time)) d[[fit$time]] <- data[[fit$time]]
   d[[mesh$xy_cols[1]]] <- data[[mesh$xy_cols[1]]]
   d[[mesh$xy_cols[2]]] <- data[[mesh$xy_cols[2]]]
-  d[["omega_s"]] <- if (all(s$omega_s_A != 0)) s$omega_s_A
-  d[["epsilon_st"]] <- if (all(s$epsilon_st_A_vec != 0)) s$epsilon_st_A_vec
-  d[["zeta_s"]] <- if (all(s$zeta_s_A != 0)) s$zeta_s_A
-  d[["mu"]] <- family$linkinv(s$eta_i)
+
+  d[["omega_s"]] <- if (sum(sigma_O) > 0) s$omega_s_A
+  d[["epsilon_st"]] <- if (sum(sigma_E) > 0) s$epsilon_st_A_vec
+  d[["zeta_s"]] <- if (sum(sigma_Z) > 0) s$zeta_s_A
+
+  # # Warnings for fields collapsing to 0
+  # info_collapse <- function(sig, vec, .par, .name) {
+  #   if (sum(sig) > 0 && all (vec == 0)) {
+  #     msg <- paste0("The ", .name, " has been returned as all zeros although ", .par,
+  #       " was specified as > 0. Try making your mesh finer, e.g., with a lower ",
+  #       "`cutoff` or a higher number of knots. Triangle edge length needs to be ",
+  #       "lower than the range size (distance correlation is effectively independent.")
+  #     cli::cli_alert_info(msg)
+  #   }
+  # }
+  # info_collapse(sigma_O, s$omega_s_A, "sigma_O", "spatial field")
+  # info_collapse(sigma_E, s$epsilon_st_A_vec, "sigma_E", "spatiotemporal field")
+  # info_collapse(sigma_Z, s$zeta_s_A, "sigma_Z", "spatially varying coefficient field")
+
+  if (any(family$family %in% c("truncated_nbinom1", "truncated_nbinom2"))) {
+    d[["mu"]] <- family$linkinv(s$eta_i, phi = phi)
+  } else {
+    d[["mu"]] <- family$linkinv(s$eta_i)
+  }
   d[["eta"]] <- s$eta_i
   d[["observed"]] <- s$y_i
   d <- do.call("data.frame", d)
@@ -407,8 +427,7 @@ simulate.sdmTMB <- function(object, nsim = 1L, seed = sample.int(1e6, 1L),
   if (!is.null(mcmc_samples)) { # we have a matrix
     for (i in seq_len(nsim)) {
       if (!silent) cli::cli_progress_update()
-      ret[[i]] <-
-        newobj$simulate(par = new_par[, i, drop = TRUE], complete = FALSE)$y_i
+      ret[[i]] <- newobj$simulate(par = new_par[, i, drop = TRUE], complete = FALSE)$y_i
     }
   } else {
     for (i in seq_len(nsim)) {