📉 CMSY++仅渔获量资源评估 • faofishingreport

fao_settings(lang = "zh")   # 本文所有图表使用中文

CMSY++（Froese et al. 2023）是CMSY（Froese et al. 2017）的后续版本，它根据渔获量序列、内禀增长率r的先验和相对生物量B/k的先验，估计MSY、Bmsy、Fmsy，以及每年的相对生物量B/Bmsy和相对捕捞强度F/Fmsy。本包内置了CMSY++的官方代码（CMSY++16.R，2021年1月版）、神经网络文件和用户手册（Froese et al. 2021），由几个函数完成从安装JAGS（Plummer 2003）到出图的全部步骤。本文展示的CMSY++结果是事先运行得到的，生成本文时不会运行CMSY++。

FishStat只有渔获量，没有CPUE。FAO捕捞产量数据不含CPUE、调查生物量等丰度指数，所以每个种群都只能走CMSY++的仅渔获量路线（CMSY，btype = "None"），不会运行贝叶斯Schaefer模型（BSM）。这时结果主要取决于B/k先验。独立检验发现，在默认先验下，仅渔获量方法对种群状态的判断精度有限（Free et al. 2020；Ovando et al. 2021；Kell et al. 2022）。因此结果只是指示性的，不能替代正式的种群评估。按照FAO（2026）的做法，状态只分两类：最大可持续捕捞（最后一年B/Bmsy ≥ 0.8）和过度捕捞（B/Bmsy < 0.8），不使用仅渔获量结果判定“低度捕捞”。

CMSY++免责说明与来源。CMSY++由Rainer Froese等人开发（Froese et al. 2021, 2023）。本包与CMSY++的作者无关联，也未得到其认可。本包内置的是作者发布的2021年3月版原始文件（GEOMAR OceanRep，记录52147，CC BY 4.0），同样的代码也发布在GitHub：https://github.com/SISTA16/cmsyPlusPlus。cmsy_run()运行的是它的一份副本，所做的修改都列在每个运行文件夹的code_changes.csv中（操作系统设置、不自动安装R包），不改变方法本身。方法上的问题请咨询CMSY++的作者；本包运行CMSY++时出现的问题请提交到faofishingreport的Issues。评估结果需由使用者自行核对和解释。

⚙️ 通用参数

下面这些参数在所有带有它们的CMSY函数中含义相同；后面各函数的参数表只列出其他参数。

参数	含义	默认值
`stocks`	要使用的种群：种群代码（如`"PRA_27"`），或物种的英文名、中文名、学名或FAO三字母代码（该物种的所有种群）	`NULL`：全部种群
`table`	是否保存表格（在`tables/cmsy/`）	`TRUE`
`lang`	图中文字、名称和表格列名的语言：`"en"`、`"zh"`或`"both"`	`NULL`：使用`fao_settings()`设定的语言（未修改时为`"en"`）
`pptx`	是否同时把图保存为可编辑的PowerPoint文件（在`figures/pptx/cmsy/`）	`FALSE`；`fao_settings(pptx = TRUE)`可对所有函数生效
`plot_show`	保存图表的同时把图画在绘图窗口中（例如RStudio的Plots窗格）	`FALSE`，或`fao_settings()`中的`plot_show`
`annotations`	在坐标轴标题和刻度之外加入的文字：`FALSE`不加，`TRUE`全部加入，或`"title"`、`"subtitle"`、`"caption"`、`"legend"`和`"labels"`（图内文字）中的任意几项	`FALSE`；`options(faofishingreport.cmsy_annotations = TRUE)`可改变本次会话的默认值
`silent`	是否不显示进度条和提示信息	`FALSE`；`fao_settings(silent = TRUE)`可改变本次会话的默认值
`data_dir`	存放下载文件的文件夹	`fao_data_dir()`

🚀 一步完成：cmsy_area()

cmsy_setup()                                                     # 第一次使用前：安装JAGS和R包
area27 <- cmsy_area("Atlantic, Northeast", top = 10)

cmsy_area()把下文第四步到第七步以及第九步合成一步：选出某个渔区近recent_years年渔获量最大的top个物种（也可以用species直接指定），依次准备输入、运行CMSY++、画全部图表、做情景预测、写表格，每一步都有进度条。无法评估的主要物种（如没有r先验）会列出原因。area可以是多个渔区，每个渔区的结果保存在out_dir下各自的文件夹中，并分类存放：

cmsy_areas/Atlantic__Eastern_Central/
├── input/      CMSY++渔获文件和ID文件（catch.csv、id.csv）、所评估的种群（stocks.csv），
│               以及未能评估的主要物种及原因（excluded.csv）
├── runs/       每个种群的CMSY++原始输出（JPEG图、cmsy_output.csv、log.txt）、修改后的代码
│               和修改清单（code_changes.csv）
├── figures/    每个种群一个文件夹，另有summary/汇总图（pptx = TRUE时还有figures_pptx/）
└── tables/     results.csv、series.csv、priors_table.csv、results_table.csv、forecast.csv、
                forecast_share.csv、settings_used.csv和sources.csv

参数	含义	默认值
`area`	FAO主要渔区名称（可以是多个），写法与`fishing_area_data()`相同（如`"Atlantic, Eastern Central"`），也可以是`fishing_area_data()$area_data`这样的FAO产量数据	必填
`top`	评估渔获量最大的几个物种	`10`
`species`	改为评估指定的物种：英文名、中文名、学名或FAO三字母代码	`NULL`
`recent_years`	按数据最后`recent_years`年的平均渔获量给物种排序	`10`
`out_dir`	在其中为每个渔区建一个文件夹	`"cmsy_areas"`
`auto_start`	起始生物量难以判断时，是否从CMSY++建议的年份开始（见第五步）	`TRUE`（`cmsy_run()`中为`FALSE`）
`seed`	随机数种子，使结果可以重复	`1`
`parallel`	同时运行的种群数；每个种群的两条MCMC链还会另外占用两个进程	`NULL`：处理器核数的一半
`settings`	替换CMSY++默认值的设置（见CMSY++自带的设置）	`list()`：不修改
`charts`	要画的图（见第六步）	全部
`classes`	状态分级：`"fao"`（两级）或`"palomares"`（五级）	`"fao"`
`forecast`	是否按`cmsy_forecast()`的默认情景做预测	`TRUE`
`annotations`、`lang`、`pptx`、`silent`	见通用参数
`...`	传给`cmsy_input()`的其他参数，如`priors`、`default_resilience`或`min_years`

cmsy_area()运行CMSY++时不做回顾性分析，所以不画回顾性分析图（见第五步）。函数返回cmsy_run()的结果，另加input、forecast、tables和dir（该渔区的文件夹）；多个渔区时返回以渔区命名的列表。在Apple芯片的Mac上，27区前10个物种从输入到表格共用时约1分钟。

💻 第一步：检查电脑并自动调整设置

CMSY++的代码是在Windows上写的。用户手册第12步要求Mac用户把所有windows()改成quartz()，把antialias = "cleartype"改成"default"。代码用grepl("win", sysname)判断是不是Windows，而macOS的系统名”Darwin”中也含有”win”，所以在Mac上会调用windows()并报错。cmsy_system()检测当前系统、R版本、处理器、JAGS和所需的R包，并列出cmsy_run()将在这个系统上使用的设置：

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.

这些修改只作用于每次运行时写入运行目录的代码副本（CMSY++16_faofishingreport.R，修改清单见同目录下的code_changes.csv），包内的原始文件始终不变。默认使用离屏绘图设备，在任何系统上都能直接保存JPEG图，不会弹出窗口；device = "screen"则改用系统自己的绘图窗口（Windows上是windows()，macOS上是quartz()，Linux上是x11()）。如果自动检测有误，可以手动指定系统：

cmsy_system(os = "windows")                       # 查看Windows下的设置
options(faofishingreport.cmsy_os = "macos")       # 本次会话内都按macOS设置

参数	含义	默认值
`os`	按哪个操作系统设置：`"auto"`（自动检测）、`"windows"`、`"macos"`或`"linux"`	`"auto"`
`device`	CMSY++保存JPEG图之前在哪里绘图：`"file"`（离屏设备，任何系统和批处理都可用）或`"screen"`（系统自己的窗口）	`"file"`
`jags`	推荐的JAGS版本：`"recommended"`（适用于当前R版本的最新JAGS 4，可配合CRAN上的rjags）、`"latest"`（JAGS 5.0.0，需要SourceForge上的rjags 5）或版本号，如`"4.3.2"`	`"recommended"`
`silent`	见通用参数

函数返回系统信息（os、os_detected、arch、r_version）、已安装的jags和rjags版本、packages表、推荐的jags_installer和rjags_source、绘图设置device和antialias，以及对CMSY++代码所做的修改changes。

🧰 第二步：安装JAGS和R包

cmsy_setup()

cmsy_setup()在安装任何东西之前都会先询问，然后依次：

从CRAN安装缺少的R包：R2jags、rjags、coda、gplots、mvtnorm、neuralnet、conicfit。原代码依赖的snpar已于2022-05-10从CRAN下架，本包用一段相同的近似游程检验代码代替了它，该代码复现自snpar 1.0（Debin Qiu，GPL ≥ 2；CRAN存档，源代码镜像https://github.com/cran/snpar）；
按当前系统和R版本，从固定链接下载JAGS安装程序，核对文件大小和MD5后打开安装程序：Windows上是.exe，macOS上是.pkg；Linux上下载源码，并显示安装命令。JAGS是独立程序，安装完成后需要重启R，再运行一次cmsy_setup()检查；
核对包内CMSY++文件。

参数	含义	默认值
`jags`	JAGS版本：`"recommended"`、`"latest"`或版本号（与`cmsy_system()`相同）	`"recommended"`
`os`	操作系统：`"auto"`、`"windows"`、`"macos"`或`"linux"`	`"auto"`
`ask`	安装前是否询问；`ask`为`TRUE`时，非交互会话中不会安装任何东西	`TRUE`
`data_dir`、`silent`	见通用参数

函数返回安装后的cmsy_system()结果。默认安装JAGS 4.3.2（较早的R版本安装4.3.1），因为CRAN上的rjags（4-17）只支持JAGS 4。最新的JAGS 5.0.0需要使用SourceForge上的rjags 5-3，可以这样安装：cmsy_setup(jags = "latest")。如果把安装程序放进源码包的inst/jags/文件夹，安装时会直接使用这些文件，不再联网下载。每个组件的版本、下载链接、文件大小、校验值、许可证和引用文献，都可以用cmsy_sources()查看：

src <- cmsy_sources()
show_table(src[, c("component", "version", "platform", "file", "license", "url")], n = 30)

component	version	platform	file	license	url
CMSY++ code	16 (January 2021)	all	CMSY++16.R	CC BY 4.0	https://oceanrep.geomar.de/id/eprint/52147/4/UserGuideCodeDataMarch2021.zip
CMSY++ neural network	16 (January 2021)	all	ffnn.bin	CC BY 4.0	https://oceanrep.geomar.de/id/eprint/52147/4/UserGuideCodeDataMarch2021.zip
CMSY++ user guide	March 2021	all	CMSYUserGuideMarch2021.pdf	CC BY 4.0	https://oceanrep.geomar.de/id/eprint/52147/4/UserGuideCodeDataMarch2021.zip
CMSY++ training data (catch)	9e	all	Train_Catch_9e.csv	CC BY 4.0	https://oceanrep.geomar.de/id/eprint/52147/4/UserGuideCodeDataMarch2021.zip
CMSY++ training data (stocks)	9e	all	Train_ID_9e.csv	CC BY 4.0	https://oceanrep.geomar.de/id/eprint/52147/4/UserGuideCodeDataMarch2021.zip
JAGS	4.3.2	windows	JAGS-4.3.2.exe	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/4.x/Windows/JAGS-4.3.2.exe
JAGS	4.3.1	windows	JAGS-4.3.1.exe	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/4.x/Windows/JAGS-4.3.1.exe
JAGS	4.3.2	macos	JAGS-4.3.2.pkg	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/4.x/Mac%20OS%20X/JAGS-4.3.2.pkg
JAGS	4.3.1	macos	JAGS-4.3.1.pkg	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/4.x/Mac%20OS%20X/JAGS-4.3.1.pkg
JAGS	4.3.2	source	JAGS-4.3.2.tar.gz	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/4.x/Source/JAGS-4.3.2.tar.gz
JAGS	5.0.0	windows	JAGS-5.0.0.exe	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/5.x/Windows/Rtools45/JAGS-5.0.0.exe
JAGS	5.0.0	macos	JAGS-5.0.0-universal.pkg	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/5.x/macOS/JAGS-5.0.0-universal.pkg
JAGS	5.0.0	source	JAGS-5.0.0.tar.gz	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/JAGS/5.x/Source/JAGS-5.0.0.tar.gz
rjags	4-17	all	rjags	GPL-2	https://cran.r-project.org/package=rjags
rjags	5-3	windows	rjags_5-3.zip	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/rjags/5/Windows/Rtools45/rjags_5-3.zip
rjags	5-3	macos-arm64	rjags_5-3-R4.6.1-arm64.tgz	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/rjags/5/macOS/rjags_5-3-R4.6.1-arm64.tgz
rjags	5-3	macos-x86_64	rjags_5-3-R4.6.1-x86_64.tgz	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/rjags/5/macOS/rjags_5-3-R4.6.1-x86_64.tgz
rjags	5-3	source	rjags_5-3.tar.gz	GPL-2	https://downloads.sourceforge.net/project/mcmc-jags/rjags/5/Source/rjags_5-3.tar.gz
snpar (runs test reproduced in cmsy_run())	1.0	all	snpar_1.0.tar.gz	GPL (>= 2)	https://cran.r-project.org/src/contrib/Archive/snpar/snpar_1.0.tar.gz
R2jags	0.8-9	all	R2jags	GPL (> 2)	https://cran.r-project.org/package=R2jags
coda	0.19-4.1	all	coda	GPL (>= 2)	https://cran.r-project.org/package=coda
gplots	3.3.0	all	gplots	GPL-2	https://cran.r-project.org/package=gplots
mvtnorm	1.4-2	all	mvtnorm	GPL-2	https://cran.r-project.org/package=mvtnorm
neuralnet	1.44.2	all	neuralnet	GPL (>= 2)	https://cran.r-project.org/package=neuralnet
conicfit	1.0.4	all	conicfit	GPL (>= 3)	https://cran.r-project.org/package=conicfit
FishBase	26.06	all	species, synonyms, stocks, estimate	CC BY-NC 4.0 (FishBase.org team, as stated for the snapshots)	https://data.source.coop/cboettig/fishbase/fb/v26.06/parquet/
SeaLifeBase	26.06	all	species, synonyms, stocks, estimate	CC BY-NC 4.0 (FishBase.org team, as stated for the snapshots)	https://data.source.coop/cboettig/fishbase/slb/v26.06/parquet/

参数	含义	默认值
`component`	只保留名称含有这段文字的组件，如`"JAGS"`或`"CMSY"`	`NULL`：全部

表中还有每个组件的page、size、md5、sha256、released、requires和citation（引用文献）。

R版本与JAGS版本必须匹配

JAGS是独立的程序，rjags包通过它运行MCMC。Windows和macOS上的JAGS安装程序都是针对一定范围的R版本编译的，装错版本时rjags会加载失败，CMSY++也就无法运行。另外，CRAN上的rjags（4-17）只支持JAGS 4。对应关系如下（默认使用JAGS 4）：

系统	R版本	JAGS安装文件	rjags
Windows	4.5、4.6	`JAGS-4.3.2.exe`	4-17（CRAN）
Windows	4.2–4.4	`JAGS-4.3.1.exe`	4-17（CRAN）
macOS	4.3及以后	`JAGS-4.3.2.pkg`	4-17（CRAN）
macOS	4.2	`JAGS-4.3.1.pkg`	4-15或以后（CRAN上R 4.2的二进制包停在4-15）
Linux	任意	系统软件包（如`sudo apt-get install jags`）或`JAGS-4.3.2.tar.gz`源码	4-17（CRAN，从源码编译）
Windows	4.5、4.6	`JAGS-5.0.0.exe`（可选，最新版）	5-3（SourceForge，尚未上CRAN）
macOS	4.3及以后	`JAGS-5.0.0-universal.pkg`（可选）	5-3（SourceForge：R 4.6有二进制包，其他版本需从源码编译）

依据是JAGS下载页面的官方说明，例如Windows版JAGS 4.3.2用Rtools45编译，适用于R 4.5；macOS版JAGS-4.3.2.pkg适用于R 4.3.0及以后，JAGS-4.3.1.pkg适用于R 4.0–4.2。此外，本包在GitHub Actions上对表中每一种组合都实际安装并运行了CMSY++。cmsy_jags_versions()给出这张表，标出本机对应的那一行，并附上每条规则的依据和官方说明页面链接：

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

参数	含义	默认值
`os`	`"all"`（全部），或某一个系统：`"auto"`（本机）、`"windows"`、`"macos"`或`"linux"`	`"all"`
`r_version`	标为`this_computer`的R版本	`getRversion()`：当前运行的R
`silent`	是否不显示本机的匹配结果	见通用参数

不需要自己查表：cmsy_download_jags()根据本机的系统、R版本和芯片自动挑出匹配的安装文件，下载后核对文件大小、MD5和SHA-256；选JAGS 5时会同时下载配套的rjags 5-3。文件保存在fao_data_dir()/cmsy/downloads，只需下载一次。也可以替其他系统或其他R版本下载，例如给实验室的离线电脑准备安装包。cmsy_install_jags()则下载并打开安装程序：

cmsy_download_jags()                                        # 本机对应的JAGS
cmsy_download_jags(version = "latest")                      # JAGS 5.0.0和rjags 5-3
cmsy_download_jags(os = "all", r_version = "4.4.1", dest = "jags_installers")   # 三个系统、R 4.4
cmsy_install_jags()                                         # 下载并打开安装程序

cmsy_download_jags()：

参数	含义	默认值
`version`	`"recommended"`（适用于当前R版本的JAGS 4）、`"latest"`（JAGS 5.0.0和rjags 5-3）或JAGS版本号，如`"4.3.1"`	`"recommended"`
`os`	`"auto"`（本机）、`"windows"`、`"macos"`、`"linux"`或`"all"`（三个系统中适用于`r_version`的安装文件）	`"auto"`
`r_version`	安装文件对应的R版本	`getRversion()`：当前运行的R
`arch`	处理器，用于选择macOS上的rjags 5二进制包：`"auto"`、`"arm64"`（Apple芯片）或`"x86_64"`（Intel）	`"auto"`
`dest`	把文件复制到哪个文件夹	`NULL`：不复制
`rjags`	是否同时下载JAGS 5所需的rjags包（JAGS 4所用的rjags 4-17由`cmsy_setup()`从CRAN安装）	`TRUE`
`data_dir`、`silent`	见通用参数

函数返回一张表，每个文件一行：系统、R版本、组件、版本、文件、本地路径、链接、SHA-256以及校验是否通过。

cmsy_install_jags()：

参数	含义	默认值
`version`	`"recommended"`、`"latest"`或版本号	`"recommended"`
`os`	`"auto"`、`"windows"`、`"macos"`或`"linux"`	`"auto"`
`open`	是否打开安装程序（Windows上是`.exe`，macOS上是`.pkg`；Linux上显示安装JAGS的命令）	`TRUE`
`data_dir`、`silent`	见通用参数

函数返回下载的安装程序路径；源码包inst/jags/中放有安装程序时，直接使用该文件，不再下载。

🐟 第三步：r先验（FishBase和SeaLifeBase）

包里内置了一张表，收录FAO捕捞统计中所有种一级的物种（3,155个）在FishBase或SeaLifeBase（26.06版）中的恢复力（resilience）和先验r，可以离线搜索。搜索支持英文名、中文名、学名或FAO三字母代码，每一行都附有来源页面的链接（open = TRUE会在浏览器中打开这些页面）。FishBase和SeaLifeBase数据怎样获取、怎样与FAO名称匹配，见🐠 FishBase与SeaLifeBase物种特征。

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
大西洋鳕	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

r的先验范围按照CMSY++用户手册确定：

情况	r范围	依据
物种在FishBase或SeaLifeBase有先验r（“Estimates based on models”）	先验r的95%置信限	用户手册表1（“更好的做法”）
只有恢复力	高0.6–1.5；中0.2–0.8；低0.05–0.5；很低0.015–0.1	用户手册表2
两者都没有	无默认值，需要由用户提供（如r ≈ 2M或r ≈ 2Fmsy）	用户手册表1

FAO学名先按有效种名匹配，再按数据库自己的同物异名表匹配。有歧义的异名，以及属、科一级的条目，都不赋值。FishBase对同一物种列出恢复力不同的几个种群时，采用全种的值，其他值写在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  2693

物种可以用英文名、中文名、中英文名（即lang = "both"时显示的形式）、学名或FAO三字母代码给出。cmsy_priors()、cmsy_input(species = )、cmsy_run(stocks = )和cmsy_plot(stocks = )都使用同一套规则：先按完整名称匹配（不区分大小写），没有完全匹配时再按名称的一部分匹配；如果匹配到多个物种，会列出候选项，请从中选择更具体的名称或代码。结果用哪种语言显示由lang决定，species_en和species_zh两列则始终同时给出英文名和中文名。

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/

参数	含义	默认值
`species`	要搜索的名称：英文名、中文名、学名、FAO三字母代码，或其中的一部分（如`"Atlantic cod"`、`"Gadus morhua"`、`"COD"`、`"cod"`）	`NULL`：全部
`open`	是否在浏览器中打开找到的物种（最多10个）的FishBase或SeaLifeBase页面	`FALSE`
`lang`	`species_name`列的语言：`"en"`、`"zh"`或`"both"`	`NULL`：使用`fao_settings()`设定的语言

每一行都附有三个来源链接：url是FishBase或SeaLifeBase的物种页面，可以在页面上核对恢复力和先验r；data_url是读取这个数值所用的FishBase或SeaLifeBase数据表（先验r来自estimate表，恢复力来自stocks表）；rule_url是CMSY++用户手册，根据手册表1、表2把这些数值换算成r的范围。cmsy_input()和cmsy_run()的结果还多一列fao_data_url，即所用FAO数据版本的下载链接。

🧾 第四步：把FAO数据转成CMSY++输入文件

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.942

cmsy_input()生成CMSY++需要的两张输入表：渔获文件（Stock、yr、ct、bt）和ID文件（每个种群一行，包含先验和设置）。它同时返回种群表（stocks）和被排除种群及原因（excluded）；table = TRUE时把这四张表保存在tables/cmsy/（<渔区>_cmsy_catch.csv、..._id.csv、..._stocks.csv和..._excluded.csv）。处理规则按照用户手册：

种群：默认为某个物种在某个渔区内所有国家渔获的合计（by = "area"）；by = "country"则按国家分开。只使用鉴定到种、以吨计的条目。
年份：年份必须连续。序列从第一个有渔获的年份开始，到最后一个有渔获的年份结束；中间缺报或为零的年份用插值补齐，并记录补齐的年数（filled_years）。
排除：少于min_years年（CMSY++的默认先验要用到序列前10年和后10年）、总渔获低于min_catch吨、补齐年份超过max_filled，或没有任何r先验的种群。

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	Resilience	r.low	r.hi	r_source
HER_27	大西洋鲱	1950	2024	Medium	0.29	0.65	FishBase prior r, 95% CL (Based on 53 full stock assessments.)
COD_27	大西洋鳕	1950	2024	Medium	0.34	0.76	FishBase prior r, 95% CL (Based on 66 full stock assessments.)
CAP_27	毛鳞鱼	1950	2024	Medium	0.41	0.94	FishBase prior r, 95% CL (Based on 5 full stock assessments.)
WHB_27	小鳍鳕	1950	2024	Medium	0.31	0.70	FishBase prior r, 95% CL (Based on 9 full stock assessments.)
MAC_27	大西洋鲭	1950	2024	Medium	0.28	0.64	FishBase prior r, 95% CL (Based on 10 full stock assessments.)
SPR_27	黍鲱	1950	2024	Medium	0.32	0.72	FishBase prior r, 95% CL (Based on 17 full stock assessments.)
POK_27	青鳕 Saithe(=Pollock)	1950	2024	Medium	0.34	0.78	FishBase prior r, 95% CL (Based on 21 full stock assessments.)
HAD_27	黑线鳕	1950	2024	Medium	0.33	0.75	FishBase prior r, 95% CL (Based on 41 full stock assessments.)
NOP_27	挪威长臀鳕	1959	2024	Medium	0.33	0.75	FishBase prior r, 95% CL (Based on 3 full stock assessments.)
PIL_27	沙丁鱼	1950	2024	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 
#>                                                                       482

如果有专家判断，可以用priors替换默认值。例如按用户手册表3设置最后一年的B/k先验，或给没有FishBase值的种群指定恢复力：

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")  # 没有FishBase值的种群统一用"中"

参数	含义	默认值
`area_data`	FAO产量数据：`fishing_area_data()$area_data`，或`load_fao_data()`的任意子集	必填
`species`	只保留这些物种：英文名、中文名、中英文名、学名或FAO三字母代码	`NULL`：全部物种
`by`	`"area"`：每个物种条目和渔区一个种群，各国合计；`"country"`：每个物种条目、渔区和国家一个种群	`"area"`
`priors`	替换默认值的先验和设置，数据框，每个种群或物种一行，用`Stock`、`alpha3`、`species_en`或`species`列识别，可包含ID文件的`Resilience`、`r.low`、`r.hi`、`stb.low`、`stb.hi`、`int.yr`、`intb.low`、`intb.hi`、`endb.low`、`endb.hi`、`StartYear`、`EndYear`和`Comment`列（来源记为”user”）	`NULL`
`default_resilience`	物种在FishBase和SeaLifeBase中既无恢复力也无先验r时使用的恢复力（`"High"`、`"Medium"`、`"Low"`或`"Very low"`）	`NULL`：排除这些种群，因为用户手册没有给出默认值
`min_years`	从首次到最后一次有渔获的最少年数	`10`
`min_catch`	种群的最低总渔获量（吨）	`1000`
`max_filled`	插值补齐的年份占全部年份的最高比例	`0.2`
`label`	表格文件名开头的名称	`NULL`：数据所属的渔区（或国家）
`table`、`lang`、`silent`	见通用参数

▶️ 第五步：运行CMSY++

result <- cmsy_run(input, stocks = "Pandalus borealis", auto_start = TRUE, seed = 1,
                   plots = c("analysis", "management", "kobe", "prior_posterior", "retrospective"))
result$summary

每个种群在cmsy/<渔区>_cmsy/<种群>/（如cmsy/Atlantic__Northeast_cmsy/PRA_27/）下有独立的文件夹，并在单独的R进程中运行。文件夹里保留了CMSY++的全部原始输出：分析图_AN.jpg、管理图_MAN.jpg、Kobe图_KOBE.jpg、cmsy_output.csv和运行日志log.txt；运行目录中还有修改后的代码、code_changes.csv、settings_used.csv和sources.csv。parallel控制同时运行几个种群；运行时有进度条，显示已完成的种群数、已用时间和预计剩余时间（见下文“进度条与silent”）。

很多FAO序列从1950年开始，那时的渔获量很低。对此CMSY++会提示：起始生物量可能偏低，也可能尚未开发，难以判断。它建议改用渔获量首次达到最大值40%的年份作为起始年。auto_start = TRUE就是采用这个建议的起始年（CMSY++自带的auto.start选项在CMSY++16.R中会报错，本包按同一规则在运行前改写起始年；如果改写后不足10年，则保留原起始年）。用户手册更推荐由用户在cmsy_input(priors = )中自己设定StartYear或首年B/k先验（stb.low、stb.hi）。

参数	含义	默认值
`input`	`cmsy_input()`的结果	必填
`stocks`	要运行的种群（见通用参数）	`NULL`：`input`中的全部种群
`out_dir`	存放CMSY++输出的文件夹	工作目录下的`"cmsy"`
`parallel`	同时运行的种群数；每个种群的两条MCMC链还会另外占用两个进程	`NULL`：处理器核数的一半
`plots`	保存为JPEG的CMSY++图：`"analysis"`、`"management"`、`"kobe"`、`"prior_posterior"`和`"retrospective"`中的任意几项（第六步中的同名图需要后两项）	`c("analysis", "management", "kobe")`
`os`、`device`	操作系统和绘图设备（与`cmsy_system()`相同）	`"auto"`、`"file"`
`seed`	随机数种子，使结果可以重复	`NULL`
`auto_start`	起始生物量难以判断时，是否从CMSY++建议的年份开始	`FALSE`，与CMSY++相同
`settings`	替换CMSY++默认值的设置（见CMSY++自带的设置）	`list()`：不修改
`timeout`	单个种群的最长运行时间（秒），超时的种群记为失败	`3600`
`table`、`lang`、`silent`	见通用参数；表格为`tables/cmsy/<渔区>_cmsy_results.csv`、`..._series.csv`和`..._sources.csv`

本文以东北大西洋（FAO 27区）2015–2024年渔获量最大的30个种群为例。按CMSY++建议的起始年运行后，结果如下（p_green为最后一年处在Kobe图绿区，即B ≥ Bmsy且F ≤ Fmsy的概率）：

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	大西洋鲱	Clupea harengus	1950	0.64	1.02	0.03	overfished
2	WHB_27	小鳍鳕	Micromesistius poutassou	1980	0.77	1.62	0.03	overfished
3	COD_27	大西洋鳕	Gadus morhua	1950	0.67	0.83	0.03	overfished
4	MAC_27	大西洋鲭	Scomber scombrus	1966	0.87	1.34	0.14	maximally sustainably fished
5	SPR_27	黍鲱	Sprattus sprattus	1972	0.66	1.02	0.03	overfished
6	POK_27	青鳕	Pollachius virens	1963	0.72	1.00	0.06	overfished
7	CAP_27	毛鳞鱼	Mallotus villosus	1970	0.26	1.05	0.00	overfished
8	HAD_27	黑线鳕	Melanogrammus aeglefinus	1955	0.71	1.09	0.05	overfished
9	HKE_27	欧洲无须鳕	Merluccius merluccius	1950	0.73	0.84	0.09	overfished
10	HOM_27	大西洋竹筴鱼	Trachurus trachurus	1971	0.27	0.99	0.00	overfished
11	PIL_27	沙丁鱼	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	扁海鲽	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	挪威长臀鳕	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	马舌鲽	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	科利鲭	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	欧洲鳀	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	紫贻贝	Mytilus edulis	1978	0.29	0.89	0.00	overfished
26	WHG_27	牙鳕	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	鮟鱇	Lophius piscatorius	1984	0.99	1.10	0.32	maximally sustainably fished
29	ALB_27	长鳍金枪鱼	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", ]

下文详细展示其中两个种群：状态最好的一个和过度捕捞最明显的一个。

状态良好：Northern prawn（Pandalus borealis，FAO和FishBase都没有中文名，PRA_27，渔获量排第14）：最后一年B/Bmsy为1.39（95%置信区间1.08–1.63），F/Fmsy为0.58，处在绿区的概率为99%。
过度捕捞：毛鳞鱼（Mallotus villosus，CAP_27，渔获量排第7）：最后一年B/Bmsy为0.26（95%置信区间0.10–0.54），F/Fmsy为1.05，B/Bmsy低于0.8的概率为100%，处在红区的概率为52%。

毛鳞鱼寿命短，资源量随补充量大起大落，部分年份曾禁止捕捞；仅基于渔获量的模型会把这种渔获量的下降当作生物量的下降。这些结果是使用神经网络B/k先验完成的CMSY++运行，不是经过核实的资源评估，使用前应与ICES的评估结果比较。

两个种群的先验表和结果表由cmsy_tables()生成（见第九步），先是Northern prawn：

show_table(read.csv(file.path(fig_dir, "green_priors_table.csv"), check.names = FALSE))

种群	物种	学名	年份	r 先验	r 先验来源	k 先验（吨）	MSY 先验（吨）	首年 B/k 先验	中间年 B/k 先验	末年 B/k 先验	FishBase / SeaLifeBase 页面
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 (神经网络)	2013: 0.0857-0.361 (神经网络)	0.466-0.837 (神经网络)	https://www.sealifebase.se/summary/25892

show_table(read.csv(file.path(fig_dir, "green_results_table.csv"), check.names = FALSE))

种群	物种	r（95% 置信区间）	k（吨，95% 置信区间）	MSY（吨，95% 置信区间）	Fmsy（95% 置信区间）	Bmsy（吨，95% 置信区间）	末年 B/Bmsy（95% 置信区间）	末年 F/Fmsy（95% 置信区间）	绿区概率	红区概率	状态（FAO）	状态（Palomares 等 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%	最大可持续捕捞	健康

show_table(read.csv(file.path(fig_dir, "overfished_priors_table.csv"), check.names = FALSE))

种群	物种	学名	年份	r 先验	r 先验来源	k 先验（吨）	MSY 先验（吨）	首年 B/k 先验	中间年 B/k 先验	末年 B/k 先验	FishBase / SeaLifeBase 页面
CAP_27	毛鳞鱼	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 (神经网络)	2019: 0.0239-0.23 (神经网络)	0.0372-0.258 (神经网络)	https://www.fishbase.se/summary/252

show_table(read.csv(file.path(fig_dir, "overfished_results_table.csv"), check.names = FALSE))

种群	物种	r（95% 置信区间）	k（吨，95% 置信区间）	MSY（吨，95% 置信区间）	Fmsy（95% 置信区间）	Bmsy（吨，95% 置信区间）	末年 B/Bmsy（95% 置信区间）	末年 F/Fmsy（95% 置信区间）	绿区概率	红区概率	状态（FAO）	状态（Palomares 等 2018）
CAP_27	毛鳞鱼	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%	过度捕捞	严重过度捕捞

result$summary每个种群一行，包括MSY、Bmsy、Fmsy、r、k及其95%置信限，最后一年的B/Bmsy、F/Fmsy、生物量及其置信限，Kobe四个象限的概率，B/Bmsy < 0.8的概率（p_overfished），两种状态分级（FAO两级的status和Palomares et al. 2018五级的status_palomares），所用的全部先验及其来源，以及CMSY++自己的提示（如建议的起始年suggested_start_year）。result$series是逐年结果，包括渔获量、预测渔获量、B/k、生物量、B/Bmsy、F和F/Fmsy及其95%置信限；result$kobe是最后一年的后验抽样；result$details保存r、k和B/k的后验抽样、先验参数和回顾性分析结果；result$settings是所用的CMSY++设置；result$changes是对代码的修改；result$sources记录本次运行所用的CMSY++、JAGS、R、各R包、FAO数据和FishBase的版本与链接，可以直接写进报告。

📊 第六步：图表

cmsy_plot(result)

根据对54篇CMSY应用文献所用图表的统计，cmsy_plot()为每个种群画出CMSY++用户手册中的全部图表，并为多个种群画出常见的汇总图。图表保存在figures/png/cmsy/<种群>/和figures/png/cmsy/summary/下，pptx = TRUE时同时导出可编辑的PowerPoint文件。函数返回这些ggplot对象。

默认只保留坐标轴标题和刻度文字，方便在报告中另配图注。标题、副标题、图注（方法、数据和先验来源）、图例和图内文字（象限概率、PPVR、Mohn’s rho）都可以用annotations打开：annotations = TRUE全部打开，annotations = c("title", "legend")只打开其中几项。options(faofishingreport.cmsy_annotations = TRUE)可以改变本次会话的默认值。

图表	内容
`analysis`	CMSY++分析图六联（手册图1）：渔获量与MSY、探索的r-k组合与先验框、可行的r-k组合、B/k与三段先验、F/Fmsy、平衡曲线
`catch`、`rk`、`biomass`、`b_bmsy`、`exploitation`、`equilibrium`	上述各图的单图版本，以及B/Bmsy（参考线Bmsy、0.8 Bmsy、0.5 Bmsy）
`management`	管理图四联（手册图2）：渔获量与MSY、B/Bmsy、F/Fmsy、Kobe图
`kobe`	Kobe图（手册图3）：首年（方块）、中间年（圆点）到末年（三角）的轨迹，末年后验的50/80/95%概率区，各象限概率
`prior_posterior`	r、k、MSY和三个年份B/k的先验与后验分布，以及PPVR（手册图4）
`retrospective`	去掉最后1–3年后的B/Bmsy、F/Fmsy及Mohn’s rho（Mohn 1999；手册图5，需`cmsy_run(plots = "retrospective")`）
`overview`、`kobe_all`	各种群末年B/Bmsy及置信区间；全部种群的Kobe图（点的大小为末年渔获量）
`status`、`status_time`、`histogram`	各状态的种群数占比和渔获量占比（末年及历年），末年B/Bmsy的分布
`b_bmsy_all`、`trend`、`aggregate`	各种群B/Bmsy小多图；B/Bmsy中位数和低于0.8 Bmsy的种群比例；总渔获量/总MSY和总生物量/总Bmsy
`map`	各FAO渔区种群末年B/Bmsy中位数的地图

参数	含义	默认值
`result`	`cmsy_run()`的结果	必填
`stocks`	要画的种群（见通用参数）	`NULL`：已运行的全部种群
`charts`	要画的图（见上表）	全部
`classes`	状态分级：`"fao"`（两级）或`"palomares"`（五级）	`"fao"`
`annotations`、`lang`、`pptx`、`silent`	见通用参数

状态默认按FAO（2026）对第3层（数据有限）评估的规定分为两级：B/Bmsy ≥ 0.8为“最大可持续捕捞”，< 0.8为“过度捕捞”。classes = "palomares"改用很多CMSY研究采用的Palomares et al. (2018)五级：健康（≥ 1）、轻度过度捕捞（0.8–1）、过度捕捞（0.5–0.8）、严重过度捕捞（0.2–0.5）、崩溃（< 0.2）。

Northern prawn和毛鳞鱼的CMSY++分析图（默认不加注释）：

管理图：

Kobe图（annotations = TRUE，打开了全部注释）：

先验与后验分布。PPVR越小，说明数据给先验带来的信息越多：

回顾性分析：

27区前30个种群的汇总图（annotations = "legend"）：

🔁 MCMC：设置与收敛诊断

mcmc <- cmsy_mcmc(result)
mcmc$settings        # 每个种群的链数、迭代次数、预烧、抽稀和保留的抽样数
mcmc$diagnostics     # 每个参数的后验统计量、Rhat和有效样本量
cmsy_mcmc(result, colors = c("grey30", "#D7191C"))   # 更换各条链的颜色

CMSY++用JAGS以马尔可夫链蒙特卡罗（MCMC）方法从后验分布中抽样。cmsy_run()保存r、k、B/k、B/Bmsy和F/Fmsy每条链的抽样（MSY按CMSY++16.R的算法由r和k计算），cmsy_mcmc()据此做诊断。Rhat（Gelman和Rubin，1992）比较链间方差和链内方差，不超过1.1时视为收敛（Gelman等，2013）。有效样本量（n_eff，用coda计算；Plummer等，2006）是自相关的抽样相当于多少个独立抽样，低于400（Vehtari等，2021）时，95%区间不够精确。未收敛的种群和有效样本量不足的种群分开报告。图中按保留顺序画出每条链的抽样（轨迹：混合良好的链互相重叠、没有趋势）及其密度（各链应当一致）。

下表是两个示例种群分别按CMSY++16.R的设置（2条链，每条迭代60,000次，前30,000次作为预烧舍弃，之后每10次保留1次）和用4条链、每条迭代300,000次（预烧150,000次，每50次保留1次；_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	内禀增长率 r	0.459	1.003	183	0.468	1.001	1,471
PRA_27	环境容纳量 k（吨）	1,070,000	1.002	164	1,050,000	1.001	1,254
PRA_27	MSY（吨）	124,000	1.000	151	123,000	1.001	1,451
PRA_27	第一年的 B/k	0.424	1.003	262	0.432	1.000	1,929
PRA_27	最后一年的 B/k	0.697	1.008	230	0.693	1.000	1,528
PRA_27	最后一年的 B/Bmsy	1.39	1.008	230	1.39	1.000	1,528
PRA_27	最后一年的 F/Fmsy	0.579	1.001	492	0.588	1.000	2,811
CAP_27	内禀增长率 r	0.386	1.055	256	0.383	1.001	2,105
CAP_27	环境容纳量 k（吨）	25,400,000	1.094	124	25,800,000	1.001	703
CAP_27	MSY（吨）	2,460,000	1.034	80	2,480,000	1.000	524
CAP_27	第一年的 B/k	0.331	1.034	216	0.326	1.000	867
CAP_27	最后一年的 B/k	0.128	1.018	142	0.13	1.007	1,437
CAP_27	最后一年的 B/Bmsy	0.255	1.018	142	0.259	1.007	1,437
CAP_27	最后一年的 F/Fmsy	1.05	1.082	262	1.03	1.022	1,926

按默认设置，各条链已经收敛（Rhat最大为1.09），但r、k和MSY的有效样本量只有80至256：r和k的相邻抽样高度相关。用更长的链时，有效样本量升至524至2105，而中位数几乎不变（Northern prawn最后一年的B/Bmsy分别为1.395和1.386，毛鳞鱼分别为0.255和0.259）；运行时间约为原来的五倍（4条链与2条链一样同时运行，但每条链长五倍）。要写进报告的结果，建议用更长、更多的链运行CMSY++，并检查收敛情况：

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)

两个示例种群（默认设置）各条链的轨迹和密度，先是Northern prawn，后是毛鳞鱼：

参数	含义	默认值
`result`	`cmsy_run()`或`cmsy_area()`的结果	必填
`stocks`	要检查的种群（见通用参数）	`NULL`：所有运行过的种群
`rhat_max`	视为收敛的最大Rhat	`1.1`
`n_eff_min`	视为足够的最小有效样本量（所有链合计）	`400`
`plot`、`table`	保存图figures/png/cmsy/<种群>/mcmc.png，以及表tables/cmsy/mcmc_settings.csv和mcmc_diagnostics.csv	`TRUE`
`colors`	各条链的颜色：`NULL`（Okabe-Ito配色）、调色板名称（如`"npg"`、`"Set1"`、`"viridis"`）或每条链一个颜色	`NULL`
`lang`、`pptx`、`silent`	见通用参数

这张图也是cmsy_plot()的"mcmc"图（颜色用mcmc_colors设置）。本包早期版本的运行结果没有保存抽样，需要重新运行CMSY++才能检查。链收敛并不说明先验或模型适合该种群。

🔮 第七步：捕捞情景预测

projection <- cmsy_forecast(result, catch = 100000)

CMSY++16.R本身没有预测功能。cmsy_forecast()用CMSY++自身的模型做预测：Schaefer剩余产量模型，生物量低于0.25 k时生产力线性下降。预测从同一组MCMC迭代的r、k和末年B/k后验抽样出发，加入标准差为sigma的对数正态过程误差（默认0.1，即CMSY++16.R给设置sigmaR的值；模型本身并不使用sigmaR，而是估计其过程误差）。与CMSY++相同，某一年的生物量是该年渔获之前的生物量，F是渔获量除以该生物量，F/Fmsy中的Fmsy在生物量低于Bmsy的一半时线性降低。数据最后一年保留其实际渔获量；从第一个预测年起，F情景的渔获量等于情景F乘以同一年的生物量，因此F/Fmsy就是情景值。默认情景是0.5、0.6、0.8、0.95倍Fmsy（与Froese et al. 2018相同），以及维持末年捕捞死亡率；也可以设定固定渔获量（catch）。数据结束年份不同的种群，各自从自己的最后一年预测到同一个结束年份。

参数	含义	默认值
`result`	`cmsy_run()`的结果	必填
`stocks`	要预测的种群（见通用参数）	`NULL`：已运行的全部种群
`years`	在各种群最晚的最后一年之后预测的年数	`10`
`f`	捕捞死亡率情景，Fmsy的倍数	`c(0.5, 0.6, 0.8, 0.95)`
`status_quo`	是否同时预测维持末年捕捞死亡率的情景	`TRUE`
`catch`	固定渔获量情景（吨/年）	`NULL`
`sigma`	对数正态过程误差的标准差	`0.1`
`seed`	随机数种子	`NULL`
`plot`	是否保存图`figures/png/cmsy/<种群>/forecast.png`和`figures/png/cmsy/summary/forecast_share.png`	`TRUE`
`table`	是否保存`tables/cmsy/forecast.csv`和`forecast_share.csv`	`TRUE`
`lang`、`annotations`、`pptx`、`silent`	见通用参数

函数返回projection（各种群、情景和年份的B/Bmsy、F/Fmsy和渔获量中位数及95%区间，以及B ≥ Bmsy的概率）和share（各情景逐年达到或高于Bmsy的种群比例及95%区间，以及种群数n_stocks；从各种群最晚的最后一年开始）。

固定渔获量取接近各种群MSY的值（Northern prawn 100,000吨，毛鳞鱼2,000,000吨），另加各F情景：

🔍 第八步：比较不同运行或其他评估

result_1950 <- cmsy_run(input, stocks = "PRA_27", seed = 1)        # 从1950年开始
cmsy_compare(`1950` = result_1950, `CMSY++` = result, stocks = "PRA_27")

cmsy_compare()把几次运行（例如不同的起始年或先验，即敏感性分析）以及其他评估的结果（external，如ICES、JABBA或文献中的B/Bmsy）画在同一张图上，并给出各自最后一年的值。起始年的影响很大：Northern prawn从1950年开始运行时，最后一年B/Bmsy为0.82，低于Bmsy；按CMSY++建议从1980年开始时为1.39，高于Bmsy。毛鳞鱼无论从哪一年开始都处于过度捕捞状态（从1950年开始为0.25，从1970年开始为0.26）：

参数	含义	默认值
`...`	`cmsy_run()`的结果，可以命名，也可以不命名（如`default = run1, start_1980 = run2`）
`stocks`	要画的种群（见通用参数）	`NULL`：已运行的全部种群
`external`	其他评估结果的数据框，包括`Stock`、`Year`、`B_Bmsy`（可另加`F_Fmsy`）和`source`（图例中显示的名称）列	`NULL`
`lang`、`annotations`、`pptx`、`silent`	见通用参数

图保存为figures/png/cmsy/compare/<种群>.png；函数返回last_year（每个种群每次运行一行）和plots。

📋 第九步：报告用的表格

tables <- cmsy_tables(result)

cmsy_tables()生成CMSY文献中常见的两张表，即上文第五步中的先验表和结果表。数值已取整，范围写成“下限-上限”，列名使用所选语言，保存为tables/cmsy/cmsy_priors_table.csv和cmsy_results_table.csv。结果表同时给出FAO（2026）两级和Palomares et al. (2018)五级的状态。

参数	含义	默认值
`result`	`cmsy_run()`的结果	必填
`stocks`	表中的种群（见通用参数）	`NULL`：已运行的全部种群
`table`、`lang`、`silent`	见通用参数

⚙️ CMSY++自带的设置

CMSY++代码里有一组“一般设置”，例如渔获量观测误差CV.C、MCMC链数和迭代次数。其中sigmaR和ni两项在CMSY++16.R中有赋值，但代码并不使用（模型的过程误差由先验itau2 ~ dgamma(4, 0.01)估计）：它们标为“not used”（不使用），不能修改，以免误以为做了过程误差的敏感性分析。cmsy_run()默认完全沿用CMSY++原来的数值，cmsy_defaults()列出每一项的原默认值、实际使用的值和含义：

show_table(cmsy_defaults()[, c("setting", "cmsy_default", "applies_to", "description")], n = 30)

setting	cmsy_default	applies_to	description
CV.C	0.15	all	渔获量观测误差的变异系数
CV.cpue	0.2	BSM	CPUE 的最小变异系数
sigmaR	0.1	not used	CMSY++16.R不使用：代码中有赋值（sigR <- sigmaR），但从未使用；模型的过程误差由先验 itau2 ~ dgamma(4, 0.01) 估计。不能修改
cor.log.rk	-0.76	charts	log r 与 log k 的相关系数，仅用于作图
rk.cor.beta	2.52, 3.37	all	r-k 相关系数加 1 的 beta 先验
nbk	3	BSM	BSM 使用的 B/k 先验个数：1（首年）、2（加中间年）、3（再加末年）
bt4pr	FALSE	BSM	用丰度数据设定 B/k 先验
ct_MSY.lim	1.21	all	渔获量与 MSY 先验之比超过此值时，B/k 先验保持不变
q.biomass.pr	0.9, 1.1	BSM	丰度数据为生物量时 q 的先验范围
n	5000	charts	分析图中 r-k 点云的点数
ni	3	not used	CMSY++16.R不使用：代码中有赋值，但从未使用。不能修改
nab	3	BSM	运行 BSM 所需丰度数据的最少年数
bw	3	all	渔获量平滑（ksmooth()）的带宽
e.creep.line	TRUE	charts	在生物量图中显示未按渔效增长（effort creep）校正的 CPUE
select.yr	NA	output	另外报告该年 B、F、B/Bmsy 和 F/Fmsy 的年份（NA 为不报告）
n.chains	2	all	MCMC 链数（JAGS）
n.iter	60000	all	每条链的 MCMC 迭代次数（JAGS）
n.burnin	30000	all	预烧（burn-in）迭代次数（JAGS）
n.thin	10	all	MCMC 抽稀间隔（JAGS）

applies_to为”BSM”的设置只在有丰度数据（CPUE或生物量）时才起作用。FAO数据只有渔获量，所以真正影响结果的是标为”all”的设置。需要改时，用settings只写要改的几项，其余保持原值；每一项都会先检查（名称、取值范围、迭代次数是否足够）：

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, ]          # 改动的设置

参数	含义	默认值
`settings`	要修改的设置，命名列表，如`list(CV.C = 0.2)`；与默认值并列显示并检查	`list()`
`lang`	`description`列的语言：`"en"`、`"zh"`或`"both"`	`NULL`：使用`fao_settings()`设定的语言

实际使用的设置保存在运行目录的settings_used.csv，也记录在result$settings和result$sources中。绘图和输出开关（mgraphs、kobe.plot、retros等）由cmsy_run()的plots参数控制，不在settings里设置。

🖥️ Windows、macOS和Linux

同一段代码在三个系统上都能直接运行。cmsy_system()自动识别系统，并做相应修改，例如macOS上不调用windows()、不使用antialias = "cleartype"，修改清单见第一步。JAGS安装程序会按系统和R版本自动选择。每次修改CMSY代码，GitHub Actions上的”CMSY-cross-platform”检验都会在三个系统上用cmsy_download_jags()下载并安装匹配的JAGS、真实运行CMSY++、运行全部CMSY测试，并保存各系统的输出，方便核对。

2026年10月的检验结果，用的是当时本文示例渔区34区的两个种群（每种组合170项测试全部通过；种子相同，cmsy_area()运行大眼裸颌鲈GRB_34和沙丁鱼PIL_34）：

系统	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芯片）	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芯片）	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芯片）	4.6.1	5.0.0	5-3	1.392	0.672	0.740	1.317

JAGS 4.3.x和rjags 4-17相同时，三个系统的结果逐位相同。换成JAGS 5或较旧的rjags后，MCMC抽样不同，数值略有差别（GRB_34的B/Bmsy在1.37–1.42之间，都在95%置信区间1.08–1.65内），但两个种群的状态判定都不变：GRB_34为最大可持续捕捞，PIL_34为过度捕捞。要完全重现某次结果，需要使用相同版本的JAGS和rjags；所用版本记录在result$sources和tables/sources.csv中。

⏱️ 进度条与silent

cmsy_run()、cmsy_plot()、cmsy_forecast()和cmsy_area()运行时会显示进度条，例如：

CMSY++ [==================            ] 18/30 ( 60%) 00:51 elapsed, 00:34 left COD_27

所有会输出信息的CMSY函数都有silent参数：silent = TRUE时不显示进度条和提示信息（cmsy_setup()需要确认时仍会询问）。fao_settings(silent = TRUE)可以设定本次会话的默认值。

📝 写入Word报告

cmsy = TRUE时，fao_report()会用CMSY++评估该渔区的主要物种（见📝 一键生成评估报告（Word、LaTeX、PDF）），并在报告中加入资源评估部分，给出每个物种的结果，以及参考点、先验及其来源的附表。它对报告的top_n个主要物种（默认10个）运行cmsy_area()，不做情景预测。cmsy_area()和cmsy_input()的设置以cmsy_为前缀传入，例如cmsy_top = 5、cmsy_default_resilience = "Medium"或cmsy_settings = list(n.iter = 100000)。CMSY++的全部结果（运行、图表和表格）保存在报告旁边的cmsy/文件夹中，其中的cmsy_result.rds可以用来重写报告而不必重新运行CMSY++：cmsy = readRDS("cmsy/<渔区>/cmsy_result.rds")，路径从工作目录算起（默认是reports/cmsy/<渔区>/cmsy_result.rds）。以前用cmsy_area()得到的结果也可以这样传入。

report <- fao_report(27, cmsy = TRUE, cmsy_top = 5)
report$cmsy_reference_points                  # r、k、MSY、B/Bmsy等及其95%置信限
fao_report(27, cmsy = readRDS("reports/cmsy/Atlantic__Northeast/cmsy_result.rds"))

📖 参考文献

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