Fix bug from last commit #465

Author: seananderson

R/project.R

@@ -343,7 +343,7 @@ insert_pars <- function(par, nm, .n, delta_model = FALSE) {
   ret
 }
 
-move_proj_to_tmbdat <- function(x, object, newdata, called_by_simulate = FALSE) {
+move_proj_to_tmbdat <- function(x, object, newdata, called_by_simulate = FALSE, size = NULL) {
   x$do_predict <- 0L
   x$year_i <- x$proj_year
   ## x$A_st <- x$proj_mesh
@@ -364,13 +364,30 @@ move_proj_to_tmbdat <- function(x, object, newdata, called_by_simulate = FALSE)
   x$y_i <- matrix(NA, ncol = n_m, nrow = nrow(x$proj_X_ij[[1]])) # fake
   x$weights_i <- rep(1, nrow(x$y_i)) # fake: FIXME: bring in?
   x$area_i <- rep(1, nrow(x$y_i)) # fake FIXME: bring in for index??
+
   if (!called_by_simulate) {
     unique_size <- unique(x$size)
     if (length(unique_size) != 1L) {
       cli_abort("This function hasn't been set up to work with binomial size specified yet.")
     }
     x$size <- rep(1, nrow(x$y_i)) # FIXME: bring in?
   }
+  if (!is.null(size)) {
+    if (length(size) != nrow(x$y_i)) cli_abort("size argument doesn't match rows of newdata")
+    x$size <- size
+  }
+  if (!all(x$size == 1L)) {
+    if (is.null(size)) {
+      msg <- c(
+        "It looks like the size/trials were specified for a binomial family with the weights argument.",
+        "When simulating with newdata all sizes/trials will be assumed to be 1 unless you specify a vector via the `size` argument."
+      )
+      cli_inform(msg)
+      x$size <- rep(1, nrow(x$y_i))
+    }
+  } else { # all 1s, but newdata; expand as needed silently
+    x$size <- rep(1, nrow(x$y_i))
+  }
   x$do_predict <- 0L
 
   # nullify large data objects that are no longer needed:

R/tmb-sim.R

@@ -345,6 +345,9 @@ sdmTMB_simulate <- function(formula,
 #'   lets you parse out whatever elements you want from the simulation.
 #'   Not usually needed.
 #' @param observation_error Logical. Simulate observation error?
+#' @param size A vector of size (trials) in the case of a binomial family with
+#'   `newdata` specified. If left `NULL` with `newdata`, will be assumed to
+#'   be a vector of 1s.
 #' @param silent Logical. Silent?
 #' @param ... Extra arguments passed to [predict.sdmTMB()]. E.g., one may wish
 #'   to pass an `offset` argument if `newdata` are supplied in a model with an
@@ -397,6 +400,7 @@ simulate.sdmTMB <- function(object, nsim = 1L, seed = sample.int(1e6, 1L),
                             mcmc_samples = NULL,
                             return_tmb_report = FALSE,
                             observation_error = TRUE,
+                            size = NULL,
                             silent = FALSE,
                             ...) {
   set.seed(seed)
@@ -426,7 +430,7 @@ simulate.sdmTMB <- function(object, nsim = 1L, seed = sample.int(1e6, 1L),
     # generate prediction TMB data list
     p <- predict(object, newdata = newdata, return_tmb_data = TRUE, ...)
     # move data elements over
-    p <- move_proj_to_tmbdat(p, object, newdata, called_by_simulate = TRUE)
+    p <- move_proj_to_tmbdat(p, object, newdata, called_by_simulate = TRUE, size = size)
     p$sim_re <- tmb_dat$sim_re
     tmb_dat <- p
   }

man/simulate.sdmTMB.Rd

@@ -195,10 +195,10 @@ test_that("simulate.sdmTMB works with sizes in binomial GLMs #465", {
   fit <- sdmTMB(prop ~ x, data = dat, weights = w, family = binomial(), spatial = "off")
   dat$X <- dat$Y <- NA # FIXME
   set.seed(1)
-  snd <- simulate(fit, nsim = 500, newdata = dat)
+  snd <- simulate(fit, nsim = 500, newdata = dat, size = w)
   set.seed(1)
-  s <- simulate(fit, nsim = 500, newdata = dat)
-  expect_identical(s, snd)
+  s <- simulate(fit, nsim = 500)
+  expect_equal(as.matrix(s), as.matrix(snd), ignore_attr = TRUE)
   expect_true(nrow(snd) == 200L)
   expect_true(ncol(snd) == 500L)
   sim_means <- rowMeans(snd)
@@ -265,4 +265,25 @@ test_that("simulate without observation error works for binomial likelihoods #43
   m <- apply(s.b, 1, mean)
   p <- predict(fit.b, newdata = qcs_grid)
   expect_gt(cor(plogis(p$est), m), 0.95)
+
+  # with size specified
+  expect_error({simulate(
+    fit.b,
+    newdata = qcs_grid,
+    nsim = 1,
+    observation_error = FALSE,
+    size = c(1, 2, 3)
+  )}, regexp = "size")
+
+  s.b1 <- simulate(
+    fit.b,
+    newdata = qcs_grid,
+    type = "mle-mvn", # fixed effects at MLE values and random effect MVN draws
+    mle_mvn_samples = "multiple", # take an MVN draw for each sample
+    nsim = 50, # increase this for more stable results
+    observation_error = FALSE, # do not include observation error
+    seed = 23859,
+    size = sample(1:9, size = nrow(qcs_grid), replace = TRUE)
+  )
+  expect_false(identical(s.b, s.b1))
 })