Add priors to logistic() and breakpt() models

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.6.0.9025
+Version: 0.6.0.9026
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca",
            role = c("aut", "cre"),

NEWS.md

@@ -70,6 +70,8 @@
 
 * Add time varying AR(1) correlation to tidy() and print() #374
 
+* Add priors for `breakpt()` and `logistic()` parameters
+
 # sdmTMB 0.6.0
 
 * Pass several arguments to `DHARMa::plotQQunif()`.

R/priors.R

@@ -62,6 +62,13 @@
 #'   parameter has support `1 < tweedie_p < 2` so choose a mean appropriately.
 #' @param b `normal()` priors for the main population-level 'beta' effects.
 #' @param sigma_G `halfnormal()` priors for the random intercept SDs.
+#' @param threshold_breakpt_slope prior for the slope of the linear (hockey stick)
+#' function
+#' @param threshold_breakpt_cut prior for the cutoff of the linear (hockey stick)
+#' function
+#' @param threshold_logistic_s50 prior for which f(x) = 0.5
+#' @param threshold_logistic_s95 prior for which f(x) = 0.95
+#' @param threshold_logistic_smax prior for the value for which f(x) is maximized
 #'
 #' @rdname priors
 #'
@@ -76,7 +83,12 @@ sdmTMBpriors <- function(
   ar1_rho = normal(NA, NA),
   tweedie_p = normal(NA, NA),
   b = normal(NA, NA),
-  sigma_G = halfnormal(NA, NA)
+  sigma_G = halfnormal(NA, NA),
+  threshold_breakpt_slope = normal(NA, NA),
+  threshold_breakpt_cut = normal(NA, NA),
+  threshold_logistic_s50 = normal(NA, NA),
+  threshold_logistic_s95 = normal(NA, NA),
+  threshold_logistic_smax = normal(NA, NA)
 ) {
   assert_that(attr(matern_s, "dist") == "pc_matern")
   assert_that(attr(matern_st, "dist") == "pc_matern")
@@ -91,7 +103,12 @@ sdmTMBpriors <- function(
     ar1_rho = ar1_rho,
     tweedie_p = tweedie_p,
     b = b,
-    sigma_G = sigma_G
+    sigma_G = sigma_G,
+    threshold_breakpt_slope = threshold_breakpt_slope,
+    threshold_breakpt_cut = threshold_breakpt_cut,
+    threshold_logistic_s50 = threshold_logistic_s50,
+    threshold_logistic_s95 = threshold_logistic_s95,
+    threshold_logistic_smax = threshold_logistic_smax
   )
 }
 

src/sdmTMB.cpp

@@ -1035,17 +1035,32 @@ Type objective_function<Type>::operator()()
       // log abs derivative = log((2 * exp(x)) / (1 + exp(x))^2)
       if (stan_flag) jnll -= log(2.) + ar1_phi(m) - 2. * log(1. + exp(ar1_phi(m)));
     }
+    if(!sdmTMB::isNA(priors(14)) && !sdmTMB::isNA(priors(15))) { // hockeystick
+      jnll -= dnorm(s_slope(m), priors(14), priors(15), true);
+    }
+    if(!sdmTMB::isNA(priors(16)) && !sdmTMB::isNA(priors(17))) { // hockeystick
+      jnll -= dnorm(s_cut(m), priors(16), priors(17), true);
+    }
+    if(!sdmTMB::isNA(priors(18)) && !sdmTMB::isNA(priors(19))) { // logistic
+      jnll -= dnorm(s50(m), priors(18), priors(19), true);
+    }
+    if(!sdmTMB::isNA(priors(20)) && !sdmTMB::isNA(priors(21))) { // logistic
+      jnll -= dnorm(s95(m), priors(20), priors(21), true);
+    }
+    if(!sdmTMB::isNA(priors(22)) && !sdmTMB::isNA(priors(23))) { // logistic
+      jnll -= dnorm(s_max(m), priors(22), priors(23), true);
+    }
     if (priors_sigma_G.rows() != sigma_G.rows())
       error("sigma_G prior dimensions are incorrect");
-    for (int m = 0; m < n_m; m++) {
-      for (int g = 0; g < sigma_G.rows(); g++) {
-        if (!sdmTMB::isNA(priors_sigma_G(g,0)) && !sdmTMB::isNA(priors_sigma_G(g,1))) {
-          jnll -= dnorm(sigma_G(g,m), priors_sigma_G(g,0), priors_sigma_G(g,1), true);
-          if (stan_flag) jnll -= log(sigma_G(g,m)); // Jacobian adjustment
+      for (int m = 0; m < n_m; m++) {
+        for (int g = 0; g < sigma_G.rows(); g++) {
+          if (!sdmTMB::isNA(priors_sigma_G(g,0)) && !sdmTMB::isNA(priors_sigma_G(g,1))) {
+            jnll -= dnorm(sigma_G(g,m), priors_sigma_G(g,0), priors_sigma_G(g,1), true);
+            if (stan_flag) jnll -= log(sigma_G(g,m)); // Jacobian adjustment
+          }
         }
       }
     }
-  }
 
   // ------------------ Predictions on new data --------------------------------
 

tests/testthat/test-priors.R

@@ -163,3 +163,54 @@ test_that("Additional priors work", {
       abs(tidy(m_mvn_na)$estimate[tidy(m_mvn_na)$term == "depth_scaled"])
   )
 })
+
+test_that("Threshold priors work", {
+  skip_on_cran()
+  d <- pcod_2011
+  pcod_spde <- pcod_mesh_2011
+
+  m <- sdmTMB(density ~ 0 + logistic(depth_scaled),
+                   data = d, mesh = pcod_spde,
+                   family = tweedie(link = "log"),
+                   priors = sdmTMBpriors(threshold_logistic_s50 = normal(-1,0.005)),
+                   #priors = sdmTMBpriors(b = normal(0, 10)),
+                   spatial = "off", spatiotemporal = "off"
+  )
+  expect_equal(tidy(m)$estimate[1], -1, tolerance = 0.01)
+
+  m <- sdmTMB(density ~ 0 + logistic(depth_scaled),
+              data = d, mesh = pcod_spde,
+              family = tweedie(link = "log"),
+              priors = sdmTMBpriors(threshold_logistic_s95 = normal(-1,0.005)),
+              #priors = sdmTMBpriors(b = normal(0, 10)),
+              spatial = "off", spatiotemporal = "off"
+  )
+  expect_equal(tidy(m)$estimate[2], -1, tolerance = 0.01)
+
+  m <- sdmTMB(density ~ 0 + logistic(depth_scaled),
+              data = d, mesh = pcod_spde,
+              family = tweedie(link = "log"),
+              priors = sdmTMBpriors(threshold_logistic_smax = normal(4.0,0.005)),
+              #priors = sdmTMBpriors(b = normal(0, 10)),
+              spatial = "off", spatiotemporal = "off"
+  )
+  expect_equal(tidy(m)$estimate[3], 4.0, tolerance = 0.01)
+
+  m <- sdmTMB(density ~ 0 + breakpt(depth_scaled),
+              data = d, mesh = pcod_spde,
+              family = tweedie(link = "log"),
+              priors = sdmTMBpriors(threshold_breakpt_slope = normal(-4,0.005)),
+              #priors = sdmTMBpriors(b = normal(0, 10)),
+              spatial = "off", spatiotemporal = "off"
+  )
+  expect_equal(tidy(m)$estimate[1], -4, tolerance = 0.01)
+
+  m <- sdmTMB(density ~ 0 + breakpt(depth_scaled),
+              data = d, mesh = pcod_spde,
+              family = tweedie(link = "log"),
+              priors = sdmTMBpriors(threshold_breakpt_cut = normal(-1,0.005)),
+              #priors = sdmTMBpriors(b = normal(0, 10)),
+              spatial = "off", spatiotemporal = "off"
+  )
+  expect_equal(tidy(m)$estimate[2], -1, tolerance = 0.01)
+})