We’ll repeat the same models used for the index standardization vignette, using the built-in data for Pacific cod.

As a summary, - We’ve included 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.

glimpse(pcod)
#> Rows: 2,143
#> Columns: 12
#> $ year          <int> 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2…
#> $ X             <dbl> 446.4752, 446.4594, 448.5987, 436.9157, 420.6101, 417.7…
#> $ Y             <dbl> 5793.426, 5800.136, 5801.687, 5802.305, 5771.055, 5772.…
#> $ depth         <dbl> 201, 212, 220, 197, 256, 293, 410, 387, 285, 270, 381, …
#> $ density       <dbl> 113.138476, 41.704922, 0.000000, 15.706138, 0.000000, 0…
#> $ present       <dbl> 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ lat           <dbl> 52.28858, 52.34890, 52.36305, 52.36738, 52.08437, 52.09…
#> $ lon           <dbl> -129.7847, -129.7860, -129.7549, -129.9265, -130.1586, …
#> $ depth_mean    <dbl> 5.155194, 5.155194, 5.155194, 5.155194, 5.155194, 5.155…
#> $ depth_sd      <dbl> 0.4448783, 0.4448783, 0.4448783, 0.4448783, 0.4448783, …
#> $ depth_scaled  <dbl> 0.3329252, 0.4526914, 0.5359529, 0.2877417, 0.8766077, …
#> $ depth_scaled2 <dbl> 0.11083919, 0.20492947, 0.28724555, 0.08279527, 0.76844…

As before, the SPDE mesh is created using 75 knots for a balance between speed and accuracy. You will likely want to use more for applied scenarios. You will want to make sure that increasing the number of knots does not change the conclusions.

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

We will extend the simple GLMM used in the index standardization vignette to include a threshold effect. For real applications, the threshold effect might be a function of habitat covariates, environmental variables, etc. Because ‘depth’ is the only external variable in the pcod data frame, we’ll use that. For any threshold model, you can specify the model as a linear threshold (breakpt()) or logistic (logistic()) as part of the formula. Note: as before, if we want to use this model for index standardization then we need to include 0 + as.factor(year) or -1 + as.factor(year) so that we have a factor predictor that represents the mean estimate for each time slice.

This first example uses the ‘depth_scaled’ covariate (depth, standardized ~ N(0,1)) and the logistic function, similar to selectivity curves used in fisheries models. The form is

\[s(x)=\tau + \psi *\quad { \left[ 1+{ \exp }^{ -ln\left( 19 \right) \quad *\quad \left( x-s50 \right) \quad /\quad \left( s95\quad -\quad s50 \right) } \right] }^{ -1 }\] where \(\psi\) is a scaling parameter (controlling the height of the y-axis for the response, and is unconstrained), \(\tau\) is an intercept, \(s50\) is a parameter controlling the point at which the function reaches 50% of the maximum (\(\psi\)), and \(s95\) is a parameter controlling the point at which the function reaches 95%. The parameter \(s50\) is unconstrained, and \(s95\) is constrained to be larger than \(s50\).

m <- sdmTMB(
  data = pcod,
  formula = density ~ 0 + logistic(depth_scaled) + as.factor(year),
  time = "year", spde = pcod_spde,
  family = tweedie(link = "log")
)

We can then look at the important coefficients from the model. Here we’re looking at ‘s50’, ‘s95’, and ‘smax’.

print(m)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: density ~ 0 + logistic(depth_scaled) + as.factor(year)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: pcod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> as.factor(year)2003    -0.52    0.63
#> as.factor(year)2004     0.06    0.62
#> as.factor(year)2005    -0.08    0.62
#> as.factor(year)2007    -1.28    0.63
#> as.factor(year)2009    -0.80    0.63
#> as.factor(year)2011    -0.33    0.62
#> as.factor(year)2013    -0.13    0.61
#> as.factor(year)2015    -0.16    0.62
#> as.factor(year)2017    -0.75    0.63
#> depth_scaled-s50       -1.48    0.12
#> depth_scaled-s95       -0.53    0.22
#> depth_scaled-smax       3.22    0.53
#> 
#> Matern range parameter: 22.09
#> Dispersion parameter: 12.33
#> Spatial SD (sigma_O): 2.02
#> Spatiotemporal SD (sigma_E): 1.40
#> ML criterion at convergence: 6435.001

We can visualize the effect:

nd <- data.frame(
  depth_scaled = seq(min(pcod$depth_scaled) + 0.2, 
    max(pcod$depth_scaled) - 0.2, length.out = 100), 
  year = 2015L # a chosen year
)
p <- predict(m, newdata = nd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
#> Warning: argument `xy_cols` is deprecated; this information is already in the
#> output of the `make_mesh(). Did you use `make_spde()`?
ggplot(p, aes(depth_scaled, exp(est), 
  ymin = exp(est - 1.96 * est_se), 
  ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4)

Next, we can try to fit the same model but use a linear breakpoint / cutpoint model:

m <- sdmTMB(
  data = pcod,
  formula = density ~ 0 + breakpt(depth_scaled) + as.factor(year),
  time = "year", spde = pcod_spde,
  family = tweedie(link = "log")
)

For this model, the important parameters are the slope, ‘depth_scaled-slope’, and breakpoint, ‘depth_scaled-breakpt’.

print(m)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: density ~ 0 + breakpt(depth_scaled) + as.factor(year)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: pcod
#> Family: tweedie(link = 'log')
#>                      coef.est coef.se
#> as.factor(year)2003      3.63    0.33
#> as.factor(year)2004      4.19    0.32
#> as.factor(year)2005      4.00    0.32
#> as.factor(year)2007      2.89    0.33
#> as.factor(year)2009      3.38    0.33
#> as.factor(year)2011      3.78    0.32
#> as.factor(year)2013      4.09    0.32
#> as.factor(year)2015      3.94    0.32
#> as.factor(year)2017      3.39    0.33
#> depth_scaled-slope       1.12    0.17
#> depth_scaled-breakpt    -1.59    0.23
#> 
#> Matern range parameter: 15.93
#> Dispersion parameter: 12.14
#> Spatial SD (sigma_O): 1.82
#> Spatiotemporal SD (sigma_E): 1.70
#> ML criterion at convergence: 6397.346

We can visualize the linear breakpoint effect:

p <- predict(m, newdata = nd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
ggplot(p, aes(depth_scaled, exp(est), 
  ymin = exp(est - 1.96 * est_se), 
  ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4)