Fix time-varying vignette example #398

Author: seananderson

vignettes/articles/basic-intro.Rmd

@@ -46,12 +46,12 @@ glimpse(pcod)
 ```
 
 The most basic model structure possible in sdmTMB replicates a GLM as can be fit with `glm()` or a GLMM as can be fit with lme4 or glmmTMB, for example.
-The spatial components in sdmTMB are included as random fields using a triangulated mesh with vertices, known as knots, used to approximate the spatial variability in observations. 
+The spatial components in sdmTMB are included as random fields using a triangulated mesh with vertices, known as knots, used to approximate the spatial variability in observations.
 Bilinear interpolation is used to approximate a continuous spatial field (Rue et al., 2009; Lindgren et al., 2011) from the estimated values of the spatial surface at these knot locations to other locations including those of actual observations.
 These spatial random effects are assumed to be drawn from Gaussian Markov random fields (e.g., Cressie & Wikle, 2011; Lindgren et al., 2011) with covariance matrices that are constrained by Matérn covariance functions (Cressie & Wikle, 2011).
 
-There are different options for creating the spatial mesh (see `sdmTMB::make_mesh()`). 
-We will start with a relatively coarse mesh for a balance between speed and accuracy (`cutoff = 10`, where cutoff is in the units of X and Y (km here) and represents the minimum distance between knots before a new mesh vertex is added). 
+There are different options for creating the spatial mesh (see `sdmTMB::make_mesh()`).
+We will start with a relatively coarse mesh for a balance between speed and accuracy (`cutoff = 10`, where cutoff is in the units of X and Y (km here) and represents the minimum distance between knots before a new mesh vertex is added).
 Smaller values create meshes with more knots.
 You will likely want to use a higher resolution mesh (more knots) in applied scenarios, but care must be taken to avoid overfitting.
 The circles represent observations and the vertices are the knot locations.
@@ -146,8 +146,8 @@ And similarly for the random effect and variance parameters:
 tidy(m3, "ran_pars", conf.int = TRUE)
 ```
 
-Note that standard errors are not reported when coefficients are in log space, but the confidence intervals are reported. These parameters are defined as follows: 
- 
+Note that standard errors are not reported when coefficients are in log space, but the confidence intervals are reported. These parameters are defined as follows:
+
 * `range`: A derived parameter that defines the distance at which 2 points are effectively independent (actually about 13% correlated). If the `share_range` argument is changed to `FALSE` then the spatial and spatiotemporal ranges will be unique, otherwise the default is for both to share the same range.
 
 * `phi`: Observation error scale parameter (e.g., SD in Gaussian).
@@ -181,11 +181,11 @@ ggplot(pcod, aes(X, Y, col = resids)) +
   coord_fixed()
 ```
 
-Those were fast to calculate but can look 'off' even when the model is consistent with the data. 
+Those were fast to calculate but can look 'off' even when the model is consistent with the data.
 MCMC-based residuals are more reliable but slow.
 We can calculate them through some help from the sdmTMBextra package. <https://github.com/pbs-assess/sdmTMBextra>.
-In practice you would like want more `mcmc_iter` and `mcmc_warmup`. 
-Total samples is `mcmc_iter - mcmc_warmup`. 
+In practice you would like want more `mcmc_iter` and `mcmc_warmup`.
+Total samples is `mcmc_iter - mcmc_warmup`.
 We will also just use the spatial model here so the vignette builds quickly.
 
 ```{r residuals-mcmc, fig.width = 8}
@@ -229,7 +229,7 @@ plot_map <- function(dat, column) {
 }
 ```
 
-There are four kinds of predictions that we get out of the model. 
+There are four kinds of predictions that we get out of the model.
 
 First, we will show the predictions that incorporate all fixed effects and random effects:
 
@@ -254,7 +254,7 @@ plot_map(predictions, exp(est_non_rf)) +
   ggtitle("Prediction (fixed effects only)")
 ```
 
-We can look at the spatial random effects that represent consistent deviations in space through time that are not accounted for by our fixed effects. 
+We can look at the spatial random effects that represent consistent deviations in space through time that are not accounted for by our fixed effects.
 In other words, these deviations represent consistent biotic and abiotic factors that are affecting biomass density but are not accounted for in the model.
 
 ```{r plot-spatial-effects, fig.width = 6}
@@ -263,7 +263,7 @@ plot_map(predictions, omega_s) +
   ggtitle("Spatial random effects only")
 ```
 
-And finally we can look at the spatiotemporal random effects that represent deviation from the fixed effect predictions and the spatial random effect deviations. 
+And finally we can look at the spatiotemporal random effects that represent deviation from the fixed effect predictions and the spatial random effect deviations.
 These represent biotic and abiotic factors that are changing through time and are not accounted for in the model.
 
 ```{r plot-spatiotemporal-effects, fig.width = 8}
@@ -273,8 +273,8 @@ plot_map(predictions, epsilon_st) +
   ggtitle("Spatiotemporal random effects only")
 ```
 
-We can also estimate the uncertainty in our spatiotemporal density predictions using simulations from the joint precision matrix by setting `nsim > 0` in the predict function. 
-Here we generate 100 estimates and use `apply()` to calculate upper and lower confidence intervals, a standard deviation, and a coefficient of variation (CV).  
+We can also estimate the uncertainty in our spatiotemporal density predictions using simulations from the joint precision matrix by setting `nsim > 0` in the predict function.
+Here we generate 100 estimates and use `apply()` to calculate upper and lower confidence intervals, a standard deviation, and a coefficient of variation (CV).
 
 ```{r sim-cv}
 sim <- predict(m3, newdata = grid_yrs, nsim = 100)
@@ -333,28 +333,32 @@ ggeffects::ggeffect(m3,  "depth [0:500 by=1]") %>% plot()
 
 ### Time-varying effects
 
-We could also let the effect of depth vary through time. 
-In this example, it helps to give each year a separate mean effect (`as.factor(year)`). 
-The `~ 0` part of the formula omits the intercept.
-For models like this that take longer to fit, we might want to set `silent = FALSE` so we can monitor progress.
+We could also let the effect of depth vary through time.
+We set up the time-varying coefficients to follow an AR1 process by setting `type_varying_type = "ar1"`.
+With `"ar1"` or `"rw0" (random walk), the fixed effects represent the starting point of the time series
+and the time-varying process represents deviations from this over time.
+We include a full length of time increments with `extra_time` to ensure we estimate time-varying coefficient values
+for each year, including any years that are missing from our data.
+For this example, we turn off the spatiotemporal random effects because we were having convergence issues with them turned on.
 
 ```{r tv-effect}
 m4 <- sdmTMB(
-  density ~ 0 + as.factor(year),
+  density ~ 1 + depth_scaled + depth_scaled2,
   data = pcod,
-  time_varying = ~ 0 + depth_scaled + depth_scaled2,
-  time_varying_type = "rw0",
+  time_varying = ~ 1 + depth_scaled + depth_scaled2,
+  time_varying_type = "ar1",
+  extra_time = seq(min(pcod$year), max(pcod$year)),
   mesh = mesh,
   family = tweedie(link = "log"),
   spatial = "on",
   time = "year",
-  spatiotemporal = "IID"
+  spatiotemporal = "off"
 )
 m4
-AIC(m4)
 ```
 
-To plot these, we make a data frame that contains all combinations of the time-varying covariate and time. This is easily created using `expand.grid()` or `tidyr::expand_grid()`.
+To plot these, we make a data frame that contains all combinations of the time-varying covariate and time.
+This is easily created using `expand.grid()` or `tidyr::expand_grid()`.
 
 ```{r tv-depth-eff, fig.width = 6}
 nd <- expand.grid(