`vignettes/spatial-trend-models.Rmd`

`spatial-trend-models.Rmd`

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/km
^{2}. - 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:

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:

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, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 (Intercept) 3.4251463 0.3897312 2.6612872 4.18900536
#> 2 scaled_year -0.2140606 0.1164730 -0.4423436 0.01422232
tidy(m1, "ran_pars", conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 range 29.7381725 NA 20.775802 42.5667765
#> 3 phi 14.0667396 NA 13.226193 14.9607048
#> 4 sigma_O 1.9167968 NA 1.489051 2.4674171
#> 5 sigma_Z 0.5442651 NA 0.398302 0.7437182
#> 6 tweedie_p 1.5952763 NA 1.575281 1.6149607
```

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_string("X", "Y", fill = column)) +
geom_raster() +
facet_wrap(~year) +
coord_fixed() +
scale_fill_viridis_c()
}
```

First, we need to predict on a grid. We also need to add a column for `scaled_year`

to match the fitting:

```
nd <- qcs_grid
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")`