Move experimental(upr) to sdmTMBcontrol()

R/fit.R

@@ -662,20 +662,8 @@ sdmTMB <- function(
     } else {
       epsilon_model <- NULL
     }
-
-    if ("lwr" %in% names(experimental) && "upr" %in% names(experimental)) {
-      lwr <- experimental$lwr
-      upr <- experimental$upr
-    } else {
-      #epsilon_predictor <- NULL
-      lwr <- 0
-      upr <- Inf
-    }
-  } else {
-    #epsilon_predictor <- NULL
-    lwr <- 0
-    upr <- Inf
   }
+
   normalize <- control$normalize
   nlminb_loops <- control$nlminb_loops
   newton_loops <- control$newton_loops
@@ -686,8 +674,9 @@ sdmTMB <- function(
   lower <- control$lower
   upper <- control$upper
   get_joint_precision <- control$get_joint_precision
+  upr <- control$censored_upper
 
-  dot_checks <- c("lower", "upper", "profile", "parallel",
+  dot_checks <- c("lower", "upper", "profile", "parallel", "censored_upper",
     "nlminb_loops", "newton_steps", "mgcv", "quadratic_roots", "multiphase",
     "newton_loops", "start", "map", "get_joint_precision", "normalize")
   .control <- control
@@ -715,13 +704,6 @@ sdmTMB <- function(
     cli_warn(c("`length(map) != length(start)`.",
       "You likely want to specify `start` values if you are setting the `map` argument."))
   }
-  if (family$family[1] == "censored_poisson") {
-    assert_that("lwr" %in% names(experimental) && "upr" %in% names(experimental),
-      msg = "`lwr` and `upr` must be specified in `experimental` as elements of a named list to use the censored Poisson likelihood.")
-    assert_that(length(lwr) == nrow(data) && length(upr) == nrow(data))
-    assert_that(length(lwr) == length(upr))
-    assert_that(mean(upr-lwr, na.rm = TRUE)>=0)
-  }
 
   if (!is.null(time)) {
     assert_that(time %in% names(data),
@@ -752,6 +734,15 @@ sdmTMB <- function(
   }
   # FIXME parallel setup here?
 
+  if (family$family[1] == "censored_poisson") {
+    if ("lwr" %in% names(experimental) || "upr" %in% names(experimental)) {
+      cli_abort("Detected `lwr` or `upr` in `experimental`. `lwr` is no longer needed and `upr` is now specified as `control = sdmTMBcontrol(censored_upper = ...)`.")
+    }
+    if (is.null(upr)) cli_abort("`censored_upper` must be defined in `control = sdmTMBcontrol()` to use the censored Poisson distribution.")
+    assert_that(length(upr) == nrow(data))
+  }
+  if (is.null(upr)) upr <- Inf
+
   # thresholds not yet enabled for delta model, where formula is a list
   if (inherits(formula, "formula")) {
     original_formula <- replicate(n_m, list(formula))
@@ -773,12 +764,11 @@ sdmTMB <- function(
   }
 
   if (!is.null(extra_time)) { # for forecasting or interpolating
-    data <- expand_time(df = data, time_slices = extra_time, time_column = time, weights = weights, offset = offset, upr = upr, lwr = lwr)
+    data <- expand_time(df = data, time_slices = extra_time, time_column = time, weights = weights, offset = offset, upr = upr)
     if (!is.null(offset)) offset <- data[["__sdmTMB_offset__"]] # expanded
     if (!is.null(weights)) weights <- data[["__weight_sdmTMB__"]] # expanded
     if (!is.null(upr)) upr <- data[["__dcens_upr__"]] # expanded
-    if (!is.null(lwr)) lwr <- data[["__dcens_lwr__"]] # expanded
-    data[["__dcens_upr__"]] <- data[["__dcens_lwr__"]] <- NULL
+    data[["__dcens_upr__"]] <- NULL
     spde$loc_xy <- as.matrix(data[,spde$xy_cols,drop=FALSE])
     spde$A_st <- INLA::inla.spde.make.A(spde$mesh, loc = spde$loc_xy)
     spde$sdm_spatial_id <- seq(1, nrow(data)) # FIXME
@@ -898,6 +888,9 @@ sdmTMB <- function(
   if (family$family[1] %in% c("Gamma", "lognormal") && min(y_i) <= 0 && !delta) {
     cli_abort("Gamma and lognormal must have response values > 0.")
   }
+  if (family$family[1] == "censored_poisson") {
+    assert_that(mean(upr-y_i, na.rm = TRUE)>=0)
+  }
 
   # This is taken from approach in glmmTMB to match how they handle binomial
   # yobs could be a factor -> treat as binary following glm
@@ -1126,7 +1119,7 @@ sdmTMB <- function(
     est_epsilon_re = as.integer(est_epsilon_re),
     has_smooths = as.integer(sm$has_smooths),
     upr = upr,
-    lwr = lwr,
+    lwr = 0L, # in case we want to reintroduce this
     poisson_link_delta = as.integer(isTRUE(family$type == "poisson_link_delta")),
     stan_flag = as.integer(bayesian),
     no_spatial = no_spatial

R/utils.R

@@ -49,6 +49,10 @@
 #'   cases) can be specified as a numeric vector. E.g.
 #'   `lower = list(b_j = c(-5, -5))`.
 #' @param upper An optional named list of upper bounds within the optimization.
+#' @param censored_upper An optional vector of upper bounds for
+#'   [sdmTMBcontrol()]. Values of `NA` indicate an unbounded right-censored to
+#'   distribution, values greater that the observation indicate and upper bound,
+#'   and values equal to the observation indicate no censoring.
 #' @param get_joint_precision Logical. Passed to `getJointPrecision` in
 #'   [TMB::sdreport()]. Must be `TRUE` to use simulation-based methods in
 #'   [predict.sdmTMB()] or `[get_index_sims()]`. If not needed, setting this
@@ -57,15 +61,6 @@
 #'   cores, as an example, use `TMB::openmp(n = 3, DLL = "sdmTMB")`. But be
 #'   careful, because it's not always faster with more cores and there is
 #'   definitely an upper limit.
-# Number of cores to use.
-# Best set with the option
-#   `options(sdmTMB.cores = n)`, where `n` is the number of cores.
-#   Currently only works on Mac/Linux. Note that the models in sdmTMB
-#   often slow down with too many cores. Ideal numbers appear to be
-#   a bit less than half the available cores or ~3-4 on the machines
-#   we have tested. Also propogates to prediction and index calculation.
-#   Can we tweaked after the fact in a fitted model by modifying
-#   `fit$control$parallel`.
 #' @param ... Anything else. See the 'Control parameters' section of
 #'   [stats::nlminb()].
 #'
@@ -75,7 +70,7 @@
 #' Usually used within [sdmTMB()]. For example:
 #'
 #' ```
-#' sdmTMB(..., control = sdmTMBcontrol(newton_loops = 1))
+#' sdmTMB(..., control = sdmTMBcontrol(newton_loops = 2))
 #' ```
 #' @examples
 #' sdmTMBcontrol()
@@ -92,6 +87,7 @@ sdmTMBcontrol <- function(
   map = NULL,
   lower = NULL,
   upper = NULL,
+  censored_upper = NULL,
   multiphase = TRUE,
   profile = FALSE,
   get_joint_precision = TRUE,
@@ -130,6 +126,7 @@ sdmTMBcontrol <- function(
     map,
     lower,
     upper,
+    censored_upper,
     multiphase,
     parallel,
     get_joint_precision
@@ -231,16 +228,14 @@ parse_threshold_formula <- function(formula, thresh_type_short = "lin_thresh",
   list(formula = formula, threshold_parameter = threshold_parameter)
 }
 
-expand_time <- function(df, time_slices, time_column, weights, offset, upr, lwr) {
+expand_time <- function(df, time_slices, time_column, weights, offset, upr) {
   if (!is.null(weights)) df[["__weight_sdmTMB__"]] <- weights
   if (!is.null(offset)) df[["__sdmTMB_offset__"]] <- offset
   if (!is.null(upr)) df[["__dcens_upr__"]] <- upr
-  if (!is.null(lwr)) df[["__dcens_lwr__"]] <- lwr
   fake_df <- df[1L, , drop = FALSE]
   if (!is.null(weights)) fake_df[["__weight_sdmTMB__"]] <- 0
   if (!is.null(offset))fake_df[["__sdmTMB_offset__"]] <- 0
   if (!is.null(upr)) fake_df[["__dcens_upr__"]] <- NA_real_
-  if (!is.null(lwr)) fake_df[["__dcens_lwr__"]] <- 0
   missing_years <- time_slices[!time_slices %in% df[[time_column]]]
   fake_df <- do.call("rbind", replicate(length(missing_years), fake_df, simplify = FALSE))
   fake_df[[time_column]] <- missing_years

tests/testthat/test-3-families.R

@@ -226,74 +226,83 @@ test_that("Censored Poisson fits", {
   m_nocens_pois <- sdmTMB(
     data = sim_dat, formula = observed ~ 1,
     mesh = mesh, family = censored_poisson(link = "log"),
-    experimental = list(upr = sim_dat$observed, lwr = sim_dat$observed)
+    control = sdmTMBcontrol(censored_upper = sim_dat$observed)
   )
-  expect_equal(m_nocens_pois$tmb_data$lwr, m_nocens_pois$tmb_data$upr)
-  expect_equal(m_nocens_pois$tmb_data$lwr, as.numeric(m_nocens_pois$tmb_data$y_i))
+  expect_equal(m_nocens_pois$tmb_data$y_i[,1], m_nocens_pois$tmb_data$upr)
   expect_equal(names(m_nocens_pois$tmb_data$family), "censored_poisson")
   expect_equal(m_pois$model, m_nocens_pois$model)
 
-  # left-censored version
-  L_1 <- 5 # zeros and ones cannot be observed directly - observed as <= L1
-  y <- sim_dat$observed
-  lwr <- ifelse(y <= L_1, 0, y)
-  upr <- ifelse(y <= L_1, L_1, y)
-  m_left_cens_pois <- sdmTMB(
-    data = sim_dat, formula = observed ~ 1,
-    mesh = mesh, family = censored_poisson(link = "log"),
-    experimental = list(lwr = lwr, upr = upr),
-    spatial = "off"
-  )
+  # # left-censored version
+  # L_1 <- 5 # zeros and ones cannot be observed directly - observed as <= L1
+  # y <- sim_dat$observed
+  # lwr <- ifelse(y <= L_1, 0, y)
+  # upr <- ifelse(y <= L_1, L_1, y)
+  # m_left_cens_pois <- sdmTMB(
+  #   data = sim_dat, formula = observed ~ 1,
+  #   mesh = mesh, family = censored_poisson(link = "log"),
+  #   control = sdmTMBcontrol(censored_lower = lwr, censored_upper = upr),
+  #   spatial = "off"
+  # )
 
   # right-censored version
   U_1 <- 8 # U_1 and above cannot be directly observed - instead we see >= U1
   y <- sim_dat$observed
   lwr <- ifelse(y >= U_1, U_1, y)
   upr <- ifelse(y >= U_1, NA, y)
-  m_right_cens_pois <- sdmTMB(
+
+  # old:
+  expect_error(m_right_cens_pois <- sdmTMB(
     data = sim_dat, formula = observed ~ 1,
-    mesh = mesh, family = censored_poisson(link = "log"),
+    family = censored_poisson(link = "log"),
     experimental = list(lwr = lwr, upr = upr),
     spatial = "off"
+  ), regexp = "upr")
+
+  # new:
+  m_right_cens_pois <- sdmTMB(
+    data = sim_dat, formula = observed ~ 1,
+    family = censored_poisson(link = "log"),
+    control = sdmTMBcontrol(censored_upper = upr),
+    spatial = "off"
   )
 
   # interval-censored tough example
   # unique bounds per observation with upper limit 500 to test numerical underflow issues
   set.seed(123)
   U_2 <- sample(c(5:9), size = length(y), replace = TRUE)
   L_2 <- sample(c(1, 2, 3, 4), size = length(y), replace = TRUE)
-  lwr <- ifelse(y >= U_2, U_2, ifelse(y <= L_2, 0, y))
+  # lwr <- ifelse(y >= U_2, U_2, ifelse(y <= L_2, 0, y))
   upr <- ifelse(y >= U_2, 500, ifelse(y <= L_2, L_2, y))
   m_interval_cens_pois <- sdmTMB(
     data = sim_dat, formula = observed ~ 1,
-    mesh = mesh, family = censored_poisson(link = "log"),
-    experimental = list(lwr = lwr, upr = upr),
+    family = censored_poisson(link = "log"),
+    control = sdmTMBcontrol(censored_upper = upr),
     spatial = "off"
   )
   expect_true(all(!is.na(summary(m_interval_cens_pois$sd_report)[, "Std. Error"])))
 
-  # reversed upr and lwr:
-  expect_error(
-    m <- sdmTMB(
-      data = sim_dat, formula = observed ~ 1,
-      mesh = mesh, family = censored_poisson(link = "log"),
-      experimental = list(lwr = upr, upr = lwr)
-    ), regexp = "lwr")
+  # # reversed upr and lwr:
+  # expect_error(
+  #   m <- sdmTMB(
+  #     data = sim_dat, formula = observed ~ 1,
+  #     mesh = mesh, family = censored_poisson(link = "log"),
+  #     experimental = list(lwr = upr, upr = lwr)
+  #   ), regexp = "lwr")
 
   # wrong length lwr and upr
   expect_error(
     m <- sdmTMB(
       data = sim_dat, formula = observed ~ 1,
       mesh = mesh, family = censored_poisson(link = "log"),
-      experimental = list(lwr = c(1, 2), upr = c(4, 5, 6))
-    ), regexp = "lwr")
+      control = sdmTMBcontrol(censored_upper = c(4, 5, 6))
+    ), regexp = "upr")
 
   # missing lwr/upr
   expect_error(
     m <- sdmTMB(
       data = sim_dat, formula = observed ~ 1,
       mesh = mesh, family = censored_poisson(link = "log"),
-    ), regexp = "lwr")
+    ), regexp = "censored_upper")
 
 })