Add table with MCMC diagnostics Rhat, ESSbulk, and ESStail

doc-sr/03_analysis.Rmd

@@ -5,7 +5,7 @@ cat("# Analysis and Response
 
 The model used to assess this stock was the Integrated Statistical Catch-at-Age Model (`r iscam`). It was tuned to four fishery-independent trawl survey series covering `r start_catch_yr`--2023, a Discard CPUE series as an index of abundance, annual estimates of commercial catch from two fleets (Freezer Trawlers and Shoreside), and age composition data from the two fleets in the commercial fishery and the four surveys. A two-sex, two-fleet base model was selected and implemented in a Bayesian context using Markov Chain Monte Carlo (MCMC) methods. Leading parameters estimated included $R_0$, initial recruitment, $h$, steepness of the stock-recruitment relationship, $\bar{R}$, average recruitment, and $q_k, k=1,2,3,4,5$, catchability of the four surveys and the Discard CPUE index. Selectivity parameters were also estimated for each sex, fleet, and survey.
 
-Parameter estimates and fixed values are given in Table \@ref(tab:param-estimates-table). As in the `r ca`, the natural mortality was fixed at 0.2 for females and 0.35 for males. Selectivity was estimated for all survey indices but fixed for the Discard CPUE ($a\\hat_7 = 9.5$ and $\\gamma_7 = 0.5$ for both sexes in Table \@ref(tab:param-estimates-table)).
+Parameter estimates and fixed values are given in Table \@ref(tab:param-estimates-table). As in the `r ca`, the natural mortality was fixed at 0.2 for females and 0.35 for males. Selectivity was estimated for all survey indices but fixed for the Discard CPUE ($\hat{a}_{7,sex,1} = 9.5$ and $\gamma_{7,sex,1} = 0.5$ for both sexes in Table \@ref(tab:param-estimates-table)).
 
 All estimated parameter estimates were close to those in Table 6 of the `r ca`. The median posterior for $B_0$ decreased from 184.16 in `r ca` to `r f(median(base_model$mcmccalcs$params$sbo), 2)` for this model. The median posterior biomass estimates were also slightly less than the estimates in the `r ca`, so there is almost no scaling effect in the relative biomass and Figure \@ref(fig:fig-base-depletion) looks almost identical to Figure 9 in the `r ca`, other than the two new points for 2022 and 2023. In both cases, relative biomass estimates from 2020 and forward were under the USR $0.4B_0$ reference line.
 
@@ -81,6 +81,8 @@ Figure \@ref(fig:fig-base-trace) shows the traceplots for the leading parameters
 
 The catchability parameters appear to show more correlation with each other in this update model than seen in the `r ca`. A comparison of the pairs plots from this update (Figure \@ref(fig:fig-base-pairs)) and Figure 41 from the `r ca` shows that there is a higher correlation between the $q_1$ and $q_3$ parameters in particular. These are the catchability parameters for the `r qcsss` and the `r hsss`.
 
+The within-chain $\hat{R}$ and Effective sample sizes (ESS) for the main part of the posterior distribution ($ESS_{bulk}$) as well as the tails ($ESS_{tail}$) were computed for each lead parameter. All values are considered acceptable for convergence of the MCMC chain. The criteria are a value of 1 of less for $\hat{R}$ and a value close to the number of posteriors for the ESS values. Table \@ref(tab:mcmc-diagnostics-rhat-table) shows the values calculated.
+
 ")
 ```
 
@@ -91,11 +93,11 @@ The catchability parameters appear to show more correlation with each other in t
 ```{r analysis-and-response-projection-assumptions-en, eval = !fr(), results = 'asis'}
 cat("## Projection assumptions
 
-Projection assumptions for this update were identical to those in the `r ca`. Projected log recruitment deviations in the years 2025-2027 were drawn randomly with replacement from the estimated 2010--2019 deviations (omitting poorly estimated 2020--2023 deviations). All other model parameterizations remained the same as in the pre-projection period.
+Projection assumptions for this update were identical to those in the `r ca`. Projected log recruitment deviations in the years 2025-2027 were drawn randomly with replacement from the estimated 2010--2019 deviations (omitting poorly estimated 2020--2024 deviations). All other model parameterizations remained the same as in the pre-projection period.
 ")
 ```
 
-```{r analysis-and-response-projection-assumptions-fr, eval = fr(), needs_trans = TRUE}
+```{r analysis-and-response-projection-assumptions-fr, eval = fr(), results = "asis", needs_trans = TRUE}
 <<analysis-and-response-projection-assumptions-en>>
 ```
 
@@ -300,6 +302,22 @@ plot_pairs_mcmc(base_model,
 
 \clearpage
 
+```{r mcmc-diagnostics-rhat-table, results = "asis"}
+cap <- paste0("Values calculated to show convergence of the MCMC chain. Numerical indices on Catchability ($q$) parameters for the survey indices are defined as: ",
+paste(base_index_gears$gear, " = ", base_index_gears$gear_name, collapse = ", "),
+".")
+if(fr()){
+cap <- paste0("")
+}
+
+csas_table(rhat_df,
+           format = "latex",
+           align = rep("r", ncol(rhat_df)),
+           cap = cap,
+           longtable = FALSE)
+```
+
+
 ```{r indicators-of-stock-status-en, eval = !fr(), results = 'asis'}
 cat("## Indicators of Stock Status
 

doc-sr/04_conclusions.Rmd

@@ -38,6 +38,6 @@ Some observations from the decision tables:
 
 9. A probability of 0.5 of biomass increasing year-to-year for each year from `r end_yr + 1` to `r end_yr + base_model$proj$num.projyrs` occurs between 7 and 8 kt of constant catch in each year (Table \@ref(tab:decision-table-decreasing-biomass)).
 
-This stock should be reviewed once again in three years, or September/October 2027 as a Science Response. In addition to new catch data, that update model must have new age data included for the three synoptic surveys and both fleets in the commercial fishery for the time period 2022-2026.
+This stock should be reviewed once again in three years, in September or October 2027 as a Science Response. In addition to new catch and survey index data, that update model must have new age data included for the three synoptic surveys and both fleets in the commercial fishery for the time period 2022-2026.
 ")
 ```

doc-sr/index.Rmd

@@ -712,6 +712,52 @@ ct_ss_2023 <- ct_df |>
 perc_dec_ss_2021_2022 <- f(100 - ct_ss_2022 / ct_ss_2021 * 100, 1)
 perc_dec_ss_2022_2023 <- f(100 - ct_ss_2023 / ct_ss_2022 * 100, 1)
 
+# create and filter an MCMC diagnostic data frame with one value column
+# (second column parameter name) for rstan Rhat and ESS values
+create_mcmc_diag <- function(df,
+                             fun = rstan::Rhat,
+                             digits = 3,
+                             value_nm = "Rhat"){
+  
+  d <- df |>
+  map_dbl(~{fun(.x)}) |>
+  enframe() |> 
+  filter(!is.na(value)) |> 
+  filter(!grepl("^sel", name)) |> 
+  filter(!grepl("msy", name)) |> 
+  filter(!grepl("ssb|beta", name)) |> 
+  filter(!grepl("^bo$", name)) |> 
+  filter(!grepl("^so$", name)) |> 
+  mutate(value = f(value, digits))
+  names(d) <- c("param", value_nm)
+  
+  d
+}
+# Rhat values for MCMC diagnostics as requested by Mazur review Oct. 2, 2024
+rhat_df <- base_model$mcmc$params |>
+  create_mcmc_diag()
+ess_bulk_df <- base_model$mcmc$params |>
+  create_mcmc_diag(fun = rstan::ess_bulk,
+                   digits = 0,
+                   value_nm = "ESS bulk")
+ess_tail_df <- base_model$mcmc$params |>
+  create_mcmc_diag(fun = rstan::ess_tail,
+                   digits = 0,
+                   value_nm = "ESS tail")
+
+lu_param_nm <- tibble(param = c("ro", "h", "rbar", "rinit", "sbo",
+                                "q1", "q2", "q3", "q4", "q5"),
+                      Parameter = c("$R_0$", "$h$", "$\\bar{R}$", "$\\bar{R}_{init}$",
+                                    "$SB_0$", "$q_1$", "$q_2$", "$q_3$", "$q_4$", "$q_5$"))
+rhat_df <- rhat_df |> 
+  full_join(ess_bulk_df, by = "param") |> 
+  full_join(ess_tail_df, by = "param") |> 
+  full_join(lu_param_nm, by = "param") |> 
+  select(-param) |> 
+  select(Parameter, everything())
+
+names(rhat_df) <- c("Parameter", "$\\hat{R}$", "$ESS_{bulk}$", "$ESS_{tail}$")
+
 ca <- "2022 stock assessment"
 ```
 

doc/bib/refs.bib

@@ -1872,7 +1872,7 @@ @article{arf2022
 @Article{         dfo2023arrowtoothSAR,
   author        = {DFO},
   year          = {2023},
-  title         = {{Arrowtooth Flounder} ({\emph{Atheresthes stomias}}) Stock Assessment for {{British Columbia} in 2021},
+  title         = {{Arrowtooth Flounder} ({\emph{Atheresthes stomias}}) Stock Assessment for British Columbia in 2021},
   journal       = {DFO Can. Sci. Advis. Sec. Sci. Advis. Rep.},
   volume        = {2023/042}
 }