+ - 0:00:00
Notes for current slide
Notes for next slide

Example: Spatiotemporal modeling of Pacific cod

DFO TESA sdmTMB workshop

January 16–19 2023

1 / 20

Extending to a spatiotemporal model

Reminder, the Pacific Cod dataset:

pcod %>% 
  select(year, X, Y, depth, density, present) %>% 
  head()
#> # A tibble: 6 × 6
#>    year     X     Y depth density present
#>   <int> <dbl> <dbl> <dbl>   <dbl>   <dbl>
#> 1  2003  446. 5793.   201   113.        1
#> 2  2003  446. 5800.   212    41.7       1
#> 3  2003  449. 5802.   220     0         0
#> 4  2003  437. 5802.   197    15.7       1
#> 5  2003  421. 5771.   256     0         0
#> 6  2003  418. 5772.   293     0         0
mesh <- make_mesh(pcod, xy_cols = c("X", "Y"), cutoff = 10)
2 / 20

Our previous spatial model fit

fit <- sdmTMB(
  present ~ poly(log(depth), 2),
  data = pcod,
  mesh = mesh,
  family = binomial(link = "logit"),
  spatial = "on"
)
3 / 20

Switching to density

fit <- sdmTMB(
  density ~ poly(log(depth), 2),
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "on"
)
4 / 20

Switching to a smoother + annual mean

fit <- sdmTMB(
  density ~ s(depth) + 0 + as.factor(year),
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "on"
)
5 / 20

Adding spatiotemporal fields

fit <- sdmTMB(
  density ~ s(depth, k = 5) + 0 + as.factor(year),
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "on",
  time = "year",
  spatiotemporal = "iid",
  silent = FALSE # show progress!
)
6 / 20

Inspecting model convergence

sanity(fit)
#> ✔ 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
7 / 20

Inspecting the model fit

fit
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: density ~ s(depth) + 0 + as.factor(year)
#> Time column: "year"
#> Mesh: mesh
#> Data: pcod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> as.factor(year)2003     1.79    0.29
#> as.factor(year)2004     2.36    0.27
#> ...
#> 
#> Dispersion parameter: 10.78
#> Tweedie p: 1.49
#> Matern range: 12.85
#> Spatial SD: 1.81
#> Spatiotemporal SD: 1.78
#> ML criterion at convergence: 6246.433
8 / 20

Model residuals

Warning: these residuals are fast but might look off even if the model is fine. Also see MCMC residuals. See the 'Residual checking' vignette and the example in the exercises.

set.seed(1)
rq_res <- residuals(fit) # randomized quantile residuals
qqnorm(rq_res);qqline(rq_res)

9 / 20

Model residuals in space

pcod$resids <- residuals(fit)
filter(pcod, year %in% c(2015, 2017)) %>% 
  ggplot(aes(X, Y, colour = resids)) +
  geom_point() +
  facet_wrap(~year) +
  scale_colour_gradient2() +
  coord_fixed()

10 / 20

Predicting on the survey grid

p <- predict(fit, newdata = qcs_grid)
#> Free parallelADFun object.
  • est: Overall estimate in link space (log here)
  • est_non_rf: Estimate of non-random-field components
  • est_rf: Estimate of random-field components
  • omega_s: Spatial random field
  • epsilon_st: Spatiotemporal random field
11 / 20

Plotting overall predictions

ggplot(p, aes(X, Y, fill = exp(est))) +
  geom_raster() +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~year) +
  coord_fixed()

12 / 20

Plotting overall predictions (truncated)

max_est <- quantile(p$est, 0.999)
p <- mutate(p, est_trim = if_else(est > max_est, max_est, est))

13 / 20

Plotting main effect contributions

ggplot(p, aes(X, Y, fill = exp(est_non_rf))) +
  geom_raster() +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~year) +
  coord_fixed()

14 / 20

Plotting spatial random effects

ggplot(p, aes(X, Y, fill = omega_s)) +
  geom_raster() +
  scale_fill_gradient2() +
  facet_wrap(~year) +
  coord_fixed()

15 / 20

Plotting spatiotemporal random effects

ggplot(p, aes(X, Y, fill = epsilon_st)) +
  geom_raster() +
  scale_fill_gradient2() +
  facet_wrap(~year) +
  coord_fixed()

16 / 20

Visualizing spatial uncertainty

By sampling from the joint precision matrix

psims <- predict(fit, newdata = qcs_grid, nsim = 500)
dim(psims)
#> [1] 65826   500
psims[1:3, 1:4]
#>           [,1]       [,2]       [,3]      [,4]
#> 2003 -2.953744 -0.6653106 -4.2355453 -1.764134
#> 2003  2.601457  3.2971739  0.7631676  2.999720
#> 2003  3.584363  4.1505816  1.7911469  4.206469
p$sd <- apply(psims, 1, sd)
17 / 20

Visualizing spatial uncertainty

filter(p, year %in% c(2013, 2015, 2017)) %>% 
  ggplot(aes(X, Y, fill = sd)) +
  geom_raster() +
  facet_wrap(~year) +
  scale_fill_viridis_c(trans = "log10") +
  coord_fixed()

18 / 20

Considering anisotropy

  • Default is isotropic: correlation decays in all directions at same rate
  • Anisotropic: directionally dependent spatial correlation
fit_aniso <- update(fit, anisotropy = TRUE)
AIC(fit, fit_aniso)
#>           df      AIC
#> fit       16 12524.87
#> fit_aniso 18 12528.46
  • Not favoured here; often important on narrow continental shelves
19 / 20

Plotting anisotropy

Plot shows range in all directions from zero in the center

plot_anisotropy(fit_aniso)

20 / 20

Extending to a spatiotemporal model

Reminder, the Pacific Cod dataset:

pcod %>% 
  select(year, X, Y, depth, density, present) %>% 
  head()
#> # A tibble: 6 × 6
#>    year     X     Y depth density present
#>   <int> <dbl> <dbl> <dbl>   <dbl>   <dbl>
#> 1  2003  446. 5793.   201   113.        1
#> 2  2003  446. 5800.   212    41.7       1
#> 3  2003  449. 5802.   220     0         0
#> 4  2003  437. 5802.   197    15.7       1
#> 5  2003  421. 5771.   256     0         0
#> 6  2003  418. 5772.   293     0         0
mesh <- make_mesh(pcod, xy_cols = c("X", "Y"), cutoff = 10)
2 / 20
Paused

Help

Keyboard shortcuts

↑, ←, Pg Up, k Go to previous slide
↓, →, Pg Dn, Space, j Go to next slide
Home Go to first slide
End Go to last slide
Number + Return Go to specific slide
b / m / f Toggle blackout / mirrored / fullscreen mode
c Clone slideshow
p Toggle presenter mode
t Restart the presentation timer
?, h Toggle this help
Esc Back to slideshow