From 9c7bd5136aa71db0aab18c308708925ad1afb317 Mon Sep 17 00:00:00 2001 From: JaneSullivan-NOAA Date: Wed, 4 Dec 2024 00:51:11 +0000 Subject: [PATCH] =?UTF-8?q?Deploying=20to=20gh-pages=20from=20@=20afsc-ass?= =?UTF-8?q?essments/rema@14cf71d9ff83902d1e83f1d017f8f9e9584f2dc5=20?= =?UTF-8?q?=F0=9F=9A=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- articles/ex2_cpue.html | 2 +- articles/rema_equations.html | 2 +- index.html | 2 +- pkgdown.yml | 2 +- search.json | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/articles/ex2_cpue.html b/articles/ex2_cpue.html index 34e46e9..4cab23b 100644 --- a/articles/ex2_cpue.html +++ b/articles/ex2_cpue.html @@ -209,7 +209,7 @@

Mode extra_cpue_cv = list(assumption = 'extra_cv')) m2 <- fit_rema(input2) -#> Model runtime: 0.2 seconds +#> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.58e-05 #> Max gradient parameter: log_q diff --git a/articles/rema_equations.html b/articles/rema_equations.html index 76228f1..7442094 100644 --- a/articles/rema_equations.html +++ b/articles/rema_equations.html @@ -222,7 +222,7 @@

The ADMB version of REMAln(Ît)=ln(q)+ln(B̂)ln(\hat{I}_{t})=ln(q)+ln(\hat{B})] occurred outside rather than inside the SEPARABLE_FUNCTION, -and as a results, violated this rule and affecting the accuracy of the +and as a result, violated this rule and affecting the accuracy of the Laplace approximation.

We explored this error in a simplified example of the REMA model in both ADMB and TMB. We were able to diff --git a/index.html b/index.html index 81a9dcb..33c3cac 100644 --- a/index.html +++ b/index.html @@ -71,7 +71,7 @@

A generalized Random Effects model for fitting survey biomass estimates with the options of including Multiple survey strata and an Additional survey index

Sullivan, J., C. Monnahan, P. Hulson, J. Ianelli, J. Thorson, and A. Havron. 2022a. REMA: a consensus version of the random effects model for ABC apportionment and Tier 4/5 assessments. Plan Team Report, Joint Groundfish Plan Teams, North Pacific Fishery Management Council. 605 W 4th Ave, Suite 306 Anchorage, AK 99501. Available through the Oct 2022 Joint GPT e-Agenda.

-

Living documentation for REMA model equations is available on the REMA website:

+

Living documentation for REMA model equations is available on the REMA website: https://afsc-assessments.github.io/rema/articles/rema_equations.html

Background diff --git a/pkgdown.yml b/pkgdown.yml index bea721f..0da1bc1 100644 --- a/pkgdown.yml +++ b/pkgdown.yml @@ -7,7 +7,7 @@ articles: ex3_zeros: ex3_zeros.html ex4_model_validation: ex4_model_validation.html rema_equations: rema_equations.html -last_built: 2024-12-04T00:29Z +last_built: 2024-12-04T00:50Z urls: reference: https://afsc-assessments.github.io/rema/reference article: https://afsc-assessments.github.io/rema/articles diff --git a/search.json b/search.json index e56d6af..6017f68 100644 --- a/search.json +++ b/search.json @@ -1 +1 @@ -[{"path":"https://afsc-assessments.github.io/rema/articles/ex1_basics.html","id":"the-rema-workflow","dir":"Articles","previous_headings":"","what":"The rema workflow:","title":"REMA basics","text":"Load rema data. user can read biomass abundance index data file (e.g. .csv file), can use rwout.rep report file ADMB version RE model using read_admb_re(). Specify model structure assumptions using prepare_rema_input(). function allows users quickly transition single two survey model, specify alternative process error structures, add likelihood penalties priors parameters, evaluate alternative assumptions zero biomass observations. Fit specified REMA model using fit_rema() determine whether model met basic convergence criteria (e.g., Hessian positive definite, maximum gradient component approximately equal zero). Extract rema model output clean, consistently formatted data frames using tidy_rema(). user can visualize model output using plot_rema(), quickly format tables report. Compare alternative REMA models conduct model selection using compare_rema_models(). Output function includes table Akaike Information Criteria (AIC) appropriate, figures, tidied data frames. function also accepts model output ADMB version RE model easy comparison past models. Taken together, functions allow R users quickly fit interrogate suite simple statistical models TMB without needing software-specific training expertise.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex1_basics.html","id":"example-1-univariate-random-effects-re-model-with-a-single-survey-and-stratum","dir":"Articles","previous_headings":"","what":"Example 1: Univariate random effects (RE) model with a single survey and stratum","title":"REMA basics","text":"example uses Aleutian Islands shortraker (aisr_rwout.rep) fits NMFS bottom trawl survey estimates. Read rwout.rep, custom report file generated ADMB version random effects model. Prepare REMA model inputs using admb_re data object. input$data, input$par, input$map, input$random used TMB fit REMA model. objects used functions process model results conduct residual analyses. Data alternatively entered prepare_rema_input() using biomass_dat argument: Fit REMA model Get tidied model output plot results Note 95% confidence intervals observations (.e., obs_lci obs_uci output$biomass_by_strata based assumption normality log-space; therefore, asymmetric arithmetic scale. Compare ADMB RE model results","code":"# ?read_admb_re admb_re <- read_admb_re(filename = 'aisr_rwout.rep', # optional label for the single biomass survey stratum biomass_strata_names = 'Aleutians Islands', model_name = 'ADMB: AI shortraker') names(admb_re) #> [1] \"biomass_dat\" \"cpue_dat\" \"model_yrs\" #> [4] \"init_log_biomass_pred\" \"admb_re_results\" kable(admb_re$biomass_dat) # ?prepare_rema_input input <- prepare_rema_input(model_name = 'TMB: AI shortraker', admb_re = admb_re) names(input) #> [1] \"data\" \"par\" \"map\" \"random\" \"model_name\" #> [6] \"osa\" \"biomass_dat\" # (not run) input <- prepare_rema_input(model_name = 'tmb_rema_aisr', biomass_dat = admb_re$biomass_dat) # ?fit_rema m <- fit_rema(input) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.90e-11 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model # ?tidy_rema output <- tidy_rema(rema_model = m) kable(output$parameter_estimates) # estimated fixed effects parameters kable(head(output$biomass_by_strata, 5)) # data.frame of predicted and observed biomass by stratum # ?plot_rema plots <- plot_rema(tidy_rema = output, biomass_ylab = 'Biomass (t)') # optional y-axis label plots$biomass_by_strata # ?compare_rema_models compare <- compare_rema_models(rema_models = list(m), admb_re = admb_re, biomass_ylab = 'Biomass (t)') compare$plots$biomass_by_strata"},{"path":"https://afsc-assessments.github.io/rema/articles/ex1_basics.html","id":"example-2-multivariate-random-effects-model-rem-with-a-single-survey-and-multiple-strata","dir":"Articles","previous_headings":"","what":"Example 2: Multivariate random effects model (REM) with a single survey and multiple strata","title":"REMA basics","text":"example uses Bering Sea Aleutian Islands shortspine thornyhead (bsaisst_rwout.rep) fits NMFS EBS slope AI bottom trawl survey estimates. original ADMB model single, shared process error (PE) three strata. example demonstrates define alternative PE structures, u perform model selection using AIC. can use ggplot2 functions modify formatting plots: One primary uses RE model apportioning catch management area. tidy_rema() function output includes table proportions predicted biomass strata: can compare TMB model results original ADMB model. Note different confidence intervals ADMB TMB versions. ADMB code uses Marlow method sum variances biomass log-space, whereas TMB version applies standard Delta method. can easily fit alternative REMA model strata-specific PE parameters, compare two models using AIC. case see single PE model lowest AIC value. Note ADMB models currently compared REMA models using AIC, therefore omit admb_re = admb_re argument compare_rema_models() function call :","code":"admb_re <- read_admb_re(filename = 'bsaisst_rwout.rep', biomass_strata_names = c('AI survey', 'EBS slope survey', 'S. Bering Sea (AI survey)'), model_name = 'ADMB: single PE') # the original ADMB model shared PE across all strata input1 <- prepare_rema_input(model_name = 'TMB: single PE', admb_re = admb_re, PE_options = list(pointer_PE_biomass = c(1, 1, 1))) m1 <- fit_rema(input1) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 9.81e-10 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model output <- tidy_rema(rema_model = m1) kable(output$parameter_estimates) # estimated fixed effects parameters kable(head(output$biomass_by_strata, 5)) # data.frame of predicted and observed biomass by stratum kable(head(output$total_predicted_biomass, 5)) # combined/total predicted biomass plots <- plot_rema(tidy_rema = output, biomass_ylab = 'Biomass (t)') plots$biomass_by_strata + facet_wrap(~strata, ncol = 1, scales = 'free_y') kable(tail(output$proportion_biomass_by_strata, 3)) # figure: # plots$proportion_biomass_by_strata compare <- compare_rema_models(rema_models = list(m1), admb_re = admb_re, biomass_ylab = 'Biomass (t)') # Note different confidence intervals between the ADMB version (Marlow method) and the TMB version (Delta method) compare$plots$biomass_by_strata + facet_wrap(~strata, ncol = 1, scales = 'free_y') compare$plots$total_predicted_biomass + ggplot2::ggtitle('BSAI Shortspine thornyhead') # REMA defaults to strata-specific parameters, which could be explicitly defined # as: PE_options = list(pointer_PE_biomass = c(1, 2, 3)) input2 <- prepare_rema_input(model_name = 'TMB: strata-specific PE', admb_re = admb_re) m2 <- fit_rema(input2) #> Model runtime: 0.3 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 2.31e-07 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model compare <- compare_rema_models(rema_models = list(m1, m2), biomass_ylab = 'Biomass (t)') compare$plots$biomass_by_strata + facet_wrap(~strata, ncol = 1, scales = 'free_y') kable(compare$aic) kable(compare$output$parameter_estimates)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-1-fit-to-the-biomass-and-cpue-survey-through-the-estimation-of-a-scaling-parameter-q","dir":"Articles","previous_headings":"","what":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","title":"Fitting to an additional CPUE survey","text":"","code":"# read in the data using the report file from the original ADMB model admb_re <- read_admb_re(filename = 'goasr_rwout.rep', biomass_strata_names = c('CGOA', 'EGOA', 'WGOA'), cpue_strata_names = c('CGOA', 'EGOA', 'WGOA'), model_name = 'M0: ADMB status quo (invalid)') input <- prepare_rema_input(model_name = 'M1: TMB status quo', # fit to both biomass and CPUE survey data multi_survey = 1, admb_re = admb_re, # the RPW index is summable across strata sum_cpue_index = TRUE, # likelihood weight wt_cpue = 0.5, # one process error parameters (log_PE) estimated PE_options = list(pointer_PE_biomass = c(1, 1, 1)), # three scaling parameters (log_q) estimated, indexed as # follows for each biomass survey stratum: q_options = list(pointer_biomass_cpue_strata = c(1, 2, 3))) m1 <- fit_rema(input) #> Model runtime: 0.2 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 4.12e-04 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model output <- tidy_rema(m1) output$parameter_estimates #> model_name parameter estimate std_err lci #> 1 M1: TMB status quo process_error 0.1696964 0.03914612 0.1079729 #> 2 M1: TMB status quo scaling_parameter_q 0.4123880 0.04474590 0.3333857 #> 3 M1: TMB status quo scaling_parameter_q 1.2127809 0.10874662 1.0173198 #> 4 M1: TMB status quo scaling_parameter_q 2.1164194 0.30469400 1.5960887 #> uci #> 1 0.2667047 #> 2 0.5101115 #> 3 1.4457965 #> 4 2.8063798 plots <- plot_rema(output, biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') cowplot::plot_grid(plots$biomass_by_strata, plots$cpue_by_strata, ncol = 1)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"a-comparison-of-the-admb-and-tmb-models","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"A comparison of the ADMB and TMB models","title":"Fitting to an additional CPUE survey","text":"compare fits status quo ADMB model identically structured correctly-specified TMB model, see large differences model predictions. particular PE variance different, ADMB model producing PE seven times PE TMB model: M0: ADMB status quo (invalid) PE = 0.0233 M1: TMB status quo PE = 0.1697","code":"compare <- compare_rema_models(rema_models = list(m1), admb_re = admb_re, biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') compare$plots$biomass_by_strata + theme(legend.position = 'top')"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-2-estimating-additional-observation-error-in-the-two-surveys","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Model 2: Estimating additional observation error in the two surveys","title":"Fitting to an additional CPUE survey","text":"outcome TMB model problematic expect shortraker rockfish, notably long-lived species, exhibit extreme inter-annual variability. part, finding attributed relatively low observation error survey data. following section explore alternative model estimates additional biomass CPUE survey observation error. use AIC conduct model selection. following figure shows fits data ; however, used tidy_extra_cv() plot_extra_cv() functions obtain total 95% confidence intervals survey biomass CPUE observations (.e., fixed + estimated observation error). Note 95% confidence intervals observations (.e., obs_lci/obs_uci tot_obs_lci/tot_obs_uci output$biomass_by_strata output$cpue_by_strata based assumption normality log-space; therefore, asymmetric arithmetic scale. figure , error bars around survey observations show 95% confidence intervals based assumed design-based estimates (bold whiskers) 95% confidence intervals based total observation error (design-based estimates CV + additional estimated CV; whiskers).","code":"input2 <- prepare_rema_input(model_name = 'M2: TMB extra obs error', # fit to both biomass and CPUE survey data multi_survey = 1, admb_re = admb_re, # the RPW index is summable across strata sum_cpue_index = TRUE, # likelihood weight wt_cpue = 0.5, # one process error parameters (log_PE) estimated PE_options = list(pointer_PE_biomass = c(1, 1, 1)), # three scaling parameters (log_q) estimated, indexed as # follows for each biomass survey stratum: q_options = list(pointer_biomass_cpue_strata = c(1, 2, 3)), # estimate additional obs error on the biomass # survey (one additional CV shared across all 3 # strata) extra_biomass_cv = list(assumption = 'extra_cv'), # estimate additional obs error on the CPUE # survey (one additional CV shared across all 3 # strata) extra_cpue_cv = list(assumption = 'extra_cv')) m2 <- fit_rema(input2) #> Model runtime: 0.2 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.58e-05 #> Max gradient parameter: log_q #> TMB:sdreport() was performed successfully for this model out2 <- tidy_rema(m2) # adds new columns with total observation error and 95% confidence intervals to # the biomass_by_strata and cpue_by_strata out2 <- tidy_extra_cv(out2) cvplots <- plot_extra_cv(out2) cowplot::plot_grid(cvplots$biomass_by_strata, cvplots$cpue_by_strata, ncol = 1)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-comparison","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Model comparison","title":"Fitting to an additional CPUE survey","text":"comparison three models shows estimation additional observation error substantially reduces PE variance, resulting smoother trajectory biomass predictions.","code":"compare <- compare_rema_models(rema_models = list(m1, m2), admb_re = admb_re, biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') compare$plots$biomass_by_strata + theme(legend.position = 'top')"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-selection","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Model selection","title":"Fitting to an additional CPUE survey","text":"Next, can compare fits biomass CPUE survey indices alternative TMB models remove admb_re object comparison function demonstrated . examination parameter estimates shows M2 (extra observation error) estimates much lower PE, approximately 50% M1 PE. Using AIC, find M2 superior fit M1 statistical support additional parameters (note example intended illustrative model fitting model selection process exhaustive).","code":"compare <- compare_rema_models(rema_models = list(m1, m2), biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') cowplot::plot_grid(compare$plots$biomass_by_strata + theme(legend.position = 'top'), compare$plots$cpue_by_strata + theme(legend.position = 'none'), ncol = 1, rel_widths = c(0.65, 0.35)) kable(compare$output$parameter_estimates) kable(compare$aic)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"apportionment-results","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Apportionment results","title":"Fitting to an additional CPUE survey","text":"proportion prediction biomass strata (management area example), often used apportionment Acceptable Biological Catch estimates. Results compare_rema_models() function include table proportioned biomass management area purpose.","code":"compare$output$proportion_biomass_by_strata %>% filter(year == max(year)) %>% kable() # compare$plots$proportion_biomass_by_strata # optional figure"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"model-1-zeros-as-nas","dir":"Articles","previous_headings":"","what":"Model 1: Zeros as NAs","title":"Strategies for handling zero biomass observations","text":"","code":"nonsst <- read.csv('ebsshelf_orox_biomass.csv') input1 <- prepare_rema_input(model_name = 'M1: Zeros as NAs', biomass_dat = nonsst, zeros = list(assumption = 'NA')) m1 <- fit_rema(input1) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 4.23e-10 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"model-2-add-a-small-constant","dir":"Articles","previous_headings":"","what":"Model 2: Add a small constant","title":"Strategies for handling zero biomass observations","text":"","code":"input2 <- prepare_rema_input(model_name = 'M2: Small constant=0.01, CV=1.5', biomass_dat = nonsst, zeros = list(assumption = 'small_constant', # values: 1) small constant, 2) assumed CV options_small_constant = c(0.01, 1.5))) m2 <- fit_rema(input2) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.33e-07 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"model-3-tweedie-distribution","dir":"Articles","previous_headings":"","what":"Model 3: Tweedie distribution","title":"Strategies for handling zero biomass observations","text":"","code":"input3 <- prepare_rema_input(model_name = 'M3: Tweedie', biomass_dat = nonsst, zeros = list(assumption = 'tweedie')) m3 <- fit_rema(input3) #> Model runtime: 0.2 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.72e-10 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"compare-model-results","dir":"Articles","previous_headings":"","what":"Compare model results","title":"Strategies for handling zero biomass observations","text":"Models 1 (Zeros NAs) 3 (Tweedie) similar results. Tweedie appears perform well case, can slow run suffer convergence issues, especially observation errors relatively small. Tweedie alternative considered experimental issues can resolved.","code":"compare <- compare_rema_models(rema_models = list(m1, m2, m3)) compare$plots$biomass_by_strata"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"background","dir":"Articles","previous_headings":"","what":"Background","title":"REMA model validation","text":"Fisheries stock assessments moving towards state-space estimation, boasts range benefits, including separation estimation observation process error, elegant framework handling missing data, high degree flexibility respect model architecture inclusion different data types, potential improved projections predictive skill (Aeberhard et al., 2018). flexibility relative ease fitting state-space models means can increase complexity dimensionality rapidly. easy fit, state-space models often challenging validate, even simple models can suffer estimation issues result biased inference (Auger‐Méthé et al., 2016). North Pacific Fishery Management Council (NPFMC), “random effects model” (REMA) far common state-space model used fishery management (Sullivan et al., 2022). REMA state-space random walk model can customized estimate multiple process errors, fit additional abundance index, estimate additional observation error. REMA used estimate biomass within Tier 5 groundfish Tier 4 crab stock assessments apportion Acceptable Biological Catches (ABCs) management area many stocks. Despite high impact model within North Pacific fishery management process, standard model validation practices REMA. goals: Apply established state-space model validation techniques real life REMA examples. Create REMA model validation guide clear descriptions method motivation use, explanation interpret model validation results, reproducible code run operational REMA models NPFMC. Provide preliminary recommendations REMA users reviewers model validation methods relevant informative REMA. Interested otherwise busy readers encouraged skip “care failed model diagnostics?” “Recommendations” sections. Table 1 Figures 6-8 show examples REMA model failing model validation criteria.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"a-testing-framework-for-rema-model-validation","dir":"Articles","previous_headings":"","what":"A testing framework for REMA model validation","title":"REMA model validation","text":"analysis, testing (1) REMA working expected (.e., code implemented correctly parameters estimated without bias), (2) assumptions made parameterizing estimating REMA valid, including assumptions related random effects error structure. use two example stocks: Aleutian Islands Pacific cod (AI Pcod; Spies et al., 2023): simple, univariate case, fits single time series AI bottom trawl survey estimates one process error Gulf Alaska Thornyhead rockfish (GOA thornyhead; Echave et al., 2022): complex, multi-strata multi-survey case, estimates multiple process errors, scaling parameter longline survey, two additional observation errors GOA bottom trawl longline surveys models appear converge per standard REMA TMB diagnostics (small maximum gradient, invertible Hessian), fit data reasonably, model fits shown Figure 1. stocks span levels complexity NPFMC REMA models; testing REMA across range complexity helps ensure model model assumptions valid variety realistic, management-relevant cases. important note model validation mean model “right” “correct.” Model validation help analyst model selection (except identify models might function properly), ensure model makes “better” predictions. Rather, model validation way ensure model operating expected, without introducing bias violating statistical assumptions upon model based (Auger‐Méthé et al., 2016; Auger‐Méthé et al., 2021), data reasonably come model (Thygesen et al., 2017). key questions aim answer model validation include following: model perform expected, introduce bias? Using simulation self-check, test model coded correctly able recover known parameters. plausible data generated model? Using one-step ahead (OSA) residuals, appropriate residual type state-space models, test model assumptions look trends residuals reveal characteristics dynamics data aren’t adequately captured model. normality assumptions made estimating random effects via Laplace approximation accurate? comparing posterior distribution fixed effects without Laplace approximation, test whether distribution random effects normal thus accuracy Laplace approximation. allows us test accuracy fixed effects estimates uncertainties. parameters unique non-redundant? Checking correlation parameters helps us identify parameters redundant .","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"prepare-data-for-each-stock-here","dir":"Articles","previous_headings":"","what":"Prepare data for each stock here","title":"REMA model validation","text":"First, run model stock. AI Pcod model (pcod_mod) uses single process error describe stock time (one parameter). GOA thornyhead rockfish model (thrn_mod) uses data two surveys (bottom trawl longline survey) describe stock time (six fixed effects parameters: three process errors, one scaling parameter, additional observation error biomass survey, additional observation error longline survey). code , show specify fit models within rema framework plot fits data. Figure 1 shows fits data AI Pcod GOA Thornyhead. particular, GOA Thornyhead fits show trade biomass longline survey data. Figure 1. Fits AI Pcod (top panel; gold) GOA Thornyhead (bottom panels; teal) models.","code":"set.seed(415) # for reproducibility, use a seed # first, get the p-cod data set up pcod_bio_dat <- read.csv(\"ai_pcod_2022_biomass_dat.csv\") pcod_input <- prepare_rema_input(model_name = \"p_cod\", biomass_dat = pcod_bio_dat, # one strata PE_options = list(pointer_PE_biomass = 1) ) # run the model pcod_mod <- fit_rema(pcod_input) # next, get the thornyhead data set up thrn_bio_dat <- read.csv(\"goa_thornyhead_2022_biomass_dat.csv\") thrn_cpue_dat <-read.csv(\"goa_thornyhead_2022_cpue_dat.csv\") thrn_input <- prepare_rema_input(model_name = 'thrnhead_rockfish', multi_survey = TRUE, biomass_dat = thrn_bio_dat, cpue_dat = thrn_cpue_dat, # RPWs are a summable/area-weighted effort index sum_cpue_index = TRUE, # three process error parameters (log_PE) estimated, # indexed as follows for each biomass survey stratum # (shared within an area across depths): PE_options = list(pointer_PE_biomass = c(1, 1, 1, 2, 2, 2, 3, 3, 3)), # scaling parameter options: q_options = list( # longline survey strata (n=3) indexed as follows for the # biomass strata (n=9) pointer_biomass_cpue_strata = c(1, 1, 1, 2, 2, 2, 3, 3, 3), # one scaling parameters (log_q) estimated, shared # over all three LLS strata pointer_q_cpue = c(1, 1, 1)), # estimate extra trawl survey observation error extra_biomass_cv = list(assumption = 'extra_cv'), # estimate extra longline survey observation error extra_cpue_cv = list(assumption = 'extra_cv') ) # run the model thrn_mod <- fit_rema(thrn_input) # tidy output and plot fitted data tidy_pcod <- tidy_rema(pcod_mod) p1 <- plot_rema(tidy_pcod)$biomass_by_strata + ggtitle(label = \"Model Fits to the AI Pcod Data\", subtitle = \"Trawl Survey Biomass Strata\") + geom_ribbon(aes(ymin = pred_lci, ymax = pred_uci), col = 'goldenrod', fill = 'goldenrod', alpha = 0.4) + geom_line() + geom_point(aes(x = year, y = obs)) + geom_errorbar(aes(x = year, ymin = obs_lci, ymax = obs_uci)) p1 ggsave(paste0(\"vignettes/ex4_pcod_fits.png\"), width = 6, height = 4, units = \"in\", bg = \"white\") tidy_thrn <- tidy_rema(thrn_mod) p2 <- plot_rema(tidy_thrn, biomass_ylab = 'Biomass (t)')$biomass_by_strata + ggtitle(label = \"Model Fits to the GOA Thornyhead Data\", subtitle = \"Trawl Survey Biomass Strata\") + geom_ribbon(aes(ymin = pred_lci, ymax = pred_uci), col = \"#21918c\", fill = \"#21918c\", alpha = 0.4) + geom_line() + geom_point(aes(x = year, y = obs)) + geom_errorbar(aes(x = year, ymin = obs_lci, ymax = obs_uci)) p2 p3 <- plot_rema(tidy_thrn, cpue_ylab = \"RPW\")$cpue_by_strata + ggtitle(label = NULL, subtitle = \"Longline Survey Relative Population Weight (RPW) Strata\") + geom_ribbon(aes(ymin = pred_lci, ymax = pred_uci), col = \"#21918c\", fill = \"#21918c\", alpha = 0.4) + geom_line() + geom_point(aes(x = year, y = obs)) + geom_errorbar(aes(x = year, ymin = obs_lci, ymax = obs_uci)) p3 cowplot::plot_grid(p2, p3, rel_heights = c(0.7, 0.3), ncol = 1) ggsave(paste0(\"vignettes/ex4_thrn_fits.png\"), width = 11, height = 10, units = \"in\", bg = \"white\")"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"simulation-self-test-can-we-recover-parameters-without-bias","dir":"Articles","previous_headings":"","what":"1. Simulation self-test: Can we recover parameters without bias?","title":"REMA model validation","text":"Simulation testing ensures model properly coded performs consistently without introducing bias (Auger‐Méthé et al., 2021; Gimenez et al., 2004). First, use REMA model estimate parameters (e.g., process error variance) real data. Next, use estimated parameters “true values” simulate new latent states (random effects) data conditioned states using REMA equations. use REMA model re-estimate model parameters (“recovered parameters”) calculate relative error (RE; .e., (true-estimated values)/true value*100)) parameters simulation replicate (N=500). AI Pcod model fitted bottom trawl data simulated data set biomass trajectory; GOA Thornyhead model fitted bottom trawl longline survey data, simulation data set includes biomass trajectories nine strata longline survey relative population weights (RPWs) three strata. Models fail simulation testing recovered parameters biased (RE ≠ 0), deviate consistently true values used simulate data (span recovered parameters large): can indicate model coding error, non-identifiable, redundant parameters, biased (.e., kind model misspecification). following function runs simulation self-test REMA model. Run self-test AI Pcod GOA Thornyhead examine results: Results simulation self-test 1 500 simulation replicates converge AI Pcod model, suggests high degree model stability. 11 500 simulation replicates GOA Thornyhead model converge. Figure 2 top row within panel shows distribution resulting recovered parameter estimates models based simulations, horizontal line median estimated value black dot true value (.e., parameters estimates real data sets).  demonstrated boxplots bottom rows Figure 2, parameters models recovered well simulation testing, low median relative error. models instances long negative tails log standard deviation process error (log_PE; Figure 2 top row within panel), suggesting PE small tended towards zero simulations (recall parameters estimated log space, negative number log space close zero natural space). also case extra observation error biomass survey (log_tau_biomass) GOA Thornyhead model (Figure 2, top row within lower panel).  outlier replicates, occur optimizer used rema (nlminb()) returns “NA/NaN function evaluation” error, help shed light replicates failed converge. hood, optimizer pushing PE towards zero (PE=0 taking global mean biomass time series), violating central assumption REMA model biomass underlying trend (PE > 0, thus PE can log-transformed). negative log-scale values extreme GOA Thornyhead model, raise concern potential model misspecification -parameterization. example, given trade-process observation error REMA (including additional observation dampens biomass trend estimated, pushing PE towards zero), possible estimation additional observation error biomass survey (log_tau_biomass) might justified. Figure 2. Simulation self-test results AI Pcod (top panels; gold) GOA Thornyhead (bottom panels; teal). upper row panel gives distribution recovered parameters, dot giving “true value” used simulate data. lower row panel gives relative error recovered parameters, zero line indicating simulation returns “true value,” box plot line giving median recovered parameter’s relative errors, whiskers dots boxplots giving spread recovered parameter’s relative errors.","code":"sim_test <- function(mod_name, replicates, cpue) { # storage things re_est <- matrix(NA, replicates, length(mod_name$par)) # parameter estimates # go through the model suppressMessages(for(i in 1:replicates) { sim <- mod_name$simulate(complete = TRUE) # simulates the data # simulated biomass observations: tmp_biomass <- matrix(data = exp(sim$log_biomass_obs), ncol = ncol(mod_name$input$data$biomass_obs)) colnames(tmp_biomass) <- colnames(mod_name$data$biomass_obs) # simulated cpue observations, when applicable: if (cpue) {tmp_cpue <- matrix(data = sim$cpue_obs, ncol = ncol(mod_name$input$data$cpue_obs)) colnames(tmp_cpue) <- colnames(mod_name$data$cpue_obs)} # set up new data for input newinput <- mod_name$input newinput$data$biomass_obs <- tmp_biomass # biomass data if (cpue) {newinput$data$cpue_obs <- tmp_cpue} # cpue data # create \"obsvec\" which is used internally in cpp file as the observation # vector for all observation (log biomass + cpue) in the likelihood # functions (and required for OSA residuals). note the transpose t() needed # to get these matrices in the correct order (by row instead of by col) -- # note obsvec is masked from users normally in prepare_rema_input() newinput$data$obsvec <- t(sim$log_biomass_obs)[!is.na(t(sim$log_biomass_obs))] if (cpue) {newinput$data$obsvec <- c(t(sim$log_biomass_obs)[!is.na(t(sim$log_biomass_obs))], t(sim$log_cpue_obs)[!is.na(t(sim$log_cpue_obs))])} # refit model mod_new <- fit_rema(newinput, do.sdrep = FALSE) # add parameter estimates to matrix if(mod_new$opt$convergence == 0) { re_est[i, ] <- mod_new$env$last.par[1:length(mod_name$par)] } else { re_est[i, ] <- rep(NA, length(mod_name$par)) } }) re_est <- as.data.frame(re_est); re_est$type <- rep(\"recovered\") return(re_est) } # run for pcod and prep data frame # run simulation testing par_ests <- sim_test(mod_name = pcod_mod, replicates = 500, cpue = FALSE) # get data frame with simulation values for each parameter n_not_converged_pcod <- length(which(is.na(par_ests[,1]))); n_pcod <- length(par_ests[,1]) prop_converged_pcod <- 1-n_not_converged_pcod/n_pcod mod_par_ests <- data.frame(\"log_PE1\" = pcod_mod$env$last.par[1], type = \"model\") names(par_ests) <- names(mod_par_ests) # rename for ease pcod_par_ests <- rbind(mod_par_ests, par_ests) # recovered and model in one data frame pcod_par_ests$sp <- rep(\"AI Pcod\") # run for thorny and prep data frame -- same process # run simulation testing par_ests <- sim_test(mod_name = thrn_mod, replicates = 500, cpue = TRUE) # note some warnings # sometimes spits out: In stats::nlminb(model$par, model$fn, model$gr, control = # list(iter.max = 1000,...: NA/NaN function evaluation n_not_converged_thrn <- length(which(is.na(par_ests[,1]))); n_thrn <- length(par_ests[,1]) prop_converged_thrn <- 1-n_not_converged_thrn/n_thrn nrow(par_ests) # get data frame with simulation values for each parameter mod_par_ests <- data.frame(\"log_PE1\" = thrn_mod$env$last.par[1], \"log_PE2\" = thrn_mod$env$last.par[2], \"log_PE3\" = thrn_mod$env$last.par[3], \"log_q\" = thrn_mod$env$last.par[4], \"log_tau_biomass\" = thrn_mod$env$last.par[5], \"log_tau_cpue\" = thrn_mod$env$last.par[6], type = \"model\") names(par_ests) <- names(mod_par_ests) # rename for ease thrn_par_ests <- rbind(mod_par_ests, par_ests) # recovered and model parameters in one data frame thrn_par_ests$sp <- rep(\"GOA Thornyhead\") # reogranize data pcod_sim <- pcod_par_ests %>% pivot_longer(1, names_to = \"parameter\") thrn_sim <- thrn_par_ests %>% pivot_longer(1:6, names_to = \"parameter\") sim_dat <- rbind(pcod_sim, thrn_sim) # Relative Error: ((om-em)/om) sim_re <- sim_dat %>% pivot_wider(id_cols = c(\"sp\", \"parameter\"), names_from = type, values_from = value) %>% unnest(cols = c(model, recovered)) %>% mutate(RE = (model-recovered)/model*100) %>% group_by(sp, parameter) %>% mutate(label = paste0(\"Median RE=\", formatC(median(RE, na.rm = TRUE), format = \"f\", digits = 1), \"%\")) %>% suppressWarnings() # plot distribution of simulated parameter estimates plot_sim <- function(sim_dat, plot_title, fill_col = \"#21918c\") { ggplot(NULL, aes(parameter, value)) + # add distribution of recovered parameters geom_violin(data = sim_dat %>% filter(type == \"recovered\"), fill = fill_col, alpha = 0.6, draw_quantiles = 0.5) + # add \"true values\" from original model geom_point(data = sim_dat %>% filter(type == \"model\"), size = 2, col = \"black\") + facet_wrap(~ parameter, scales = \"free\", nrow = 1) + labs(x = NULL, y = \"Parameter estimate\", title = plot_title, subtitle = \"Distribution of parameters estimates (median=horizontal line, true value=point)\") + scale_x_discrete(labels = NULL, breaks = NULL) } # plot relative error plot_re <- function(sim_re, fill_col = \"#21918c\") { ggplot(NULL, aes(parameter, RE)) + # add distribution of recovered parameters geom_boxplot(data = sim_re, alpha = 0.6, fill = fill_col, na.rm = TRUE, outlier.size = 0.8) + geom_hline(yintercept = 0) + # separate by parameter facet_wrap(~ parameter+label, scales = \"free\", nrow = 1) + #, space = \"free\") + labs(x = NULL, y = \"Relative error (%)\", subtitle = \"Distribution of relative error (RE; i.e., (true-estimated values)/true value*100)\") + scale_x_discrete(labels = NULL, breaks = NULL) } p1pcod <- plot_sim(pcod_sim, \"AI Pcod Simulation\", fill_col = \"goldenrod\") p1thrn <- plot_sim(thrn_sim, \"GOA Thornyhead Simulation\") p2pcod <- plot_re(sim_re %>% filter(sp == \"AI Pcod\"), fill_col = \"goldenrod\") p2thrn <- plot_re(sim_re %>% filter(sp == \"GOA Thornyhead\")) cowplot::plot_grid(p1pcod, p2pcod, ncol = 1) ggsave(paste0(\"vignettes/ex4_sim_pcod.png\"), width = 6.5, height = 7, units = \"in\", bg = \"white\") cowplot::plot_grid(p1thrn, p2thrn, ncol = 1) ggsave(paste0(\"vignettes/ex4_sim_thrn.png\"), width = 11, height = 7, units = \"in\", bg = \"white\")"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"residual-analysis","dir":"Articles","previous_headings":"","what":"2. Residual analysis","title":"REMA model validation","text":"Residuals used test underlying assumptions structure error distributions model. example, linear regression, residuals independent, normal, constant variance. Traditional residuals (e.g., Pearson’s residuals) inappropriate state-space models like REMA, random effects induce correlations predicted data residuals longer independent. Additionally, process error variance may overestimated cases model mis-specified, thus leading artificially small residuals (see Section 3 Thygesen et al. 2017 example). appropriate residual type validation state-space models one-step ahead (OSA) residuals. Instead comparing observed expected values time step, OSA residuals use forecasted values based previous observations (.e., exclude observed value current time step prediction). way, OSA residuals account non-normality correlation residuals among years.  correctly specified model, resultant OSA residuals independent identically distributed (..d.) standard normal distribution N(0,1)N(0,1). can tested normal QQ plot, theoretical quantile values standard normal distribution plotted x-axis, corresponding empirical quantile values OSA residuals plotted y-axis. OSA residuals standard normally distributed, points fall 0/1 reference line, standard deviation normalized residuals (SDNR) statistic close 1. residuals standard normally distributed (SDNR far 1), points deviate reference line. case model includes multiple strata surveys, OSA residuals normally distributed within across strata surveys, however small sample sizes (e.g., within single stratum) may sufficient statistical power identify model misfit.  addition normality OSA residuals, residuals random independent (.e., correlated year). Standard approaches autocorrelation function (ACF) plots residual runs tests (e.g., Wald-Wolfowitz test) inappropriate many REMA applications missing years data survey time series, instead, directly plot residuals years survey time series. Visual patterns extreme outliers (greater 3 less -3 N(0,1) distribution) residuals can indicate structure data (.e., temporal correlation) adequately captured model.  Methods calculating OSA residuals REMA implemented rema::get_osa_residuals() using TMB::oneStepPredict() function fullGaussian method using lognormal error distributions. rema::get_osa_residuals() returns tidied dataframes observations, REMA model predictions, OSA residuals biomass CPUE survey data appropriate, along QQ residual-time series diagnostic plots.  Results residual analysis normal QQ plot AI Pcod (Figure 3, top panel) suggests slight negative skewness residuals (.e., majority points fall 0/1 line), though small sample size makes interpretation challenging. SDNR 0.99 (recall perfect model SDNR=1.00), indicating  assumptions likely met residuals. plot OSA residuals year (Figure 3, bottom panel) shows obvious patterns residuals, suggesting independent; evidence outliers residuals warrant inspection. combined normal QQ plot GOA Thornyhead OSA residuals (Figure 4, top panel) suggests follow normal distribution (points fall along 0/1 line), though evidence slight positive, right skewness (majority points fall 0/1 line). SDNR approximately 1, indicating variance assumptions likely met residuals. QQ plots biomass strata (Figure 4, middle panels) highlight positive skewness may coming (e.g., EGOA 701-1000m, WGOA 0-500 m); however, small sample sizes stratum make interpretation QQ plots difficult. QQ plots CPUE strata (Figure 4, bottom panel) mostly normal, though evidence light tails WGOA stratum. plot OSA residuals year GOA Thornyhead OSA residuals biomass strata (Figure 5, top panels) show patterns residuals, suggesting independent. However, runs residuals CPUE strata (Figure 5, bottom panels), especially CGOA WGOA, might indicate model misspecification misfit. evidence outliers. Figure 3. One-step ahead (OSA) residual plots AI Pcod. top panel gives QQ plot, comparing normal distribution quantiles (x-axis) OSA residual distribution quantiles (y-axis). diagonal line 0-1 line, two quantiles equal, SDNR statistic given upper left plot. bottom panel gives residuals (y-axis) plotted year (x-axis) used check independence. Figure 4. One-step ahead (OSA) QQ plots across (top panel) within strata (middle panel, bottom trawl survey strata; bottom panel, longline survey strata) GOA Thornyhead. panels, colors correspond bottom trawl survey strata. SDNR statistic given across strata QQ plot; within strata plots useful check extreme skew lack statistical power. Figure 5. One-step ahead (OSA) residual plots GOA Thornyhead. residuals (y-axis) annual time step (random effect; x-axis) survey stratum.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"laplace-approximation-are-the-model-assumptions-related-to-random-effects-estimation-reasonable","dir":"Articles","previous_headings":"","what":"3. Laplace approximation: Are the model assumptions related to random effects estimation reasonable?","title":"REMA model validation","text":"rema library developed Template Model Builder (TMB), uses maximum marginal likelihood estimation Laplace approximation efficiently estimate high dimensional, non-linear mixed effects models frequentist framework (Skaug Fournier, 2006; Kristensen et al., 2016). primary assumption models using Laplace approximation random effects follow normal distribution. assumption simplifies complex integrals make likelihood function. Laplace approximation fast accurate normality assumption met; however, true distribution random effects normal, Laplace approximation may introduce bias parameter estimates. can test validity Laplace approximation using Markov chain Monte Carlo (MCMC) sampling tmbstan library, comparing distributions fixed effects parameters MCMC-sampled models without Laplace approximation (Monnahan Kristensen, 2018). , compare distributions using QQ plots two model cases. Laplace approximation reasonably accurate sampling quantiles two models (without Laplace approximation) similar, .e., fall 1:1 line. can also help us identify bias introduced Laplace approximation, example, fixed effects estimates associated estimates uncertainty differ base model (run marginal maximum likelihood estimation using Laplace approximation) MCMC versions without Laplace approximation. MCMC models run 1000 iterations, assuming warmup 500. discussed Results section, iterations increased 5000 warmup 2500 GOA Thornyhead effort improve diagnostics model convergence. function used run two model cases: following code uses output mcmc_comp() make QQ plot compare posterior distributions parameters parameter estimates MCMC models without Laplace approximation: Results testing validity Laplace approximation Figure 6 gives compared quantiles full MCMC testing (x-axis) Laplace approximation (y-axis) AI Pcod GOA Thornyhead models. simpler AI Pcod model (Figure 6, left panel), Laplace approximation full MCMC testing show similar distribution, indicating Laplace approximation reasonable. Additionally estimates fixed effects nearly identical base model MCMC models without Laplace approximation (Table 1). complex GOA Thornyhead model, preliminary model results based 1000 iterations exhibited poor diagnostics process error parameters additional observation error biomass survey (log_tau_biomass). increased number MCMC iterations 5000 effort improve convergence diagnostics; however, log_tau_biomass parameter continued show significant deviations models without Laplace approximation (Figure 6, right panels). comparison fixed effects estimates base models show large differences parameter estimates associated estimates uncertainty (Table 1). indicates Laplace approximation likely inaccurate, investigation model structure might necessary. Table 1. AI Pcod GOA Thornyhead fixed effect estimates 95% confidence credible intervals maximum likelihood MCMC models, respectively (CI). Fixed effects estimates compared base model run marginal maximum likelihood estimation assuming Laplace approximation (“Base_MMLE_Laplace”) MCMC models without Laplace approximation (“MCMC_Laplace” “MCMC_withoutLaplace”, respectively). Rhat statistic provided MCMC models convergence metric (convergence = Rhat = 1). Nonconverged parameters highlighted bold. parameter estimate MCMC models median posterior samples. Figure 6. Results Laplace approximation testing. compare results running model without assuming normality random effects (MCMC, x-axis) running model assuming normality random effects (Laplace approximation, y-axis), quantiles plotted fixed effect model. gray line 1:1 line; points fall gray line indicate quantile value two cases. top panel gives plot AI Pcod model (one fixed effect), bottom panels give plots GOA Thornyhead model (six fixed effects).","code":"# function to (1) run models with and without laplace approximation for # comparison and (2) return posterior data frames and models # function input: model (e.g., AI Pcod or GOA Thornyhead), number of iterations # (samples), number of chains mcmc_comp <- function(mod_name, it_numb, chain_numb) { # set up MCMC chain information it_num <- it_numb chain_num <- chain_numb # mod_name = pcod_mod; it_num = 1000; chain_num = 4 # run model with laplace approximation mod_la <- tmbstan(obj = mod_name, chains = chain_num, init = mod_name$par, laplace = TRUE, iter = it_num) # run model without laplace approximation, i.e., all parameters fully estimated without assumptions of normality mod_mcmc <- tmbstan(obj = mod_name, chains = chain_num, init = mod_name$par, laplace = FALSE, iter = it_num) # posteriors as data frame post_la <- as.data.frame(mod_la); post_la$type <- (\"la\") post_mcmc <- as.data.frame(mod_mcmc); post_mcmc <- post_mcmc[, c(1:length(mod_name$par), dim(post_mcmc)[2])]; post_mcmc$type <- rep(\"mcmc\") # informational things... this is for getting the posterior draws, i.e., to # test traceplots and such post_la$chain <- rep(1:chain_num, each = it_num/2) post_la$iter_num <- rep(1:(it_num/2), chain_num) post_mcmc$chain <- rep(1:chain_num, each = it_num/2) post_mcmc$iter_num <- rep(1:(it_num/2), chain_num) post_draws <- rbind(post_la, post_mcmc) # get quantiles qv <- seq(from = 0, to = 1, by = 0.01) quant_dat <- data.frame(quant = NULL, la = NULL, mcmc = NULL, par = NULL) for (i in 1:(dim(post_la)[2]-4)) { # post_la has type, chains, lp, iteration number column that don\"t count tmp <- data.frame(quant = qv, la = quantile(unlist(post_la[i]), probs = qv), mcmc = quantile(unlist(post_mcmc[i]), probs = qv), par = rep(paste0(\"V\", i))) quant_dat <- rbind(quant_dat, tmp) } return(list(quant_dat, post_draws, mod_la, mod_mcmc)) } # run models pcod_comp <- mcmc_comp(mod_name = pcod_mod, it_numb = 1000, chain_numb = 3) thrn_comp_log_tau <- mcmc_comp(mod_name = thrn_mod, it_numb = 5000, chain_numb = 3) # clean up data frame names by renaming things pcod_qq <- pcod_comp[[1]] pcod_qq$par_name <- rep(\"log_PE1\"); pcod_qq$sp <- rep(\"AI Pcod\") thrn_qq <- thrn_comp_log_tau[[1]] thrn_qq$par_name <- recode(thrn_qq$par, V1 = \"log_PE1\", V2 = \"log_PE2\", V3 = \"log_PE3\", V4 = \"log_q\", V5 = \"log_tau_biomass\", V6 = \"log_tau_cpue\") thrn_qq$sp <- rep(\"GOA Thornyhead\") qq_dat <- rbind(pcod_qq, thrn_qq) # plot plot_laplace_mcmc <- function(qq_dat, plot_title) { ggplot(qq_dat, aes(mcmc, la)) + geom_abline(intercept = 0, slope = 1, col = \"lightgray\") + geom_point() + facet_wrap(~par_name, scales = \"free\", nrow = 2, dir = \"v\") + labs(x = \"MCMC\", y = \"Laplace approx.\", title = plot_title) + theme(plot.title = element_text(size = 12), axis.text = element_text(size = 7), axis.title = element_text(size = 11), strip.text = element_text(size = 11)) } p1 <- plot_laplace_mcmc(pcod_qq, \"AI Pcod\") p2 <- plot_laplace_mcmc(thrn_qq, \"GOA Thornyhead\") cowplot::plot_grid(p1, p2, rel_widths = c(0.4,0.6), ncol = 2) ggsave(paste0(\"vignettes/ex4_mcmc_qq.png\"), width = 11, height = 7, units = \"in\", bg = \"white\") # parameter summary table - compare the marginal max likelihood estimates (MMLE) # with the mcmc estimates. Rhat is the potential scale reduction factor on split # chains (at convergence, Rhat=1) m <- summary(pcod_comp[[3]])$summary[1, c(6,4,8,10), drop = FALSE] psum <- data.frame(Stock = 'AI Pcod', Model = 'MCMC_Laplace', FixedEffect = row.names(m), m, row.names = NULL) m <- summary(pcod_comp[[4]])$summary[1, c(6,4,8,10), drop = FALSE] psum <- psum %>% bind_rows(data.frame(Stock = 'AI Pcod', Model = 'MCMC_withoutLaplace', FixedEffect = row.names(m), m, row.names = NULL)) m <- summary(thrn_comp_log_tau[[3]])$summary[1:6, c(6,4,8,10), drop = FALSE] psum <- psum %>% bind_rows(data.frame(Stock = 'GOA Thornyhead', Model = 'MCMC_Laplace', FixedEffect = row.names(m), m, row.names = NULL)) m <- summary(thrn_comp_log_tau[[4]])$summary[1:6, c(6,4,8,10), drop = FALSE] psum <- psum %>% bind_rows(data.frame(Stock = 'GOA Thornyhead', Model = 'MCMC_withoutLaplace', FixedEffect = row.names(m), m, row.names = NULL)) pars <- as.list(pcod_mod$sdrep, \"Est\")$log_PE sd <- as.list(pcod_mod$sdrep, \"Std\")$log_PE psum <- psum %>% bind_rows(data.frame(Stock = 'AI Pcod', Model = 'Base_MMLE_Laplace', FixedEffect = \"log_PE\", X50. = pars, X2.5. = pars-1.96*sd, X97.5. = pars+1.96*sd, Rhat = NA)) pars <- unlist(as.list(thrn_mod$sdrep, \"Est\")) pars <- pars[names(pars)[grepl(\"log_PE|log_q|log_tau\", names(pars))]] sd <- unlist(as.list(thrn_mod$sdrep, \"Std\")) sd <- sd[names(sd)[grepl(\"log_PE|log_q|log_tau\", names(sd))]] psum <- psum %>% bind_rows(data.frame(Stock = 'GOA Thornyhead', Model = 'Base_MMLE_Laplace', FixedEffect = rownames(summary(thrn_comp_log_tau[[3]])$summary[1:6,,drop=FALSE]), X50. = pars, sd = sd, row.names = NULL) %>% mutate(X2.5. = X50.-1.96*sd, X97.5. = X50.+1.96*sd, Rhat = NA) %>% select(-sd)) # write.csv(psum, 'vignettes/ex4_param_table.csv')"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"parameter-correlation-are-the-model-parameters-identifiable-and-non-redundant","dir":"Articles","previous_headings":"","what":"4. Parameter correlation: Are the model parameters identifiable and non-redundant?","title":"REMA model validation","text":"Parameter redundancy refers idea multiple parameters contribute model way. intuitive case y∼β1+β2+αxy \\sim \\beta_1 + \\beta_2 + \\alpha x, since model estimate many combinations β1\\beta_1 β2\\beta_2 minimize log-likelihood, .e., sampled parameters correlated. reduce redundancy, β1+β2\\beta_1 + \\beta_2 can redefined β0\\beta_0.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"results-of-parameter-correlation-analysis","dir":"Articles","previous_headings":"4. Parameter correlation: Are the model parameters identifiable and non-redundant?","what":"Results of parameter correlation analysis","title":"REMA model validation","text":"Figure 7 shows pairwise correlation matrix posterior draws GOA Thornyhead model (MCMC model Laplace approximation), histograms univariate marginal distributions diagonal scatterplot bivariate distributions diagonal. significant unresolved sampling problems log_tau_biomass parameter, shown non-normal heavily skewed distribution log_tau_biomass. pairwise scatterplots -diagonals show heavy left tails log_tau_biomass causing interactions (.e., correlation) almost parameters model. diagnostic clear indication model mis-specification warrants investigation. Figure 7. Pairs plot GOA Thornyhead MCMC model assuming Laplace approximation. diagonal plots give distribution fixed effect 2500 MCMC draws. -diagonal elements give parameter--parameter correlation, point gives parameter values estimated MCMC simulation.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"is-the-goa-thornyhead-model-really-converged","dir":"Articles","previous_headings":"","what":"5. Is the GOA Thornyhead model really converged?","title":"REMA model validation","text":"Despite poor diagnostics shown vignette GOA Thornyhead model, standard output TMB rema suggests model converged (.e., low maximum gradient component standard error estimates appear reasonable). Within MCMC framework, can look mixing MCMC across multiple chains using diagnostic called traceplot assess convergence. following code shows within rstan library (Stan Development Team, 2024): traceplots models figures , see AI Pcod model MCMC chains fully mixed model converged (top panel). However, GOA Thornyhead model (bottom panel), see estimation additional observation error (log_tau_biomass) causing convergence issues MCMC chains well mixed. Figure 8.  MCMC traceplots AI Pcod (top panel) GOA Thornyhead (bottom panels). subpanel individual parameter, plot gives value parameter (y-axis) simulation draw (x-axis) model chain (color).","code":"p1_mcmc_pcod <- rstan::traceplot(pcod_comp[[3]]) + ggtitle(label = \"Traceplots to assess mixing across Markov chains and convergence\", subtitle = \"AI Pcod\") p1_mcmc_thrn <- rstan::traceplot(thrn_comp_log_tau[[3]], ncol = 3) + ggtitle(label = NULL, subtitle = \"GOA Thornyhead\") cowplot::plot_grid(p1_mcmc_pcod, p1_mcmc_thrn, ncol = 1) ggsave(paste0(\"vignettes/ex4_traceplots.png\"), width = 11, height = 8, units = \"in\", bg = \"white\")"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"when-should-we-care-about-failed-model-diagnostics","dir":"Articles","previous_headings":"","what":"When should we care about failed model diagnostics?","title":"REMA model validation","text":"REMA model used within NPFMC obtain exploitable biomass estimates Tier 4 crab Tier 5 groundfish method apportion Acceptable Biological Catches among management regions (Sullivan et al., 2022). , primary purpose smooth noisy survey biomass estimates (.e., using fancy average), need care potential red flags raised model validation?  Yes: case GOA Thornyhead model, example, appears issue estimation additional observation error, pointing potential -parameterization model /parameter redundancy. terminal year predictions REMA model quantities directly used management, estimation additional observation error definition increases “smoothness” model predictions, answer likely yes particular case.  Maybe : hypothetical case estimation year predictions show diagnostic concerns, uncertainty estimations show possible diagnostic problems, model validation red flags might less important. uncertainty results REMA generally used fisheries management NPFMC. words, failed model validation criteria primarily concern impact quantities interest model (case, biomass predictions). Model validation, subsequent interpretation, considered light goal model, rather binary pass/fail testing procedure.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"preliminary-recommendations","dir":"Articles","previous_headings":"","what":"Preliminary recommendations","title":"REMA model validation","text":"model validation methods presented tools stock assessment scientists evaluate existing models guide development future models. analysis, found REMA model can recover parameters across range complexity relevant NPFMC Tier 4 Tier 5 stocks. analysis suggests MCMC diagnostics OSA residuals useful diagnostics REMA model. models fail diagnostics, may indicate model complex simplifications considered. highlighted several places document diagnostic features can help make choices appropriate model simplifications (pool process errors, remove additional observation errors, etc.) light REMA model assumptions trade-offs model parameters. results indicate diagnostics important explore existing models using multivariate configurations multiple process errors, multi-survey versions using CPUE indices like longline survey, models estimating additional observation error. model validation may needed existing REMA models, recommend used recommending new models management.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"references","dir":"Articles","previous_headings":"","what":"References","title":"REMA model validation","text":"Aeberhard, W.H., Flemming, J.M., Nielsen, . 2018. Review State-Space Models Fisheries Science. Annual Review Statistics Application, 5(1), 215-235. https://doi.org/10.1146/annurev-statistics-031017-100427 Auger-Méthé, M., Field, C., Albertsen, C.M., Derocher, .E., Lewis, M.., Jonsen, .D., Flemming, J.M. 2016. State-space models’ dirty little secrets: even simple linear Gaussian models can estimation problems. Scientific reports, 6(1), 26677. Auger‐Méthé, M., Newman, K., Cole, D., Empacher, F., Gryba, R., King, . Leos-Barajas, V., Flemming, J.M., Nielsen, , Petris, G., Thomas, L. 2021. guide state–space modeling ecological time series. Ecological Monographs, 91(4), e01470. Campbell, D., & Lele, S. 2014. ANOVA test parameter estimability using data cloning application statistical inference dynamic systems. Computational Statistics & Data Analysis, 70, 257-267. Cole, D. J. 2019. Parameter redundancy identifiability hidden Markov models. Metron, 77(2), 105-118. Echave, K.B., Siwicke, K.., Sullivan, J., Ferriss, B., Hulson, P.J.F. 2022. Assessment Thornyhead stock complex Gulf Alaska. Stock assessment fishery evaluation report groundfish fisheries Gulf Alaska. North Pacific Fishery Management Council, 605 W. 4th. Avenue, Suite 306, Anchorage, AK 99501. https://apps-afsc.fisheries.noaa.gov/Plan_Team/2022/GOAthorny.pdf Gabry, J. Mahr, T. 2024. bayesplot: Plotting Bayesian Models. R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry, J., Simpson, D., Vehtari, ., Betancourt, M., Gelman, . 2019. Visualization Bayesian workflow. J. R. Stat. Soc. . (182) 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378. Kristensen, K., Nielsen, ., Berg, C.W., Skaug, H., Bell, B.M. 2016. TMB: Automatic differentiation Laplace approximation. J Stat Softw. 2016; 70(5):21. Epub 2016-04-04. https://doi.org/10.18637/jss.v070.i05 Monnahan, C.C., & Kristensen, K. 2018. -U-turn sampling fast Bayesian inference ADMB TMB: Introducing adnuts tmbstan R packages. PloS one, 13(5), e0197954. Skaug, H.J., Fournier, D.. 2006. Automatic approximation marginal likelihood non-Gaussian hierarchical models. Computational Statistics & Data Analysis, 51(2), 699-709. Spies, ., Barbeaux, S., Hulson, H., Ortiz, . 2023. Assessment Pacifc cod stock Aleutian Islands. Stock assessment fishery evaluation report groundfish fisheries Bering Sea Aleutian Islands. North Pacific Fishery Management Council, 605 W. 4th. Avenue, Suite 306, Anchorage, AK 99501. https://apps-afsc.fisheries.noaa.gov/Plan_Team/2023/AIpcod.pdf Stan Development Team. 2024. RStan: R interface Stan. R package version 2.32.6, https://mc-stan.org/. Sullivan, J., Monnahan, C. Hulson, P. Ianelli, J. Thorson, J. Havron, . 2022. REMA: consensus version random effects model ABC apportionment Tier 4/5 assessments. Plan Team Report, Joint Groundfish Plan Teams, North Pacific Fishery Management Council. 605 W 4th Ave, Suite 306 Anchorage, AK 99501. Available Oct 2022 Joint GPT e-Agenda. Thygesen, U.H., Albertsen, C.M., Berg, C.W., Kristensen, K., & Nielsen, . 2017. Validation ecological state space models using Laplace approximation. Environmental Ecological Statistics, 24(2), 317-339. https://doi.org/10.1007/s10651-017-0372-4","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"base-model-structure-for-a-single-survey-and-stratum","dir":"Articles","previous_headings":"","what":"Base model structure for a single survey and stratum","title":"REMA model equations","text":"base REMA model can represented simple state-space random walk model added noise. observation model comprised index log-transformed annual survey biomass estimates ln(Bt)ln(B_t) associated standard deviations σln(Bt)\\sigma_{ln(B_t)}, σln(Bt)\\sigma_{ln(B_t)} approximated using coefficient variation BtB_t (σBt/Bt\\sigma_{B_t}/B_t) σln(Bt)=ln((σBtBt)2+1). \\sigma_{ln(B_t)}=\\sqrt{ln((\\frac{\\sigma_{B_t}}{B_t})^2+1)}. measurement observation equation, describes relationship observed survey biomass ln(Bt)ln(B_t) latent state variable, population biomass ln(B̂t)ln(\\hat{B}_t), expressed ln(Bt)=ln(B̂t)+ϵB,ϵB∼N(0,σln(Bt)2). ln(B_t) = ln(\\hat{B}_t) + \\epsilon_{B}, \\text{} \\epsilon_{B} \\sim N(0, \\sigma_{ln(B_t)}^2). state equation associated process error variance σPE2\\sigma_{PE}^2 defined ln(B̂t)=ln(B̂t−1)+ηt−1,ηt∼N(0,σPE2), ln(\\hat{B}_{t}) = ln(\\hat{B}_{t-1})+\\eta_{t-1}, \\text{} \\eta_t \\sim N(0, \\sigma_{PE}^2), initial lnB̂1ln\\hat{B}_{1} constrained random walk process. base model, process error variance σPE2\\sigma_{PE}^2 fixed effect parameter estimated, unobserved population biomass ln(B̂t)ln(\\hat{B}_t) estimated series random effects. model fit using marginal maximum likelihood estimation, marginal negative log-likelihood obtained via Laplace approximation (Skaug Fournier 2006). reproducible example fitting univariate (.e., single stratum, single survey) version TMB model comparing results ADMB available REMA basics vignette.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"extending-to-multiple-biomass-survey-strata","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Extending to multiple biomass survey strata","title":"REMA model equations","text":"single survey, single stratum version model can extended include one additional strata jj survey. inclusion multiple strata model advantageous scenarios apportionment biomass among areas desired result. extension assumes correlation observation errors among survey strata, though PE variance can shared estimated independently among strata. ADMB TMB versions REMA model utilize different methods estimating variance total predicted biomass. Therefore, strata-specific estimates predicted biomass associated confidence intervals close identical ADMB TMB versions model, confidence intervals total predicted biomass differ slightly. original ADMB code, Marlow (1967) method applied, total variance σJ2\\sigma_J^2 approximated σJ2=ln(∑e2B̂j+σB̂j2(eσB̂j2−1)(∑eB̂j+σB̂j2/2)2+1). \\sigma_J^2 = ln\\Bigg(\\frac{\\sum{e^{2 \\hat{B}_j + \\sigma_{\\hat{B}_j}^2} (e^{\\sigma_{\\hat{B}_j}^2 - 1})}} {(\\sum{e^{\\hat{B}_j + \\sigma_{\\hat{B}_j}^2/2}})^2} + 1\\Bigg). updated rema package, total variance estimated using standard Delta method can replicated ADMB using sdreport_number TMB using ADREPORT macro. described Monnahan et al. (2021), exploration methods summing log-normal variables research topic potential impacts beyond scope package. reproducible example fitting multivariate version REMA model TMB comparison results ADMB available REMA basics vignette. example, also show Akaike Information Criteria (AIC) can used explore inclusion strata-specific versus single, shared process error variance.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"addition-of-an-auxiliary-catch-per-unit-effort-cpue-survey","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Addition of an auxiliary catch per unit effort (CPUE) survey","title":"REMA model equations","text":"situations auxiliary biomass catch per unit effort (CPUE) survey ItI_{t} associated variance σIt\\sigma_{I_{t}} available, additional scaling parameter qq (estimated log-space) can estimated facilitate inclusion new information biomass predictions. predicted annual CPUE survey index ̂t\\hat{}_{t} calculated ̂t=eq*eB̂t \\hat{}_{t} = e^{q} * e^{\\hat{B}_{t}} CPUE survey observations additional measurement equation likelihood component similar biomass survey: ln()=ln(̂t)+ϵI,ϵI∼N(0,σln()2), ln(I_{t}) = ln(\\hat{}_{t}) + \\epsilon_{}, \\text{} \\epsilon_{} \\sim N(0, \\sigma_{ln(I_t)}^2), default, rema estimates single qq stratum; however, users can alternatively specify shared qq parameters among strata. strata definitions t,biomass CPUE surveys (e.g. biomass estimated finer geographic resolution CPUE), user can define relationship two surveys’ strata using q_options argument prepare_rema_input() function (see q_options pointer_biomass_cpue_strata). However, since auxiliary CPUE survey index related unobserved population biomass level biomass survey strata, REMA model can accommodate scenarios CPUE survey strata coarser resolution equivalent biomass survey strata. reproducible example fitting two-survey version REMA model provided Fitting additional CPUE survey vignette. example also describes error previously used ADMB version REMA model, presents several alternative models estimate additional observation error. topics described detail following two sections.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"the-admb-version-of-rema","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum > Addition of an auxiliary catch per unit effort (CPUE) survey","what":"The ADMB version of REMA","title":"REMA model equations","text":"Separability, implemented SEPARABLE_FUNCTION ADMB automatic detection TMB, increases computational efficiency Laplace approximation breaking complex, multivariate integrals product simpler, univariate integrals (Fournier et al. 2012). use separability particularly relevant state-space models relates efficiency Laplace approximation marginal negative log-likelihood integration random effects joint negative log-likelihood. random effects implementation ADMB, parameters defined inside PARAMETER_SECTION template file used within SEPARABLE_FUNCTION unless passed arguments function (Skaug Fournier 2013). ADMB version REMA model, calculation predicted CPUE [.e. ln(̂t)=ln(q)+ln(B̂)ln(\\hat{}_{t})=ln(q)+ln(\\hat{B})] occurred outside rather inside SEPARABLE_FUNCTION, results, violated rule affecting accuracy Laplace approximation. explored error simplified example REMA model ADMB TMB. able reproduce results TMB version REMA ADMB passing ln(q)ln(q) argument SEPARABLE_FUNCTION performing ln(̂t)ln(\\hat{}_{t}) calculation inside SEPARABLE_FUNCTION. used Markov Chain Monte Carlo (MCMC) methods statistical inference instead MLE, results ADMB versions closely matched correct version REMA. Finally, comparison individual negative log-likelihood (NLL) components, along joint marginal NLLs, revealed TMB model, correct “inside” ADMB version, incorrect “outside” ADMB version, except marginal NLL. Taken together, analysis confirms bug “outside” ADMB version, affects accuracy Laplace approximation model results.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"estimation-of-additional-observation-error","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Estimation of additional observation error","title":"REMA model equations","text":"Based experience gained using alternative observed index estimates (e.g., relative CPUE indices), appears cases estimates observation error variances biomass /CPUE survey low (e.g., Echave et al. 2020). , mismatch biologically reasonable inter-annual variability precision index estimates. instances, model estimates (σBt,j2+σIt,j2)/σPE2(\\sigma^2_{B_{t,j}}+\\sigma^2_{I_{t,j}})/\\sigma^2_{PE} may lower expected based individual species life history traits. example, ratio observation PE variation low, model predictions population biomass may exhibit high inter-annual variability. behavior unexpected low productivity species, exhibit low inter-annual variation biomass (.e. low PE variance), especially situations fishing exploitation low. One approach addressing issue estimate additional observation error. method commonly implemented Stock Synthesis, implemented Alaskan crab stock assessments (Zheng Siddeek 2020), explored several groundfish assessments well. Using biomass survey variance example, extra estimated error (στ\\sigma_{\\tau}) specified : σln(Bt)=ln((σBtBt)2+στ2+1). \\sigma_{ln(B_t)}=\\sqrt{ln((\\frac{\\sigma_{B_t}}{B_t})^2+{\\sigma_{\\tau}}^2+1)}. approach new method Tier 5 stock assessments apportionment methods AFSC. reproducible example GOA shortraker outlined Fitting additional CPUE survey. example, estimation additional observation error biomass CPUE survey resulted better fit AIC status quo approaches biologically realistic trend predicted population biomass.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"exploration-of-the-tweedie-distribution-for-zero-biomass-observations","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Exploration of the Tweedie distribution for zero biomass observations","title":"REMA model equations","text":"Tweedie distributions family probability distributions can generalized include Gaussian, inverse Gaussian, gamma, Poisson, compound Poisson-gamma distributions. Tweedie distribution can defined three parameters, mean (μ\\mu), power parameter (ρ\\rho), dispersion (ϕ\\phi), relationship variance parameters defined σ2=ϕμρ\\sigma^2=\\phi\\mu^\\rho. ρ\\rho bounded 1 2, Tweedie positive, continuous distribution can handle zero values, thus allowing naturally handle survey time series one zero observations. Values ρ\\rho bounds special cases Tweedie distribution, ρ\\rho=1 equivalent Poisson distribution ρ\\rho=2 gamma distribution. rema, Tweedie ρ\\rho constrained 1 2 using logit-transformation. Similar normal distribution, observation error variances treated known Tweedie distribution REMA model. Using biomass survey observation example, dispersion biomass observation strata jj year tt derived ϕBt,j=σBt,j2(B̂t,j)ρ. \\phi_{B_{t,j}}=\\frac{\\sigma_{B_{t,j}}^2}{{(\\hat{B}_{t,j})}^\\rho}. result one additional parameter (ρ\\rho) estimated applying alternative distribution. measurement equation Tweedie becomes Bt,j=B̂t,j+ϵB,ϵB∼Twρ(0,ϕBt,j). B_{t,j} = \\hat{B}_{t,j} + \\epsilon_{B}, \\text{} \\epsilon_{B} \\sim Tw_{\\rho}(0, \\phi_{B_{t,j}}). σBt,j\\sigma_{B_{t,j}} undefined Bt,jB_{t,j}=0, zeros assumed CV=1.5. assumption can explored user rema, example, user wanted define CV=2.0 zeros, define zeros = list(assumption = 'tweedie', options_tweedie = list(zeros_cv = 2.0)) within prepare_rema_input() function. reproducible example using Tweedie distribution observation error outlined Strategies handling zero biomass observations vignette. use non-shortspine thornyhead component Bering Sea/Aleutian Islands rockfish stock, unique 13 38 bottom trawl survey biomass estimates zeros. Tweedie performs well case, exploration revealed Tweedie models can slow run often converge, especially instances observation error variance estimates low. possible development estimating additional observation error variance (e.g. τ2{\\tau}^2 previous section) alleviate convergence errors. However, interim, Tweedie distribution rema considered experimental, recommend viable alternative tactical assessments time.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"experimental-one-step-ahead-osa-residuals","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Experimental One-Step Ahead (OSA) residuals","title":"REMA model equations","text":"use one-step ahead (OSA) residuals, also referred forecast residuals prediction errors, crucial validation state-space models like REMA (Thygesen et al. 2017). Instead comparing observed expected values time step, OSA residuals use forecasted values based previous observations (.e. exclude observed value current time step prediction). Traditional residuals (e.g. Pearson’s residuals) inappropriate state-space models, process error variance may -estimated cases model mispecified, thus leading artificially small residuals (see Section 3 Thygesen et al. 2017 example). Methods calculating OSA residuals implemented TMB TMB::oneStepPredict() function. methods straight forward implement TMB, computationally demanding, validity (accuracy) OSA residuals may vary situation method used. Methods implement OSA residuals rema development considered experimental. Currently, OSA residuals implemented using method = \"cdf\" option TMB::oneStepPredict(), benefit speeding residual calculations saving copies one-step cumulative density function data point thus reducing number calculations function evaluation. OSA residuals appear calculated correctly REMA models, occasionally NaN values residuals returned, especially cases small measurement errors. One potential cause error may lie initialization state process within REMA model explored future versions package.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"references","dir":"Articles","previous_headings":"","what":"References","title":"REMA model equations","text":"Echave, K. B., P. J. F. Hulson, K. . Siwicke. 2020. Assessment Thornyhead stock complex Gulf Alaska. : Stock assessment fishery evaluation report groundfish resources Gulf Alaska projected 2021. North Pacific Fishery Management Council, 605 W. 4th. Avenue, Suite 306, Anchorage, AK 99501. Fournier D. ., H. J. Skaug , J. Ancheta , J. Ianelli , . Magnusson, M. N. Maunder , . Nielsen, J. Sibert. 2012. AD Model Builder: using automatic differentiation statistical inference highly parameterized complex nonlinear models, Optimization Methods Software, 27:2, 33-249, DOI: 10.1080/10556788.2011.597854 Marlow, N. . 1967. normal limit theorem power sums independent normal random variables. Bell System Technical Journal. 46 (9): 2081–2089. doi:10.1002/j.1538-7305.1967.tb04244.x. Methot, R. D. Wetzel, C. R. 2013. Stock Synthesis: biological statistical framework fish stock assessment fishery management. Fisheries Research, 142: 86-99. https://doi.org/10.1016/j.fishres.2012.10.012 Skaug H. D. Fournier. 2013. Random effects AD Model Builder: ADMB-RE User Guide. http://ftp.admb-project.org/admb-11.2pre/admbre-11.2pre.pdf Skaug, H. J. D. . Fournier. 2006. Automatic approximation marginal likelihood non-Gaussian hierarchical models. Computational Statistics & Data Analysis, 51(2): 699–709. https://doi.org/10.1016/j.csda.2006.03.005 Thygesen, U. H., Albertsen, C. M., Berg, C. W., Kristensen, K., Nielsen, . 2017. Validation ecological state space models using Laplace approximation. Environ Ecol Stat 24, 317–339. https://doi.org/10.1007/s10651-017-0372-4 Zheng, J., M. S. M. Siddeek. 2020. Bristol Bay red king crab stock assessment fall 2020. : Stock assessment fishery evaluation report king tanner crab fisheries Bering Sea Aleutian Islands regions. North Pacific Fishery Management Council, 605 W 4th Ave, Suite 306 Anchorage, AK 99501.","code":""},{"path":"https://afsc-assessments.github.io/rema/authors.html","id":null,"dir":"","previous_headings":"","what":"Authors","title":"Authors and Citation","text":"Jane Sullivan. Author, maintainer. Laurinne Balstad. Author, contributor. Cole Monnahan. Contributor. Pete Hulson. Contributor.","code":""},{"path":"https://afsc-assessments.github.io/rema/authors.html","id":"citation","dir":"","previous_headings":"","what":"Citation","title":"Authors and Citation","text":"Sullivan J, Balstad L (2024). rema: generalized framework fit random effects (RE) model, state-space random walk model developed Alaska Fisheries Science Center (AFSC) apportionment biomass estimation groundfish crab stocks.. R package version 1.2.0, https://afsc-assessments.github.io/rema/, https://github.com/JaneSullivan-NOAA/rema.","code":"@Manual{, title = {rema: A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.}, author = {Jane Sullivan and Laurinne Balstad}, year = {2024}, note = {R package version 1.2.0, https://afsc-assessments.github.io/rema/}, url = {https://github.com/JaneSullivan-NOAA/rema}, }"},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"a-generalized-random-effects-model-for-fitting-survey-biomass-estimates-with-the-options-of-including-multiple-survey-strata-and-an-additional-survey-index","dir":"","previous_headings":"","what":"A generalized Random Effects model for fitting survey biomass estimates with the options of including Multiple survey strata and an Additional survey index","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"Sullivan, J., C. Monnahan, P. Hulson, J. Ianelli, J. Thorson, . Havron. 2022a. REMA: consensus version random effects model ABC apportionment Tier 4/5 assessments. Plan Team Report, Joint Groundfish Plan Teams, North Pacific Fishery Management Council. 605 W 4th Ave, Suite 306 Anchorage, AK 99501. Available Oct 2022 Joint GPT e-Agenda. Living documentation REMA model equations available REMA website:","code":""},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"background","dir":"","previous_headings":"","what":"Background","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"random effects (RE) model developed North Pacific Fisheries Management Council (NPFMC) Groundfish Plan Team’s Survey Averaging working group used Alaska Fisheries Science Center (AFSC) since 2013 estimate biomass data-limited groundfish crab stock assessments, apportion catch among management areas. RE model estimates biomass series random effects underlying state dynamics modeled random walk (Oct 2013 Joint GPT minutes). several versions original single-strata (univariate) RE single-survey, multiple-strata (multivariate; REM) models currently use AFSC. RE REM models extended 2017 include addition second relative index abundance (e.g., NMFS longline survey Relative Population Numbers) inform biomass trend estimation additional scaling parameters (REMA; Hulson et al. 2021). Although RE, REM, REMA models share underlying state-space random walk dynamics, Monnahan et al. (2021) found inconsistencies treatment zero biomass observations, use prior penalty process error parameter, weighting observation error likelihood components. purpose rema R library provide unified random effects model flexible enough accommodate multitude use cases AFSC. presented , REMA model generalized include RE REM allows users quickly develop evaluate complex models alternative parameterizations. variable names updated original versions increase interpretability transparency, model rewritten Template Model Builder (TMB; Kristensen et al. 2016). Additionally, rema provides flexible framework following analyses: Compare bridge existing ADMB models currently used Tiers 4 5 stock assessments Tier 3 apportionment REMA; Visualize multiple model results simultaneously conduct model selection appropriate using Akaike Information Criteria (AIC); Analyze alternative approaches handling zero biomass observations, including treating zeros missing values failed surveys, adding small constant zero values manually defining associated coefficient variation values, exploring experimental option model observations using Tweedie distribution, positive, continuous distribution accommodates zeros; Perform model validation using appropriate diagnostic tools latent variable models [e.g., reporting convergence criteria conducting one-step-ahead (OSA) residual analysis]; Evaluate use priors process error scaling parameter parameters, along alternative likelihood penalties. structure, naming conventions, functions, documentation rema inspired modeled Woods Hole Assessment Model (WHAM; Miller Stock 2020), open-source, state-space age-structured assessment model R package.","code":""},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"installation","dir":"","previous_headings":"","what":"Installation","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"Use devtools package install rema package github. devtools installed, must first.","code":"# install.packages(\"devtools\") devtools::install_github(\"afsc-assessments/rema\", dependencies = TRUE, build_vignettes = TRUE) # Example R scripts are downloaded when `rema` is installed. Locate them on your computer by running the following commands: (rema_path <- find.package('rema')) (rema_examples <- file.path(rema_path, 'example_scripts')) list.files(rema_examples) # Vignettes library(rema) browseVignettes(\"rema\") vignette(topic = \"rema_equations\") # view technical details offline vignette(topic = \"ex1_basics\") vignette(topic = \"ex2_cpue\") vignette(topic = \"ex3_zeros\") vignette(topic = \"ex4_model_validation\")"},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"the-rema-worflow","dir":"","previous_headings":"","what":"The rema worflow","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"Load rema data. user can read biomass abundance index data file (e.g. csv files), can use rwout.rep report file ADMB version RE model using read_admb_re(). Define model structure assumptions using prepare_rema_input(). function allows users quickly transition single two survey model, specify alternative process error structures, add likelihood penalties priors parameters, evaluate alternative assumptions zero biomass observations. Fit specified REMA model using fit_rema() determine whether model met basic convergence criteria. Extract rema model output clean, consistently formatted data frames using tidy_rema(). user can visualize model output using plot_rema(), quickly format tables report. Compare alternative REMA models conduct model selection using compare_rema_models(). Output function includes table Akaike Information Criteria (AIC) appropriate, figures, tidied data frames. function also accepts model output ADMB version RE model easy comparison past models. Taken together, functions allow R users quickly fit interrogate suite simple statistical models TMB without needing software-specific training expertise.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":null,"dir":"Reference","previous_headings":"","what":"Check convergence of REMA model — check_convergence","title":"Check convergence of REMA model — check_convergence","text":"Access quick convergence checks `TMB` `nlminb`. Function modified wham::check_convergence.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Check convergence of REMA model — check_convergence","text":"","code":"check_convergence(mod, ret = FALSE, f = \"\")"},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Check convergence of REMA model — check_convergence","text":"mod output fit_rema ret T/F, return list? Default = FALSE, just prints console","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Check convergence of REMA model — check_convergence","text":"list least first three components: $convergence stats::nlminb, \"0 indicates successful convergence nlminb\" $maxgr Max absolute gradient value, `max(abs(mod$gr(mod$opt$par)))` $maxgr_par Name parameter max gradient $is_sdrep TMB::sdreport performed model, indicates whether performed without error $na_sdrep TMB::sdreport performed without error model, indicates () components diagonal inverted hessian returned NA","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Check convergence of REMA model — check_convergence","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":null,"dir":"Reference","previous_headings":"","what":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"check_estimability calculates matrix second-derivatives marginal likelihood w.r.t. fixed effects, see linear combinations estimable (.e. uniquely estimated conditional upon model structure available data, e.g., resulting likelihood ridge singular, non-invertable Hessian matrix)","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"","code":"check_estimability(obj, h)"},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"obj compiled object h optional argument containing pre-computed Hessian matrix","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"tagged list hessian message","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"Takes list REMA models fit_rema, returns list ggplot2 objects plotted saved, list tidy_rema data.frames, AIC values.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"","code":"compare_rema_models( rema_models, admb_re = NULL, save = FALSE, filetype = \"png\", path = NULL, xlab = NULL, biomass_ylab = \"Biomass\", cpue_ylab = \"CPUE\" )"},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"rema_models list REMA models compared. REMA model list list object output fit_rema admb_re list ADMB RE model input/output read_admb_re. Accepts single list, list multiple ADMB RE models. admb_re provided, AIC calculations conducted. save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. filetype (optional) character string; type figure file. Default = 'png'. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. xlab (optional) label x-axis biomass CPUE plots (e.g. 'Year'). Default = NULL. biomass_ylab (optional) label y-axis biomass plots (e.g. 'Biomass (t)'). Default = 'Biomass'. cpue_ylab (optional) label y-axis CPUE plots (e.g. 'Relative Population Number'). Default = 'CPUE'.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"list following items: $output list tidied dataframes include parameter estimates, biomass optional CPUE data, REMA model predictions model compared. Results given variable included applicable comparison models. example, CPUE fit one model another, compare$output$cpue_by_strata) return informational message instead dataframe. See tidy_rema information. $plots ggplot2 figure objects compare$output data. $aic dataframe Akaike Information Criteria (AIC) values. output underlying models fit data.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":null,"dir":"Reference","previous_headings":"","what":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"extract_fixed extracts best previous value fixed effects, way works mixed fixed effect models","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"","code":"extract_fixed(obj)"},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"obj, compiled object","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"vector fixed-effect estimates","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":null,"dir":"Reference","previous_headings":"","what":"Fit REMA model — fit_rema","title":"Fit REMA model — fit_rema","text":"Fits compiled REMA model using TMB::MakeADFun stats::nlminb. Source code documentation modified wham::fit_wham.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Fit REMA model — fit_rema","text":"","code":"fit_rema( input, n.newton = 0, do.sdrep = TRUE, model = NULL, do.check = FALSE, MakeADFun.silent = TRUE, do.fit = TRUE, save.sdrep = TRUE )"},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Fit REMA model — fit_rema","text":"input Named list output prepare_rema_input, includes following components needed fit model using TMB::MakeADFun: $data Data, list data objects model fitting specification (e.g., user-defined pentalties, index pointers, etc.). required input MakeADFun. $par Parameters, list random fixed effects parameter objects. required input MakeADFun. $map Map, mechanism collecting fixing parameters TMB. input MakeADFun. $random Character vector defining parameters treat random effects. input MakeADFun. $model_name Character, name model, e.g. \"GOA shortraker LLS depth strata\". Useful model comparison. n.newton integer, number additional Newton steps optimization. option currently needed, passed fit_tmb. Default = 0. .sdrep T/F, calculate standard deviations model parameters? See sdreport. Default = TRUE. model (optional), previously fit rema model. .check T/F, check model parameters identifiable? Passed fit_tmb. Runs internal function check_estimability, originally provided https://github.com/kaskr/TMB_contrib_R/TMBhelper. Default = TRUE. MakeADFun.silent T/F, Passed silent argument TMB::MakeADFun. Default = TRUE. .fit T/F, fit model using fit_tmb. Default = TRUE. save.sdrep T/F, save full TMB::sdreport object? FALSE, save summary.sdreport reduce model object file size. Default = TRUE.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Fit REMA model — fit_rema","text":"fit TMB model additional output specified: $rep List derived quantity estimates (e.g. estimated biomass) $sdrep Parameter estimates (standard errors .sdrep = TRUE)","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Fit REMA model — fit_rema","text":"Future development: Implement one-step-ahead (OSA) residuals evaluating model goodness--fit TMB::oneStepPredict). OSA residuals appropriate standard residuals models random effects (Thygeson et al. (2017). See wham example OSA implementation additional OSA residual options (e.g. full Gaussian approximation instead (default) generic method using osa.opts=list(method=\"fullGaussian\").","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Fit REMA model — fit_rema","text":"","code":"if (FALSE) { # \\dontrun{ # place holder for example code } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":null,"dir":"Reference","previous_headings":"","what":"Fit TMB model using nlminb — fit_tmb","title":"Fit TMB model using nlminb — fit_tmb","text":"Runs optimization TMB model using stats::nlminb. specified, takes additional Newton steps calculates standard deviations. Internal function called fit_rema. Source code documentation modified wham::fit_tmb.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Fit TMB model using nlminb — fit_tmb","text":"","code":"fit_tmb( model, n.newton = 0, do.sdrep = TRUE, do.check = FALSE, save.sdrep = FALSE )"},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Fit TMB model using nlminb — fit_tmb","text":"model Output TMB::MakeADFun. n.newton Integer, number additional Newton steps optimization. Default = 0. .sdrep T/F, calculate standard deviations model parameters? See TMB::sdreport. Default = TRUE. .check T/F, check model parameters identifiable? Runs internal check_estimability, originally provided https://github.com/kaskr/TMB_contrib_R/TMBhelper. Default = TRUE. save.sdrep T/F, save full TMB::sdreport object? FALSE, save summary.sdreport) reduce model object file size. Default = FALSE.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Fit TMB model using nlminb — fit_tmb","text":"model, appends following: model$opt Output stats::nlminb model$date System date model$dir Current working directory model$rep model$report() model$TMB_version Version TMB installed model$parList List parameters, model$env$parList() model$final_gradient Final gradient, model$gr() model$sdrep Estimated standard deviations model parameters, TMB::sdreport summary.sdreport)","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":null,"dir":"Reference","previous_headings":"","what":"Get one-step-head (OSA) — get_osa_residuals","title":"Get one-step-head (OSA) — get_osa_residuals","text":"Takes rema model output fit_rema returns OSA residuals calculated using TMB::oneStepPredict accompanying residual analysis plots. IMPORTANT: OSA residuals work users implementing Tweedie distribution.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Get one-step-head (OSA) — get_osa_residuals","text":"","code":"get_osa_residuals( rema_model, options = list(method = \"fullGaussian\", parallel = TRUE) )"},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Get one-step-head (OSA) — get_osa_residuals","text":"rema_model list output fit_rema, includes model results also inputs. note OSA residual calculations rema_model$input$osa object, data.frame containing data observations fit model residuals associated . options list options calculating OSA residuals, passed TMB::oneStepPredict. Default: options = list(method = \"fullGaussian\", parallel = TRUE). Alternative methods include \"cdf\", \"oneStepGeneric\", \"oneStepGaussianOffMode\", \"oneStepGaussian\".","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Get one-step-head (OSA) — get_osa_residuals","text":"list tidied data.frames containing biomass CPUE survey residuals accompanying data, well QQ-plot, histogram residuals, plots residuals~year residuals~fitted values strata biomass CPUE survey.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Get one-step-head (OSA) — get_osa_residuals","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/pipe.html","id":null,"dir":"Reference","previous_headings":"","what":"Pipe function — %>%","title":"Pipe function — %>%","text":"Allows use pipe function, %>%","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"Takes list output tidy_rema returns list ggplot2 objects plotted saved.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"","code":"plot_extra_cv( tidy_rema, save = FALSE, filetype = \"png\", path = NULL, xlab = NULL, biomass_ylab = \"Biomass\", cpue_ylab = \"CPUE\" )"},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"tidy_rema list output tidy_extra_cv, includes inputs, model results, confidence intervals total observation error (fixed + estimated) save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. filetype (optional) character string; type figure file. Default = 'png'. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. xlab (optional) label x-axis biomass CPUE plots (e.g. 'Year'). Default = NULL. biomass_ylab (optional) label y-axis biomass plots (e.g. 'Biomass (t)'). Default = 'Biomass'. cpue_ylab (optional) label y-axis CPUE plots (e.g. 'Relative Population Number'). Default = 'CPUE'.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"list ggplot2 plots character string messages data. Except parameter estimates, objects output tidy_rema outputted function.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot survey data and model output — plot_rema","title":"Plot survey data and model output — plot_rema","text":"Takes list output tidy_rema returns list ggplot2 objects plotted saved.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot survey data and model output — plot_rema","text":"","code":"plot_rema( tidy_rema, save = FALSE, filetype = \"png\", path = NULL, xlab = NULL, biomass_ylab = \"Biomass\", cpue_ylab = \"CPUE\" )"},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot survey data and model output — plot_rema","text":"tidy_rema list output tidy_rema, includes model results also inputs save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. filetype (optional) character string; type figure file. Default = 'png'. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. xlab (optional) label x-axis biomass CPUE plots (e.g. 'Year'). Default = NULL. biomass_ylab (optional) label y-axis biomass plots (e.g. 'Biomass (t)'). Default = 'Biomass'. cpue_ylab (optional) label y-axis CPUE plots (e.g. 'Relative Population Number'). Default = 'CPUE'.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot survey data and model output — plot_rema","text":"list ggplot2 plots character string messages data. Except parameter estimates, objects output tidy_rema outputted function.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot survey data and model output — plot_rema","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":null,"dir":"Reference","previous_headings":"","what":"Prepare input data and parameters for REMA model — prepare_rema_input","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"data read R (either manually .csv data file using read_admb_re), function prepares data parameter settings fit_rema. model can set run single survey mode one strata, multi-survey mode, uses additional relative abundance index (.e. cpue) inform predicted biomass. optional inputs described related CPUE survey data scaling parameter q, cpue_dat options_q used multi_survey = 1. function structure documentation modeled wham::prepare_wham_input.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"","code":"prepare_rema_input( model_name = \"REMA for unnamed stock\", multi_survey = 0, admb_re = NULL, biomass_dat = NULL, cpue_dat = NULL, sum_cpue_index = FALSE, start_year = NULL, end_year = NULL, wt_biomass = NULL, wt_cpue = NULL, PE_options = NULL, q_options = NULL, zeros = NULL, extra_biomass_cv = NULL, extra_cpue_cv = NULL )"},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"model_name name stock identifier REMA model multi_survey switch run model single multi-survey mode. 0 (default) = single survey, 1 = multi-survey. admb_re list object returned read_admb_re.R, includes biomass survey data (admb_re$biomass_dat), optional cpue survey data (admb_re$cpue_dat), years model predictions (admb_re$model_yrs), model predictions log biomass strata correct format input REMA (admb_re$init_log_biomass_pred). supplied, user need enter biomass_dat cpue_dat. biomass_dat data.frame biomass survey data long format following columns: strata character; survey name, survey region, management unit, depth strata. Note user must include column even one survey strata year integer; survey year. Note user needs include years observations (.e. need supply NULL NA values missing survey years) biomass numeric; biomass estimate/observation (e.g. bottom trawl survey biomass mt). default, biomass == 0, value treated NA (.e., failed survey). user wants make assumptions zeros (e.g. adding small constant), must define data manually. cv numeric; coefficient variation (CV) biomass estimate (.e. sd(biomass)/biomass) cpue_dat (optional) data.frame relative abundance index (.e. cpue) data long format following columns: strata character; survey name, survey region, management unit, depth strata (note user must include column even one survey strata) year integer; survey year. Note user needs include years observations (.e. need supply NULL NA values missing survey years) cpue numeric; cpue estimate/observation (e.g. longline survey cpue relative population number); default, cpue == 0, value treated NA (.e., failed survey). user wants make assumptions zeros (e.g. adding small constant), must define data manually. cv numeric; coefficient variation (CV) cpue estimate (.e. sd(cpue)/cpue) sum_cpue_index T/F 1/0, CPUE survey index able summed across strata get total CPUE survey index? example, Longline survey relative population numbers (RPNs) summable longline survey numbers per hachi (CPUE) . Default = FALSE. start_year (optional) integer value specifying start year estimation model; admb_re supplied, value defaults start_year = min(admb_re$model_yrs); admb_re supplied, value defaults first year either biomass_dat cpue_dat end_year (optional) integer value specifying last year estimation model; admb_re supplied, value defaults end_year = max(admb_re$model_yrs); admb_re supplied, value defaults last year either biomass_dat cpue_dat wt_biomass (optional) multiplier biomass survey data component negative log likelihood. example, nll = wt_biomass * nll. Defaults wt_biomass = 1 wt_cpue (optional) multiplier CPUE survey data component negative log likelihood. example, nll = wt_cpue * nll. Defaults wt_cpue = 1 PE_options (optional) customize implementation process error (PE) parameters, including options share PE across biomass survey strata, change starting values, fix parameters, add penalties priors (see details) q_options (optional) customize implementation scaling parameters (q), including options define q biomass cpue survey cpue strata, change starting values, fix parameters, add penalties priors (see details). used multi_survey = 1 zeros (optional) define assumptions treat zero biomass CPUE observations, including treating zeros NAs, changing zeros small constants fixed CVs, modeling zeros using Tweedie distribution (see details). extra_biomass_cv (optional) estimate additional observation error biomass survey data (see details). default, assumption = \"extra_cv\" estimate one extra CV parameter, regardless number biomass survey strata. extra_cpue_cv (optional) estimate additional observation error CPUE survey data (see details). default, assumption = \"extra_cv\" estimate one extra CV parameter, regardless number CPUE survey strata.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"function returns named list following components: data Named list data, passed TMB::MakeADFun par Named list parameters, passed TMB::MakeADFun map Named list defining optionally collect fix parameters, passed TMB::MakeADFun random Character vector parameters treat random effects, passed TMB::MakeADFun model_name Name stock identifier REMA model biomass_dat tidied long format data.frame biomass survey observations associated CVs strata. data.frame 'complete' include modeled years, missing values treated NAs. Note data.frame differ admb_re$biomass_dat input biomass assumptions zero biomass observations different ADMB model user specifies REMA. user can change assumptions zeros using zeros argument. cpue_dat optional CPUE survey data provided multi_survey = 1, tidied long-format data.frame CPUE survey observations associated CVs strata. data.frame 'complete' include modeled years, missing values treated NAs. Note data.frame differ admb_re$biomass_dat input biomass assumptions zero CPUE observations different ADMB model user specifies REMA. user can change assumptions zeros using zeros argument. optional CPUE survey data provided multi_survey = 0, object NULL.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"PE_options allows user specify options process error (PE) parameters. NULL, default PE specifications used: one PE parameter estimated biomass survey strata, initial values log_PE set 1, penalties priors added. user can modify default PE_options using following list entries: $pointer_PE_biomass index customize assignment PE parameters individual biomass strata. Vector length = number biomass strata, starting index 1 ending number unique PE estimated. example, three biomass survey strata user wants estimate one PE, specify pointer_PE_biomass = c(1, 1, 1). default one unique log_PE estimated unique biomass survey stratum $initial_pars vector initial values log_PE. default initial value log_PE 1. $fix_pars Option fix PE parameters, user specifies index value PE parameter like fix initial value. example, three biomass survey strata, user wants fix log_PE second stratum estimate log_PE first third strata specify fix_pars = c(2) Note option recommended. $penalty_options Warning: following options experimental well-tested. Options penalizing PE likelihood adding prior log_PE include following: \"none\" (default) penalty prior used \"wt\" multiplier PE random effects component negative log likelihood. example, nll = wt * nll, wt = 1.5 specified single value penalty_values argument \"squared_penalty\" implemented earlier version RE.tpl, penalty prevents PE shrinking zero. example, nll = nll + (log_PE + squared_penalty)^2, squared_penalty = 1.5. vector squared_penalty values specified PE penalty_values argument \"normal_prior\" Normal prior log space, nll = nll - dnorm(log_PE, pmu_log_PE, psig_log_PE, 1) pmu_log_PE psig_log_PE specified PE parameter penalty_values argument penalty_values user-defined values penalty_options. penalty type entered follows: \"none\" (default) NULL example, penalty_values = NULL \"wt\" single numeric value. example, penalty_values = 1.5 \"squared_penalty\" vector numeric values length = number PE parameters. example, three PE parameters estimated user wants penalty one, use penalty_values = c(1.5, 1.5, 1.5) \"normal_prior\" vector paired values PE parameter, vector pair prior mean log_PE pmu_log_PE associated standard deviation psig_log_PE. example, three PE parameters estimated user wants normal prior log_PE ~ N(1.0, 0.08), penalty_values = c(c(1.0, 0.08), c(1.0, 0.08), c(1.0, 0.08)) q_options allows user specify options CPUE survey scaling parameters (q). multi_survey = 0 (default), q parameters estimated regardless user defines q_options. multi_survey = 0 q_options = NULL, default q specifications used: one q parameter estimated CPUE survey strata, biomass CPUE surveys assumed share strata definitions (.e., biomass_dat cpue_dat number columns columns represent strata), initial values log_q set 1, penalties priors added. user can modify default q_options using following list entries: $pointer_q_cpue index customize assignment q parameters individual CPUE survey strata. Vector length = number CPUE strata, starting index 1 ending number unique q parameters estimated. example, three CPUE survey strata user wanted estimate one q, specify pointer_q_cpue = c(1, 1, 1). recommended model configuration estimate one log_q CPUE survey stratum. $pointer_biomass_cpue_strata index customize assignment biomass predictions individual CPUE survey strata. Vector length = number biomass survey strata, starting index 1 ending number unique CPUE survey strata. pointer needs defined number biomass CPUE strata equal. pointer_biomass_cpue_strata option allows user calculate predicted biomass CPUE survey strata level scenario biomass survey strata higher resolution CPUE survey strata. example, 3 biomass survey strata represented 2 CPUE survey strata, user may specify pointer_biomass_cpue_strata = c(1, 1, 2). specification assign first 2 biomass strata first CPUE strata, third biomass stratum second CPUE stratum. CPUE data compliment specific biomass stratum, user can populate NAs. example pointer_biomass_cpue_strata = c(1, NA, 3), means CPUE data biomass strata 1 3 2. NOTE: scenario CPUE survey strata biomass survey strata CPUE survey used inform biomass survey trend. error thrown q_options$pointer_biomass_cpue_strata defined biomass CPUE survey strata definitions . $initial_pars vector initial values log_q. default initial value log_q 1. $fix_pars Option fix q parameters, user specifies index value q parameter like fix initial value. example, three CPUE survey strata, user wants fix log_q second stratum estimate log_q first third strata specify fix_pars = c(2) $penalty_options Options penalizing q likelihood adding prior log_q include following: \"none\" (default) penalty prior used \"normal_prior\" Warning, experimental well-tested. Normal prior log space, nll = nll - dnorm(log_q, pmu_log_q, psig_log_q, 1) pmu_log_q psig_log_q specified q parameter penalty_values argument penalty_values user-defined values penalty_options. penalty type entered follows: \"none\" (default) NULL example, penalty_values = NULL \"normal_prior\" vector paired values q parameter, vector pair prior mean log_q pmu_log_q associated standard deviation psig_log_q. example, 2 q parameters estimated user wants normal prior log_q ~ N(1.0, 0.05), penalty_values = c(c(1.0, 0.05), c(1.0, 0.05)) zeros allows user specify options treat zero biomass CPUE survey observations. default zero observations treated NAs warning msg effect returned console. zeros allows user specify non-default zero assumptions using following list entries: $assumption character, name assumption using. three alternatives currently implemented, zeros = list(assumption = c(\"NA\", \"small_constant\", \"tweedie\"). \"NA\" default; option assumes zero estimates failed surveys removes . \"small_constant\" ad hoc method fixed value added zero assumed CV. default, small constant = 0.0001 CV value entered user data. user can change assumed value CV using options_small_constant. \"tweedie\" uses Tweedie assumed error distribution survey data, allows zeros. alternative estimates one additional power parameter. assumed CV zero biomass zero cpue survey observations defaults 1.5. user can change assumed CV, change initial values inverse logit transformed power parameter, fix initial values using options_tweedie. $options_small_constant vector length two numeric values. first value small constant add zero observation, second user-defined coefficient value. user can specify small value use input CV specifying NA second value. E.g., 'options_small_constant = c(0.0001, NA)'. $options_tweedie list entries control initial fixed values Tweedie parameters. Currently, argument accepts following entries: extra_biomass_cv allows user specify options estimating additional CV parameter (log_tau_biomass source code, estimated log-space) biomass survey observations. extra_biomass_cv = NULL (default), extra CV estimated. user can modify default extra_biomass_cv options using following list entries: $assumption string identifying assumption used biomass survey observations. Options include \"none\" (default extra CV estimated) \"extra_cv\". assumption = \"extra_cv\", default one extra CV estimated, regardless many biomass strata defined. extra_biomass_cv NULL, user must define appropriate assumption. $pointer_extra_biomass_cv index customize assignment extra CV parameters individual biomass survey strata. Vector length = number biomass strata, starting index 1 ending number unique extra CV parameters estimated. three biomass survey strata user wanted estimate extra CV per stratum, specify pointer_extra_biomass_cv = c(1, 2, 3). default, one additional parameter estimated, regardless many strata defined (.e. pointer_extra_biomass_cv = c(1, 1, 1)). $initial_pars vector initial values extra biomass log_tau_biomass. default initial value log_tau_biomass log(1e-7) (approximately 0 arithmetic scale). $fix_pars Option fix extra biomass CV parameters, user specifies index value parameter like fix initial value. example, three biomass survey strata defined pointer_extra_biomass_cv, user wants fix log_tau_biomass second stratum estimate log_tau_biomass first third strata specify fix_pars = c(2). extra_cpue_cv allows user specify options estimating additional CV parameter (log_tau_cpue source code, estimated log-space) cpue survey observations. extra_cpue_cv = NULL (default), extra CV estimated. user can modify default extra_cpue_cv options using following list entries: $assumption string identifying assumption used cpue survey observations. Options include \"none\" (default extra CV estimated) \"extra_cv\". assumption = \"extra_cv\", default one extra CV estimated, regardless many cpue strata defined. extra_cpue_cv NULL, user must define appropriate assumption. $pointer_extra_cpue_cv index customize assignment extra CV parameters individual cpue survey strata. Vector length = number cpue strata, starting index 1 ending number unique extra CV parameters estimated. three cpue survey strata user wanted estimate extra CV per stratum, specify pointer_extra_cpue_cv = c(1, 2, 3). default, one additional parameter estimated, regardless many strata defined (.e. pointer_extra_cpue_cv = c(1, 1, 1)). $initial_pars vector initial values extra cpue log_tau_cpue. default initial value log_tau_cpue log(1e-7) (approximately 0 arithmetic scale). $fix_pars Option fix extra cpue CV parameters, user specifies index value parameter like fix initial value. example, three cpue survey strata defined pointer_extra_cpue_cv, user wants fix log_tau_cpue second stratum estimate log_tau_cpue first third strata specify fix_pars = c(2).","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"","code":"if (FALSE) { # \\dontrun{ # place holder for example code } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":null,"dir":"Reference","previous_headings":"","what":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"Read report file ADMB version RE model (rwout.rep) convert long format survey data estimates CVs input REMA.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"","code":"read_admb_re( filename, model_name = \"Unnamed ADMB RE model\", biomass_strata_names = NULL, cpue_strata_names = NULL )"},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"filename name ADMB output file read (e.g. rwout.rep) model_name (optional) Name stock identifier ADMB version RE model. Defaults 'ADMB RE' biomass_strata_names (optional) vector character names corresponding names biomass survey strata. Vector order columns srv_est rwout.rep cpue_strata_names (optional) vector character names corresponding names CPUE survey strata. Vector order columns srv_est_LL rwout.rep version ADMB RE model accepts additional survey index","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"object type \"list\" biomass optional cpue survey data long format, initial parameter values log_biomass_pred (random effects matrix), ready input REMA list following items: $biomass_dat dataframe biomass survey data strata, year, biomass estimates, CVs. Note CVs back-transformed natural space. $cpue_dat Optional dataframe CPUE survey data strata, year, CPUE estimates, CVs. Note CVs back-transformed natural space. $model_yrs Vector prediction years. $init_log_biomass_pred Matrix initial parameter values log_biomass_pred (random effects matrix), ready input REMA. $admb_re_results list ADMB RE model results ready comparison REMA models using compare_rema_models(). User beware... many, many versions RE.tpl existence individual variances may cause errors output.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"","code":"if (FALSE) { # \\dontrun{ # place holder for example code } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":null,"dir":"Reference","previous_headings":"","what":"Read ADMB .rep file and return an R object of type 'list' — read_rep","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"Code modified original function provided Steve Martell, D'Arcy N. Webber called read_admb_re","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"","code":"read_rep(fn)"},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"fn full path name ADMB output file read","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"object type \"list\" ADMB outputs therein","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"","code":"if (FALSE) { # \\dontrun{ read_rep(fn = 'inst/example_data/goasr.rep') } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":null,"dir":"Reference","previous_headings":"","what":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"Takes list output tidy_rema returns list enhanced versions biomass_by_strata cpue_by_strata appropriate. enhanced dataframes include three new columns, tot_log_obs_cv, tot_obs_lci, tot_obs_uci, represent combined log-space standard error associated confidence intervals include assumed estimated additional observation error.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"","code":"tidy_extra_cv(tidy_rema, save = FALSE, path = NULL, alpha_ci = 0.05)"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"tidy_rema list output tidy_rema, includes model results also inputs save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. alpha_ci (optional) significance level generating confidence intervals. Default = 0.05","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"list following items: $parameter_estimates data.frame fixed effects parameters REMA (e.g. log_PE log_q) standard errors confidence intervals transformed log space natural space ease interpretation. $biomass_by_strata tidy, long format data.frame model predicted observed biomass biomass survey strata. data.frame now enhanced new columns include log-space standard error associated confidence intervals account additional estimated observation error. $cpue_by_strata tidy, long format data.frame model predicted observed CPUE CPUE survey strata. data.frame now enhanced new columns include log-space standard error associated confidence intervals account additional estimated observation error. REMA run multi-survey mode, CPUE data provided, explanatory character string instructions fitting CPUE data returned. $biomass_by_cpue_strata tidy, long format data.frame model predicted biomass CPUE survey strata. Note observed/summed biomass observations returned case missing values one stratum another within given year. output reserved instances number biomass strata exceeds CPUE survey strata, user wants visualize predicted biomass resolution CPUE predictions. scenarios, character string returned explaining special use case object. $total_predicted_biomass tidy, long format data.frame total model predicted biomass summed across biomass survey strata. one stratum used (.e. univariate RE), predicted values output$biomass_by_strata. $total_predicted_cpue tidy, long format data.frame total model predicted CPUE summed across CPUE survey strata. one stratum used (.e. univariate RE), predicted values output$cpue_by_strata. CPUE survey index provided defined summable prepare_rema_input(), character string returned explaining change using 'sum_cpue_index' ?prepare_rema_input appropriate.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":null,"dir":"Reference","previous_headings":"","what":"Tidy REMA model output — tidy_rema","title":"Tidy REMA model output — tidy_rema","text":"Takes outputs fit_rema, returns named list tidied data.frames include parameter estimates standard errors, derived variables model. information \"tidy\" data, please see Wickham 2014. code modified wham::par_tables_fun.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Tidy REMA model output — tidy_rema","text":"","code":"tidy_rema(rema_model, save = FALSE, path = NULL, alpha_ci = 0.05)"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Tidy REMA model output — tidy_rema","text":"rema_model list output fit_rema, includes model results also inputs save (optional) logical (T/F) save output data.frames csvs path. Default = FALSE. YET IMPLEMENTED. path (optional) directory path location csvs saved save = TRUE. YET IMPLEMENTED. alpha_ci (optional) significance level generating confidence intervals. Default = 0.05","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Tidy REMA model output — tidy_rema","text":"list following items: $parameter_estimates data.frame fixed effects parameters REMA (e.g. log_PE log_q) standard errors confidence intervals transformed log space natural space ease interpretation. $biomass_by_strata tidy, long format data.frame model predicted observed biomass biomass survey strata. $cpue_by_strata tidy, long format data.frame model predicted observed CPUE CPUE survey strata. REMA run multi-survey mode, CPUE data provided, explanatory character string instructions fitting CPUE data returned. $biomass_by_cpue_strata tidy, long format data.frame model predicted biomass CPUE survey strata. Note observed/summed biomass observations returned case missing values one stratum another within given year. output reserved instances number biomass strata exceeds CPUE survey strata, user wants visualize predicted biomass resolution CPUE predictions. scenarios, character string returned explaining special use case object. $total_predicted_biomass tidy, long format data.frame total model predicted biomass summed across biomass survey strata. one stratum used (.e. univariate RE), predicted values output$biomass_by_strata. $total_predicted_cpue tidy, long format data.frame total model predicted CPUE summed across CPUE survey strata. one stratum used (.e. univariate RE), predicted values output$cpue_by_strata. CPUE survey index provided defined summable prepare_rema_input(), character string returned explaining change using 'sum_cpue_index' ?prepare_rema_input appropriate.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Tidy REMA model output — tidy_rema","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"}] +[{"path":"https://afsc-assessments.github.io/rema/articles/ex1_basics.html","id":"the-rema-workflow","dir":"Articles","previous_headings":"","what":"The rema workflow:","title":"REMA basics","text":"Load rema data. user can read biomass abundance index data file (e.g. .csv file), can use rwout.rep report file ADMB version RE model using read_admb_re(). Specify model structure assumptions using prepare_rema_input(). function allows users quickly transition single two survey model, specify alternative process error structures, add likelihood penalties priors parameters, evaluate alternative assumptions zero biomass observations. Fit specified REMA model using fit_rema() determine whether model met basic convergence criteria (e.g., Hessian positive definite, maximum gradient component approximately equal zero). Extract rema model output clean, consistently formatted data frames using tidy_rema(). user can visualize model output using plot_rema(), quickly format tables report. Compare alternative REMA models conduct model selection using compare_rema_models(). Output function includes table Akaike Information Criteria (AIC) appropriate, figures, tidied data frames. function also accepts model output ADMB version RE model easy comparison past models. Taken together, functions allow R users quickly fit interrogate suite simple statistical models TMB without needing software-specific training expertise.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex1_basics.html","id":"example-1-univariate-random-effects-re-model-with-a-single-survey-and-stratum","dir":"Articles","previous_headings":"","what":"Example 1: Univariate random effects (RE) model with a single survey and stratum","title":"REMA basics","text":"example uses Aleutian Islands shortraker (aisr_rwout.rep) fits NMFS bottom trawl survey estimates. Read rwout.rep, custom report file generated ADMB version random effects model. Prepare REMA model inputs using admb_re data object. input$data, input$par, input$map, input$random used TMB fit REMA model. objects used functions process model results conduct residual analyses. Data alternatively entered prepare_rema_input() using biomass_dat argument: Fit REMA model Get tidied model output plot results Note 95% confidence intervals observations (.e., obs_lci obs_uci output$biomass_by_strata based assumption normality log-space; therefore, asymmetric arithmetic scale. Compare ADMB RE model results","code":"# ?read_admb_re admb_re <- read_admb_re(filename = 'aisr_rwout.rep', # optional label for the single biomass survey stratum biomass_strata_names = 'Aleutians Islands', model_name = 'ADMB: AI shortraker') names(admb_re) #> [1] \"biomass_dat\" \"cpue_dat\" \"model_yrs\" #> [4] \"init_log_biomass_pred\" \"admb_re_results\" kable(admb_re$biomass_dat) # ?prepare_rema_input input <- prepare_rema_input(model_name = 'TMB: AI shortraker', admb_re = admb_re) names(input) #> [1] \"data\" \"par\" \"map\" \"random\" \"model_name\" #> [6] \"osa\" \"biomass_dat\" # (not run) input <- prepare_rema_input(model_name = 'tmb_rema_aisr', biomass_dat = admb_re$biomass_dat) # ?fit_rema m <- fit_rema(input) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.90e-11 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model # ?tidy_rema output <- tidy_rema(rema_model = m) kable(output$parameter_estimates) # estimated fixed effects parameters kable(head(output$biomass_by_strata, 5)) # data.frame of predicted and observed biomass by stratum # ?plot_rema plots <- plot_rema(tidy_rema = output, biomass_ylab = 'Biomass (t)') # optional y-axis label plots$biomass_by_strata # ?compare_rema_models compare <- compare_rema_models(rema_models = list(m), admb_re = admb_re, biomass_ylab = 'Biomass (t)') compare$plots$biomass_by_strata"},{"path":"https://afsc-assessments.github.io/rema/articles/ex1_basics.html","id":"example-2-multivariate-random-effects-model-rem-with-a-single-survey-and-multiple-strata","dir":"Articles","previous_headings":"","what":"Example 2: Multivariate random effects model (REM) with a single survey and multiple strata","title":"REMA basics","text":"example uses Bering Sea Aleutian Islands shortspine thornyhead (bsaisst_rwout.rep) fits NMFS EBS slope AI bottom trawl survey estimates. original ADMB model single, shared process error (PE) three strata. example demonstrates define alternative PE structures, u perform model selection using AIC. can use ggplot2 functions modify formatting plots: One primary uses RE model apportioning catch management area. tidy_rema() function output includes table proportions predicted biomass strata: can compare TMB model results original ADMB model. Note different confidence intervals ADMB TMB versions. ADMB code uses Marlow method sum variances biomass log-space, whereas TMB version applies standard Delta method. can easily fit alternative REMA model strata-specific PE parameters, compare two models using AIC. case see single PE model lowest AIC value. Note ADMB models currently compared REMA models using AIC, therefore omit admb_re = admb_re argument compare_rema_models() function call :","code":"admb_re <- read_admb_re(filename = 'bsaisst_rwout.rep', biomass_strata_names = c('AI survey', 'EBS slope survey', 'S. Bering Sea (AI survey)'), model_name = 'ADMB: single PE') # the original ADMB model shared PE across all strata input1 <- prepare_rema_input(model_name = 'TMB: single PE', admb_re = admb_re, PE_options = list(pointer_PE_biomass = c(1, 1, 1))) m1 <- fit_rema(input1) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 9.81e-10 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model output <- tidy_rema(rema_model = m1) kable(output$parameter_estimates) # estimated fixed effects parameters kable(head(output$biomass_by_strata, 5)) # data.frame of predicted and observed biomass by stratum kable(head(output$total_predicted_biomass, 5)) # combined/total predicted biomass plots <- plot_rema(tidy_rema = output, biomass_ylab = 'Biomass (t)') plots$biomass_by_strata + facet_wrap(~strata, ncol = 1, scales = 'free_y') kable(tail(output$proportion_biomass_by_strata, 3)) # figure: # plots$proportion_biomass_by_strata compare <- compare_rema_models(rema_models = list(m1), admb_re = admb_re, biomass_ylab = 'Biomass (t)') # Note different confidence intervals between the ADMB version (Marlow method) and the TMB version (Delta method) compare$plots$biomass_by_strata + facet_wrap(~strata, ncol = 1, scales = 'free_y') compare$plots$total_predicted_biomass + ggplot2::ggtitle('BSAI Shortspine thornyhead') # REMA defaults to strata-specific parameters, which could be explicitly defined # as: PE_options = list(pointer_PE_biomass = c(1, 2, 3)) input2 <- prepare_rema_input(model_name = 'TMB: strata-specific PE', admb_re = admb_re) m2 <- fit_rema(input2) #> Model runtime: 0.3 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 2.31e-07 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model compare <- compare_rema_models(rema_models = list(m1, m2), biomass_ylab = 'Biomass (t)') compare$plots$biomass_by_strata + facet_wrap(~strata, ncol = 1, scales = 'free_y') kable(compare$aic) kable(compare$output$parameter_estimates)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-1-fit-to-the-biomass-and-cpue-survey-through-the-estimation-of-a-scaling-parameter-q","dir":"Articles","previous_headings":"","what":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","title":"Fitting to an additional CPUE survey","text":"","code":"# read in the data using the report file from the original ADMB model admb_re <- read_admb_re(filename = 'goasr_rwout.rep', biomass_strata_names = c('CGOA', 'EGOA', 'WGOA'), cpue_strata_names = c('CGOA', 'EGOA', 'WGOA'), model_name = 'M0: ADMB status quo (invalid)') input <- prepare_rema_input(model_name = 'M1: TMB status quo', # fit to both biomass and CPUE survey data multi_survey = 1, admb_re = admb_re, # the RPW index is summable across strata sum_cpue_index = TRUE, # likelihood weight wt_cpue = 0.5, # one process error parameters (log_PE) estimated PE_options = list(pointer_PE_biomass = c(1, 1, 1)), # three scaling parameters (log_q) estimated, indexed as # follows for each biomass survey stratum: q_options = list(pointer_biomass_cpue_strata = c(1, 2, 3))) m1 <- fit_rema(input) #> Model runtime: 0.2 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 4.12e-04 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model output <- tidy_rema(m1) output$parameter_estimates #> model_name parameter estimate std_err lci #> 1 M1: TMB status quo process_error 0.1696964 0.03914612 0.1079729 #> 2 M1: TMB status quo scaling_parameter_q 0.4123880 0.04474590 0.3333857 #> 3 M1: TMB status quo scaling_parameter_q 1.2127809 0.10874662 1.0173198 #> 4 M1: TMB status quo scaling_parameter_q 2.1164194 0.30469400 1.5960887 #> uci #> 1 0.2667047 #> 2 0.5101115 #> 3 1.4457965 #> 4 2.8063798 plots <- plot_rema(output, biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') cowplot::plot_grid(plots$biomass_by_strata, plots$cpue_by_strata, ncol = 1)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"a-comparison-of-the-admb-and-tmb-models","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"A comparison of the ADMB and TMB models","title":"Fitting to an additional CPUE survey","text":"compare fits status quo ADMB model identically structured correctly-specified TMB model, see large differences model predictions. particular PE variance different, ADMB model producing PE seven times PE TMB model: M0: ADMB status quo (invalid) PE = 0.0233 M1: TMB status quo PE = 0.1697","code":"compare <- compare_rema_models(rema_models = list(m1), admb_re = admb_re, biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') compare$plots$biomass_by_strata + theme(legend.position = 'top')"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-2-estimating-additional-observation-error-in-the-two-surveys","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Model 2: Estimating additional observation error in the two surveys","title":"Fitting to an additional CPUE survey","text":"outcome TMB model problematic expect shortraker rockfish, notably long-lived species, exhibit extreme inter-annual variability. part, finding attributed relatively low observation error survey data. following section explore alternative model estimates additional biomass CPUE survey observation error. use AIC conduct model selection. following figure shows fits data ; however, used tidy_extra_cv() plot_extra_cv() functions obtain total 95% confidence intervals survey biomass CPUE observations (.e., fixed + estimated observation error). Note 95% confidence intervals observations (.e., obs_lci/obs_uci tot_obs_lci/tot_obs_uci output$biomass_by_strata output$cpue_by_strata based assumption normality log-space; therefore, asymmetric arithmetic scale. figure , error bars around survey observations show 95% confidence intervals based assumed design-based estimates (bold whiskers) 95% confidence intervals based total observation error (design-based estimates CV + additional estimated CV; whiskers).","code":"input2 <- prepare_rema_input(model_name = 'M2: TMB extra obs error', # fit to both biomass and CPUE survey data multi_survey = 1, admb_re = admb_re, # the RPW index is summable across strata sum_cpue_index = TRUE, # likelihood weight wt_cpue = 0.5, # one process error parameters (log_PE) estimated PE_options = list(pointer_PE_biomass = c(1, 1, 1)), # three scaling parameters (log_q) estimated, indexed as # follows for each biomass survey stratum: q_options = list(pointer_biomass_cpue_strata = c(1, 2, 3)), # estimate additional obs error on the biomass # survey (one additional CV shared across all 3 # strata) extra_biomass_cv = list(assumption = 'extra_cv'), # estimate additional obs error on the CPUE # survey (one additional CV shared across all 3 # strata) extra_cpue_cv = list(assumption = 'extra_cv')) m2 <- fit_rema(input2) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.58e-05 #> Max gradient parameter: log_q #> TMB:sdreport() was performed successfully for this model out2 <- tidy_rema(m2) # adds new columns with total observation error and 95% confidence intervals to # the biomass_by_strata and cpue_by_strata out2 <- tidy_extra_cv(out2) cvplots <- plot_extra_cv(out2) cowplot::plot_grid(cvplots$biomass_by_strata, cvplots$cpue_by_strata, ncol = 1)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-comparison","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Model comparison","title":"Fitting to an additional CPUE survey","text":"comparison three models shows estimation additional observation error substantially reduces PE variance, resulting smoother trajectory biomass predictions.","code":"compare <- compare_rema_models(rema_models = list(m1, m2), admb_re = admb_re, biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') compare$plots$biomass_by_strata + theme(legend.position = 'top')"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"model-selection","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Model selection","title":"Fitting to an additional CPUE survey","text":"Next, can compare fits biomass CPUE survey indices alternative TMB models remove admb_re object comparison function demonstrated . examination parameter estimates shows M2 (extra observation error) estimates much lower PE, approximately 50% M1 PE. Using AIC, find M2 superior fit M1 statistical support additional parameters (note example intended illustrative model fitting model selection process exhaustive).","code":"compare <- compare_rema_models(rema_models = list(m1, m2), biomass_ylab = 'Biomass (t)', cpue_ylab = 'Relative Population Weights') cowplot::plot_grid(compare$plots$biomass_by_strata + theme(legend.position = 'top'), compare$plots$cpue_by_strata + theme(legend.position = 'none'), ncol = 1, rel_widths = c(0.65, 0.35)) kable(compare$output$parameter_estimates) kable(compare$aic)"},{"path":"https://afsc-assessments.github.io/rema/articles/ex2_cpue.html","id":"apportionment-results","dir":"Articles","previous_headings":"Model 1: Fit to the biomass and CPUE survey through the estimation of a scaling parameter qq","what":"Apportionment results","title":"Fitting to an additional CPUE survey","text":"proportion prediction biomass strata (management area example), often used apportionment Acceptable Biological Catch estimates. Results compare_rema_models() function include table proportioned biomass management area purpose.","code":"compare$output$proportion_biomass_by_strata %>% filter(year == max(year)) %>% kable() # compare$plots$proportion_biomass_by_strata # optional figure"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"model-1-zeros-as-nas","dir":"Articles","previous_headings":"","what":"Model 1: Zeros as NAs","title":"Strategies for handling zero biomass observations","text":"","code":"nonsst <- read.csv('ebsshelf_orox_biomass.csv') input1 <- prepare_rema_input(model_name = 'M1: Zeros as NAs', biomass_dat = nonsst, zeros = list(assumption = 'NA')) m1 <- fit_rema(input1) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 4.23e-10 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"model-2-add-a-small-constant","dir":"Articles","previous_headings":"","what":"Model 2: Add a small constant","title":"Strategies for handling zero biomass observations","text":"","code":"input2 <- prepare_rema_input(model_name = 'M2: Small constant=0.01, CV=1.5', biomass_dat = nonsst, zeros = list(assumption = 'small_constant', # values: 1) small constant, 2) assumed CV options_small_constant = c(0.01, 1.5))) m2 <- fit_rema(input2) #> Model runtime: 0.1 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.33e-07 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"model-3-tweedie-distribution","dir":"Articles","previous_headings":"","what":"Model 3: Tweedie distribution","title":"Strategies for handling zero biomass observations","text":"","code":"input3 <- prepare_rema_input(model_name = 'M3: Tweedie', biomass_dat = nonsst, zeros = list(assumption = 'tweedie')) m3 <- fit_rema(input3) #> Model runtime: 0.2 seconds #> stats::nlminb thinks the model has converged: mod$opt$convergence == 0 #> Maximum gradient component: 1.72e-10 #> Max gradient parameter: log_PE #> TMB:sdreport() was performed successfully for this model"},{"path":"https://afsc-assessments.github.io/rema/articles/ex3_zeros.html","id":"compare-model-results","dir":"Articles","previous_headings":"","what":"Compare model results","title":"Strategies for handling zero biomass observations","text":"Models 1 (Zeros NAs) 3 (Tweedie) similar results. Tweedie appears perform well case, can slow run suffer convergence issues, especially observation errors relatively small. Tweedie alternative considered experimental issues can resolved.","code":"compare <- compare_rema_models(rema_models = list(m1, m2, m3)) compare$plots$biomass_by_strata"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"background","dir":"Articles","previous_headings":"","what":"Background","title":"REMA model validation","text":"Fisheries stock assessments moving towards state-space estimation, boasts range benefits, including separation estimation observation process error, elegant framework handling missing data, high degree flexibility respect model architecture inclusion different data types, potential improved projections predictive skill (Aeberhard et al., 2018). flexibility relative ease fitting state-space models means can increase complexity dimensionality rapidly. easy fit, state-space models often challenging validate, even simple models can suffer estimation issues result biased inference (Auger‐Méthé et al., 2016). North Pacific Fishery Management Council (NPFMC), “random effects model” (REMA) far common state-space model used fishery management (Sullivan et al., 2022). REMA state-space random walk model can customized estimate multiple process errors, fit additional abundance index, estimate additional observation error. REMA used estimate biomass within Tier 5 groundfish Tier 4 crab stock assessments apportion Acceptable Biological Catches (ABCs) management area many stocks. Despite high impact model within North Pacific fishery management process, standard model validation practices REMA. goals: Apply established state-space model validation techniques real life REMA examples. Create REMA model validation guide clear descriptions method motivation use, explanation interpret model validation results, reproducible code run operational REMA models NPFMC. Provide preliminary recommendations REMA users reviewers model validation methods relevant informative REMA. Interested otherwise busy readers encouraged skip “care failed model diagnostics?” “Recommendations” sections. Table 1 Figures 6-8 show examples REMA model failing model validation criteria.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"a-testing-framework-for-rema-model-validation","dir":"Articles","previous_headings":"","what":"A testing framework for REMA model validation","title":"REMA model validation","text":"analysis, testing (1) REMA working expected (.e., code implemented correctly parameters estimated without bias), (2) assumptions made parameterizing estimating REMA valid, including assumptions related random effects error structure. use two example stocks: Aleutian Islands Pacific cod (AI Pcod; Spies et al., 2023): simple, univariate case, fits single time series AI bottom trawl survey estimates one process error Gulf Alaska Thornyhead rockfish (GOA thornyhead; Echave et al., 2022): complex, multi-strata multi-survey case, estimates multiple process errors, scaling parameter longline survey, two additional observation errors GOA bottom trawl longline surveys models appear converge per standard REMA TMB diagnostics (small maximum gradient, invertible Hessian), fit data reasonably, model fits shown Figure 1. stocks span levels complexity NPFMC REMA models; testing REMA across range complexity helps ensure model model assumptions valid variety realistic, management-relevant cases. important note model validation mean model “right” “correct.” Model validation help analyst model selection (except identify models might function properly), ensure model makes “better” predictions. Rather, model validation way ensure model operating expected, without introducing bias violating statistical assumptions upon model based (Auger‐Méthé et al., 2016; Auger‐Méthé et al., 2021), data reasonably come model (Thygesen et al., 2017). key questions aim answer model validation include following: model perform expected, introduce bias? Using simulation self-check, test model coded correctly able recover known parameters. plausible data generated model? Using one-step ahead (OSA) residuals, appropriate residual type state-space models, test model assumptions look trends residuals reveal characteristics dynamics data aren’t adequately captured model. normality assumptions made estimating random effects via Laplace approximation accurate? comparing posterior distribution fixed effects without Laplace approximation, test whether distribution random effects normal thus accuracy Laplace approximation. allows us test accuracy fixed effects estimates uncertainties. parameters unique non-redundant? Checking correlation parameters helps us identify parameters redundant .","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"prepare-data-for-each-stock-here","dir":"Articles","previous_headings":"","what":"Prepare data for each stock here","title":"REMA model validation","text":"First, run model stock. AI Pcod model (pcod_mod) uses single process error describe stock time (one parameter). GOA thornyhead rockfish model (thrn_mod) uses data two surveys (bottom trawl longline survey) describe stock time (six fixed effects parameters: three process errors, one scaling parameter, additional observation error biomass survey, additional observation error longline survey). code , show specify fit models within rema framework plot fits data. Figure 1 shows fits data AI Pcod GOA Thornyhead. particular, GOA Thornyhead fits show trade biomass longline survey data. Figure 1. Fits AI Pcod (top panel; gold) GOA Thornyhead (bottom panels; teal) models.","code":"set.seed(415) # for reproducibility, use a seed # first, get the p-cod data set up pcod_bio_dat <- read.csv(\"ai_pcod_2022_biomass_dat.csv\") pcod_input <- prepare_rema_input(model_name = \"p_cod\", biomass_dat = pcod_bio_dat, # one strata PE_options = list(pointer_PE_biomass = 1) ) # run the model pcod_mod <- fit_rema(pcod_input) # next, get the thornyhead data set up thrn_bio_dat <- read.csv(\"goa_thornyhead_2022_biomass_dat.csv\") thrn_cpue_dat <-read.csv(\"goa_thornyhead_2022_cpue_dat.csv\") thrn_input <- prepare_rema_input(model_name = 'thrnhead_rockfish', multi_survey = TRUE, biomass_dat = thrn_bio_dat, cpue_dat = thrn_cpue_dat, # RPWs are a summable/area-weighted effort index sum_cpue_index = TRUE, # three process error parameters (log_PE) estimated, # indexed as follows for each biomass survey stratum # (shared within an area across depths): PE_options = list(pointer_PE_biomass = c(1, 1, 1, 2, 2, 2, 3, 3, 3)), # scaling parameter options: q_options = list( # longline survey strata (n=3) indexed as follows for the # biomass strata (n=9) pointer_biomass_cpue_strata = c(1, 1, 1, 2, 2, 2, 3, 3, 3), # one scaling parameters (log_q) estimated, shared # over all three LLS strata pointer_q_cpue = c(1, 1, 1)), # estimate extra trawl survey observation error extra_biomass_cv = list(assumption = 'extra_cv'), # estimate extra longline survey observation error extra_cpue_cv = list(assumption = 'extra_cv') ) # run the model thrn_mod <- fit_rema(thrn_input) # tidy output and plot fitted data tidy_pcod <- tidy_rema(pcod_mod) p1 <- plot_rema(tidy_pcod)$biomass_by_strata + ggtitle(label = \"Model Fits to the AI Pcod Data\", subtitle = \"Trawl Survey Biomass Strata\") + geom_ribbon(aes(ymin = pred_lci, ymax = pred_uci), col = 'goldenrod', fill = 'goldenrod', alpha = 0.4) + geom_line() + geom_point(aes(x = year, y = obs)) + geom_errorbar(aes(x = year, ymin = obs_lci, ymax = obs_uci)) p1 ggsave(paste0(\"vignettes/ex4_pcod_fits.png\"), width = 6, height = 4, units = \"in\", bg = \"white\") tidy_thrn <- tidy_rema(thrn_mod) p2 <- plot_rema(tidy_thrn, biomass_ylab = 'Biomass (t)')$biomass_by_strata + ggtitle(label = \"Model Fits to the GOA Thornyhead Data\", subtitle = \"Trawl Survey Biomass Strata\") + geom_ribbon(aes(ymin = pred_lci, ymax = pred_uci), col = \"#21918c\", fill = \"#21918c\", alpha = 0.4) + geom_line() + geom_point(aes(x = year, y = obs)) + geom_errorbar(aes(x = year, ymin = obs_lci, ymax = obs_uci)) p2 p3 <- plot_rema(tidy_thrn, cpue_ylab = \"RPW\")$cpue_by_strata + ggtitle(label = NULL, subtitle = \"Longline Survey Relative Population Weight (RPW) Strata\") + geom_ribbon(aes(ymin = pred_lci, ymax = pred_uci), col = \"#21918c\", fill = \"#21918c\", alpha = 0.4) + geom_line() + geom_point(aes(x = year, y = obs)) + geom_errorbar(aes(x = year, ymin = obs_lci, ymax = obs_uci)) p3 cowplot::plot_grid(p2, p3, rel_heights = c(0.7, 0.3), ncol = 1) ggsave(paste0(\"vignettes/ex4_thrn_fits.png\"), width = 11, height = 10, units = \"in\", bg = \"white\")"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"simulation-self-test-can-we-recover-parameters-without-bias","dir":"Articles","previous_headings":"","what":"1. Simulation self-test: Can we recover parameters without bias?","title":"REMA model validation","text":"Simulation testing ensures model properly coded performs consistently without introducing bias (Auger‐Méthé et al., 2021; Gimenez et al., 2004). First, use REMA model estimate parameters (e.g., process error variance) real data. Next, use estimated parameters “true values” simulate new latent states (random effects) data conditioned states using REMA equations. use REMA model re-estimate model parameters (“recovered parameters”) calculate relative error (RE; .e., (true-estimated values)/true value*100)) parameters simulation replicate (N=500). AI Pcod model fitted bottom trawl data simulated data set biomass trajectory; GOA Thornyhead model fitted bottom trawl longline survey data, simulation data set includes biomass trajectories nine strata longline survey relative population weights (RPWs) three strata. Models fail simulation testing recovered parameters biased (RE ≠ 0), deviate consistently true values used simulate data (span recovered parameters large): can indicate model coding error, non-identifiable, redundant parameters, biased (.e., kind model misspecification). following function runs simulation self-test REMA model. Run self-test AI Pcod GOA Thornyhead examine results: Results simulation self-test 1 500 simulation replicates converge AI Pcod model, suggests high degree model stability. 11 500 simulation replicates GOA Thornyhead model converge. Figure 2 top row within panel shows distribution resulting recovered parameter estimates models based simulations, horizontal line median estimated value black dot true value (.e., parameters estimates real data sets).  demonstrated boxplots bottom rows Figure 2, parameters models recovered well simulation testing, low median relative error. models instances long negative tails log standard deviation process error (log_PE; Figure 2 top row within panel), suggesting PE small tended towards zero simulations (recall parameters estimated log space, negative number log space close zero natural space). also case extra observation error biomass survey (log_tau_biomass) GOA Thornyhead model (Figure 2, top row within lower panel).  outlier replicates, occur optimizer used rema (nlminb()) returns “NA/NaN function evaluation” error, help shed light replicates failed converge. hood, optimizer pushing PE towards zero (PE=0 taking global mean biomass time series), violating central assumption REMA model biomass underlying trend (PE > 0, thus PE can log-transformed). negative log-scale values extreme GOA Thornyhead model, raise concern potential model misspecification -parameterization. example, given trade-process observation error REMA (including additional observation dampens biomass trend estimated, pushing PE towards zero), possible estimation additional observation error biomass survey (log_tau_biomass) might justified. Figure 2. Simulation self-test results AI Pcod (top panels; gold) GOA Thornyhead (bottom panels; teal). upper row panel gives distribution recovered parameters, dot giving “true value” used simulate data. lower row panel gives relative error recovered parameters, zero line indicating simulation returns “true value,” box plot line giving median recovered parameter’s relative errors, whiskers dots boxplots giving spread recovered parameter’s relative errors.","code":"sim_test <- function(mod_name, replicates, cpue) { # storage things re_est <- matrix(NA, replicates, length(mod_name$par)) # parameter estimates # go through the model suppressMessages(for(i in 1:replicates) { sim <- mod_name$simulate(complete = TRUE) # simulates the data # simulated biomass observations: tmp_biomass <- matrix(data = exp(sim$log_biomass_obs), ncol = ncol(mod_name$input$data$biomass_obs)) colnames(tmp_biomass) <- colnames(mod_name$data$biomass_obs) # simulated cpue observations, when applicable: if (cpue) {tmp_cpue <- matrix(data = sim$cpue_obs, ncol = ncol(mod_name$input$data$cpue_obs)) colnames(tmp_cpue) <- colnames(mod_name$data$cpue_obs)} # set up new data for input newinput <- mod_name$input newinput$data$biomass_obs <- tmp_biomass # biomass data if (cpue) {newinput$data$cpue_obs <- tmp_cpue} # cpue data # create \"obsvec\" which is used internally in cpp file as the observation # vector for all observation (log biomass + cpue) in the likelihood # functions (and required for OSA residuals). note the transpose t() needed # to get these matrices in the correct order (by row instead of by col) -- # note obsvec is masked from users normally in prepare_rema_input() newinput$data$obsvec <- t(sim$log_biomass_obs)[!is.na(t(sim$log_biomass_obs))] if (cpue) {newinput$data$obsvec <- c(t(sim$log_biomass_obs)[!is.na(t(sim$log_biomass_obs))], t(sim$log_cpue_obs)[!is.na(t(sim$log_cpue_obs))])} # refit model mod_new <- fit_rema(newinput, do.sdrep = FALSE) # add parameter estimates to matrix if(mod_new$opt$convergence == 0) { re_est[i, ] <- mod_new$env$last.par[1:length(mod_name$par)] } else { re_est[i, ] <- rep(NA, length(mod_name$par)) } }) re_est <- as.data.frame(re_est); re_est$type <- rep(\"recovered\") return(re_est) } # run for pcod and prep data frame # run simulation testing par_ests <- sim_test(mod_name = pcod_mod, replicates = 500, cpue = FALSE) # get data frame with simulation values for each parameter n_not_converged_pcod <- length(which(is.na(par_ests[,1]))); n_pcod <- length(par_ests[,1]) prop_converged_pcod <- 1-n_not_converged_pcod/n_pcod mod_par_ests <- data.frame(\"log_PE1\" = pcod_mod$env$last.par[1], type = \"model\") names(par_ests) <- names(mod_par_ests) # rename for ease pcod_par_ests <- rbind(mod_par_ests, par_ests) # recovered and model in one data frame pcod_par_ests$sp <- rep(\"AI Pcod\") # run for thorny and prep data frame -- same process # run simulation testing par_ests <- sim_test(mod_name = thrn_mod, replicates = 500, cpue = TRUE) # note some warnings # sometimes spits out: In stats::nlminb(model$par, model$fn, model$gr, control = # list(iter.max = 1000,...: NA/NaN function evaluation n_not_converged_thrn <- length(which(is.na(par_ests[,1]))); n_thrn <- length(par_ests[,1]) prop_converged_thrn <- 1-n_not_converged_thrn/n_thrn nrow(par_ests) # get data frame with simulation values for each parameter mod_par_ests <- data.frame(\"log_PE1\" = thrn_mod$env$last.par[1], \"log_PE2\" = thrn_mod$env$last.par[2], \"log_PE3\" = thrn_mod$env$last.par[3], \"log_q\" = thrn_mod$env$last.par[4], \"log_tau_biomass\" = thrn_mod$env$last.par[5], \"log_tau_cpue\" = thrn_mod$env$last.par[6], type = \"model\") names(par_ests) <- names(mod_par_ests) # rename for ease thrn_par_ests <- rbind(mod_par_ests, par_ests) # recovered and model parameters in one data frame thrn_par_ests$sp <- rep(\"GOA Thornyhead\") # reogranize data pcod_sim <- pcod_par_ests %>% pivot_longer(1, names_to = \"parameter\") thrn_sim <- thrn_par_ests %>% pivot_longer(1:6, names_to = \"parameter\") sim_dat <- rbind(pcod_sim, thrn_sim) # Relative Error: ((om-em)/om) sim_re <- sim_dat %>% pivot_wider(id_cols = c(\"sp\", \"parameter\"), names_from = type, values_from = value) %>% unnest(cols = c(model, recovered)) %>% mutate(RE = (model-recovered)/model*100) %>% group_by(sp, parameter) %>% mutate(label = paste0(\"Median RE=\", formatC(median(RE, na.rm = TRUE), format = \"f\", digits = 1), \"%\")) %>% suppressWarnings() # plot distribution of simulated parameter estimates plot_sim <- function(sim_dat, plot_title, fill_col = \"#21918c\") { ggplot(NULL, aes(parameter, value)) + # add distribution of recovered parameters geom_violin(data = sim_dat %>% filter(type == \"recovered\"), fill = fill_col, alpha = 0.6, draw_quantiles = 0.5) + # add \"true values\" from original model geom_point(data = sim_dat %>% filter(type == \"model\"), size = 2, col = \"black\") + facet_wrap(~ parameter, scales = \"free\", nrow = 1) + labs(x = NULL, y = \"Parameter estimate\", title = plot_title, subtitle = \"Distribution of parameters estimates (median=horizontal line, true value=point)\") + scale_x_discrete(labels = NULL, breaks = NULL) } # plot relative error plot_re <- function(sim_re, fill_col = \"#21918c\") { ggplot(NULL, aes(parameter, RE)) + # add distribution of recovered parameters geom_boxplot(data = sim_re, alpha = 0.6, fill = fill_col, na.rm = TRUE, outlier.size = 0.8) + geom_hline(yintercept = 0) + # separate by parameter facet_wrap(~ parameter+label, scales = \"free\", nrow = 1) + #, space = \"free\") + labs(x = NULL, y = \"Relative error (%)\", subtitle = \"Distribution of relative error (RE; i.e., (true-estimated values)/true value*100)\") + scale_x_discrete(labels = NULL, breaks = NULL) } p1pcod <- plot_sim(pcod_sim, \"AI Pcod Simulation\", fill_col = \"goldenrod\") p1thrn <- plot_sim(thrn_sim, \"GOA Thornyhead Simulation\") p2pcod <- plot_re(sim_re %>% filter(sp == \"AI Pcod\"), fill_col = \"goldenrod\") p2thrn <- plot_re(sim_re %>% filter(sp == \"GOA Thornyhead\")) cowplot::plot_grid(p1pcod, p2pcod, ncol = 1) ggsave(paste0(\"vignettes/ex4_sim_pcod.png\"), width = 6.5, height = 7, units = \"in\", bg = \"white\") cowplot::plot_grid(p1thrn, p2thrn, ncol = 1) ggsave(paste0(\"vignettes/ex4_sim_thrn.png\"), width = 11, height = 7, units = \"in\", bg = \"white\")"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"residual-analysis","dir":"Articles","previous_headings":"","what":"2. Residual analysis","title":"REMA model validation","text":"Residuals used test underlying assumptions structure error distributions model. example, linear regression, residuals independent, normal, constant variance. Traditional residuals (e.g., Pearson’s residuals) inappropriate state-space models like REMA, random effects induce correlations predicted data residuals longer independent. Additionally, process error variance may overestimated cases model mis-specified, thus leading artificially small residuals (see Section 3 Thygesen et al. 2017 example). appropriate residual type validation state-space models one-step ahead (OSA) residuals. Instead comparing observed expected values time step, OSA residuals use forecasted values based previous observations (.e., exclude observed value current time step prediction). way, OSA residuals account non-normality correlation residuals among years.  correctly specified model, resultant OSA residuals independent identically distributed (..d.) standard normal distribution N(0,1)N(0,1). can tested normal QQ plot, theoretical quantile values standard normal distribution plotted x-axis, corresponding empirical quantile values OSA residuals plotted y-axis. OSA residuals standard normally distributed, points fall 0/1 reference line, standard deviation normalized residuals (SDNR) statistic close 1. residuals standard normally distributed (SDNR far 1), points deviate reference line. case model includes multiple strata surveys, OSA residuals normally distributed within across strata surveys, however small sample sizes (e.g., within single stratum) may sufficient statistical power identify model misfit.  addition normality OSA residuals, residuals random independent (.e., correlated year). Standard approaches autocorrelation function (ACF) plots residual runs tests (e.g., Wald-Wolfowitz test) inappropriate many REMA applications missing years data survey time series, instead, directly plot residuals years survey time series. Visual patterns extreme outliers (greater 3 less -3 N(0,1) distribution) residuals can indicate structure data (.e., temporal correlation) adequately captured model.  Methods calculating OSA residuals REMA implemented rema::get_osa_residuals() using TMB::oneStepPredict() function fullGaussian method using lognormal error distributions. rema::get_osa_residuals() returns tidied dataframes observations, REMA model predictions, OSA residuals biomass CPUE survey data appropriate, along QQ residual-time series diagnostic plots.  Results residual analysis normal QQ plot AI Pcod (Figure 3, top panel) suggests slight negative skewness residuals (.e., majority points fall 0/1 line), though small sample size makes interpretation challenging. SDNR 0.99 (recall perfect model SDNR=1.00), indicating  assumptions likely met residuals. plot OSA residuals year (Figure 3, bottom panel) shows obvious patterns residuals, suggesting independent; evidence outliers residuals warrant inspection. combined normal QQ plot GOA Thornyhead OSA residuals (Figure 4, top panel) suggests follow normal distribution (points fall along 0/1 line), though evidence slight positive, right skewness (majority points fall 0/1 line). SDNR approximately 1, indicating variance assumptions likely met residuals. QQ plots biomass strata (Figure 4, middle panels) highlight positive skewness may coming (e.g., EGOA 701-1000m, WGOA 0-500 m); however, small sample sizes stratum make interpretation QQ plots difficult. QQ plots CPUE strata (Figure 4, bottom panel) mostly normal, though evidence light tails WGOA stratum. plot OSA residuals year GOA Thornyhead OSA residuals biomass strata (Figure 5, top panels) show patterns residuals, suggesting independent. However, runs residuals CPUE strata (Figure 5, bottom panels), especially CGOA WGOA, might indicate model misspecification misfit. evidence outliers. Figure 3. One-step ahead (OSA) residual plots AI Pcod. top panel gives QQ plot, comparing normal distribution quantiles (x-axis) OSA residual distribution quantiles (y-axis). diagonal line 0-1 line, two quantiles equal, SDNR statistic given upper left plot. bottom panel gives residuals (y-axis) plotted year (x-axis) used check independence. Figure 4. One-step ahead (OSA) QQ plots across (top panel) within strata (middle panel, bottom trawl survey strata; bottom panel, longline survey strata) GOA Thornyhead. panels, colors correspond bottom trawl survey strata. SDNR statistic given across strata QQ plot; within strata plots useful check extreme skew lack statistical power. Figure 5. One-step ahead (OSA) residual plots GOA Thornyhead. residuals (y-axis) annual time step (random effect; x-axis) survey stratum.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"laplace-approximation-are-the-model-assumptions-related-to-random-effects-estimation-reasonable","dir":"Articles","previous_headings":"","what":"3. Laplace approximation: Are the model assumptions related to random effects estimation reasonable?","title":"REMA model validation","text":"rema library developed Template Model Builder (TMB), uses maximum marginal likelihood estimation Laplace approximation efficiently estimate high dimensional, non-linear mixed effects models frequentist framework (Skaug Fournier, 2006; Kristensen et al., 2016). primary assumption models using Laplace approximation random effects follow normal distribution. assumption simplifies complex integrals make likelihood function. Laplace approximation fast accurate normality assumption met; however, true distribution random effects normal, Laplace approximation may introduce bias parameter estimates. can test validity Laplace approximation using Markov chain Monte Carlo (MCMC) sampling tmbstan library, comparing distributions fixed effects parameters MCMC-sampled models without Laplace approximation (Monnahan Kristensen, 2018). , compare distributions using QQ plots two model cases. Laplace approximation reasonably accurate sampling quantiles two models (without Laplace approximation) similar, .e., fall 1:1 line. can also help us identify bias introduced Laplace approximation, example, fixed effects estimates associated estimates uncertainty differ base model (run marginal maximum likelihood estimation using Laplace approximation) MCMC versions without Laplace approximation. MCMC models run 1000 iterations, assuming warmup 500. discussed Results section, iterations increased 5000 warmup 2500 GOA Thornyhead effort improve diagnostics model convergence. function used run two model cases: following code uses output mcmc_comp() make QQ plot compare posterior distributions parameters parameter estimates MCMC models without Laplace approximation: Results testing validity Laplace approximation Figure 6 gives compared quantiles full MCMC testing (x-axis) Laplace approximation (y-axis) AI Pcod GOA Thornyhead models. simpler AI Pcod model (Figure 6, left panel), Laplace approximation full MCMC testing show similar distribution, indicating Laplace approximation reasonable. Additionally estimates fixed effects nearly identical base model MCMC models without Laplace approximation (Table 1). complex GOA Thornyhead model, preliminary model results based 1000 iterations exhibited poor diagnostics process error parameters additional observation error biomass survey (log_tau_biomass). increased number MCMC iterations 5000 effort improve convergence diagnostics; however, log_tau_biomass parameter continued show significant deviations models without Laplace approximation (Figure 6, right panels). comparison fixed effects estimates base models show large differences parameter estimates associated estimates uncertainty (Table 1). indicates Laplace approximation likely inaccurate, investigation model structure might necessary. Table 1. AI Pcod GOA Thornyhead fixed effect estimates 95% confidence credible intervals maximum likelihood MCMC models, respectively (CI). Fixed effects estimates compared base model run marginal maximum likelihood estimation assuming Laplace approximation (“Base_MMLE_Laplace”) MCMC models without Laplace approximation (“MCMC_Laplace” “MCMC_withoutLaplace”, respectively). Rhat statistic provided MCMC models convergence metric (convergence = Rhat = 1). Nonconverged parameters highlighted bold. parameter estimate MCMC models median posterior samples. Figure 6. Results Laplace approximation testing. compare results running model without assuming normality random effects (MCMC, x-axis) running model assuming normality random effects (Laplace approximation, y-axis), quantiles plotted fixed effect model. gray line 1:1 line; points fall gray line indicate quantile value two cases. top panel gives plot AI Pcod model (one fixed effect), bottom panels give plots GOA Thornyhead model (six fixed effects).","code":"# function to (1) run models with and without laplace approximation for # comparison and (2) return posterior data frames and models # function input: model (e.g., AI Pcod or GOA Thornyhead), number of iterations # (samples), number of chains mcmc_comp <- function(mod_name, it_numb, chain_numb) { # set up MCMC chain information it_num <- it_numb chain_num <- chain_numb # mod_name = pcod_mod; it_num = 1000; chain_num = 4 # run model with laplace approximation mod_la <- tmbstan(obj = mod_name, chains = chain_num, init = mod_name$par, laplace = TRUE, iter = it_num) # run model without laplace approximation, i.e., all parameters fully estimated without assumptions of normality mod_mcmc <- tmbstan(obj = mod_name, chains = chain_num, init = mod_name$par, laplace = FALSE, iter = it_num) # posteriors as data frame post_la <- as.data.frame(mod_la); post_la$type <- (\"la\") post_mcmc <- as.data.frame(mod_mcmc); post_mcmc <- post_mcmc[, c(1:length(mod_name$par), dim(post_mcmc)[2])]; post_mcmc$type <- rep(\"mcmc\") # informational things... this is for getting the posterior draws, i.e., to # test traceplots and such post_la$chain <- rep(1:chain_num, each = it_num/2) post_la$iter_num <- rep(1:(it_num/2), chain_num) post_mcmc$chain <- rep(1:chain_num, each = it_num/2) post_mcmc$iter_num <- rep(1:(it_num/2), chain_num) post_draws <- rbind(post_la, post_mcmc) # get quantiles qv <- seq(from = 0, to = 1, by = 0.01) quant_dat <- data.frame(quant = NULL, la = NULL, mcmc = NULL, par = NULL) for (i in 1:(dim(post_la)[2]-4)) { # post_la has type, chains, lp, iteration number column that don\"t count tmp <- data.frame(quant = qv, la = quantile(unlist(post_la[i]), probs = qv), mcmc = quantile(unlist(post_mcmc[i]), probs = qv), par = rep(paste0(\"V\", i))) quant_dat <- rbind(quant_dat, tmp) } return(list(quant_dat, post_draws, mod_la, mod_mcmc)) } # run models pcod_comp <- mcmc_comp(mod_name = pcod_mod, it_numb = 1000, chain_numb = 3) thrn_comp_log_tau <- mcmc_comp(mod_name = thrn_mod, it_numb = 5000, chain_numb = 3) # clean up data frame names by renaming things pcod_qq <- pcod_comp[[1]] pcod_qq$par_name <- rep(\"log_PE1\"); pcod_qq$sp <- rep(\"AI Pcod\") thrn_qq <- thrn_comp_log_tau[[1]] thrn_qq$par_name <- recode(thrn_qq$par, V1 = \"log_PE1\", V2 = \"log_PE2\", V3 = \"log_PE3\", V4 = \"log_q\", V5 = \"log_tau_biomass\", V6 = \"log_tau_cpue\") thrn_qq$sp <- rep(\"GOA Thornyhead\") qq_dat <- rbind(pcod_qq, thrn_qq) # plot plot_laplace_mcmc <- function(qq_dat, plot_title) { ggplot(qq_dat, aes(mcmc, la)) + geom_abline(intercept = 0, slope = 1, col = \"lightgray\") + geom_point() + facet_wrap(~par_name, scales = \"free\", nrow = 2, dir = \"v\") + labs(x = \"MCMC\", y = \"Laplace approx.\", title = plot_title) + theme(plot.title = element_text(size = 12), axis.text = element_text(size = 7), axis.title = element_text(size = 11), strip.text = element_text(size = 11)) } p1 <- plot_laplace_mcmc(pcod_qq, \"AI Pcod\") p2 <- plot_laplace_mcmc(thrn_qq, \"GOA Thornyhead\") cowplot::plot_grid(p1, p2, rel_widths = c(0.4,0.6), ncol = 2) ggsave(paste0(\"vignettes/ex4_mcmc_qq.png\"), width = 11, height = 7, units = \"in\", bg = \"white\") # parameter summary table - compare the marginal max likelihood estimates (MMLE) # with the mcmc estimates. Rhat is the potential scale reduction factor on split # chains (at convergence, Rhat=1) m <- summary(pcod_comp[[3]])$summary[1, c(6,4,8,10), drop = FALSE] psum <- data.frame(Stock = 'AI Pcod', Model = 'MCMC_Laplace', FixedEffect = row.names(m), m, row.names = NULL) m <- summary(pcod_comp[[4]])$summary[1, c(6,4,8,10), drop = FALSE] psum <- psum %>% bind_rows(data.frame(Stock = 'AI Pcod', Model = 'MCMC_withoutLaplace', FixedEffect = row.names(m), m, row.names = NULL)) m <- summary(thrn_comp_log_tau[[3]])$summary[1:6, c(6,4,8,10), drop = FALSE] psum <- psum %>% bind_rows(data.frame(Stock = 'GOA Thornyhead', Model = 'MCMC_Laplace', FixedEffect = row.names(m), m, row.names = NULL)) m <- summary(thrn_comp_log_tau[[4]])$summary[1:6, c(6,4,8,10), drop = FALSE] psum <- psum %>% bind_rows(data.frame(Stock = 'GOA Thornyhead', Model = 'MCMC_withoutLaplace', FixedEffect = row.names(m), m, row.names = NULL)) pars <- as.list(pcod_mod$sdrep, \"Est\")$log_PE sd <- as.list(pcod_mod$sdrep, \"Std\")$log_PE psum <- psum %>% bind_rows(data.frame(Stock = 'AI Pcod', Model = 'Base_MMLE_Laplace', FixedEffect = \"log_PE\", X50. = pars, X2.5. = pars-1.96*sd, X97.5. = pars+1.96*sd, Rhat = NA)) pars <- unlist(as.list(thrn_mod$sdrep, \"Est\")) pars <- pars[names(pars)[grepl(\"log_PE|log_q|log_tau\", names(pars))]] sd <- unlist(as.list(thrn_mod$sdrep, \"Std\")) sd <- sd[names(sd)[grepl(\"log_PE|log_q|log_tau\", names(sd))]] psum <- psum %>% bind_rows(data.frame(Stock = 'GOA Thornyhead', Model = 'Base_MMLE_Laplace', FixedEffect = rownames(summary(thrn_comp_log_tau[[3]])$summary[1:6,,drop=FALSE]), X50. = pars, sd = sd, row.names = NULL) %>% mutate(X2.5. = X50.-1.96*sd, X97.5. = X50.+1.96*sd, Rhat = NA) %>% select(-sd)) # write.csv(psum, 'vignettes/ex4_param_table.csv')"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"parameter-correlation-are-the-model-parameters-identifiable-and-non-redundant","dir":"Articles","previous_headings":"","what":"4. Parameter correlation: Are the model parameters identifiable and non-redundant?","title":"REMA model validation","text":"Parameter redundancy refers idea multiple parameters contribute model way. intuitive case y∼β1+β2+αxy \\sim \\beta_1 + \\beta_2 + \\alpha x, since model estimate many combinations β1\\beta_1 β2\\beta_2 minimize log-likelihood, .e., sampled parameters correlated. reduce redundancy, β1+β2\\beta_1 + \\beta_2 can redefined β0\\beta_0.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"results-of-parameter-correlation-analysis","dir":"Articles","previous_headings":"4. Parameter correlation: Are the model parameters identifiable and non-redundant?","what":"Results of parameter correlation analysis","title":"REMA model validation","text":"Figure 7 shows pairwise correlation matrix posterior draws GOA Thornyhead model (MCMC model Laplace approximation), histograms univariate marginal distributions diagonal scatterplot bivariate distributions diagonal. significant unresolved sampling problems log_tau_biomass parameter, shown non-normal heavily skewed distribution log_tau_biomass. pairwise scatterplots -diagonals show heavy left tails log_tau_biomass causing interactions (.e., correlation) almost parameters model. diagnostic clear indication model mis-specification warrants investigation. Figure 7. Pairs plot GOA Thornyhead MCMC model assuming Laplace approximation. diagonal plots give distribution fixed effect 2500 MCMC draws. -diagonal elements give parameter--parameter correlation, point gives parameter values estimated MCMC simulation.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"is-the-goa-thornyhead-model-really-converged","dir":"Articles","previous_headings":"","what":"5. Is the GOA Thornyhead model really converged?","title":"REMA model validation","text":"Despite poor diagnostics shown vignette GOA Thornyhead model, standard output TMB rema suggests model converged (.e., low maximum gradient component standard error estimates appear reasonable). Within MCMC framework, can look mixing MCMC across multiple chains using diagnostic called traceplot assess convergence. following code shows within rstan library (Stan Development Team, 2024): traceplots models figures , see AI Pcod model MCMC chains fully mixed model converged (top panel). However, GOA Thornyhead model (bottom panel), see estimation additional observation error (log_tau_biomass) causing convergence issues MCMC chains well mixed. Figure 8.  MCMC traceplots AI Pcod (top panel) GOA Thornyhead (bottom panels). subpanel individual parameter, plot gives value parameter (y-axis) simulation draw (x-axis) model chain (color).","code":"p1_mcmc_pcod <- rstan::traceplot(pcod_comp[[3]]) + ggtitle(label = \"Traceplots to assess mixing across Markov chains and convergence\", subtitle = \"AI Pcod\") p1_mcmc_thrn <- rstan::traceplot(thrn_comp_log_tau[[3]], ncol = 3) + ggtitle(label = NULL, subtitle = \"GOA Thornyhead\") cowplot::plot_grid(p1_mcmc_pcod, p1_mcmc_thrn, ncol = 1) ggsave(paste0(\"vignettes/ex4_traceplots.png\"), width = 11, height = 8, units = \"in\", bg = \"white\")"},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"when-should-we-care-about-failed-model-diagnostics","dir":"Articles","previous_headings":"","what":"When should we care about failed model diagnostics?","title":"REMA model validation","text":"REMA model used within NPFMC obtain exploitable biomass estimates Tier 4 crab Tier 5 groundfish method apportion Acceptable Biological Catches among management regions (Sullivan et al., 2022). , primary purpose smooth noisy survey biomass estimates (.e., using fancy average), need care potential red flags raised model validation?  Yes: case GOA Thornyhead model, example, appears issue estimation additional observation error, pointing potential -parameterization model /parameter redundancy. terminal year predictions REMA model quantities directly used management, estimation additional observation error definition increases “smoothness” model predictions, answer likely yes particular case.  Maybe : hypothetical case estimation year predictions show diagnostic concerns, uncertainty estimations show possible diagnostic problems, model validation red flags might less important. uncertainty results REMA generally used fisheries management NPFMC. words, failed model validation criteria primarily concern impact quantities interest model (case, biomass predictions). Model validation, subsequent interpretation, considered light goal model, rather binary pass/fail testing procedure.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"preliminary-recommendations","dir":"Articles","previous_headings":"","what":"Preliminary recommendations","title":"REMA model validation","text":"model validation methods presented tools stock assessment scientists evaluate existing models guide development future models. analysis, found REMA model can recover parameters across range complexity relevant NPFMC Tier 4 Tier 5 stocks. analysis suggests MCMC diagnostics OSA residuals useful diagnostics REMA model. models fail diagnostics, may indicate model complex simplifications considered. highlighted several places document diagnostic features can help make choices appropriate model simplifications (pool process errors, remove additional observation errors, etc.) light REMA model assumptions trade-offs model parameters. results indicate diagnostics important explore existing models using multivariate configurations multiple process errors, multi-survey versions using CPUE indices like longline survey, models estimating additional observation error. model validation may needed existing REMA models, recommend used recommending new models management.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/ex4_model_validation.html","id":"references","dir":"Articles","previous_headings":"","what":"References","title":"REMA model validation","text":"Aeberhard, W.H., Flemming, J.M., Nielsen, . 2018. Review State-Space Models Fisheries Science. Annual Review Statistics Application, 5(1), 215-235. https://doi.org/10.1146/annurev-statistics-031017-100427 Auger-Méthé, M., Field, C., Albertsen, C.M., Derocher, .E., Lewis, M.., Jonsen, .D., Flemming, J.M. 2016. State-space models’ dirty little secrets: even simple linear Gaussian models can estimation problems. Scientific reports, 6(1), 26677. Auger‐Méthé, M., Newman, K., Cole, D., Empacher, F., Gryba, R., King, . Leos-Barajas, V., Flemming, J.M., Nielsen, , Petris, G., Thomas, L. 2021. guide state–space modeling ecological time series. Ecological Monographs, 91(4), e01470. Campbell, D., & Lele, S. 2014. ANOVA test parameter estimability using data cloning application statistical inference dynamic systems. Computational Statistics & Data Analysis, 70, 257-267. Cole, D. J. 2019. Parameter redundancy identifiability hidden Markov models. Metron, 77(2), 105-118. Echave, K.B., Siwicke, K.., Sullivan, J., Ferriss, B., Hulson, P.J.F. 2022. Assessment Thornyhead stock complex Gulf Alaska. Stock assessment fishery evaluation report groundfish fisheries Gulf Alaska. North Pacific Fishery Management Council, 605 W. 4th. Avenue, Suite 306, Anchorage, AK 99501. https://apps-afsc.fisheries.noaa.gov/Plan_Team/2022/GOAthorny.pdf Gabry, J. Mahr, T. 2024. bayesplot: Plotting Bayesian Models. R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry, J., Simpson, D., Vehtari, ., Betancourt, M., Gelman, . 2019. Visualization Bayesian workflow. J. R. Stat. Soc. . (182) 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378. Kristensen, K., Nielsen, ., Berg, C.W., Skaug, H., Bell, B.M. 2016. TMB: Automatic differentiation Laplace approximation. J Stat Softw. 2016; 70(5):21. Epub 2016-04-04. https://doi.org/10.18637/jss.v070.i05 Monnahan, C.C., & Kristensen, K. 2018. -U-turn sampling fast Bayesian inference ADMB TMB: Introducing adnuts tmbstan R packages. PloS one, 13(5), e0197954. Skaug, H.J., Fournier, D.. 2006. Automatic approximation marginal likelihood non-Gaussian hierarchical models. Computational Statistics & Data Analysis, 51(2), 699-709. Spies, ., Barbeaux, S., Hulson, H., Ortiz, . 2023. Assessment Pacifc cod stock Aleutian Islands. Stock assessment fishery evaluation report groundfish fisheries Bering Sea Aleutian Islands. North Pacific Fishery Management Council, 605 W. 4th. Avenue, Suite 306, Anchorage, AK 99501. https://apps-afsc.fisheries.noaa.gov/Plan_Team/2023/AIpcod.pdf Stan Development Team. 2024. RStan: R interface Stan. R package version 2.32.6, https://mc-stan.org/. Sullivan, J., Monnahan, C. Hulson, P. Ianelli, J. Thorson, J. Havron, . 2022. REMA: consensus version random effects model ABC apportionment Tier 4/5 assessments. Plan Team Report, Joint Groundfish Plan Teams, North Pacific Fishery Management Council. 605 W 4th Ave, Suite 306 Anchorage, AK 99501. Available Oct 2022 Joint GPT e-Agenda. Thygesen, U.H., Albertsen, C.M., Berg, C.W., Kristensen, K., & Nielsen, . 2017. Validation ecological state space models using Laplace approximation. Environmental Ecological Statistics, 24(2), 317-339. https://doi.org/10.1007/s10651-017-0372-4","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"base-model-structure-for-a-single-survey-and-stratum","dir":"Articles","previous_headings":"","what":"Base model structure for a single survey and stratum","title":"REMA model equations","text":"base REMA model can represented simple state-space random walk model added noise. observation model comprised index log-transformed annual survey biomass estimates ln(Bt)ln(B_t) associated standard deviations σln(Bt)\\sigma_{ln(B_t)}, σln(Bt)\\sigma_{ln(B_t)} approximated using coefficient variation BtB_t (σBt/Bt\\sigma_{B_t}/B_t) σln(Bt)=ln((σBtBt)2+1). \\sigma_{ln(B_t)}=\\sqrt{ln((\\frac{\\sigma_{B_t}}{B_t})^2+1)}. measurement observation equation, describes relationship observed survey biomass ln(Bt)ln(B_t) latent state variable, population biomass ln(B̂t)ln(\\hat{B}_t), expressed ln(Bt)=ln(B̂t)+ϵB,ϵB∼N(0,σln(Bt)2). ln(B_t) = ln(\\hat{B}_t) + \\epsilon_{B}, \\text{} \\epsilon_{B} \\sim N(0, \\sigma_{ln(B_t)}^2). state equation associated process error variance σPE2\\sigma_{PE}^2 defined ln(B̂t)=ln(B̂t−1)+ηt−1,ηt∼N(0,σPE2), ln(\\hat{B}_{t}) = ln(\\hat{B}_{t-1})+\\eta_{t-1}, \\text{} \\eta_t \\sim N(0, \\sigma_{PE}^2), initial lnB̂1ln\\hat{B}_{1} constrained random walk process. base model, process error variance σPE2\\sigma_{PE}^2 fixed effect parameter estimated, unobserved population biomass ln(B̂t)ln(\\hat{B}_t) estimated series random effects. model fit using marginal maximum likelihood estimation, marginal negative log-likelihood obtained via Laplace approximation (Skaug Fournier 2006). reproducible example fitting univariate (.e., single stratum, single survey) version TMB model comparing results ADMB available REMA basics vignette.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"extending-to-multiple-biomass-survey-strata","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Extending to multiple biomass survey strata","title":"REMA model equations","text":"single survey, single stratum version model can extended include one additional strata jj survey. inclusion multiple strata model advantageous scenarios apportionment biomass among areas desired result. extension assumes correlation observation errors among survey strata, though PE variance can shared estimated independently among strata. ADMB TMB versions REMA model utilize different methods estimating variance total predicted biomass. Therefore, strata-specific estimates predicted biomass associated confidence intervals close identical ADMB TMB versions model, confidence intervals total predicted biomass differ slightly. original ADMB code, Marlow (1967) method applied, total variance σJ2\\sigma_J^2 approximated σJ2=ln(∑e2B̂j+σB̂j2(eσB̂j2−1)(∑eB̂j+σB̂j2/2)2+1). \\sigma_J^2 = ln\\Bigg(\\frac{\\sum{e^{2 \\hat{B}_j + \\sigma_{\\hat{B}_j}^2} (e^{\\sigma_{\\hat{B}_j}^2 - 1})}} {(\\sum{e^{\\hat{B}_j + \\sigma_{\\hat{B}_j}^2/2}})^2} + 1\\Bigg). updated rema package, total variance estimated using standard Delta method can replicated ADMB using sdreport_number TMB using ADREPORT macro. described Monnahan et al. (2021), exploration methods summing log-normal variables research topic potential impacts beyond scope package. reproducible example fitting multivariate version REMA model TMB comparison results ADMB available REMA basics vignette. example, also show Akaike Information Criteria (AIC) can used explore inclusion strata-specific versus single, shared process error variance.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"addition-of-an-auxiliary-catch-per-unit-effort-cpue-survey","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Addition of an auxiliary catch per unit effort (CPUE) survey","title":"REMA model equations","text":"situations auxiliary biomass catch per unit effort (CPUE) survey ItI_{t} associated variance σIt\\sigma_{I_{t}} available, additional scaling parameter qq (estimated log-space) can estimated facilitate inclusion new information biomass predictions. predicted annual CPUE survey index ̂t\\hat{}_{t} calculated ̂t=eq*eB̂t \\hat{}_{t} = e^{q} * e^{\\hat{B}_{t}} CPUE survey observations additional measurement equation likelihood component similar biomass survey: ln()=ln(̂t)+ϵI,ϵI∼N(0,σln()2), ln(I_{t}) = ln(\\hat{}_{t}) + \\epsilon_{}, \\text{} \\epsilon_{} \\sim N(0, \\sigma_{ln(I_t)}^2), default, rema estimates single qq stratum; however, users can alternatively specify shared qq parameters among strata. strata definitions t,biomass CPUE surveys (e.g. biomass estimated finer geographic resolution CPUE), user can define relationship two surveys’ strata using q_options argument prepare_rema_input() function (see q_options pointer_biomass_cpue_strata). However, since auxiliary CPUE survey index related unobserved population biomass level biomass survey strata, REMA model can accommodate scenarios CPUE survey strata coarser resolution equivalent biomass survey strata. reproducible example fitting two-survey version REMA model provided Fitting additional CPUE survey vignette. example also describes error previously used ADMB version REMA model, presents several alternative models estimate additional observation error. topics described detail following two sections.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"the-admb-version-of-rema","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum > Addition of an auxiliary catch per unit effort (CPUE) survey","what":"The ADMB version of REMA","title":"REMA model equations","text":"Separability, implemented SEPARABLE_FUNCTION ADMB automatic detection TMB, increases computational efficiency Laplace approximation breaking complex, multivariate integrals product simpler, univariate integrals (Fournier et al. 2012). use separability particularly relevant state-space models relates efficiency Laplace approximation marginal negative log-likelihood integration random effects joint negative log-likelihood. random effects implementation ADMB, parameters defined inside PARAMETER_SECTION template file used within SEPARABLE_FUNCTION unless passed arguments function (Skaug Fournier 2013). ADMB version REMA model, calculation predicted CPUE [.e. ln(̂t)=ln(q)+ln(B̂)ln(\\hat{}_{t})=ln(q)+ln(\\hat{B})] occurred outside rather inside SEPARABLE_FUNCTION, result, violated rule affecting accuracy Laplace approximation. explored error simplified example REMA model ADMB TMB. able reproduce results TMB version REMA ADMB passing ln(q)ln(q) argument SEPARABLE_FUNCTION performing ln(̂t)ln(\\hat{}_{t}) calculation inside SEPARABLE_FUNCTION. used Markov Chain Monte Carlo (MCMC) methods statistical inference instead MLE, results ADMB versions closely matched correct version REMA. Finally, comparison individual negative log-likelihood (NLL) components, along joint marginal NLLs, revealed TMB model, correct “inside” ADMB version, incorrect “outside” ADMB version, except marginal NLL. Taken together, analysis confirms bug “outside” ADMB version, affects accuracy Laplace approximation model results.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"estimation-of-additional-observation-error","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Estimation of additional observation error","title":"REMA model equations","text":"Based experience gained using alternative observed index estimates (e.g., relative CPUE indices), appears cases estimates observation error variances biomass /CPUE survey low (e.g., Echave et al. 2020). , mismatch biologically reasonable inter-annual variability precision index estimates. instances, model estimates (σBt,j2+σIt,j2)/σPE2(\\sigma^2_{B_{t,j}}+\\sigma^2_{I_{t,j}})/\\sigma^2_{PE} may lower expected based individual species life history traits. example, ratio observation PE variation low, model predictions population biomass may exhibit high inter-annual variability. behavior unexpected low productivity species, exhibit low inter-annual variation biomass (.e. low PE variance), especially situations fishing exploitation low. One approach addressing issue estimate additional observation error. method commonly implemented Stock Synthesis, implemented Alaskan crab stock assessments (Zheng Siddeek 2020), explored several groundfish assessments well. Using biomass survey variance example, extra estimated error (στ\\sigma_{\\tau}) specified : σln(Bt)=ln((σBtBt)2+στ2+1). \\sigma_{ln(B_t)}=\\sqrt{ln((\\frac{\\sigma_{B_t}}{B_t})^2+{\\sigma_{\\tau}}^2+1)}. approach new method Tier 5 stock assessments apportionment methods AFSC. reproducible example GOA shortraker outlined Fitting additional CPUE survey. example, estimation additional observation error biomass CPUE survey resulted better fit AIC status quo approaches biologically realistic trend predicted population biomass.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"exploration-of-the-tweedie-distribution-for-zero-biomass-observations","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Exploration of the Tweedie distribution for zero biomass observations","title":"REMA model equations","text":"Tweedie distributions family probability distributions can generalized include Gaussian, inverse Gaussian, gamma, Poisson, compound Poisson-gamma distributions. Tweedie distribution can defined three parameters, mean (μ\\mu), power parameter (ρ\\rho), dispersion (ϕ\\phi), relationship variance parameters defined σ2=ϕμρ\\sigma^2=\\phi\\mu^\\rho. ρ\\rho bounded 1 2, Tweedie positive, continuous distribution can handle zero values, thus allowing naturally handle survey time series one zero observations. Values ρ\\rho bounds special cases Tweedie distribution, ρ\\rho=1 equivalent Poisson distribution ρ\\rho=2 gamma distribution. rema, Tweedie ρ\\rho constrained 1 2 using logit-transformation. Similar normal distribution, observation error variances treated known Tweedie distribution REMA model. Using biomass survey observation example, dispersion biomass observation strata jj year tt derived ϕBt,j=σBt,j2(B̂t,j)ρ. \\phi_{B_{t,j}}=\\frac{\\sigma_{B_{t,j}}^2}{{(\\hat{B}_{t,j})}^\\rho}. result one additional parameter (ρ\\rho) estimated applying alternative distribution. measurement equation Tweedie becomes Bt,j=B̂t,j+ϵB,ϵB∼Twρ(0,ϕBt,j). B_{t,j} = \\hat{B}_{t,j} + \\epsilon_{B}, \\text{} \\epsilon_{B} \\sim Tw_{\\rho}(0, \\phi_{B_{t,j}}). σBt,j\\sigma_{B_{t,j}} undefined Bt,jB_{t,j}=0, zeros assumed CV=1.5. assumption can explored user rema, example, user wanted define CV=2.0 zeros, define zeros = list(assumption = 'tweedie', options_tweedie = list(zeros_cv = 2.0)) within prepare_rema_input() function. reproducible example using Tweedie distribution observation error outlined Strategies handling zero biomass observations vignette. use non-shortspine thornyhead component Bering Sea/Aleutian Islands rockfish stock, unique 13 38 bottom trawl survey biomass estimates zeros. Tweedie performs well case, exploration revealed Tweedie models can slow run often converge, especially instances observation error variance estimates low. possible development estimating additional observation error variance (e.g. τ2{\\tau}^2 previous section) alleviate convergence errors. However, interim, Tweedie distribution rema considered experimental, recommend viable alternative tactical assessments time.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"experimental-one-step-ahead-osa-residuals","dir":"Articles","previous_headings":"Base model structure for a single survey and stratum","what":"Experimental One-Step Ahead (OSA) residuals","title":"REMA model equations","text":"use one-step ahead (OSA) residuals, also referred forecast residuals prediction errors, crucial validation state-space models like REMA (Thygesen et al. 2017). Instead comparing observed expected values time step, OSA residuals use forecasted values based previous observations (.e. exclude observed value current time step prediction). Traditional residuals (e.g. Pearson’s residuals) inappropriate state-space models, process error variance may -estimated cases model mispecified, thus leading artificially small residuals (see Section 3 Thygesen et al. 2017 example). Methods calculating OSA residuals implemented TMB TMB::oneStepPredict() function. methods straight forward implement TMB, computationally demanding, validity (accuracy) OSA residuals may vary situation method used. Methods implement OSA residuals rema development considered experimental. Currently, OSA residuals implemented using method = \"cdf\" option TMB::oneStepPredict(), benefit speeding residual calculations saving copies one-step cumulative density function data point thus reducing number calculations function evaluation. OSA residuals appear calculated correctly REMA models, occasionally NaN values residuals returned, especially cases small measurement errors. One potential cause error may lie initialization state process within REMA model explored future versions package.","code":""},{"path":"https://afsc-assessments.github.io/rema/articles/rema_equations.html","id":"references","dir":"Articles","previous_headings":"","what":"References","title":"REMA model equations","text":"Echave, K. B., P. J. F. Hulson, K. . Siwicke. 2020. Assessment Thornyhead stock complex Gulf Alaska. : Stock assessment fishery evaluation report groundfish resources Gulf Alaska projected 2021. North Pacific Fishery Management Council, 605 W. 4th. Avenue, Suite 306, Anchorage, AK 99501. Fournier D. ., H. J. Skaug , J. Ancheta , J. Ianelli , . Magnusson, M. N. Maunder , . Nielsen, J. Sibert. 2012. AD Model Builder: using automatic differentiation statistical inference highly parameterized complex nonlinear models, Optimization Methods Software, 27:2, 33-249, DOI: 10.1080/10556788.2011.597854 Marlow, N. . 1967. normal limit theorem power sums independent normal random variables. Bell System Technical Journal. 46 (9): 2081–2089. doi:10.1002/j.1538-7305.1967.tb04244.x. Methot, R. D. Wetzel, C. R. 2013. Stock Synthesis: biological statistical framework fish stock assessment fishery management. Fisheries Research, 142: 86-99. https://doi.org/10.1016/j.fishres.2012.10.012 Skaug H. D. Fournier. 2013. Random effects AD Model Builder: ADMB-RE User Guide. http://ftp.admb-project.org/admb-11.2pre/admbre-11.2pre.pdf Skaug, H. J. D. . Fournier. 2006. Automatic approximation marginal likelihood non-Gaussian hierarchical models. Computational Statistics & Data Analysis, 51(2): 699–709. https://doi.org/10.1016/j.csda.2006.03.005 Thygesen, U. H., Albertsen, C. M., Berg, C. W., Kristensen, K., Nielsen, . 2017. Validation ecological state space models using Laplace approximation. Environ Ecol Stat 24, 317–339. https://doi.org/10.1007/s10651-017-0372-4 Zheng, J., M. S. M. Siddeek. 2020. Bristol Bay red king crab stock assessment fall 2020. : Stock assessment fishery evaluation report king tanner crab fisheries Bering Sea Aleutian Islands regions. North Pacific Fishery Management Council, 605 W 4th Ave, Suite 306 Anchorage, AK 99501.","code":""},{"path":"https://afsc-assessments.github.io/rema/authors.html","id":null,"dir":"","previous_headings":"","what":"Authors","title":"Authors and Citation","text":"Jane Sullivan. Author, maintainer. Laurinne Balstad. Author, contributor. Cole Monnahan. Contributor. Pete Hulson. Contributor.","code":""},{"path":"https://afsc-assessments.github.io/rema/authors.html","id":"citation","dir":"","previous_headings":"","what":"Citation","title":"Authors and Citation","text":"Sullivan J, Balstad L (2024). rema: generalized framework fit random effects (RE) model, state-space random walk model developed Alaska Fisheries Science Center (AFSC) apportionment biomass estimation groundfish crab stocks.. R package version 1.2.0, https://afsc-assessments.github.io/rema/, https://github.com/JaneSullivan-NOAA/rema.","code":"@Manual{, title = {rema: A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.}, author = {Jane Sullivan and Laurinne Balstad}, year = {2024}, note = {R package version 1.2.0, https://afsc-assessments.github.io/rema/}, url = {https://github.com/JaneSullivan-NOAA/rema}, }"},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"a-generalized-random-effects-model-for-fitting-survey-biomass-estimates-with-the-options-of-including-multiple-survey-strata-and-an-additional-survey-index","dir":"","previous_headings":"","what":"A generalized Random Effects model for fitting survey biomass estimates with the options of including Multiple survey strata and an Additional survey index","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"Sullivan, J., C. Monnahan, P. Hulson, J. Ianelli, J. Thorson, . Havron. 2022a. REMA: consensus version random effects model ABC apportionment Tier 4/5 assessments. Plan Team Report, Joint Groundfish Plan Teams, North Pacific Fishery Management Council. 605 W 4th Ave, Suite 306 Anchorage, AK 99501. Available Oct 2022 Joint GPT e-Agenda. Living documentation REMA model equations available REMA website: https://afsc-assessments.github.io/rema/articles/rema_equations.html","code":""},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"background","dir":"","previous_headings":"","what":"Background","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"random effects (RE) model developed North Pacific Fisheries Management Council (NPFMC) Groundfish Plan Team’s Survey Averaging working group used Alaska Fisheries Science Center (AFSC) since 2013 estimate biomass data-limited groundfish crab stock assessments, apportion catch among management areas. RE model estimates biomass series random effects underlying state dynamics modeled random walk (Oct 2013 Joint GPT minutes). several versions original single-strata (univariate) RE single-survey, multiple-strata (multivariate; REM) models currently use AFSC. RE REM models extended 2017 include addition second relative index abundance (e.g., NMFS longline survey Relative Population Numbers) inform biomass trend estimation additional scaling parameters (REMA; Hulson et al. 2021). Although RE, REM, REMA models share underlying state-space random walk dynamics, Monnahan et al. (2021) found inconsistencies treatment zero biomass observations, use prior penalty process error parameter, weighting observation error likelihood components. purpose rema R library provide unified random effects model flexible enough accommodate multitude use cases AFSC. presented , REMA model generalized include RE REM allows users quickly develop evaluate complex models alternative parameterizations. variable names updated original versions increase interpretability transparency, model rewritten Template Model Builder (TMB; Kristensen et al. 2016). Additionally, rema provides flexible framework following analyses: Compare bridge existing ADMB models currently used Tiers 4 5 stock assessments Tier 3 apportionment REMA; Visualize multiple model results simultaneously conduct model selection appropriate using Akaike Information Criteria (AIC); Analyze alternative approaches handling zero biomass observations, including treating zeros missing values failed surveys, adding small constant zero values manually defining associated coefficient variation values, exploring experimental option model observations using Tweedie distribution, positive, continuous distribution accommodates zeros; Perform model validation using appropriate diagnostic tools latent variable models [e.g., reporting convergence criteria conducting one-step-ahead (OSA) residual analysis]; Evaluate use priors process error scaling parameter parameters, along alternative likelihood penalties. structure, naming conventions, functions, documentation rema inspired modeled Woods Hole Assessment Model (WHAM; Miller Stock 2020), open-source, state-space age-structured assessment model R package.","code":""},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"installation","dir":"","previous_headings":"","what":"Installation","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"Use devtools package install rema package github. devtools installed, must first.","code":"# install.packages(\"devtools\") devtools::install_github(\"afsc-assessments/rema\", dependencies = TRUE, build_vignettes = TRUE) # Example R scripts are downloaded when `rema` is installed. Locate them on your computer by running the following commands: (rema_path <- find.package('rema')) (rema_examples <- file.path(rema_path, 'example_scripts')) list.files(rema_examples) # Vignettes library(rema) browseVignettes(\"rema\") vignette(topic = \"rema_equations\") # view technical details offline vignette(topic = \"ex1_basics\") vignette(topic = \"ex2_cpue\") vignette(topic = \"ex3_zeros\") vignette(topic = \"ex4_model_validation\")"},{"path":"https://afsc-assessments.github.io/rema/index.html","id":"the-rema-worflow","dir":"","previous_headings":"","what":"The rema worflow","title":"A generalized framework to fit the random effects (RE) model, a state-space random walk model developed at the Alaska Fisheries Science Center (AFSC) for apportionment and biomass estimation of groundfish and crab stocks.","text":"Load rema data. user can read biomass abundance index data file (e.g. csv files), can use rwout.rep report file ADMB version RE model using read_admb_re(). Define model structure assumptions using prepare_rema_input(). function allows users quickly transition single two survey model, specify alternative process error structures, add likelihood penalties priors parameters, evaluate alternative assumptions zero biomass observations. Fit specified REMA model using fit_rema() determine whether model met basic convergence criteria. Extract rema model output clean, consistently formatted data frames using tidy_rema(). user can visualize model output using plot_rema(), quickly format tables report. Compare alternative REMA models conduct model selection using compare_rema_models(). Output function includes table Akaike Information Criteria (AIC) appropriate, figures, tidied data frames. function also accepts model output ADMB version RE model easy comparison past models. Taken together, functions allow R users quickly fit interrogate suite simple statistical models TMB without needing software-specific training expertise.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":null,"dir":"Reference","previous_headings":"","what":"Check convergence of REMA model — check_convergence","title":"Check convergence of REMA model — check_convergence","text":"Access quick convergence checks `TMB` `nlminb`. Function modified wham::check_convergence.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Check convergence of REMA model — check_convergence","text":"","code":"check_convergence(mod, ret = FALSE, f = \"\")"},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Check convergence of REMA model — check_convergence","text":"mod output fit_rema ret T/F, return list? Default = FALSE, just prints console","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Check convergence of REMA model — check_convergence","text":"list least first three components: $convergence stats::nlminb, \"0 indicates successful convergence nlminb\" $maxgr Max absolute gradient value, `max(abs(mod$gr(mod$opt$par)))` $maxgr_par Name parameter max gradient $is_sdrep TMB::sdreport performed model, indicates whether performed without error $na_sdrep TMB::sdreport performed without error model, indicates () components diagonal inverted hessian returned NA","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/check_convergence.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Check convergence of REMA model — check_convergence","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":null,"dir":"Reference","previous_headings":"","what":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"check_estimability calculates matrix second-derivatives marginal likelihood w.r.t. fixed effects, see linear combinations estimable (.e. uniquely estimated conditional upon model structure available data, e.g., resulting likelihood ridge singular, non-invertable Hessian matrix)","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"","code":"check_estimability(obj, h)"},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"obj compiled object h optional argument containing pre-computed Hessian matrix","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/check_estimability.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Check for identifiability of fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by fit_tmb. — check_estimability","text":"tagged list hessian message","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"Takes list REMA models fit_rema, returns list ggplot2 objects plotted saved, list tidy_rema data.frames, AIC values.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"","code":"compare_rema_models( rema_models, admb_re = NULL, save = FALSE, filetype = \"png\", path = NULL, xlab = NULL, biomass_ylab = \"Biomass\", cpue_ylab = \"CPUE\" )"},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"rema_models list REMA models compared. REMA model list list object output fit_rema admb_re list ADMB RE model input/output read_admb_re. Accepts single list, list multiple ADMB RE models. admb_re provided, AIC calculations conducted. save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. filetype (optional) character string; type figure file. Default = 'png'. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. xlab (optional) label x-axis biomass CPUE plots (e.g. 'Year'). Default = NULL. biomass_ylab (optional) label y-axis biomass plots (e.g. 'Biomass (t)'). Default = 'Biomass'. cpue_ylab (optional) label y-axis CPUE plots (e.g. 'Relative Population Number'). Default = 'CPUE'.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"list following items: $output list tidied dataframes include parameter estimates, biomass optional CPUE data, REMA model predictions model compared. Results given variable included applicable comparison models. example, CPUE fit one model another, compare$output$cpue_by_strata) return informational message instead dataframe. See tidy_rema information. $plots ggplot2 figure objects compare$output data. $aic dataframe Akaike Information Criteria (AIC) values. output underlying models fit data.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/compare_rema_models.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot REMA model comparisons and return AIC values when appropriate — compare_rema_models","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":null,"dir":"Reference","previous_headings":"","what":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"extract_fixed extracts best previous value fixed effects, way works mixed fixed effect models","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"","code":"extract_fixed(obj)"},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"obj, compiled object","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/extract_fixed.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Extract fixed effects Originally provided by https://github.com/kaskr/TMB_contrib_R/TMBhelper Internal function called by check_estimability. — extract_fixed","text":"vector fixed-effect estimates","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":null,"dir":"Reference","previous_headings":"","what":"Fit REMA model — fit_rema","title":"Fit REMA model — fit_rema","text":"Fits compiled REMA model using TMB::MakeADFun stats::nlminb. Source code documentation modified wham::fit_wham.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Fit REMA model — fit_rema","text":"","code":"fit_rema( input, n.newton = 0, do.sdrep = TRUE, model = NULL, do.check = FALSE, MakeADFun.silent = TRUE, do.fit = TRUE, save.sdrep = TRUE )"},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Fit REMA model — fit_rema","text":"input Named list output prepare_rema_input, includes following components needed fit model using TMB::MakeADFun: $data Data, list data objects model fitting specification (e.g., user-defined pentalties, index pointers, etc.). required input MakeADFun. $par Parameters, list random fixed effects parameter objects. required input MakeADFun. $map Map, mechanism collecting fixing parameters TMB. input MakeADFun. $random Character vector defining parameters treat random effects. input MakeADFun. $model_name Character, name model, e.g. \"GOA shortraker LLS depth strata\". Useful model comparison. n.newton integer, number additional Newton steps optimization. option currently needed, passed fit_tmb. Default = 0. .sdrep T/F, calculate standard deviations model parameters? See sdreport. Default = TRUE. model (optional), previously fit rema model. .check T/F, check model parameters identifiable? Passed fit_tmb. Runs internal function check_estimability, originally provided https://github.com/kaskr/TMB_contrib_R/TMBhelper. Default = TRUE. MakeADFun.silent T/F, Passed silent argument TMB::MakeADFun. Default = TRUE. .fit T/F, fit model using fit_tmb. Default = TRUE. save.sdrep T/F, save full TMB::sdreport object? FALSE, save summary.sdreport reduce model object file size. Default = TRUE.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Fit REMA model — fit_rema","text":"fit TMB model additional output specified: $rep List derived quantity estimates (e.g. estimated biomass) $sdrep Parameter estimates (standard errors .sdrep = TRUE)","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Fit REMA model — fit_rema","text":"Future development: Implement one-step-ahead (OSA) residuals evaluating model goodness--fit TMB::oneStepPredict). OSA residuals appropriate standard residuals models random effects (Thygeson et al. (2017). See wham example OSA implementation additional OSA residual options (e.g. full Gaussian approximation instead (default) generic method using osa.opts=list(method=\"fullGaussian\").","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/fit_rema.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Fit REMA model — fit_rema","text":"","code":"if (FALSE) { # \\dontrun{ # place holder for example code } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":null,"dir":"Reference","previous_headings":"","what":"Fit TMB model using nlminb — fit_tmb","title":"Fit TMB model using nlminb — fit_tmb","text":"Runs optimization TMB model using stats::nlminb. specified, takes additional Newton steps calculates standard deviations. Internal function called fit_rema. Source code documentation modified wham::fit_tmb.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Fit TMB model using nlminb — fit_tmb","text":"","code":"fit_tmb( model, n.newton = 0, do.sdrep = TRUE, do.check = FALSE, save.sdrep = FALSE )"},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Fit TMB model using nlminb — fit_tmb","text":"model Output TMB::MakeADFun. n.newton Integer, number additional Newton steps optimization. Default = 0. .sdrep T/F, calculate standard deviations model parameters? See TMB::sdreport. Default = TRUE. .check T/F, check model parameters identifiable? Runs internal check_estimability, originally provided https://github.com/kaskr/TMB_contrib_R/TMBhelper. Default = TRUE. save.sdrep T/F, save full TMB::sdreport object? FALSE, save summary.sdreport) reduce model object file size. Default = FALSE.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/fit_tmb.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Fit TMB model using nlminb — fit_tmb","text":"model, appends following: model$opt Output stats::nlminb model$date System date model$dir Current working directory model$rep model$report() model$TMB_version Version TMB installed model$parList List parameters, model$env$parList() model$final_gradient Final gradient, model$gr() model$sdrep Estimated standard deviations model parameters, TMB::sdreport summary.sdreport)","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":null,"dir":"Reference","previous_headings":"","what":"Get one-step-head (OSA) — get_osa_residuals","title":"Get one-step-head (OSA) — get_osa_residuals","text":"Takes rema model output fit_rema returns OSA residuals calculated using TMB::oneStepPredict accompanying residual analysis plots. IMPORTANT: OSA residuals work users implementing Tweedie distribution.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Get one-step-head (OSA) — get_osa_residuals","text":"","code":"get_osa_residuals( rema_model, options = list(method = \"fullGaussian\", parallel = TRUE) )"},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Get one-step-head (OSA) — get_osa_residuals","text":"rema_model list output fit_rema, includes model results also inputs. note OSA residual calculations rema_model$input$osa object, data.frame containing data observations fit model residuals associated . options list options calculating OSA residuals, passed TMB::oneStepPredict. Default: options = list(method = \"fullGaussian\", parallel = TRUE). Alternative methods include \"cdf\", \"oneStepGeneric\", \"oneStepGaussianOffMode\", \"oneStepGaussian\".","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Get one-step-head (OSA) — get_osa_residuals","text":"list tidied data.frames containing biomass CPUE survey residuals accompanying data, well QQ-plot, histogram residuals, plots residuals~year residuals~fitted values strata biomass CPUE survey.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/get_osa_residuals.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Get one-step-head (OSA) — get_osa_residuals","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/pipe.html","id":null,"dir":"Reference","previous_headings":"","what":"Pipe function — %>%","title":"Pipe function — %>%","text":"Allows use pipe function, %>%","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"Takes list output tidy_rema returns list ggplot2 objects plotted saved.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"","code":"plot_extra_cv( tidy_rema, save = FALSE, filetype = \"png\", path = NULL, xlab = NULL, biomass_ylab = \"Biomass\", cpue_ylab = \"CPUE\" )"},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"tidy_rema list output tidy_extra_cv, includes inputs, model results, confidence intervals total observation error (fixed + estimated) save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. filetype (optional) character string; type figure file. Default = 'png'. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. xlab (optional) label x-axis biomass CPUE plots (e.g. 'Year'). Default = NULL. biomass_ylab (optional) label y-axis biomass plots (e.g. 'Biomass (t)'). Default = 'Biomass'. cpue_ylab (optional) label y-axis CPUE plots (e.g. 'Relative Population Number'). Default = 'CPUE'.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"list ggplot2 plots character string messages data. Except parameter estimates, objects output tidy_rema outputted function.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/plot_extra_cv.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot the additional estimated observation error for biomass by strata and/or cpue by strata — plot_extra_cv","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot survey data and model output — plot_rema","title":"Plot survey data and model output — plot_rema","text":"Takes list output tidy_rema returns list ggplot2 objects plotted saved.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot survey data and model output — plot_rema","text":"","code":"plot_rema( tidy_rema, save = FALSE, filetype = \"png\", path = NULL, xlab = NULL, biomass_ylab = \"Biomass\", cpue_ylab = \"CPUE\" )"},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot survey data and model output — plot_rema","text":"tidy_rema list output tidy_rema, includes model results also inputs save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. filetype (optional) character string; type figure file. Default = 'png'. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. xlab (optional) label x-axis biomass CPUE plots (e.g. 'Year'). Default = NULL. biomass_ylab (optional) label y-axis biomass plots (e.g. 'Biomass (t)'). Default = 'Biomass'. cpue_ylab (optional) label y-axis CPUE plots (e.g. 'Relative Population Number'). Default = 'CPUE'.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot survey data and model output — plot_rema","text":"list ggplot2 plots character string messages data. Except parameter estimates, objects output tidy_rema outputted function.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/plot_rema.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot survey data and model output — plot_rema","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":null,"dir":"Reference","previous_headings":"","what":"Prepare input data and parameters for REMA model — prepare_rema_input","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"data read R (either manually .csv data file using read_admb_re), function prepares data parameter settings fit_rema. model can set run single survey mode one strata, multi-survey mode, uses additional relative abundance index (.e. cpue) inform predicted biomass. optional inputs described related CPUE survey data scaling parameter q, cpue_dat options_q used multi_survey = 1. function structure documentation modeled wham::prepare_wham_input.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"","code":"prepare_rema_input( model_name = \"REMA for unnamed stock\", multi_survey = 0, admb_re = NULL, biomass_dat = NULL, cpue_dat = NULL, sum_cpue_index = FALSE, start_year = NULL, end_year = NULL, wt_biomass = NULL, wt_cpue = NULL, PE_options = NULL, q_options = NULL, zeros = NULL, extra_biomass_cv = NULL, extra_cpue_cv = NULL )"},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"model_name name stock identifier REMA model multi_survey switch run model single multi-survey mode. 0 (default) = single survey, 1 = multi-survey. admb_re list object returned read_admb_re.R, includes biomass survey data (admb_re$biomass_dat), optional cpue survey data (admb_re$cpue_dat), years model predictions (admb_re$model_yrs), model predictions log biomass strata correct format input REMA (admb_re$init_log_biomass_pred). supplied, user need enter biomass_dat cpue_dat. biomass_dat data.frame biomass survey data long format following columns: strata character; survey name, survey region, management unit, depth strata. Note user must include column even one survey strata year integer; survey year. Note user needs include years observations (.e. need supply NULL NA values missing survey years) biomass numeric; biomass estimate/observation (e.g. bottom trawl survey biomass mt). default, biomass == 0, value treated NA (.e., failed survey). user wants make assumptions zeros (e.g. adding small constant), must define data manually. cv numeric; coefficient variation (CV) biomass estimate (.e. sd(biomass)/biomass) cpue_dat (optional) data.frame relative abundance index (.e. cpue) data long format following columns: strata character; survey name, survey region, management unit, depth strata (note user must include column even one survey strata) year integer; survey year. Note user needs include years observations (.e. need supply NULL NA values missing survey years) cpue numeric; cpue estimate/observation (e.g. longline survey cpue relative population number); default, cpue == 0, value treated NA (.e., failed survey). user wants make assumptions zeros (e.g. adding small constant), must define data manually. cv numeric; coefficient variation (CV) cpue estimate (.e. sd(cpue)/cpue) sum_cpue_index T/F 1/0, CPUE survey index able summed across strata get total CPUE survey index? example, Longline survey relative population numbers (RPNs) summable longline survey numbers per hachi (CPUE) . Default = FALSE. start_year (optional) integer value specifying start year estimation model; admb_re supplied, value defaults start_year = min(admb_re$model_yrs); admb_re supplied, value defaults first year either biomass_dat cpue_dat end_year (optional) integer value specifying last year estimation model; admb_re supplied, value defaults end_year = max(admb_re$model_yrs); admb_re supplied, value defaults last year either biomass_dat cpue_dat wt_biomass (optional) multiplier biomass survey data component negative log likelihood. example, nll = wt_biomass * nll. Defaults wt_biomass = 1 wt_cpue (optional) multiplier CPUE survey data component negative log likelihood. example, nll = wt_cpue * nll. Defaults wt_cpue = 1 PE_options (optional) customize implementation process error (PE) parameters, including options share PE across biomass survey strata, change starting values, fix parameters, add penalties priors (see details) q_options (optional) customize implementation scaling parameters (q), including options define q biomass cpue survey cpue strata, change starting values, fix parameters, add penalties priors (see details). used multi_survey = 1 zeros (optional) define assumptions treat zero biomass CPUE observations, including treating zeros NAs, changing zeros small constants fixed CVs, modeling zeros using Tweedie distribution (see details). extra_biomass_cv (optional) estimate additional observation error biomass survey data (see details). default, assumption = \"extra_cv\" estimate one extra CV parameter, regardless number biomass survey strata. extra_cpue_cv (optional) estimate additional observation error CPUE survey data (see details). default, assumption = \"extra_cv\" estimate one extra CV parameter, regardless number CPUE survey strata.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"function returns named list following components: data Named list data, passed TMB::MakeADFun par Named list parameters, passed TMB::MakeADFun map Named list defining optionally collect fix parameters, passed TMB::MakeADFun random Character vector parameters treat random effects, passed TMB::MakeADFun model_name Name stock identifier REMA model biomass_dat tidied long format data.frame biomass survey observations associated CVs strata. data.frame 'complete' include modeled years, missing values treated NAs. Note data.frame differ admb_re$biomass_dat input biomass assumptions zero biomass observations different ADMB model user specifies REMA. user can change assumptions zeros using zeros argument. cpue_dat optional CPUE survey data provided multi_survey = 1, tidied long-format data.frame CPUE survey observations associated CVs strata. data.frame 'complete' include modeled years, missing values treated NAs. Note data.frame differ admb_re$biomass_dat input biomass assumptions zero CPUE observations different ADMB model user specifies REMA. user can change assumptions zeros using zeros argument. optional CPUE survey data provided multi_survey = 0, object NULL.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"PE_options allows user specify options process error (PE) parameters. NULL, default PE specifications used: one PE parameter estimated biomass survey strata, initial values log_PE set 1, penalties priors added. user can modify default PE_options using following list entries: $pointer_PE_biomass index customize assignment PE parameters individual biomass strata. Vector length = number biomass strata, starting index 1 ending number unique PE estimated. example, three biomass survey strata user wants estimate one PE, specify pointer_PE_biomass = c(1, 1, 1). default one unique log_PE estimated unique biomass survey stratum $initial_pars vector initial values log_PE. default initial value log_PE 1. $fix_pars Option fix PE parameters, user specifies index value PE parameter like fix initial value. example, three biomass survey strata, user wants fix log_PE second stratum estimate log_PE first third strata specify fix_pars = c(2) Note option recommended. $penalty_options Warning: following options experimental well-tested. Options penalizing PE likelihood adding prior log_PE include following: \"none\" (default) penalty prior used \"wt\" multiplier PE random effects component negative log likelihood. example, nll = wt * nll, wt = 1.5 specified single value penalty_values argument \"squared_penalty\" implemented earlier version RE.tpl, penalty prevents PE shrinking zero. example, nll = nll + (log_PE + squared_penalty)^2, squared_penalty = 1.5. vector squared_penalty values specified PE penalty_values argument \"normal_prior\" Normal prior log space, nll = nll - dnorm(log_PE, pmu_log_PE, psig_log_PE, 1) pmu_log_PE psig_log_PE specified PE parameter penalty_values argument penalty_values user-defined values penalty_options. penalty type entered follows: \"none\" (default) NULL example, penalty_values = NULL \"wt\" single numeric value. example, penalty_values = 1.5 \"squared_penalty\" vector numeric values length = number PE parameters. example, three PE parameters estimated user wants penalty one, use penalty_values = c(1.5, 1.5, 1.5) \"normal_prior\" vector paired values PE parameter, vector pair prior mean log_PE pmu_log_PE associated standard deviation psig_log_PE. example, three PE parameters estimated user wants normal prior log_PE ~ N(1.0, 0.08), penalty_values = c(c(1.0, 0.08), c(1.0, 0.08), c(1.0, 0.08)) q_options allows user specify options CPUE survey scaling parameters (q). multi_survey = 0 (default), q parameters estimated regardless user defines q_options. multi_survey = 0 q_options = NULL, default q specifications used: one q parameter estimated CPUE survey strata, biomass CPUE surveys assumed share strata definitions (.e., biomass_dat cpue_dat number columns columns represent strata), initial values log_q set 1, penalties priors added. user can modify default q_options using following list entries: $pointer_q_cpue index customize assignment q parameters individual CPUE survey strata. Vector length = number CPUE strata, starting index 1 ending number unique q parameters estimated. example, three CPUE survey strata user wanted estimate one q, specify pointer_q_cpue = c(1, 1, 1). recommended model configuration estimate one log_q CPUE survey stratum. $pointer_biomass_cpue_strata index customize assignment biomass predictions individual CPUE survey strata. Vector length = number biomass survey strata, starting index 1 ending number unique CPUE survey strata. pointer needs defined number biomass CPUE strata equal. pointer_biomass_cpue_strata option allows user calculate predicted biomass CPUE survey strata level scenario biomass survey strata higher resolution CPUE survey strata. example, 3 biomass survey strata represented 2 CPUE survey strata, user may specify pointer_biomass_cpue_strata = c(1, 1, 2). specification assign first 2 biomass strata first CPUE strata, third biomass stratum second CPUE stratum. CPUE data compliment specific biomass stratum, user can populate NAs. example pointer_biomass_cpue_strata = c(1, NA, 3), means CPUE data biomass strata 1 3 2. NOTE: scenario CPUE survey strata biomass survey strata CPUE survey used inform biomass survey trend. error thrown q_options$pointer_biomass_cpue_strata defined biomass CPUE survey strata definitions . $initial_pars vector initial values log_q. default initial value log_q 1. $fix_pars Option fix q parameters, user specifies index value q parameter like fix initial value. example, three CPUE survey strata, user wants fix log_q second stratum estimate log_q first third strata specify fix_pars = c(2) $penalty_options Options penalizing q likelihood adding prior log_q include following: \"none\" (default) penalty prior used \"normal_prior\" Warning, experimental well-tested. Normal prior log space, nll = nll - dnorm(log_q, pmu_log_q, psig_log_q, 1) pmu_log_q psig_log_q specified q parameter penalty_values argument penalty_values user-defined values penalty_options. penalty type entered follows: \"none\" (default) NULL example, penalty_values = NULL \"normal_prior\" vector paired values q parameter, vector pair prior mean log_q pmu_log_q associated standard deviation psig_log_q. example, 2 q parameters estimated user wants normal prior log_q ~ N(1.0, 0.05), penalty_values = c(c(1.0, 0.05), c(1.0, 0.05)) zeros allows user specify options treat zero biomass CPUE survey observations. default zero observations treated NAs warning msg effect returned console. zeros allows user specify non-default zero assumptions using following list entries: $assumption character, name assumption using. three alternatives currently implemented, zeros = list(assumption = c(\"NA\", \"small_constant\", \"tweedie\"). \"NA\" default; option assumes zero estimates failed surveys removes . \"small_constant\" ad hoc method fixed value added zero assumed CV. default, small constant = 0.0001 CV value entered user data. user can change assumed value CV using options_small_constant. \"tweedie\" uses Tweedie assumed error distribution survey data, allows zeros. alternative estimates one additional power parameter. assumed CV zero biomass zero cpue survey observations defaults 1.5. user can change assumed CV, change initial values inverse logit transformed power parameter, fix initial values using options_tweedie. $options_small_constant vector length two numeric values. first value small constant add zero observation, second user-defined coefficient value. user can specify small value use input CV specifying NA second value. E.g., 'options_small_constant = c(0.0001, NA)'. $options_tweedie list entries control initial fixed values Tweedie parameters. Currently, argument accepts following entries: extra_biomass_cv allows user specify options estimating additional CV parameter (log_tau_biomass source code, estimated log-space) biomass survey observations. extra_biomass_cv = NULL (default), extra CV estimated. user can modify default extra_biomass_cv options using following list entries: $assumption string identifying assumption used biomass survey observations. Options include \"none\" (default extra CV estimated) \"extra_cv\". assumption = \"extra_cv\", default one extra CV estimated, regardless many biomass strata defined. extra_biomass_cv NULL, user must define appropriate assumption. $pointer_extra_biomass_cv index customize assignment extra CV parameters individual biomass survey strata. Vector length = number biomass strata, starting index 1 ending number unique extra CV parameters estimated. three biomass survey strata user wanted estimate extra CV per stratum, specify pointer_extra_biomass_cv = c(1, 2, 3). default, one additional parameter estimated, regardless many strata defined (.e. pointer_extra_biomass_cv = c(1, 1, 1)). $initial_pars vector initial values extra biomass log_tau_biomass. default initial value log_tau_biomass log(1e-7) (approximately 0 arithmetic scale). $fix_pars Option fix extra biomass CV parameters, user specifies index value parameter like fix initial value. example, three biomass survey strata defined pointer_extra_biomass_cv, user wants fix log_tau_biomass second stratum estimate log_tau_biomass first third strata specify fix_pars = c(2). extra_cpue_cv allows user specify options estimating additional CV parameter (log_tau_cpue source code, estimated log-space) cpue survey observations. extra_cpue_cv = NULL (default), extra CV estimated. user can modify default extra_cpue_cv options using following list entries: $assumption string identifying assumption used cpue survey observations. Options include \"none\" (default extra CV estimated) \"extra_cv\". assumption = \"extra_cv\", default one extra CV estimated, regardless many cpue strata defined. extra_cpue_cv NULL, user must define appropriate assumption. $pointer_extra_cpue_cv index customize assignment extra CV parameters individual cpue survey strata. Vector length = number cpue strata, starting index 1 ending number unique extra CV parameters estimated. three cpue survey strata user wanted estimate extra CV per stratum, specify pointer_extra_cpue_cv = c(1, 2, 3). default, one additional parameter estimated, regardless many strata defined (.e. pointer_extra_cpue_cv = c(1, 1, 1)). $initial_pars vector initial values extra cpue log_tau_cpue. default initial value log_tau_cpue log(1e-7) (approximately 0 arithmetic scale). $fix_pars Option fix extra cpue CV parameters, user specifies index value parameter like fix initial value. example, three cpue survey strata defined pointer_extra_cpue_cv, user wants fix log_tau_cpue second stratum estimate log_tau_cpue first third strata specify fix_pars = c(2).","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/prepare_rema_input.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Prepare input data and parameters for REMA model — prepare_rema_input","text":"","code":"if (FALSE) { # \\dontrun{ # place holder for example code } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":null,"dir":"Reference","previous_headings":"","what":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"Read report file ADMB version RE model (rwout.rep) convert long format survey data estimates CVs input REMA.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"","code":"read_admb_re( filename, model_name = \"Unnamed ADMB RE model\", biomass_strata_names = NULL, cpue_strata_names = NULL )"},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"filename name ADMB output file read (e.g. rwout.rep) model_name (optional) Name stock identifier ADMB version RE model. Defaults 'ADMB RE' biomass_strata_names (optional) vector character names corresponding names biomass survey strata. Vector order columns srv_est rwout.rep cpue_strata_names (optional) vector character names corresponding names CPUE survey strata. Vector order columns srv_est_LL rwout.rep version ADMB RE model accepts additional survey index","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"object type \"list\" biomass optional cpue survey data long format, initial parameter values log_biomass_pred (random effects matrix), ready input REMA list following items: $biomass_dat dataframe biomass survey data strata, year, biomass estimates, CVs. Note CVs back-transformed natural space. $cpue_dat Optional dataframe CPUE survey data strata, year, CPUE estimates, CVs. Note CVs back-transformed natural space. $model_yrs Vector prediction years. $init_log_biomass_pred Matrix initial parameter values log_biomass_pred (random effects matrix), ready input REMA. $admb_re_results list ADMB RE model results ready comparison REMA models using compare_rema_models(). User beware... many, many versions RE.tpl existence individual variances may cause errors output.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_admb_re.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Convert ADMB version of the RE model data and output to REMA inputs — read_admb_re","text":"","code":"if (FALSE) { # \\dontrun{ # place holder for example code } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":null,"dir":"Reference","previous_headings":"","what":"Read ADMB .rep file and return an R object of type 'list' — read_rep","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"Code modified original function provided Steve Martell, D'Arcy N. Webber called read_admb_re","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"","code":"read_rep(fn)"},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"fn full path name ADMB output file read","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"object type \"list\" ADMB outputs therein","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/read_rep.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Read ADMB .rep file and return an R object of type 'list' — read_rep","text":"","code":"if (FALSE) { # \\dontrun{ read_rep(fn = 'inst/example_data/goasr.rep') } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":null,"dir":"Reference","previous_headings":"","what":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"Takes list output tidy_rema returns list enhanced versions biomass_by_strata cpue_by_strata appropriate. enhanced dataframes include three new columns, tot_log_obs_cv, tot_obs_lci, tot_obs_uci, represent combined log-space standard error associated confidence intervals include assumed estimated additional observation error.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"","code":"tidy_extra_cv(tidy_rema, save = FALSE, path = NULL, alpha_ci = 0.05)"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"tidy_rema list output tidy_rema, includes model results also inputs save (optional) logical (T/F) save figures filetype path. Default = FALSE. YET IMPLEMENTED. path (optional) directory path location figure files saved save = TRUE. YET IMPLEMENTED. alpha_ci (optional) significance level generating confidence intervals. Default = 0.05","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"list following items: $parameter_estimates data.frame fixed effects parameters REMA (e.g. log_PE log_q) standard errors confidence intervals transformed log space natural space ease interpretation. $biomass_by_strata tidy, long format data.frame model predicted observed biomass biomass survey strata. data.frame now enhanced new columns include log-space standard error associated confidence intervals account additional estimated observation error. $cpue_by_strata tidy, long format data.frame model predicted observed CPUE CPUE survey strata. data.frame now enhanced new columns include log-space standard error associated confidence intervals account additional estimated observation error. REMA run multi-survey mode, CPUE data provided, explanatory character string instructions fitting CPUE data returned. $biomass_by_cpue_strata tidy, long format data.frame model predicted biomass CPUE survey strata. Note observed/summed biomass observations returned case missing values one stratum another within given year. output reserved instances number biomass strata exceeds CPUE survey strata, user wants visualize predicted biomass resolution CPUE predictions. scenarios, character string returned explaining special use case object. $total_predicted_biomass tidy, long format data.frame total model predicted biomass summed across biomass survey strata. one stratum used (.e. univariate RE), predicted values output$biomass_by_strata. $total_predicted_cpue tidy, long format data.frame total model predicted CPUE summed across CPUE survey strata. one stratum used (.e. univariate RE), predicted values output$cpue_by_strata. CPUE survey index provided defined summable prepare_rema_input(), character string returned explaining change using 'sum_cpue_index' ?prepare_rema_input appropriate.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_extra_cv.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Tidy estimates of extra biomass or CPUE index CV — tidy_extra_cv","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":null,"dir":"Reference","previous_headings":"","what":"Tidy REMA model output — tidy_rema","title":"Tidy REMA model output — tidy_rema","text":"Takes outputs fit_rema, returns named list tidied data.frames include parameter estimates standard errors, derived variables model. information \"tidy\" data, please see Wickham 2014. code modified wham::par_tables_fun.","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Tidy REMA model output — tidy_rema","text":"","code":"tidy_rema(rema_model, save = FALSE, path = NULL, alpha_ci = 0.05)"},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Tidy REMA model output — tidy_rema","text":"rema_model list output fit_rema, includes model results also inputs save (optional) logical (T/F) save output data.frames csvs path. Default = FALSE. YET IMPLEMENTED. path (optional) directory path location csvs saved save = TRUE. YET IMPLEMENTED. alpha_ci (optional) significance level generating confidence intervals. Default = 0.05","code":""},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Tidy REMA model output — tidy_rema","text":"list following items: $parameter_estimates data.frame fixed effects parameters REMA (e.g. log_PE log_q) standard errors confidence intervals transformed log space natural space ease interpretation. $biomass_by_strata tidy, long format data.frame model predicted observed biomass biomass survey strata. $cpue_by_strata tidy, long format data.frame model predicted observed CPUE CPUE survey strata. REMA run multi-survey mode, CPUE data provided, explanatory character string instructions fitting CPUE data returned. $biomass_by_cpue_strata tidy, long format data.frame model predicted biomass CPUE survey strata. Note observed/summed biomass observations returned case missing values one stratum another within given year. output reserved instances number biomass strata exceeds CPUE survey strata, user wants visualize predicted biomass resolution CPUE predictions. scenarios, character string returned explaining special use case object. $total_predicted_biomass tidy, long format data.frame total model predicted biomass summed across biomass survey strata. one stratum used (.e. univariate RE), predicted values output$biomass_by_strata. $total_predicted_cpue tidy, long format data.frame total model predicted CPUE summed across CPUE survey strata. one stratum used (.e. univariate RE), predicted values output$cpue_by_strata. CPUE survey index provided defined summable prepare_rema_input(), character string returned explaining change using 'sum_cpue_index' ?prepare_rema_input appropriate.","code":""},{"path":[]},{"path":"https://afsc-assessments.github.io/rema/reference/tidy_rema.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Tidy REMA model output — tidy_rema","text":"","code":"if (FALSE) { # \\dontrun{ # placeholder for example } # }"}]