From 82e54cfd708dd4ac43075cc6854812b44f1542d7 Mon Sep 17 00:00:00 2001
From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com>
Date: Mon, 19 Aug 2024 08:36:13 -0700
Subject: [PATCH 1/7] Multispecies / multiresponse vignette

I talked to Lewis last week about this. Currently includes examples of:
- main effects by species
- spatial fields by species, with and without shared variance
- spatiotemporal fields by species
---
 vignettes/web_only/multispecies.Rmd | 247 ++++++++++++++++++++++++++++
 1 file changed, 247 insertions(+)
 create mode 100644 vignettes/web_only/multispecies.Rmd

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
new file mode 100644
index 000000000..10bf8c21e
--- /dev/null
+++ b/vignettes/web_only/multispecies.Rmd
@@ -0,0 +1,247 @@
+---
+title: "Fitting multispecies models with sdmTMB"
+date: "`r Sys.Date()`"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Fitting multispecies models with sdmTMB}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+**If the code in this vignette has not been evaluated, a rendered version is available on the [documentation site](https://pbs-assess.github.io/sdmTMB/index.html) under 'Articles'.**
+
+```{r setup, include = FALSE, cache=FALSE}
+dplyr_installed <- require("dplyr", quietly = TRUE)
+ggplot_installed <- require("ggplot2", quietly = TRUE)
+pkgs <- dplyr_installed && ggplot_installed
+EVAL <- identical(Sys.getenv("NOT_CRAN"), "true") && pkgs
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  warning=FALSE,
+  message=FALSE,
+  comment = "#>",
+  fig.width = 7,
+  fig.asp = 0.618,
+  eval = EVAL,
+  purl = EVAL
+)
+```
+
+```{r packages, message=FALSE, warning=TRUE}
+library(ggplot2)
+library(dplyr)
+library(sdmTMB)
+```
+
+For some applications, we might be interested in fitting a model that includes multiple responses (such as 2+ species). The most important step in fitting these models is understanding which parameters are shared, and which parameters are species specific.
+
+Below, we illustrate a series of models that may be used to answer different questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial ranges and variances, as well as different year effects.
+
+```{r sim_dat}
+  set.seed(123)
+
+  # make fake predictor(s) (a1) and sampling locations:
+  predictor_dat <- data.frame(
+    X = runif(500), Y = runif(500),
+    year = rep(1:5, each = 100)
+  )
+  predictor_dat$fyear <- as.factor(predictor_dat$year)
+  mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
+
+  sim_dat_A <- sdmTMB_simulate(
+    formula = ~ 0 + fyear,
+    data = predictor_dat,
+    time = "year",
+    mesh = mesh,
+    range = 0.5,
+    family = tweedie(),
+    seed = 42,
+    sigma_O = 0.2,
+    phi = 0.1,
+    sigma_E = 0.1,
+    B = runif(5, min = 5, max = 8) # 5 random year effects
+  )
+  sim_dat_A$species <- "A"
+  
+  sim_dat_B <- sdmTMB_simulate(
+    formula = ~ 0 + fyear,
+    data = predictor_dat,
+    time = "year",
+    mesh = mesh,
+    range = 0.2,
+    family = gaussian(),
+    seed = 42,
+    sigma_O = 0.3,
+    phi = 0.1,
+    sigma_E = 0.1,
+    B = runif(5, min = 5, max = 8) # 5 random year effects
+  )
+  sim_dat_B$species <- "B"
+  
+  sim_dat <- rbind(sim_dat_A, sim_dat_B)
+  sim_dat$fyear <- as.factor(sim_dat$year)
+```
+
+We'll start by making the mesh across the full dataset
+
+```{r mesh_fig, fig.asp=1}
+spde <- make_mesh(sim_dat, c("X", "Y"), cutoff = 0.1)
+plot(spde)
+```
+
+### Model 1: species specific intercepts
+
+As a first model, we can include species specific year effects. This can be done in a couple ways -- one option would be to estimate the `species*year` interaction, letting the year effects for each species be independent. All other parameters (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared.
+
+```{r}
+fit <- sdmTMB(
+  formula = observed ~ fyear*species,
+  data = sim_dat, 
+  time = "year", 
+  spatiotemporal = "iid",
+  mesh = spde, family = gaussian())
+```
+
+This model doesn't appear to converge, and the results from `sanity(fit)` indicate problems with the spatial and spatiotemporal variances. This is likely a sign that our assumptions in the estimation model are wrong for this particular dataset (not surprising, as our simulated data has different parameters by species).
+
+Our simulated dataset only has 5 years of data -- if we had a longer time series, an alternative formulation for our year effects could be to fit smooth terms, like in a GAM with `mgcv`. 
+
+```{r eval=FALSE}
+fit <- sdmTMB(
+  formula = observed ~ s(year,by="species"),
+  data = sim_dat, 
+  time = "year", 
+  spatiotemporal = "iid",
+  mesh = spde, family = gaussian())
+```
+
+### Model 2: species specific fields
+
+We may be interested in fitting a model that lets the spatial patterning differ by species. These kinds of models can be expressed using spatially varying coefficients. Note that we use `spatial = off` because this represents a global spatial intercept -- turning this off is equivalent to using a `-1` of `0` in a main formula including a factor. 
+
+```{r}
+fit <- sdmTMB(
+  observed ~ 0 + fyear, # shared year effects  
+  data = sim_dat,
+  mesh = spde, 
+  family = gaussian(),
+  spatial = "off",
+  time = "year",
+  spatiotemporal = "iid",
+  spatial_varying = ~ 0 + as.factor(species)
+)
+```
+
+This model also appears to be problematic. If we look at the estimated parameters, you'll notice that there are two row entries for `sigma_Z`. This means that our model is trying to estimate separate species specific variance terms for the species specific spatial fields. 
+
+```{r}
+fit$sd_report
+```
+
+If we wanted to estimate species specific spatial fields with a single shared variance (meaning the magnitude of the peaks and valleys in the field was similar), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`
+
+```{r}
+map_list = list(ln_tau_Z = factor(c(1, 1)))
+
+fit <- sdmTMB(
+  observed ~ 0 + fyear, # shared year effects  
+  data = sim_dat,
+  mesh = spde, 
+  family = gaussian(),
+  spatial = "off",
+  time = "year",
+  spatiotemporal = "iid",
+  spatial_varying = ~ 0 + as.factor(species),
+  control = sdmTMBcontrol(map = map_list)
+)
+```
+
+### Model 3: species specific spatiotemporal fields
+
+As a last example, we can set up a model that also includes spatiotemporal fields unique to each species. In all of the examples above, spatiotemporal fields are included -- but shared among species.
+
+One approach to including spatiotemporal fields is creating a new variable combining species and year. We can then include this spatiotemporal variation by changing the `time` argument to be `time = "species_year"`
+
+```{r}
+# This needs to be a factor
+sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
+
+map_list = list(ln_tau_Z = factor(c(1, 1)))
+
+fit <- sdmTMB(
+  observed ~ 0 + fyear, # shared year effects  
+  data = sim_dat,
+  mesh = spde, 
+  family = gaussian(),
+  spatial = "off",
+  time = "species_year",
+  spatiotemporal = "iid",
+  spatial_varying = ~ 0 + as.factor(species),
+  control = sdmTMBcontrol(map = map_list)
+)
+
+```
+
+This is getting a little better -- there is a gradient warning, but the Hessian appears to be ok.
+
+### Putting it all together
+
+Our spatiotemporal model above came close to passing all sanity() checks, and could be modified to also include the species-specific year effects that we started with
+
+```{r}
+fit <- sdmTMB(
+  observed ~ 0 + fyear*species, # shared year effects  
+  data = sim_dat,
+  mesh = spde, 
+  family = gaussian(),
+  spatial = "off",
+  time = "species_year",
+  spatiotemporal = "iid",
+  spatial_varying = ~ 0 + as.factor(species),
+  control = sdmTMBcontrol(map = map_list)
+)
+```
+
+These examples illustrate a number of ways that species - specific effects can be included in `sdmTMB` models, and can be extended to other types of groupings. A brief summary of these can be summarized as:
+
+```{r echo=FALSE}
+desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year by species interactions, or smooths", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
+knitr::kable(desc)
+```
+
+### Extensions
+
+Including spatiotemporal effects like species-year combinations in the `time` argument is relatively straightforward, but these models could get more complicated if we wanted to also include other grouping factors (such as cohort effects for each species). There may be some identifiability issues that would have to be resolved, but even if the `time` argument is being used, it's possible to add other factors through the `spatial_varying` formula. As an example, if our dataframe included a cohort variable for each species, we could first create a new combination  
+
+```{r eval=FALSE}
+sim_dat$species_cohort <- as.factor(paste(sim_dat$species, sim_dat$cohort))
+```
+
+If there were ~ 10 levels of `species_cohort`, and we were including our original spatial effects in `spatial_varying`, we could create a new map list that would indicate we wanted all species-cohort fields to have the same spatial variance, 
+
+```{r eval=FALSE}
+sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
+
+# The first 2 elements correspond to species, remaining 10 to species-cohort
+map_list = list(ln_tau_Z = factor(c(1, 1, rep(2, 10))))
+
+fit <- sdmTMB(
+  observed ~ 0 + fyear*species,
+  data = sim_dat,
+  mesh = spde, 
+  family = gaussian(),
+  spatial = "off",
+  time = "species_year",
+  spatiotemporal = "iid",
+  spatial_varying = ~ 0 + as.factor(species) + species_cohort,
+  control = sdmTMBcontrol(map = map_list)
+)
+```
+
+
+
+
+
+
+
+

From a677156a7f491e7407fdfe6916c10f25b8142c73 Mon Sep 17 00:00:00 2001
From: Lewis Barnett <lewis.barnett@noaa.gov>
Date: Fri, 22 Nov 2024 12:46:29 -0800
Subject: [PATCH 2/7] edit and show more of output

---
 vignettes/web_only/multispecies.Rmd | 82 +++++++++++++++++------------
 1 file changed, 48 insertions(+), 34 deletions(-)

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
index 10bf8c21e..7da13f0cd 100644
--- a/vignettes/web_only/multispecies.Rmd
+++ b/vignettes/web_only/multispecies.Rmd
@@ -10,7 +10,7 @@ vignette: >
 
 **If the code in this vignette has not been evaluated, a rendered version is available on the [documentation site](https://pbs-assess.github.io/sdmTMB/index.html) under 'Articles'.**
 
-```{r setup, include = FALSE, cache=FALSE}
+```{r setup, include = FALSE, cache = FALSE}
 dplyr_installed <- require("dplyr", quietly = TRUE)
 ggplot_installed <- require("ggplot2", quietly = TRUE)
 pkgs <- dplyr_installed && ggplot_installed
@@ -27,15 +27,15 @@ knitr::opts_chunk$set(
 )
 ```
 
-```{r packages, message=FALSE, warning=TRUE}
+```{r packages, message = FALSE, warning = TRUE}
 library(ggplot2)
 library(dplyr)
 library(sdmTMB)
 ```
 
-For some applications, we might be interested in fitting a model that includes multiple responses (such as 2+ species). The most important step in fitting these models is understanding which parameters are shared, and which parameters are species specific.
+For some applications, we might be interested in fitting a model that includes multiple responses (such as 2+ species, or multiple size or age classes within a species). The most important step in fitting these models is understanding which parameters are shared, and which parameters are species-specific.
 
-Below, we illustrate a series of models that may be used to answer different questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial ranges and variances, as well as different year effects.
+Below, we illustrate a series of models that may be used to answer a few different common questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial ranges and variances, as well as different year effects.
 
 ```{r sim_dat}
   set.seed(123)
@@ -82,42 +82,50 @@ Below, we illustrate a series of models that may be used to answer different que
   sim_dat$fyear <- as.factor(sim_dat$year)
 ```
 
-We'll start by making the mesh across the full dataset
+We'll start by making the mesh across the full dataset.
 
 ```{r mesh_fig, fig.asp=1}
 spde <- make_mesh(sim_dat, c("X", "Y"), cutoff = 0.1)
 plot(spde)
 ```
 
-### Model 1: species specific intercepts
+### Model 1: species-specific intercepts
 
-As a first model, we can include species specific year effects. This can be done in a couple ways -- one option would be to estimate the `species*year` interaction, letting the year effects for each species be independent. All other parameters (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared.
+As a first model, we can include species-specific year effects. This can be done in a couple ways -- one option would be to estimate the `species*year` interaction, letting the year effects for each species be independent. All other parameters (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared.
 
 ```{r}
 fit <- sdmTMB(
-  formula = observed ~ fyear*species,
+  formula = observed ~ fyear * species,
   data = sim_dat, 
   time = "year", 
   spatiotemporal = "iid",
-  mesh = spde, family = gaussian())
+  mesh = spde, 
+  family = gaussian()
+  )
+
+fit
+
+print(sanity(fit))
 ```
 
 This model doesn't appear to converge, and the results from `sanity(fit)` indicate problems with the spatial and spatiotemporal variances. This is likely a sign that our assumptions in the estimation model are wrong for this particular dataset (not surprising, as our simulated data has different parameters by species).
 
-Our simulated dataset only has 5 years of data -- if we had a longer time series, an alternative formulation for our year effects could be to fit smooth terms, like in a GAM with `mgcv`. 
+Our simulated dataset only has 5 years of data -- if we had a longer time series, an alternative formulation for our year effects could be to fit smooth terms, like in a GAM with `mgcv`.
 
 ```{r eval=FALSE}
 fit <- sdmTMB(
-  formula = observed ~ s(year,by="species"),
+  formula = observed ~ s(year, by = "species"),
   data = sim_dat, 
   time = "year", 
   spatiotemporal = "iid",
-  mesh = spde, family = gaussian())
+  mesh = spde, 
+  family = gaussian()
+  )
 ```
 
-### Model 2: species specific fields
+### Model 2: species-specific fields
 
-We may be interested in fitting a model that lets the spatial patterning differ by species. These kinds of models can be expressed using spatially varying coefficients. Note that we use `spatial = off` because this represents a global spatial intercept -- turning this off is equivalent to using a `-1` of `0` in a main formula including a factor. 
+We may be interested in fitting a model that lets the spatial patterning differ by species. These kinds of models can be expressed using spatially varying coefficients. Note that we use `spatial = off` because this represents a global spatial intercept -- turning this off is akin to using a `-1` of `0` in a main formula including a factor.
 
 ```{r}
 fit <- sdmTMB(
@@ -130,15 +138,17 @@ fit <- sdmTMB(
   spatiotemporal = "iid",
   spatial_varying = ~ 0 + as.factor(species)
 )
+
+fit
 ```
 
-This model also appears to be problematic. If we look at the estimated parameters, you'll notice that there are two row entries for `sigma_Z`. This means that our model is trying to estimate separate species specific variance terms for the species specific spatial fields. 
+This model also appears to be problematic. If we look at the estimated parameters, you'll notice that there are two rows of entries for `sigma_Z`. This means that our model is trying to estimate separate species-specific variance terms for the species-specific spatial fields (say *that* 10 times fast!).
 
 ```{r}
 fit$sd_report
 ```
 
-If we wanted to estimate species specific spatial fields with a single shared variance (meaning the magnitude of the peaks and valleys in the field was similar), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`
+If we wanted to estimate species-specific spatial fields with a single shared variance (meaning the net magnitude of the peaks and valleys in the fields were similar), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`
 
 ```{r}
 map_list = list(ln_tau_Z = factor(c(1, 1)))
@@ -154,13 +164,15 @@ fit <- sdmTMB(
   spatial_varying = ~ 0 + as.factor(species),
   control = sdmTMBcontrol(map = map_list)
 )
+
+fit
 ```
 
-### Model 3: species specific spatiotemporal fields
+### Model 3: species-specific spatiotemporal fields
 
-As a last example, we can set up a model that also includes spatiotemporal fields unique to each species. In all of the examples above, spatiotemporal fields are included -- but shared among species.
+In all of the examples above, spatiotemporal fields are included -- but shared among species. As a last example, we can extend the above approaches to set up a model that includes spatiotemporal fields unique to each species.
 
-One approach to including spatiotemporal fields is creating a new variable combining species and year. We can then include this spatiotemporal variation by changing the `time` argument to be `time = "species_year"`
+One approach to including separate spatiotemporal fields by species is creating a new variable that is a concatenation of species and year (or any given time step factor). For example, we can then implement this form of species-specific spatiotemporal variation by changing the `time` argument to be `time = "species_year"`
 
 ```{r}
 # This needs to be a factor
@@ -180,9 +192,12 @@ fit <- sdmTMB(
   control = sdmTMBcontrol(map = map_list)
 )
 
+fit
+
+print(sanity(fit))
 ```
 
-This is getting a little better -- there is a gradient warning, but the Hessian appears to be ok.
+This is getting a little better -- there is a gradient warning, but the Hessian appears to be invertible (`hessian_ok == TRUE)`, which is critical to convergence.
 
 ### Putting it all together
 
@@ -200,33 +215,40 @@ fit <- sdmTMB(
   spatial_varying = ~ 0 + as.factor(species),
   control = sdmTMBcontrol(map = map_list)
 )
+
+
+fit
+
+print(sanity(fit))
+
+fit$sd_report
 ```
 
-These examples illustrate a number of ways that species - specific effects can be included in `sdmTMB` models, and can be extended to other types of groupings. A brief summary of these can be summarized as:
+These examples illustrate a number of ways that species-specific effects can be included in `sdmTMB` models, and can be extended to other categories/groups/cohorts within a species for which one wants to control the amount of information shared between groups (e.g., age-, size-, or stage-specific estimates). A brief summary of these approaches can be summarized as:
 
 ```{r echo=FALSE}
 desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year by species interactions, or smooths", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
 knitr::kable(desc)
 ```
 
-### Extensions
+### Further extensions
 
-Including spatiotemporal effects like species-year combinations in the `time` argument is relatively straightforward, but these models could get more complicated if we wanted to also include other grouping factors (such as cohort effects for each species). There may be some identifiability issues that would have to be resolved, but even if the `time` argument is being used, it's possible to add other factors through the `spatial_varying` formula. As an example, if our dataframe included a cohort variable for each species, we could first create a new combination  
+Including spatiotemporal effects like species-year combinations in the `time` argument is relatively straightforward, but these models could get more complicated if we wanted to also include other grouping factors (such as cohort effects for each species). There may be some identifiability issues that would have to be resolved, but even if the `time` argument is being used, it's possible to add other factors through the `spatial_varying` formula. As an example, if our dataframe included a cohort variable for each species, we could first create a new factor that is a concatenation of species and cohort
 
 ```{r eval=FALSE}
 sim_dat$species_cohort <- as.factor(paste(sim_dat$species, sim_dat$cohort))
 ```
 
-If there were ~ 10 levels of `species_cohort`, and we were including our original spatial effects in `spatial_varying`, we could create a new map list that would indicate we wanted all species-cohort fields to have the same spatial variance, 
+If there were \~ 10 levels of `species_cohort`, and we were including our original spatial effects in the `spatial_varying` term, we could create a new map list that would indicate we wanted all species-cohort fields to have the same spatial variance as below.
 
 ```{r eval=FALSE}
 sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
 
-# The first 2 elements correspond to species, remaining 10 to species-cohort
+# The first 2 elements correspond to species, the remaining 10 to species-cohort
 map_list = list(ln_tau_Z = factor(c(1, 1, rep(2, 10))))
 
 fit <- sdmTMB(
-  observed ~ 0 + fyear*species,
+  observed ~ 0 + fyear * species,
   data = sim_dat,
   mesh = spde, 
   family = gaussian(),
@@ -237,11 +259,3 @@ fit <- sdmTMB(
   control = sdmTMBcontrol(map = map_list)
 )
 ```
-
-
-
-
-
-
-
-

From 6f493848e74b17836b7e07682763aa5580fba2c4 Mon Sep 17 00:00:00 2001
From: Sean Anderson <sean@seananderson.ca>
Date: Tue, 3 Dec 2024 15:05:55 -0800
Subject: [PATCH 3/7] Run styler + copy edit

---
 vignettes/web_only/multispecies.Rmd | 130 ++++++++++++++--------------
 1 file changed, 65 insertions(+), 65 deletions(-)

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
index 7da13f0cd..b0d0d6be0 100644
--- a/vignettes/web_only/multispecies.Rmd
+++ b/vignettes/web_only/multispecies.Rmd
@@ -17,8 +17,8 @@ pkgs <- dplyr_installed && ggplot_installed
 EVAL <- identical(Sys.getenv("NOT_CRAN"), "true") && pkgs
 knitr::opts_chunk$set(
   collapse = TRUE,
-  warning=FALSE,
-  message=FALSE,
+  warning = FALSE,
+  message = FALSE,
   comment = "#>",
   fig.width = 7,
   fig.asp = 0.618,
@@ -33,53 +33,53 @@ library(dplyr)
 library(sdmTMB)
 ```
 
-For some applications, we might be interested in fitting a model that includes multiple responses (such as 2+ species, or multiple size or age classes within a species). The most important step in fitting these models is understanding which parameters are shared, and which parameters are species-specific.
+For some applications, we might be interested in fitting a model that includes multiple responses such as 2+ species, or multiple size or age classes within a species. The most important step in fitting these models is understanding which parameters are shared, and which parameters are species-specific.
 
 Below, we illustrate a series of models that may be used to answer a few different common questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial ranges and variances, as well as different year effects.
 
 ```{r sim_dat}
-  set.seed(123)
-
-  # make fake predictor(s) (a1) and sampling locations:
-  predictor_dat <- data.frame(
-    X = runif(500), Y = runif(500),
-    year = rep(1:5, each = 100)
-  )
-  predictor_dat$fyear <- as.factor(predictor_dat$year)
-  mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
-
-  sim_dat_A <- sdmTMB_simulate(
-    formula = ~ 0 + fyear,
-    data = predictor_dat,
-    time = "year",
-    mesh = mesh,
-    range = 0.5,
-    family = tweedie(),
-    seed = 42,
-    sigma_O = 0.2,
-    phi = 0.1,
-    sigma_E = 0.1,
-    B = runif(5, min = 5, max = 8) # 5 random year effects
-  )
-  sim_dat_A$species <- "A"
-  
-  sim_dat_B <- sdmTMB_simulate(
-    formula = ~ 0 + fyear,
-    data = predictor_dat,
-    time = "year",
-    mesh = mesh,
-    range = 0.2,
-    family = gaussian(),
-    seed = 42,
-    sigma_O = 0.3,
-    phi = 0.1,
-    sigma_E = 0.1,
-    B = runif(5, min = 5, max = 8) # 5 random year effects
-  )
-  sim_dat_B$species <- "B"
-  
-  sim_dat <- rbind(sim_dat_A, sim_dat_B)
-  sim_dat$fyear <- as.factor(sim_dat$year)
+set.seed(123)
+
+# make fake predictor(s) (a1) and sampling locations:
+predictor_dat <- data.frame(
+  X = runif(500), Y = runif(500),
+  year = rep(1:5, each = 100)
+)
+predictor_dat$fyear <- as.factor(predictor_dat$year)
+mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
+
+sim_dat_A <- sdmTMB_simulate(
+  formula = ~ 0 + fyear,
+  data = predictor_dat,
+  time = "year",
+  mesh = mesh,
+  range = 0.5,
+  family = tweedie(),
+  seed = 42,
+  sigma_O = 0.2,
+  phi = 0.1,
+  sigma_E = 0.1,
+  B = runif(5, min = 5, max = 8) # 5 random year effects
+)
+sim_dat_A$species <- "A"
+
+sim_dat_B <- sdmTMB_simulate(
+  formula = ~ 0 + fyear,
+  data = predictor_dat,
+  time = "year",
+  mesh = mesh,
+  range = 0.2,
+  family = gaussian(),
+  seed = 42,
+  sigma_O = 0.3,
+  phi = 0.1,
+  sigma_E = 0.1,
+  B = runif(5, min = 5, max = 8) # 5 random year effects
+)
+sim_dat_B$species <- "B"
+
+sim_dat <- rbind(sim_dat_A, sim_dat_B)
+sim_dat$fyear <- as.factor(sim_dat$year)
 ```
 
 We'll start by making the mesh across the full dataset.
@@ -96,12 +96,12 @@ As a first model, we can include species-specific year effects. This can be done
 ```{r}
 fit <- sdmTMB(
   formula = observed ~ fyear * species,
-  data = sim_dat, 
-  time = "year", 
+  data = sim_dat,
+  time = "year",
   spatiotemporal = "iid",
-  mesh = spde, 
+  mesh = spde,
   family = gaussian()
-  )
+)
 
 fit
 
@@ -115,12 +115,12 @@ Our simulated dataset only has 5 years of data -- if we had a longer time series
 ```{r eval=FALSE}
 fit <- sdmTMB(
   formula = observed ~ s(year, by = "species"),
-  data = sim_dat, 
-  time = "year", 
+  data = sim_dat,
+  time = "year",
   spatiotemporal = "iid",
-  mesh = spde, 
+  mesh = spde,
   family = gaussian()
-  )
+)
 ```
 
 ### Model 2: species-specific fields
@@ -129,9 +129,9 @@ We may be interested in fitting a model that lets the spatial patterning differ
 
 ```{r}
 fit <- sdmTMB(
-  observed ~ 0 + fyear, # shared year effects  
+  observed ~ 0 + fyear, # shared year effects
   data = sim_dat,
-  mesh = spde, 
+  mesh = spde,
   family = gaussian(),
   spatial = "off",
   time = "year",
@@ -151,12 +151,12 @@ fit$sd_report
 If we wanted to estimate species-specific spatial fields with a single shared variance (meaning the net magnitude of the peaks and valleys in the fields were similar), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`
 
 ```{r}
-map_list = list(ln_tau_Z = factor(c(1, 1)))
+map_list <- list(ln_tau_Z = factor(c(1, 1)))
 
 fit <- sdmTMB(
-  observed ~ 0 + fyear, # shared year effects  
+  observed ~ 0 + fyear, # shared year effects
   data = sim_dat,
-  mesh = spde, 
+  mesh = spde,
   family = gaussian(),
   spatial = "off",
   time = "year",
@@ -178,12 +178,12 @@ One approach to including separate spatiotemporal fields by species is creating
 # This needs to be a factor
 sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
 
-map_list = list(ln_tau_Z = factor(c(1, 1)))
+map_list <- list(ln_tau_Z = factor(c(1, 1)))
 
 fit <- sdmTMB(
-  observed ~ 0 + fyear, # shared year effects  
+  observed ~ 0 + fyear, # shared year effects
   data = sim_dat,
-  mesh = spde, 
+  mesh = spde,
   family = gaussian(),
   spatial = "off",
   time = "species_year",
@@ -205,9 +205,9 @@ Our spatiotemporal model above came close to passing all sanity() checks, and co
 
 ```{r}
 fit <- sdmTMB(
-  observed ~ 0 + fyear*species, # shared year effects  
+  observed ~ 0 + fyear * species, # shared year effects
   data = sim_dat,
-  mesh = spde, 
+  mesh = spde,
   family = gaussian(),
   spatial = "off",
   time = "species_year",
@@ -245,12 +245,12 @@ If there were \~ 10 levels of `species_cohort`, and we were including our origin
 sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
 
 # The first 2 elements correspond to species, the remaining 10 to species-cohort
-map_list = list(ln_tau_Z = factor(c(1, 1, rep(2, 10))))
+map_list <- list(ln_tau_Z = factor(c(1, 1, rep(2, 10))))
 
 fit <- sdmTMB(
   observed ~ 0 + fyear * species,
   data = sim_dat,
-  mesh = spde, 
+  mesh = spde,
   family = gaussian(),
   spatial = "off",
   time = "species_year",

From 3f97c07a9d98b94f19e07582d9aa10ee5d49fc1e Mon Sep 17 00:00:00 2001
From: Sean Anderson <sean@seananderson.ca>
Date: Tue, 3 Dec 2024 15:51:29 -0800
Subject: [PATCH 4/7] Multispecies vignette edits

For one, make it converge!
---
 vignettes/web_only/multispecies.Rmd | 176 +++++++++++++++-------------
 1 file changed, 93 insertions(+), 83 deletions(-)

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
index b0d0d6be0..6154d8538 100644
--- a/vignettes/web_only/multispecies.Rmd
+++ b/vignettes/web_only/multispecies.Rmd
@@ -35,34 +35,30 @@ library(sdmTMB)
 
 For some applications, we might be interested in fitting a model that includes multiple responses such as 2+ species, or multiple size or age classes within a species. The most important step in fitting these models is understanding which parameters are shared, and which parameters are species-specific.
 
-Below, we illustrate a series of models that may be used to answer a few different common questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial ranges and variances, as well as different year effects.
+Below, we illustrate a series of models that may be used to answer a few different common questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial standard deviations (`sigma_O`) as well as different year effects.
 
 ```{r sim_dat}
-set.seed(123)
-
-# make fake predictor(s) (a1) and sampling locations:
+set.seed(1)
 predictor_dat <- data.frame(
-  X = runif(500), Y = runif(500),
-  year = rep(1:5, each = 100)
+  X = runif(1000), Y = runif(1000),
+  year = rep(1:5, each = 200)
 )
 predictor_dat$fyear <- as.factor(predictor_dat$year)
 mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
-
 sim_dat_A <- sdmTMB_simulate(
   formula = ~ 0 + fyear,
   data = predictor_dat,
   time = "year",
   mesh = mesh,
-  range = 0.5,
-  family = tweedie(),
+  range = 0.2,
+  family = gaussian(),
   seed = 42,
   sigma_O = 0.2,
   phi = 0.1,
-  sigma_E = 0.1,
+  sigma_E = 0.3,
   B = runif(5, min = 5, max = 8) # 5 random year effects
 )
 sim_dat_A$species <- "A"
-
 sim_dat_B <- sdmTMB_simulate(
   formula = ~ 0 + fyear,
   data = predictor_dat,
@@ -73,189 +69,203 @@ sim_dat_B <- sdmTMB_simulate(
   seed = 42,
   sigma_O = 0.3,
   phi = 0.1,
-  sigma_E = 0.1,
+  sigma_E = 0.3,
   B = runif(5, min = 5, max = 8) # 5 random year effects
 )
 sim_dat_B$species <- "B"
-
 sim_dat <- rbind(sim_dat_A, sim_dat_B)
-sim_dat$fyear <- as.factor(sim_dat$year)
+sim_dat$fyear <- factor(sim_dat$year)
 ```
 
 We'll start by making the mesh across the full dataset.
 
 ```{r mesh_fig, fig.asp=1}
-spde <- make_mesh(sim_dat, c("X", "Y"), cutoff = 0.1)
-plot(spde)
+mesh <- make_mesh(sim_dat, c("X", "Y"), cutoff = 0.1)
+plot(mesh)
 ```
 
 ### Model 1: species-specific intercepts
 
-As a first model, we can include species-specific year effects. This can be done in a couple ways -- one option would be to estimate the `species*year` interaction, letting the year effects for each species be independent. All other parameters (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared.
+As a first model, we can include species-specific year effects. This can be done in a couple ways. One option would be to estimate the `species * year` interaction, letting the year effects for each species be independent. Here, all other parameters (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared. Also, both species share a single spatial random field.
 
 ```{r}
 fit <- sdmTMB(
-  formula = observed ~ fyear * species,
+  observed ~ fyear * species,
   data = sim_dat,
   time = "year",
   spatiotemporal = "iid",
-  mesh = spde,
+  mesh = mesh,
   family = gaussian()
 )
-
 fit
-
-print(sanity(fit))
 ```
 
-This model doesn't appear to converge, and the results from `sanity(fit)` indicate problems with the spatial and spatiotemporal variances. This is likely a sign that our assumptions in the estimation model are wrong for this particular dataset (not surprising, as our simulated data has different parameters by species).
+### Model 2: species-specific spatial fields
 
-Our simulated dataset only has 5 years of data -- if we had a longer time series, an alternative formulation for our year effects could be to fit smooth terms, like in a GAM with `mgcv`.
+We may be interested in fitting a model that lets the spatial patterning differ by species. These kinds of models can be expressed using spatially varying coefficients. Note that we use `spatial = off` because this represents a global spatial intercept---turning this off is akin to using a `-1` of `0` in a main formula including a factor. Both species take their spatial fields from the `spatial_varying` field here.
 
-```{r eval=FALSE}
+```{r}
 fit <- sdmTMB(
-  formula = observed ~ s(year, by = "species"),
+  observed ~ fyear * species,
   data = sim_dat,
+  mesh = mesh,
+  family = gaussian(),
+  spatial = "off",
   time = "year",
   spatiotemporal = "iid",
-  mesh = spde,
-  family = gaussian()
+  spatial_varying = ~ 0 + factor(species)
 )
+fit
 ```
 
-### Model 2: species-specific fields
+Here we had some convergence issues: one of the spatially varying SDs has collapsed to zero.
 
-We may be interested in fitting a model that lets the spatial patterning differ by species. These kinds of models can be expressed using spatially varying coefficients. Note that we use `spatial = off` because this represents a global spatial intercept -- turning this off is akin to using a `-1` of `0` in a main formula including a factor.
+You'll notice that there are two rows of entries for `sigma_Z`:
 
 ```{r}
-fit <- sdmTMB(
-  observed ~ 0 + fyear, # shared year effects
-  data = sim_dat,
-  mesh = spde,
-  family = gaussian(),
-  spatial = "off",
-  time = "year",
-  spatiotemporal = "iid",
-  spatial_varying = ~ 0 + as.factor(species)
-)
-
-fit
+tidy(fit, "ran_pars")
 ```
 
-This model also appears to be problematic. If we look at the estimated parameters, you'll notice that there are two rows of entries for `sigma_Z`. This means that our model is trying to estimate separate species-specific variance terms for the species-specific spatial fields (say *that* 10 times fast!).
+This means that our model is trying to estimate separate species-specific variance terms for the species-specific spatial fields (say *that* 10 times fast!). It's struggling.
+
+If we wanted to estimate species-specific spatial fields with a single shared variance (meaning the net magnitude of the peaks and valleys in the fields were similar but the wiggles themselves were species specific), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`. Any shared factor levels in the `map` are gathered to have 'mirrored' or shared parameter values. We assign these to `ln_tau_Z` because, internally, this is the parameter that gets converted into the spatially-varying field variances (the SDs of those fields are `sigma_Z`).
+
+This case is pretty simple, but for more complicated cases we could figure out the structure of our needed `map` vector as follows:
 
 ```{r}
-fit$sd_report
+colnames(model.matrix(~ 0 + factor(species), data = sim_dat))
 ```
 
-If we wanted to estimate species-specific spatial fields with a single shared variance (meaning the net magnitude of the peaks and valleys in the fields were similar), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`
+So, we need a vector of length two with shared factor values:
 
 ```{r}
 map_list <- list(ln_tau_Z = factor(c(1, 1)))
+```
 
+Then, we can use our map list to share the spatially varying coefficient SDs:
+
+```{r}
 fit <- sdmTMB(
-  observed ~ 0 + fyear, # shared year effects
+  observed ~ fyear * factor(species),
   data = sim_dat,
-  mesh = spde,
+  mesh = mesh,
   family = gaussian(),
   spatial = "off",
   time = "year",
   spatiotemporal = "iid",
-  spatial_varying = ~ 0 + as.factor(species),
+  spatial_varying = ~ 0 + factor(species),
   control = sdmTMBcontrol(map = map_list)
 )
-
 fit
 ```
 
+Notice the spatially varying coefficient SD is now shared.
+
 ### Model 3: species-specific spatiotemporal fields
 
-In all of the examples above, spatiotemporal fields are included -- but shared among species. As a last example, we can extend the above approaches to set up a model that includes spatiotemporal fields unique to each species.
+In all of the examples above, spatiotemporal fields are included, but shared among species. As another example, we can extend the above approaches to set up a model that includes spatiotemporal fields unique to each species.
 
-One approach to including separate spatiotemporal fields by species is creating a new variable that is a concatenation of species and year (or any given time step factor). For example, we can then implement this form of species-specific spatiotemporal variation by changing the `time` argument to be `time = "species_year"`
+One approach to including separate spatiotemporal fields by species is creating a new variable that is a concatenation of species and year (or any given time step factor). For example, we can then implement this form of species-specific spatiotemporal variation by changing the `time` argument to be `time = "species_year"`.
 
 ```{r}
-# This needs to be a factor
-sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
-
+sim_dat$species_year <- factor(paste(sim_dat$species, sim_dat$year))
 map_list <- list(ln_tau_Z = factor(c(1, 1)))
-
 fit <- sdmTMB(
-  observed ~ 0 + fyear, # shared year effects
+  observed ~ fyear * factor(species),
   data = sim_dat,
-  mesh = spde,
+  mesh = mesh,
   family = gaussian(),
   spatial = "off",
   time = "species_year",
   spatiotemporal = "iid",
-  spatial_varying = ~ 0 + as.factor(species),
+  spatial_varying = ~ 0 + factor(species),
   control = sdmTMBcontrol(map = map_list)
 )
-
 fit
-
-print(sanity(fit))
 ```
 
-This is getting a little better -- there is a gradient warning, but the Hessian appears to be invertible (`hessian_ok == TRUE)`, which is critical to convergence.
-
-### Putting it all together
+### Model 4: hack species into the time element for spatial models
 
-Our spatiotemporal model above came close to passing all sanity() checks, and could be modified to also include the species-specific year effects that we started with
+If we only wanted to fit a spatial model but had several species (or other groups), one approach is to pretend our species (or other group) is the time element.
 
 ```{r}
-fit <- sdmTMB(
-  observed ~ 0 + fyear * species, # shared year effects
+sim_dat$numeric_species <- as.numeric(factor(sim_dat$species)) # needs to be numeric
+fit_fake_time <- sdmTMB(
+  observed ~ 0 + factor(species),
   data = sim_dat,
-  mesh = spde,
+  mesh = mesh,
   family = gaussian(),
   spatial = "off",
-  time = "species_year",
-  spatiotemporal = "iid",
-  spatial_varying = ~ 0 + as.factor(species),
-  control = sdmTMBcontrol(map = map_list)
+  time = "numeric_species", #< hack
+  spatiotemporal = "iid" #< 'AR1' or 'RW' probably wouldn't make sense here
 )
+fit_fake_time
+```
 
+This is just a convenience though. We could instead do the same thing using the `spatial_varying` argument making sure to 'map' the field variances to be shared to match the above:
 
-fit
+```{r}
+fit_svc <- sdmTMB(
+  observed ~ 0 + factor(species),
+  data = sim_dat,
+  mesh = mesh,
+  family = gaussian(),
+  spatial = "off",
+  spatial_varying = ~ 0 + factor(species),
+  control = sdmTMBcontrol(map = list(ln_tau_Z = factor(c(1, 1))))
+)
+fit_svc
+```
 
-print(sanity(fit))
+We can prove they're identical:
 
-fit$sd_report
+```{r}
+logLik(fit_fake_time)
+logLik(fit_svc)
 ```
 
+### Putting it all together
+
 These examples illustrate a number of ways that species-specific effects can be included in `sdmTMB` models, and can be extended to other categories/groups/cohorts within a species for which one wants to control the amount of information shared between groups (e.g., age-, size-, or stage-specific estimates). A brief summary of these approaches can be summarized as:
 
 ```{r echo=FALSE}
-desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year by species interactions, or smooths", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
+desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year-by-species interactions or smooths", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
 knitr::kable(desc)
 ```
 
 ### Further extensions
 
-Including spatiotemporal effects like species-year combinations in the `time` argument is relatively straightforward, but these models could get more complicated if we wanted to also include other grouping factors (such as cohort effects for each species). There may be some identifiability issues that would have to be resolved, but even if the `time` argument is being used, it's possible to add other factors through the `spatial_varying` formula. As an example, if our dataframe included a cohort variable for each species, we could first create a new factor that is a concatenation of species and cohort
+Including spatiotemporal effects like species-year combinations in the `time` argument is relatively straightforward, but these models could get more complicated if we wanted to also include other grouping factors (such as cohort effects for each species). There may be some identifiability issues that would have to be resolved, but even if the `time` argument is being used, it's possible to add other factors through the `spatial_varying` formula. As an example, if our dataframe included a cohort variable for each species, we could first create a new factor that is a concatenation of species and cohort.
 
-```{r eval=FALSE}
+```{r eval=TRUE}
+sim_dat$cohort <- rep(seq_len(10), 200)
 sim_dat$species_cohort <- as.factor(paste(sim_dat$species, sim_dat$cohort))
 ```
 
-If there were \~ 10 levels of `species_cohort`, and we were including our original spatial effects in the `spatial_varying` term, we could create a new map list that would indicate we wanted all species-cohort fields to have the same spatial variance as below.
+If there were 10 levels of `species_cohort`, and we were including our original spatial effects in the `spatial_varying` term, we could create a new map list that would indicate we wanted all species-cohort fields to have the same spatial variance as below.
+
+Figure out how our spatial-varying formula will be parsed:
+
+```{r, eval=TRUE}
+colnames(model.matrix(~0 + as.factor(species) + species_cohort, data = sim_dat))
+```
 
 ```{r eval=FALSE}
-sim_dat$species_year <- as.factor(paste(sim_dat$species, sim_dat$year))
+sim_dat$species_year <- factor(paste(sim_dat$species, sim_dat$year))
 
 # The first 2 elements correspond to species, the remaining 10 to species-cohort
 map_list <- list(ln_tau_Z = factor(c(1, 1, rep(2, 10))))
 
 fit <- sdmTMB(
-  observed ~ 0 + fyear * species,
+  observed ~ 0 + fyear * species_cohort,
   data = sim_dat,
-  mesh = spde,
+  mesh = mesh,
   family = gaussian(),
   spatial = "off",
   time = "species_year",
   spatiotemporal = "iid",
   spatial_varying = ~ 0 + as.factor(species) + species_cohort,
-  control = sdmTMBcontrol(map = map_list)
+  do_fit = FALSE
+  # control = sdmTMBcontrol(map = map_list)
 )
 ```

From b94dc68d4b2ea5967b842e36be83e780a145627d Mon Sep 17 00:00:00 2001
From: Sean Anderson <sean@seananderson.ca>
Date: Tue, 3 Dec 2024 16:17:08 -0800
Subject: [PATCH 5/7] More edits

---
 vignettes/web_only/multispecies.Rmd | 40 +++--------------------------
 1 file changed, 3 insertions(+), 37 deletions(-)

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
index 6154d8538..aa2108fdb 100644
--- a/vignettes/web_only/multispecies.Rmd
+++ b/vignettes/web_only/multispecies.Rmd
@@ -66,7 +66,7 @@ sim_dat_B <- sdmTMB_simulate(
   mesh = mesh,
   range = 0.2,
   family = gaussian(),
-  seed = 42,
+  seed = 43,
   sigma_O = 0.3,
   phi = 0.1,
   sigma_E = 0.3,
@@ -118,8 +118,6 @@ fit <- sdmTMB(
 fit
 ```
 
-Here we had some convergence issues: one of the spatially varying SDs has collapsed to zero.
-
 You'll notice that there are two rows of entries for `sigma_Z`:
 
 ```{r}
@@ -229,43 +227,11 @@ logLik(fit_svc)
 These examples illustrate a number of ways that species-specific effects can be included in `sdmTMB` models, and can be extended to other categories/groups/cohorts within a species for which one wants to control the amount of information shared between groups (e.g., age-, size-, or stage-specific estimates). A brief summary of these approaches can be summarized as:
 
 ```{r echo=FALSE}
-desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year-by-species interactions or smooths", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
+desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year-by-species interactions or smooths by year", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
 knitr::kable(desc)
 ```
 
 ### Further extensions
 
-Including spatiotemporal effects like species-year combinations in the `time` argument is relatively straightforward, but these models could get more complicated if we wanted to also include other grouping factors (such as cohort effects for each species). There may be some identifiability issues that would have to be resolved, but even if the `time` argument is being used, it's possible to add other factors through the `spatial_varying` formula. As an example, if our dataframe included a cohort variable for each species, we could first create a new factor that is a concatenation of species and cohort.
-
-```{r eval=TRUE}
-sim_dat$cohort <- rep(seq_len(10), 200)
-sim_dat$species_cohort <- as.factor(paste(sim_dat$species, sim_dat$cohort))
-```
-
-If there were 10 levels of `species_cohort`, and we were including our original spatial effects in the `spatial_varying` term, we could create a new map list that would indicate we wanted all species-cohort fields to have the same spatial variance as below.
-
-Figure out how our spatial-varying formula will be parsed:
-
-```{r, eval=TRUE}
-colnames(model.matrix(~0 + as.factor(species) + species_cohort, data = sim_dat))
-```
-
-```{r eval=FALSE}
-sim_dat$species_year <- factor(paste(sim_dat$species, sim_dat$year))
-
-# The first 2 elements correspond to species, the remaining 10 to species-cohort
-map_list <- list(ln_tau_Z = factor(c(1, 1, rep(2, 10))))
+As long as you're willing to treat spatiotemporal and group-level fields (e.g., for different species or age cohorts) as independent, sdmTMB can be used to fit models to these data. For example, this allows sdmTMB to be used for standardization of age or length composition data as in [Thorson and Haltuch (2018) CJFAS](https://doi.org/10.1139/cjfas-2018-0015). The approach is to similar to the above and we plan to write a separate vignette on the topic.
 
-fit <- sdmTMB(
-  observed ~ 0 + fyear * species_cohort,
-  data = sim_dat,
-  mesh = mesh,
-  family = gaussian(),
-  spatial = "off",
-  time = "species_year",
-  spatiotemporal = "iid",
-  spatial_varying = ~ 0 + as.factor(species) + species_cohort,
-  do_fit = FALSE
-  # control = sdmTMBcontrol(map = map_list)
-)
-```

From 0f2ee5a181db0b2a9c20be452ecaefe0936ae8b7 Mon Sep 17 00:00:00 2001
From: Sean Anderson <sean@seananderson.ca>
Date: Tue, 3 Dec 2024 16:21:22 -0800
Subject: [PATCH 6/7] Edits to article

---
 vignettes/web_only/multispecies.Rmd | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
index aa2108fdb..0be19fbbd 100644
--- a/vignettes/web_only/multispecies.Rmd
+++ b/vignettes/web_only/multispecies.Rmd
@@ -35,7 +35,7 @@ library(sdmTMB)
 
 For some applications, we might be interested in fitting a model that includes multiple responses such as 2+ species, or multiple size or age classes within a species. The most important step in fitting these models is understanding which parameters are shared, and which parameters are species-specific.
 
-Below, we illustrate a series of models that may be used to answer a few different common questions. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial standard deviations (`sigma_O`) as well as different year effects.
+Below, we illustrate a series of models. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial standard deviations (`sigma_O`) as well as different year effects.
 
 ```{r sim_dat}
 set.seed(1)
@@ -77,7 +77,7 @@ sim_dat <- rbind(sim_dat_A, sim_dat_B)
 sim_dat$fyear <- factor(sim_dat$year)
 ```
 
-We'll start by making the mesh across the full dataset.
+We'll start by making an SPDE mesh across the full dataset.
 
 ```{r mesh_fig, fig.asp=1}
 mesh <- make_mesh(sim_dat, c("X", "Y"), cutoff = 0.1)
@@ -86,7 +86,7 @@ plot(mesh)
 
 ### Model 1: species-specific intercepts
 
-As a first model, we can include species-specific year effects. This can be done in a couple ways. One option would be to estimate the `species * year` interaction, letting the year effects for each species be independent. Here, all other parameters (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared. Also, both species share a single spatial random field.
+As a first model, we can include species-specific year effects. This can be done in a couple ways. One option would be to estimate the `species * year` interaction, letting the year effects for each species be independent. Here, all other parameters and random effect values (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared.
 
 ```{r}
 fit <- sdmTMB(
@@ -118,13 +118,13 @@ fit <- sdmTMB(
 fit
 ```
 
-You'll notice that there are two rows of entries for `sigma_Z`:
+You'll notice that there are two rows of entries for `sigma_Z` our spatially varying random field standard deviation:
 
 ```{r}
 tidy(fit, "ran_pars")
 ```
 
-This means that our model is trying to estimate separate species-specific variance terms for the species-specific spatial fields (say *that* 10 times fast!). It's struggling.
+This means that our model is trying to estimate separate species-specific variance terms for the species-specific spatial fields (say *that* 10 times fast!). Here, that matches how we simulated the data. In other contexts, especially if we ran into estimation issues, we might want to share those SDs.
 
 If we wanted to estimate species-specific spatial fields with a single shared variance (meaning the net magnitude of the peaks and valleys in the fields were similar but the wiggles themselves were species specific), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`. Any shared factor levels in the `map` are gathered to have 'mirrored' or shared parameter values. We assign these to `ln_tau_Z` because, internally, this is the parameter that gets converted into the spatially-varying field variances (the SDs of those fields are `sigma_Z`).
 

From ae837b6e81c90f598a19e637512b27336d90bc99 Mon Sep 17 00:00:00 2001
From: Sean Anderson <sean@seananderson.ca>
Date: Tue, 3 Dec 2024 16:24:49 -0800
Subject: [PATCH 7/7] Simplify article setup

---
 vignettes/web_only/multispecies.Rmd | 28 ++++++++++++++++++----------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd
index 0be19fbbd..c61f4e331 100644
--- a/vignettes/web_only/multispecies.Rmd
+++ b/vignettes/web_only/multispecies.Rmd
@@ -11,10 +11,7 @@ vignette: >
 **If the code in this vignette has not been evaluated, a rendered version is available on the [documentation site](https://pbs-assess.github.io/sdmTMB/index.html) under 'Articles'.**
 
 ```{r setup, include = FALSE, cache = FALSE}
-dplyr_installed <- require("dplyr", quietly = TRUE)
-ggplot_installed <- require("ggplot2", quietly = TRUE)
-pkgs <- dplyr_installed && ggplot_installed
-EVAL <- identical(Sys.getenv("NOT_CRAN"), "true") && pkgs
+EVAL <- identical(Sys.getenv("NOT_CRAN"), "true")
 knitr::opts_chunk$set(
   collapse = TRUE,
   warning = FALSE,
@@ -27,9 +24,7 @@ knitr::opts_chunk$set(
 )
 ```
 
-```{r packages, message = FALSE, warning = TRUE}
-library(ggplot2)
-library(dplyr)
+```{r packages, message=FALSE}
 library(sdmTMB)
 ```
 
@@ -227,11 +222,24 @@ logLik(fit_svc)
 These examples illustrate a number of ways that species-specific effects can be included in `sdmTMB` models, and can be extended to other categories/groups/cohorts within a species for which one wants to control the amount of information shared between groups (e.g., age-, size-, or stage-specific estimates). A brief summary of these approaches can be summarized as:
 
 ```{r echo=FALSE}
-desc <- data.frame("Form" = c("Main effects", "Spatial effects", "Spatial effects w/shared variance", "Spatiotemporal effects"), "Implementation" = c("Year-by-species interactions or smooths by year", "Spatially varying coefficients", "Spatially varying coefficients + map argument", "Species-year factor as time variable"))
-knitr::kable(desc)
+desc <- data.frame(
+  Form = c(
+    "Main effects", 
+    "Spatial effects", 
+    "Spatial effects w/shared variance", 
+    "Spatiotemporal effects"),
+  Implementation = c(
+    "Year-by-species interactions or smooths by year", 
+    "Spatially varying coefficients",
+    "Spatially varying coefficients + map argument",
+    "Species-year factor as time variable")
+)
+if (require("knitr", quietly = TRUE)) {
+  knitr::kable(desc)
+} else 
+  print(desc)
 ```
 
 ### Further extensions
 
 As long as you're willing to treat spatiotemporal and group-level fields (e.g., for different species or age cohorts) as independent, sdmTMB can be used to fit models to these data. For example, this allows sdmTMB to be used for standardization of age or length composition data as in [Thorson and Haltuch (2018) CJFAS](https://doi.org/10.1139/cjfas-2018-0015). The approach is to similar to the above and we plan to write a separate vignette on the topic.
-