-
Notifications
You must be signed in to change notification settings - Fork 4
Example cat 1 stock assessment
See also: Getting set up, Running a TAF analysis, Creating TAF datasets.
This page was based on using the icesTAF package version 4.2.0 dated
2023-03-21.
This example requires packages to be installed that is not on CRAN.
Please ensure you have the package stockasessment by running
# Download and install stockassessment in R
install.packages("stockassessment", repos = "https://fishfollower.r-universe.dev")and the FLfse package by running
# Download and install stockassessment in R
install.packages("FLfse", repos = "https://ices-tools-prod.r-universe.dev")The resulting TAF analysis created from this example is on GitHub: github.com/ices-taf-dev/wiki-example-stockassessment, if you dont want to build up the code yourself by following the example, you can skip to the end and run on your computer by downloading the complete code and running it like this:
# get code
run_dir <- download.analysis("ices-taf-dev/wiki-example-stockassessment")
# view downloaded code
browseURL(run_dir)
# run analysis
run.analysis(run_dir)
# view result
browseURL(run_dir)- Creating an empty stockassessment.org TAF project
- Preprocess the data
- Running a model
- Writing TAF tables
- Formatted output for reporting
First we create an empty TAF project set up with a boot section for a
stock assessment created on stock
stockassessment.org/, and in this case
we will call it example-4 and choose the run called:
WBCod_2021_cand01, and then moving our working directory to this
folder. We do this by running
taf.skeleton.sa.org("example-4", "Cod_27_5b1")## To run this template please ensure you have the package:
## stockasessment
## by running:
## install.packages("stockassessment", repos = "https://fishfollower.r-universe.dev")
setwd("example-4")resulting in the following:
example-4
¦--boot
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data.R
¦--model.R
¦--output.R
°--report.R and after running
taf.boot()## [09:16:56] Boot procedure running...
## Processing DATA.bib
## [09:16:56] * sam_data
## [09:16:58] * sam_config
## [09:16:58] Boot procedure done
your project should now look like this:
example-4
¦--boot
¦ ¦--data
¦ ¦ ¦--sam_config
¦ ¦ ¦ °--model.cfg
¦ ¦ °--sam_data
¦ ¦ ¦--cn.dat
¦ ¦ ¦--cw.dat
¦ ¦ ¦--dw.dat
¦ ¦ ¦--lf.dat
¦ ¦ ¦--lw.dat
¦ ¦ ¦--mo.dat
¦ ¦ ¦--nm.dat
¦ ¦ ¦--pf.dat
¦ ¦ ¦--pm.dat
¦ ¦ ¦--survey.dat
¦ ¦ °--sw.dat
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data.R
¦--model.R
¦--output.R
°--report.R Now we have a input data and the fitted model object that we can use to create some standard TAF files.
Now that you have created the file structure and uploaded your data,
it’s time to turn your attention to the data.R script. The purpose of
this script is to do some data processing. This could be (among other
things):
- calculating a plus group,
- taking averages of survey indices,
- calculating fill-in values and removing outliers,
- calculating smoothed time series
The data.R script should also write out input data into flat .csv
files, so they are readable and reviewable by others.
In this section we will build the script up in sections. First we need to load the libraries and get the data required for the script. It is good practice to write your analysis in several scripts with each script being as independent as possible. This will allow you to return to your analysis and work on a specific section, and not worry too much about whether you have loaded the correct libraries and data. So the first step is always to load libraries, then load the data required.
It is also considered good practice to create (or ensure) that the output directory for the TAF section you are working in exists. Hence the start of your script will be structured like this:
## Preprocess data, write TAF data tables
## Before:
## After:
## load libraries
library(icesTAF)
library(stockassessment)
# ensure directory
mkdir("data")
# Read underlying data from boot/data
# ...For the example we are working with here this becomes
## Preprocess data, write TAF data tables
## Before:
## After:
## load libraries
library(icesTAF)
library(stockassessment)
# ensure directory
mkdir("data")
# Read underlying data from boot/data
# quick utility function
read.ices.taf <- function(...) {
read.ices(taf.data.path("sam_data", ...))
}
# ## Catch-numbers-at-age ##
catage <- read.ices.taf("cn.dat")
# ## Catch-weight-at-age ##
wcatch <- read.ices.taf("cw.dat")
wdiscards <- read.ices.taf("dw.dat")
wlandings <- read.ices.taf("lw.dat")
# ## Natural-mortality ##
natmort <- read.ices.taf("nm.dat")
# maturity
maturity <- read.ices.taf("mo.dat")
# ## Proportion of F before spawning ##
propf <- read.ices.taf("pf.dat")
# ## Proportion of M before spawning ##
propm <- read.ices.taf("pm.dat")
# ## Stock-weight-at-age ##
wstock <- read.ices.taf("sw.dat")
# Landing fraction in catch at age
landfrac <- read.ices.taf("lf.dat")
# survey data
surveys <- read.ices.taf("survey.dat")The next step is to do any preprocessing or calculations that you may need to do. In this simple example we don’t need to do much
## 2 Preprocess data
# landings
latage <- catage * landfrac
datage <- catage * (1 - landfrac)
# put surveys in seperate matrices (for writing out)
survey_summer <- surveys[[1]]
survey_spring <- surveys[[2]]And finally when we are done, it is recommended that you write out your input data in a human (easily) readable form. CSV is always a good option as it is easily read into R and Excel etc.
## 3 Write TAF tables to data directory
write.taf(
c(
"catage", "latage", "datage", "wstock", "wcatch",
"wdiscards", "wlandings", "natmort", "propf", "propm",
"landfrac", "survey_summer", "survey_spring"
),
dir = "data"
)
## write model files
dat <- setup.sam.data(
surveys = surveys,
residual.fleet = catage,
prop.mature = maturity,
stock.mean.weight = wstock,
catch.mean.weight = wcatch,
dis.mean.weight = wdiscards,
land.mean.weight = wlandings,
prop.f = propf,
prop.m = propm,
natural.mortality = natmort,
land.frac = landfrac
)
save(dat, file = "data/data.RData")And this concludes the data script. To test the script you can run
sourceTAF("data")
sourceTAF("data")## [09:16:58] data.R running...
## [09:16:59] data.R done
your project should now look like this:
example-4
¦--boot
¦ ¦--data
¦ ¦ ¦--sam_config
¦ ¦ ¦ °--model.cfg
¦ ¦ °--sam_data
¦ ¦ ¦--cn.dat
¦ ¦ ¦--cw.dat
¦ ¦ ¦--dw.dat
¦ ¦ ¦--lf.dat
¦ ¦ ¦--lw.dat
¦ ¦ ¦--mo.dat
¦ ¦ ¦--nm.dat
¦ ¦ ¦--pf.dat
¦ ¦ ¦--pm.dat
¦ ¦ ¦--survey.dat
¦ ¦ °--sw.dat
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data
¦ ¦--catage.csv
¦ ¦--data.RData
¦ ¦--datage.csv
¦ ¦--landfrac.csv
¦ ¦--latage.csv
¦ ¦--natmort.csv
¦ ¦--propf.csv
¦ ¦--propm.csv
¦ ¦--survey_spring.csv
¦ ¦--survey_summer.csv
¦ ¦--wcatch.csv
¦ ¦--wdiscards.csv
¦ ¦--wlandings.csv
¦ °--wstock.csv
¦--data.R
¦--model.R
¦--output.R
°--report.R Although in this example the model has already been run elsewhere, we
can still choose to rerun the model or even conduct further analyses.
Here we load the model from the boot/data folder and run a
retrospective analysis. then save the results in the model folder.
## Run analysis, write model results
## Before:
## After:
# load libraries
library(icesTAF)
library(stockassessment)
# ensure model dir exists
mkdir("model")
# load data required
load("data/data.RData", verbose = TRUE)
# setup configuration
conf <-
loadConf(
dat,
"boot/data/sam_config/model.cfg",
patch = TRUE
)
# define parameters
par <- defpar(dat, conf)
# fit model
fit <- sam.fit(dat, conf, par)
# retrospective fits
retro_fit <- retro(fit, year = 5)
# save model fits
save(fit, file = "model/fit.RData")
save(retro_fit, file = "model/retro_fit.RData")And this concludes the model script. To test the script you can run
sourceTAF("model")
your project should now look like this:
example-4
¦--boot
¦ ¦--data
¦ ¦ ¦--sam_config
¦ ¦ ¦ °--model.cfg
¦ ¦ °--sam_data
¦ ¦ ¦--cn.dat
¦ ¦ ¦--cw.dat
¦ ¦ ¦--dw.dat
¦ ¦ ¦--lf.dat
¦ ¦ ¦--lw.dat
¦ ¦ ¦--mo.dat
¦ ¦ ¦--nm.dat
¦ ¦ ¦--pf.dat
¦ ¦ ¦--pm.dat
¦ ¦ ¦--survey.dat
¦ ¦ °--sw.dat
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data
¦ ¦--catage.csv
¦ ¦--data.RData
¦ ¦--datage.csv
¦ ¦--landfrac.csv
¦ ¦--latage.csv
¦ ¦--natmort.csv
¦ ¦--propf.csv
¦ ¦--propm.csv
¦ ¦--survey_spring.csv
¦ ¦--survey_summer.csv
¦ ¦--wcatch.csv
¦ ¦--wdiscards.csv
¦ ¦--wlandings.csv
¦ °--wstock.csv
¦--data.R
¦--model
¦--model.R
¦--output.R
°--report.R In the output section we take the model and any aditional analyses and create human readable, but full precision, outputs. This is much like the data section, but we are dealing with model outputs.
## Extract results of interest, write TAF output tables
## Before:
## After: csv tables of assessment output
library(icesTAF)
library(stockassessment)
mkdir("output")
# load fit
load("model/fit.rData", verbose = TRUE)
# Model Parameters
partab <- partable(fit)
# Fs
fatage <- faytable(fit)
fatage <- fatage[, -1]
fatage <- as.data.frame(fatage)
# Ns
natage <- as.data.frame(ntable(fit))
# Catch
catab <- as.data.frame(catchtable(fit))
colnames(catab) <- c("Catch", "Low", "High")
# TSB
tsb <- as.data.frame(tsbtable(fit))
colnames(tsb) <- c("TSB", "Low", "High")
# Summary Table
tab.summary <- cbind(as.data.frame(summary(fit)), tsb)
tab.summary <- cbind(tab.summary, catab)
# should probably make Low and High column names unique R_Low etc.
mohns_rho <- mohn(retro_fit)
mohns_rho <- as.data.frame(t(mohns_rho))
## Write tables to output directory
write.taf(
c("partab", "tab.summary", "natage", "fatage", "mohns_rho"),
dir = "output"
)This concludes the ouput script. To test the script you can run
sourceTAF("output")
sourceTAF("output")## [09:16:59] output.R running...
## Warning: replacing previous import 'FLCore::weighted.mean' by 'stats::weighted.mean' when loading 'FLfse'
## Warning in readChar(con, 5L, useBytes = TRUE): cannot open compressed file 'model/fit.rData', probable reason 'No such file or directory'
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## [09:17:00] output.R failed
your project should now look like this:
example-4
¦--boot
¦ ¦--data
¦ ¦ ¦--sam_config
¦ ¦ ¦ °--model.cfg
¦ ¦ °--sam_data
¦ ¦ ¦--cn.dat
¦ ¦ ¦--cw.dat
¦ ¦ ¦--dw.dat
¦ ¦ ¦--lf.dat
¦ ¦ ¦--lw.dat
¦ ¦ ¦--mo.dat
¦ ¦ ¦--nm.dat
¦ ¦ ¦--pf.dat
¦ ¦ ¦--pm.dat
¦ ¦ ¦--survey.dat
¦ ¦ °--sw.dat
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data
¦ ¦--catage.csv
¦ ¦--data.RData
¦ ¦--datage.csv
¦ ¦--landfrac.csv
¦ ¦--latage.csv
¦ ¦--natmort.csv
¦ ¦--propf.csv
¦ ¦--propm.csv
¦ ¦--survey_spring.csv
¦ ¦--survey_summer.csv
¦ ¦--wcatch.csv
¦ ¦--wdiscards.csv
¦ ¦--wlandings.csv
¦ °--wstock.csv
¦--data.R
¦--model
¦--model.R
¦--output
¦--output.R
°--report.R It is often of use to split up the reporting into several scripts. This
can be done in any section, but is most common in the report section. Do
do this, you create scripts called report_*.R and your report.R
script becomes a like a recipe, that sources each relavent script in
turn. In this example we are going to make three seperate scripts, and
an R markdown script to create a report:
- report_plots.R
- report_tables.R
This means that the central report.R script will look like this:
## Prepare plots and tables for report
## Before:
## After:
library(icesTAF)
mkdir("report")
sourceTAF("report_plots.R")
sourceTAF("report_tables.R")
sourceTAF("report_doc.R")The contents of these scripts are very much incomplete and are here just
to serve as an example. A simple example of the report_plots.R script
is:
library(icesTAF)
library(stockassessment)
mkdir("report")
load("model/fit.rData")
load("model/retro_fit.rData")
## input data plots
## ....
## model output plots ##
taf.png("summary", width = 1600, height = 2000)
plot(fit)
dev.off()
taf.png("SSB")
ssbplot(fit, addCI = TRUE)
dev.off()
taf.png("Fbar")
fbarplot(fit, xlab = "", partial = FALSE)
dev.off()
taf.png("Rec")
recplot(fit, xlab = "")
dev.off()
taf.png("Landings")
catchplot(fit, xlab = "")
dev.off()
taf.png("retrospective", width = 1600, height = 2000)
plot(retro_fit)
dev.off()and this can be run and tested using sourceTAF("report_plots"), or
simply source("report_plots.R")
A simple report_tables.R script is
mkdir("report")
load("model/fit.RData", verbose = TRUE)
years <- unique(fit$data$aux[, "year"])
## catage
catage <- read.taf("data/catage.csv")
# row.names(catage) <- years[1:nrow(catage)]
catage <- cbind(catage, total = rowSums(catage))
catage <- rbind(catage, mean = colMeans(catage))
write.taf(catage, "report/catage.csv")In completing this section we will create three more scripts:
- report_doc.R
- report.Rmd
- utilities.R
And add an extra record to the DATA.bib file, which will get a
template word document that will be used in the report.Rmd to produce
a work document with preset formatting.
And the report_doc.R script will generally look like this
library(rmarkdown)
source("utilities.R")
mkdir("report")
# combine into a word and html document
render("report.Rmd",
output_format = c("word_document", "html_document"),
encoding = "UTF-8"
)
# move to report folder
cp("report.html", "report", move = TRUE)
cp("report.docx", "report", move = TRUE)The utilities script is a place you can keep useful functions. In this exmaple it is a function to help make nice tables, and contains just one function
# utilities.R
style_table1 <- function(tab) {
# Capitalize first letter of column, make header, last column and second
# last row in boldface and make last row italic
names(tab) <- pandoc.strong.return(names(tab))
emphasize.strong.cols(ncol(tab))
emphasize.strong.rows(nrow(tab))
set.alignment("right")
return(tab)
}The report.Rmd script is very much an example and not a suggestion.
(more to come later).
The code is found here: report.Rmd
and this requires an additional file in the boot data folder, a word document that acts like a template for formats and styles. To add this run:
draft.data(
data.scripts = NULL,
data.files = "reportTemplate.docx",
originator = "ICES",
year = 2019,
title = "ICES TAF word template for report automation",
source = "https://github.com/ices-taf/doc/raw/master/reportTemplate.docx",
append = TRUE,
file = TRUE
)
setwd("..")and then run taf.boot() again to fetch the word template into the
boot/data folder
taf.boot()## [09:17:00] Boot procedure running...
## Processing DATA.bib
## [09:17:00] * sam_data
## Skipping script 'sam_data.R' (directory contains files already).
## [09:17:00] * sam_config
## Skipping script 'sam_config.R' (directory contains files already).
## [09:17:00] * reportTemplate.docx
## [09:17:01] Boot procedure done
your project should now look like this:
example-4
¦--boot
¦ ¦--data
¦ ¦ ¦--reportTemplate.docx
¦ ¦ ¦--sam_config
¦ ¦ ¦ °--model.cfg
¦ ¦ °--sam_data
¦ ¦ ¦--cn.dat
¦ ¦ ¦--cw.dat
¦ ¦ ¦--dw.dat
¦ ¦ ¦--lf.dat
¦ ¦ ¦--lw.dat
¦ ¦ ¦--mo.dat
¦ ¦ ¦--nm.dat
¦ ¦ ¦--pf.dat
¦ ¦ ¦--pm.dat
¦ ¦ ¦--survey.dat
¦ ¦ °--sw.dat
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data
¦ ¦--catage.csv
¦ ¦--data.RData
¦ ¦--datage.csv
¦ ¦--landfrac.csv
¦ ¦--latage.csv
¦ ¦--natmort.csv
¦ ¦--propf.csv
¦ ¦--propm.csv
¦ ¦--survey_spring.csv
¦ ¦--survey_summer.csv
¦ ¦--wcatch.csv
¦ ¦--wdiscards.csv
¦ ¦--wlandings.csv
¦ °--wstock.csv
¦--data.R
¦--model
¦--model.R
¦--output
¦--output.R
¦--report.R
¦--report.Rmd
¦--report_doc.R
¦--report_plots.R
¦--report_tables.R
°--utilities.R Once all this code is in place the whole report section can be run
sourceTAF("report")
sourceTAF("report")## [09:17:01] report.R running...
## [09:17:01] report_plots.R running...
## Warning in readChar(con, 5L, useBytes = TRUE): cannot open compressed file 'model/fit.rData', probable reason 'No such file or directory'
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## [09:17:01] report_plots.R failed
## [09:17:01] report_tables.R running...
## Warning in readChar(con, 5L, useBytes = TRUE): cannot open compressed file 'model/fit.RData', probable reason 'No such file or directory'
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## [09:17:01] report_tables.R failed
## [09:17:01] report_doc.R running...
##
##
## processing file: report.Rmd
## 1/17
## 2/17 [libraries]
## 3/17
## 4/17 [chunk_setup]
## 5/17
## 6/17 [pander_settings]
## 7/17
## 8/17 [caption_counters]
## 9/17
## 10/17 [catch_table]
##
## Quitting from lines 61-68 [catch_table] (report.Rmd)
## Error in file(file, "rt", encoding = fileEncoding) :
## cannot open the connection
## [09:17:01] report_doc.R failed
## [09:17:01] report.R done
your completed project should now look like this:
example-4
¦--boot
¦ ¦--data
¦ ¦ ¦--reportTemplate.docx
¦ ¦ ¦--sam_config
¦ ¦ ¦ °--model.cfg
¦ ¦ °--sam_data
¦ ¦ ¦--cn.dat
¦ ¦ ¦--cw.dat
¦ ¦ ¦--dw.dat
¦ ¦ ¦--lf.dat
¦ ¦ ¦--lw.dat
¦ ¦ ¦--mo.dat
¦ ¦ ¦--nm.dat
¦ ¦ ¦--pf.dat
¦ ¦ ¦--pm.dat
¦ ¦ ¦--survey.dat
¦ ¦ °--sw.dat
¦ ¦--DATA.bib
¦ ¦--initial
¦ ¦ °--data
¦ ¦--sam_config.R
¦ °--sam_data.R
¦--data
¦ ¦--catage.csv
¦ ¦--data.RData
¦ ¦--datage.csv
¦ ¦--landfrac.csv
¦ ¦--latage.csv
¦ ¦--natmort.csv
¦ ¦--propf.csv
¦ ¦--propm.csv
¦ ¦--survey_spring.csv
¦ ¦--survey_summer.csv
¦ ¦--wcatch.csv
¦ ¦--wdiscards.csv
¦ ¦--wlandings.csv
¦ °--wstock.csv
¦--data.R
¦--model
¦--model.R
¦--output
¦--output.R
¦--report
¦--report.R
¦--report.Rmd
¦--report_doc.R
¦--report_plots.R
¦--report_tables.R
°--utilities.R The final check that everything is working from start to finish is to run:
# re run boot, downloading all data again
taf.boot(taf = TRUE)
# clean all folders and rerun from scratch
sourceAll()