Using sdmTMB for spatiotemporal modelling
Understanding spatial and spatiotemporal models:
- Understand the spatial range.
- Become familiar with spatial random fields through simulation.
- Become familiar with spatiotemporal random fields and how they can combine with spatial random fields.
- Understand how correlation in spatiotemporal random fields looks.
Understanding the spatial range
The spatial range is the distance at which two data points are effectively independent. Technically, this is the distance at which correlation declines to approximately 0.13.
We can visualize the Matérn correlation function for various range values:
range <- 0.5 #< try changing this
kappa <- sqrt(8) / range
distance <- seq(0.001, 1, length.out = 100)
correlation <- kappa * distance * besselK(kappa * distance, nu = 1)
plot(distance, correlation, type = "l", ylim = c(0, 1))
abline(h = 0.13, lty = 2) # 0.13 correlation
abline(v = range, lty = 2)
Here we have a unitless distance from 0 to 1. In real situations these would be in something like m or km and so the range would be units of m or km.
Exercise:
- Try running the previous chunk over and over with various
rangevalues. Note how the correlation function changes.
Spatial simulation
Next we will simulate some spatial random fields to understand what they look like and how their parameters control them.
First, this chunk creates a set of X and Y spatial coordinates and constructs a mesh using a minimum triangle edge length cutoff distance of 0.05. We need this for the internal random field calculations with the SPDE.
# A grid of X and Y values:
predictor_dat <- expand.grid(
X = seq(0, 1, length.out = 100), Y = seq(0, 1, length.out = 100)
)
mesh <- make_mesh(predictor_dat, c("X", "Y"), cutoff = 0.05)
plot(mesh)
The basic function for simulating random fields in sdmTMB is sdmTMB_simulate().
The formula here contains only an intercept but could include covariates.
The family is Gaussian to specify that we simulate a normally distributed response.
The range (see above) is 0.6.
The two important standard deviation parameters are
sigma_O(controlling the marginal standard deviation ofomega, the spatial random field) andphi(controlling the observation error standard deviation).We are plotting
observed(observed values: the expected value plus observation error).
sim_dat <- sdmTMB_simulate(
formula = ~ 1,
data = predictor_dat,
mesh = mesh,
family = gaussian(link = "identity"),
range = 0.6, #< try changing this (spatial range)
sigma_O = 0.2, #< try changing this (spatial standard deviation)
phi = 0.001, #< try changing this observation error standard deviation
B = 0 # intercept
)Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
TMB was built with Matrix version 1.5.4
Current Matrix version is 1.5.3
Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' packageWarning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
glmmTMB was built with TMB version 1.9.3
Current TMB version is 1.9.4
Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)# head(sim_dat)
ggplot(sim_dat, aes(X, Y, fill = observed)) +
geom_raster()
Exercise:
- Try running the previous chunk over and over. Try changing the
rangeandsigma_Ovalues. What do you observe? - Try increasing the observation error standard deviation
phiand re-running it many times. Try making it equal tosigma_Oand twice its size. What happens when observation error is much larger thansigma_O?
Note that the colour scale may be changing between runs!
Spatiotemporal simulation
predictor_dat2 <- expand.grid(
X = seq(0, 1, length.out = 100), Y = seq(0, 1, length.out = 100), year = 1:6
)
mesh2 <- make_mesh(predictor_dat2, c("X", "Y"), cutoff = 0.1)sim_dat2 <- sdmTMB_simulate(
formula = ~ 1,
data = predictor_dat2,
time = "year",
mesh = mesh2,
family = gaussian(),
range = 0.5, #< spatial/spatiotemporal range
sigma_O = 0.2, #< spatial SD
sigma_E = 0.1, #< spatiotemporal SD
phi = 0.001,
B = 0 # intercept
)
# head(sim_dat2)ggplot(sim_dat2, aes(X, Y, fill = omega_s)) +
geom_raster() +
facet_wrap(~year) +
scale_fill_gradient2()
ggplot(sim_dat2, aes(X, Y, fill = epsilon_st)) +
geom_raster() +
facet_wrap(~year) +
scale_fill_gradient2()
ggplot(sim_dat2, aes(X, Y, fill = observed)) +
geom_raster() +
facet_wrap(~year) +
scale_fill_gradient2()
Exercise:
- Try running the previous 4 chunks over and over and observe the various components.
- What do you notice about
omega_s(spatial random field) across the time slices? Conversely, what do you notice aboutepsilon_st(spatiotemporal random fields) across time? - Try increasing
sigma_E(spatiotemporal SD) to be larger thansigma_O(spatial SD). How do the observations change?
Correlated spatiotemporal fields
The above spatiotemporal fields were left as their default: independent across years (spatiotemporal = "iid"). They can also be correlated. Here we will simulate and plot spatiotemporal fields with various levels of correlation.
sim_dat3 <- sdmTMB_simulate(
formula = ~ 1,
data = predictor_dat2,
time = "year",
mesh = mesh2,
family = gaussian(),
range = 0.5,
sigma_O = 0,
sigma_E = 0.1,
spatiotemporal = "ar1", #< field type
rho = 0.6, #< AR(1) field correlation
phi = 0.001,
B = 0 # intercept
)
ggplot(sim_dat3, aes(X, Y, fill = epsilon_st)) +
geom_raster() +
facet_wrap(~year) +
scale_fill_gradient2()
Exercise:
- In the previous chunk, try setting the AR(1) random field correlation (
rho) to values ranging from 0 to near 1 (e.g., 0.99). What do you observe? - When might you observe a process that resembles the above?
- What happens when
rhois < 0? Isrhomore likely to be > or < than 0 in ecological systems? Can you think of a situation whererho < 0?
Fitting a geostatistical spatial model
Let’s fit a basic spatial model to some simulated data. We’ll simulate on a grid and then down-sample from that grid. (We’ll do it this way so we can plot a smooth known random field in the end for comparison.)
predictor_dat <- expand.grid(
X = seq(0, 1, length.out = 100),
Y = seq(0, 1, length.out = 100)
)
set.seed(19201) # for consistency
predictor_dat$covariate <- rnorm(n = nrow(predictor_dat))
mesh <- make_mesh(predictor_dat, c("X", "Y"), cutoff = 0.05)
plot(mesh)
sim_dat4 <- sdmTMB_simulate(
formula = ~ 1 + covariate,
data = predictor_dat,
mesh = mesh,
family = gaussian(link = "identity"),
range = 0.5,
sigma_O = 0.4,
phi = 0.3,
B = c(0.4, 0.2), # intercept and covariate coefficient
seed = 98259 # for consistency
)
head(sim_dat4) X Y omega_s mu eta observed (Intercept) covariate
1 0.00000000 0 0.20377670 0.6369297 0.6369297 0.4566973 1 0.1657648
2 0.01010101 0 0.17837392 0.6613120 0.6613120 0.9819005 1 0.4146902
3 0.02020202 0 0.15297114 0.3625643 0.3625643 0.4087907 1 -0.9520340
4 0.03030303 0 0.12756836 0.3930963 0.3930963 0.9683649 1 -0.6723601
5 0.04040404 0 0.10216557 0.6080553 0.6080553 0.9007784 1 0.5294488
6 0.05050505 0 0.07676279 0.3461807 0.3461807 0.2452402 1 -0.6529106ggplot(sim_dat4, aes(X, Y, fill = omega_s)) +
geom_raster()
ggplot(sim_dat4, aes(X, Y, fill = mu)) +
geom_raster()
ggplot(sim_dat4, aes(X, Y, fill = observed)) +
geom_raster()
Exercise:
- What did we just simulate above? What is the one addition from the first spatial simulation?
- What do the last 3 plots represent? Why do they get increasingly noisy?
Now we’ll down-sample to simulate sampling/observing only 300 locations:
set.seed(1098)
to_sample <- sample(seq_len(nrow(sim_dat4)), size = 300)
obs_dat <- sim_dat4[to_sample, ]
ggplot(obs_dat, aes(X, Y, colour = observed)) +
geom_point()
Start by fitting a regular GLM:
fit_glm <- glm(observed ~ 1 + covariate, family = gaussian(), data = obs_dat)
summary(fit_glm)
Call:
glm(formula = observed ~ 1 + covariate, family = gaussian(),
data = obs_dat)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.48113 -0.36757 0.00013 0.31924 1.33303
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.81382 0.02850 28.55 < 2e-16 ***
covariate 0.19316 0.02721 7.10 9.2e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.243673)
Null deviance: 84.897 on 299 degrees of freedom
Residual deviance: 72.615 on 298 degrees of freedom
AIC: 431.78
Number of Fisher Scoring iterations: 2We can calculate and plot the residuals spatially:
obs_dat$resid <- residuals(fit_glm)
ggplot(obs_dat, aes(X, Y, colour = resid)) +
geom_point() +
scale_colour_gradient2()
And plot predictions spatially on the full grid:
sim_dat4$pred_glm <- predict(fit_glm, newdata = sim_dat4)
lims <- range(c(sim_dat4$mu, sim_dat4$pred_glm))
g1 <- ggplot(sim_dat4, aes(X, Y, fill = mu)) +
geom_raster() + ggtitle("Truth") +
scale_fill_viridis_c(limits = lims)
g2 <- ggplot(sim_dat4, aes(X, Y, fill = pred_glm)) +
geom_raster() + ggtitle("Estimated") +
scale_fill_viridis_c(limits = lims)
cowplot::plot_grid(g1, g2, nrow = 2)
Exercise:
- Do you see spatial correlation in the residuals?
- What are at least 2 problems this creates?
- How good do those predictions look compared to the known truth?
- In the context of this workshop, what is a solution to this?
First, let’s convince ourselves that we can fit the exact same model with sdmTMB. All we have to do is turn off the spatial random fields:
fit_glm2 <- sdmTMB(observed ~ 1 + covariate, family = gaussian(),
data = obs_dat, spatial = "off") # default is spatial = "on"
summary(fit_glm2)Model fit by ML ['sdmTMB']
Formula: observed ~ 1 + covariate
Mesh: NULL (isotropic covariance)
Data: obs_dat
Family: gaussian(link = 'identity')
coef.est coef.se
(Intercept) 0.81 0.03
covariate 0.19 0.03
Dispersion parameter: 0.49
ML criterion at convergence: 212.889
See ?tidy.sdmTMB to extract these values as a data frame.Exercise:
- Check if
fit_glmandfit_glm2match and represent the same model.
Now, try fitting a GLM with Gaussian random fields:
fitting_mesh <- make_mesh(obs_dat, c("X", "Y"), cutoff = 0.05)
plot(fitting_mesh)
fit_sdmTMB <- sdmTMB(
observed ~ 1 + covariate,
family = gaussian(),
data = obs_dat,
spatial = "on",
mesh = fitting_mesh
)
sanity(fit_sdmTMB)✔ 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 largesummary(fit_sdmTMB)Spatial model fit by ML ['sdmTMB']
Formula: observed ~ 1 + covariate
Mesh: fitting_mesh (isotropic covariance)
Data: obs_dat
Family: gaussian(link = 'identity')
coef.est coef.se
(Intercept) 0.81 0.19
covariate 0.20 0.02
Dispersion parameter: 0.29
Matérn range: 0.54
Spatial SD: 0.34
ML criterion at convergence: 104.334
See ?tidy.sdmTMB to extract these values as a data frame.tidy(fit_sdmTMB, conf.int = TRUE)# A tibble: 2 × 5
term estimate std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.811 0.192 0.435 1.19
2 covariate 0.199 0.0176 0.164 0.233tidy(fit_sdmTMB, "ran_pars", conf.int = TRUE)Standard errors intentionally omitted because they have been calculated in log
space.# A tibble: 3 × 5
term estimate std.error conf.low conf.high
<chr> <dbl> <lgl> <dbl> <dbl>
1 range 0.540 NA 0.309 0.943
2 phi 0.295 NA 0.268 0.325
3 sigma_O 0.344 NA 0.245 0.483summary(fit_glm2)Model fit by ML ['sdmTMB']
Formula: observed ~ 1 + covariate
Mesh: NULL (isotropic covariance)
Data: obs_dat
Family: gaussian(link = 'identity')
coef.est coef.se
(Intercept) 0.81 0.03
covariate 0.19 0.03
Dispersion parameter: 0.49
ML criterion at convergence: 212.889
See ?tidy.sdmTMB to extract these values as a data frame.summary(fit_sdmTMB)Spatial model fit by ML ['sdmTMB']
Formula: observed ~ 1 + covariate
Mesh: fitting_mesh (isotropic covariance)
Data: obs_dat
Family: gaussian(link = 'identity')
coef.est coef.se
(Intercept) 0.81 0.19
covariate 0.20 0.02
Dispersion parameter: 0.29
Matérn range: 0.54
Spatial SD: 0.34
ML criterion at convergence: 104.334
See ?tidy.sdmTMB to extract these values as a data frame.Exercise:
- How did the dispersion parameter (Gaussian SD) differ between
fit_glm2andfit_sdmTMB? Why? - How close are the estimated values to the truth? Is that about what you expected?
Now let’s calculate and plot residuals for the random field model:
obs_dat$resid_sdmTMB <- residuals(fit_sdmTMB)
ggplot(obs_dat, aes(X, Y, colour = resid_sdmTMB)) +
geom_point() +
scale_colour_gradient2()
And make predictions on the full grid:
p <- predict(fit_sdmTMB, newdata = sim_dat4)lims <- range(c(sim_dat4$mu, p$est))
g1 <- ggplot(sim_dat4, aes(X, Y, fill = mu)) +
geom_raster() + ggtitle("True mean") +
scale_fill_viridis_c(limits = lims)
g2 <- ggplot(p, aes(X, Y, fill = est)) +
geom_raster() + ggtitle("Estimated value") +
scale_fill_viridis_c(limits = lims)
cowplot::plot_grid(g1, g2, nrow = 2)
lims <- range(c(sim_dat4$omega_s, p$omega_s))
g1 <- ggplot(sim_dat4, aes(X, Y, fill = omega_s)) +
geom_raster() + ggtitle("True random field") +
scale_fill_gradient2(limits = lims)
g2 <- ggplot(p, aes(X, Y, fill = omega_s)) +
geom_raster() + ggtitle("Estimated random field") +
scale_fill_gradient2(limits = lims)
cowplot::plot_grid(g1, g2, nrow = 2)
Predicted vs. true for the two models:
par(mfrow = c(2, 1))
plot(sim_dat4$mu, sim_dat4$pred_glm, main = "GLM",
xlab = "True", ylab = "Predicted")
plot(sim_dat4$mu, p$est, main = "GLM with random fields",
xlab = "True", ylab = "Predicted")
Exercise:
- How do the residuals look now? Why?
- How do the predictions look compared to before? Why?
- How do the random fields look? Is this as good as you expected? What would make these easier/harder to estimate?
Goals:
- Practice fitting a basic spatiotemporal model.
- Gain experience using cross validation for model selection.
- Understand how to inspect the model output.
- Practice predicting from the model on new data and making visualizations of those predictions.
- Gain familiarity with fitting and interpreting different random field structures.
The data
We will work with data representing Pacific Cod in the West Coast Vancouver Island synoptic trawl survey.
dat <- pcodhead(dat)# A tibble: 6 × 12
year X Y depth density present lat lon depth_mean depth_sd
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2003 446. 5793. 201 113. 1 52.3 -130. 5.16 0.445
2 2003 446. 5800. 212 41.7 1 52.3 -130. 5.16 0.445
3 2003 449. 5802. 220 0 0 52.4 -130. 5.16 0.445
4 2003 437. 5802. 197 15.7 1 52.4 -130. 5.16 0.445
5 2003 421. 5771. 256 0 0 52.1 -130. 5.16 0.445
6 2003 418. 5772. 293 0 0 52.1 -130. 5.16 0.445
# ℹ 2 more variables: depth_scaled <dbl>, depth_scaled2 <dbl>The dataset contains sampling locations (X and Y) and year. It also contains sampling depth in meters and sample density density (CPUE) in units of tonnes/km2.
ggplot(dat, aes(X,Y, size = density)) + geom_point()
And check to make sure that looks right:
ggplot(dat, aes(X, Y, size = density)) + geom_point(shape = 21) + coord_fixed()
We can also plot the data by year:
ggplot(dat, aes(X, Y, size = density, colour = log(density + 1))) +
geom_point(alpha = 0.3) +
facet_wrap(~year) + coord_fixed()
We will create these new columns to use later:
dat$log_depth <- log(dat$depth)
dat$year_factor <- as.factor(dat$year)Constructing a mesh
We start by constructing an SPDE mesh with INLA. This creates some matrices that are used internally when fitting the model. We will use the shortcut function make_mesh() and use a cutoff (minimum triangle length of 10 km). Note: this is the only parameter that can be changed in make_mesh() but more control over mesh arguments can be done with INLA (inla.mesh.2d()).
mesh <- make_mesh(dat, xy_cols = c("X", "Y"), cutoff = 10)
plot(mesh)
mesh$mesh$n # number of vertices or knots[1] 205# ggplot alternative:
ggplot() + inlabru::gg(mesh$mesh) + coord_fixed() +
geom_point(aes(X, Y), data = dat, alpha = 0.2, size = 0.5)
Exercise:
- Try adjusting the size of the cutoff distance to explore the effects of decreasing the mesh resolution. Make sure to reset the
cutoffvalue to10in the end so the rest of the exercise behaves as intended, because model convergence issues can be caused by meshes that are either too fine or too coarse.
Fitting a spatial model
The most basic model we could fit would be a model with a single spatial random field, and no covariates. Using silent = FALSE lets us see what is happening, but it’s awkward when running code in Rmd/Quarto chunks (with the default RStudio setting to ‘Show output inline for all R Markdown document’) so we will comment it out in most cases here. But it is a good idea to use it if models are running slowly or not converging to monitor progress.
fit_spatial <- sdmTMB(
density ~ 1, # intercept only
data = dat,
family = tweedie(link = "log"),
mesh = mesh,
spatial = "on",
# silent = FALSE
)sanity(fit_spatial)✔ 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 largeDid it have trouble? If so, try an extra optimization run and see if that’s sufficient:
fit_spatial <- run_extra_optimization(fit_spatial)
sanity(fit_spatial)✔ 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 largefit_spatialSpatial model fit by ML ['sdmTMB']
Formula: density ~ 1
Mesh: mesh (isotropic covariance)
Data: dat
Family: tweedie(link = 'log')
coef.est coef.se
(Intercept) 2.74 0.35
Dispersion parameter: 13.72
Tweedie p: 1.60
Matérn range: 20.65
Spatial SD: 2.65
ML criterion at convergence: 6566.570
See ?tidy.sdmTMB to extract these values as a data frame.Exercise:
A refresher:
- What is the ‘Matern range’ in the output?
- What is the ‘Spatial SD’ in the output?
- What is the ‘Dispersion parameter’ in the output?
Fitting a model with spatial + spatiotemporal fields
This first model includes a quadratic effect of log depth (poly(log_depth, 2)), a factor effect for each year, and models total density using a Tweedie distribution and a log link. The spatial field and spatiotemporal fields are estimated.
The year factors give each year its own mean (this is generally the approach used in fisheries stock assessment). The 0 + omits the intercept such that each year’s estimate represents a mean as opposed to a difference from the intercept. This part is arbitrary and chosen for the sake of this exercise.
The choice to use log depth helps the model fit because we have fewer samples from the deeper depths.
fit <- sdmTMB(
density ~ 0 + year_factor + poly(log_depth, 2),
data = dat,
family = tweedie(link = "log"),
mesh = mesh,
spatial = "on",
spatiotemporal = "iid", #< new
time = "year", #< new
# silent = FALSE
)Exercise:
- Skim the help file (
?sdmTMB). What are the other options forspatiotemporal?
Print the model:
fitSpatiotemporal model fit by ML ['sdmTMB']
Formula: density ~ 0 + year_factor + poly(log_depth, 2)
Mesh: mesh (isotropic covariance)
Time column: year
Data: dat
Family: tweedie(link = 'log')
coef.est coef.se
year_factor2003 1.83 0.29
year_factor2004 2.39 0.27
year_factor2005 2.17 0.27
year_factor2007 0.98 0.29
year_factor2009 1.49 0.28
year_factor2011 1.96 0.28
year_factor2013 2.31 0.27
year_factor2015 2.12 0.27
year_factor2017 1.39 0.29
poly(log_depth, 2)1 -63.66 6.10
poly(log_depth, 2)2 -95.54 5.87
Dispersion parameter: 11.01
Tweedie p: 1.50
Matérn range: 14.70
Spatial SD: 1.54
Spatiotemporal SD: 1.65
ML criterion at convergence: 6262.261
See ?tidy.sdmTMB to extract these values as a data frame.Exercise:
- What is new in the output of
fitcompared tofit_spatial?
Run a basic sanity check on the model:
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 largeWe can use the tidy function to get the coefficients for the fixed effects:
tidy(fit, conf.int = TRUE)# A tibble: 11 × 5
term estimate std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 year_factor2003 1.83 0.285 1.28 2.39
2 year_factor2004 2.39 0.267 1.87 2.92
3 year_factor2005 2.17 0.270 1.64 2.70
4 year_factor2007 0.983 0.289 0.416 1.55
5 year_factor2009 1.49 0.283 0.932 2.04
6 year_factor2011 1.96 0.276 1.42 2.50
7 year_factor2013 2.31 0.266 1.79 2.83
8 year_factor2015 2.12 0.274 1.59 2.66
9 year_factor2017 1.39 0.287 0.823 1.95
10 poly(log_depth, 2)1 -63.7 6.10 -75.6 -51.7
11 poly(log_depth, 2)2 -95.5 5.87 -107. -84.0 and by setting effects = "ran_pars", we can extract the variance and random effect components:
tidy(fit, effects = "ran_pars", conf.int = 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 14.7 NA 10.0 21.5
2 phi 11.0 NA 10.3 11.8
3 sigma_O 1.54 NA 1.18 2.01
4 sigma_E 1.65 NA 1.31 2.07
5 tweedie_p 1.50 NA 1.48 1.52There are multiple ways we can plot the depth effect on density. First, we can create a new data frame of all potential values, setting the other predictors (here year) to a fixed value.
nd <- data.frame(log_depth = seq(log(50), log(700), length.out = 100), year = 2004)
# (picking any one year)
nd$year_factor <- as.factor(nd$year)
p <- predict(
fit,
newdata = nd,
re_form = ~ 0, # means only include the fixed effects (not the default)
se_fit = TRUE # means calculate standard errors (not the default)
)
ggplot(p, aes(log_depth, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4)
The second approach is to pass the sdmTMB object to the visreg package. This shows the conditional effect, where all values other than depth are held at a particular value. Note the default visreg plot is in link (here, log) space and the dots are randomized quantile residuals.
visreg::visreg(fit, xvar = "log_depth")
visreg::visreg(fit, xvar = "log_depth", scale = "response")
Third, we could use the ggeffects package, which can be used to show the marginal effects of predictors (averaging over all other covariates rather than using a single fixed value). For more details see the visualizing marginal effects vignette. Note that won’t yet work with smoother s() terms, but will soon work with the similar ggeffects::ggpredict() .
g <- ggeffects::ggeffect(fit, "log_depth [3.5:6.7 by=0.05]")
plot(g)
Prediction
Let’s now predict on a grid that covers the entire survey (wcvi_grid).
grid <- replicate_df(qcs_grid, "year", unique(pcod$year))
grid$log_depth <- log(grid$depth)
grid$year_factor <- as.factor(grid$year)We can predict on the original data:
p0 <- predict(fit)To predict on a new data frame, we can specify newdata. Here, we will predict on the survey grid. (1) This makes it easy to make visualizations. (2) This will be useful if we wanted to generate an area-weighted standardized population index later.
p <- predict(fit, newdata = grid)We can plot each of the components of the prediction data frame spatially:
# Depth and year effect contribution:
# (Everything not a random field)
ggplot(p, aes(X, Y, fill = exp(est_non_rf))) +
facet_wrap(~year) +
geom_raster() +
coord_fixed()
# Spatial random field:
ggplot(p, aes(X, Y, fill = omega_s)) +
facet_wrap(~year) +
geom_raster() +
scale_fill_gradient2() +
coord_fixed()
# Spatial-temporal random field:
ggplot(p, aes(X, Y, fill = epsilon_st)) +
facet_wrap(~year) +
geom_raster() +
scale_fill_gradient2() +
coord_fixed()
# Overall estimate of density in link (log) space:
ggplot(p, aes(X, Y, fill = est)) +
facet_wrap(~year) +
geom_raster() +
coord_fixed()
# Overall estimate of density: (with log-distributed colour)
ggplot(p, aes(X, Y, fill = exp(est))) +
facet_wrap(~year) +
geom_raster() +
coord_fixed() +
scale_fill_viridis_c(trans = "log10")
Exercise:
- Step through each of the plots above and discuss what each represents.
Residual checking
We can calculate randomized quantile residuals with the residuals.sdmTMB() function.
dat$resid <- residuals(fit)We can plot those residuals spatially:
ggplot(dat, aes(X, Y, colour = resid)) +
facet_wrap(~year) +
geom_point(size = 0.5) +
coord_fixed() +
scale_colour_gradient2()
Exercise:
- What are you looking for in the above? Does it look OK?
- Usually visual inspection is fine but what kind of formal test could you use for spatial correlation?
We can check the distribution of the residuals with a QQ plot:
qqnorm(dat$resid)
qqline(dat$resid)
These don’t look great, but, this is largely a function of error from the Laplace approximation on the random effects (see Thygesen et al. 2017). MCMC-based are a better approach, but are slower. The following line of code will take a while to run and will require rstan and tmbstan. You may not want to run it. There’s a whole vignette on residual checking with sdmTMB here: https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html It’s an active area of research.
# warning: will take a couple minutes and requires rstan!
# see ?residuals.sdmTMB
set.seed(1)
dat$resid_mcmc <- residuals(
fit,
type = "mle-mcmc",
mcmc_iter = 201,
mcmc_warmup = 200
)
qqnorm(dat$resid_mcmc)
qqline(dat$resid_mcmc)We can simulate new observations from the model and check properties of the simulation compared to the observed data:
s <- simulate(fit, nsim = 100)mean(s == 0)[1] 0.555063mean(dat$density == 0)[1] 0.5380308hist(s, xlim = c(0, 30), breaks = 200, freq = FALSE)
hist(dat$density, xlim = c(0, 30), breaks = 200, freq = FALSE)
Exercise:
- What is the difference between
sdmTMB_simulate()andsimulate.sdmTMB()?
Fitting an AR(1) fields model
An alternative model would use AR(1) (first order autoregressive) correlated fields with or without spatial fields. We’ll fit a version without spatial fields.
fit_ar1 <- update(fit, spatial = "off", spatiotemporal = "ar1")sanity(fit_ar1)✔ 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 largefit_ar1Spatiotemporal model fit by ML ['sdmTMB']
Formula: density ~ 0 + year_factor + poly(log_depth, 2)
Mesh: mesh (isotropic covariance)
Time column: year
Data: dat
Family: tweedie(link = 'log')
coef.est coef.se
year_factor2003 1.97 0.28
year_factor2004 2.53 0.26
year_factor2005 2.29 0.27
year_factor2007 1.07 0.29
year_factor2009 1.56 0.28
year_factor2011 2.04 0.27
year_factor2013 2.32 0.26
year_factor2015 2.13 0.27
year_factor2017 1.42 0.29
poly(log_depth, 2)1 -60.21 5.60
poly(log_depth, 2)2 -93.85 5.73
Dispersion parameter: 11.24
Tweedie p: 1.50
Spatiotemporal AR1 correlation (rho): 0.62
Matérn range: 15.53
Spatiotemporal SD: 2.11
ML criterion at convergence: 6284.092
See ?tidy.sdmTMB to extract these values as a data frame.AIC(fit, fit_ar1) df AIC
fit 16 12556.52
fit_ar1 16 12600.18tidy(fit_ar1, "ran_pars", conf.int = 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 15.5 NA 10.7 22.6
2 phi 11.2 NA 10.5 12.0
3 sigma_E 2.11 NA 1.74 2.56
4 tweedie_p 1.50 NA 1.48 1.53
5 rho 0.618 NA 0.456 0.740Exercise:
- What is the estimated random field AR(1) correlation value (
rho)? (Hint: seeprint(fit_ar1)and thetidy()line above.) - What does this mean? Does this make sense ecologically?
- What is a subtle reason why you might want to include spatial fields if you’re using AR(1) fields? (Hint: think about the first time step and what the random field SD represents in the AR(1) process.)