--- title: "Defining Priors" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{k_defining_priors} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` In some cases, prior information might be available for use to inform model estimation. Priors can either serve as informative priors that have a strong influence on model estimation (e.g., natural mortality, catchability), but can also serve as regularizing priors that help stabilize model estimation. These regularizing priors can be helpful particularly for stabalizing tag reporting, movement, and selectivity parameters because some parameters are not well informed by data (e.g., slope of selectivity when it approaches infinity) such that the gradient is zero, with respect to the likelihood surface. In `SPoRC`, priors can be placed on several parameters. These include: 1. Steepness, 2. Natural Mortality, 3. Catchability, 4. Selectivity, 5. Movement, 6. Tag Reporting Rates, and 7. Recruitment Proportions In the following, we will demonstrate how priors can be specified for each of these parameters. For further mathematical details on the formulation of these priors, refer to the vignette describing model equations. Let us first load in any required packages, data we may use, and define model dimensions. ```{r setup} # Load in packages library(SPoRC) data("sgl_rg_sable_data") # load in data # Define model dimensions input_list <- Setup_Mod_Dim(years = 1:length(sgl_rg_sable_data$years), # vector of years # (corresponds to year 1960 - 2024) ages = 1:length(sgl_rg_sable_data$ages), # vector of ages lens = seq(41,99,2), # number of lengths n_regions = 1, # number of regions n_sexes = sgl_rg_sable_data$n_sexes, # number of sexes == 1, # female, == 2 male n_fish_fleets = sgl_rg_sable_data$n_fish_fleets, # number of fishery # fleet == 1, fixed gear, == 2 trawl gear n_srv_fleets = sgl_rg_sable_data$n_srv_fleets, # number of survey fleets verbose = TRUE ) ``` # Steepness To specify priors on the steepness of a Beverton-Holt stock-recruitment relationship, we use the `Setup_Mod_Rec` function with default settings for simplicity. First, we define a dataframe that includes the region where the prior applies, the mean (`mu`), and standard deviation (`sd`) of the steepness prior. These values are specified in normal space but are internally transformed to a beta distribution bounded between 0.2 and 1. To activate steepness priors, set `Use_h_prior = 1` and provide the dataframe to the `h_prior` argument: ```{r} # set up dataframe specifying the steepness prior and the dimensions it is applied to steepness_prior <- data.frame( region = 1, # apply steepness prior to region 1 mu = 0.6, # mean of steepness sd = 0.2 # sd of steepness ) # setup recruitment options input_list <- Setup_Mod_Rec( input_list = input_list, # input data list from above # Model options rec_model = "bh_rec", # recruitment model type Use_h_prior = 1, # whether steepness priors are used h_prior = steepness_prior # Steepness prior dataframe ) ``` In a spatial model where region-specific steepness priors are required, the prior dataframe can include multiple entries. In particular, each row corresponds to a region-specific steepness parameter, along with its associated mean (`mu`) and standard deviation (`sd`). For instance: ```{r} steepness_prior <- data.frame( region = c(1,2,3), # apply steepness prior to region 1, 2 and, 3 mu = c(0.6, 0.8, 0.3), # mean of steepness prior corresponding to the regions defined sd = c(0.2, 0.2, 0.2) # sd of steepness corresponding to the regions defined ) steepness_prior ``` # Natural Mortality To define natural mortality priors, we use the `Setup_Mod_Biologicals` function to specify a natural mortality prior. To do so, we will set `Use_M_prior = 1` and provide a dataframe of values corresponding the region, year, age, and sex blocks specified for natural mortality, along with their associated mean and standard deviation. ```{r} # Define natural mortality prior M_prior <- data.frame( regionblk = 1, # region block yearblk = 1, # year block ageblk = c(1), # age block sexblk = c(1,2), # sex block mu = 0.085, # mean sd = 0.1 #sd ) input_list <- Setup_Mod_Biologicals(input_list = input_list, # Data inputs WAA = sgl_rg_sable_data$WAA, # weight-at-age MatAA = sgl_rg_sable_data$MatAA, # maturity at age AgeingError = as.matrix(sgl_rg_sable_data$age_error), # ageing error SizeAgeTrans = sgl_rg_sable_data$SizeAgeTrans, # size age transition matrix # Model options Use_M_prior = 1, # use natural mortality prior M_prior = M_prior, # mean and sd for M prior M_spec = "est_ln_M", # estmating M M_sexblk_spec = list(1,2) # specifying M as sex-specific ) ``` # Catchability Defining catchability priors for the fishery and survey follows the same process. For demonstration purposes, the following section will illustrate how catchability priors can be defined for the survey using the `Setup_Mod_Srvsel_and_Q` function. Similar to the steepness section, a dataframe of catchability priors is first defined that determines the region, fleet, and block the catchability prior is applied to. The dataframe also includes columns defining the the mean (`mu`), and standard deviation (`sd`) of the catchability prior. Again, these priors are provided in normal space, but are internally transformed to be a lognormal density. The `Use_srv_q_prior = 1` argument is utilized to activate the prior, and the dataframe is provided to the `srv_q_prior` argument. ```{r} # set up dataframe specifying the survey catchability prior and the dimensions it is applied to srv_q_prior <- data.frame( region = 1, # apply catchability prior to region 1 fleet = c(1, 2, 3), # apply catchability prior to region 1, fleets 1 - 3 block = 1, # apply catchability prior to region 1, fleets 1 - 3, block 1 mu = c(6, 0.5, 4), # corresponding mean catchability prior values sd = c(0.1, 0.35, 0.1) # corresponding standard deviation prior values ) srv_q_prior # Setup survey selectivity and catchability input_list <- Setup_Mod_Srvsel_and_Q(input_list = input_list, # Model options # survey selectivity, whether continuous time-varying cont_tv_srv_sel = c("none_Fleet_1", "none_Fleet_2", "none_Fleet_3"), # survey selectivity blocks srv_sel_blocks = c("none_Fleet_1", "none_Fleet_2", "none_Fleet_3"), # survey selectivity form srv_sel_model = c("logist1_Fleet_1", "exponential_Fleet_2", "logist1_Fleet_3"), # survey catchability blocks srv_q_blocks = c("none_Fleet_1", "none_Fleet_2", "none_Fleet_3"), # whether to estiamte all fixed effects for survey selectivity srv_fixed_sel_pars_spec = c("est_all", "est_all", "est_all"), # whether to estiamte all fixed effects for survey catchability srv_q_spec = c("est_all", "est_all", "est_all"), # Using catchability priors Use_srv_q_prior = 1, # Input survey catchability prior dataframe srv_q_prior = srv_q_prior ) ``` # Selectivity Defining selectivity priors for the fishery and survey follows the same process. Here, we will focus on defining selectivity priors for the survey for brevity, using the `Setup_Mod_Srvsel_and_Q` function. We first define a dataframe with the following columns: `region`, `par`, `block`, `sex`, `fleet`, `mu`, `sd`, which defines the region, parameter, time block, sex, and fleet combination to apply the prior to. Given the complexity of the number of selectivity parameters that can be estimated across model partitions (regions, sexes, blocks, and fleets), users will need to be cautious how the selectivity prior dataframe is defined to ensure they are correctly applied. In the following, we will apply priors to the logistic selectivity parameters for both sexes in fleet 1, and priors to the power parameter describing exponential selectivity in fleet 2. ```{r} srv_selex_prior <- data.frame( region = 1, # only a single region par = c(1,2,1,2,1,1), # parameter to apply prior to sex = c(1,1,2,2,1,2), # sexes to apply prior to fleet = c(1,1,1,1,2,2), # fleets to apply prior to block = 1, # only a single time block mu = rep(0, 6), # mean specified at 0 sd = rep(5, 6) # sd specified to be wide for all parameters ) srv_selex_prior ``` In the example provided above, the first row of the dataframe indicates that a prior is pplied to region 1, parameter 1, sex 1, fleet 1, and block 1, which corresponds to the $b_50$ parameter for logistic selectivity specified for the first survey fleet. The last row of the datframe indicates that a prior is applied to region 1, parameter 1, sex 2, fleet 2, and block 1, which corresponds to the male $\phi$ power paramter for exponential selectivity specified for the second survey fleet. Next, we can activate selectivity priors by setting `Use_srv_selex_prior = 1` and providing the associated selectivity prior dataframe to `srv_selex_prior`. ```{r} # Setup survey selectivity and catchability input_list <- Setup_Mod_Srvsel_and_Q(input_list = input_list, # Model options # survey selectivity, whether continuous time-varying cont_tv_srv_sel = c("none_Fleet_1", "none_Fleet_2", "none_Fleet_3"), # survey selectivity blocks srv_sel_blocks = c("none_Fleet_1", "none_Fleet_2", "none_Fleet_3"), # survey selectivity form srv_sel_model = c("logist1_Fleet_1", "exponential_Fleet_2", "logist1_Fleet_3"), # survey catchability blocks srv_q_blocks = c("none_Fleet_1", "none_Fleet_2", "none_Fleet_3"), # whether to estiamte all fixed effects for survey selectivity srv_fixed_sel_pars_spec = c("est_all", "est_all", "est_all"), # whether to estiamte all fixed effects for survey catchability srv_q_spec = c("est_all", "est_all", "est_all"), # Using selex priors Use_srv_selex_prior = 1, # Input survey selex priors srv_selex_prior = srv_selex_prior ) ``` # Movement Defining priors for movement, tag reporting rates, and recruitment proportions (see section below) is best illustrated with setting up a spatial model. Here, we will load in data from the multi-region sablefish case study to demonstrate defining priors for these processes. ```{r} data(mlt_rg_sable_data) # load in data # Initialize model dimensions and data list input_list <- Setup_Mod_Dim(years = 1:length(mlt_rg_sable_data$years), # vector of years (1 - 62) ages = 1:length(mlt_rg_sable_data$ages), # vector of ages (1 - 30) lens = mlt_rg_sable_data$lens, # number of lengths (41 - 99) n_regions = mlt_rg_sable_data$n_regions, # number of regions (5) n_sexes = mlt_rg_sable_data$n_sexes, # number of sexes (2) n_fish_fleets = mlt_rg_sable_data$n_fish_fleets, # number of fishery fleet (2) n_srv_fleets = mlt_rg_sable_data$n_srv_fleets, # number of survey fleets (2) verbose = TRUE ) ``` Because movement parameters are correlated and must sum to 1, `SPoRC` uses a Dirichlet prior to model movement probabilities. To illustrate, we set up the movement specification using the `Setup_Mod_Movement` function. Movement priors are enabled by setting `Use_Movement_Prior = 1`. The `Movement_prior` argument can be provided as a scalar (as in the example below), a vector, or a fully specified array. A scalar value, such as 1.5, implies a vague, symmetric Dirichlet prior that mildly penalizes values near the extremes (i.e., near 0 or 1). ```{r} # setting up movement parameterization input_list <- Setup_Mod_Movement(input_list = input_list, # Model options Movement_ageblk_spec = "constant", # estimating movement in 3 age blocks # (ages 1-6, ages 7-15, ages 16-30) Movement_yearblk_spec = "constant", # time-invariant movement Movement_sexblk_spec = "constant", # sex-invariant movement Use_Movement_Prior = 1, # priors used for movement Movement_prior = 1.5 # vague prior to penalize movement away from the extremes ) ``` If a vector or array is supplied to `Movement_prior`, it will be reshaped to match the full set of movement dimensions: ```{r} # array(Movement_prior, dim = c(input_list$data$n_regions, input_list$data$n_regions, # length(input_list$data$years), # length(input_list$data$ages), # input_list$data$n_sexes)) ``` It is important to ensure that the shape and values of the prior align with the structure of the movement parameters being estimated. Note that priors are only applied to the unique movement parameters. Thus, if movement is specified as constant across time, age, and sex, only one prior will be applied, regardless of the dimensionality of the array provided. # Tag Reporting Rates Defining priors for tag reporting rates is done by constructing a dataframe and passing it to the `Setup_Mod_Tagging` function. The dataframe must contain the following columns: `region`, `block`, `mu`, `sd`, and `type.` These columns define the region and time block the prior applies to, the mean and standard deviation of the prior distribution, and the type of beta prior to use. Two prior types are supported: a symmetric beta prior (`type = 0`), which penalizes extreme values near 0 and 1, and a standard beta prior (`type = 1`), which allows users to specify both the mean and standard deviation directly. In the example below, symmetric beta priors are applied to two time blocks in region 1. Because the symmetric form does not require a mean, the `mu` column is set to NA: ```{r} # setup tagging priors tag_prior <- data.frame( region = 1, # apply to region 1 block = c(1, 2), # two time blocks mu = NA, # not needed for symmetric beta sd = 5, # standard deviation type = 0 # symmetric beta prior ) tag_prior ``` In this example, the first row of the dataframe indicates that a symmetric beta prior is applied to tag reporting rates in region 1 and block 1, with a standard deviation of 5. The second row applies the same prior structure to block 2. Because `type = 0`, the prior is symmetric around 0.5 and discourages extreme reporting rate estimates. If `type = 1` were used instead, the mu column would need to be specified with the prior mean. Tag reporting rate priors can then be applied by setting `Use_TagRep_Prior = 1` and inputting the dataframe into the `TagRep_Prior` argument. ```{r} # setting up tagging parameterization input_list <- Setup_Mod_Tagging(input_list = input_list, UseTagging = 1, # using tagging data max_tag_liberty = 15, # maximum number of years to track a cohort # Data Inputs tag_release_indicator = mlt_rg_sable_data$tag_release_indicator, Tagged_Fish = mlt_rg_sable_data$Tagged_Fish, # Released fish Obs_Tag_Recap = mlt_rg_sable_data$Obs_Tag_Recap, # Model options Use_TagRep_Prior = 1, # tag reporting rate priors are used TagRep_Prior = tag_prior, # sex-specific data when fitting Init_Tag_Mort_spec = "fix", # fixing initial tag mortality Tag_Shed_spec = "fix", # fixing chronic shedding TagRep_spec = "est_shared_r", # tag reporting rates are # not region specific # Time blocks for tag reporting rates Tag_Reporting_blocks = c( paste("Block_1_Year_1-35_Region_", c(1:input_list$data$n_regions), sep = ''), paste("Block_2_Year_36-terminal_Region_", c(1:input_list$data$n_regions), sep = '') ), # Specify starting values or fixing values ln_Init_Tag_Mort = log(0.1), # fixing initial tag mortality ln_Tag_Shed = log(0.02) # fixing tag shedding ) ``` # Recruitment Proportions In the context of a spatial model, Dirichlet priors can also be placed on recruitment proportions. This process is relatively straightforward and uses the `Setup_Mod_Rec` functions to do so. In particular, users will need to define whether the prior is being used by setting `Use_Rec_prop_Prior = 1`. Following that, a scalar or vector of (`n_regions`) can be defined for `Rec_prop_prior`. In the case where `Rec_prop_prior` is a scalar, a uniform Dirichlet prior will be applied. If a vector of different values for each region are provided, then the Dirichlet prior will be applied according to those concentration parameters. ```{r} # Setup recruitment stuff (using defaults for other stuff) input_list <- Setup_Mod_Rec(input_list = input_list, # input data list from above do_rec_bias_ramp = 0, # Model options rec_model = "mean_rec", # recruitment model Use_Rec_prop_Prior = 1, # using recruitment proportion prior Rec_prop_prior = 5 # prior value for recruitment proportions ) ```