FRA-SA2024-ABCWG02-05

令和 6(2024)年 度再生産関係の決定に関するガイドライン
2024 年 6 月 21 日 1
資源評価高度化作業部会(文責 宮川光代・市野川桃子）
再生産関係は将来予測や管理基準値の推定に影響する最も重要な要因である一方，再生産関係のパラメータは資源評価モデルの中で推定が最も難しいパラメータの一つである (Walters 2004)．また環境要因な
どの要因が再生産関係に大きな影響を及ぼすことが明らかになっているが，その詳細なプロセスについてはほとんど解明されていない．そのため，利用可能なデータのみでは十分な確度をもって再生産関係の形を選択したりパラメータ推定ができない場合が多い．その場合には様々な要因を考慮して再生産関係を選択する必要がある．ここでは，資源評価モデルで推定された加入尾数と親魚量の関係をもとに MSY
管理基準値やそれに準じた管理基準値を計算する資源（1A，1B ルール，水産研究･教育機構 2024）において，再生産関係を決定する際の基本的な考え方と手続きについて述べる．
1． 用いるデータ範囲と基本的な考え方
a. 基本的には，過去に推定された加入尾数と親魚量の利用可能な時系列データを全て用いる．
b. ただし，レジームシフトなど明らかな根拠がある場合は，再生産関係の推定に用いるデータを特定の期間のみのものとするような方法も考えられる．その場合は，抽出した期間に明確な科学的根拠があること，それを報告書に詳述するか適切な引用を明記すること，関係科学機関や有識者間で合意することが必要である．明確な科学的根拠とは，再生産関係のシフトを示唆する既往研究だけでなく，
AIC 等の情報量規準をもとにした統計的根拠もあわせて示す．さらに，不確実性に留意し，簡易的なシミュレーション，感度試験などで不確実性に対する頑健性を調べることが推奨される．
c. 資源評価最終年またはそれに準じる年の加入尾数については，統計解析によって不確実性が非常に高いことが明らかになっている場合，再生産関係の推定から除くことも考えられる．その場合，シミュレーションや感度試験により，除いた場合と除かない場合のリスクの評価を行うことが推奨される．
d. 再生産関係は利用可能なデータと既存の科学的知見を最大限活用して最善のものを選択することを基本とし，想定外のデータの更新や持続性を脅かすほどの極端な加入が起こった場合（付録．管理期
間内に管理基準値や漁獲管理規則を変更するためのガイドライン）以外には，原則的に 5 年間の管理期間内で同じものを用いる．
e. MSY 管理基準値は再生産関係の仮定に大きく影響されるため，その決定には科学的に十分な根拠と客観性が求められる．複数の再生産関係の候補のうち一つを選択するための客観的な根拠がなく，また，その選択によって MSY 管理基準値が大きく変わるような場合には，どの選択をしても管理を大1

令和 5 年度版からの変更点：1) 1A ルール適用事例集に令和 4 年度の事例で抜けていたキンメダイ太平洋系群とカタクチイワシ瀬戸内海系群を追加，2) 引用文献の体裁の調整，3) 前書きにガイドラインの目的を加筆．
1

FRA-SA2024-ABCWG02-05
きく失敗することがない頑健な管理が実現されるように再生産関係を選択するか（3f， 3g， 3i）や管理基準値そのものを再生産関係によらずに決定する（1B ルール）方法もある．再生産関係によらずに管理基準値を決定する 1B ルールの適用範囲については「4.1B ルールの適用範囲」に詳述した．
2． 候補となる再生産関係式と注意点
a. ホッケー・スティック型再生産関係（HS，Clark et al. 1985），ベバートン・ホルト型再生産関係（BH,
Beverton and Holt 1957），リッカー型再生産関係（RI，Ricker 1954）など，適切な引用がある再生産
関係式を候補とする．
b. 基本的には，密度効果が見られない資源や，極端な密度効果がある資源においても現実的な再生産関係が得られる HS の利用を推奨する（Ichinokawa et al．2017）．しかし，「3 再生産関係を選択するときの基準」で解説する基準に照らしたときに，使用の利点が明らかである場合にはその他の再生産関係式やそれらのモデル平均も使用できる．
c. 加入変動に影響を及ぼす環境変数やその他の要因の導入によって管理上の利点が大きい場合には，適切な科学的根拠・引用とともに，再生産関係モデルにこれらの要因も導入する．
d. 再生産関係のパラメータ推定の際には，補足資料 1 から 3 を参考に，一連のモデル診断手法をあて
はめ，パラメータ推定が適切かどうかを確認する（図 1）．また，この結果を研究機関会議資料の補足として掲載する．
資源の生産性を測る指標として，SB0（漁獲がないときの平均親魚量），R0（SB0 のときの平均加入尾数），h（スティープネス）を計算し、これらの推定値と不確実性を示す．h は BH と RI を使う場合
は 0.2SB0 のときの平均加入尾数を R0 で割った値（Hilborn and Walters, 1992），HS の場合には 1SBhs/SB0（SBhs は HS の折れ点）（Punt et al．2014）として計算する．h は小さいほど密度補償効果（資源量が少なくなると親魚あたりの加入尾数が多くなる）が小さくなることを示す指標で，資源の漁獲に対する抵抗性の強さを示している（Hilborn and Walters，1992）．BH を用いたときは，条鰭綱における h の平均値が 0.74（標準偏差 0.23）程度になることがメタ解析の結果から知られている（Thorson2020）．また，h の値は解析的に Fmsy に直接影響することが知られている（Mangel et al．2013）ため，h の値が小さすぎる・大きすぎる場合には非現実的な Fmsy の推定につながり，h の不確実性が非常に大きいことは Fmsy の不確実性が大きいことを示唆しているため，h の値とその不確実性を，加入量の残差のブートストラップ解析等を実施することにより示すことが重要である．
3． 再生産関係を選択するときの基準仮定する再生産関係は，過去に見られた再生産関係を記述するためのものであると同時に，将来の親魚量の変化に対して加入尾数がどのように変化するかを予測するために用いられ，管理基準値に大きな影響を与えることになる．そこで，再生産関係の決定は，予測力に加えて，以下に挙げる他のさまざまな要
因も勘案した上で行うことが望ましい．
a. 予測力・

赤池情報量規準（AIC），または，小標本サイズによるバイアスを補正した bias corrected AIC (AICc)（Hurvich and Tsai 1989）などをモデルの予測力を測る規準として用いる．たとえば AICc または2

FRA-SA2024-ABCWG02-05
AIC を使う場合では，おおまかに，AIC 最小モデルに対して候補モデルの AIC の差（ΔAIC）が４以上であれば，候補モデルの確からしさは小さいと言われている（Burnham and Anderson
2002）．
・

代表的な事例：H31 年度研究機関会議ホッケ道北系群（森田ほか 2019）b. 生物学的妥当性 or 便宜的仮定
・

RI は，親魚が多すぎる場合に加入尾数が減少するような強い密度効果を持つ再生産関係を表現できる．RI において特に強い密度効果（たとえば h が 1 を大きく超えるような場合）が推定された
場合には，対象種の生活史を考慮して，推定されたような強い密度効果が起こりうるかを検討することが重要である．

・

加入尾数と親魚量が直線関係あるいは加入尾数と親魚量が無相関の場合，BH や RI ではパラメー
タ a，b に高い相関が見られ，解が一意に定まらなくなる．そのような場合でも，HS では最大親魚量よりも大きい親魚量の範囲において加入を一定と仮定する，あるいは最小親魚量よりも小さい親魚量の範囲において直線的に加入が減少すると仮定することで，現実的な範囲での加入の予測値を与えることが可能である（Ichinokawa et al．2017）．ただし，これらの仮定には生物学的な背景がなく，便宜的な仮定であること (Walters and Martel

2004)，また，ブートストラップ等を用いて推定された h の不確実性は（本来なら過去最小親魚量以下，過去最大親魚量以上のところに真の HS の折れ点があるかもしれないので）過小評価されていることに注意する必要がある．
・

代表的な事例：H31 年度研究機関会議スケトウダラ太平洋系群（HS の利用）（境ほか 2019），R3
年度マダラ本州太平洋北部系群（RI の利用（成松ほか 2021））

c. 外れ値に対する頑健性・

最小二乗法（正規分布を仮定した最尤法に対応， 以下 L2 と表記）を用いたパラメータ推定は
外れ値の影響を受けやすいことが知られている．外れ値がパラメータ推定に大きな影響を与えていると考えられるケース等，特に管理基準値の頑健性を重視するような場合は，外れ値の影響を受けにくい最小絶対値法（ラプラス分布を仮定した最尤法に対応， L1 と表記）によるパラメータ推定を利用したほうが良い．

・

代表的な事例：H31 年度研究機関会議ホッケ道北系群（森田ほか 2019）

d. 観察された最大親魚量以上で，加入尾数が過去に観察された最大尾数以上の極端な外挿値になるような場合の回避・

親魚量の観察範囲内で密度効果がほとんど認められない場合には，再生産関係が直線に近い形になる．そのため，BH や RI では最大親魚量以上の範囲で，加入尾数の期待値が過去に観察された最大加入尾数よりも極端に大きくなるような再生産関係が推定されることがある（図 2）．
それにより MSY 管理基準値も非現実的に過大な値になってしまうため，このような場合には，
BH や RI の使用は避け，加入尾数の期待値を過去最大加入以上に外挿しない性質を持つ HS を選ぶことが望ましい．

・

代表的な事例：H31 年度研究機関会議スケトウダラ日本海系群（山下ほか 2019）3

FRA-SA2024-ABCWG02-05

e. 観察された最小親魚量以下で加入尾数が保守的でない外挿値になるような場合の回避・観察された範囲内で親魚量と加入尾数のあいだに明瞭な正の相関関係が見られない場合，または，負の相関関係が見られるような場合には，過去最低親魚量以下の範囲において，親が減っても加入が減らない，または，増加するような予測値が得られる場合がある（図 3）．過去に経験したことがないくらい低い親魚量において，加入尾数が保守的でない外挿値になる場合にはリスクが高いため，予防的な観点から，このような再生産関係の使用は避けるほうが望ましい．

・

代表的な事例：H31 年度研究機関会議スケトウダラ太平洋系群（山下ほか 2019）
・マサバ対馬暖流系群（安田ほか 2019）
・ゴマサバ太平洋系群（由上ほか 2019）f. 推定された管理基準値の頑健性
・複数の再生産式，最適化手法を用いた再生産関係を比較する際には，あわせて MSY 管理基準値の計算もおこない，MSY 管理基準値の再生産関係式の選択に対する頑健性（どの再生産関係を選んでも比較的同程度の MSY 管理基準値が得られる）を確認することが望ましい（図 4）．

・

一方で，尤度の差が非常に小さいにもかかわらず MSY 管理基準値が大きく異なる場合には，再生産関係を両者の尤度によって重みづけ平均するモデル平均（複数の再生産関係のモデル平均を用いた管理基準値推定，FRA-SA2020-BRP01-08）や（図 5），1B ルールの利用（4 節），g) リスクの非対称性の検討といった方法を用い，頑健な管理方策を与えるような方法を考慮する．

・

代表的な事例：H31 年度研究機関会議ゴマサバ東シナ海系群・スケトウダラ太平洋系群，R2 年度
研究機関会議マアジ太平洋系群（モデル平均，FRA-SA2020-BRP01-08）

g. 異なる再生産関係を用いた場合のリスクの非対称性
・

予測力やその他の点で同等と考えられる再生産関係が複数ある場合には，異なる再生産関係を用いて計算した管理基準値・漁獲管理方策をもとに管理を実施した場合のリスクに非対称性（再生産関係 A， B がどちらも確からしい場合，誤って A を使った場合と誤って B を使った場合
で，将来予測におけるリスクに大きな違いがある，など）があるかどうかを検討することで，再生産関係を選択する方法も考えられる（図 6）．

・

代表的な事例：H31 年度研究機関会議ゴマサバ東シナ海系群（林ほか 2019）・ゴマサバ太平洋系
群（由上ほか 2019）h. 自己相関・

選択した再生産関係からの予測値と観測値との残差において，有意な自己相関関係が認められるかどうかを検討する．自己相関が有意である場合には，管理基準値計算や ABC 計算における将来予測において，加入プロセスに自己相関を考慮することが望ましい．これは，自己相関が有意である場合，悪い加入が続けて発生するような事態がより高い確率で発生するため，基準値以下に資源が下がるリスクがそれだけ高くなり，より保守的な管理が必要になる場合が想定されるためである（図 7）．

・

自己相関を考慮しない場合には AR0，過去 1 年の残差のみを考慮する場合には AR1 と表記する
4

FRA-SA2024-ABCWG02-05（AR は auto-regression の頭文字）．
・

再生産関係パラメータと自己相関係数は時系列解析の方法（正確尤度など）を使って同時推定することも可能であるが，計算の安定性を考慮し，AR を考慮しないでパラメータ推定したのちに残差に対して自己相関係数を推定する方法（2 段階推定）を基本的に用いることとする．これにより，データが 1 年増えただけでパラメータの推定値が大きく変わってしまうようなことが防げる．また，
現行の漁獲管理規則のもととなった MSE においても 2 段階推定が用いられている．
但し，2 段階推定と同時推定で推定される自己相関係数が大きく異なり，結果として管理基準値も大きく異なる場合には，同時推定の利用も考慮する．FRA-SA2020-BRP01-06 は，同時推定のほうが推定のバイアスが小さい利点がある一方で，分散が大きい欠点があることを示している．
同時推定を用いる場合には，最適化手法が最小二乗法であること，パラメータが十分収束していることを確認し，また，自己相関パラメータの推定値の不確実性が管理基準値に与える影響を十分考慮する必要がある．

・

自己相関を再生産関係の推定の際に考慮した場合，将来予測においても過去年の残差を引き継
ぐような加入尾数を推定することで，傾向を持った残差が将来も起こりうる状況を表現できる．
実際の数式については「再生産関係の推定・管理基準値計算・将来予測シミュレーションに関する技術ノート（令和 6 年度）（FRA-SA2024-ABCWG02-04）」を参照のこと．

・

自己相関が非常に大きい場合，推定される再生産関係が親子関係のデータに全く当てはまっていないように見える場合がある．このような場合には，自己相関分を観測データから差し引いたプロットを作成することで，再生産関係の当てはまりを視覚的に示すことが望ましい（図 8）．

・

代表的な事例：H31 年度研究機関会議マサバ太平洋系群（西嶋ほか 2019），R2 年度研究機関会議
マアジ太平洋系群（井須ほか 2020）

i. 管理方策の頑健性（MSE による評価）
・

自己相関の考慮だけでは加入の予測値からのずれに対処できず，また，その問題が短期的な将来予測や ABC 推定に影響を与えることが懸念される場合には，再生産関係に複数のシナリオをあてはめ，どのシナリオが真だったとしても頑健な管理方策・管理基準値を選択する方法をとる．
これによって，再生産関係には不確実性があるものの，管理方策においては頑健なものを選択することができる．

・

代表的な事例：H31 年度研究機関会議ホッケ道北系群（森田ほか 2019）（図 9），R3 年度研究機関
会議ブリ，簡易的 MSE を用いた複数の管理基準値の頑健性の比較・HCR の検討（古川ほか 2021））

4．1B ルールの適用範囲
1B ルール（詳細は「令和 6 年度漁獲管理規則および ABC 算定のための基本指針」FRA-SA2024ABCWG0201 を参照）の適用を検討する場合の例を以下に示す．
1.

再生産関係のパラメータを推定するための親魚量・加入尾数のデータ数が充分にない.

2.

信頼できる再生産関係のパラメータを推定できず，結果として MSY 管理基準値の不確実性が非常に大きくなる場合．信頼できる再生産関係のパラメータを推定できない場合とは，例えば，①モデル選
5

FRA-SA2024-ABCWG02-05
択の不確実性が高い，②再生産関係の推定パラメータの信頼区間が広すぎる，③パラメータ間の相関が強く，一意なパラメータ推定が困難，④h の不確実性が高すぎるなど．ただし①～④のどれかがあてはまっていたとしても，得られる MSY 管理基準値の推定値が頑健であるような場合においては，便宜的仮定（たとえば HS の利用やパラメータの固定）のもとで再生産関係を一つに決めて MSY 管
理基準値を得る方法もある．
3.

過去に得られた親魚量または資源量の範囲が SB0 またはそれに準じる値と比べたときに非常に狭い
範囲となる

4.

対象とした系群内でみられる加入が他海域からの移入に大きく依存している可能性が高い

5.

資源評価モデル内で再生産関係を推定できるモデル（例：Stock Synthesis，SS や State-Space Assessment
Model，SAM など）を使っているが，再生産関係のパラメータ推定が不安定なためにモデル内で再生
産関係のパラメータ（h など）を仮定しており，その仮定によって MSY 管理基準値が大きく変わる
場合．

6.

加入尾数の推定は比較的頑健であるが，それに対応する親魚量の推定精度が悪いと考えられる場合（例えば，卵から加入に至るまでの期間が不明なため，加入がどの親魚量と対応しているか不明（ズワイガニなど），年齢データの不確実性により親魚量の推定精度が加入尾数の推定精度よりも低いと考えられる場合）．

ただし，上記のような例以外でも，何等かの理由によって MSY 管理基準値を計算するための再生産関係の決定が困難である場合，それが客観的に説明され，関係研究機関の合意があれば 1B ルールを適用する．
令和 4 年度以降には 3 系群で 1B ルールが適用された（6. 1B ルール適用事例集）．
5． 1A ルール適用事例集
平成 31 年度研究機関会議（7 系群）
系群

再生産

選択基準

説明

引用文献

3e

HS は AICc 最小モデル（RI, ΔAICc=8.1)よりも予測力

境ほか（2019）

関係
スケトウ

HS,

ダラ太平

L2，

が低いが，親魚が少ない範囲で RI が加入を過大評価

洋系群

AR0

するリスクを避けた．
3b，3f

RI の場合 2010～2014 年の 5 年間のデータのみが曲線の形を決めるため，生物学的背景の吟味および関係式間での MSY 管理基準値算出値の比較の上で RI の使用を避けた．
HS，RI，BH ともに直線的な再生産関係になり AICc

山下ほか

L2，

はほぼ変わらないため，観察された最大親魚量以上（2019）海系群

AR0で極端な加入の外挿を避ける HS を選択．

ホッケ道

HS，

北系群

スケトウ

HS，

ダラ日本

3d

3a，3c

AICc が最小のモデルを選択．結果として，L1 が選択

森田ほか

L1，

されたため，外れ値に対しても頑健なモデルとなっ（2019）

AR0

た．
6

FRA-SA2024-ABCWG02-05
3i

簡易的な MSE を行い，近年の低い加入を想定した場合のリスクを漁獲管理規則で推奨するβにおいて考慮した（図 9）．

マサバ太

HS，

平洋系群

3h

残差の自己相関が有意だったため，将来予測では加

西嶋ほか

L2，

入の自己相関構造も考慮した．これにより加入の自（2019）

AR1（2己相関から生じる潜在的なリスクを管理基準値に取段階）

り入れることができた（図 7）．また，近年の残差は正であったため，直近年の将来予測は加入が通常よりも良い状況となることが仮定された．
HS は AICc 最小モデル（RI，ΔAICc=5.9）よりも予測

安田ほか

L2，

力が低いが，親魚が少ない範囲で RI が加入を過大評（2019）

群

AR0

価するリスクを避けた（図 10）．

ゴマサバ

HS，

太平洋系

マサバ対

HS，

馬暖流系

3e

3e, 3g

RI は AICc 最小モデル（BH or HS，ΔAICc=6.4 or

由上ほか

L1，

6.5）よりも予測力が低いが，親魚が少ない範囲で RI（2019）

群

AR0が加入を過大評価するリスクを避けた．

ゴマサバ

HS，

3f

図 4 の解説を参照．

林ほか

東シナ海

L1，

3g

図 6 の解説を参照．

（2019）系群

AR0(注) 先行 7 魚種における管理基準値計算のための将来予測の加入の残差分布としては，すべての系群に対して対数正規分布が仮定された．
令和 2 年度研究機関会議（8 系群）
系群

再生産

選択基準

説明

引用文献

3a

最小二乗法を用いた場合で AICc が最小となる HS を古市ほか
（2020）

関係
マイワシ

HS，

太平洋系

L2，

選択．また，期間を分けない場合（HS，L2，AR1）

群

AR0，

の AICc は 108.30，期間を分けた場合は 96.04 で，10レジー以上の差が見られた．

ム
マイワシ

HS，

対馬暖流

L2，

系群

AR0，

3a

期間を分けない場合の AICc（HS,

L2,

AR1）は

向ほか

124，分けた場合の AICc は 119．

（2020）

AICc がほとんど変わらない 2 つの再生産関係間で管安田ほか（2020）

レジー
ム
3f

マアジ太

モデル

平洋系群

平均

理基準値が大きく異なるため，モデル平均した再生

RI・

産関係を用いた（図 5）．選択された管理基準値の頑BH，健性は簡易 MSE にて確認した．
7

FRA-SA2024-ABCWG02-05
L2，3h

同時推定の場合に自己相関が非常に高く推定される

AR1

ため，その不確実性を考慮した上で，同時推定によ（同るパラメータ推定を行った（図 8）
．

時）
マアジ対

HS，

馬暖流系

L1，

群

AR0

ズワイガ

RI，

ニ日本海
系群 A 海3a， 3cAICc が最小となるモデルを選択した．

依田ほか（2020）

3a， 3b

AICc が最小となるモデルを選択した．その際，RI で

上田ほか

L1，

は BH や HS に比べて強い密度効果が推定されたが，（2020）

AR0

本種では大型個体による小型個体の共食いが見られていることから、RI のあてはめは生物学的にも妥当域と考えられた．
ズワイガ

HS，

3a，

AICc が最小の HS，最小二乗法，自己相関ありを選

森川ほか

ニ太平洋

L2，

3d， 3e択．それにより，3d，3e も考慮した．

（2020）

系群

AR1（2
3b

産卵後に死亡するなどの生態学的特徴を考慮して，

久保田ほ

RI は再生産関係の候補から外した．

か

尤度が似た HS と BH の間で簡易 MSE を実施し（報（2020）

段階）
スルメイ

HS，

カ秋季発

L2，

生系群

AR0

3i

告書詳細版，補足資料 5），頑健な管理基準値（HS）
を与える再生産関係を選択した．
スルメイ

BH，

カ冬季発

L1，

生系群

AR0

3b
3i

産卵後に死亡するなどの生態学的特徴を考慮して，

加賀ほか

RI は再生産関係の候補から外した．

（2020）

尤度が似た HS と BH の間で簡易 MSE を実施し（報告書詳細版，補足資料 3），頑健な管理基準値（BH）
を与える再生産関係を選択した．

令和 3 年度研究機関会議（12 系群）
系群

再生産

選択基準

説明

引用文献

RI，

3a，

AICc が最小となるモデル(RI)を選択 (3a)．しかし，

古川ほか

L2，

3g，3i

➀RI のような強い密度効果が生じるのかは明らかで（2021）

関係ブリ

AR0はない，②BH および HS との AICc の差は 4 未満のため，簡易 MSE を行った (3g，3i)．それにより，真の再生産関係が HS もしくは BH のときに誤って RIと想定して管理を行っても，短期的にも，中長期的にも漁獲量および資源量の減少リスクはほとんどな
いため RI を選択．
8

FRA-SA2024-ABCWG02-05
期間を分けない場合（BH，L2，AR1）の AICc は

木下ほか

L2，

32．8 期間を分けた場合は 26.11 であり (3a)，通常加（2021）

AR0，

入期と高加入期で分割するモデル (1b) を採用．

カタクチ

BH，

イワシ太
平洋系群

3a，1b

レジー
ム
カタクチ

HS，

2d，

HS と BH における自己相関 (3h) の同時推定の AICc

日野ほか

イワシ対

L2，

3a，

は同じ (3a) であったが，当てはめる上での個々のデ（2021）

馬暖流系

AR1 (同3b，3hータの影響をジャックナイフ法により検討したとこ

群時)

ろ，HS では推定値の頑健性に大きな問題がなかった
が，BH では推定パラメータが大きく変化した (2d)．
また，BH では親魚量が低水準の際の観察データが少ないことによるパラメータの推定の不安定性 (3b, HSの利用) がみられたため，HS を採用．

ウルメイ

HS，

3a，

HS で自己相関の同時推定 (3h)，RI で自己相関の二

依田ほか

ワシ対馬

L2，

3h，

段階推定，BH で自己相関の二段階推定の三つのモデ（2021）

暖流系群

AR1(同

3g，3i

ルの予測力に差は認められないため (3a)，不確実性に対する頑健性と，それぞれの再生産関係を誤って時)選択した場合のリスク(3g) を簡易 MSE で評価した
(3i)．
ヤナギム

HS，

3a，

HS， RI， BH における同時推定 (3h)の AICc に大き

三澤ほか

シガレイ

L2，

3b，

な差はなかったが (3a)，ジャックナイフ法による検（2021）

太平洋北

AR1(同

3h，2d

討を行ったところ (2d)，BH では影響が大きかった

部系群

時)

が，HS および RI では影響が BH 小さく安定していたため，HS を採用 (3b)．

マダイ日

HS，

3b，

RI およびを仮定した場合は，親魚量が過去最低値よ

下瀬ほか

本海西

L2，

3e，3h

り低い場合に加入尾数が保守的でない外挿値となる（2021）

部・東シ

AR1(同

ナ海系群

時)

マダラ本

HS，

3a，

HS と BH で AICc は同値となり RI より低かった

佐久間ほ

州日本海

L2，

3b，3e

(3a)．観測範囲内で明瞭な親子関係が見られず BH で

か

北部系群

AR0

はパラメータが一意に定まらなかったため過去最小（2021）

ことから HS を選択 (3b，3e)．

親魚を変曲点とする HS を選択 (3b，3e)．
ムシガレ

HS，

3a，

BH および RI では，観察された親魚量の範囲では右

八木ほか

イ日本海

L2，

3b，

肩上がりの直線として推定され密度効果が認められ（2021）

南西部系

AR1(二

3d，3h

なかったため過去最大親魚量を折れ点とする HS を選

群

段階)

ソウハチ

HS，

3a，

BH，RI，HS の AICc に大きな差はなく(3a)，親子関

飯田ほか

日本海南

L2，

3b，3h

係に明瞭な密度効果が認められないため，そのよう（2021）択 (3b, 3d)．

9

FRA-SA2024-ABCWG02-05
西部系群

AR1(二

な資源でも現実的な再生産関係が得られる HS を採用

段階)

(3b)．

ヒラメ瀬

HS，

3a，

HS では，折れ点を観測範囲から求めない場合に最も

山下ほか

戸内海系

L2，

3b，

尤度が小さくなった（最小親魚量を変曲点とする）（2021）群

AR1(同3e， 3h(3a，3b， 3e)．RI では，観測された最小親魚量より少ない状況において最大の加入量が予測．BH では観時)

測値から変曲点が推定されなかった．
マダイ瀬

RI，

3a，

BH，RI，HS の AICc に大きな差はないが (3a)，HS

山本ほか

戸内海

L2，

3b，3h

と BH では推定の不確実性が高く，特に BH では再生

（2021）

中・西部

AR1(同

産関係のパラメータ間に高い相関があり解が一意に系群

時)計算できなかった．また，従来の資源評価結果から
親魚量と天然加入量の増減に反比例の関係が観察されており，さらに，生活史のいくつかの段階で密度過多の影響が示唆されていることより RI を採用
(3b)．

マダラ本

RI，

3a，

BH，RI，HS の AICc に大きな差はないが (3a)，ジャ

成松ほか

州太平洋

L2，

3b，

ックナイフ解析より (2d) BH，HS では影響が大きか（2021）

北部系群

AR1(同

3h，2d

ったが，RI では影響が小さく安定．さらに，RI のよ
うな密度効果を示唆するような要因が確認されてい時)るため RI を採用 (3b)．
令和 4 年度研究機関会議（8 系群）
系群

再生産

選択基準

説明

引用文献

関係
ヒラメ太

HS, L1,3a, 3b, 3g,

HS, BH, RI の AICc はそれぞれ 32.4, 32.7, 34.3 と大

富樫ほか

平洋北部

AR0

3h

差はなかった (3a)．また，異なる再生産関係を用い（2022）た場合のリスクの非対称性の評価においても，将来予測のリスクに大きな違いはみられなかった (3g)．
また MSY 管理基準値は HS のほうがより現実的な目標となっているため HS を採用 (3b)．
HS, BH, RI の AICc はそれぞれ 12.7, 12.7, 2.0 であ

八木ほか

AR1(同

り， RI を当てはめた場合で最も低かった．しか（2022）

時)

し，RI を仮定した場合は，親魚量が観測範囲より低

ヒラメ日HS, L2,

本海北部

3b, 3e, 3h

い状態で高い加入量が予測され，加入尾数が保守的でない外挿値となるため，RI を選択するリスクは高く回避することが望ましい (3e)．また，RI で推定された強い密度効果が生物学的に妥当であることを
10

FRA-SA2024-ABCWG02-05
示す情報は得られていない (3b)．HS では折れ点を観測範囲から求めない場合に最も尤度が小さいた
め，観測範囲の最小親魚量を変曲点とした (3b)．
BH では観測値から変曲点が推定されなかった．
サワラ瀬HS, L2,戸内海

AR03a, 3b，3d

RI と BH では，親魚量が増加するほど加入量が増え

安田ほか

続け，明確な密度効果はみつからず，管理基準値の（2022）

過大な外挿につながるため候補から外した (3b, 3d).
また，HS で L2 が最も AICc が低かった (3a).
BH, HS, RI の順で AICc が低く，また L2, AR1 (同時)

山本ほか

AR1(同

が最も低かった．RI は AICc の差も大きいことから（2022）

時)

除外．予測力(3a), プロファイル尤度において BH がマダイ瀬

BH, L2,戸内海東

3a, 3hHS より優位（推定値が安定）であった.
BH は変曲点が一意に求まらないため除外．RI と

白川ほか

AR1(同

HS と AICc の差は小さかった (3a).そのため，リス（2022）

時)

ク非対称性の評価を行い，将来予測におけるリスク

アカガレ

HS, L2,イ日本海

3a, 3b, 3hに大きな違いはないことを確認した．また，％SPRを考慮すると HS のほうが現実的な目標であると判断した(3b).
ホッケ道

HS, L1,北系群

AR03a, 3c, 3i

2018 年～2020 年の最新のデータが加わり，平成 31

森田ほか

年度の評価が更新されたが，選択された再生産関係（2022）

は同じであった．具体的には，AICc が最小の HS で
L1 が選択され(3a)，RI との差は 2 程度とわずかであったが， HS は RI と比べて推定の頑健性が高かった(3c)．また，真の再生産関係が RI だったとしても，漁獲や持続性の点で重大な損失が生じるリスク
は小さいことを簡易 MSE により確認した(3i).
キンメダ

HS, L2,イ太平洋

AR03a, 3b, 3d

同一の推定方法，自己相関の考慮の有無において

亘ほか

は，HS,BH,RI の当てはまりは同程度であった(3a).

（2022）

BH および RI は，推定された親魚量の極大値が過去

系群の親魚量の範囲外に位置し，その生物学的な妥当性を検証する材料は少ないため，便宜的に HS とした(3b, 3d)．加入量の自己相関を考慮する場合のほうが当てはまりはよいが，自己相関係数が非常に高く推定され，生物過程によりもたらされる加入量の自己相関の程度はそれよりもはるかに小さいと考えられるため，自己相関なしを採用．
カタクチ

HS, L2,イワシ瀬

AR03a, 3e, 3h

HS, BH, RI で AICc はほぼ同じ(3a)．その中で，HS河野ほか

では親魚量が少ない場合に RI と BH よりも予防的(2022)

11

FRA-SA2024-ABCWG02-05
戸内海系

な加入量となっていること(3e)から HS を採用．自群

己相関はみとめられなかった(3h)．

6． 1B ルール適用事例集
令和 4 年度研究機関会議（3 系群）
ヒラメ日本海中西部・東シナ海系群(管理基準値：SBlimit→SBmin, Fmsy_proxy→F25%SPR, SBban→0.2×SBlimit)



1A 適用時の再生産関係

RI と HS が候補となったが，両者共に，自己相関係数が非常に高く(0.94-0.97)，かつ直近年の加入の残差が負になっていたため，自己相関を考慮した 1A の将来予測で期待される加入量は，負の自己相関が 10 年以上に渡って継続することによって低い水準となる．そのため，水産庁による目標管理基準値の考え方（現在の環境下において持続的な最大の漁獲量を達成するための資源量）と齟齬が生じる．このような自己相関が生じる理由は生息環境あるいは資源評価精度の問題である可能性が指摘されたが現時点では不明である．
1B 適用理由

上記の理由に加えて，RI と HS のどちらの再生産関係を選ぶかによって管理基準値
案が大きく変わり，再生産関係の仮定に不確実性が大きい状況下では頑健な MSY 管
理基準値を得ることが難しい（4 節：1B ルール適用範囲の例 2①に該当）

Fmsy_proxy の 根

Fmax,F0.1,F25%SPR,F30%SPR,F40%SPR を検討．平衡状態における平均親魚量は，いずれの場合

拠でも SBmax を大きく上回っていた．また，現状の漁獲圧（F2018-2020）は，いずれの生物学的管理基準値よりも高い水準にあった． その中で，想定される加入の有効利用による MSY を達成する YPR 基準である Fmax に相当する F25％SPR を Fmsy の代替値
（Fmsy_proxy）の候補とした．これは，ヒラメの他系群で算出されている Fmsy（25%SPR 相当）と同程度であり，F25％SPR を Fmsy_proxy として管理した場合であっても，10 年間で 1 回でも限界管理基準値案(SBmin)を下回る確率は，βに寄らず 0％であった．このことから F25％
SPR でも少なくとも資源を崩壊させるような漁獲圧ではない

1B で用いた

SBtarget 計算時：直近過去 5 年間の天然由来の加入尾数に対して対数正規分布をあては

加入量の仮定め、その分布に基づく加入を仮定
調整係数(β)計算時：同上

引用文献

増渕ほか（2022）

12

FRA-SA2024-ABCWG02-05
トラフグ日本海・東シナ海・瀬戸内海(管理基準値：SBlimit→SBmin, Fmsy_proxy→F30%SPR, SBban→0)

1A 適用時の再生産関係

近年の天然加入が過去と比較して極端に低く，過去最低を更新し続けているものの，加入の減少要因については不明．
1B 適用理由

現状利用可能なデータでは広い範囲の親魚量に対して十分な精度で加入を予測するための十分なデータが得られていないことに加え，HS と RI で管理基準値に大きな違いが生じることから，再生産関係に基づいた頑健な管理基準値を提案することが困難と判断した．適用条件としては，過去に得られた親魚量または資源量の範囲がSB0 またはそれに準じる値と比べたときに非常に狭い範囲となること（4 節の例 3 に該当），再生産関係のパラメータの推定の不安定性のために頑健な MSY 管理基準値を得ることが出来ない（4 節の例 2①に該当）場合に該当する．

Fmsy_proxy の 根

Fmax,F0.1,F20％SPR,F30％SPR を検討．Fmax(=0.32)での漁獲や Fmax と％SPR が近似する F20％SPR

拠での漁獲では，さらなる資源量低下を引き起こす懸念がある．管理基準値に対応する総資源量で比較した場合，F0.1 での MSY 達成時の資源量を超えた資源量は 2006 年漁期のみであり，現在の経年的な加入尾数の減少下では，現実的な目標としては困難である．一方，F30％SPR（＝0.24）での MSY 達成時の総資源量は過去最大の総資源量よりも少なく，現実的な目標として F30％SPR を選択．

1B で用いた

SBtarget 計算時：2002-2020 年（全年）の天然当歳魚の加入尾数に対して対数正規分布

加入量の仮定をあてはめ、その分布に基づく加入調整係数(β)計算時：2002-2020 年（全年）の天然当歳魚の加入尾数に対して対数正規分布をあてはめ、その平均値を予測値として、過去の観測値との残差を 3 年を 1 ブロ
ックとし、ブロックバックワードサンプリング

引用文献


平井ほか（2022）

トラフグ伊勢・三河(管理基準値：SBlimit→SBmin, Fmsy_proxy→F15%SPR, SBban→0)
1A 適用時の再生産関係

近年の加入の減少が親魚の増加による極端な密度効果なのか，環境変化に起因するものなのかによって管理基準値についての考え方が大きく変わるが，現状ではどちらか
13

FRA-SA2024-ABCWG02-05
を正しいと判断するための十分な科学的知見が得られていない．また，RI によって推定された再生産関係は高いスティープネスを導出する密度効果を仮定しており，過去に観測された高い親魚量において極端に低い加入量を予測するほか，過去最低に近い
親魚量において最大の加入を期待するなど，リスク管理の面からも懸念すべき点がみられた．
1B 適用理由

密度効果が弱く仮定される HS をあてはめた場合の SBmsy は RI のときの推定値と大きく乖離したことから，現状では再生産関係に基づいた頑健な管理基準値を提案することが困難と判断した．適用条件としては，再生産関係のパラメータの推定の不安定性のために頑健な MSY 管理基準値を得ることが出来ない（4 節の例 2①に該当）場合に該当する．また，4 節の例 3（過去に得られた親魚量または資源量の範囲が SB0 またはそれに準じる値と比べたときに非常に狭い範囲となる）にも該当する．

Fmsy_proxy の 根

Fmax,F0.1,F20％SPR,F30％SPR を検討．Fmax（F15%SPR に相当）以上の漁獲圧がかかり親魚資源

拠が極端に少なくなった場合でも加入乱獲の兆候が見られなかったという過去の傾向を考慮し，Fmax で漁獲しても過度な加入乱獲を引き起こす懸念は少ないと考え，Fmax
を Fmsy_proxy として提案．

1B で用いた

SBtarget 計算時：2009-2020 年（低加入の時期）の天然当歳魚の加入尾数に対して対数正加入量の仮定規分布をあてはめ、その分布に基づく加入調整係数(β)計算時：同上

引用文献

真鍋ほか（2022）

引用文献
Beverton, R. J. H. and Holt, S. J. (1957) On the Dynamics of Exploited Fish Populations. Her Majesty’s Stationary
Office, London.
Burnham, K. P. and Anderson, A. D. (2002) Model selection and multimodel inference: a practical informationtheoretic approach. Springer, New York, New York, USA.
Clark, C. W., Charles, A. T., Beddington, J. R. and Mangel, M. (1985) Optimal capacity decisions in a developing
fishery. Marine Resource Economics 2, 25–53.
古市 生、由上龍嗣、上村泰洋、林 晃、井須小羊子、渡部亮介 (2020) 令和 2（2020）年度マイワシ太平
洋系群の管理基準値等に関する研究機関会議報告. FRA-SA2020-BRP01-1, 水産研究・教育機構, 横浜, 50pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_maiwashi_p.pdf
古川誠志郎、加賀敏樹、久保田 洋 (2021) 令和会 議 資 料 .

FRA-SA2021-BRP07-01,

3(2021) 年度ブリの管理基準値等に関する研究機関

水 産 研 究 ・ 教 育 機 構 ,

横 浜 ,

74pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20211208/FRA-SA2021-BRP0701.pdf
林 晃、黒田啓行、安田十也 (2019) 平成 31（2019）年度ゴマサバ東シナ海系群の管理基準値等に関する
研 究 機 関 会 議 報 告 .

水 産 研 究 ・ 教 育 機 構 ,

横 浜 ,

39pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_gomasaba_e_r.pdf
Hilborn, R., and Walters, C. J. (1992) Quantitative fisheries stock assessment: Choice, dynamics, and uncertainty.
14

FRA-SA2024-ABCWG02-05
New York: Chapman and Hall.
日野晴彦、黒田啓行、向草世香、佐々千由紀、国松翔太 (2021) 令和 3(2021) 年度カタクチイワシ対

馬暖流系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP03-2, 水産研究・教育機構,
横 浜 ,

56pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20210906/FRA-

SA2021-BRP03-2.pdf
平井慈恵、片町太輔、真鍋明弘、 (2022) 令和 4 年度（2022）年度トラフグ日本海・東シナ海・瀬戸内海系群の管理基準値等に関する研究機関会議. FRA-SA2022-BRP18-01, 水産研究・教育機構, 横浜, 91pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20221027/FRA-SA2022-BRP18-01.pdf
Hurvich, C. M., and Tsai, C. L. (1989) Regression and time series model selection in small samples. Biometrika. 76(2) 297-307.
Ichinokawa，M.， Okamura，H. and Kurota，H. (2017) The status of Japanese fisheries relative to fisheries aroundthe world. ICES Journal of Marine Science 74， 1277–1287.
飯田真也、吉川 茜、八木佑太、藤原邦浩 (2021) 令和

3 (2021) 年度ソウハチ日本海南西部系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP12-3, 水産研究・教育機構, 横浜, 26pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20211026/FRA-SA2021-BRP123.pdf
井須小羊子、由上龍嗣、中神正康、渡邊千夏子、高橋紀夫、上村泰洋、 古市 生、渡部亮介 (2020) 令和
2（2020）年度マアジ太平洋系群の管理基準値等に関する研究機関会議報告 . FRA-SA2020-BRP01-3,
水産研究・教育機構, 横浜, 46pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/20191/detail_maaji_p.pdf
加賀敏樹、岡本 俊、久保田洋、宮原寿恵、西嶋翔太 (2020) 令和 2（2020）年度スルメイカ冬季発生系
群の管理基準値等に関する研究機関会議報告. FRA-SA2020-BRP04-1, 水産研究・教育機構, 横浜,
81pp,https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/20191/detail_surume_w_20201014.pdf
木下順二、安田十也、渡邊千夏子、上村泰洋 (2021) 令和 3(2021) 年度カタクチイワシ太平洋系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP03-1, 水産研究・教育機構, 横浜, 54pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20210906/FRA-SA2021-BRP031.pdf
河野悌昌、高橋正知、安田十也、渡邊千夏子、渡井幹雄、井元順一、木下順二（2022）令和 4（2022）年
度カタクチイワシ瀬戸内海系群の管理基準等に関する研究機関会議資料．FRA-SA2022-BRP05-01, 水産研究・教育機構，横浜，55pp，
久保田洋、宮原寿恵、西嶋翔太、加賀敏樹、岡本 俊 (2020) 令和 2（2020）年度スルメイカ秋季発生系
群の管理基準値等に関する研究機関会議報告. FRA-SA2020-BRP04-2, 水産研究・教育機構, 横浜,87pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_surume-a_20201216.pdf
Mangel，M.，MacCall，A. D.， Brodziak，J.，Dick，E. J.，Forrest，R. E.，Pourzand，R.， & Ralston，S.
(2013). A perspective on steepness，reference points，and stock assessment. Canadian Journal of Fisheries andAquatic Sciences，70(6)，930–940.
15

FRA-SA2024-ABCWG02-05
増渕隆仁、下瀬 環、井関智明 (2022) 令和 4 年度（2022）年度ヒラメ日本海中西部・東シナ海系群の管理基準値等に関する研究機関会議資料. FRA-SA2022-BRP16-01, 水産研究・教育機構, 横浜, 130pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20220930/FRA-SA2022-BRP16-01.pdf
真鍋明弘、平井慈恵、片町太輔、西嶋翔太、澤山周平、青木一弘 (2022) 令和 4 年度（2022）年度トラフ
グ伊勢・三河系群の管理基準値等に関する研究機関会議資料. FRA-SA2022-BRP17-01, 水産研究・教育機構, 横浜, 115pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20221028/FRASA2022-BRP17-01.pdf
三澤 遼、成松庸二、鈴木勇人、森川英祐、時岡 駿、金森由妃、冨樫博幸、 永尾次郎 (2021) 令和

3(2021) 年度ヤナギムシガレイ太平洋北部の管理基準値等に関する研究機関会議資料. FRA-SA2021BRP06-1,

水

産

研

究

・

教

育

機

構

,

横

浜

,

38pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20220825/FRA-SA2022-BRP06-01.pdf
森川英祐、成松庸二、柴田泰宙、鈴木勇人、時岡駿、永尾次郎 (2020) 令和 2（2020）年度ズワイガニ太
平洋北部系群の管理基準値等に関する研究機関会議報告. FRA-SA2020-BRP02-5, 水産研究・教育機
構

,

横

浜

,

27pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-

1/detail_zuwai_pacific_north_r.pdf
森田晶子、境 磨、千葉悟、濱津友紀、山下夕帆、市野川桃子、岡村寛 (2022) 令和 4 年度（2022）年度
ホッケ道北系群の管理基準値等に関する研究機関会議. FRA-SA2022-BRP02-01, 水産研究・教育機構,
横 浜 ,

64pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20220715/FRA-

SA2022-BRP02-01.pdf
森田晶子、山下夕帆、境 磨、磯野岳臣、服部薫、市野川桃子 (2019) 平成 31（2019）年度ホッケ道北系
群 の 管 理 基 準 値 等 に 関 す る 研 究 機 関 会 議 報 告 . 水 産 研 究 ・ 教 育 機 構 , 横 浜 , 50pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_hokke_h.pdf
成松庸二・鈴木勇人・森川英祐・時岡 駿・三澤 遼・金森由妃・冨樫博幸・ 柴田泰宙 (2021) 令和

3（2021）年度マダラ本州太平洋北部系群の管理基準値等に関する研究機関会議資料. FRA-SA2021BRP02-1,

水

産

研

究

・

教

育

機

構

,

横

浜

,

50pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20210824/FRA-SA2021-BRP02-1.pdf
西嶋翔太、由上龍嗣、井須小羊子、上村泰洋、古市生 (2019) 平成 31（2019）年度マサバ太平洋系群の管理基準値等に関する研究機関会議報告.

水産研究・教育機構,

横浜,

21pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_masaba_p.pdf
Punt，A.E.，Smith，A.D.M.，Smith，D.C.，Tuck，G.N.，and Klaer，N.L. (2014). Selecting relative abundanceproxies for BMSY and BMEY. ICES Journal of Marine Science, 71(3)，469–483.
Ricker，W. E. (1954) Stock and recruitment. Journal of Fisheries Research Board of Canada 11， 559–623.
境 磨、山下夕帆、千村昌之、石野光弘 (2019) 平成 31（2019）年度スケトウダラ太平洋系群の管理基準
値 等 に 関 す る 研 究 機 関 会 議 報 告 .

水 産 研 究 ・ 教 育 機 構 ,

横 浜 ,

46pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_suketou_p.pdf
佐久間啓、藤原邦浩、八木佑太、吉川 茜、飯田真也、白川北斗 (2021) 令和

3 （2021）年度マダラ本

州日本海北部系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP-12-1, 水産研究・教育

16

FRA-SA2024-ABCWG02-05
機構, 横浜, 26pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20211026/FRASA2021-BRP-12-1.pdf
下瀬 環、増渕隆仁、中川雅弘 (2021) 令和

3 （2021）年度マダイ日本海西部・東シナ海系群の管理基

準値等に関する研究機関会議資料. FRA-SA2021-BRP05-001, 水産研究・教育機構, 横浜, 39pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20211124/FRA-SA2021-BRP05001.pdf
白川北斗、内藤大河、八木佑太、吉川 茜、佐久間啓、藤原邦浩 (2022) 令和 4 年度(2022)年度アカガレ
イ日本海系群の管理基準値等に関する研究機関会議資料. FRA-SA2022-BRP13-01, 水産研究・教育機構,
横浜, 40pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20221031/FRA-SA2022BRP13-01.pdf
水産研究・教育機構（2024）令和 6 年度(2024)年度漁獲管理規則および ABC 算定のための基本指針．FRASA2024-ABCWG02-01.
高橋 素光、黒田 啓行、依田 真里、向 草世香、安田 十也 (2020) 令和 2（2020）年度マイワシ対馬暖流系群の管理基準値等に関する研究機関会議報告. FRA-SA2020-BRP01-2, 水産研究・教育機構, 横浜,
64pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_maiwashi_tc.pdf
冨樫博幸、成松庸二、鈴木勇人、森川英祐、時岡 駿、三澤 遼、金森由妃、永尾次郎、櫻井慎大 (2022)
令和 4 年度(2022)年度ヒラメ太平洋北部系群の管理基準値等に関する研究機関会議資料. FRA-SA2022BRP07-1,

水

産

研

究

・

教

育

機

構

,

横

浜

,

45pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20220930/FRA-SA2022-BRP07-01.pdf
Thorson，J. T. (2020) Predicting recruitment density dependence and intrinsic growth rate for all fishes worldwide
using a data-integrated life-history model. Fish and Fisheries， 21: 237-251. https://doi.org/10.1111/faf.12427
上田祐司、佐久間啓、藤原邦浩、八木佑太、吉川 茜、松倉隆一、山本岳男 (2020) 令和 2（2020）年度
ズワイガニ日本海系群 A 海域の管理基準値等に関する研究機関会議報告. FRA-SA2020-BRP02-01, 水産研究・教育機構, 横浜, 29pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/20191/detail_zuwai_jsba.pdf
Walters，C. J. and Martell，S. J. D. (2004) Fisheries Ecology and Management. Princeton University Press，Princeton(New Jersey)．
亘真吾、川内洋平、青木一弘、竹村紫苑、竹茂愛吾、半沢祐大（2022）令和 4（2022）年度キンメダイ太
平洋系群の管理基準値等に関する研究機関会議資料．FRA-SA2022-BRP04-01, 水産研究・教育機構，
横浜，48pp,
八木佑太、藤原邦浩、飯田真也、佐久間啓、吉川 茜、白川北斗 (2021) 令和

3 （2021）年度ムシガレ

イ日本海南西部系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP12-2, 水産研究・教育機構, 横浜, 22pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20211026/FRASA2021-BRP12-2.pdf
八木佑太、藤原邦浩、飯田真也、白川北斗 (2022) 令和 4 年度(2022)年度ヒラメ日本海北部系群の管理基準値等に関する研究機関会議資料. FRA-SA2022-BRP07-02, 水産研究・教育機構, 横浜, 42pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20220930/FRA-SA2022-BRP07-02.pdf

17

FRA-SA2024-ABCWG02-05
山本圭介、金谷彩友美、片町太輔、山下夕帆 (2022) 令和 4 年度(2022)年度マダイ瀬戸内海東部系群の管理基準値等に関する研究機関会議資料. FRA-SA2022-BRP14-02, 水産研究・教育機構, 横浜, 54pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20221027/FRA-SA2022-BRP14-02.pdf
山本圭介、片町太輔、山下夕帆、鈴木重則 (2021) 令和 3 （2021）年度マダイ瀬戸内海中・西部系群の管理基準値案等に関する研究機関会議資料. FRA-SA2021-BRP04-002, 水産研究・教育機構, 横浜, 50pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20210928/FRA-SA2021-BRP04002.pdf
山下夕帆、境 磨、千村昌之、石野光弘 (2019) 平成 31（2019）年度スケトウダラ日本海系群の管理基準
値 等 に 関 す る 研 究 機 関 会 議 報 告 .

水 産 研 究 ・ 教 育 機 構 ,

横 浜 ,

31pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_suketou_n.pdf
山下夕帆、山田徹生、真鍋明弘、金谷彩友美、阪地英男 (2021) 令和 3(2021) 年度ヒラメ瀬戸内海系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP04-001, 水産研究・教育機構, 横浜,40pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20210928/FRA-SA2021-

BRP04-001.pdf
安田十也、片町太輔、河野悌昌、高橋正知、渡邊千夏子、渡井幹夫、木下順二、井元順一 (2022) 令和 4
年度(2022)年度サワラ瀬戸内海系群の管理基準値等に関する研究機関会議資料. FRA-SA2022-BRP12-02,
水産研究・教育機構, 横浜, 40pp, https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/20191/20221026/FRA-SA2022-BRP12-02.pdf
安田十也、黒田啓行、林 晃 (2019) 平成 31（2019）年度マサバ対馬暖流系群の管理基準値等に関する研究 機 関 会 議 報 告 .

水 産 研 究 ・ 教 育 機 構 ,

横 浜 ,

19pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_masaba_t.pdf
依田真里、黒田啓行、向草世香、佐々千由紀 (2021) 令和 3(2021) 年度ウルメイワシ対馬暖流系群の管理基準値等に関する研究機関会議資料. FRA-SA2021-BRP03-4, 水産研究・教育機構, 横浜, 43pp,
https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/20210906/FRA-SA2021-BRP03-4.pdf
依田真里、黒田啓行、髙橋素光 (2020) 令和 2（2020）年度マアジ対馬暖流系群の管理基準値等に関する研究機関会議報告.

FRA-SA2020-BRP01-4,

水産研究・教育機構,

横浜,

28pp,

https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_maaji_tc.pdf
由上龍嗣、井須小羊子、上村泰洋、古市 生 (2019) 平成 31（2019）年度ゴマサバ太平洋系群の管理基準
値 等 に 関 す る 研 究 機 関 会 議 報 告.

水 産 研 究 ・ 教 育 機 構 ,横 浜 ,https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/detail_gomasaba_p.pdf

18

30pp,

FRA-SA2024-ABCWG02-05

図 1． モデル診断により推定パラメータの不確実性が明らかになった例（「平成 31（2019）年度ゴ
マサバ東シナ海系群の管理基準値等に関する研究機関会議報告」より）．再生産関係を HS とし，最適化手法に最小絶対値法を用いた時の尤度プロファイル（左）とジャックナイフ解析の結果（右）．
尤度プロファイルでは 2 つの局所解があることが確認でき，ジャックナイフ解析の結果からはデー
タの取り除きによって推定パラメータが大きく変化することが示されている．

図 2． 加入尾数の期待値が過去最大加入尾数以上の極端な外挿値になるために BH や RI の使
用は勧められないケース（Ichinokawa et al． 2017， Figure S3）
19

FRA-SA2024-ABCWG02-05

図 3． 極端な外挿がおこるために RI（左），BH（右）の使用は勧められないケース（Ichinokawa
et al．2017，Figure S3）

図 4． ゴマサバ東シナ海系群における 6 つの再生産関係のあてはめの結果（左）と，推定された管理基準値（右）．右図では，再生産関係式・最適化法ごとに推定された管理基準値の値をプロットし，丸の大きさは Akaike Weight に比例させた（最大の Akaike Weight は BH-L2 で 0．289，最小は HS-L2
で 0．07）．HS の L1 と L2 間での AICc の差は大きくないが，HS-L1 から推定される管理基準値は，より AICc が小さい他の再生産モデルから推定される管理基準値に近く，再生産曲線の選択に対して管理基準値が頑健と考えられる．

20

FRA-SA2024-ABCWG02-05

図 5． マアジ太平洋系群で実施されたモデル平均の例（FRA-SA2020-BRP01-08）．同時推定・最小二乗法において RI（左図，緑点線）の AICc が 27．2，BH（右図，赤点線）が 29．0 と，AIC の差はわずかだが，Fcurrent に対する Fmsy の比が大きく異なる（右図）．右図の横軸は Fmsy/Fcurrent．モデル平均をした結果，右図の青線のような余剰生産曲線が推定され，その余剰生産曲線をもとに管理基準値が決定された．

図 6． ゴマサバ東シナ海系群における感度分析結果．(左) 真の再生産関係が HS-L1（緑），HS-L2（青）であったときに HS-L1 から計算された管理基準値 (hcr1) を使って管理したときの将来予測．(右) 真の再生産関係が HS-L1（緑），HS-L2（青）であったときに HS-L2 から計算された管理基準値 (hcr2)を使って管理したときの将来予測．HS-L1 だと思って管理したが本当は HS-L2 だった場合（左図，青）よりも，HS-L2 だと思って管理したが本当は HS-L1 だった場合（右図，青）の方が，親魚量が過去最低親魚量を下回るリスクが高い．

21

FRA-SA2024-ABCWG02-05

図 7． 将来予測（上）と漁獲量曲線（下）で自己相関を考慮する場合と考慮しない場合の比較．マ
サバ太平洋系群の時系列データを例として用いた．同じ漁獲圧で漁獲しても，自己相関がある場合には将来の親魚量の信頼区間が広く（上），漁獲量曲線（横軸：平衡状態における親魚量の平均値，
縦軸：平衡状態における漁獲量の平均値）は全体的に小さくなり，管理基準値が若干保守的になる．

図 8． マアジ太平洋系群で AR1， L2， RI で推定された再生産関係と観測された加入尾数（灰色），
また，観測値から，再生産関係で予測された自己相関分を引いたときの加入尾数（赤色）．灰色の点は再生産関係にほとんどフィットしていないように見えるが，モデル内で考慮されている自己相関
分を挽くと，再生産関係の周りに分布するようになる．マアジ太平洋系群研究機関会議資料補足図
1-3b より．

22

FRA-SA2024-ABCWG02-05

図 9． ホッケ道北系群で実施した簡易的な MSE の例．上図：真の再生産関係と ABC 算定に用いる再生産関係が同じ（再生産関係の間違いなし）だが，赤（0． 標準再生産関係(MSE)）のケ
ースでは，ABC 算定年までに加入が確率的に変動するとしたもとで，ABC 算定年の漁獲量はABC どおりに漁獲すると仮定している．緑（0． 標準再生産関係(non-MSE)）は通常の将来予測の結果．下図：真の再生産関係が ABC 算定に用いる再生産関係とは異なっていた場合の MSEの結果．赤：再生産関係の間違いなし，緑：より悲観的な再生産関係が真だった場合の MSE の結果．

23

FRA-SA2024-ABCWG02-05

図 10． 過去最低親魚量付近での加入尾数の予測値が楽観的になるという理由で BH，RI の使用は避け，HS を選択したマサバ対馬暖流系群の例（「平成 31（2019）年度マサバ対馬暖流系群の管理基準値等に関する研究機関会議報告」より）．（このケースでは，HS の AICc の値は最適モデルよりも 5．9 大きかっ

たが，HS が選択された．予測力を犠牲にしてこの程度の外挿を避けるべきかは，判断が難しい．なお，この系群では RI や BH だと限界管理基準値が過去最低親魚量以下になること等も考慮してい
る．
）

24

補足資料１

再⽣産関係の診断の⼿順例
最初にfrasyr_toolの流れに沿って複数の再⽣産関係候補に対して⼀度にモデル診断を⾏うスクリプトを紹介します（後述の[#診断内容説明]で項⽬の⼀つ⼀つを説明します）
再⽣産関係の診断の⼿順の⼀例をマイワシ対⾺暖流系群を使って紹介します
診断の⼿法としては収束しているかや⼤域解かのチェック
残差の確率分布のチェック
残差の⾃⼰相関のチェック
スティープネスブートストラップ
ジャックナイフ
プロファイル尤度
パラメータ間相関の8種類を実装しています
通常の再⽣産関係にも、レジーム分けをした再⽣産関係にも適⽤できます

事前準備
frasyrのインストール・読み込み# devtools::load_all('~/git/frasyr') otherwise:
devtools::install_github("ichimomo/frasyr@dev")

Skipping install of 'frasyr' from a github remote, the SHA1 (c704b83e) has notchanged since last install.
Use `force = TRUE` to force installation

library(frasyr)library(tidyverse)

VPAデータの読み込み# setwd('~/git/frasyr_tool/SRcheck') # SRcheckのフォルダがある場所を指定
res_vpa <- get(load("../data_SC2020/res_vpa_TMI2019assess.rda"))フィットに使う設定を1do_MSYest.rとあわせる
1get.SRdataの引数weight.yearを1do_MSYest.rの設定(year_for_srfit)と同じにします

year_for_srfit <- 0
data_SR = get.SRdata(res_vpa, weight.year = year_for_srfit)

steepnessを計算するときの⽣物パラメータ設定
derive_biopar関数を使い、引数にVPAの結果オブジェクトか、将来予測の結果オブジェクトを⼊れます
VPA結果の単純な平均を使う場合は以下のような感じになります

bio_par1517 <- derive_biopar(res_vpa, derive_year = 2015:2017)

密度によって体重や成熟率が変化する設定を⽤いる場合には、将来予測で平衡状態に達したときの⽣物パラ
メータの平均をとってくることになります# res_future_0.8HCRが密度によって体重などが変化する将来予測結果の場合
bio_par1517 <- derive_biopar(res_future_0.8HCR, derive_year = 2050)

再⽣産関係の候補と名前を決定して診断(全体の流れ)
複数の再⽣産関係の候補モデルについてまとめてモデル診断を⾏うための例コードを⽰します
再⽣産関係モデルごとにフォルダを⽣成してそこに結果を出⼒します
fit.SRやfit.SRregimeのかわりにfit.SR_tolを使うことで、解の収束などを同時に診断しています

SRname_list = c("HSL2outer", "HSL2inner", "HSL2regime", "HSL1regime")for (i in 1:length(SRname_list)) {cat(SRname_list[i], "\n")SRname = SRname_list[i]dir.create(SRname)# 繰り返し計算の設定。時間短縮のため、ここの例コードでは少なめに設定しています
n_check <- 10 # 100 # 初期値を変えて繰り返し計算するときの回数。100以上が良い。
n_boot <- 10 # 200 # ブートストラップするときの回数。200以上が良い。
n_grid <- 10 # 50 # 尤度プロファイルするときのgridの数の細かさ。デフォルトは50。
# 各モデルの設定をここで変えます (SRname_listの1番⽬から順番に)# bio_parを与えるとh関連のパラメータも計算されます# fit.SRやfit.SRregimeのかわりにfit.SR_tolを使うと、パラメータ推定時に初期値を変えた計算を100回おこない、# ⼤域解に達しているか・同じ尤度で異なるパラメータ範囲があるかなどを確認します。
# ⼤域解に達していない場合、⼤域解と思われる初期値で計算したパラメータで置き換えます。
# 同じ尤度で異なるパラメータ範囲がある場合、中央値をとります
resSR = switch(i, fit.SR_tol(data_SR, SR = "HS", AR = 1, out.AR = TRUE, length2

= 50,
bio_par = bio_par1517, plus_group = res_vpa$input$plus.group, n_check =n_check),fit.SR_tol(data_SR, SR = "HS", AR = 1, out.AR = FALSE, length = 50,
bio_par = bio_par1517,
plus_group = res_vpa$input$plus.group, n_check = n_check),
fit.SR_tol(data_SR,
SR = "HS", method = "L2", regime.year = c(1976, 1988), regime.par =c("a","b", "sd")[1:3], regime.key = c(0, 1, 0), length = 10, use.fit.SR= FALSE,
bio_par = bio_par1517, plus_group = res_vpa$input$plus.group,
is_regime = TRUE,
n_check = n_check), fit.SR_tol(data_SR, SR = "HS", method = "L1",
regime.year = c(1976,1988), regime.par = c("a", "b", "sd")[1:3], regime.key = c(0, 1, 0),
length = 10, use.fit.SR = FALSE, bio_par = bio_par1517, plus_group =res_vpa$input$plus.group,
is_regime = TRUE, n_check = n_check))# 収束＆⼤域解かどうかを診断# fit.SR_tolを使う場合、以下のプロセスはfit.SR_tolに包含されるため使⽤しなくて良いです# check = check.SRfit(resSR,output=TRUE,filename =# paste0(SRname,'/checkSRfit_',SRname)) if (!is.null(check$optimum)) {# #⼤域解に達していない場合、結果を更新して再チェック resSR = check$optimum check# = check.SRfit(resSR,output=TRUE,filename =# paste0(SRname,'/checkSRfit_',SRname)) }# 結果をtxtファイルに出⼒
out.SR(resSR, filename = paste0(SRname, "/resSR_", SRname))# 再⽣産関係の図を出⼒if (class(resSR) == "fit.SR") {g1 = SRplot_gg(resSR)
ggsave_SH(g1, file = paste0(SRname, "/SRcurve_", SRname, ".png"))} else {g1 = SRregime_plot(resSR, labeling.year = c(min(data_SR$year),
data_SR$year[data_SR$year %in%seq(1900, 2100, by = 10)], max(data_SR$year)))
ggplot2::ggsave(g1, file = paste0(SRname, "/SRcurve_", SRname, ".png"),
unit = "mm",width = 240, height = 120, dpi = 600)}
# プロファイル尤度を計算 a_range,# b_rangeはデフォルトではNULL（関数内で適当に設定）ですが、⼗分な範囲出ない場合にはa_rangeとb_rangeを適当な範囲に変更してください
proflik = prof.likSR(resSR, output = TRUE, a_range = NULL, b_range = NULL,
filename = paste0(SRname,"/ProfLik_", SRname), length = n_grid)# 残差の分布チェック
check.SRdist(resSR, output = TRUE, filename = paste0(SRname, "/ResidDist_",SRname))# 残差(deviance)の⾃⼰相関をチェック
3autocor.plot(resSR, output = TRUE, filename = paste0(SRname, "/devianceSR_",
SRname))# ⾃⼰相関を推定した場合、⾃⼰相関の推定後の残差(residual)に対してもチェックを⾏う# devianceのプロットでは⾃⼰相関があるが、residualのプロットではだいたい⾃⼰相関が消えていることを確認する
if (class(resSR) == "fit.SR" && resSR$input$AR) {autocor.plot(resSR, use.resid = 2, output = TRUE, filename =paste0(SRname,"/residualSR_", SRname))}# ジャックナイフ解析
jack1 = jackknife.SR(resSR, output = TRUE, is.plot = TRUE, filename =paste0(SRname,"/jackknife_", SRname))# ブートストラップ解析 パラメトリック（⾃⼰相関ある場合はパラメトリックを推奨）
boot.res1 = boot.SR(resSR, n = n_boot, method = "p")bootSR.plot(boot.res1, filename = paste0(SRname, "/para_boot_", SRname),
output = TRUE)# ノンパラメトリック
boot.res2 = boot.SR(resSR, n = n_boot, method = "n")bootSR.plot(boot.res2, filename = paste0(SRname, "/nonpara_boot_", SRname),
output = TRUE)# 結果の保存SRdiagnosis = list(resSR = resSR, ProfLik = proflik, ParaBoot = boot.res1,
NonParaBoot = boot.res2,
Jackknife = jack1)
save(SRdiagnosis, file = paste0(SRname, "/SRdiagnosis.rda"))# resSRだけも単独に保存
save(resSR, file = paste0(SRname, "/resSR.rda"))}HSL2outer...check.SRfit 1 回⽬....
1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. パラメータが唯⼀の解として推定されているのでOKです (OK)HSL2inner...check.SRfit 1 回⽬....
1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. パラメータが唯⼀の解として推定されているのでOKです (OK)HSL2regime...check.SRfit 1 回⽬....
4

1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. パラメータが唯⼀の解として推定されているのでOKです (OK)Joining, by = c("Regime", "SSB", "R", "Category")
Standardized residuals are used for 'fit.SRregime'
HSL1regime...check.SRfit 1 回⽬....
1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. 同じ最⼤尤度(1e-06よりも⼩さい違い)を持つ複数のパラメータが⾒つかりました（L1かつHSでよく⾒られます）。
同じ尤度を持つパラメータの範囲 ( 10 回試⾏),
regime2
a1
a2
b1
b2
sd1
sd2
Min.
1 0.02508404 0.04534170 504707.9 2336658 0.6896823 0.4273242
1st Qu.
1 0.02525135 0.04534170 544176.0 2336658 0.6899769 0.4273242
Median1 0.02595267 0.04534171 557915.5 2336658 0.6915280 0.4273242
Mean
1 0.02624902 0.04534171 552561.7 2336658 0.6917446 0.4273242
3rd Qu.
1 0.02659644 0.04534171 573167.8 2336659 0.6934601 0.4273243
Max.
1 0.02867613 0.04534171 576983.3 2336659 0.6941651 0.4273243
中央値に最も近いパラメータセットを持つ推定結果をoptimumに出⼒します(そのときの初期値は-1.65729209742478--3.36255325827061--9.0518671379257--1.69383836850377 です)尤度が同じパラメータの範囲で中央値をとる結果に置き換えます
Joining, by = c("Regime", "SSB", "R", "Category")
Standardized residuals are used for 'fit.SRregime'

ここまでが複数の再⽣産候補を⼀度に診断するスクリプトです。個々の診断内容の説明は以下をご参照ください
レジームなしの場合 →補⾜資料２

レジームありの場合 →

補⾜資料３

5

補足資料２

個々の診断内容の説明（レジームなしの場合）
準備
frasyrのインストール・読み込み、vpaオブジェクトの読み込み

devtools::install_github("ichimomo/frasyr@dev")

Skipping install of 'frasyr' from a github remote, the SHA1 (c704b83e) has notchanged since last install.
Use `force = TRUE` to force installation

library(frasyr)library(tidyverse)

Registered S3 method overwritten by 'httr':
methodfromprint.response rmutil
─ Attaching packages ───────────────────────────────────── tidyverse 1.3.0 ─✔ ggplot2 3.3.5✔ tibble 3.1.3✔ tidyr
1.1.3✔ readr
2.0.0✔ purrr
0.3.4✔ dplyr
1.0.7✔ stringr 1.4.0✔ forcats 0.5.1

─ Conflicts ─────────────────────────────────────── tidyverse_conflicts() ─✖ dplyr::filter() masks stats::filter()✖ dplyr::lag()masks stats::lag()

res_vpa <- get(load("../data_SC2020/res_vpa_TMI2019assess.rda"))フィットに使うSRdata整形
vpaの結果オブジェクトを準備したらget.SRdata関数の引数に⼊れてfit.SR関数に必要なSRdataを作成します
get.SRdataではweight.yearで再⽣産関係を推定するために使う（尤度の計算に含める）年のデータを指定
できます

6

0を指定するとvpaの結果オブジェクトに含まれる全年のデータが使われます
1do_MSYest.rではyear_for_srfitで設定しているので、同じように設定します

year_for_srfit <- 0
SRdata = get.SRdata(res_vpa, weight.year = 0)

⽣物パラメータの設定
再⽣産関係の推定をするfit.SR関数にある引数bio_parに⽣物パラメータをいれるため、まずderive_biopar関数でbio_parオブジェクトを作成します
これは後述のスティープネスを算出するために必須です
vpaの結果のオブジェクトか、将来予測の結果のオブジェクトを⼀つ⽬の引数に⼊れます
ここでは⽣物パラメータとしてres_vpaに収められている⼊⼒の最終年をつかいますh_year <- max(res_vpa$input$rec.year)
bio_par <- derive_biopar(res_vpa, derive_year = h_year)

上述の診断⼿順例に記述されているように、密度によって体重や成熟率が変化する場合、将来予測で平衡状態に達したときの⽣物パラメータの平均を使います
そのときの⼀つ⽬の引数は将来予測の結果のオブジェクトになります

通常の再⽣産関係（レジームなし）
まず、レジームを考慮しない通常の再⽣産関係についての診断について説明します

再⽣産関係の推定とプロット
ここではHS再⽣産関係の①L2で⾃⼰相関ナシ、②L1で⾃⼰相関あり、③L1で⾃⼰相関を外側推定、④L2で⾃⼰相関を内側推定の4つのモデルを検討します

resL2 = fit.SR(SRdata, SR = "HS", method = "L2", out.AR = FALSE, AR = 0, bio_par =bio_par)resL1 = fit.SR(SRdata, SR = "HS", method = "L1", out.AR = FALSE, AR = 0, bio_par =bio_par)resL1outer = fit.SR(SRdata, SR = "HS", method = "L1", out.AR = TRUE, AR = 1,
bio_par = bio_par)resL2inner = fit.SR(SRdata, SR = "HS", method = "L2", out.AR = FALSE, AR = 1,
bio_par = bio_par)

結果の出⼒と図⽰
7

out.SRで結果をtxtファイルに出⼒できます
例として④を使⽤します

out.SR(resL2inner, filename = "L2inner")

プロットにはSRplot_ggが使えます(g1 = SRplot_gg(resL2inner))ggsave_SH(g1, file = "SRcurve.png")

推定値の収束の有無や最適解に達しているかをチェック
check.SRfitという関数を使って、適切な推定値が得られているかを診断できます
fit.SRやfit.SRregimeのかわりにfit.SR_tol(is_regime=FALSEでレジームなし、is_regime=TRUEでレジームあり)でパラメータを推定すると、以下のモデル診断と⼤域解に達していない場合に解を置き換える作業を⾃動でおこないます
診断の結果(⾚字でメッセージが表⽰された場合、問題がある可能性があります)①収束しているか (convergence)②ヘッセ⾏列 (Hessian) の対⾓成分が正定値 (positive definite) になっているか③境界条件（HSの場合SSBの両端、それ以外の場合直線関係や極端に強い密度依存関係）に達しているか④対数尤度が⼤域解 (global optimum) に 達しているか⑤⼤域的最適解をもつパラメータがひとつに決まるかどうか（④で⼤域解が出られている場合
のみ）
8

HSでL1を使う場合、⼤域解にたどり着いていなかったり、尤度が全く同じ解が複数得られることがあります
output = TRUEで結果を出⼒

check_L2 = check.SRfit(resL2, output = TRUE, filename = "check_L2")

1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. パラメータが唯⼀の解として推定されているのでOKです (OK)check_L1 = check.SRfit(resL1, output = TRUE, filename = "check_L1")

1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっていません (HSかつ3以降のチェックがOKであれば問題ありません)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)Maximum difference of log-likelihood is 0.00371
4. パラメータが⼤域解に達していません (fit.SRtolで再計算をおこなうか、本関数の返り値のoptimumに結果を置き換えてください。⼤域解を得るための初期値は,-0.627443985964221, -0.851863010270145です)

最⼩絶対値法 (L1)のときに’4. パラメータが⼤域解に達していません…’というメッセージが出ます。これは、初期値を変えて推定を繰り返したとき、resL1を計算したときに与えた初期値から推定されたパラメータのもとでの対数尤度よりも⼤きい対数尤度を持つパラメータが推定されたという意味です。
このメッセージが出た場合、check_L1$optimumに最適解の結果のオブジェクトが格納されています
ので、以下のように、対数尤度やAICcやパラメータ値がどのくらい違っているのかを⽐較できますc(resL1$loglik, check_L1$optimum$loglik)[1] -69.26809 -69.26438

c(resL1$AICc, check_L1$optimum$AICc)
9

[1] 144.9725 144.9651

rbind(resL1$pars, check_L1$optimum$pars)ab
sd rho1 0.03092958 2331717 0.7699262
0
2 0.03092876 2573979 0.7791758
0check_L1で得られた最適解でresL1を上書きして、4. の項⽬が(OK)になるか再度チェックします

resL1 = check_L1$optimum
check_L1_rev = check.SRfit(resL1, output = TRUE, filename = "check_L1")

1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. 同じ最⼤尤度(1e-06よりも⼩さい違い)を持つ複数のパラメータが⾒つかりました（L1かつHSでよく⾒られます）。
同じ尤度を持つパラメータの範囲 ( 100 回試⾏),
a
b
sd
Min.
0.03018238 2573910 0.7772454
1st Qu. 0.03036667 2583689 0.7776661
Median 0.03062992 2599092 0.7783313
Mean
0.03060973 2600950 0.7783052
3rd Qu. 0.03081253 2621623 0.7788363
Max.
0.03092958 2637631 0.7791782
中央値に最も近いパラメータセットを持つ推定結果をoptimumに出⼒します(そのときの初期値は-10.3982031433377-9.35999964262816 です)

その結果4.の項⽬はOKになりましたが、こんどは’5.同じ最⼤尤度を持つ複数のパラメータが⾒つかりました（L1かつHSでよく⾒られます）’というメッセージが出ます。これは、ほとんど同じ尤度だけど異なる複数のパラメータがあるという警告です。この例の場合には、b=2573910から2637631の範囲
でほぼ尤度が変わらないパラメータセットが出現していることを⽰しています。このようなケースの場合には、このパラメータの範囲のmedianに最も近いパラメータセットを持つ結果が$optimumに返
されます

10

どのくらい尤度やパラメータが異なっているかを確認してから、’check_L1_rev$optimum’の結果を
resL1に⼊れます。

# どのくらい尤度やパラメータが違うか？
c(resL1$loglik, check_L1_rev$optimum$loglik)[1] -69.26438 -69.26438

c(resL1$AICc, check_L1_rev$optimum$AICc)[1] 144.9651 144.9651

rbind(resL1$pars, check_L1_rev$optimum$pars)ab
sd rho1 0.03092876 2573979 0.7791758
0
2 0.03062992 2599092 0.7783313
0# check_L1_rev$optimをresL1に上書き
resL1 = check_L1_rev$optimum

上記のプロセスを⼀括して⾏う関数としてfit.SR_tolを⽤意しています。以下の例では、...check.SRfit 1 回⽬....で初期値を⼤域解に近いものに置き換えてます。...check.SRfit
2 回⽬....で同じ尤度を持つパラメータのうち中央値のものを選んで返しています

resL1_tol = fit.SR_tol(SRdata, SR = "HS", method = "L1", out.AR = FALSE, AR = 0,
bio_par = bio_par)...check.SRfit 1 回⽬....
1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっていません (HSかつ3以降のチェックがOKであれば問題ありません)113. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)Maximum difference of log-likelihood is 0.00371
4. パラメータが⼤域解に達していません (fit.SRtolで再計算をおこなうか、本関数の返り値のoptimumに結果を置き換えてください。⼤域解を得るための初期値は,-0.627443985964221, -0.851863010270145です)
⼤域解を得るための初期値に変えたres_SRの結果に置き換えます...check.SRfit 2回⽬....
1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. 同じ最⼤尤度(1e-06よりも⼩さい違い)を持つ複数のパラメータが⾒つかりました（L1かつHSでよく⾒られます）。
同じ尤度を持つパラメータの範囲 ( 100 回試⾏),
a
b
sd
Min.
0.03018238 2573910 0.7772454
1st Qu. 0.03036667 2583689 0.7776661
Median 0.03062992 2599092 0.7783313
Mean
0.03060973 2600950 0.7783052
3rd Qu. 0.03081253 2621623 0.7788363
Max.
0.03092958 2637631 0.7791782
中央値に最も近いパラメータセットを持つ推定結果をoptimumに出⼒します(そのときの初期値は-10.3982031433377-9.35999964262816 です)尤度が同じパラメータの範囲で中央値をとる結果に置き換えます

⼿動で計算したresL1とfit.SR_tolを使ったresL1_tolで結果は同じになります
script/1do_MSYest.rのdo_check_SRfitをTRUEにするとfit.SR_tolを使った結果が返されます

rbind(resL1$pars, resL1_tol$pars)fit.SR_tolが包含しているので# resL1とresL1_tolは同じ。resL1を計算するプロセスを

ab
sd rho1 0.03062992 2599092 0.7783313
0
2 0.03062992 2599092 0.77833130残差の分布のチェック
check.SRdistという関数を使って、残差の分布が仮定にあっているかをチェックできます
図が3つ出⼒され、左が標準化された（標準偏差で割られた）残差のヒストグラムで⾚い線が標準正規分布を表します。右上のSWはShapiro-Wilk検定により、残差が正規分布に従っているかを検定したときのP値です。右上のKSはKolmogorov-Smirnov検定により、残差が正規分布に従っているかを検定したときのP値です。どちらも、帰無仮説は「正規分布に従っている」です12真ん中は残差の累積確率密度のヒストグラムを⽰しており、理論的には0~1の⼀様分布に従うので、このヒストグラムが⼤きく凸凹しているようだと仮定した確率分布と合致していないことになります
右も累積確率密度を0~1の⼀様分布と仮定したときのQQプロットで、線は理論値を表しており、この線から外れていないほうがよいということになります。注意すべき点は、このQQプロットは⼀様分布を仮定しているので、端は線に載りやすく、中⼼付近が外れやすいことです。つまり、qqnorm等で出⼒される正規分布を仮定したQQプロットとは単純な⽐較はできません。これも、L1と⽐較するために
こういう処理をしていますcheck.SRdist(resL2)

output=TRUEとすれば、pngファイルが出⼒され、保存されます

check.SRdist(resL2, output = TRUE, filename = "ResidDistCheck_L2")

L1を使う場合も同様の図が描かれますが、左の図に⻘い点線がプロットされます。これはラプラス分
布を表しており、L1の場合の尤度はラプラス分布を使って計算されるため⽰してあります。右上のKS(Kolmogorov-Smirnov検定) のP値もラプラス分布との⼀致性を検定しています。SWの⽅は正規分布にのみ適⽤できる検定なので、正規分布との⼀致性を検定しています。L1推定の場合は正規分布を仮定しているわけではないので、この値は参考程度でいいと思います
真ん中と右の図は、ラプラス分布を仮定したときの累積確率密度のヒストグラムとQQプロットになりますcheck.SRdist(resL1)

13

⾃⼰相関を外側推定した場合にはファイルが2つ出⼒され、1つ⽬が再⽣産関係の推定、2つ⽬が⾃⼰相関の推定における残差です。しがたって、1枚⽬の出⼒は上と同じになります。⾃⼰相関の推定はmethodがL1/L2によらず正規分布を仮定していますcheck.SRdist(resL1outer)

⾃⼰相関を内側で推定した場合には再⽣産関係と⾃⼰相関係数を同時に推定するので、⾃⼰相関を除いた残差に対する図が出⼒されます14check.SRdist(resL2inner)

残差の⾃⼰相関のチェック
残差の⾃⼰相関を事後的に推定するための関数calc.residARを作りました。fit.SRのout.ARオプシ
ョンとやっていることは同じですが、AICの差だけではなく尤度の値やAICc等の値も取り出せます
後述しますが、fit.SRregimeでも外側で⾃⼰相関を推定できます
output=TRUEで結果をテキストファイルに出⼒します ※このチェックは必須ではないですouter1 = calc.residAR(resL1, output = TRUE, filename = "residARouter")rbind(outer1$pars, resL1outer$pars)

a
b
sd
rho
1 0.03062992 2599092 0.6592486 0.5243604
2 0.03092958 2331717 0.6563346 0.5157942

c(diff(outer1$AIC), diff(resL1outer$AIC.ar))AR(1)1-17.27257 -16.52632

outer1$AICc

AR(0)
AR(1)
139.9338 122.8053

15autocor.plotは残差の⾃⼰相関に関するプロットを3つ出⼒します。左は再⽣産関係との残差の時系
列を表しており、⾚い線は平滑化された曲線です。この線が0から外れていると残差がホワイトノイズ
ではなく、トレンドがあると⾔えます
真ん中はラグを1,2,3,…年と増やした場合の⾃⼰相関係数を表しています。⻘い線は95％信頼区間を⽰しています。1年のラグで有意なときは⾃⼰相関の推定を検討したほうがよいと考えられます
右側はLjung-Box検定におけるP値です。Ljung-Box検定の帰無仮説は「ラグが1~mまでの全ての⾃⼰相関が0である」であり、対⽴仮説は「ラグ1から mまでの⾃⼰相関のうち、少なくとも⼀つが0でない」となります (参考: https://to-kei.net/time-series-analysis/hypothesis-testing/ )。この場合、全ての点で有意（⾚点で表⽰される）となっているのは、ラグ1の相関係数が⾼いからだと考えられます

autocor.plot(resL1outer)上は*再⽣産関係との逸脱度(deviance)*についてチェックをしましたが、*⾃⼰相関を推定した後での残差(residual)*についても、use.resid=2とすれば同様の図を描くことができます (ディフォルトは
use.resid=1)。この場合、有意な⾃⼰相関が消えていることが確認できますので、⾃⼰相関モデルを⽤いる
場合にはdevianceでは有意な⾃⼰相関構造があるのに対して、residualには⾃⼰相関構造がなくなっていることを確認してください。

autocor.plot(resL1outer, use.resid = 2)16output=TRUEとすればpngファイルが出⼒されます

autocor.plot(resL1outer, output = TRUE, filename = "devianceAR")
autocor.plot(resL1outer, use.resid = 2, output = TRUE, filename = "residAR")

スティープネス
SR=BH,RIにおいて、F=0の時の初期資源量から理論上求められる加⼊量R0に対する、0.2B0の時の加⼊量Rの⽐であるスティープネスhが計算されます。
スティープネスは再⽣産モデルの密度補償効果の程度を⽰します。0.2から1の間の値を取り、0.2の時は密度補償効果がないことを、1の時は親⿂量に関わらず⼀定の加⼊量が得られることを意味します。
derive_biopar関数で作成したbio_parオブジェクトをfit.SR関数の引数に⼊れてあれば、fit.SR関数の戻り値オブジェクトに格納されています
スティープネスhのみならず、SPR0、SB0、R0、B0も出⼒できます

resL2$steepness$h#スティープネスのみ[1] 0.8010001

resL2$steepness

SPR0
SB0
R0
B0
h
1 166.5378 7477300 44898.51 9412160 0.8010001

ブートストラップ17

boot.SRには3種類のブートストラップ⽅法を実装しており、うち2つが残差のブートストラップで、残り⼀つがデータのブートストラップです。残差ブートストラップには、パラメトリックとノンパラ
メトリックの2つあります。⾃⼰相関を推定している場合には、ノンパラメトリックは対応していない
ので、パラメトリックブートストラップを使⽤するようにしてください。
残差のパラメトリックブートストラップは、fit.SRのpars$sdを標準偏差とする正規分布からランダ
ムに乱数を発⽣させ、予測値からのずれを加えて加⼊量のブートストラップデータを⽣成し、再推定しています
a,bの信頼区間は数値を⾒ただけで狭い・広いの判断がつきにくいので、h(0.2SB0のときの平均加⼊尾
数/R0)やR0(SB0のときの平均加⼊尾数)の信頼区間を⾒て、パラメータの不確実性が⼤きいかどうかを確認してください
残差のパラメトリックブートストラップはmethod="p"で実⾏可能で、図のプロットにはbootSR.plot
を使⽤しますここでは100回でやっていますが、最終的に使⽤する予定の再⽣産関係には数を増やしたほうがよいと思われます（500回とか1000回）
⾃⼰相関を推定していない場合は最後のrhoの図は表⽰されません
fit.SRの引数にbio_parを⼊れてスティープネスを計算した場合、SB0、R0、B0、hの図も出⼒されます
図の出⼒bootSR.plotのオプションoutput=Tとすると、各パラメータのヒストグラムが出⼒されます
図の出⼒bootSR.plotのオプションggplt=Tとすると、ヒストグラムはggplotで描画されます

boot.res1 = boot.SR(resL2inner, n = 100, method = "p")# bootSR.plot(boot.res1)
bootSR.plot(boot.res1, ggplt = T)

Joining, by = "name"`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

18

19

bootSR.plot(boot.res1, output = T, ggplt = T)

Joining, by = "name"`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

20

png
2

残差のノンパラメトリックブートストラップは残差の確率分布を仮定せず、残差を重複ありでリサン
プリングして、加⼊量のブートストラップデータを⽣成します
確率分布を仮定しないので、L1とL2を⽐較する場合にはパラメトリックよりもこちらの⽅がよいかもしれません
ノンパラメトリックブートストラップはmethod="n"で⾏います

boot.res2 = boot.SR(resL2inner, n = 100, method = "n")bootSR.plot(boot.res2, ggplt = T)21

Joining, by = "name"`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

22

データのブートストラップはデータを重複ありでリサンプリングしたデータを使⽤して、再⽣産関係の再推定を⾏います
親⿂量データもリサンプリングにより変化するため、親⿂量の不確実性も考慮されることになります
親⿂量データに偏りがあったり、データ数が少なかったり、あるデータ点に推定値が⼤きく依存している場合はバイアスや不確実性が⼤きくなりやすいと思われます
データブートストラップはmethod="d"で⾏います
残差ブートストラップを基本として、データのブートストラップはオプションで参考程度でよいと思います
この場合データブートストラップの⽅が残差ブートストラップよりもバイアスが⼩さいのは親⿂量データのバイアスは⼩さく、加⼊変動が⼤きいからでしょうかね…

boot.res3 = boot.SR(resL2inner, n = 100, method = "d")bootSR.plot(boot.res3, ggplt = T)
23

Joining, by = "name"`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

24

図を画像ファイルとして保存する場合はoutput=TRUEとします

bootSR.plot(boot.res1, filename = "para_boot", output = TRUE)
bootSR.plot(boot.res2, filename = "nonpara_boot", output = TRUE)# bootSR.plot(boot.res3,filename='data_boot',output=TRUE)

ジャックナイフ
データを⼀点ずつ除いて再推定し、各データが与える影響を評価します
こちらについても、a,bのぶれの⼤きさは判断がつきにくいので、hやR0がどのくらいずれるかに着⽬
してくださいjack1 = jackknife.SR(resL2inner, is.plot = TRUE)
25

26

27

28

図を画像ファイルとして保存する場合︓

jack1 = jackknife.SR(resL2inner, output = TRUE, is.plot = TRUE)

プロファイル尤度
再⽣産関係のパラメータa,bを変化させたときの尤度を計算し、点推定値以外に、尤度が⾼くなっている⼭が無いかをチェックします
⾃⼰相関を内部推定している場合はa,bを固定して尤度が最も⾼くなる⾃⼰相関係数rhoを推定しています
HSで⾃⼰相関を内部推定すると複数の⼭が現れることがあるが、この場合は⼤丈夫そう
a_range, b_rangeで描画するパラメータの範囲を、lengthでgridの細かさを指定できます。デフォルトはa_rangeとb_rangeはNULL(内部で適当に範囲を決定)ですが、できた図を⾒てみていい感じでない場
合ではa_rangeとb_rangeを適当な範囲に変更してください

prof1 = prof.likSR(resL2inner, a_range = NULL, b_range = NULL, length = 50)図を保存する場合:

prof1 = prof.likSR(resL2inner, output = TRUE)

パラメータ間相関
パラメータa, bの値がどれだけ相関を持つか、相関係数を算出します
パラメータa, b間の相関が強いと推定が不安定になっている可能性があります29check.SRfitで収束が安定しなかったり、尤度プロファイルを出⼒した時に推定値が⼭の頂点になっていない場合に相関をチェックしてください
SR関係をプロットして、原点からmin(SSB)まで直線的な関係があると不安定な傾向があります

res_corSR <- corSR(resSR = resL1)
res_corSR$cor[1, 2][1] -0.7009081

30補足資料３

個々の診断内容の説明（レジームあり）
準備
frasyrのインストール・読み込み、vpaオブジェクトの読み込みはうえに同じ

devtools::install_github("ichimomo/frasyr@dev")

Skipping install of 'frasyr' from a github remote, the SHA1 (c704b83e) has notchanged since last install.
Use `force = TRUE` to force installation

library(frasyr)library(tidyverse)

Registered S3 method overwritten by 'httr':
methodfromprint.response rmutil
─ Attaching packages ─────────────────────────────────────────────── tidyverse
1.3.0 ─✔ ggplot2 3.3.5✔ tibble 3.1.3✔ tidyr
1.1.3✔ readr
2.0.0✔ purrr
0.3.4✔ dplyr
1.0.7✔ stringr 1.4.0✔ forcats 0.5.1

─ Conflicts ─────────────────────────────────────────────────
tidyverse_conflicts() ─✖ dplyr::filter() masks stats::filter()✖ dplyr::lag()masks stats::lag()

res_vpa <- get(load("../data_SC2020/res_vpa_TMI2019assess.rda"))フィットに使うSRdata整形
vpaの結果オブジェクトを準備したらget.SRdata関数の引数に⼊れてfit.SR関数に必要なSRdataを作成します31get.SRdataではweight.yearで再⽣産関係を推定するために使う（尤度の計算に含める）年のデータを指定できます0を指定するとvpaの結果オブジェクトに含まれる全年のデータが使われます
1do_MSYest.rではyear_for_srfitで設定しているので、同じように設定します

year_for_srfit <- 0
SRdata = get.SRdata(res_vpa, weight.year = 0)

⽣物パラメータの設定
再⽣産関係の推定をするfit.SR関数にある引数bio_parに⽣物パラメータをいれるため、まずderive_biopar関数でbio_parオブジェクトを作成します
これは後述のスティープネスを算出するために必須です
vpaの結果のオブジェクトか、将来予測の結果のオブジェクトを⼀つ⽬の引数に⼊れます
ここでは⽣物パラメータとしてres_vpaに収められている⼊⼒の最終年をつかいますh_year <- max(res_vpa$input$rec.year)
bio_par <- derive_biopar(res_vpa, derive_year = h_year)

上述の診断⼿順例に記述されているように、密度によって体重や成熟率が変化する場合、将来予測で平衡状態に達したときの⽣物パラメータの平均をつかいます
そのときの⼀つ⽬の引数は将来予測の結果のオブジェクトになります。

再⽣産関係の推定
推定にはfit.SRregimeを使います。ここではHS, L2で、1976年と1988年にレジームシフトが⽣じた場
合を候補としますregime.key = c(0,1,0)というのは、レジームがA->B->Aのように変化することを想定しています
regime.parはレジームごとに異なるパラメータを設定します。今はすべてのパラメータが異なる場合とaは共通でbとsdが異なる場合を使います
AICc等の結果はfit.SRと⽐較できます
詳しくは?fit.SRregimeをご覧ください

resR1 <- fit.SRregime(SRdata, SR = "HS", method = "L2", regime.year = c(1976,1988),
regime.par = c("a", "b", "sd")[1:3], use.fit.SR = TRUE, regime.key = c(0, 1,0), bio_par = bio_par, plus_group = res_vpa$input$plus.group)
resR2 <- fit.SRregime(SRdata, SR = "HS", method = "L2", regime.year = c(1976,1988),
regime.par = c("a", "b", "sd")[2:3], use.fit.SR = TRUE, regime.key = c(0, 1,0), bio_par = bio_par, plus_group = res_vpa$input$plus.group)c(resR1$AICc, resR2$AICc)

32

[1] 120.6456 131.1173

結果の出⼒と図⽰
fit.SRのときと同様にout.SRで結果をtxtファイルに出⼒できます

out.SR(resR1, filename = "Regime1")

プロットにはSRregime_plotが使えます(g2 = SRregime_plot(resR1, regime.name = c("Usual", "High")))Joining, by = c("Regime", "SSB", "R", "Category")

ggplot2::ggsave(g2, file = "SRregime_curve.png", unit = "mm", width = 240, height= 120,
dpi = 600)

推定値の収束の有無や最適解に達しているかをチェック
fit.SRregimeの場合でもcheck.SRfitという関数を使って、適切な推定値が得られているかを診断できます33

（fit.SRregimeのかわりにfit.SR_tol(fit.SRregimeへの引数, is_regime=TRUE)とすると、以下の診断の結果得られるパラメータが⾃動的に得られます）
以下の例の場合resR2の⽅でHSの折れ点が最⼤値に達しており、（境界付近の）異なるパラメータ値
でも同じ対数尤度が得られるという診断結果となっていますcheck_R1 = check.SRfit(resR1, output = TRUE, filename = "check_R1")

1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっています (OK)
3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)
4. パラメータが⼤域解に達しているのでOKです (OK)
5. パラメータが唯⼀の解として推定されているのでOKです (OK)check_R2 = check.SRfit(resR2, output = TRUE, filename = "check_R2")

1. 収束しています (OK)
2. Hessian⾏列の対⾓成分が正定値になっていません (HSかつ3以降のチェックがOKであれば問題ありません)
3. パラメータは壁にあたっています(HSで折れ点が過去最⼩・最⼤親⿂量になっているときにそうなります。HSでない
場合は推定の不安定性を⽰唆します。)Parameter b of regime 1 reaching the maximum SSB
4. パラメータが⼤域解に達しているのでOKです (OK)
5. 同じ最⼤尤度(1e-06よりも⼩さい違い)を持つ複数のパラメータが⾒つかりました（L1かつHSでよく⾒られます）。
同じ尤度を持つパラメータの範囲 ( 100 回試⾏),
regime2
a1
a2
b1
b2
sd1
sd2
Min.
1 0.03162750 0.03162750 642690.9 3592760 0.6888562 0.5960809
1st Qu.
1 0.03162946 0.03162946 642866.8 3592767 0.6888669 0.5961052
Median1 0.03163016 0.03163016 642956.5 3592769 0.6888708 0.5961176
Mean
1 0.03163016 0.03163016 642944.4 3592768 0.6888708 0.5961175
3rd Qu.
1 0.03163083 0.03163083 643023.3 3592769 0.6888745 0.5961307
Max.
1 0.03163213 0.03163213 643148.5 3592769 0.6888816 0.5961672
中央値に最も近いパラメータセットを持つ推定結果をoptimumに出⼒します(そのときの初期値は0.0311954100844516-0.841070178834585-13.5729167838772 です)

残差の分布のチェック
残差の確率分布のプロットはfit.SRと同様に出⼒できます34check.SRdist(resR1)

残差の⾃⼰相関のチェック
レジームを推定してから⾃⼰相関を事後的に推定することができます
calc.residARを使い、per_regime = TRUEでレジームごとに⾃⼰相関を推定、per_regime =FALSEでレジーム間で共通の⾃⼰相関係数を推定できます
ここで⾃⼰相関係数が⼤きかったりAICcが⼤きく下がったからと⾔って将来予測および管理基準値計
算に⾃⼰相関を考慮すべきかはわかりません。個⼈的には、⾃⼰相関はレジームの代替⼿段として考慮するものですし、⾃⼰相関とレジームを両⽅考慮した再⽣産関係を⾒たことがないので、⾃⼰相関は考慮しないでいいのではないかと思います

outer_R1 = calc.residAR(resR1, per_regime = TRUE)

Joining, by = "regime"

outer_R1$regime_pars# A tibble: 2 × 5
regimeab
sd
rho<dbl> <dbl><dbl> <dbl> <dbl>
1
0 0.0275 740401. 0.644 0.313
2
1 0.0633 1460893. 0.346 0.389

outer_R1$AICc
35

AR(0)
AR(1)107.0302 104.7579

outer_R0 = calc.residAR(resR1, per_regime = FALSE)

Joining, by = "regime"

outer_R0$regime_pars# A tibble: 2 × 5
regimeab
sd
rho<dbl> <dbl><dbl> <dbl> <dbl>
1
0 0.0275 740401. 0.644 0.328
2
1 0.0633 1460893. 0.347 0.328

outer_R0$AICc

AR(0)
AR(1)
107.0302 102.6727

もしどうしてもレジームありで⾃⼰相関を内部推定したい場合には、fit.SRの引数にいれるSRdataのwightを変更することで、レジームを完全に分けての推定はできます

regime.year = c(1976)
regime.year2 = c(1988)# regime0を取り出したい場合
SRdata_regime0 <- SRdata
SRdata_regime0$weight[which(SRdata$year < regime.year)] <- 0
SRdata_regime0$weight[which(SRdata$year >= regime.year2)] <- 0
SRdata_regime0$weight <- as.numeric(!SRdata_regime0$weight)
resR3 <- fit.SR(SRdata_regime0, SR = "HS", method = "L2", AR = 1, out.AR = FALSE,
bio_par = bio_par)resR3$pars

36

a
b
sd
rho
1 0.02756965 741716.6 0.6448769 0.2997622# regime1を取り出したい場合
SRdata_regime1 <- SRdata
SRdata_regime1$weight[which(SRdata$year < regime.year)] <- 0
SRdata_regime1$weight[which(SRdata$year >= regime.year2)] <- 0
resR3_1 <- fit.SR(SRdata_regime1, SR = "HS", method = "L2", AR = 1, out.AR =FALSE,
bio_par = bio_par)resR3_1$pars

a
b
sd
rho
1 0.07682997 943941.6 0.3315344 0.6648866

⾃⼰相関に関するプロットはautocor.plotで出⼒されますが、標準偏差が年代によって変わっているので標準化残差を使っています

autocor.plot(resR1)

Standardized residuals are used for 'fit.SRregime'

スティープネスfit.SRregime関数の引数にbio_parオブジェクトを代⼊してあればレジームごとに結果が返されます
レジームによって⽣物パラメータが異なるような状況にはまだ対応していません
37

resR2$steepness

SPR0
SB0
R0
B0
h regime
1 166.5378 3386593 20335.28 4262924 0.8101652
0
2 166.5378 18925763 113642.43 23823078 0.8101652
1

ブートストラップ
fit.SRregimeでもboot.SRを使ってfit.SRと同様に、①残差のパラメトリックブートストラップ、②残差のノンパラメトリックブートストラップ、③データのブートストラップを⾏うことができます
残差の標準偏差がレジームによって変わるので、標準化残差に直したうえで残差ブートストラップをやり、元の標準偏差の⼤きさに戻しています
ここでは例として残差のノンパラメトリックブートストラップを⽰します。レジーム毎に各パラメー
タのヒストグラムが表⽰されます

boot.res4 = boot.SR(resR1, n = 100, method = "n")bootSR.plot(boot.res4, ggplt = T)

Joining, by = "name"`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Joining, by = "name"38`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

39

40

レジーム間で共通のパラメータは各レジームで同じ図となります

boot.res5 = boot.SR(resR2, n = 100, method = "n")bootSR.plot(boot.res5, ggplt = T)

Joining, by = "name"`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Joining, by = "name"41`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

42

43

ジャックナイフジャックナイフも同様にレジームごとに結果が表⽰されます
レジームシフトの位置を⻘点線で⽰しています
レジーム間で共通のパラメータがある場合はレジーム間で同じ図が表⽰されますjack2 = jackknife.SR(resR1, is.plot = TRUE)

44

45

jack3 = jackknife.SR(resR2, is.plot = TRUE)

46

47

プロファイル尤度
プロファイル尤度もfit.SRのときと同様の関数prof.likSRで計算・プロットできます
レジームごとにa,bのパラメータを動かし、別のレジームで推定が必要なパラメータは尤度が最⼤になるように推定しています（そのため計算に時間がかかります）
図はレジームごとに出⼒されますprof2 = prof.likSR(resR1, length = 10)
ゃ時間かかります# ほんとはlength=50くらいのほうが良い？でもめちゃくちパラメータ間相関48

レジームがある場合でもパラメータ間の相関係数が計算できます
パラメータはa,bのセットが２つ⽣成されますので、それぞれのセットの相関係数が出⼒されます

res_corSR_R1 <- corSR(resSR = resR1)
res_corSR_R1$cor[1, 3][1] -0.3859249

res_corSR_R1$cor[2, 4][1] -0.7071037

49