Merge branch 'main' into cAIC

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.6.0.9013
+Version: 0.6.0.9015
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "[email protected]",
            role = c("aut", "cre"),
@@ -64,7 +64,6 @@ Depends:
 Imports:
     assertthat,
     abind,
-    clisymbols,
     cli,
     fmesher,
     fishMod,

NEWS.md

@@ -1,5 +1,13 @@
 # sdmTMB (development version)
 
+* Fix bug in `est` column when predicting on new data with Poisson-link
+  delta models with `type = "link"` and `re_form = NA`. #389
+
+* Fix bug in `s95` param reporting from the `tidy()` method. `s95` is present
+  in the logistic threshold models. The model itself was fine but the `s95`
+  parameter was supposed to be reported by `tidy()` as a combination of two 
+  other parameters. This also affected the output in `print()`/`summary()`.
+
 * Add `return_tmb_report` to `simulate.sdmTMB()`.
 
 * Add `newdata` argument to `simulate.sdmTMB()`. This enables simulating on

R/mesh.R

@@ -207,12 +207,12 @@ binary_search_knots <- function(loc_xy,
     cat("cutoff =", pretty_cutoff, "| knots =", realized_knots)
     if (realized_knots > n_knots) {
       L <- m + 1
-      cat(" |", clisymbols::symbol$arrow_down, "\n")
+      cat(" |", cli::symbol$arrow_down, "\n")
     } else if (realized_knots < n_knots) {
       R <- m - 1
-      cat(" |", clisymbols::symbol$arrow_up, "\n")
+      cat(" |", cli::symbol$arrow_up, "\n")
     } else {
-      cat(" |", clisymbols::symbol$tick, "\n")
+      cat(" |", cli::symbol$tick, "\n")
       return(mesh)
     }
   }

R/predict.R

@@ -785,14 +785,28 @@ predict.sdmTMB <- function(object, newdata = NULL,
         if (type == "response") {
           nd$est1 <- object$family[[1]]$linkinv(r$proj_fe[,1])
           nd$est2 <- object$family[[2]]$linkinv(r$proj_fe[,2])
-          nd$est <- nd$est1 * nd$est2
+          if (object$tmb_data$poisson_link_delta) {
+            .n <- nd$est1 # expected group density (already exp())
+            .p <- 1 - exp(-.n) # expected encounter rate
+            .w <- nd$est2 # expected biomass per group (already exp())
+            .r <- (.n * .w) / .p # (n * w)/p # positive expectation
+            nd$est1 <- .p # expected encounter rate
+            nd$est2 <- .r # positive expectation
+            nd$est <- .n * .w # expected combined value
+          } else {
+            nd$est <- nd$est1 * nd$est2
+          }
         } else {
           nd$est1 <- r$proj_fe[,1]
           nd$est2 <- r$proj_fe[,2]
           if (is.na(model)) {
             p1 <- object$family[[1]]$linkinv(r$proj_fe[,1])
             p2 <- object$family[[2]]$linkinv(r$proj_fe[,2])
-            nd$est <- object$family[[2]]$linkfun(p1 * p1)
+            if (object$tmb_data$poisson_link_delta) {
+              nd$est <- nd$est1 + nd$est2
+            } else {
+              nd$est <- object$family[[2]]$linkfun(p1 * p1)
+            }
             if (se_fit) {
               nd$est <- sr_est_rep$proj_rf_delta
               nd$est_se <- sr_se_rep$proj_rf_delta

R/tidy.R

@@ -135,14 +135,18 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
   if (x$tmb_data$threshold_func > 0) {
     if (x$threshold_function == 1L) {
       par_name <- paste0(x$threshold_parameter, c("-slope", "-breakpt"))
+      estimates <- est$b_threshold[,model,drop=TRUE]
+      ses <- se$b_threshold[,model,drop=TRUE]
     } else {
       par_name <- paste0(x$threshold_parameter, c("-s50", "-s95", "-smax"))
+      estimates <- c(est$s50[model], est$s95[model], est$s_max[model])
+      ses <- c(se$s50[model], se$s95[model], se$s_max[model])
     }
     out <- rbind(
       out,
       data.frame(
-        term = par_name, estimate = est$b_threshold[,model,drop=TRUE],
-        std.error = se$b_threshold[,model,drop=TRUE], stringsAsFactors = FALSE
+        term = par_name, estimate = estimates,
+        std.error = ses, stringsAsFactors = FALSE
       )
     )
   }

README.md

@@ -1,7 +1,7 @@
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
-# sdmTMB <a href='https://github.com/pbs-assess/sdmTMB'><img src='man/figures/logo-sdmTMB.png' align="right" style="height:139px;"/></a>
+# sdmTMB <a href='https://github.com/pbs-assess/sdmTMB'><img src='man/figures/logo.png' align="right" style="height:139px;"/></a>
 
 > Spatial and spatiotemporal GLMMs with TMB
 

header.md

@@ -1,6 +1,6 @@
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
-# sdmTMB <a href='https://github.com/pbs-assess/sdmTMB'><img src='man/figures/logo-sdmTMB.png' align="right" style="height:139px;"/></a>
+# sdmTMB <a href='https://github.com/pbs-assess/sdmTMB'><img src='man/figures/logo.png' align="right" style="height:139px;"/></a>
 
 > Spatial and spatiotemporal GLMMs with TMB
 

pkgdown/favicon/site.webmanifest

@@ -0,0 +1,21 @@
+{
+  "name": "",
+  "short_name": "",
+  "icons": [
+    {
+      "src": "/web-app-manifest-192x192.png",
+      "sizes": "192x192",
+      "type": "image/png",
+      "purpose": "maskable"
+    },
+    {
+      "src": "/web-app-manifest-512x512.png",
+      "sizes": "512x512",
+      "type": "image/png",
+      "purpose": "maskable"
+    }
+  ],
+  "theme_color": "#ffffff",
+  "background_color": "#ffffff",
+  "display": "standalone"
+}

scratch/binary-search.R

@@ -47,12 +47,12 @@ binary_search_knots <- function(x, y,
     cat("cutoff =", pretty_cutoff, "| knots =", realized_knots)
     if (realized_knots > n_knots) {
       L <- m + 1
-      cat(" |", clisymbols::symbol$arrow_down, "\n")
+      cat(" |", cli::symbol$arrow_down, "\n")
     } else if (realized_knots < n_knots) {
       R <- m - 1
-      cat(" |", clisymbols::symbol$arrow_up, "\n")
+      cat(" |", cli::symbol$arrow_up, "\n")
     } else {
-      cat(" |", clisymbols::symbol$tick, "\n")
+      cat(" |", cli::symbol$tick, "\n")
       return(mesh)
     }
   }

tests/testthat/test-poisson-link.R

@@ -0,0 +1,23 @@
+test_that("Poisson-link prediction issues from GitHub issues #389,", {
+  skip_on_cran()
+  dogfish$log_depth <- log(dogfish$depth)
+  mesh <- make_mesh(dogfish, c("X", "Y"), cutoff = 30)
+  fit_dpg <- sdmTMB(catch_weight ~ 0 + as.factor(year) + s(log_depth),
+    family = delta_gamma(type = "poisson-link"),
+    spatial = "on",
+    mesh = mesh,
+    data = dogfish,
+    offset = log(dogfish$area_swept)
+  )
+  p <- predict(fit_dpg, re_form = NA)
+  p$est_est1est2 <- p$est1 + p$est2
+  expect_equal(p$est, p$est_est1est2)
+  expect_equal(p$est_est1est2[1:5],
+    c(7.01028301430728, 2.34143881755487, 6.96979232578834, 6.99973559970208,
+      7.03187981132451), tolerance = 1e-3)
+
+  pp <- predict(fit_dpg, type = "response")
+  expect_equal(pp$est[1:5],
+    c(693.502617933497, 107.803298324919, 8414.54507288536, 5770.52404422525,
+    6545.06096568627), tolerance = 1e-3)
+})

tests/testthat/test-threshold-models.R

@@ -10,7 +10,7 @@ test_that("A logistic threshold model fits", {
   expect_true("depth_scaled-s50" %in% tidy(m)$term)
   expect_true("depth_scaled-s95" %in% tidy(m)$term)
   expect_true("depth_scaled-smax" %in% tidy(m)$term)
-  expect_equal(tidy(m)[,"estimate",drop=TRUE], c(1.555 , 1.655 , 1.718 , 1.138, -0.979, -3.173 , 1.760), tolerance = 1e-3)
+  expect_equal(tidy(m)[,"estimate",drop=TRUE], c(1.555 , 1.655 , 1.718 , 1.138, -0.979, -0.937 , 1.760), tolerance = 1e-3)
 })
 
 test_that("A linear threshold model fits", {
@@ -114,3 +114,4 @@ test_that("A linear threshold *delta* model fits", {
   expect_equal(t1$std.error, td1$std.error, tolerance = 1e-5)
   expect_equal(t2$std.error, td2$std.error, tolerance = 1e-5)
 })
+

vignettes/articles/ggeffects.Rmd

@@ -81,7 +81,7 @@ plot(g2)
 We can add in data points
 
 ```{r}
-plot(g, add.data = TRUE)
+plot(g, show_data = TRUE)
 ```
 
 We can also use `ggeffect` to plot multiple variables by listing them in `terms = c()`, with the first term listed indicating the variable to be plotted on the x-axis, and the remaining listed terms (up to four total) indicating the groups. 
@@ -120,7 +120,7 @@ g2 <- ggeffect(fit2, "depth_scaled [-3:2.7 by=0.05]")
 plot(g2)
 
 # note the high density values dwarf the fitted curve here
-plot(g2, add.dat = TRUE)
+plot(g2, show_data = TRUE)
 ```
 
 We can fit a model with an interaction of two continuous variables: