2024-02-21

Problem

Spatially explicit models

  • Spatially explicit models are often not spatially explicit
  • We often just transform it to tabular data
  • We have to go out of our wa

Spatially explicit GLMM

data("pcod")
bnd <- INLA::inla.nonconvex.hull(cbind(pcod$X, pcod$Y), convex = -0.1)
mesh_inla <- INLA::inla.mesh.2d(boundary = bnd, max.edge = c(25, 50))
mesh <- make_mesh(pcod, c("X", "Y"), mesh = mesh_inla)
plot(mesh)

First model

m <- sdmTMB(
  data = pcod,
  formula = present ~ depth_scaled + + I(depth_scaled^2),
  mesh = mesh, # can be omitted for a non-spatial model
  family = binomial(link = "logit"),
  spatial = "off"
)

Spatial residuals

coords <- pcod[, c("X", "Y")]
coords$Resid <- residuals(m)

Spatial on

m1 <- sdmTMB(data = pcod, formula = present ~ depth_scaled + +I(depth_scaled^2),
    mesh = mesh, family = binomial(link = "logit"), spatial = "on")

Spatial residuals

coords$ResidSpatial <- residuals(m1)

Lets look at the parameters

t <- tidy(m, conf.int = T) |>
    mutate(Spatial = "off")
t1 <- tidy(m1, conf.int = T) |>
    mutate(Spatial = "on")

tidies <- dplyr::bind_rows(t, t1)
term estimate std.error conf.low conf.high Spatial
(Intercept) 0.5659870 0.0597859 0.4488087 0.6831653 off
depth_scaled -1.0359044 0.0726568 -1.1783092 -0.8934996 off
I(depth_scaled^2) -0.9925942 0.0606616 -1.1114888 -0.8736996 off
(Intercept) 1.2340958 0.4845461 0.2844028 2.1837888 on
depth_scaled -2.2558559 0.2110190 -2.6694455 -1.8422664 on
I(depth_scaled^2) -1.6291480 0.1313962 -1.8866798 -1.3716162 on

Add Spatiotemporal random fields

m2 <- sdmTMB(
  data = pcod,
  formula = present ~ depth_scaled + I(depth_scaled^2),
  mesh = mesh,
  family = binomial(link = "logit"),
  spatial = "on",
  time = "year",
  # first-order autoregressive 
  spatiotemporal = "AR1",
  # Show that there are missing years
  extra_time = c(2006, 2008, 2010, 2012, 2014, 2016)
)

Check for issues

sanity(m2)
## ✔ Non-linear minimizer suggests successful convergence
## ✔ Hessian matrix is positive definite
## ✔ No extreme or very small eigenvalues detected
## ✔ No gradients with respect to fixed effects are >= 0.001
## ✔ No fixed-effect standard errors are NA
## ✔ No standard errors look unreasonably large
## ✔ No sigma parameters are < 0.01
## ✔ No sigma parameters are > 100
## ✔ Range parameter doesn't look unreasonably large

Calculate response

NewDf <- data.frame(depth_scaled = seq(from = min(pcod$depth_scaled), to = max(pcod$depth_scaled),
    length.out = 100), year = 2015)

p <- predict(m2, newdata = NewDf, se_fit = TRUE, re_form = NA, type = "response")

p <- p |>
    mutate(prob = 1/(1 + exp(-est)), upper = 1/(1 + exp(-(est + 1.96 * est_se))),
        lower = 1/(1 + exp(-(est - 1.96 * est_se))))

Plot response

Time varying response

m3 <- sdmTMB(
    data = pcod,
    formula = present ~ 0 + as.factor(year),
    time_varying = ~ 0 + depth_scaled + I(depth_scaled^2),
    mesh = mesh,
    family = binomial(link = "logit"),
    spatial = "on",
    time = "year",
    # first-order autoregressive 
    spatiotemporal = "AR1",
    # Show that there are missing years
    extra_time = c(2006, 2008, 2010, 2012, 2014, 2016)
)

Prediction time varying

NewDf <- expand.grid(depth_scaled = seq(from = min(pcod$depth_scaled), to = max(pcod$depth_scaled),
    length.out = 50), year = unique(pcod$year))

p <- predict(m3, newdata = NewDf, se_fit = TRUE, re_form = NA, type = "response")

p <- p |>
    mutate(prob = 1/(1 + exp(-est)), upper = 1/(1 + exp(-(est + 1.96 * est_se))),
        lower = 1/(1 + exp(-(est - 1.96 * est_se))))

Plot time varying response

Quick model comparison

Sel <- MuMIn::model.sel(list(m, m1, m2, m3)) |>
    as.data.frame() |>
    dplyr::select("spatial", "time", "spatiotemporal", "time_varying", "df", "AICc",
        "delta", "weight")
spatial time spatiotemporal time_varying df AICc delta weight
4 on year AR1 Yes 21 2036.88 0.00 7.411335e-01
3 on year AR1 No 7 2038.98 2.10 2.588665e-01
2 on No 5 2078.71 41.83 6.105728e-10
1 off No 3 2392.08 355.21 5.468264e-78