Skip to contents

If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under ‘Articles’.

Using the built-in British Columbia Queen Charlotte Sound Pacific Cod dataset, we might be interested in fitting a model that describes spatially varying trends through time. The data are as follows:

  • There are columns for depth and depth squared.
  • Depth was centred and scaled by its standard deviation and we’ve included those in the data frame so that they could be used to similarly scale the prediction grid.
  • The density units should be kg/km2.
  • Here, X and Y are coordinates in UTM zone 9.

We will set up our SPDE mesh with a relatively coarse resolution so that this vignette builds quickly:

pcod_spde <- make_mesh(pcod, c("X", "Y"), cutoff = 12)

We will fit a model that includes a slope for ‘year’, an intercept spatial random field, and another random field for spatially varying slopes the represent trends over time in space (spatial_varying argument). Our model just estimates an intercept and accounts for all other variation through the random effects.

First, we will set up a column for time that is Normal(0, 1) to help with estimation:

d <- pcod
d$scaled_year <- (pcod$year - mean(pcod$year)) / sd(pcod$year)

Now fit a model using spatial_varying ~ 0 + scaled_year:

(The 0 + drops the intercept, although sdmTMB would take care of that anyways here.)

m1 <- sdmTMB(density ~ scaled_year, data = d,
  mesh = pcod_spde, family = tweedie(link = "log"),
  spatial_varying = ~ 0 + scaled_year, time = "year",
  spatiotemporal = "off")

We have turned off spatiotemporal random fields for this example for simplicity, but they also could be IID or AR1.

Let’s extract some parameter estimates. Look for sigma_Z:

tidy(m1, = TRUE)
#> # A tibble: 2 × 5
#>   term        estimate std.error conf.low conf.high
#>   <chr>          <dbl>     <dbl>    <dbl>     <dbl>
#> 1 (Intercept)    2.85      0.346    2.17      3.52 
#> 2 scaled_year   -0.126     0.117   -0.356     0.104
tidy(m1, "ran_pars", = TRUE)
#> Standard errors intentionally omitted because they have been calculated in log
#> space.
#> # A tibble: 5 × 5
#>   term      estimate std.error conf.low conf.high
#>   <chr>        <dbl> <lgl>        <dbl>     <dbl>
#> 1 range       26.2   NA          18.2      37.8  
#> 2 phi         14.1   NA          13.3      15.0  
#> 3 sigma_O      2.13  NA           1.75      2.60 
#> 4 sigma_Z      0.625 NA           0.459     0.850
#> 5 tweedie_p    1.59  NA           1.57      1.61

Let’s look at the predictions and estimates of the spatially varying coefficients on a grid:

plot_map_raster <- function(dat, column = est) {
  ggplot(dat, aes(X, Y, fill = {{ column }})) +
    geom_raster() +
    facet_wrap(~year) +
    coord_fixed() +

First, we need to predict on a grid. We also need to add a column for scaled_year to match the fitting:

nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
nd$scaled_year <- (nd$year - mean(pcod$year)) / sd(pcod$year)
p1 <- predict(m1, newdata = nd)

First let’s look at the spatial trends.

We will just pick out a single year to plot since they should all be the same for the slopes. Note that these are in log space. zeta_s are the spatially varying coefficients.

plot_map_raster(filter(p1, year == 2003), zeta_s_scaled_year)

This is the spatially varying intercept:

plot_map_raster(filter(p1, year == 2003), omega_s) + scale_fill_gradient2()

These are the predictions including all fixed and random effects plotted in log space.

plot_map_raster(filter(p1, year == 2003), est)

And we can look at just the spatiotemporal random effects for models 2 and 3 (intercept + slope combined):

plot_map_raster(filter(p1, year == 2003), est_rf)