Merge branch 'zeta-intercept'

R/fit.R

@@ -81,9 +81,15 @@ NULL
 #'   used, it will represent a local-time-trend model. See
 #'   \doi{10.1111/ecog.05176} and the [spatial trends
 #'   vignette](https://pbs-assess.github.io/sdmTMB/articles/spatial-trend-models.html).
-#'   Note this predictor should be centered to have mean zero and have a
+#'   Note this predictor should usually be centered to have mean zero and have a
 #'   standard deviation of approximately 1 and should likely also be included as
-#'   a main effect. Structure must currently be shared in delta models.
+#'   a main effect.
+#'   **The spatial intercept is controlled by the `spatial` argument**; therefore,
+#'   include or exclude the spatial intercept by setting `spatial = 'on'` or
+#'   `'off'`. The only time when it matters whether `spatial_varying` excludes
+#'   an intercept is in the case of factor predictors. In this case, if
+#'   `spatial_varying` excludes the intercept (`~ 0` or `~ -1`), you should set
+#'   `spatial = 'off'` to match.  Structure must be shared in delta models.
 #' @param weights A numeric vector representing optional likelihood weights for
 #'   the conditional model. Implemented as in \pkg{glmmTMB}: weights do not have
 #'   to sum to one and are not internally modified. Can also be used for trials
@@ -573,7 +579,7 @@ sdmTMB <- function(
   ) {
 
   data <- droplevels(data) # if data was subset, strips absent factors
-  
+
   delta <- isTRUE(family$delta)
   n_m <- if (delta) 2L else 1L
 
@@ -616,6 +622,15 @@ sdmTMB <- function(
   }
   include_spatial <- "on" %in% spatial
 
+  if (!include_spatial && !is.null(spatial_varying)) {
+    # move intercept into spatial_varying
+    omit_spatial_intercept <- TRUE
+    include_spatial <- TRUE
+    spatial <- "on" # checked later
+  } else {
+    omit_spatial_intercept <- FALSE
+  }
+
   if (!include_spatial && all(spatiotemporal == "off") || !include_spatial && all(spatial_only)) {
     # message("Both spatial and spatiotemporal fields are set to 'off'.")
     # control$map_rf <- TRUE
@@ -771,6 +786,27 @@ sdmTMB <- function(
   spatial_varying_formula <- spatial_varying # save it
   if (!is.null(spatial_varying)) {
     z_i <- model.matrix(spatial_varying, data)
+    .int <- sum(grep("(Intercept)", colnames(z_i)) > 0)
+    if (length(attr(z_i, "contrasts")) && !.int && !omit_spatial_intercept) { # factors with ~ 0 or ~ -1
+      msg <- c("Detected predictors with factor levels in `spatial_varying` with the intercept omitted from the `spatial_varying` formula.",
+        "You likely want to set `spatial = 'off'` since the constant spatial field (`omega_s`) also represents a spatial intercept.`")
+        # "As of version 0.3.1, sdmTMB turns off the constant spatial field `omega_s` when `spatial_varying` is specified so that the intercept or factor-level means are fully described by the spatially varying random fields `zeta_s`.")
+      cli_inform(paste(msg, collapse = " "))
+    }
+    # if (!omit_spatial_intercept & !length(attr(z_i, "contrasts"))) {
+    #   msg <- c("The spatial intercept is now in the first element of the spatially varying random fields `zeta_s` instead of the constant spatial random field `omega_s`. This change in the output format occurred in version 0.3.1.")
+    #   cli_inform(msg)
+    # }
+#    .int <- sum(grep("(Intercept)", colnames(z_i)) > 0)
+#    if (.int && !omit_spatial_intercept) {
+#      # msg <- c("Detected an intercept in `spatial_varying`.",
+#      #   "Make sure you have `spatial = 'off'` set since this also represents a spatial intercept.")
+#      # cli_warn(msg)
+#      # actually, just do it!
+#      omit_spatial_intercept <- TRUE
+#      include_spatial <- TRUE
+#      spatial <- "on"
+#    }
     .int <- grep("(Intercept)", colnames(z_i))
     if (sum(.int) > 0) z_i <- z_i[,-.int,drop=FALSE]
     spatial_varying <- colnames(z_i)
@@ -1044,6 +1080,7 @@ sdmTMB <- function(
     priors = as.numeric(unlist(.priors)),
     share_range = as.integer(if (length(share_range) == 1L) rep(share_range, 2L) else share_range),
     include_spatial = as.integer(include_spatial),
+    omit_spatial_intercept = as.integer(omit_spatial_intercept),
     proj_mesh  = Matrix::Matrix(c(0,0,2:0), 3, 5), # dummy
     proj_X_ij  = list(matrix(0, ncol = 1, nrow = 1)), # dummy
     proj_X_rw_ik = matrix(0, ncol = 1, nrow = 1), # dummy
@@ -1105,7 +1142,7 @@ sdmTMB <- function(
     ln_tau_G   = matrix(0, ncol(RE_indexes), n_m),
     RE         = matrix(0, sum(nobs_RE), n_m),
     b_rw_t     = array(0, dim = c(tmb_data$n_t, ncol(X_rw_ik), n_m)),
-    omega_s    = matrix(0, n_s, n_m),
+    omega_s    = matrix(0, if (!omit_spatial_intercept) n_s else 0L, n_m),
     zeta_s    = array(0, dim = c(n_s, n_z, n_m)),
     epsilon_st = array(0, dim = c(n_s, tmb_data$n_t, n_m)),
     b_threshold = if(thresh[[1]]$threshold_func == 2L) matrix(0, 3L, n_m) else matrix(0, 2L, n_m),
@@ -1194,7 +1231,7 @@ sdmTMB <- function(
   }
 
   tmb_random <- c()
-  if (any(spatial == "on")) {
+  if (any(spatial == "on") && !omit_spatial_intercept) {
     tmb_random <- c(tmb_random, "omega_s")
     tmb_map <- unmap(tmb_map, c("omega_s", "ln_tau_O"))
   }

R/predict.R

@@ -292,7 +292,7 @@ predict.sdmTMB <- function(object, newdata = object$data,
   type <- match.arg(type)
   # FIXME parallel setup here?
 
-  if (!is.null(re_form) && isTRUE(se_fit)) {
+  if (is.null(re_form) && isTRUE(se_fit)) {
     msg <- paste0("Prediction can be slow when `se_fit = TRUE` and random fields ",
       "are included (i.e., `re_form = NA`). Consider using the `nsim` argument ",
       "to take draws from the joint precision matrix and summarizing the standard ",

R/print.R

@@ -277,8 +277,14 @@ print_other_parameters <- function(x, m = 1L) {
   rho <- get_term_text("rho", "Spatiotemporal AR1 correlation (rho)")
 
   if ("sigma_Z" %in% b$term) {
-    a <- mround(b$estimate[b$term == "sigma_Z"], 2L)
+    # tidy() takes sigma_Z from the sdreport,
+    # which condenses them to unique values if mapped, so:
+    sigma_Z <- x$tmb_obj$report(x$tmb_obj$env$last.par.best)$sigma_Z
+    sigma_Z <- sigma_Z[,m,drop=TRUE]
+    a <- mround(sigma_Z, 2L)
     sigma_Z <- paste0("Spatially varying coefficient SD (", x$spatial_varying,  "): ", a, "\n")
+    sigma_Z <- gsub("\\(\\(", "\\(", sigma_Z) # ((Intercept))
+    sigma_Z <- gsub("\\)\\)", "\\)", sigma_Z) # ((Intercept))
     sigma_Z <- paste(sigma_Z, collapse = "")
   } else {
     sigma_Z <- ""

src/sdmTMB.cpp

@@ -158,6 +158,7 @@ Type objective_function<Type>::operator()()
   DATA_IVECTOR(ar1_fields);
   DATA_IVECTOR(rw_fields);
   DATA_INTEGER(include_spatial);
+  DATA_INTEGER(omit_spatial_intercept);
   DATA_INTEGER(random_walk);
   DATA_INTEGER(ar1_time);
   DATA_IVECTOR(exclude_RE); // DELTA TODO currently shared...
@@ -403,12 +404,14 @@ Type objective_function<Type>::operator()()
       } else {
         Q_temp = Q_s2;
       }
-      PARALLEL_REGION jnll += SCALE(GMRF(Q_temp, s), 1. / exp(ln_tau_O(m)))(omega_s.col(m));
-      if (sim_re(0)) {
-        vector<Type> omega_s_tmp(omega_s.rows());
-        SIMULATE {
-          GMRF(Q_temp, s).simulate(omega_s_tmp);
-          omega_s.col(m) = omega_s_tmp / exp(ln_tau_O(m));
+      if (!omit_spatial_intercept) {
+        PARALLEL_REGION jnll += SCALE(GMRF(Q_temp, s), 1. / exp(ln_tau_O(m)))(omega_s.col(m));
+        if (sim_re(0)) {
+          vector<Type> omega_s_tmp(omega_s.rows());
+          SIMULATE {
+            GMRF(Q_temp, s).simulate(omega_s_tmp);
+            omega_s.col(m) = omega_s_tmp / exp(ln_tau_O(m));
+          }
         }
       }
       if (spatial_covariate) {
@@ -590,7 +593,7 @@ Type objective_function<Type>::operator()()
       for (int t = 0; t < n_t; t++)
         if (!spatial_only(m)) epsilon_st_A.col(m).col(t) =
           A_st * vector<Type>(epsilon_st.col(m).col(t));
-      omega_s_A.col(m) = A_st * vector<Type>(omega_s.col(m));
+      if (!omit_spatial_intercept) omega_s_A.col(m) = A_st * vector<Type>(omega_s.col(m));
       for (int z = 0; z < n_z; z++)
         zeta_s_A.col(m).col(z) = A_st * vector<Type>(zeta_s.col(m).col(z));
     }
@@ -665,7 +668,8 @@ Type objective_function<Type>::operator()()
 
       // Spatially varying effects:
       if (include_spatial) {
-        eta_i(i,m) += omega_s_A(i,m);  // spatial omega
+        if (!omit_spatial_intercept)
+          eta_i(i,m) += omega_s_A(i,m);  // spatial omega
         if (spatial_covariate)
           for (int z = 0; z < n_z; z++)
             eta_i(i,m) += zeta_s_A(i,z,m) * z_i(i,z); // spatially varying covariate DELTA
@@ -1016,7 +1020,7 @@ Type objective_function<Type>::operator()()
         for (int t = 0; t < n_t; t++) {
           proj_epsilon_st_A_unique.col(m).col(t) = proj_mesh * vector<Type>(epsilon_st.col(m).col(t));
         }
-        proj_omega_s_A_unique.col(m) = proj_mesh * vector<Type>(omega_s.col(m));
+        if (!omit_spatial_intercept) proj_omega_s_A_unique.col(m) = proj_mesh * vector<Type>(omega_s.col(m));
       }
 
       // Spatially varying coefficients:

tests/testthat/test-3-spatial-varying.R

@@ -0,0 +1,85 @@
+test_that("Spatially-varying coefficients are estimated correctly for binomial and delta models", {
+  skip_on_cran()
+  skip_on_ci()
+  skip_if_not_installed("INLA")
+  local_edition(2)
+  d <- pcod
+  d$year_scaled <- as.numeric(scale(d$year))
+  mesh10 <- make_mesh(d, c("X", "Y"), cutoff = 10)
+  m1 <- sdmTMB(
+    data = d,
+    formula = present ~ 1 + year_scaled,
+    spatial_varying = ~ 0 + year_scaled,
+    mesh = mesh10,
+    family = binomial()
+  )
+  # m1.1 <- sdmTMB(
+  #   data = d,
+  #   formula = present ~ 1 + year_scaled,
+  #   spatial_varying = ~ 1 + year_scaled, #<
+  #   spatial = "off", #<
+  #   mesh = mesh10,
+  #   family = binomial()
+  # )
+  # expect_equal(m1$model$objective, m1.1$model$objective)
+
+  b1 <- tidy(m1, effects = "ran_pars", conf.int = TRUE)
+  expect_equal(b1$estimate[3], 0.312, tolerance = 0.1)
+
+  m1.2 <- sdmTMB(
+    data = d,
+    formula = present ~ 1 + year_scaled,
+    spatial_varying = ~ 1 + year_scaled,
+    spatial = "on",
+    mesh = mesh10,
+    family = binomial()
+  )
+
+  expect_equal(m1$model$objective, m1.2$model$objective)
+
+  # warn: probably don't want to do this!
+  expect_message({
+    m1.3 <- sdmTMB(
+      data = d,
+      formula = present ~ 1 + year_scaled,
+      spatial_varying = ~ 0 + as.factor(year), #<
+      spatial = "on", #<
+      mesh = mesh10,
+      family = binomial(),
+      do_fit = FALSE
+    )
+  }, regexp = "intercept")
+
+  # better:
+  m1.4 <- sdmTMB(
+    data = d,
+    formula = present ~ 1 + year_scaled,
+    spatial_varying = ~ 0 + as.factor(year),
+    spatial = "off",
+    mesh = mesh10,
+    family = binomial()
+  )
+  m1.4
+
+  m1.5 <- sdmTMB(
+    data = d,
+    formula = present ~ 1 + year_scaled,
+    spatial_varying = ~ 1 + as.factor(year),
+    spatial = "on",
+    mesh = mesh10,
+    family = binomial()
+  )
+  m1.5
+  p <- predict(m1.5, newdata = d)
+
+  # also check that binomial portion of delta model matches the above
+  m2 <- sdmTMB(
+    data = d,
+    formula = density ~ 1 + year_scaled,
+    spatial_varying = ~ 0 + year_scaled,
+    mesh = mesh10,
+    family = delta_gamma()
+  )
+  b2 <- tidy(m2, effects = "ran_pars", conf.int = TRUE)
+  expect_equal(b2$estimate[3], b1$estimate[3], tolerance = 1e-3)
+})