CMSY++ (Froese et al. 2023), the successor of CMSY (Froese et
al. 2017), estimates MSY, Bmsy and Fmsy, and the relative stock size
(B/Bmsy) and exploitation (F/Fmsy) of every year, from a catch series, a
prior for the intrinsic growth rate r and priors for the relative
biomass B/k. The package bundles the official CMSY++ code
(CMSY++16.R, January 2021), its neural network and its user
guide (Froese et al. 2021), and a few functions take care of everything
from installing JAGS (Plummer 2003) to drawing the charts. The CMSY++
runs shown in this article were made beforehand; building the article
does not run CMSY++.
FishStat has catches only, no CPUE. The FAO capture
statistics have no abundance index (CPUE or survey biomass), so every
stock is assessed with the catch-only part of CMSY++ (CMSY,
btype = "None") and the Bayesian Schaefer model (BSM) does
not run. The results then depend mostly on the B/k priors, and
independent tests found that catch-only methods classify stock status
poorly with default priors (Free et al. 2020; Ovando et al. 2021; Kell
et al. 2022). The results are indicative and do not
replace stock assessments. As in FAO (2026), the status has two classes
only: maximally sustainably fished (B/Bmsy at least 0.8
in the last year) or overfished (below 0.8); catch-only
results do not make a stock “underfished”.
Disclaimer and source of CMSY++. CMSY++ was
developed by Rainer Froese and colleagues (Froese et al. 2021, 2023).
faofishingreport is not affiliated with or endorsed by the authors of
CMSY++. The package bundles the authors’ unmodified March 2021 release
(GEOMAR OceanRep, record 52147, CC BY
4.0), whose code is also published on GitHub at https://github.com/SISTA16/cmsyPlusPlus.
cmsy_run() runs a copy of it with only the changes listed
in code_changes.csv of each run folder (operating-system
settings, no automatic package installation), which leave the method
unchanged. Questions on the method go to the authors of CMSY++; problems
with how this package runs it go to the faofishingreport
issues. Users are responsible for checking and interpreting the
results.
⚙️ Common arguments
These arguments have the same meaning in all CMSY functions that have them; the tables of each function below list its other arguments.
| Argument | Meaning | Default |
|---|---|---|
stocks |
Stocks to use: Stock identifiers
(e.g. "PRA_27") or species, by English or Chinese FAO name,
scientific name or FAO 3-alpha code (all stocks of the species) |
NULL: all stocks |
table |
Save the tables (in tables/cmsy/) |
TRUE |
lang |
Language of the chart text, names and table columns:
"en", "zh" or "both"
|
NULL: the language of fao_settings()
("en" unless changed) |
pptx |
Also save the charts as editable PowerPoint files (in
figures/pptx/cmsy/) |
FALSE; fao_settings(pptx = TRUE) changes
it for all functions |
plot_show |
Also draw each chart in the plot window (e.g. the Plots pane of RStudio) as it is saved |
FALSE, or plot_show of
fao_settings()
|
annotations |
Text added to the axis titles and labels: FALSE none,
TRUE all, or any of "title",
"subtitle", "caption", "legend"
and "labels" (text inside the chart) |
FALSE;
options(faofishingreport.cmsy_annotations = TRUE) changes
it for the session |
silent |
Print no progress bars and messages |
FALSE; fao_settings(silent = TRUE) changes
it for the session |
data_dir |
Folder where downloads are kept | fao_data_dir() |
🚀 One step: cmsy_area()
cmsy_setup() # once: JAGS and the R packages
area27 <- cmsy_area("Atlantic, Northeast", top = 10)cmsy_area() does Steps 4 to 7 and Step 9 below in one
call: it takes the top species with the largest catch of
the last recent_years years of a fishing area (or the
species given in species), prepares the input, runs CMSY++,
draws all charts, projects the stocks and writes the tables, with a
progress bar for each step. Main species that cannot be assessed
(e.g. without an r prior) are listed with the reason. area
can name several areas; each gets its own folder in
out_dir, sorted into:
cmsy_areas/Atlantic__Eastern_Central/
├── input/ CMSY++ catch and ID files (catch.csv, id.csv), the stocks assessed (stocks.csv)
│ and the main species left out with the reason (excluded.csv)
├── runs/ CMSY++ output of each stock (JPEG charts, cmsy_output.csv, log.txt), the changed
│ code and the list of changes (code_changes.csv)
├── figures/ a folder of charts per stock and summary/ (figures_pptx/ with pptx = TRUE)
└── tables/ results.csv, series.csv, priors_table.csv, results_table.csv, forecast.csv,
forecast_share.csv, settings_used.csv and sources.csv| Argument | Meaning | Default |
|---|---|---|
area |
FAO major fishing area name(s), as in
fishing_area_data()
(e.g. "Atlantic, Eastern Central"), or FAO catch data such
as fishing_area_data()$area_data
|
required |
top |
Number of species with the largest catch to assess | 10 |
species |
Species to assess instead of the top ones: English or
Chinese name, scientific name or FAO 3-alpha code |
NULL |
recent_years |
Species are ranked by their mean catch over the last
recent_years years of the data |
10 |
out_dir |
Folder in which a folder per area is made | "cmsy_areas" |
auto_start |
Start a stock in the year suggested by CMSY++ when its start biomass is ambiguous (see Step 5) |
TRUE (FALSE in
cmsy_run()) |
seed |
Random seed, so runs can be repeated with the same results | 1 |
parallel |
Number of stocks run at the same time; each uses two more processes for its two MCMC chains |
NULL: half the processor cores |
settings |
CMSY++ settings that replace its defaults (see The settings of CMSY++) |
list(): none |
charts |
Charts to draw (see Step 6) | all charts |
classes |
Status classes: "fao" (two) or "palomares"
(five) |
"fao" |
forecast |
Also project the stocks under the default scenarios of
cmsy_forecast()
|
TRUE |
annotations, lang, pptx,
silent
|
See Common arguments | |
... |
Further arguments of cmsy_input(),
e.g. priors, default_resilience or
min_years
|
cmsy_area() runs CMSY++ without the retrospective
analysis, so its retrospective chart is skipped (see Step 5). It returns
the result of cmsy_run() with input,
forecast, tables and dir (the
folder of the area); for several areas, a list of these named by area.
On a Mac with Apple silicon, the ten main species of area 27 take about
a minute from input to tables.
💻 Step 1: check the computer and adapt the settings
The CMSY++ code was written on Windows. Its user guide (step 12)
tells Mac users to change every windows() to
quartz() and antialias = "cleartype" to
"default". The code tests for Windows with
grepl("win", sysname), which is also true for macOS
(“Darwin”), so on a Mac it calls windows() and stops.
cmsy_system() finds the operating system, R version,
processor, JAGS and R packages, and shows the settings
cmsy_run() will use:
sys <- cmsy_system()
show_table(sys$changes)| change | applies_to | what |
|---|---|---|
| Package installation | all | Lines that install missing packages on every run are removed (cmsy_setup() installs them). |
| snpar | all | library(snpar) is removed and snpar::runs.test() is replaced by the same approximate runs test, reproduced from snpar 1.0 (Debin Qiu, GPL >= 2; archived on CRAN on 2022-05-10: https://cran.r-project.org/src/contrib/Archive/snpar/; source mirror: https://github.com/cran/snpar). |
| Working directory | all | setwd(rstudioapi::getActiveDocumentContext()$path) is removed, so RStudio is not needed; each stock runs in its own folder. |
| Settings | all | File names, the stock, the chart switches and the general settings (CMSY++ defaults unless changed with cmsy_run(settings = ), including the MCMC iterations of the JAGS call) are set by cmsy_run() after the original settings section. |
| Chart windows | linux | The test grepl(“win”, sysname), which is also true on macOS (“Darwin”), is replaced; windows(): changed to an off-screen device, so charts are drawn and saved without opening windows. |
| JPEG antialias | linux | antialias = “cleartype” (Windows only) is changed to “default” (user guide, step 12). |
| Results for charts | all | After the main analysis of a stock, its priors, the posterior draws of r, k and B/k, the yearly estimates with 95% limits and the last-year Kobe draws are saved to cmsy_extra.rds (retrospective runs to cmsy_retro.rds) for cmsy_plot(); no calculation is changed. |
The changes are made to a copy of the code written into each run
folder (CMSY++16_faofishingreport.R, with the list in
code_changes.csv); the bundled original never changes.
Charts are drawn on an off-screen device by default, which saves the
JPEG files on every system without opening windows;
device = "screen" uses the system’s own window instead
(windows(), quartz() or x11()).
If the detection is wrong, set the system:
cmsy_system(os = "windows") # settings used on Windows
options(faofishingreport.cmsy_os = "macos") # for the whole session| Argument | Meaning | Default |
|---|---|---|
os |
Operating system whose settings are used: "auto"
(detect it), "windows", "macos" or
"linux"
|
"auto" |
device |
Where CMSY++ draws its charts before saving them as JPEG files:
"file" (an off-screen device that works on every system and
in batch runs) or "screen" (the system’s own window) |
"file" |
jags |
JAGS version to recommend: "recommended" (the newest
JAGS 4 for this R version, which works with rjags from CRAN),
"latest" (JAGS 5.0.0, which needs rjags 5 from SourceForge)
or a version number such as "4.3.2"
|
"recommended" |
silent |
See Common arguments |
It returns the system (os, os_detected,
arch, r_version), the installed
jags and rjags versions, the
packages table, the recommended jags_installer
and rjags_source, the chart device and
antialias, and the changes made to the CMSY++
code.
🧰 Step 2: install JAGS and the R packages
cmsy_setup() asks before installing anything, then:
- installs the missing R packages from CRAN (R2jags, rjags, coda, gplots, mvtnorm, neuralnet, conicfit); snpar, archived on CRAN on 2022-05-10, is replaced by the same approximate runs test, reproduced from snpar 1.0 (Debin Qiu, GPL ≥ 2; CRAN archive, source mirror https://github.com/cran/snpar);
- downloads the JAGS installer for the system and R version from its
pinned link, checks its size and MD5 checksum and opens it
(
.exeon Windows,.pkgon macOS; on Linux it downloads the source and shows the commands). JAGS is a separate program: restart R after installing it and runcmsy_setup()again to check; - checks the bundled CMSY++ files.
| Argument | Meaning | Default |
|---|---|---|
jags |
JAGS version: "recommended", "latest" or a
version number (as in cmsy_system()) |
"recommended" |
os |
Operating system: "auto", "windows",
"macos" or "linux"
|
"auto" |
ask |
Ask before installing anything; in a non-interactive session nothing
is installed when ask is TRUE
|
TRUE |
data_dir, silent
|
See Common arguments |
It returns the result of cmsy_system() after the
installation. JAGS 4.3.2 (4.3.1 for older R versions) is the default
because rjags on CRAN (4-17) needs JAGS 4. The newest JAGS, 5.0.0, needs
rjags 5-3 from SourceForge: cmsy_setup(jags = "latest").
Installers copied into inst/jags/ of the package source are
used instead of downloading. Versions, links, sizes, checksums, licences
and references of every component are listed by
cmsy_sources():
src <- cmsy_sources()
show_table(src[, c("component", "version", "platform", "file", "license", "url")], n = 30)| Argument | Meaning | Default |
|---|---|---|
component |
Text to keep only matching components, e.g. "JAGS" or
"CMSY"
|
NULL: all |
The table also gives page, size,
md5, sha256, released,
requires and citation (the reference to cite)
of each component.
The R version and the JAGS version must match
JAGS is a separate program that rjags uses for the MCMC. On Windows and macOS each JAGS installer is built for a range of R versions; with the wrong one rjags cannot be loaded and CMSY++ cannot run. rjags on CRAN (4-17) needs JAGS 4. The matches (JAGS 4 is the default):
| System | R version | JAGS installer | rjags |
|---|---|---|---|
| Windows | 4.5, 4.6 | JAGS-4.3.2.exe |
4-17 (CRAN) |
| Windows | 4.2 to 4.4 | JAGS-4.3.1.exe |
4-17 (CRAN) |
| macOS | 4.3 and later | JAGS-4.3.2.pkg |
4-17 (CRAN) |
| macOS | 4.2 | JAGS-4.3.1.pkg |
4-15 or later (the CRAN binary for R 4.2 stays at 4-15) |
| Linux | any | the system package (e.g. sudo apt-get install jags) or
JAGS-4.3.2.tar.gz
|
4-17 (CRAN, built from source) |
| Windows | 4.5, 4.6 |
JAGS-5.0.0.exe (optional, newest) |
5-3 (SourceForge, not yet on CRAN) |
| macOS | 4.3 and later |
JAGS-5.0.0-universal.pkg (optional) |
5-3 (SourceForge: binary for R 4.6, else from source) |
The rules come from the README of the JAGS downloads (e.g. the
Windows build of JAGS 4.3.2 “was built with Rtools45 and is suitable for
interfacing with R 4.5.x”; “JAGS-4.3.2.pkg for compatibility with R
version 4.3.0 and later”; “JAGS-4.3.1.pkg for compatibility with R
versions 4.0.x to 4.2.x”), and each combination is installed and run
with CMSY++ by the package’s GitHub Actions workflow.
cmsy_jags_versions() gives the table, marks the row for
this computer and gives the evidence and README link of each rule:
jv <- cmsy_jags_versions()
show_table(jv[, c("system", "r_versions", "jags", "installer", "rjags", "default", "this_computer")], n = 10)| system | r_versions | jags | installer | rjags | default | this_computer |
|---|---|---|---|---|---|---|
| Windows | 4.5 to 4.6.x | 4.3.2 | JAGS-4.3.2.exe | 4-17 (CRAN) | TRUE | FALSE |
| Windows | 4.2 to 4.4.x | 4.3.1 | JAGS-4.3.1.exe | 4-17 (CRAN) | TRUE | FALSE |
| macOS | >= 4.3 | 4.3.2 | JAGS-4.3.2.pkg | 4-17 (CRAN) | TRUE | FALSE |
| macOS | 4.0 to 4.2.x | 4.3.1 | JAGS-4.3.1.pkg | 4-15 or later (CRAN; its binary for R 4.2 is 4-15) | TRUE | FALSE |
| Linux | >= 4.2 | 4.3.2 | JAGS-4.3.2.tar.gz | 4-17 (CRAN, built from source) | TRUE | TRUE |
| Windows | 4.5 to 4.6.x | 5.0.0 | JAGS-5.0.0.exe | 5-3 (SourceForge, Windows binary) | FALSE | FALSE |
| macOS | >= 4.3 | 5.0.0 | JAGS-5.0.0-universal.pkg | 5-3 (SourceForge: binary for R 4.6, else source) | FALSE | FALSE |
| Linux | >= 4.2 | 5.0.0 | JAGS-5.0.0.tar.gz | 5-3 (SourceForge, source) | FALSE | FALSE |
| Argument | Meaning | Default |
|---|---|---|
os |
"all", or one system: "auto" (this
computer), "windows", "macos" or
"linux"
|
"all" |
r_version |
R version marked as this_computer
|
getRversion(): the running R |
silent |
Do not print the match for this computer | see Common arguments |
No need to look it up: cmsy_download_jags() picks the
installer for the system, R version and processor, downloads it and
checks its size, MD5 and SHA-256; for JAGS 5 it also downloads rjags
5-3. The files are kept in fao_data_dir()/cmsy/downloads,
so they are downloaded only once. It can also download for other systems
or R versions, e.g. to prepare installers for computers without
internet. cmsy_install_jags() downloads the installer and
opens it:
cmsy_download_jags() # the JAGS for this computer
cmsy_download_jags(version = "latest") # JAGS 5.0.0 and rjags 5-3
cmsy_download_jags(os = "all", r_version = "4.4.1", dest = "jags_installers") # all systems, R 4.4
cmsy_install_jags() # download and open the installer| Argument | Meaning | Default |
|---|---|---|
version |
"recommended" (JAGS 4 for this R version),
"latest" (JAGS 5.0.0 with rjags 5-3) or a JAGS version such
as "4.3.1"
|
"recommended" |
os |
"auto" (this computer), "windows",
"macos", "linux" or "all" (the
installers of the three systems for r_version) |
"auto" |
r_version |
R version the installer is for |
getRversion(): the running R |
arch |
Processor, for the rjags 5 binary on macOS: "auto",
"arm64" (Apple silicon) or "x86_64"
(Intel) |
"auto" |
dest |
Folder to copy the files to |
NULL: none |
rjags |
Also download the rjags package needed by JAGS 5 (rjags 4-17 for
JAGS 4 comes from CRAN through cmsy_setup()) |
TRUE |
data_dir, silent
|
See Common arguments |
It returns a table with one row per file: system, R version, component, version, file, local path, link, SHA-256 and whether the checksums matched.
| Argument | Meaning | Default |
|---|---|---|
version |
"recommended", "latest" or a version
number |
"recommended" |
os |
"auto", "windows", "macos" or
"linux"
|
"auto" |
open |
Open the installer (.exe on Windows, .pkg
on macOS; on Linux the commands to install JAGS are shown) |
TRUE |
data_dir, silent
|
See Common arguments |
It returns the path of the downloaded installer; an installer placed
in inst/jags/ of the package source is used instead of
downloading.
🐟 Step 3: priors for r from FishBase and SeaLifeBase
A built-in table gives the resilience and prior r of the 3,155
species-level items of the FAO capture statistics in FishBase or
SeaLifeBase (release 26.06). It can be searched offline by English,
Chinese or scientific name or by 3-alpha code, and every row links to
its source page (open = TRUE opens the pages). How FishBase
and SeaLifeBase data are obtained and matched to FAO names is explained
in 🐠 FishBase and SeaLifeBase traits.
cod <- cmsy_priors("Gadus morhua")
show_table(cod[, c("species_name", "scientific_name", "SpecCode", "resilience", "prior_r",
"r_low", "r_high", "r_source", "url")])| species_name | scientific_name | SpecCode | resilience | prior_r | r_low | r_high | r_source | url |
|---|---|---|---|---|---|---|---|---|
| Atlantic cod | Gadus morhua | 69 | Medium | 0.51 | 0.34 | 0.76 | FishBase prior r, 95% CL (Based on 66 full stock assessments.) | https://www.fishbase.se/summary/69 |
The r range follows the CMSY++ user guide:
| Case | r range | Source |
|---|---|---|
| The species has a prior r in FishBase or SeaLifeBase (“Estimates based on models”) | The 95% confidence limits of the prior r | User guide Table 1 (“alternatively and better”) |
| Only the resilience is known | High 0.6–1.5; Medium 0.2–0.8; Low 0.05–0.5; Very low 0.015–0.1 | User guide Table 2 |
| Neither | No default: the user gives it (e.g. r close to 2M or 2Fmsy) | User guide Table 1 |
FAO names are matched as valid species names, then through each
database’s own synonym table; ambiguous synonyms and genus or family
items get no values. When FishBase lists several stocks of a species
with different resilience, the species-wide value is used and the others
are given in resilience_note.
priors <- cmsy_priors()
table(priors$match_type)
#>
#> ambiguous synonym exact not a species not found
#> 3 3071 2 43
#> synonym
#> 36
table(!is.na(priors$r_source))
#>
#> FALSE TRUE
#> 462 2693Species can be named by English or Chinese FAO name, both together
(as shown with lang = "both"), scientific name or FAO
3-alpha code, in cmsy_priors(),
cmsy_input(species = ), cmsy_run(stocks = )
and cmsy_plot(stocks = ): full names are matched first
(ignoring case), then parts of names; if several species match, they are
listed so a more precise name or code can be given. lang
sets the language of the results; species_en and
species_zh always give both names.
p <- cmsy_priors(c("COD", "Clupea harengus", "Atlantic mackerel", "小鳍鳕"), lang = "both")
show_table(p[, c("species_name", "scientific_name", "alpha3", "r_low", "r_high", "url", "data_url", "rule_url")])| species_name | scientific_name | alpha3 | r_low | r_high | url | data_url | rule_url |
|---|---|---|---|---|---|---|---|
| 大西洋鳕 Atlantic cod | Gadus morhua | COD | 0.34 | 0.76 | https://www.fishbase.se/summary/69 | https://data.source.coop/cboettig/fishbase/fb/v26.06/parquet/estimate.parquet | https://oceanrep.geomar.de/id/eprint/52147/ |
| 大西洋鲱 Atlantic herring | Clupea harengus | HER | 0.29 | 0.65 | https://www.fishbase.se/summary/24 | https://data.source.coop/cboettig/fishbase/fb/v26.06/parquet/estimate.parquet | https://oceanrep.geomar.de/id/eprint/52147/ |
| 大西洋鲭 Atlantic mackerel | Scomber scombrus | MAC | 0.28 | 0.64 | https://www.fishbase.se/summary/118 | https://data.source.coop/cboettig/fishbase/fb/v26.06/parquet/estimate.parquet | https://oceanrep.geomar.de/id/eprint/52147/ |
| 小鳍鳕 Blue whiting(=Poutassou) | Micromesistius poutassou | WHB | 0.31 | 0.70 | https://www.fishbase.se/summary/31 | https://data.source.coop/cboettig/fishbase/fb/v26.06/parquet/estimate.parquet | https://oceanrep.geomar.de/id/eprint/52147/ |
| Argument | Meaning | Default |
|---|---|---|
species |
Names to search for: English, Chinese or scientific names, FAO
3-alpha codes, or parts of them (e.g. "Atlantic cod",
"Gadus morhua", "COD",
"cod") |
NULL: all rows |
open |
Open the FishBase or SeaLifeBase pages of the rows found (at most 10) in the browser | FALSE |
lang |
Language of the species_name column: "en",
"zh" or "both"
|
NULL: the language of fao_settings()
|
Every row links to its sources: url, the FishBase or
SeaLifeBase page of the species (where resilience and prior r can be
checked); data_url, the FishBase or SeaLifeBase table the
value was read from (estimate for prior r,
stocks for resilience); and rule_url, the
CMSY++ user guide, whose Tables 1 and 2 turn the values into the r
range. The results of cmsy_input() and
cmsy_run() add fao_data_url, the FAO data
release used.
🧾 Step 4: FAO data to CMSY++ input files
input <- cmsy_input(FAO27$area_data)
input
#> CMSY++ input from FAO capture 2026.1.0: 220 stocks ready, 534 left out
#> Stock species_en area first_year last_year
#> HER_27 Atlantic herring Atlantic, Northeast 1950 2024
#> COD_27 Atlantic cod Atlantic, Northeast 1950 2024
#> CAP_27 Capelin Atlantic, Northeast 1950 2024
#> WHB_27 Blue whiting(=Poutassou) Atlantic, Northeast 1950 2024
#> MAC_27 Atlantic mackerel Atlantic, Northeast 1950 2024
#> SPR_27 European sprat Atlantic, Northeast 1950 2024
#> POK_27 Saithe(=Pollock) Atlantic, Northeast 1950 2024
#> HAD_27 Haddock Atlantic, Northeast 1950 2024
#> NOP_27 Norway pout Atlantic, Northeast 1959 2024
#> PIL_27 European pilchard(=Sardine) Atlantic, Northeast 1950 2024
#> Resilience r.low r.hi
#> Medium 0.288 0.654
#> Medium 0.336 0.764
#> Medium 0.412 0.937
#> Medium 0.307 0.698
#> Medium 0.283 0.644
#> Medium 0.318 0.722
#> Medium 0.344 0.781
#> Medium 0.328 0.746
#> Medium 0.329 0.747
#> Medium 0.415 0.942cmsy_input() makes the two input tables of CMSY++: the
catch file (Stock, yr, ct,
bt) and the ID file (one row per stock, with its priors and
settings). It also returns the table of stocks (stocks) and
of stocks left out with the reason (excluded), and with
table = TRUE saves the four tables in
tables/cmsy/ (<area>_cmsy_catch.csv,
..._id.csv, ..._stocks.csv and
..._excluded.csv). Following the user guide:
-
Stock: the catch of one species item in one fishing
area summed over countries (
by = "area");by = "country"keeps countries apart. Only items identified to species and measured in tonnes are used. -
Years: years are consecutive. A stock runs from its
first to its last year with catch, and missing or zero years in between
are interpolated (and counted in
filled_years). -
Left out: stocks with fewer than
min_yearsyears (CMSY++ uses the first and last ten years for its default priors), less thanmin_catchtonnes, more thanmax_filledyears filled, or no r prior.
show_table(input$stocks[, c("Stock", "species_name", "first_year", "last_year", "filled_years",
"Resilience", "r.low", "r.hi", "r_source")])| Stock | species_name | first_year | last_year | filled_years | Resilience | r.low | r.hi | r_source |
|---|---|---|---|---|---|---|---|---|
| HER_27 | Atlantic herring | 1950 | 2024 | 0 | Medium | 0.29 | 0.65 | FishBase prior r, 95% CL (Based on 53 full stock assessments.) |
| COD_27 | Atlantic cod | 1950 | 2024 | 0 | Medium | 0.34 | 0.76 | FishBase prior r, 95% CL (Based on 66 full stock assessments.) |
| CAP_27 | Capelin | 1950 | 2024 | 0 | Medium | 0.41 | 0.94 | FishBase prior r, 95% CL (Based on 5 full stock assessments.) |
| WHB_27 | Blue whiting(=Poutassou) | 1950 | 2024 | 0 | Medium | 0.31 | 0.70 | FishBase prior r, 95% CL (Based on 9 full stock assessments.) |
| MAC_27 | Atlantic mackerel | 1950 | 2024 | 0 | Medium | 0.28 | 0.64 | FishBase prior r, 95% CL (Based on 10 full stock assessments.) |
| SPR_27 | European sprat | 1950 | 2024 | 0 | Medium | 0.32 | 0.72 | FishBase prior r, 95% CL (Based on 17 full stock assessments.) |
| POK_27 | Saithe(=Pollock) | 1950 | 2024 | 0 | Medium | 0.34 | 0.78 | FishBase prior r, 95% CL (Based on 21 full stock assessments.) |
| HAD_27 | Haddock | 1950 | 2024 | 0 | Medium | 0.33 | 0.75 | FishBase prior r, 95% CL (Based on 41 full stock assessments.) |
| NOP_27 | Norway pout | 1959 | 2024 | 0 | Medium | 0.33 | 0.75 | FishBase prior r, 95% CL (Based on 3 full stock assessments.) |
| PIL_27 | European pilchard(=Sardine) | 1950 | 2024 | 0 | Medium | 0.42 | 0.94 | FishBase prior r, 95% CL (Based on 26 full stock assessments.) |
table(input$excluded$reason)
#>
#> fewer than 10 years
#> 10
#> more than 20% of years filled
#> 20
#> no resilience or r prior: give it in `priors` or set `default_resilience`
#> 22
#> total catch below 1000 t
#> 482Expert priors replace the defaults, e.g. a B/k prior for the last year as in Table 3 of the user guide, or a resilience for species without a FishBase value:
input <- cmsy_input(FAO27$area_data, priors = data.frame(
alpha3 = c("HER", "MAC"),
endb.low = c(0.2, 0.01),
endb.hi = c(0.6, 0.4)
))
input <- cmsy_input(FAO27$area_data, default_resilience = "Medium") # stocks without a FishBase value| Argument | Meaning | Default |
|---|---|---|
area_data |
FAO catch data: fishing_area_data()$area_data, or any
subset of load_fao_data()
|
required |
species |
Species to keep: English or Chinese FAO name, both together, scientific name or FAO 3-alpha code |
NULL: all species |
by |
"area": one stock per species item and fishing area,
summed over countries; "country": one stock per species
item, area and country |
"area" |
priors |
Data frame of priors and settings that replace the defaults, one row
per stock or species, identified by a Stock,
alpha3, species_en or species
column, with any of the ID file columns Resilience,
r.low, r.hi, stb.low,
stb.hi, int.yr, intb.low,
intb.hi, endb.low, endb.hi,
StartYear, EndYear and Comment
(recorded as “user”) |
NULL |
default_resilience |
Resilience ("High", "Medium",
"Low" or "Very low") for stocks whose species
has no resilience or prior r |
NULL: such stocks are left out, since the user guide
gives no default |
min_years |
Minimum number of years from the first to the last year with catch | 10 |
min_catch |
Minimum total catch of a stock (tonnes) | 1000 |
max_filled |
Maximum share of the years of a stock filled by interpolation | 0.2 |
label |
Name at the start of the table file names |
NULL: the fishing area (or country) of the data |
table, lang, silent
|
See Common arguments |
▶️ Step 5: run CMSY++
result <- cmsy_run(input, stocks = "Pandalus borealis", auto_start = TRUE, seed = 1,
plots = c("analysis", "management", "kobe", "prior_posterior", "retrospective"))
result$summaryEach stock runs in its own folder,
cmsy/<area>_cmsy/<stock>/
(e.g. cmsy/Atlantic__Northeast_cmsy/PRA_27/), and R
process; the folders keep CMSY++’s own output (_AN.jpg,
_MAN.jpg, _KOBE.jpg,
cmsy_output.csv, log.txt), and the folder of
the run has the changed code, code_changes.csv,
settings_used.csv and sources.csv.
parallel sets how many stocks run at the same time, and a
progress bar shows the stocks finished, the elapsed and the remaining
time (see “Progress bars and silent” below).
Many FAO series start in 1950 with low catches, which CMSY++ flags as
ambiguous (depleted or unexploited) and for which it suggests a start
year: the first year with 40% of the maximum catch.
auto_start = TRUE uses that year (CMSY++’s own
auto.start stops with an error in CMSY++16.R,
so the same rule is applied to the start year before the run, and not
used if fewer than ten years would remain). The user guide prefers
setting StartYear or the start B/k prior
(stb.low, stb.hi) yourself in
cmsy_input(priors = ).
| Argument | Meaning | Default |
|---|---|---|
input |
The result of cmsy_input()
|
required |
stocks |
Stocks to run (see Common arguments) |
NULL: all stocks of input
|
out_dir |
Folder for the CMSY++ output |
"cmsy" in the working directory |
parallel |
Number of stocks run at the same time; each uses two more processes for its two MCMC chains |
NULL: half the processor cores |
plots |
CMSY++ charts saved as JPEG files: any of "analysis",
"management", "kobe",
"prior_posterior" and "retrospective" (the
last two are needed by the charts of the same name in Step 6) |
c("analysis", "management", "kobe") |
os, device
|
Operating system and chart device (as in
cmsy_system()) |
"auto", "file"
|
seed |
Random seed, so runs can be repeated with the same results | NULL |
auto_start |
Start each stock in the year suggested by CMSY++ when its start biomass is ambiguous |
FALSE, as in CMSY++ |
settings |
CMSY++ settings that replace its defaults (see The settings of CMSY++) |
list(): none |
timeout |
Maximum time for one stock (seconds); a stock still running after it is reported as failed | 3600 |
table, lang, silent
|
See Common arguments; the tables are
tables/cmsy/<area>_cmsy_results.csv,
..._series.csv and ..._sources.csv
|
The example is the 30 stocks with the largest catch in 2015–2024 in
FAO area 27 (Northeast Atlantic), run with the start years suggested by
CMSY++ (p_green: probability of the green Kobe quadrant, B
≥ Bmsy and F ≤ Fmsy, in the last year):
top30 <- read.csv(file.path(fig_dir, "..", "top30_area27.csv"), check.names = FALSE)
top30$species <- fao_translate(top30$species_en, "species", getOption("faofishingreport.lang", "en"))
show_table(top30[, c("rank", "Stock", "species", "scientific_name", "first_year", "B_Bmsy", "F_Fmsy", "p_green", "status")],
n = 30, digits = 2)| rank | Stock | species | scientific_name | first_year | B_Bmsy | F_Fmsy | p_green | status |
|---|---|---|---|---|---|---|---|---|
| 1 | HER_27 | Atlantic herring | Clupea harengus | 1950 | 0.64 | 1.02 | 0.03 | overfished |
| 2 | WHB_27 | Blue whiting(=Poutassou) | Micromesistius poutassou | 1980 | 0.77 | 1.62 | 0.03 | overfished |
| 3 | COD_27 | Atlantic cod | Gadus morhua | 1950 | 0.67 | 0.83 | 0.03 | overfished |
| 4 | MAC_27 | Atlantic mackerel | Scomber scombrus | 1966 | 0.87 | 1.34 | 0.14 | maximally sustainably fished |
| 5 | SPR_27 | European sprat | Sprattus sprattus | 1972 | 0.66 | 1.02 | 0.03 | overfished |
| 6 | POK_27 | Saithe(=Pollock) | Pollachius virens | 1963 | 0.72 | 1.00 | 0.06 | overfished |
| 7 | CAP_27 | Capelin | Mallotus villosus | 1970 | 0.26 | 1.05 | 0.00 | overfished |
| 8 | HAD_27 | Haddock | Melanogrammus aeglefinus | 1955 | 0.71 | 1.09 | 0.05 | overfished |
| 9 | HKE_27 | European hake | Merluccius merluccius | 1950 | 0.73 | 0.84 | 0.09 | overfished |
| 10 | HOM_27 | Atlantic horse mackerel | Trachurus trachurus | 1971 | 0.27 | 0.99 | 0.00 | overfished |
| 11 | PIL_27 | European pilchard(=Sardine) | Sardina pilchardus | 1950 | 1.35 | 0.48 | 0.98 | maximally sustainably fished |
| 12 | REB_27 | Beaked redfish | Sebastes mentella | 2000 | 0.94 | 1.22 | 0.23 | maximally sustainably fished |
| 13 | PLE_27 | European plaice | Pleuronectes platessa | 1950 | 0.42 | 0.94 | 0.01 | overfished |
| 14 | PRA_27 | Northern prawn | Pandalus borealis | 1980 | 1.39 | 0.58 | 0.99 | maximally sustainably fished |
| 15 | SCE_27 | Great Atlantic scallop | Pecten maximus | 1973 | 1.05 | 1.11 | 0.36 | maximally sustainably fished |
| 16 | NOP_27 | Norway pout | Trisopterus esmarkii | 1968 | 0.20 | 0.51 | 0.00 | overfished |
| 17 | REG_27 | Golden redfish | Sebastes norvegicus | 1998 | 0.93 | 0.99 | 0.33 | maximally sustainably fished |
| 18 | GHL_27 | Greenland halibut | Reinhardtius hippoglossoides | 1965 | 1.28 | 0.60 | 0.94 | maximally sustainably fished |
| 19 | NEP_27 | Norway lobster | Nephrops norvegicus | 1966 | 0.98 | 0.92 | 0.42 | maximally sustainably fished |
| 20 | VMA_27 | Atlantic chub mackerel | Scomber colias | 2011 | 0.64 | 0.79 | 0.01 | overfished |
| 21 | CRE_27 | Edible crab | Cancer pagurus | 1970 | 0.81 | 1.06 | 0.17 | maximally sustainably fished |
| 22 | LIN_27 | Ling | Molva molva | 1950 | 0.92 | 0.82 | 0.36 | maximally sustainably fished |
| 23 | ANE_27 | European anchovy | Engraulis encrasicolus | 1950 | 0.64 | 1.03 | 0.01 | overfished |
| 24 | WHE_27 | Whelk | Buccinum undatum | 1996 | 0.92 | 0.94 | 0.32 | maximally sustainably fished |
| 25 | MUS_27 | Blue mussel | Mytilus edulis | 1978 | 0.29 | 0.89 | 0.00 | overfished |
| 26 | WHG_27 | Whiting | Merlangius merlangus | 1950 | 0.26 | 0.96 | 0.00 | overfished |
| 27 | CSH_27 | Common shrimp | Crangon crangon | 1950 | 0.57 | 0.73 | 0.02 | overfished |
| 28 | MON_27 | Angler(=Monk) | Lophius piscatorius | 1984 | 0.99 | 1.10 | 0.32 | maximally sustainably fished |
| 29 | ALB_27 | Albacore | Thunnus alalunga | 1950 | 0.67 | 0.93 | 0.02 | overfished |
| 30 | COC_27 | Common edible cockle | Cerastoderma edule | 1981 | 0.34 | 1.14 | 0.00 | overfished |
ex <- read.csv(file.path(fig_dir, "..", "example_values.csv"), check.names = FALSE)
green <- ex[ex$example == "green", ]
over <- ex[ex$example == "overfished", ]Two stocks are shown in detail below: the healthiest and the most clearly overfished.
-
Healthy: northern prawn (Pandalus
borealis,
PRA_27, 14th by catch): last-year B/Bmsy 1.39 (95% CI 1.08–1.63), F/Fmsy 0.58, a 99% probability of the green quadrant. -
Overfished: capelin (Mallotus villosus,
CAP_27, 7th by catch): last-year B/Bmsy 0.26 (95% CI 0.10–0.54), F/Fmsy 1.05, a 100% probability that B/Bmsy is below 0.8 and 52% of the red quadrant.
Capelin is a short-lived fish whose stocks rise and fall with recruitment, and fisheries on it have been closed in some years; a catch-only model reads such falls in catch as falls in biomass. These are completed CMSY++ runs with the B/k priors of the neural network, not validated stock assessments; compare them with the assessments of ICES before using them.
Their priors and results tables, from cmsy_tables()
(Step 9), northern prawn first:
| Stock | Species | Scientific name | Years | r prior | Source of the r prior | k prior (t) | MSY prior (t) | B/k prior, first year | B/k prior, intermediate year | B/k prior, last year | FishBase / SeaLifeBase page |
|---|---|---|---|---|---|---|---|---|---|---|---|
| PRA_27 | Northern prawn | Pandalus borealis | 1980-2024 | 0.408-0.927 | SeaLifeBase prior r, 95% CL (Based on 7 full stock assessments.) | 354,482-1,504,879 | 112,898 | 0.202-0.608 (neural network) | 2013: 0.0857-0.361 (neural network) | 0.466-0.837 (neural network) | https://www.sealifebase.se/summary/25892 |
| Stock | Species | r (95% CI) | k (t, 95% CI) | MSY (t, 95% CI) | Fmsy (95% CI) | Bmsy (t, 95% CI) | Last-year B/Bmsy (95% CI) | Last-year F/Fmsy (95% CI) | P(green quadrant) | P(red quadrant) | Status (FAO) | Status (Palomares et al. 2018) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| PRA_27 | Northern prawn | 0.502 (0.356-0.615) | 985,962 (786,871-1,437,661) | 123,757 (110,512-143,108) | 0.251 (0.178-0.308) | 492,981 (393,435-718,830) | 1.39 (1.08-1.63) | 0.579 (0.395-0.866) | 99% | 0% | Maximally sustainably fished | Healthy |
| Stock | Species | Scientific name | Years | r prior | Source of the r prior | k prior (t) | MSY prior (t) | B/k prior, first year | B/k prior, intermediate year | B/k prior, last year | FishBase / SeaLifeBase page |
|---|---|---|---|---|---|---|---|---|---|---|---|
| CAP_27 | Capelin | Mallotus villosus | 1970-2024 | 0.412-0.937 | FishBase prior r, 95% CL (Based on 5 full stock assessments.) | 7,217,484-30,658,265 | 2,322,843 | 0.156-0.509 (neural network) | 2019: 0.0239-0.23 (neural network) | 0.0372-0.258 (neural network) | https://www.fishbase.se/summary/252 |
| Stock | Species | r (95% CI) | k (t, 95% CI) | MSY (t, 95% CI) | Fmsy (95% CI) | Bmsy (t, 95% CI) | Last-year B/Bmsy (95% CI) | Last-year F/Fmsy (95% CI) | P(green quadrant) | P(red quadrant) | Status (FAO) | Status (Palomares et al. 2018) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CAP_27 | Capelin | 0.423 (0.31-0.517) | 23,255,085 (18,563,960-35,085,235) | 2,456,458 (2,070,092-3,026,716) | 0.211 (0.155-0.258) | 11,627,543 (9,281,980-17,542,618) | 0.255 (0.0975-0.541) | 1.05 (0.238-7.44) | 0% | 52% | Overfished | Grossly overfished |
result$summary has one row per stock: MSY, Bmsy, Fmsy, r
and k with their 95% limits, the last-year B/Bmsy, F/Fmsy and biomass
with their limits, the probability of each Kobe quadrant and of B/Bmsy
below 0.8 (p_overfished), the status in two schemes
(status, FAO; status_palomares, Palomares et
al. 2018), all priors with their sources and CMSY++’s own notes (such as
suggested_start_year). result$series gives
catch, predicted catch, B/k, biomass, B/Bmsy, F and F/Fmsy by year with
95% limits; result$kobe the last-year draws;
result$details the posterior draws of r, k and B/k, the
prior parameters and the retrospective runs;
result$settings the CMSY++ settings used;
result$changes the changes to the code; and
result$sources the versions and links of CMSY++, JAGS, R,
the R packages, the FAO data and FishBase used.
📊 Step 6: charts
cmsy_plot(result)Based on a survey of the charts of 54 CMSY application papers,
cmsy_plot() draws, for each stock, every chart of the
CMSY++ user guide, and the usual summaries of several stocks, in
figures/png/cmsy/<stock>/ and
figures/png/cmsy/summary/ (and as PowerPoint with
pptx = TRUE). It returns the ggplot objects.
By default only axis titles and labels are shown, so
charts can take report captions of their own. Titles, subtitles,
captions (method, data and prior sources), legends and text in the
charts (quadrant probabilities, PPVR, Mohn’s rho) are switched on with
annotations: TRUE for all, or
e.g. c("title", "legend");
options(faofishingreport.cmsy_annotations = TRUE) changes
the session default.
| Chart | Content |
|---|---|
analysis |
The six panels of the CMSY++ analysis chart (user guide Figure 1): catch and MSY, explored r-k pairs with the prior box, viable r-k pairs, B/k with its three priors, F/Fmsy, equilibrium curve |
catch, rk, biomass,
b_bmsy, exploitation,
equilibrium
|
The same panels as single charts, and B/Bmsy (lines at Bmsy, 0.8 Bmsy, 0.5 Bmsy) |
management |
The four management panels (Figure 2): catch and MSY, B/Bmsy, F/Fmsy, Kobe plot |
kobe |
Kobe plot (Figure 3): path from the first (square) through the intermediate (circle) to the last year (triangle), 50/80/95% regions of the last-year draws, quadrant probabilities |
prior_posterior |
Prior and posterior of r, k, MSY and B/k in three years, with PPVR (Figure 4) |
retrospective |
B/Bmsy and F/Fmsy without the last 1–3 years, with Mohn’s rho (Mohn
1999; Figure 5; needs
cmsy_run(plots = "retrospective")) |
overview, kobe_all
|
Last-year B/Bmsy of each stock with limits; Kobe plot of all stocks (sized by catch) |
status, status_time,
histogram
|
Share of stocks and of catch by status, in the last year and every year; distribution of last-year B/Bmsy |
b_bmsy_all, trend,
aggregate
|
Small multiples of B/Bmsy; median B/Bmsy and share below 0.8 Bmsy; total catch / total MSY and total biomass / total Bmsy |
map |
Median last-year B/Bmsy by FAO fishing area |
| Argument | Meaning | Default |
|---|---|---|
result |
The result of cmsy_run()
|
required |
stocks |
Stocks to draw (see Common arguments) |
NULL: all stocks that ran |
charts |
Charts to draw (table above) | all charts |
classes |
Status classes: "fao" (two classes) or
"palomares" (five classes) |
"fao" |
annotations, lang, pptx,
silent
|
See Common arguments |
Status follows FAO (2026) for Tier 3 (data-limited) assessments by
default: maximally sustainably fished (B/Bmsy ≥ 0.8) or overfished.
classes = "palomares" uses the five classes of Palomares et
al. (2018), common in CMSY studies: healthy (≥ 1), slightly overfished
(0.8–1), overfished (0.5–0.8), grossly overfished (0.2–0.5), collapsed
(< 0.2).
The CMSY++ analysis chart (no annotations, the default), northern prawn and capelin:
Management charts:
Kobe plots with annotations = TRUE:
Prior and posterior distributions; the smaller the PPVR, the more the data added to the prior:
Retrospective analyses:
Summaries of the 30 stocks of area 27
(annotations = "legend"):
🔁 MCMC: settings and convergence
mcmc <- cmsy_mcmc(result)
mcmc$settings # chains, iterations, burn-in, thinning and draws kept of each stock
mcmc$diagnostics # posterior summary, Rhat and effective sample size of each parameter
cmsy_mcmc(result, colors = c("grey30", "#D7191C")) # other colours for the chainsCMSY++ samples the posterior with JAGS by Markov chain Monte Carlo
(MCMC). cmsy_run() keeps the draws of each chain of r, k,
B/k, B/Bmsy and F/Fmsy (MSY is computed from r and k, as in
CMSY++16.R), and cmsy_mcmc() checks them.
Rhat (Gelman and Rubin 1992) compares the variance
between chains with the variance within them; at most 1.1 the chains are
taken as converged (Gelman et al. 2013). The effective sample
size (n_eff, computed with coda; Plummer et
al. 2006) is the number of independent draws that the autocorrelated
draws are worth; below 400 (Vehtari et al. 2021) the 95% limits are
imprecise. Stocks that did not converge and stocks with too few
effective draws are reported apart. The chart shows the draws of each
chain in the order they were kept (trace: chains that mix well overlap
with no trend) and their densities (which should agree).
The two example stocks, with the settings of CMSY++16.R
(2 chains of 60,000 iterations, the first 30,000 discarded as burn-in,
one draw in 10 kept) and with 4 chains of 300,000 iterations (150,000
burn-in, one draw in 50 kept; columns _long):
mc <- rbind(cbind(stock = "PRA_27", read.csv(file.path(fig_dir, "green_mcmc.csv"), check.names = FALSE)),
cbind(stock = "CAP_27", read.csv(file.path(fig_dir, "overfished_mcmc.csv"), check.names = FALSE)))
mc <- mc[mc$parameter != "deviance", ]
rk <- mc[mc$parameter %in% c("r", "k", "msy"), ]
shown <- mc[, c("stock", "label", "median", "rhat", "n_eff", "median_long", "rhat_long", "n_eff_long")]
shown[c("median", "median_long")] <- lapply(shown[c("median", "median_long")], function(x)
vapply(signif(x, 3), format, character(1), big.mark = ",", scientific = FALSE))
shown[c("rhat", "rhat_long")] <- lapply(shown[c("rhat", "rhat_long")], sprintf, fmt = "%.3f")
shown[c("n_eff", "n_eff_long")] <- lapply(shown[c("n_eff", "n_eff_long")], function(x) format(round(x), big.mark = ","))
show_table(data.frame(shown, row.names = NULL), n = 20, digits = 3)| stock | label | median | rhat | n_eff | median_long | rhat_long | n_eff_long |
|---|---|---|---|---|---|---|---|
| PRA_27 | Intrinsic growth rate r | 0.459 | 1.003 | 183 | 0.468 | 1.001 | 1,471 |
| PRA_27 | Carrying capacity k (t) | 1,070,000 | 1.002 | 164 | 1,050,000 | 1.001 | 1,254 |
| PRA_27 | MSY (t) | 124,000 | 1.000 | 151 | 123,000 | 1.001 | 1,451 |
| PRA_27 | B/k in the first year | 0.424 | 1.003 | 262 | 0.432 | 1.000 | 1,929 |
| PRA_27 | B/k in the last year | 0.697 | 1.008 | 230 | 0.693 | 1.000 | 1,528 |
| PRA_27 | B/Bmsy in the last year | 1.39 | 1.008 | 230 | 1.39 | 1.000 | 1,528 |
| PRA_27 | F/Fmsy in the last year | 0.579 | 1.001 | 492 | 0.588 | 1.000 | 2,811 |
| CAP_27 | Intrinsic growth rate r | 0.386 | 1.055 | 256 | 0.383 | 1.001 | 2,105 |
| CAP_27 | Carrying capacity k (t) | 25,400,000 | 1.094 | 124 | 25,800,000 | 1.001 | 703 |
| CAP_27 | MSY (t) | 2,460,000 | 1.034 | 80 | 2,480,000 | 1.000 | 524 |
| CAP_27 | B/k in the first year | 0.331 | 1.034 | 216 | 0.326 | 1.000 | 867 |
| CAP_27 | B/k in the last year | 0.128 | 1.018 | 142 | 0.13 | 1.007 | 1,437 |
| CAP_27 | B/Bmsy in the last year | 0.255 | 1.018 | 142 | 0.259 | 1.007 | 1,437 |
| CAP_27 | F/Fmsy in the last year | 1.05 | 1.082 | 262 | 1.03 | 1.022 | 1,926 |
With the default settings the chains converged (Rhat at most 1.09), but the effective sample size of r, k and MSY was only 80 to 256: successive draws of r and k are strongly correlated. With the longer chains it rose to 524 to 2105, while the medians hardly changed (last-year B/Bmsy of northern prawn 1.395 and 1.386, of capelin 0.255 and 0.259), and the runs took about five times as long (4 chains instead of 2 run side by side, but each five times longer). For results to report, run CMSY++ with longer and more chains and check them:
result <- cmsy_run(input, seed = 1, auto_start = TRUE,
settings = list(n.chains = 4, n.iter = 300000, n.burnin = 150000, n.thin = 50))
cmsy_mcmc(result)Trace and density of each chain, northern prawn and capelin (default settings):
| Argument | Meaning | Default |
|---|---|---|
result |
The result of cmsy_run() or
cmsy_area()
|
required |
stocks |
Stocks to check (see Common arguments) |
NULL: all stocks that ran |
rhat_max |
Largest Rhat taken as converged | 1.1 |
n_eff_min |
Smallest effective sample size (all chains together) taken as enough | 400 |
plot, table
|
Save the chart figures/png/cmsy/<stock>/mcmc.png
and the tables tables/cmsy/mcmc_settings.csv and
mcmc_diagnostics.csv
|
TRUE |
colors |
Colours of the chains: NULL (Okabe-Ito colours), a
palette name (e.g. "npg", "Set1",
"viridis") or one colour per chain |
NULL |
lang, pptx, silent
|
See Common arguments |
The same chart is the "mcmc" chart of
cmsy_plot() (colours with mcmc_colors).
Results of earlier versions of the package have no draws: run CMSY++
again to check them. Converged chains do not mean that the priors or the
model fit the stock.
🔮 Step 7: projections under fishing scenarios
projection <- cmsy_forecast(result, catch = 100000)CMSY++16.R has no projection, so
cmsy_forecast() projects with the CMSY++ model itself
(Schaefer surplus production with productivity declining below 0.25 k),
from the posterior draws of r, k and last-year B/k of the same MCMC
iterations, with lognormal process error of standard deviation
sigma (0.1 by default, the value CMSY++16.R
gives to its setting sigmaR, which the model itself does
not use: it estimates its process error). As in CMSY++, the biomass of a
year is the biomass before that year’s catch, F is the catch divided by
that biomass, and F/Fmsy uses an Fmsy reduced linearly below half of
Bmsy. The last year of the data keeps its catch; from the first
projected year on, the catch of an F scenario is the scenario F times
the biomass of the same year, so F/Fmsy is the scenario value. The
default scenarios are 0.5, 0.6, 0.8 and 0.95 Fmsy (as in Froese et
al. 2018) and the last-year fishing mortality; constant catches can be
added (catch). Stocks whose data end in different years are
each projected from their own last year to the same final year.
| Argument | Meaning | Default |
|---|---|---|
result |
The result of cmsy_run()
|
required |
stocks |
Stocks to project (see Common arguments) |
NULL: all stocks that ran |
years |
Number of years to project after the latest last year of the stocks | 10 |
f |
Fishing mortality scenarios as multiples of Fmsy | c(0.5, 0.6, 0.8, 0.95) |
status_quo |
Also project the last-year fishing mortality | TRUE |
catch |
Constant catches to project (tonnes per year) | NULL |
sigma |
Standard deviation of the lognormal process error | 0.1 |
seed |
Random seed | NULL |
plot |
Save the charts
figures/png/cmsy/<stock>/forecast.png and
figures/png/cmsy/summary/forecast_share.png
|
TRUE |
table |
Save tables/cmsy/forecast.csv and
forecast_share.csv
|
TRUE |
lang, annotations, pptx,
silent
|
See Common arguments |
It returns projection (by stock, scenario and year:
median B/Bmsy, F/Fmsy and catch with 95% limits, and the probability
that B is at least Bmsy) and share (by scenario and year:
share of stocks at or above Bmsy, with 95% limits, and the number of
stocks n_stocks; from the latest last year of the stocks
on).
With constant catches near the MSY of each stock (100,000 t for northern prawn, 2,000,000 t for capelin) and the F scenarios:
🔍 Step 8: compare runs or other assessments
result_1950 <- cmsy_run(input, stocks = "PRA_27", seed = 1) # starting in 1950
cmsy_compare(`1950` = result_1950, `CMSY++` = result, stocks = "PRA_27")cmsy_compare() draws several runs (e.g. different start
years or priors, a sensitivity analysis) and other estimates
(external, e.g. B/Bmsy from ICES, JABBA or a paper)
together and lists their last-year values. The start year matters.
Started in 1950, northern prawn has a last-year B/Bmsy of 0.82, below
Bmsy; started in 1980 as CMSY++ suggests, 1.39, above it. Capelin is
overfished from either start year (0.25 from 1950, 0.26 from 1970):
| Argument | Meaning | Default |
|---|---|---|
... |
Results of cmsy_run(), named or not
(e.g. default = run1, start_1980 = run2) |
|
stocks |
Stocks to draw (see Common arguments) |
NULL: all stocks that ran |
external |
Data frame of other estimates, with columns Stock,
Year, B_Bmsy (optionally F_Fmsy)
and source (the name shown in the legend) |
NULL |
lang, annotations, pptx,
silent
|
See Common arguments |
The charts are saved as
figures/png/cmsy/compare/<stock>.png; the function
returns last_year (one row per stock and run) and
plots.
📋 Step 9: tables for reports
tables <- cmsy_tables(result)cmsy_tables() makes the two tables of CMSY papers, the
priors table and the results table shown in Step 5, rounded, with ranges
as “low-high” and column names in the chosen language, saved as
tables/cmsy/cmsy_priors_table.csv and
cmsy_results_table.csv. The results give the status in the
two classes of FAO (2026) and in the five classes of Palomares et
al. (2018).
| Argument | Meaning | Default |
|---|---|---|
result |
The result of cmsy_run()
|
required |
stocks |
Stocks in the tables (see Common arguments) |
NULL: all stocks that ran |
table, lang, silent
|
See Common arguments |
⚙️ The settings of CMSY++
CMSY++ has a set of general settings, such as the catch observation
error CV.C and the number of MCMC chains and iterations.
Two of them, sigmaR and ni, are set in
CMSY++16.R but not used by its code (the model estimates
its process error with the prior itau2 ~ dgamma(4, 0.01)):
they are listed as “not used” and cannot be changed, so that a run is
not mistaken for a sensitivity analysis of the process error.
cmsy_run() keeps the CMSY++ values unless they are changed;
cmsy_defaults() lists each one with its CMSY++ default, the
value used and what it does:
show_table(cmsy_defaults()[, c("setting", "cmsy_default", "applies_to", "description")], n = 30)| setting | cmsy_default | applies_to | description |
|---|---|---|---|
| CV.C | 0.15 | all | Coefficient of variation of the catch (observation error) |
| CV.cpue | 0.2 | BSM | Minimum realistic coefficient of variation of CPUE |
| sigmaR | 0.1 | not used | Not used by CMSY++16.R: the code sets it (sigR <- sigmaR) but never uses it; the model estimates its process error with the prior itau2 ~ dgamma(4, 0.01). It cannot be changed |
| cor.log.rk | -0.76 | charts | Correlation of log r and log k, used only in charts |
| rk.cor.beta | 2.52, 3.37 | all | Beta prior of the r-k correlation plus 1 |
| nbk | 3 | BSM | Number of B/k priors used by BSM: 1 (first year), 2 (and intermediate year), 3 (and last year) |
| bt4pr | FALSE | BSM | Use the abundance data to set the B/k priors |
| ct_MSY.lim | 1.21 | all | Ratio of catch to the MSY prior above which the B/k prior is kept constant |
| q.biomass.pr | 0.9, 1.1 | BSM | Prior range of q when the abundance data are biomass |
| n | 5000 | charts | Number of points of the r-k cloud in the analysis chart |
| ni | 3 | not used | Not used by CMSY++16.R: the code sets it but never uses it. It cannot be changed |
| nab | 3 | BSM | Minimum number of years with abundance data to run BSM |
| bw | 3 | all | Bandwidth of ksmooth() for smoothing the catch |
| e.creep.line | TRUE | charts | Show the CPUE not corrected for effort creep in the biomass chart |
| select.yr | NA | output | Year for which B, F, B/Bmsy and F/Fmsy are also reported (NA: none) |
| n.chains | 2 | all | Number of MCMC chains (JAGS) |
| n.iter | 60000 | all | Number of MCMC iterations per chain (JAGS) |
| n.burnin | 30000 | all | Number of burn-in iterations (JAGS) |
| n.thin | 10 | all | Thinning of the MCMC chains (JAGS) |
Settings marked “BSM” only matter with abundance data (CPUE or
biomass); with FAO catch data only those marked “all” affect the
results. To change some, give only those in settings; the
others keep their CMSY++ values, and every value is checked (name,
range, enough MCMC draws):
cmsy_defaults(list(CV.C = 0.2, n.iter = 100000, n.burnin = 50000))[c(1, 17, 18), c("setting", "cmsy_default", "used", "changed")]
#> setting cmsy_default used changed
#> 1 CV.C 0.15 0.2 TRUE
#> 17 n.iter 60000 100000 TRUE
#> 18 n.burnin 30000 50000 TRUE
result <- cmsy_run(input, stocks = "PRA_27", seed = 1, auto_start = TRUE,
settings = list(CV.C = 0.2, n.iter = 100000, n.burnin = 50000))
result$settings[result$settings$changed, ] # the changed settings| Argument | Meaning | Default |
|---|---|---|
settings |
Named list of changed settings, e.g. list(CV.C = 0.2),
shown next to the defaults and checked |
list() |
lang |
Language of the description column: "en",
"zh" or "both"
|
NULL: the language of fao_settings()
|
The settings used are saved in settings_used.csv of the
run folder and recorded in result$settings and
result$sources. The chart and output switches of CMSY++
(mgraphs, kobe.plot, retros, …)
are set by the plots argument of cmsy_run(),
not in settings.
🖥️ Windows, macOS and Linux
The same code runs on all three systems: cmsy_system()
detects the system and makes the changes it needs (on macOS, no
windows() and no antialias = "cleartype"; see
Step 1), and the JAGS installer is chosen for the system and R version.
Each change to the CMSY code is checked by the GitHub Actions workflow
“CMSY-cross-platform” on Windows, macOS and Linux: it downloads and
installs the matching JAGS with cmsy_download_jags(), runs
CMSY++ for real, runs all CMSY tests and keeps the output of each
system.
Results of the check in October 2026, made on two stocks of area 34,
the example area of this article at the time (all 170 tests pass for
each combination; same seed; cmsy_area() on bigeye grunt
GRB_34 and sardine PIL_34):
| System | R | JAGS | rjags | GRB_34 B/Bmsy | GRB_34 F/Fmsy | PIL_34 B/Bmsy | PIL_34 F/Fmsy |
|---|---|---|---|---|---|---|---|
| Windows | 4.6.1 | 4.3.2 | 4-17 | 1.422 | 0.639 | 0.706 | 1.391 |
| Windows | 4.4.3 | 4.3.1 | 4-17 | 1.422 | 0.639 | 0.706 | 1.391 |
| Windows | 4.3.3 | 4.3.1 | 4-17 | 1.422 | 0.639 | 0.706 | 1.391 |
| Windows | 4.2.3 | 4.3.1 | 4-17 | 1.422 | 0.639 | 0.706 | 1.391 |
| macOS (Apple silicon) | 4.6.1 | 4.3.2 | 4-17 | 1.422 | 0.639 | 0.706 | 1.391 |
| Ubuntu | 4.6.1 | 4.3.2 (apt) | 4-17 | 1.422 | 0.639 | 0.706 | 1.391 |
| macOS (Apple silicon) | 4.2.3 | 4.3.1 | 4-15 | 1.366 | 0.663 | 0.745 | 1.306 |
| Windows | 4.6.1 | 5.0.0 | 5-3 | 1.392 | 0.672 | 0.740 | 1.317 |
| macOS (Apple silicon) | 4.6.1 | 5.0.0 | 5-3 | 1.392 | 0.672 | 0.740 | 1.317 |
With the same JAGS 4.3.x and rjags 4-17, the results are identical on
all three systems. JAGS 5 or an older rjags draws different MCMC
samples, so the numbers differ slightly (B/Bmsy of GRB_34 from 1.37 to
1.42, all inside its 95% interval 1.08–1.65), but the status of both
stocks is the same: GRB_34 maximally sustainably fished, PIL_34
overfished. To repeat a result exactly, use the same JAGS and rjags
versions; they are recorded in result$sources and
tables/sources.csv.
⏱️ Progress bars and silent
cmsy_run(), cmsy_plot(),
cmsy_forecast() and cmsy_area() show progress
bars such as
CMSY++ [================== ] 18/30 ( 60%) 00:51 elapsed, 00:34 left COD_27All CMSY functions that print have a silent argument:
with silent = TRUE they print no progress bars or messages
(cmsy_setup() still asks before installing).
fao_settings(silent = TRUE) sets the session default.
📝 In the Word report
fao_report() assesses the main species of the area with
CMSY++ when cmsy = TRUE (see 📝
Assessment report: Word, LaTeX and PDF), and adds a part on the
assessment with the results of every species and appendix tables of the
reference points, the priors and their sources. It runs
cmsy_area() on the top_n main species of the
report (10 by default), without projections. The settings of
cmsy_area() and cmsy_input() are given with
the prefix cmsy_, e.g. cmsy_top = 5,
cmsy_default_resilience = "Medium" or
cmsy_settings = list(n.iter = 100000). All CMSY++ results
(runs, charts and tables) are kept in the folder cmsy/ next
to the report, with cmsy_result.rds, so the report can be
written again without running CMSY++:
cmsy = readRDS("cmsy/<area>/cmsy_result.rds"), with
the path from the working directory (by default
reports/cmsy/<area>/cmsy_result.rds). The result of
an earlier cmsy_area() run can be given in the same
way.
report <- fao_report(27, cmsy = TRUE, cmsy_top = 5)
report$cmsy_reference_points # r, k, MSY, B/Bmsy, ... with 95% limits
fao_report(27, cmsy = readRDS("reports/cmsy/Atlantic__Northeast/cmsy_result.rds"))📖 References
Boettiger, C., Lang, D.T. and Wainwright, P.C. (2012). rfishbase: exploring, manipulating and visualizing FishBase data from R. Journal of Fish Biology 81(6): 2030–2039. doi:10.1111/j.1095-8649.2012.03464.x. Source code: https://github.com/ropensci/rfishbase
FAO (2026). Global Capture Production 1950–2024, release 2026.1.0. In: Fisheries and Aquaculture. Rome: FAO. https://www.fao.org/fishery/collection/capture/en. Licence: CC BY 4.0. Release file: https://www.fao.org/fishery/static/Data/Capture_2026.1.0.zip
FAO (2026). The State of World Fisheries and Aquaculture 2026 – Blue Transformation: turning vision into impact. Rome: FAO. doi:10.4060/cd8357en. Section “The status of fishery resources”, Box 1.3 (assessment tiers and status classes): https://openknowledge.fao.org/server/api/core/bitstreams/09d03e2f-7097-4d59-bf94-db65e95f44fd/content/sofia-2026/the-status-fishery-resources.html
Free, C.M., Jensen, O.P., Anderson, S.C., Gutierrez, N.L., Kleisner, K.M., Longo, C., Minto, C., Osio, G.C. and Walsh, J.C. (2020). Blood from a stone: performance of catch-only methods in estimating stock biomass status. Fisheries Research 223: 105452. doi:10.1016/j.fishres.2019.105452
Froese, R., Demirel, N., Coro, G., Kleisner, K.M. and Winker, H. (2017). Estimating fisheries reference points from catch and resilience. Fish and Fisheries 18: 506–526. doi:10.1111/faf.12190
Froese, R., Demirel, N., Coro, G. and Winker, H. (2021). User
Guide for CMSY++. GEOMAR, Germany, 17 pp. https://oceanrep.geomar.de/id/eprint/52147/ (bundled
with the package:
system.file("cmsy", "CMSYUserGuideMarch2021.pdf", package = "faofishingreport"))
Froese, R. and Pauly, D. (Eds.) (2026). FishBase. World Wide Web electronic publication. Version (06/2026). https://www.fishbase.se
Froese, R., Winker, H., Coro, G., Demirel, N., Tsikliras, A.C., Dimarchopoulou, D., Scarcella, G., Quaas, M. and Matz-Lück, N. (2018). Status and rebuilding of European fisheries. Marine Policy 93: 159–170. doi:10.1016/j.marpol.2018.04.018
Froese, R., Winker, H., Coro, G., Palomares, M.L.D., Tsikliras, A.C., Dimarchopoulou, D., Touloumis, K., Demirel, N., Vianna, G.M.S., Scarcella, G., Schijns, R., Liang, C. and Pauly, D. (2023). New developments in the analysis of catch time series as the basis for fish stock assessments: the CMSY++ method. Acta Ichthyologica et Piscatoria 53: 173–189. doi:10.3897/aiep.53.105910
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B. (2013). Bayesian Data Analysis, 3rd edition. Boca Raton: CRC Press. doi:10.1201/b16018
Gelman, A. and Rubin, D.B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science 7(4): 457–472. doi:10.1214/ss/1177011136
Kell, L.T., Sharma, R. and Winker, H. (2022). Artefact and artifice: evaluation of the skill of catch-only methods for classifying stock status. Frontiers in Marine Science 9: 762203. doi:10.3389/fmars.2022.762203
Mohn, R. (1999). The retrospective problem in sequential population analysis: an investigation using cod fishery and simulated data. ICES Journal of Marine Science 56: 473–488. doi:10.1006/jmsc.1999.0481
Ovando, D., Hilborn, R., Monnahan, C., Rudd, M., Sharma, R., Thorson, J.T., Rousseau, Y. and Ye, Y. (2021). Improving estimates of the state of global fisheries depends on better data. Fish and Fisheries 22: 1377–1391. doi:10.1111/faf.12593
Palomares, M.L.D., Froese, R., Derrick, B., Nöel, S.-L., Tsui, G., Woroniak, J. and Pauly, D. (2018). A preliminary global assessment of the status of exploited marine fish and invertebrate populations. A report prepared by the Sea Around Us for OCEANA. University of British Columbia and GEOMAR. https://www.fishbase.de/rfroese/OceanaFinalReport_V7a.pdf
Palomares, M.L.D. and Pauly, D. (Eds.) (2026). SeaLifeBase. World Wide Web electronic publication. Version (06/2026). https://www.sealifebase.se
Plummer, M. (2003). JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. In: Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003), Vienna, 1–10. https://www.r-project.org/conferences/DSC-2003/Proceedings/Plummer.pdf; JAGS: https://mcmc-jags.sourceforge.io/
Plummer, M., Best, N., Cowles, K. and Vines, K. (2006). CODA: convergence diagnosis and output analysis for MCMC. R News 6(1): 7–11. https://journal.r-project.org/articles/RN-2006-002/
Plummer, M. (2025). rjags: Bayesian Graphical Models using MCMC. R package version 4-17. doi:10.32614/CRAN.package.rjags
Su, Y. and Yajima, M. (2024). R2jags: Using R to Run ‘JAGS’. R package version 0.8-9. doi:10.32614/CRAN.package.R2jags
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved Rhat for assessing convergence of MCMC. Bayesian Analysis 16(2): 667–718. doi:10.1214/20-BA1221