simulate() fails for binomial spatiotemporal model with weights

[akshayanand93]

I have fit an AR(1) spatiotemporal model with the following specification:

model <- smdTMB(y ~ poly(x1, 2) + poly(x2, 2) + poly(x3, 2),
data = dat_train,
weights = dat_train$weights,
extra_time = 18:19,
mesh = mesh,
spatial = "off",
spatiotemporal = "ar1",
family = binomial(link = "cloglog"),
time = "time_step")

I am now trying to use the simulate() function to get epsilon corrected estimates like this:

sims <- simulate(model,
newdata = dat,
type = "mle-mvn",
mle_mvn_samples = "multiple",
nsim = 500,
observation_error = FALSE)

My "dat" here contains the training data that I used to fit the model along with some holdout data with the extra time steps specified during model fitting. However, the simulate() function returns the error "This function hasn't been set up to work with binomial size specified yet."

[seananderson]

I can bypass that error here and pass the size vector through in simulate() (it's there to block some untested parts of the project() function that uses a shared helper function). However, I wanted to check that you're using weights to represent the number of trials ('size') for the binomial likelihood here. sdmTMB following glm() in using weights for size and it turns out that creates all sorts of pain points. If you're intending to use weights as 1s and 0s to drop data for validation, that works for all likelihoods except binomial.

[seananderson]

OK, try the main branch now. Changes and a unit test are shown here: ea08d86

This should hopefully be on CRAN soon.

[akshayanand93]

Thank you! I will give it a try and let you know. To answer your question, I did used the weights to represent the binomial trials in these models.

[seananderson]

Sorry, package checking showed this 'fix' made some other configurations crash when simulating. I'm looking into it.

[seananderson]

That was more complicated than I thought. Should be fixed now. There's a new size argument in simulate.sdmTMB now that can optionally be used to pass in a size vector when newdata are supplied with the binomial. This shouldn't actually do anything in the case of observation_error = FALSE, because size only affects the sampling from rbinom(), which isn't getting used. If not supplied with newdata, the size vector will be assumed to be a vector of 1s.

[akshayanand93]

Thank you very much! Just to clarify, I can use the size argument to pass in my size vector even if observation_error = FALSE is specified?

[seananderson]

It gets passed through but should currently do nothing if observation_error = FALSE. If observation_error = FALSE, then the mean should be returned but rbinom(1, size = size, prob = mu) never gets applied, so what you get back is mu. I'm not sure if this should be a special case in the code or if the user should just apply the size multiplier vector afterwords. If you or anyone is aware of a precedent for how this interface should work, let me know. I can see the argument that mu * size is most consistent with what simulation with observation error returns given a supplied size argument.

[seananderson]

I convinced myself that we should be multiplying mu by size if size is supplied with new data. That should now be happening and is tested in this unit test:

sdmTMB/tests/testthat/test-tmb-simulation.R

Lines 280 to 293 in 7ae4cd0

	w <- sample(1:9, size = nrow(qcs_grid), replace = TRUE)
	s.b1 <- simulate(
	fit.b,
	newdata = qcs_grid,
	type = "mle-mvn",
	mle_mvn_samples = "multiple",
	nsim = 50,
	observation_error = FALSE,
	seed = 23859,
	size = w
	)
	expect_true(max(s.b1) > 1)
	expect_equal(mean(s.b1[1,]), m[1] * w[1])
	expect_equal(mean(s.b1[51,]), m[51] * w[51])

Code change here: 9ba405e

[akshayanand93]

Thank you very much for fixing this! I have tested this with my data and it works well.