--- title: "Estimating ageing error using the TMB model written by André Punt" author: "Paul Burch and Ian Taylor based on documentation written by André Punt" date: "2025-03-21" output: html_document: default word_document: default vignette: > %\VignetteEngine{litedown::vignette} %\VignetteIndexEntry{Estimating ageing error using the TMB model written by André Punt} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` *Note: this vignette has not yet been updated to describe the new function `write_files()` which makes it easier to write the correctly-formatted input files, and `run()` which makes it easy to run the model and process the results. See the help pages for those functions for more detail.* # Introduction Uncertainty in age estimates needs to be accounted for when undertaking stock assessments. Punt et al. (2008) provides methods for estimating ageing error matrices that account for both ageing bias and imprecision. We provide some documentation and examples of the Template Model Builder (TMB) implementation these methods, incorporating documentation written by André Punt. This vignette is illustrated using data and settings associated with three stocks in the Australian Southern and Eastern Scalefish and Shark Fishery (SESSF) fishery (Blue Grenadier, Bight Redfish, and School Whiting), as well as sablefish from the U.S. west coast. Note: in version 2.1.0, released in March 2025, several function names were updated in the package: `CreateData()` was renamed to `load_data()`, `CreateSpecs()` was renamed to `load_specs()`. # R package The current implementation of the software was written by André Punt in 2020 and subsequently updated in 2022. There are several slightly different versions of the source code in circulation. This vignette describes the version hosted on the `pfmc-assessments` GitHub organization maintained by the assessment team at the NOAA Northwest Fisheries Science Center (NWFSC). The newer TMB version of the software replaced the ADMB version in the main branch of the repository in May 2024. The ADMB version is still available via the git history. We have consulted with André Punt and agreed that future development of the ageing error code should be done through the NWFSC GitHub repository. The easiest installation method is to use a pre-compiled version on R-universe via ```{r, echo=TRUE, eval=FALSE} install.packages("AgeingError", repos = c("https://noaa-fisheries-integrated-toolbox.r-universe.dev", "https://cloud.r-project.org")) ``` You can also install from GitHub using the code below which will compile the package locally (which takes longer and requires [Rtools](https://cran.r-project.org/bin/windows/Rtools/) for Windows users). ```{r, echo=TRUE, eval=FALSE} remotes::install_github("pfmc-assessments/AgeingError") ``` The package can then be loaded. ```{r} library(AgeingError) ``` The software requires the data and model are specified using input files which are described below. # Input files There are two input files needed to run the software the data file (.dat extension) and the specifications file (.spc extension). These are described below using Blue Grenadier (BG2022) as the example. ## Data file The data file (.dat extension) contains the age comparison data by pair of readers and specifies the range of ages permitted in the data, the minus and plus groups and the reference age. The data file comprises two or more sections (the age range then one section for each data pair of reader comparisons) The age range is simply specified by providing the minimum and maximum ages. The data can't fall outside the specified age range or the code will generate an error. The input files are parsed by the software using keywords that need to match exactly what is expected, so the text "Range_of_ages" is part of the input and not just a comment to help the user keep track of the format. The `write_data_file()` function makes it easier to create a file which matches the required format. ````{verbatim} Range_of_ages 0 25 ```` The data can be entered with a separate data set for each pair of readers, or a single data set for all pairs of readers. Both these options are described below. ### Separate Data Sets For each pair of readers (often the first and second reads by the same reader) there is one section with the header `Data_set_x` (x=1, 2, etc). The four lines immediately below the header are * The number of rows of paired age determinations (described below) * The number of readers in this data set (always 2 for the SESSF) * The minus group, plus group, and reference age * The reader numbers The remainder of the data set comprises one row of data for each pair of age determinations with the following format. The columns for each reader (1-6) represent an age, while the left most column represents the number of times that pair of ages were assigned by the two readers. The columns for the other readers (if any) need to be filled with `-999`. ````{verbatim} Data_set_1 141 # number of lines 2 # Number of readers 1 20 5 # minus group; plus group; reference age 1 2 3 4 # Which readers 9 0 0 -999 -999 1 0 10 -999 -999 153 1 1 -999 -999 25 1 2 -999 -999 ... 2 23 23 -999 -999 1 24 23 -999 -999 2 24 24 -999 -999 1 25 24 -999 -999 ```` A second (and potentially a third) data set follow the same format as the first data set, with the data being associated with the next pair of readers (see below). The headers for each data set need to have headers numbered from 1 to `NDataSet` which is an input to `load_data()`. That is, if `NDataSet=3`, the data file needs to have sections with header `Data_set_1`, `Data_set_2`, and `Data_set_3`. ````{verbatim} Data_set_2 123 # number of lines 2 # Number of readers 1 20 5 # minus group; plus group; reference age 3 4 1 2 # which readers 27 1 1 -999 -999 13 1 2 -999 -999 3 2 1 -999 -999 227 2 2 -999 -999 ... 2 22 20 -999 -999 1 22 21 -999 -999 1 23 21 -999 -999 1 23 22 -999 -999 ```` ### Combined Data Set The alternative to having separate data sets for each pair of readers is to combine multiple pairs of readers in a single data set as shown below. For this format the number of readers is now 4 and the data for the second pair of readers follows on from the end of the first pair. ````{verbatim} Data_set_1 264 # number of lines 4 # Number of readers 0 25 5 # minus group; plus group; reference age 1 2 3 4 # Which readers 9 0 0 -999 -999 1 0 10 -999 -999 153 1 1 -999 -999 30 1 2 -999 -999 ... 1 24 23 -999 -999 2 24 24 -999 -999 1 25 24 -999 -999 27 -999 -999 1 1 13 -999 -999 1 2 3 -999 -999 2 1 ... 2 -999 -999 22 20 1 -999 -999 22 21 1 -999 -999 23 21 1 -999 -999 23 22 ```` ### Other Considerations Note that the ages need to be whole numbers in the data file. The C++ code adds 0.5 to each age so the model it fitting to the midpoint of an age bin. Different agencies have created R code specific to the format of their ageing comparison data to format it for this package. For instance CSIRO uses the R package `SESSFdataproc` and hosts examples are in the assessment repository `SESSFassessments`. A version of this code is available from CSIRO. They usually create the data file in Excel (using the summarised output from the R code) and then paste into a text file. The data file is read into R using the function `load_data()` which is described in the examples below. ## Specifications file The age-estimates for at least one of the age-readers must be assumed to be unbiased (sensitivity to the choice of this age-reader should be examined if there is uncertainty about which age-reader is most likely to be unbiased). The options for bias and precision are described below and entered into the specifications file (.spc extension), an example of which is provided after the description of the model options. ### Bias options The available options for bias are: * `-x` Assume that the relationship between expected age and true age is the same as that for reader x. [Note that the ''x'' must be lower than the number of the reader for which bias is being defined. Reader 2 can have `BiasOpt=-1` to match reader 1]. * `0` Age-estimates are unbiased. * `1` The expected age of an animal of age $a$, $E_a$, is a linear function of its true age, i.e. $E_a =\alpha a$ (constant coefficient of variation). * `2` The expected age of an animal of age $a$, $E_a$, is a given by: $$E_a = E_L + (E_H - E_L)\frac{1-\text{exp}(-\beta(a-1))}{1-\text{exp}(-\beta(a_{max}-1))}$$. The C++ code (`AgeingError.cpp`) also has conditional (if) statements for Bias options 4 and 5 but these are identical to option 0 (unbiased) and therefore presumably placeholders for some future options. In the SESSF we typically assume that all readers are unbiased (`BiasOpt=0`), however, it is important to check this assumption using the residual plots (André Punt pers. comm.) ### Sigma options The available options for random ageing precision are: * `-x`. Assume that the relationship between the variance of age-reading error and true age is the same as that for reader x [note that ''x'' must be lower than the number of the reader for which bias is being defined]. * `1` This option assumes a constant coefficient of variation, it has one parameter that needs to be specified for each pair of independent readers in the specification (.spc) file. The standard deviation of random age-reading error, $\sigma_a$, is a linear function of true age, i.e.: $$ \sigma_a = \begin{cases} \gamma & \text{if } a=0 \\ \gamma a & \text{otherwise} \\ \end{cases} $$ * `2` The parameters relate to the standard deviation (Michaelis-Menten equation). This option has three parameters that need to be specified for each pair of independent readers in the specification (.spc) file. The standard deviation of random age-reading error, $\sigma_a$, is given by: $$ \sigma_a = \begin{cases} \sigma_L + (\sigma_H - \sigma_L)\frac{1}{1-exp(-\delta(a_{max}-1))} & \text{if } a=0 \\ \sigma_L + (\sigma_H - \sigma_L)\frac{1-exp(-\delta(a-1))}{1-exp(-\delta(a_{max}-1))} & \text{otherwise} \\ \end{cases} $$ * `3` The parameters relate to the coefficient of variation (Michaelis-Menten equation). This option has three parameters that need to be specified for each pair of independent readers in the specification (.spc) file. The coefficient of variation of random age-reading error, $CV_a$, is given by: $$ CV_a = \begin{cases} CV_L + (CV_H - CV_L)\frac{1}{1-\text{exp}(-\delta(a_{max}-1))} & \text{if } a=0 \\ CV_L + (CV_H - CV_L)\frac{1-\text{exp}(-\delta(a-1))}{1-\text{exp}(-\delta(a_{max}-1))} & \text{otherwise} \\ \end{cases} $$ * `4` The estimates of age are exact (use this option for ''known age'' individuals). This option returns a vector of zeroes (or perhaps just a single zero). * `5` The standard deviation of random age-reading error, $\sigma_a$, is a spline function of age (a Forsythe, Malcolm, and Moler type spline as implemented by `tmbutils::splinefun()`). The number and location (ages) of the knots must be specified, with the number of parameters equal to the number of knots. Bounds and initial parameters need to be specified in log space. * `6` The standard deviation of random age-reading error, $\sigma_a$, is a piecewise linear function of age. The number and location (ages) of the knots must be specified, with the number of parameters equal to the number of knots. Bounds and initial parameters need to be specified in log space. * `7` A linear change in the standard deviation of random age-reading error, $\sigma_a$, with age. This option has two parameters that need to be specified for each pair of independent readers in the specifications (.spc) file. * `8` A linear change in the coefficient of variation of random age-reading error, $CV_a$ , with age. This option has two parameters that need to be specified for each pair of independent readers in the specifications (.spc) file. ### Example Specifications file The .spc file contains three or four sections: 1. reader specification, 2. knots for splines or breakpoints for piecewise linear functions (if used), 3. bias parameters, and 4. the sigma parameters. The first section specifies the Bias and Sigma options for each reader. For each reader, one line specifies their Bias and Sigma options. `SigmaOpt=-x` is used to assume the relationship between expected age and true age is the same as that for reader x [note that the `x` must be lower than the number of the reader for which bias is being defined]. The example below (from Blue Grenadier) shows two pairs of readers (1-2 & 3-4) who are all specified to be unbiased (`BiasOpt=0`) with readers 1 and 2 having the same ageing precision `SigmaOpt=2`. This model requires three parameters to be specified (presumably $\sigma_L$, $\sigma_H$, and $\delta$) for each pair of readers (in the case of Blue Grenadier the first and second reads from two readers). ````{verbatim} # reader BiasOpt SigmaOpt 1 0 2 2 0 -1 3 0 2 4 0 -3 ```` The second (optional) section species the knots for the spline and linear interpolation models. The example below is taken from School Whiting (`WHS2.spc`). It specifies five knots with the first and last knots being the minimum and maximum ages. Note that the text "# Spline specifications" or "# Linear specifications" are required by the code if SigmaOpt options 5 or 6 are used. ````{verbatim} # Spline specifications 5 0 3 5 7 9 # Linear specifications 5 0 3 5 7 9 ```` The next section specifies conditions under which the Bias parameters are estimated. We need to provide the bounds, initial values and whether the parameter should be estimated (1) or pre-specified (0) at the initial value. If `Bias=0` (the case for SESSF stocks) this section is blank (NULL), although we retain the heading (not sure if the heading is needed). The Sablefish example (Sable.spc) provides an example of how to implement bias within ageing error estimation. ````{verbatim} Bias_Pars (low high init, on/off) ```` The final section specifies conditions under which the Sigma parameters are estimated. Like above, we need to provide the bounds, initial values and whether the parameter should be estimated (1). One set of parameter bounds, initial values and on/off needs to be provide for each group of readers that assume the same ageing uncertainty (group of `x/-x`). In the example below (`SigmaOpt=2`) there are three parameters for each pair of readers. ````{verbatim} Sigma_Pars (low high init, on/off) 0 1 0.2 1 # Readers 1 & 2 0.01 1 0.2 1 0 2 0.5 1 0 1 0.2 1 # Readers 3 & 4 0.01 1 0.2 1 0 2 0.5 1 ```` # Examples For the examples we focus on SESSF species, which all specify (`BiasOpt=0`). If bias in the age readings is suspected, an example of how to implement bias within ageing error estimation is provided in the Sablefish example (Sable.spc). ## Blue Grenadier Blue Grenadier (_Macruronus novaezelandiae_) is deep water species caught by trawl in south eastern Australia (Tuck and Bessell-Browne 2022). It has age readings from two separate readers who have each re-read their own reads (i.e. there are no inter-reader comparisons in the production data). In 2022, estimating ageing error for Blue Grenadier was challenging, requiring removal of problematic data, different models for each reader and improvements to the model source code. We work through this example below. Note we are using the data file with a single data set so we set `NDataSet=1`. ```{r} data_dir <- system.file("extdata", package = "AgeingError") BG2022_dat <- AgeingError::load_data(file.path(data_dir, "BG2022.dat"), NDataSet = 1, verbose = TRUE, EchoFile = "BG2022echo.out" ) ``` When we set `verbose=TRUE` a summary of the loaded data is printed to the console. Check that the this matches with the input file, some things that are worth checking include, * The last line of each data set. * The total number reads for each reader pair. * The number of readers and their order. * The plus, minus and reference ages. There's a warning about potentially missing data, this occurs when there are ages below the minimum age and above the maximum age, however, the warning does not appear to be correct in this case. Possibly the `-999` values in the data are causing the warning because they are <0. Next we load the model specifications using the `load_specs` function. ```{r} BG2022_spc <- AgeingError::load_specs(file.path(data_dir, "BG2022.spc"), DataSpecs = BG2022_dat, verbose = TRUE ) ``` The function creates a nested list with one list element for each reader in the data. By setting `verbose=TRUE`, the list will print to the console and you can check that the inputs have been read in correctly. Note a reader that has a `-x` for the Sigma option (i.e. the relationship between expected age and true age is the same as that for reader x) will not have parameter values. The output from this function varies with the model that is specified (i.e. the knots will all be zero unless a spline or piecewise linear model is specified - sigma options 5 and 6). Ageing error is the estimated using the function `DoApplyAgeError`. We need to provide a name (`Species` argument), the data, model specs, specify the weights (`AprobWght` and `SlopeWght`) and the name of a subdirectory to save the output (it gets created if it doesn't exist). We set `verbose=FALSE` in this example as a large amount of output is printed to the screen. ```{r, results="hide"} BG2022_mod <- AgeingError::DoApplyAgeError( Species = "BG2022", DataSpecs = BG2022_dat, ModelSpecsInp = BG2022_spc, AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results", verbose = FALSE ) ``` The model is stored as a list in the object `BG2022_mod` and this information along with the data and model specifications is saved in the output directory as an lda file with name `Species` (in this case "BG2022.lda"). Examine the structure of the model object ```{r} str(BG2022_mod) ``` We then save the model output using the `ProcessResults` function. We want to specify that the effective sample size is estimated (`CalcEff = TRUE`) as this produces the fits to the individual data points which can be used to identify outliers that may be impacting model convergence. The model creates several output files that contain the ageing error estimates and the information we need to assess the fit so we set `verbose = FALSE`. This following files in the output subdirectory (in this case `Results`). * The report file (.rpt extension) which contains the convergence criteria, model specifications, parameter estimates with their uncertainty and ageing error matrices for each reader (note readers with bias or sigma option `-x` will be identical to reader `x`). * An LDA file (.lda extension) containing the data and model specifications. * One csv file (.csv extension) for each reader that contains the estimated standard deviation and coefficient of variation for each age (note readers with bias or sigma option `-x` will be identical to reader `x`). * A series png file (.png extension) that contain plots of the data and model residuals. ```{r} BG2022_out <- AgeingError::ProcessResults(Species = "BG2022", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE) ``` An examination of the report file (`BG2022.rpt`) shows the maximum gradient is -93.96, indicating the model has not converged. Ideally the gradient should be between $\pm \text{1e-4}$ and zero, however, in practice achieving this for some models can be very time consuming and involves removing progressively more data so we sometimes accept a gradient ~$\pm \text{1e-3}$. The residual plots show a large number of outliers exceeding the 95% confidence interval and the table of fits to the individual data points (bottom of the report file (only present when `CalcEff = TRUE`) shows a number of data points with predicted values of 0 (or very close to zero e.g. 10e-10), see line 2 below (extracted from the report file). ````{verbatim} Data set: 1 Data_set Group Group Line Readers Obs Obs_Numbers Pred_Numbers Data Point: 1 1 1 1 0 0 -1 -1 9 15.6581407 9.241041710732 Data Point: 1 1 1 2 0 10 -1 -1 1 1.7397934 0 Data Point: 1 1 1 3 1 1 -1 -1 153 266.1883915 186.217154197499 ```` To get closer to a converged model it was necessary to undertake data trimming (removing observations values with predicted values of 0), specify a different ageing error model for the second reader and prespecify (i.e. fix) one of the parameters for the second reader. This was done in a stepwise manner, with small changes being made, the model run and the results examined before more changes were made until we achieved a final model we considered acceptable. We refit the final model (`BG2022_trim_8_1`). Load the data file. ```{r} BG2022final_dat <- AgeingError::load_data(file.path(data_dir, "BG2022_trim_8_1.dat"), NDataSet = 1, verbose = TRUE, EchoFile = "BG2022finalecho.out" ) ``` Note the number of lines of data has been reduced from 264 to 252. Next we load the model specifications. For the second reader (3 and 4 in the data) we have specified a linear change in the standard deviation of random age-reading error (`SigmaOpt=7`) and fixed the first parameter at 0.2. ```{r} BG2022final_spc <- AgeingError::load_specs(file.path(data_dir, "BG2022_trim_8_1.spc"), DataSpecs = BG2022final_dat, verbose = TRUE ) ``` We run the final model. ```{r, results="hide"} BG2022final_mod <- AgeingError::DoApplyAgeError( Species = "BG2022_trim_8_1", DataSpecs = BG2022final_dat, ModelSpecsInp = BG2022final_spc, AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results", verbose = FALSE ) ``` As with the earlier model we save the results. ```{r} BG2022final_out <- AgeingError::ProcessResults(Species = "BG2022_trim_8_1", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE) ``` This model is now close to convergence (gradient = -0.00143) and residual plots do not show any concerning features and the estimates of ageing error standard deviation are increasing with age. ## Bight Redfish Bight Redfish (_Centroberyx gerrardi_) is a long-lived species caught by trawl in the Great Australian Bight (Curin-Osorio and Burch 2022). It has age readings from three separate readers who have each re-read their own reads (i.e. there are no inter-reader comparisons in the data). ```{r} REB2022_dat <- AgeingError::load_data(file.path(data_dir, "REB2022.dat"), NDataSet = 1, verbose = TRUE, EchoFile = "REB2022echo.out" ) ``` We check the summary above against the information in the data file. There don't appear to be any problems, however, we are still getting the warning about missing data. Next we load the model specifications using the `load_specs` function. ```{r} REB2022_spc <- AgeingError::load_specs(file.path(data_dir, "REB2022.spc"), DataSpecs = REB2022_dat, verbose = TRUE ) ``` For Bight Redfish we fit a one parameter model assuming the expected age is a linear function of the true age. We assume there is no inter-reader variability (all readers have the same ageing uncertainty as reader 1). Note that because we have specified all readers have the same ageing uncertainty (readers 2--6 have `SigmaOpt=-1`), we only need to provide one line of sigma parameters. We estimate ageing error for Bight Redfish. ```{r, results="hide"} REB2022_mod <- AgeingError::DoApplyAgeError( Species = "REB2022", DataSpecs = REB2022_dat, ModelSpecsInp = REB2022_spc, AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results", verbose = FALSE ) ``` We then save the Bight Redfish model output. ```{r} REB2022_out <- AgeingError::ProcessResults(Species = "REB2022", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE) ``` An examination of the report file (`REB2022.rpt`) shows the maximum gradient is 0.0221, indicating the model has not quite converged. We remove four outliers with probability <1e-5 and also start the model at the estimated value of the CV (0.0413288). We save the revised data and model specifications as (`REB2022_trim.dat` and `REB2022_trim.spc`) and rerun the model (note we suppress the creation of output in this document). ```{r, results="hide"} REB2022_dat2 <- AgeingError::load_data(file.path(data_dir, "REB2022_trim.dat"), NDataSet = 1, verbose = FALSE, EchoFile = "REB2022_trimecho.out" ) REB2022_spc2 <- AgeingError::load_specs(file.path(data_dir, "REB2022_trim.spc"), DataSpecs = REB2022_dat2, verbose = FALSE ) REB2022_mod2 <- AgeingError::DoApplyAgeError( Species = "REB2022_final", DataSpecs = REB2022_dat2, ModelSpecsInp = REB2022_spc2, AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results", verbose = FALSE ) REB2022_out2 <- AgeingError::ProcessResults( Species = "REB2022_final", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE ) ``` The final gradient is -0.001973, close to convergence. ## School Whiting School Whiting (_Sillago flindersi_) is a short-lived species caught by Danish seine and trawl in south eastern Australia (Day et al. 2020). This species provides an example of a data file with three separate data sets and a spline model. Load the School Whiting data file, make sure to set the number of data sets to three (`NDataSet=3`). ```{r} WHS2_dat <- AgeingError::load_data(file.path(data_dir, "WHS2.dat"), NDataSet = 3, verbose = TRUE, EchoFile = "WHS2echo.out" ) ``` Make sure to check the summary above against the values in the data file. Load the specifications for School Whiting. ```{r} WHS2_spc <- AgeingError::load_specs(file.path(data_dir, "WHS2.spc"), DataSpecs = WHS2_dat, verbose = TRUE ) ``` For School Whiting we fit a spline function (`SigmaOption=5`) with five knots (0, 3, 5, 7, and 9). Note the initial parameter values are in log space. Fit the ageing error model and save the results ```{r, results="hide"} WHS2_mod <- AgeingError::DoApplyAgeError( Species = "WHS2", DataSpecs = WHS2_dat, ModelSpecsInp = WHS2_spc, AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results", verbose = FALSE ) WHS2_out <- AgeingError::ProcessResults( Species = "WHS2", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE ) ``` The maximum gradient is ~2.3-05 indicating the model has converged. The estimated standard deviation by age shows the ageing uncertainty increasing from age 0 to age 4, then declining for ages 5 and 6 before increasing substantially for ages 7--9. The fits at the far right hand side of the age distribution are poor and there is little data to inform the model in this region. The appropriateness of this model should be discussed with the ageing technicians. ## Sablefish Sablefish (_Anoplopoma fimbria_) provides an example of using the bias option. ```{r} Sable_dat <- AgeingError::load_data(file.path(data_dir, "Sable.dat"), NDataSet = 1, verbose = TRUE, EchoFile = "SableEcho.out" ) ``` We check the summary above against the information in the data file. This data set has several inter-reader comparisons. We load the model specifications using the `load_specs` function. ```{r} Sable_spc <- AgeingError::load_specs(file.path(data_dir, "Sable.spc"), DataSpecs = Sable_dat, verbose = TRUE ) ``` The Sablefish example fits three separate series (1-2, 3-4 & 5-6). The model assumes the following: Bias Options: Reader 5 is unbiased, readers 1 & 2 have the same constant bias, reader 6 has a separate constant bias while the bias of readers 3 & 4 is the same and follows a Michaelis-Menten relationship. In the specifications file the first list relates to the constant bias of readers 1 & 2, the next three lines readers 3 & 4, while the final line specifies the constant bias of reader 6. Reader 5 is assumed to be unbiased and therefore doesn't require parameters to be specified. Sigma Options: Readers 1 & 2 assume the same linear relationship (constant CV), while readers 5 & 6 have a separate linear relationship. Readers 3 & 4 assume ageing uncertainty follows a Michaelis-Menten relationship. The parameters are specified in the same manner as for the bias option (above). The code below can be used to fit the Sablefish ageing error model. We don't run it because this particular model is very slow. ```{r, eval=FALSE} ## run the Sablefish model Sable_mod <- AgeingError::DoApplyAgeError( Species = "Sable", DataSpecs = Sable_dat, ModelSpecsInp = Sable_spc, AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results", verbose = FALSE ) ## save the model results Sable_out <- AgeingError::ProcessResults(Species = "Sable", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE) ``` # Trouble Shooting Some suggestions for getting to a converged model are provided below. * Vary the initial values of the parameters in the specifications file. * For models that are close to convergence, re-run the model using the estimated parameters as the initial values. * Problematic data can be identified from the residual plots and at the bottom of the report file with prediction probabilities <1e-5. Removing data with low probability can improve the fit. * There may not be sufficient data to estimate ageing error for multiple combinations of readers, the data can be pooled to estimate fewer parameters (bias / sigma option `-x`) * The model may not be appropriate, select a different model based on the residual plots. * Prespecify (fix) one or more parameters for a pair of readers. # Acknowledgements Fish Ageing Services https://www.fishageingservices.com/ provided the data for Blue Grenadier, Bight Redfish, and School Whiting and the Australian Fisheries Management Authority (AFMA) provided funding for the data collection and otolith reading. The sablefish age data are provided by the NWFSC. This document is prepared from documentation and examples written by André Punt. Pia Bessell-Browne provided helpful comments that improved this document. # References Day, J., Hall, K., Bessell-Browne, P., and Sporcic, M. (2020) School Whiting (_Sillago flindersi_) stock assessment based on data up to 2019. For discussion at SERAG, December 2020. Curin-Osorio, S., and Burch, P. (2022). Bight Redfish (_Centroberyx gerrardi_) stock assessment based on data up to 2021-22. Technical paper presented to the GABRAG, 22 November 2022, Hobart, Tasmania. Punt, A.E., Smith, D.C., KrusicGolub, K., and Robertson, S. (2008). Quantifying age-reading error for use in fisheries stock assessments, with application to species in Australia's southern and eastern scalefish and shark fishery. _Canadian Journal of Fisheries and Aquatic Sciences_ 65: 1991--2005. https://doi.org/10.1139/F08-111 Tuck, G.N., Bessell-Browne, P. (2022). Blue Grenadier (_Macruronus novaezelandiae_) stock assessment based on data up to 2021. Technical paper presented to theSERAG2, 29--30$^{th}$ November 2022, Melbourne, Victoria. 99pp.