spatially varying coefficient uncertainty

[lysoifer]

Hi,

I am modelling proportion data with a binomial distribution using the observed proportions as the response and the number of observations (N) as the weights. The model includes a spatially varying coefficient as well as a spatial random field. I am trying to estimate the uncertainty of the spatially varying coefficient based on the owls example in Anderson et al. (2022) sdmTMB: an R package for fast, flexible, and user-friendly generalized linear mixed effects models with spatial and spatiotemporal random fields. However, when I run spread_sims(fit, nsim = n), I receive the following error: "Error in am.fit(G, control = control, gamma = gamma, ...) : y must be univariate unless binomial". Since the model is binomial, I am not sure how to solve this. Is it possible to use run spread_sims for a binomial model to estimate the SVC uncertainty?

Thanks,
Lydia

[seananderson]

Hi Lydia, can you provide a reproducible example?

I tried to reproduce the issue, but this seems to work without a problem:

library(sdmTMB)

# simulate some data:
set.seed(123)
predictor_dat <- data.frame(
  X = runif(300), Y = runif(300),
  a1 = rnorm(300), year = rep(-2:3, each = 50)
)
mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
sim_dat <- sdmTMB_simulate(
  formula = ~ 1,
  data = predictor_dat,
  time = "year",
  mesh = mesh,
  family = gaussian(),
  range = 0.5,
  phi = 0.1,
  sigma_O = 0.2,
  seed = 42,
  B = 0
)
sim_dat$mu <- (sim_dat$omega_s + 0.2) * sim_dat$year
sim_dat$size <- rpois(nrow(sim_dat), lambda = 5) + 1
sim_dat$observed <- rbinom(nrow(sim_dat), size = sim_dat$size, prob = plogis(sim_dat$mu))
sim_dat$observed_prop <- sim_dat$observed / sim_dat$size

# fit SVC model:
fit <- sdmTMB(
  observed_prop ~ 1 + year, 
  spatial = "off", 
  spatiotemporal = "off", 
  time = "year", 
  mesh = mesh, 
  spatial_varying = ~ 0 + year, 
  data = sim_dat, 
  family = binomial(), 
  weights = sim_dat$size
)

s <- spread_sims(fit)
head(s)
#>   .iteration X.Intercept.      year   ln_tau_Z     range
#> 1          1  -0.02623886 0.4164971 -2.4003479 0.2061976
#> 2          2  -0.11554581 0.2832342 -0.8288195 0.7894908
#> 3          3  -0.10153308 0.2936959 -2.2770656 0.3015632
#> 4          4  -0.06927436 0.2752159 -2.0803629 0.2422635
#> 5          5   0.08345718 0.2833141 -2.1405401 0.2118434
#> 6          6  -0.06760870 0.1898790 -0.6503099 0.5226862

hist(s$year)

hist(s$range)

^{Created on 2024-09-10 with reprex v2.1.1}