diff --git a/.Rbuildignore b/.Rbuildignore
index 41f8b9ec0..cac907582 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -21,3 +21,4 @@
 ^CRAN-SUBMISSION$
 ^vignettes/web_only$
 revdep/
+^vignettes/articles$
diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
index 3d1d9f05d..a47219023 100644
--- a/.github/workflows/pkgdown.yaml
+++ b/.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}
diff --git a/DESCRIPTION b/DESCRIPTION
index 03cd6b9f2..1784f9892 100644
--- a/DESCRIPTION
+++ b/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", , "sean@seananderson.ca",
            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
diff --git a/NAMESPACE b/NAMESPACE
index cab038703..91f8e64da 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
index 39fc01af7..a1a24b087 100644
--- a/NEWS.md
+++ b/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.
 
diff --git a/R/cross-val.R b/R/cross-val.R
index 06abe9af1..0264d5c6f 100644
--- a/R/cross-val.R
+++ b/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")
+}
diff --git a/R/dharma.R b/R/dharma.R
index a3f29323e..ca66b2103 100644
--- a/R/dharma.R
+++ b/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
diff --git a/R/families.R b/R/families.R
index 7b2b8ae4b..c69946524 100644
--- a/R/families.R
+++ b/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.
diff --git a/R/fit.R b/R/fit.R
index 538486fd5..c53eeb686 100644
--- a/R/fit.R
+++ b/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,
diff --git a/R/mesh.R b/R/mesh.R
index 8b97c69d4..0a071ed1b 100644
--- a/R/mesh.R
+++ b/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")
diff --git a/R/methods.R b/R/methods.R
index 58ec21638..314577abe 100644
--- a/R/methods.R
+++ b/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
diff --git a/R/predict.R b/R/predict.R
index c28b96a5f..1981545be 100644
--- a/R/predict.R
+++ b/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.")
diff --git a/R/tidy.R b/R/tidy.R
index e1dff0112..71fe69e06 100644
--- a/R/tidy.R
+++ b/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")
diff --git a/R/tmb-sim.R b/R/tmb-sim.R
index 8aa4a4ea0..d5a29c1b0 100644
--- a/R/tmb-sim.R
+++ b/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)) {
diff --git a/README.Rmd b/README.Rmd
index ca07128ca..95180c7cb 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -509,7 +509,7 @@ See [`?sdmTMBpriors`](https://pbs-assess.github.io/sdmTMB/reference/priors.html)
 
 ### Bayesian MCMC sampling with Stan
 
-The fitted model can be passed to the tmbstan package to sample from the posterior with Stan. See the [Bayesian vignette](https://pbs-assess.github.io/sdmTMB/articles/web_only/bayesian.html).
+The fitted model can be passed to the tmbstan package to sample from the posterior with Stan. See the [Bayesian vignette](https://pbs-assess.github.io/sdmTMB/articles/bayesian.html).
 
 ### Turning off random fields
 
diff --git a/README.md b/README.md
index 341bcc1dd..7147a42e9 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@
 [![Documentation](https://img.shields.io/badge/documentation-sdmTMB-orange.svg?colorB=E91E63)](https://pbs-assess.github.io/sdmTMB/)
 [![R-CMD-check](https://github.com/pbs-assess/sdmTMB/workflows/R-CMD-check/badge.svg)](https://github.com/pbs-assess/sdmTMB/actions)
 [![Codecov test coverage](https://codecov.io/gh/pbs-assess/sdmTMB/branch/main/graph/badge.svg)](https://app.codecov.io/gh/pbs-assess/sdmTMB?branch=main)
-[![downloads](http://cranlogs.r-pkg.org/badges/sdmTMB)](https://cranlogs.r-pkg.org/)
+[![downloads](https://cranlogs.r-pkg.org/badges/sdmTMB)](https://cranlogs.r-pkg.org/)
 <!-- badges: end -->
 
 sdmTMB is an R package that fits spatial and spatiotemporal GLMMs (Generalized Linear Mixed Effects Models) using Template Model Builder ([TMB](https://github.com/kaskr/adcomp)), [R-INLA](https://www.r-inla.org/), and Gaussian Markov random fields. One common application is for species distribution models (SDMs). See the [documentation site](https://pbs-assess.github.io/sdmTMB/) and a preprint:
@@ -729,7 +729,7 @@ for more details.
 
 The fitted model can be passed to the tmbstan package to sample from the
 posterior with Stan. See the [Bayesian
-vignette](https://pbs-assess.github.io/sdmTMB/articles/web_only/bayesian.html).
+vignette](https://pbs-assess.github.io/sdmTMB/articles/bayesian.html).
 
 ### Turning off random fields
 
diff --git a/_pkgdown.yml b/_pkgdown.yml
index beae68bc5..6adb1f359 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -32,15 +32,16 @@ reference:
     Core tools for model fitting, prediction, and inspection.
   contents:
   - sdmTMB
-  - predict.sdmTMB
-  - tidy.sdmTMB
   - sanity
-  - sdmTMBcontrol
-  - run_extra_optimization
+  - tidy.sdmTMB
+  - predict.sdmTMB
   - residuals.sdmTMB
   - dharma_residuals
+  - sdmTMBcontrol
+  - run_extra_optimization
   - replicate_df
   - set_delta_model
+  - coef.sdmTMB
 
 - title: 'Families'
   desc: |
diff --git a/header.md b/header.md
index 0d9224a54..a24995118 100644
--- a/header.md
+++ b/header.md
@@ -9,7 +9,7 @@
 [![Documentation](https://img.shields.io/badge/documentation-sdmTMB-orange.svg?colorB=E91E63)](https://pbs-assess.github.io/sdmTMB/)
 [![R-CMD-check](https://github.com/pbs-assess/sdmTMB/workflows/R-CMD-check/badge.svg)](https://github.com/pbs-assess/sdmTMB/actions)
 [![Codecov test coverage](https://codecov.io/gh/pbs-assess/sdmTMB/branch/main/graph/badge.svg)](https://app.codecov.io/gh/pbs-assess/sdmTMB?branch=main)
-[![downloads](http://cranlogs.r-pkg.org/badges/sdmTMB)](https://cranlogs.r-pkg.org/)
+[![downloads](https://cranlogs.r-pkg.org/badges/sdmTMB)](https://cranlogs.r-pkg.org/)
 <!-- badges: end -->
 
 sdmTMB is an R package that fits spatial and spatiotemporal GLMMs (Generalized Linear Mixed Effects Models) using Template Model Builder ([TMB](https://github.com/kaskr/adcomp)), [R-INLA](https://www.r-inla.org/), and Gaussian Markov random fields. One common application is for species distribution models (SDMs). See the [documentation site](https://pbs-assess.github.io/sdmTMB/) and a preprint:
diff --git a/man/coef.sdmTMB.Rd b/man/coef.sdmTMB.Rd
new file mode 100644
index 000000000..d0060e660
--- /dev/null
+++ b/man/coef.sdmTMB.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/methods.R
+\name{coef.sdmTMB}
+\alias{coef.sdmTMB}
+\title{Get fixed-effect coefficients}
+\usage{
+\method{coef}{sdmTMB}(object, complete = FALSE, model = 1, ...)
+}
+\arguments{
+\item{object}{The fitted sdmTMB model object}
+
+\item{complete}{Currently ignored}
+
+\item{model}{Linear predictor for delta models. Defaults to the first
+linear predictor.}
+
+\item{...}{Currently ignored}
+}
+\description{
+Get fixed-effect coefficients
+}
diff --git a/man/dharma_residuals.Rd b/man/dharma_residuals.Rd
index d531b460b..91523c805 100644
--- a/man/dharma_residuals.Rd
+++ b/man/dharma_residuals.Rd
@@ -50,7 +50,7 @@ around \code{\link[DHARMa:createDHARMa]{DHARMa::createDHARMa()}} to facilitate i
 expected distribution. This is \emph{not} the default.
 }
 \details{
-See the \href{https://pbs-assess.github.io/sdmTMB/articles/web_only/residual-checking.html}{residuals vignette}.
+See the \href{https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html}{residuals vignette}.
 
 Advantages to these residuals over the ones from the \code{\link[=residuals.sdmTMB]{residuals.sdmTMB()}}
 method are (1) they work with delta/hurdle models for the combined
diff --git a/man/make_mesh.Rd b/man/make_mesh.Rd
index 9ea51aea0..d993b5350 100644
--- a/man/make_mesh.Rd
+++ b/man/make_mesh.Rd
@@ -61,10 +61,8 @@ Distance to extend non-convex hull from data.}
 from \code{fmesher_func} (default is \code{\link[fmesher:fm_mesh_2d]{fmesher::fm_mesh_2d_inla()}}). See
 \code{mesh$mesh$n} for the number of vertices.
 
-\code{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 \code{your_mesh$mesh} to \code{inlabru::gg()}.
+\code{plot.sdmTMBmesh()}: A plot of the mesh and data points. To make your
+own \pkg{ggplot2} version, pass \code{your_mesh$mesh} to \code{inlabru::gg()}.
 }
 \description{
 Construct an SPDE mesh for use with sdmTMB.
diff --git a/man/predict.sdmTMB.Rd b/man/predict.sdmTMB.Rd
index c33cb8dd0..1745cd037 100644
--- a/man/predict.sdmTMB.Rd
+++ b/man/predict.sdmTMB.Rd
@@ -78,7 +78,7 @@ both sets of predictions are returned.}
 \item{mcmc_samples}{See \code{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 \code{nsim > 0} but representing Bayesian posterior samples from the Stan
 model.}
diff --git a/man/tidy.sdmTMB.Rd b/man/tidy.sdmTMB.Rd
index a9ae62a67..d7d03ce3e 100644
--- a/man/tidy.sdmTMB.Rd
+++ b/man/tidy.sdmTMB.Rd
@@ -8,7 +8,7 @@
   x,
   effects = c("fixed", "ran_pars", "ran_vals"),
   model = 1,
-  conf.int = FALSE,
+  conf.int = TRUE,
   conf.level = 0.95,
   exponentiate = FALSE,
   silent = FALSE,
diff --git a/src/sdmTMB.cpp b/src/sdmTMB.cpp
index 7c3891f01..3adee6db2 100644
--- a/src/sdmTMB.cpp
+++ b/src/sdmTMB.cpp
@@ -98,6 +98,30 @@ Type Link(Type eta, int link)
   return out;
 }
 
+// Modified from glmmTMB:
+/* log-prob of non-zero value in conditional distribution  */
+template<class Type>
+Type calc_log_nzprob(Type mu, Type phi, int family) {
+  Type ans, s1, s2;
+  switch (family) {
+  case truncated_nbinom1_family:
+    s2 = logspace_add(Type(0), log(phi));      // log(1. + phi(i)
+    ans = logspace_sub(Type(0), -mu / phi * s2); // 1-prob(0)
+    break;
+  case truncated_nbinom2_family:
+    s1 = log(mu);
+    // s2 := log( 1. + mu(i) / phi(i) )
+    s2 = logspace_add(Type(0), s1 - log(phi));
+    ans = logspace_sub(Type(0), -phi * s2);
+    break;
+  // case truncated_poisson_family:
+  //   ans = logspace_sub(Type(0), -mu);  // log(1-exp(-mu(i))) = P(x>0)
+  //   break;
+  default: ans = Type(0);
+  }
+  return ans;
+}
+
 // ------------------ Main TMB template ----------------------------------------
 
 template <class Type>
@@ -755,21 +779,21 @@ Type objective_function<Type>::operator()()
   Type s1, s2, s3, lognzprob, tmp_ll, ll_1, ll_2, p_mix, mix_ratio, tweedie_p, s1_large, s2_large;
 
   // calcs for mix distr. first:
-  int mix_model;
+  int pos_model;
   if (n_m > 1) {
-    mix_model = 1;
+    pos_model = 1;
   } else {
-    mix_model = 0;
+    pos_model = 0;
   }
   vector<Type> mu_i_large(n_i);
-  switch (family(mix_model)) {
+  switch (family(pos_model)) {
   case gamma_mix_family:
   case lognormal_mix_family:
   case nbinom2_mix_family: {
     p_mix = invlogit(logit_p_mix); // probability of larger event
     mix_ratio = exp(log_ratio_mix) + Type(1.); // ratio of large:small values, constrained > 1.0
     for (int i = 0; i < n_i; i++) {
-      mu_i_large(i) = exp(log(mu_i(i, mix_model)) + log(mix_ratio));  // mean of large component = mean of smaller * ratio
+      mu_i_large(i) = exp(log(mu_i(i, pos_model)) + log(mix_ratio));  // mean of large component = mean of smaller * ratio
     }
     ADREPORT(logit_p_mix);
     ADREPORT(log_ratio_mix);
@@ -1186,12 +1210,12 @@ Type objective_function<Type>::operator()()
     // for families that implement mixture models, adjust proj_eta by
     // proportion and ratio of means
     // (1 - p_mix) * mu_i(i,m) + p_mix * (mu(i,m) * mix_ratio);
-    switch (family(mix_model)) {
+    switch (family(pos_model)) {
       case gamma_mix_family:
       case lognormal_mix_family:
       case nbinom2_mix_family:
-        proj_eta.col(mix_model) = log((1. - p_mix) * exp(proj_eta.col(mix_model)) + // regular part
-               p_mix * exp(proj_eta.col(mix_model)) * mix_ratio); //large part
+        proj_eta.col(pos_model) = log((1. - p_mix) * exp(proj_eta.col(pos_model)) + // regular part
+               p_mix * exp(proj_eta.col(pos_model)) * mix_ratio); //large part
         break;
       default:
         break;
@@ -1237,6 +1261,17 @@ Type objective_function<Type>::operator()()
     vector<Type> mu_combined(n_p);
     mu_combined.setZero();
 
+    int truncated_dist;
+    switch (family(pos_model)) {
+      case truncated_nbinom1_family:
+      case truncated_nbinom2_family:
+        truncated_dist = 1;
+        break;
+      default:
+        truncated_dist = 0;
+        break;
+    }
+
     if (calc_index_totals || calc_cog || calc_eao) {
       // ------------------ Derived quantities ---------------------------------
       Type t1;
@@ -1249,14 +1284,26 @@ Type objective_function<Type>::operator()()
             // Type R1 = Type(1.) - exp(-exp(proj_eta(i,0)));
             // Type R2 = exp(proj_eta(i,0)) / R1 * exp(proj_eta(i,1))
             mu_combined(i) = exp(proj_eta(i,0) + proj_eta(i,1)); // prevent numerical issues
-          } else {
+          } else if (truncated_dist) {
+            t1 = InverseLink(proj_eta(i,0), link(0));
+            // convert from mean of *un-truncated* to mean of *truncated* distribution
+            Type log_nzprob = calc_log_nzprob(exp(proj_eta(i,1)), phi(1), family(1));
+            t2 = exp(proj_eta(i,1)) / exp(log_nzprob);
+            mu_combined(i) = t1 * t2;
+          } else  {
             t1 = InverseLink(proj_eta(i,0), link(0));
             t2 = InverseLink(proj_eta(i,1), link(1));
             mu_combined(i) = t1 * t2;
           }
           total(proj_year(i)) += mu_combined(i) * area_i(i);
         } else { // non-delta model
-          mu_combined(i) = InverseLink(proj_eta(i,0), link(0));
+          if (truncated_dist) {
+            // convert from mean of *un-truncated* to mean of *truncated* distribution
+            Type log_nzprob = calc_log_nzprob(exp(proj_eta(i,0)), phi(0), family(0));
+            mu_combined(i) = exp(proj_eta(i,0)) / exp(log_nzprob);
+          } else  {
+            mu_combined(i) = InverseLink(proj_eta(i,0), link(0));
+          }
           total(proj_year(i)) += mu_combined(i) * area_i(i);
         }
       }
diff --git a/tests/testthat/test-cross-validation.R b/tests/testthat/test-cross-validation.R
index e7f5abfe3..5b1be6639 100644
--- a/tests/testthat/test-cross-validation.R
+++ b/tests/testthat/test-cross-validation.R
@@ -12,6 +12,7 @@ test_that("Basic cross validation works", {
     data = d, mesh = spde,
     family = tweedie(link = "log"), time = "year", k_folds = 2
   )
+  print(x)
   expect_equal(class(x$sum_loglik), "numeric")
   expect_equal(x$sum_loglik, sum(x$data$cv_loglik))
   expect_equal(x$sum_loglik, sum(x$fold_loglik))
diff --git a/tests/testthat/test-mesh.R b/tests/testthat/test-mesh.R
new file mode 100644
index 000000000..b3d069adc
--- /dev/null
+++ b/tests/testthat/test-mesh.R
@@ -0,0 +1,13 @@
+test_that("make_mesh returns informative error if some coordinates are NA", {
+  d <- data.frame(x = c(NA, runif(10)), y = c(NA, runif(10)))
+  expect_error({make_mesh(d, xy_cols = c("x", "y"), cutoff = 1)}, regexp = "NA")
+
+  d <- data.frame(x = c(1, runif(10)), y = c(NA, runif(10)))
+  expect_error({make_mesh(d, xy_cols = c("x", "y"), cutoff = 1)}, regexp = "NA")
+
+  d <- data.frame(x = c(NA, runif(10)), y = c(1, runif(10)))
+  expect_error({make_mesh(d, xy_cols = c("x", "y"), cutoff = 1)}, regexp = "NA")
+
+  d <- data.frame(x = c(1, runif(10)), y = c(1, runif(10)))
+  mesh <- make_mesh(d, xy_cols = c("x", "y"), cutoff = 1)
+})
diff --git a/tests/testthat/test-methods.R b/tests/testthat/test-methods.R
index 6b31ba3bf..99a86ffe9 100644
--- a/tests/testthat/test-methods.R
+++ b/tests/testthat/test-methods.R
@@ -27,6 +27,18 @@ test_that("coef and vcov and confint work", {
   expect_true(grepl("Estimate", colnames(x))[3])
 })
 
+test_that("coef works with delta models and informs as needed", {
+  skip_on_cran()
+  fit <- sdmTMB(
+    density ~ depth,
+    data = pcod_2011, spatial = "off",
+    family = delta_gamma()
+  )
+  expect_message(x <- coef(fit), regexp = "model")
+  expect_message(x <- coef(fit, model = 1), regexp = "model")
+  expect_message(x <- coef(fit, model = 2), regexp = "model")
+})
+
 test_that("various methods work", {
   skip_on_cran()
   fit <- sdmTMB(
diff --git a/tests/testthat/test-offset.R b/tests/testthat/test-offset.R
index aa6a636f2..bb3e91202 100644
--- a/tests/testthat/test-offset.R
+++ b/tests/testthat/test-offset.R
@@ -138,6 +138,46 @@ test_that("Offset prediction matches glm()", {
   # p_glmmTMB <- predict(fit_glmmTMB, newdata = dat)
   # expect_equal(p$est, unname(p_glmmTMB))
 })
+
+test_that("offset gets passed through cross validation as expected #372", {
+  skip_on_cran()
+  dat <- subset(dogfish, catch_weight > 0)
+  expect_error(
+    x <- sdmTMB_cv(catch_weight ~ 1,
+      data = dat,
+      family = Gamma("log"), offset = log(dat$area_swept), spatial = "off"
+    ),
+    "offset"
+  )
+  set.seed(1)
+  x <- sdmTMB_cv(catch_weight ~ 1,
+    data = dat, family = Gamma("log"),
+    offset = "area_swept", spatial = "off",
+    mesh = make_mesh(dat, c("X", "Y"), cutoff = 10), k_folds = 2
+  )
+  y <- x$data[, c("catch_weight", "cv_predicted")]
+  # plot(y$catch_weight, y$cv_predicted)
+  # if offset is applied, will have unique values because an intercept-only model:
+  expect_true(length(unique(y$cv_predicted)) == 684L)
+})
+
+test_that("predicting on newdata with a non-null offset in fit but a null offset in predict informs the user appropriately", {
+  skip_on_cran()
+  dat <- subset(dogfish, catch_weight > 0)
+  fit <- sdmTMB(
+    catch_weight ~ 1,
+    data = dat,
+    family = Gamma("log"),
+    offset = "area_swept",
+    spatial = "off"
+  )
+  pred <- predict(fit)
+  pred <- predict(fit, offset = rep(0, nrow(dat)))
+  pred <- predict(fit, newdata = qcs_grid, offset = rep(0, nrow(qcs_grid)))
+  pred <- predict(fit, newdata = qcs_grid)
+  expect_message({pred <- predict(fit, newdata = qcs_grid)}, regexp = "offset")
+})
+
 # #
 # # offset/prediction setting checks:
 #
diff --git a/tests/testthat/test-tmb-simulation.R b/tests/testthat/test-tmb-simulation.R
index fe8851f5b..d6c418746 100644
--- a/tests/testthat/test-tmb-simulation.R
+++ b/tests/testthat/test-tmb-simulation.R
@@ -183,3 +183,28 @@ test_that("simulate.sdmTMB returns the right length", {
   s <- simulate(m, nsim = 2)
   expect_equal(nrow(s), nrow(pcod))
 })
+
+test_that("coarse meshes with zeros in simulation still return fields #370", {
+  set.seed(123)
+  predictor_dat <- data.frame(
+    X = runif(100), Y = runif(100),
+    a1 = rnorm(100), year = rep(1:2, each = 50))
+  mesh <- sdmTMB::make_mesh(predictor_dat, xy_cols = c("X", "Y"), n_knots = 30)
+  sim_dat <- sdmTMB::sdmTMB_simulate(
+    formula = ~ 1 + a1,
+    data = predictor_dat,
+    time = "year",
+    mesh = mesh,
+    family = gaussian(),
+    range = 0.5,
+    sigma_E = 0.1,
+    phi = 0.1,
+    sigma_O = 0.2,
+    seed = 42,
+    B = c(0.2, -0.4) # B0 = intercept, B1 = a1 slope
+  )
+  nm <- names(sim_dat)
+  expect_true("omega_s" %in% nm)
+  expect_true("epsilon_st" %in% nm)
+  expect_false("zeta_s" %in% nm)
+})
diff --git a/tests/testthat/test-truncated-dists.R b/tests/testthat/test-truncated-dists.R
new file mode 100644
index 000000000..2039e41a4
--- /dev/null
+++ b/tests/testthat/test-truncated-dists.R
@@ -0,0 +1,191 @@
+test_that("Delta-truncated NB2 works with various types of predictions", {
+  skip_on_cran()
+  skip_if_not_installed("glmmTMB")
+
+  set.seed(1)
+  d0 <- data.frame(y = rnbinom(n = 1e3, mu = 3, size = 5))
+
+  fit <- sdmTMB(
+    y ~ 1,
+    data = d0,
+    spatial = FALSE,
+    family = delta_truncated_nbinom2()
+  )
+  fit$sd_report
+
+  fit2 <- glmmTMB::glmmTMB(
+    y ~ 1,
+    ziformula = ~ 1,
+    data = d0,
+    family = glmmTMB::truncated_nbinom2()
+  )
+
+  p <- predict(fit)
+  p2 <- predict(fit2)
+  expect_equal(p2, p$est2, tolerance = 1e-4)
+
+  p <- predict(fit, type = "response")
+  p2 <- predict(fit2, type = "response")
+  expect_equal(p2, p$est, tolerance = 1e-4)
+
+  p <- predict(fit)
+  p2 <- predict(fit2, type = "zlink")
+  expect_equal(-p2, p$est1, tolerance = 1e-4)
+
+  p <- predict(fit, type = "response")
+  p2 <- predict(fit2, type = "conditional")
+  expect_equal(p2, p$est2, tolerance = 1e-4)
+})
+
+test_that("Delta-truncated NB1 works with various types of predictions", {
+  skip_on_cran()
+  skip_if_not_installed("glmmTMB")
+
+  set.seed(1)
+  d0 <- data.frame(y = rnbinom(n = 1e3, mu = 3, size = 5))
+
+  fit <- sdmTMB(
+    y ~ 1,
+    data = d0,
+    spatial = FALSE,
+    family = delta_truncated_nbinom1()
+  )
+  fit$sd_report
+
+  fit2 <- glmmTMB::glmmTMB(
+    y ~ 1,
+    ziformula = ~ 1,
+    data = d0,
+    family = glmmTMB::truncated_nbinom1()
+  )
+
+  p <- predict(fit)
+  p2 <- predict(fit2)
+  expect_equal(p2, p$est2, tolerance = 1e-4)
+
+  p <- predict(fit, type = "response")
+  p2 <- predict(fit2, type = "response")
+  expect_equal(p2, p$est, tolerance = 1e-4)
+
+  p <- predict(fit)
+  p2 <- predict(fit2, type = "zlink")
+  expect_equal(-p2, p$est1, tolerance = 1e-4)
+
+  p <- predict(fit, type = "response")
+  p2 <- predict(fit2, type = "conditional")
+  expect_equal(p2, p$est2, tolerance = 1e-4)
+})
+
+test_that("Truncated NB2 works with various types of predictions", {
+  skip_on_cran()
+  skip_if_not_installed("glmmTMB")
+  set.seed(1)
+  d0 <- data.frame(y = rnbinom(n = 1e3, mu = 3, size = 5))
+  d0 <- subset(d0, y > 0)
+  fit <- sdmTMB(
+    y ~ 1,
+    data = d0,
+    spatial = FALSE,
+    family = truncated_nbinom2()
+  )
+  fit2 <- glmmTMB::glmmTMB(
+    y ~ 1,
+    data = d0,
+    family = glmmTMB::truncated_nbinom2()
+  )
+  p <- predict(fit)
+  p2 <- predict(fit2)
+  expect_equal(p2, p$est, tolerance = 1e-4)
+  p <- predict(fit, type = "response")
+  p2 <- predict(fit2, type = "response")
+  expect_equal(p2, p$est, tolerance = 1e-4)
+})
+
+test_that("Truncated NB1 works with various types of predictions", {
+  skip_on_cran()
+  skip_if_not_installed("glmmTMB")
+  set.seed(1)
+  d0 <- data.frame(y = rnbinom(n = 1e3, mu = 3, size = 5))
+  d0 <- subset(d0, y > 0)
+  fit <- sdmTMB(
+    y ~ 1,
+    data = d0,
+    spatial = FALSE,
+    family = truncated_nbinom1()
+  )
+  fit2 <- glmmTMB::glmmTMB(
+    y ~ 1,
+    data = d0,
+    family = glmmTMB::truncated_nbinom1()
+  )
+  p <- predict(fit)
+  p2 <- predict(fit2)
+  expect_equal(p2, p$est, tolerance = 1e-4)
+  p <- predict(fit, type = "response")
+  p2 <- predict(fit2, type = "response")
+  expect_equal(p2, p$est, tolerance = 1e-4)
+})
+
+test_that("Truncated NB1/2 indexes are right", {
+  skip_on_cran()
+  set.seed(1)
+  d0 <- data.frame(y = rnbinom(n = 1e3, mu = 3, size = 5))
+  d0$year <- 1
+
+  # D-TNB2
+  fit <- sdmTMB(
+    y ~ 1,
+    data = d0,
+    time = "year",
+    spatiotemporal = "off",
+    spatial = FALSE,
+    family = delta_truncated_nbinom2()
+  )
+  p <- predict(fit, newdata = d0, return_tmb_object = TRUE)
+  x <- get_index(p, area = 1/nrow(d0))
+  pp <- predict(fit, type = "response")
+  expect_equal(x$est, pp$est[1], tolerance = 1e-3)
+
+  # D-TNB1
+  fit <- sdmTMB(
+    y ~ 1,
+    data = d0,
+    time = "year",
+    spatiotemporal = "off",
+    spatial = FALSE,
+    family = delta_truncated_nbinom1()
+  )
+  p <- predict(fit, newdata = d0, return_tmb_object = TRUE)
+  x <- get_index(p, area = 1/nrow(d0))
+  pp <- predict(fit, type = "response")
+  expect_equal(x$est, pp$est[1], tolerance = 1e-3)
+
+  # TNB2
+  dd <- subset(d0, y > 0)
+  fit <- sdmTMB(
+    y ~ 1,
+    data = dd,
+    time = "year",
+    spatiotemporal = "off",
+    spatial = FALSE,
+    family = truncated_nbinom2()
+  )
+  p <- predict(fit, newdata = d0, return_tmb_object = TRUE)
+  x <- get_index(p, area = 1/nrow(d0))
+  pp <- predict(fit, type = "response")
+  expect_equal(x$est, pp$est[1], tolerance = 1e-3)
+
+  # TNB1
+  fit <- sdmTMB(
+    y ~ 1,
+    data = dd,
+    time = "year",
+    spatiotemporal = "off",
+    spatial = FALSE,
+    family = truncated_nbinom2()
+  )
+  p <- predict(fit, newdata = d0, return_tmb_object = TRUE)
+  x <- get_index(p, area = 1/nrow(d0))
+  pp <- predict(fit, type = "response")
+  expect_equal(x$est, pp$est[1], tolerance = 1e-3)
+})
diff --git a/vignettes/web_only/basic-intro.Rmd b/vignettes/articles/basic-intro.Rmd
similarity index 100%
rename from vignettes/web_only/basic-intro.Rmd
rename to vignettes/articles/basic-intro.Rmd
diff --git a/vignettes/web_only/bayesian.Rmd b/vignettes/articles/bayesian.Rmd
similarity index 100%
rename from vignettes/web_only/bayesian.Rmd
rename to vignettes/articles/bayesian.Rmd
diff --git a/vignettes/web_only/cross-validation.Rmd b/vignettes/articles/cross-validation.Rmd
similarity index 74%
rename from vignettes/web_only/cross-validation.Rmd
rename to vignettes/articles/cross-validation.Rmd
index 2eb4669e3..c453bd271 100644
--- a/vignettes/web_only/cross-validation.Rmd
+++ b/vignettes/articles/cross-validation.Rmd
@@ -221,4 +221,82 @@ If we had just wanted to use the predictions from the first fold onto the 10% te
 weights <- sdmTMB_stacking(model_list, include_folds = 1)
 ```
 
+# Calculating measures of predictive skill for binary data
+
+For delta models, or models of presence-absence data, several measures of predictive ability are available.
+These are applicable to cross validation, although we demonstrate them here first in a non-cross validation context for simplicity.
+
+A first commonly used diagnostic is the AUC (Area Under the Curve), which quantifies the ability of a model to discriminate between the two classes; this is done from the Receiver Operating Characteristic (ROC) curve, which plots the true positive rate vs. false positive rate.
+There are several packages to calculate AUC in R, but this can be done with the `pROC` package, where inputs are a vector of 0s and 1s (or factor equivalents) in the raw data, and a vector of estimated probabilities (generated from a call to `predict()`, as shown below).
+The `plogis()` function is needed to convert estimated values in logit space to probabilities in natural (zero to one) space. 
+
+```{r roc}
+mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 10)
+fit <- sdmTMB(present ~ s(depth), data = pcod, mesh = mesh)
+pred <- predict(fit) # presence-absence model
+roc <- pROC::roc(pcod$present, plogis(pred$est))
+auc <- pROC::auc(roc)
+auc
+```
+
+With a delta model, two estimated values are returned, so only the first would be used. E.g.,
+
+```{r}
+fit <- sdmTMB(density ~ 1, data = pcod, 
+  mesh = mesh, family = delta_gamma())
+pred <- predict(fit)
+
+# the first linear predictor is the binomial component (est1):
+roc <- pROC::roc(pcod$present, plogis(pred$est1))
+auc <- pROC::auc(roc)
+auc
+```
+
+If we wanted to apply this in the context of cross validation, we could do it like this:
+
+```{r, eval=FALSE}
+x <- sdmTMB_cv(
+  present ~ s(depth), data = pcod, spatial = "off",
+  mesh = mesh, family = binomial(), k_folds = 2
+)
+roc <- pROC::roc(x$data$present, plogis(x$data$cv_predicted))
+auc <- pROC::auc(roc)
+auc
+```
+
+AUC may be sensitive to imbalances in the data, however, and alternative metrics may better approximate skill.
+Here we highlight an example of using true skill score (implemented in packages such as SDMtune):
+
+```{r}
+mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 10)
+fit <- sdmTMB(present ~ 1, data = pcod, 
+  mesh = mesh, family = binomial())
+```
+
+Next, we can generate predicted probabilities and classes using a threshold of 0.5 as an example:
+
+```{r}
+pred <- predict(fit)
+pred$p <- plogis(pred$est)
+pred$pred_01 <- ifelse(pred$p < 0.5, 0, 1)
+```
+
+Next we create a confusion matrix and calculate the true skill score:
+
+```{r}
+conmat <- table(pred$pred_01, pred$present)
+true_neg <- conmat[1, 1]
+false_neg <- conmat[1, 2]
+false_pos <- conmat[2, 1]
+true_pos <- conmat[2, 2]
+
+# Calculate TSS:
+true_pos_rate <- true_pos / (true_pos + false_neg)
+true_neg_rate <- true_neg / (true_neg + false_pos)
+TSS <- true_pos_rate + true_neg_rate - 1
+TSS
+```
+
+In some cases, reporting the true negative or true positive rate might be of interest in addition to TSS.
+
 # References
diff --git a/vignettes/web_only/delta-models.Rmd b/vignettes/articles/delta-models.Rmd
similarity index 99%
rename from vignettes/web_only/delta-models.Rmd
rename to vignettes/articles/delta-models.Rmd
index bfb4006a8..f15b21dc6 100644
--- a/vignettes/web_only/delta-models.Rmd
+++ b/vignettes/articles/delta-models.Rmd
@@ -47,7 +47,7 @@ Here and with other delta models, the `link1` and `link2` can be omitted and lef
 2. Delta-lognormal: `family = delta_lognormal()`.
 This fits a binomial presence-absence model (i.e., `binomial(link = "logit")`) and then a model for the positive catches only with a lognormal observation distribution and a log link (i.e., `lognormal(link = "log")`
 
-3. Poisson-link delta-Gamma or delta-lognormal. See the [Poisson-link delta model vignette](https://pbs-assess.github.io/sdmTMB/articles/web_only/poisson-link.html).
+3. Poisson-link delta-Gamma or delta-lognormal. See the [Poisson-link delta model vignette](https://pbs-assess.github.io/sdmTMB/articles/poisson-link.html).
 
 4. Delta-truncated-negative-binomial: `family = delta_truncated_nbinom1()` or `family = delta_truncated_nbinom2()`.
 This fits a binomial presence-absence model (`binomial(link = "logit")`) and a `truncated_nbinom1(link = "log")` or `truncated_nbinom1(link = "log")` distribution for positive catches.
diff --git a/vignettes/web_only/ggeffects.Rmd b/vignettes/articles/ggeffects.Rmd
similarity index 100%
rename from vignettes/web_only/ggeffects.Rmd
rename to vignettes/articles/ggeffects.Rmd
diff --git a/vignettes/web_only/index-standardization.Rmd b/vignettes/articles/index-standardization.Rmd
similarity index 100%
rename from vignettes/web_only/index-standardization.Rmd
rename to vignettes/articles/index-standardization.Rmd
diff --git a/vignettes/web_only/poisson-link.Rmd b/vignettes/articles/poisson-link.Rmd
similarity index 100%
rename from vignettes/web_only/poisson-link.Rmd
rename to vignettes/articles/poisson-link.Rmd
diff --git a/vignettes/web_only/pretty-plots.Rmd b/vignettes/articles/pretty-plots.Rmd
similarity index 100%
rename from vignettes/web_only/pretty-plots.Rmd
rename to vignettes/articles/pretty-plots.Rmd
diff --git a/vignettes/web_only/residual-checking.Rmd b/vignettes/articles/residual-checking.Rmd
similarity index 100%
rename from vignettes/web_only/residual-checking.Rmd
rename to vignettes/articles/residual-checking.Rmd
diff --git a/vignettes/web_only/spatial-trend-models.Rmd b/vignettes/articles/spatial-trend-models.Rmd
similarity index 100%
rename from vignettes/web_only/spatial-trend-models.Rmd
rename to vignettes/articles/spatial-trend-models.Rmd
diff --git a/vignettes/web_only/threshold-models.Rmd b/vignettes/articles/threshold-models.Rmd
similarity index 100%
rename from vignettes/web_only/threshold-models.Rmd
rename to vignettes/articles/threshold-models.Rmd
diff --git a/vignettes/web_only/visreg.Rmd b/vignettes/articles/visreg.Rmd
similarity index 100%
rename from vignettes/web_only/visreg.Rmd
rename to vignettes/articles/visreg.Rmd