管理基準値・ABC計算チュートリアル
1. 事前準備
- データの読み込み，RVPA関数の読み込みなど
# 関数の読み込み →
# warningまたは「警告」が出るかもしれませんが，その後動いていれば問題ありません
source("../../rvpa1.9.2.r")
source("../../future2.1.r")
source("../../utilities.r", encoding = "UTF-8") # ggplotを使ったグラフ作成用の関数
# データの読み込み
caa <- read.csv("caa_pma.csv", row.names = 1)
waa <- read.csv("waa_pma.csv", row.names = 1)
maa <- read.csv("maa_pma.csv", row.names = 1)
dat <- data.handler(caa = caa, waa = waa, maa = maa, M = 0.5)
names(dat)[1] "caa" "maa" "waa" "index" "M"
[6] "maa.tune" "waa.catch" "catch.prop"2. VPAによる資源量推定
今後はvpa関数の返り値，res.pmaを使って将来予測計算をおこなっていくので，そのためにvpaを実施します．(この辺はあまり詳しく解説しません．)
# VPAによる資源量推定
res.pma <- vpa(dat, fc.year = 2009:2011, rec = 585, rec.year = 2011, tf.year = 2008:2010,
term.F = "max", stat.tf = "mean", Pope = TRUE, tune = FALSE, p.init = 1)res.pma$Fc.at.age # 将来予測やMSY計算で使うcurrent F (fc.yearのオプションでいつのFの平均かが指定される) 0 1 2 3
0.5436309 1.1766833 1.3059521 1.3059519 plot(res.pma$Fc.at.age, type = "b", xlab = "Age", ylab = "F", ylim = c(0, max(res.pma$Fc.at.age)))
2.5 VPA結果を外部から読み込む場合
- read.vpa関数を使って読み込めます
- out.vpa関数による出力と同じ書式になりますので、out.vpa関数でひな形ファイルを作成してから、エクセルでそれを編集し、read.vpa関数で読むと良いと思います
out.vpa(res.pma) # vpa.csvというファイルが作成されます。VPAの結果のグラフ出力となるvpa.pdfも出力されます。
res.pma2 <- read.vpa("vpa.csv") # vpa.csvを編集後、read.vpa関数で読み込みますPope is TRUE... OK?3. 再生産関係を仮定しない管理基準値の計算
- ref.F関数を使います
- %SPRやFmaxなど、再生産関係を仮定しない管理基準値を計算します
- 計算結果はrres.pmaに格納されます
- YPR, SPR曲線とFcurrent (
res.pma$Fc.at.aに入っている値です), Fmax, Fmed, F0.1などの位置が表示されます
byear <- 2009:2011 # 生物パラメータを平均する期間を2009年から2011年とする
rres.pma <- ref.F(res.pma, # VPAの計算結果
waa.year=byear, maa.year=byear, M.year=byear, # weight at age, maturity at age, Mは2009から2011年までの平均とする
rps.year=2000:2011, # Fmedを計算するときに用いるRPSの範囲
max.age=Inf, # SPR計算で仮定する年齢の最大値
pSPR=c(10,20,30,35,40), # F_%SPRを計算するときに，何パーセントのSPRを計算するか
Fspr.init=1)
図：YPR, SPR曲線
- 結果のサマリーは
rres.pma$summaryによって見れます - max: F at ageの最大値，mean: F at ageの平均値，Fref/Fcur: Fcurrentを分母にしたときのF管理基準値の比
- この結果から，現状のF（Fcurrent）はFmedとほぼ同等（Fref/Fcur=0.96なので），F％SRP=10％くらいであることがわかります
rres.pma$summary Fcurrent Fmed Flow Fhigh Fmax F0.1
max 1.305952 1.2545205 1.482739 1.0427937 0.7069380 0.3956762
mean 1.083055 1.0404012 1.229668 0.8648116 0.5862791 0.3281429
Fref/Fcur 1.000000 0.9606176 1.135370 0.7984931 0.5413200 0.3029791
Fmean FpSPR.10.SPR FpSPR.20.SPR FpSPR.30.SPR FpSPR.35.SPR
max 1.2552763 1.434605 0.8363638 0.5648693 0.4739020
mean 1.0410280 1.189749 0.6936147 0.4684585 0.3930172
Fref/Fcur 0.9611963 1.098512 0.6404245 0.4325345 0.3628785
FpSPR.40.SPR
max 0.4000316
mean 0.3317549
Fref/Fcur 0.3063142# 横軸や縦線で示す管理基準値を調整する場合、plot.Fref関数を使う x.labelは
# res.pma$summaryの行名、vline.textは
# res.pma$summaryの列名前に対応させて指定する
plot.Fref(rres.pma, xlabel = "Fref/Fcur", vline.text = c("FpSPR.20.SPR", "FpSPR.30.SPR",
"FpSPR.40.SPR"))
図：YPR, SPR曲線 (x軸などを変更した場合)
4. 再生産関係の推定
データの作成
- get.SRdataを使って再生産関係のフィット用のデータを作る
- get.SRdata関数では，
rownames(res.pma$naa)を参照し、必要な年齢分のSSBをずらしたデータを作成する - yearは加入年
# VPA結果を使って再生産データを作る
SRdata <- get.SRdata(res.pma)
head(SRdata)$year
[1] 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995
[15] 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009
[29] 2010 2011
$SSB
[1] 12199.02 15266.68 15072.03 19114.22 23544.42 28769.36 34764.44
[8] 38219.49 48535.10 61891.08 63966.56 38839.78 53404.29 47322.39
[15] 54485.00 54385.04 47917.14 46090.77 59847.97 53370.62 48781.21
[22] 42719.43 40095.38 39311.72 38335.98 37564.69 27604.88 23079.33
[29] 22902.91 21427.90
$R
[1] 406.0086 498.9652 544.3007 469.6025 1106.8877 1043.4237 696.7049
[8] 923.9567 1353.1790 1698.8457 1117.5454 2381.1352 1669.1381 1818.3638
[15] 1858.0043 1458.9524 1334.9288 1116.9433 1100.4598 1693.9770 1090.7181
[22] 1081.8377 1265.0559 1024.0418 753.5207 764.9720 879.6573 513.1128
[29] 570.9195 585.0000# SSBとRのデータだけを持っている場合
SRdata0 <- get.SRdata(R.dat = exp(rnorm(10)), SSB.dat = exp(rnorm(10)))
# 特定の期間のデータだけを使う場合
SRdata0 <- get.SRdata(res.pma, years = 1990:2000)モデルのフィット
- HS,BH,RIをフィットし，再生産関係のパラメータを推定する
- 結果のオブジェクトのAICcにAICcの値が入っているので，それを比較し，再生産関係を決定する
- SR.fitオプション
- SR:再生産関係のタイプ： “HS”（ホッケー・スティック）、“BH”（べバートン・ホルト）、“RI”（リッカー）
- AR: 自己相関の考慮なし(AR=1)、過去１年分の自己相関を考慮(AR=1) （１年分しか対応していない）
- method: 最小二乗法（“L2”)か最小絶対値法（“L1”）
- 自己相関あり・なしでAICcを比較し、自己相関を入れたほうがいいかどうか判断する
- \(\log(R_t)=\log(HS(SSB_t))+\rho \times {\log(R_{t-1})-\log(HS(SSB_{t-1}))}\)
- \(\log(R_t)~N(\log(R_t),\sigma^2)\)
- 自己相関パラメータrhoの推定については不安定な部分があります。計算方法の改善により今後値が変わる可能性があります
- この例の場合はHSでARなしで最もAICcが小さい→MSY計算ではHS.par0の結果を使う
HS.par0 <- fit.SR(SRdata, SR = "HS", method = "L2", AR = 0, hessian = FALSE)
HS.par1 <- fit.SR(SRdata, SR = "HS", method = "L2", AR = 1, hessian = FALSE)
BH.par0 <- fit.SR(SRdata, SR = "BH", method = "L2", AR = 0, hessian = FALSE)
BH.par1 <- fit.SR(SRdata, SR = "BH", method = "L2", AR = 1, hessian = FALSE)
RI.par0 <- fit.SR(SRdata, SR = "RI", method = "L2", AR = 0, hessian = FALSE)
RI.par1 <- fit.SR(SRdata, SR = "RI", method = "L2", AR = 1, hessian = FALSE)
c(HS.par0$AICc, HS.par1$AICc, BH.par0$AICc, BH.par1$AICc, RI.par0$AICc, RI.par1$AICc)[1] 10.66778 13.27672 11.44803 14.11029 11.40670 14.05799- 結果の図示
plot.SRdata(SRdata)
lines(HS.par0$pred, col = 2, type = "l", lwd = 3)
lines(BH.par0$pred, col = 3, type = "l", lwd = 3)
lines(RI.par0$pred, col = 4, type = "l", lwd = 3)
図：観測値（○）に対する再生産関係式．plot=赤がHS，緑と青がBH, RIだが両者はほとんど重なっていて見えない
- TMBオプション(
TMB=TRUE)も使えます（ちょっと不安定です。使いたい場合はお問い合わせください）
autoregressiveSR2.cppをダウンロードして，作業フォルダに置く
# install.packages('TMB') #TMBがインストールされてなければ
library(TMB)
compile("autoregressiveSR2.cpp")
dyn.load(dynlib("autoregressiveSR2"))
HS.par11 <- fit.SR(SRdata, SR = "HS", method = "L2", AR = 1, TMB = TRUE) #marginal likelihoodモデル診断
再生産関係のあてはめのあとは、推定されたパラメータの信頼区間や頑健性などをチェックする必要があります。そのための関数群なども用意しています。詳しくは SRRガイドライン へ
5. 将来予測
future.vpa関数を使います
- recfuncの引数に再生産関係の関数を，rec.argにrecfuncに対する引数（再生産関係のパラメータ）を入れる
- 利用可能な再生産関数
- HS.recAR: ホッケー・スティック＋加入は対数正規分布＋自己相関ありの場合も対応
- RI.recAR・BH.recAR：HS.recARのリッカー・べバートンホルトバージョン
- HS.rec, BH.rec, RI.rec : 残差リサンプリング用
fres.HS <- future.vpa(res.pma,
multi=1, # res.pma$Fc.at.ageに掛ける乗数
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
is.plot=TRUE, # 結果をプロットするかどうか
seed=1,
silent=TRUE,
recfunc=HS.recAR, # 再生産関係の関数
# recfuncに対する引数
rec.arg=list(a=HS.par0$pars$a,b=HS.par0$pars$b,
rho=HS.par0$pars$rho, # ここではrho=0なので指定しなくてもOK
sd=HS.par0$pars$sd,resid=HS.par0$resid))
図：is.plot=TRUEで表示される図．資源量(Biomass)，親魚資源量(SSB), 漁獲量(Catch)の時系列．決定論的将来予測（Deterministic），平均値（Mean），中央値(Median)，80％信頼区間を表示
将来予測で自己相関を考慮する場合
fres.HS.AR <- future.vpa(res.pma,
multi=1,
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,is.plot=FALSE,
seed=1, silent=TRUE,recfunc=HS.recAR,
# recfuncに対する引数 => 自己相関ありのオプションで計算した結果を入れる
rec.arg=list(a=HS.par1$pars$a,b=HS.par1$pars$b,
rho=HS.par1$pars$rho, # 自己相関が高い場合、この値が>0となる
sd=HS.par1$pars$sd,
resid=HS.par1$resid, # 再生産関係における残差の時系列
resid.year=NULL # 近年の残差を何年分平均して将来予測に使うか？NULLの場合は、最後の年の残差を使う
)
)Beverton-Holtを仮定する場合
fres.BH <- future.vpa(res.pma,
multi=1,
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
is.plot=TRUE, # 結果をプロットするかどうか
seed=1,
silent=TRUE,
recfunc=BH.recAR, # 再生産関係の関数
# recfuncに対する引数
rec.arg=list(a=BH.par0$pars$a,b=BH.par0$pars$b,rho=BH.par0$rho,
sd=BH.par0$pars$sd,resid=BH.par0$resid))
図：is.plot=TRUEで表示される図．資源量(Biomass)，親魚資源量(SSB), 漁獲量(Catch)の時系列．決定論的将来予測（Deterministic），平均値（Mean），中央値(Median)，80％信頼区間を表示
同じ引数を使ってもう一度将来予測をする
fres.HS$inputに、将来予測で使った引数が入っているので、それにdo.call(関数、引数)すると同じ計算を繰り返せる
fres.HS2 <- do.call(future.vpa, fres.HS$input)
- fres.HS$inputを上書きすることで，同じ引数を使いながら設定を少しだけ変更した将来予測が実行できる
- 引数
multiがcurrent Fへの乗数になる - たとえばmulti=1からmulti=0.5に変更する例は以下のとおり
# 引数をinput.tmpに代入．
input.tmp <- fres.HS2$input
# 引数の一部を変える
input.tmp$multi <- 0.5 # current Fの1/2で漁獲
fres.HS3 <- do.call(future.vpa, input.tmp)
plot.futures関数を使って複数の結果を比較
par(mfrow = c(2, 2))
plot.futures(list(fres.HS, fres.HS3), legend.text = c("F=Fcurrent", "F=0.5Fcurrent"),
target = "SSB")
plot.futures(list(fres.HS, fres.HS3), legend.text = c("F=Fcurrent", "F=0.5Fcurrent"),
target = "Catch")
plot.futures(list(fres.HS, fres.HS3), legend.text = c("F=Fcurrent", "F=0.5Fcurrent"),
target = "Biomass")
図：plot.futures関数の結果
(5-1) Fの設定やFrec
将来予測における漁獲のシナリオ
- future.vpaの引数
ABC.yearで指定した年から，Fcurrent × multiによるFで漁獲される - ABC.year-1年まではFcurrentによる漁獲
- Frecに引数を与えることで，任意の資源量に任意の確率で回復させるような将来予測ができます．
Frecのオプション
| オプション | 説明 |
|---|---|
| stochastic | 確率的将来予測をもとにFrecを計算するかどうか |
| future.year | 条件を満たしているかどうかを判断する年 |
| Blimit | 条件として使われる閾値 |
| scenario | =“blimit”: Blimitを下回る確率をtarget.probsにする |
| =“catch.mean”: future.year年の平均漁獲量をBlimitの値と一致させる | |
| =“ssb.mean”: future.year年の平均親魚量をBlimitの値と一致させる | |
| target.probs | scenario=“blimit”のときに目的とする確率（パーセントで指定） |
| Frange | 探索するFの範囲．指定しない場合，c(0.01,multi*2)の範囲で探索しますので，うまく推定できない場合はfuture.vpaの引数multiを変えるか，このオプションでそれらしいFの値に限定してください| |
# たとえば現状の資源量に維持するシナリオ
fres.currentSSB <- future.vpa(res.pma,
multi=0.8,
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,seed=1,
is.plot=TRUE, # 結果をプロットするかどうか
Frec=list(stochastic=TRUE,future.year=2023,Blimit=rev(colSums(res.pma$ssb))[1],scenario="blimit",target.probs=50),
recfunc=HS.recAR, # 再生産関係の関数
# recfuncに対する引数
rec.arg=list(a=HS.par0$pars$a,b=HS.par0$pars$b,
rho=HS.par0$pars$rho,
sd=HS.par0$pars$sd,bias.corrected=TRUE))$ABC.year
[1] 2013
$Blim
[1] 0
$F.sigma
[1] 0
$Frec
$Frec$stochastic
[1] TRUE
$Frec$future.year
[1] 2023
$Frec$Blimit
2011
21427.9
$Frec$scenario
[1] "blimit"
$Frec$target.probs
[1] 50
$HCR
NULL
$M
NULL
$M.year
[1] 2009 2010 2011
$N
[1] 100
$Pope
[1] TRUE
$add.year
[1] 0
$currentF
NULL
$delta
NULL
$det.run
[1] TRUE
$eaa0
NULL
$faa0
NULL
$is.plot
[1] TRUE
$maa
NULL
$maa.year
[1] 2009 2010 2011
$multi
[1] 0.8
$multi.year
[1] 1
$naa0
NULL
$nyear
[1] 50
$outtype
[1] "FULL"
$plus.group
[1] TRUE
$pre.catch
NULL
$random.select
NULL
$rec.arg
$rec.arg$a
[1] 0.02861539
$rec.arg$b
[1] 51935.41
$rec.arg$rho
[1] 0
$rec.arg$sd
[1] 0.2575537
$rec.arg$bias.corrected
[1] TRUE
$rec.new
NULL
$recfunc
function (ssb, vpares, rec.resample = NULL, rec.arg = list(a = 1000,
b = 1000, sd = 0.1, rho = 0, resid = 0))
{
rec0 <- ifelse(ssb > rec.arg$b, rec.arg$a * rec.arg$b, rec.arg$a *
ssb)
rec <- rec0 * exp(rec.arg$rho * rec.arg$resid)
rec <- rec * exp(rnorm(length(ssb), -0.5 * rec.arg$sd2^2,
rec.arg$sd))
new.resid <- log(rec/rec0) + 0.5 * rec.arg$sd2^2
return(list(rec = rec, rec.resample = new.resid))
}
<bytecode: 0x26627c8>
$replace.rec.year
[1] 2012
$seed
[1] 1
$silent
[1] FALSE
$start.year
[1] 2012
$use.MSE
[1] FALSE
$waa
NULL
$waa.catch
NULL
$waa.fun
[1] FALSE
$waa.year
[1] 2009 2010 2011
F multiplier= 1.03432 
Frecオプションを使った場合は、結果の図に目的とする年・資源量のところに赤線が入ります。これが将来予測の結果と一致しているか確かめてください。もし一致していない場合、multi（初期値）かFrecのオプションのFrangeを指定してやり直してください
(5-2) 再生産関係
- 残差リサンプリングで将来予測をする場合→refuncとしてHS.rec（ホッケー・スティック）、BH.rec（べバートン・ホルト）、RI.rec（リッカー）を使う
- rec.argの引数で、必ず resample=TRUEとしてください。
- rho>0の場合には対応しておりません
# 残差リサンプリングによる将来予測
fres.HS4 <- future.vpa(res.pma,
multi=1,
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
is.plot=TRUE, # 結果をプロットするかどうか
seed=1,
recfunc=HS.rec, # 再生産関係の関数（HS.rec=Hockey-stick)
rec.arg=list(a=HS.par0$pars$a,b=HS.par0$pars$b,
rho=HS.par0$pars$rho,
sd=HS.par0$pars$sd,bias.correction=TRUE,
resample=TRUE,resid=HS.par0$resid))$ABC.year
[1] 2013
$Blim
[1] 0
$F.sigma
[1] 0
$Frec
NULL
$HCR
NULL
$M
NULL
$M.year
[1] 2009 2010 2011
$N
[1] 100
$Pope
[1] TRUE
$add.year
[1] 0
$currentF
NULL
$delta
NULL
$det.run
[1] TRUE
$eaa0
NULL
$faa0
NULL
$is.plot
[1] TRUE
$maa
NULL
$maa.year
[1] 2009 2010 2011
$multi
[1] 1
$multi.year
[1] 1
$naa0
NULL
$nyear
[1] 50
$outtype
[1] "FULL"
$plus.group
[1] TRUE
$pre.catch
NULL
$random.select
NULL
$rec.arg
$rec.arg$a
[1] 0.02861539
$rec.arg$b
[1] 51935.41
$rec.arg$rho
[1] 0
$rec.arg$sd
[1] 0.2575537
$rec.arg$bias.correction
[1] TRUE
$rec.arg$resample
[1] TRUE
$rec.arg$resid
[1] 0.15107386 0.13291912 0.23271699 -0.15249042 0.49647422
[6] 0.23700706 -0.35617750 -0.16862488 -0.02601987 0.13375888
[11] -0.28505548 0.76194318 0.11611724 0.29476430 0.22331281
[16] -0.01847154 -0.02678452 -0.16620607 -0.30046205 0.13088884
[21] -0.24669812 -0.12218191 0.09766607 -0.09395451 -0.37557687
[26] -0.34016983 0.10759521 -0.25238643 -0.13796039 -0.04702647
$rec.new
NULL
$recfunc
function (ssb, vpares, rec.resample = NULL, rec.arg = list(a = 1000,
b = 1000, sd = 0.1, resample = FALSE, resid = 0, bias.correction = TRUE))
{
rec0 <- ifelse(ssb > rec.arg$b, rec.arg$a * rec.arg$b, rec.arg$a *
ssb)
if (!isTRUE(rec.arg$resample)) {
if (isTRUE(rec.arg$bias.correction)) {
rec <- rec0 * exp(rnorm(length(ssb), -0.5 * (rec.arg$sd)^2,
rec.arg$sd))
}
else {
rec <- rec0 * exp(rnorm(length(ssb), 0, rec.arg$sd))
}
}
else {
if (isTRUE(rec.arg$bias.correction)) {
rec <- c(rec0[1], exp(log(rec0[-1]) + sample(rec.arg$resid,
length(ssb) - 1, replace = TRUE))/mean(exp(rec.arg$resid)))
}
else {
rec <- c(rec0[1], exp(log(rec0[-1]) + sample(rec.arg$resid,
length(ssb) - 1, replace = TRUE)))
}
}
return(list(rec = rec, rec.resample = rec.arg$resid))
}
$replace.rec.year
[1] 2012
$seed
[1] 1
$silent
[1] FALSE
$start.year
[1] 2012
$use.MSE
[1] FALSE
$waa
NULL
$waa.catch
NULL
$waa.fun
[1] FALSE
$waa.year
[1] 2009 2010 2011
残差リサンプリングか対数正規分布かの違いを比較
par(mfrow = c(2, 2))
plot(fres.HS$vssb[, -1], fres.HS$naa[1, , -1], xlab = "SSB", ylab = "Recruits")
plot(fres.HS4$vssb[, -1], fres.HS4$naa[1, , -1], xlab = "SSB", ylab = "Recruits")
plot.futures(list(fres.HS, fres.HS4)) # 両者の比較(5-3) 年齢別体重が資源尾数に影響される場合の将来予測
- future.vpaで，waa.fun = TRUEとすれば、年齢別資源重量が資源尾数（log(体重)~log(資源尾数)の回帰を関数内部で実行）の関数から予測されます
- 不確実性も考慮されます
- 30系群であてはめた例はこちら (データは1年分古いです)
- 太平洋マイワシ，対馬マイワシ，太平洋マサバ，ホッケ，瀬戸内サワラでは年齢別体重と年齢別資源尾数に関係がありそうなかんじです
lm.res <- plot.waa(res.pma) # weight at ageが資源尾数の関数になっているかどうか，確認してみる．この例の場合は特に有意な関係はない
# lm.resの中に回帰した結果が年齢分だけ入っています
fres.HS6 <- fres.HS
fres.HS6$input$waa.fun <- TRUE
fres.HS6$input$N <- 1000
fres.HS6 <- do.call(future.vpa, fres.HS6$input)
(5-4) その他のオプション
HCR=list(Blim=154500, Bban=49400,beta=1)のように与えると、BlimからBbanまで直線的にFが減少するようなHCRで漁獲します。betaは、全体のFに掛ける乗数です。pre.catch=list(year=2013,wcatch=13000)のように与えると特定の年の漁獲量をwcatchで与えた漁獲量で固定しますrec.new=list(year=2012,rec=1500)のように与えると、対応する年の加入を置き換えます。
# 残差リサンプリングによる将来予測
fres.HS5 <- future.vpa(res.pma,
multi=1,
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
recfunc=HS.rec,
rec.arg=list(a=HS.par0$pars$a,b=HS.par0$pars$b,
rho=HS.par0$pars$rho,
sd=HS.par0$pars$sd,bias.correction=TRUE,
resample=TRUE,resid=HS.par0$resid),
rec.new=list(year=2012,rec=100),
pre.catch=list(year=2013,wcatch=100)
)$ABC.year
[1] 2013
$Blim
[1] 0
$F.sigma
[1] 0
$Frec
NULL
$HCR
NULL
$M
NULL
$M.year
[1] 2009 2010 2011
$N
[1] 100
$Pope
[1] TRUE
$add.year
[1] 0
$currentF
NULL
$delta
NULL
$det.run
[1] TRUE
$eaa0
NULL
$faa0
NULL
$is.plot
[1] TRUE
$maa
NULL
$maa.year
[1] 2009 2010 2011
$multi
[1] 1
$multi.year
[1] 1
$naa0
NULL
$nyear
[1] 50
$outtype
[1] "FULL"
$plus.group
[1] TRUE
$pre.catch
$pre.catch$year
[1] 2013
$pre.catch$wcatch
[1] 100
$random.select
NULL
$rec.arg
$rec.arg$a
[1] 0.02861539
$rec.arg$b
[1] 51935.41
$rec.arg$rho
[1] 0
$rec.arg$sd
[1] 0.2575537
$rec.arg$bias.correction
[1] TRUE
$rec.arg$resample
[1] TRUE
$rec.arg$resid
[1] 0.15107386 0.13291912 0.23271699 -0.15249042 0.49647422
[6] 0.23700706 -0.35617750 -0.16862488 -0.02601987 0.13375888
[11] -0.28505548 0.76194318 0.11611724 0.29476430 0.22331281
[16] -0.01847154 -0.02678452 -0.16620607 -0.30046205 0.13088884
[21] -0.24669812 -0.12218191 0.09766607 -0.09395451 -0.37557687
[26] -0.34016983 0.10759521 -0.25238643 -0.13796039 -0.04702647
$rec.new
$rec.new$year
[1] 2012
$rec.new$rec
[1] 100
$recfunc
function (ssb, vpares, rec.resample = NULL, rec.arg = list(a = 1000,
b = 1000, sd = 0.1, resample = FALSE, resid = 0, bias.correction = TRUE))
{
rec0 <- ifelse(ssb > rec.arg$b, rec.arg$a * rec.arg$b, rec.arg$a *
ssb)
if (!isTRUE(rec.arg$resample)) {
if (isTRUE(rec.arg$bias.correction)) {
rec <- rec0 * exp(rnorm(length(ssb), -0.5 * (rec.arg$sd)^2,
rec.arg$sd))
}
else {
rec <- rec0 * exp(rnorm(length(ssb), 0, rec.arg$sd))
}
}
else {
if (isTRUE(rec.arg$bias.correction)) {
rec <- c(rec0[1], exp(log(rec0[-1]) + sample(rec.arg$resid,
length(ssb) - 1, replace = TRUE))/mean(exp(rec.arg$resid)))
}
else {
rec <- c(rec0[1], exp(log(rec0[-1]) + sample(rec.arg$resid,
length(ssb) - 1, replace = TRUE)))
}
}
return(list(rec = rec, rec.resample = rec.arg$resid))
}
<bytecode: 0x9be7378>
$replace.rec.year
[1] 2012
$seed
NULL
$silent
[1] FALSE
$start.year
[1] 2012
$use.MSE
[1] FALSE
$waa
NULL
$waa.catch
NULL
$waa.fun
[1] FALSE
$waa.year
[1] 2009 2010 2011
6. MSY管理基準値の計算
- MSY管理基準値計算では，上記の将来予測において，Fの値を様々に変えたときの平衡状態（世代時間×20年を
nyearで指定します）における資源量やそれに対応するF等を管理基準値として算出します - なので、ここまでのプロセスで、ABC計算のためにきちんとしたオプションを設定したfuture.vpaを実行しておいてください。その返り値
fres.HSをMSY計算では使っていきます
est.MSYの説明
- この関数で計算できる管理基準値は以下のようなものになります
- どの管理基準値がtarget, limit, banになるかは関数内では評価されません
| est.MSYの引数（重要なもののみ） | 説明 |
|---|---|
| vpares | VPA結果のオブジェクト |
| farg | 将来予測関数future.vpaへの引数のリスト。通常はfuture.vpaを実施したあとのオブジェクトがfoutだった場合、fout$inputと指定すれば良いです |
| trace.multi | Fmsyを探索したり、Yield curveを書くときにグリッドサーチをするときのFの刻み。Fcurrentに対する乗数。Fが異常に大きい場合、親魚=0になって加入＝NAとなるので現実的な値を入れる |
| long.term | 世代時間の何倍年後の状態を平衡状態と仮定するか |
| mY | 自己相関を考慮して管理基準値を計算する場合、平衡状態から何年分進めたときの資源量を管理基準値として使うか |
| resid.year | ARありの場合、最近何年分の残差を平均するかをここで指定する。ARありの設定を反映させたい場合必ずここを１以上とすること（とりあえず１としておいてください）。 |
| current.resid | 残差の値を直接入れる場合。上の年数が設定されていてもこちらが設定されたらこの値を使う |
| PGY | MSYの何割の平均漁獲量を得るときの親魚資源量を管理基準値として推定するか。比率で入力する |
| B0percent | F=0のときの平均親魚資源量の一定割合に相当する親魚量を管理基準値として推定するか。比率で入力する |
| Bempirical | 現行の親魚資源量、HSの折れ点など、特定の親魚資源量で平衡状態になるとしたときの管理基準値。親魚資源量の絶対値で与える。 |
| 管理基準値 | 説明 |
|---|---|
| SSB_MSY | 平衡状態において平均最大漁獲量が最大になるときの親魚量 |
| SSB_0 (XX%) | F=0で将来予測したときの平衡状態における親魚量(\(B_0\))に対する割合（引数B0percentでc(0.4, 0.5)のように指定します） |
| SSB_PGY (LXX%) (HXX%) | SS_MSYで達成される漁獲量のXX%を達成するときの親魚量の下限または上限（引数PGYでc(0.9, 0.95)のように指定します） |
| SSB_empirical | 過去最大・最小漁獲量、など、経験的な親魚資源量の値 |
| 関数の返り値 | 説明 |
|---|---|
| summay | 平衡状態における代表的な各種統計量（SSB・総資源量・漁獲量等の平均値やFの値）*1 |
| summayAR | 直近の加入の残差を考慮した場合に、平衡状態のmY年後における各種統計量（SSB・総資源量・漁獲量等の平均値やFの値）*1 |
| summay_tb | 上記の2つのsummaryを縦に結合したもの。tibble形式になっている。 |
| all.stat | 平衡状態における各種統計量（summaryよりも詳しい) |
| all.statAR | 直近の加入の残差を考慮した場合の各種統計量 |
| all.stat_tb | 上記の2つのall.statを縦に結合したもの。tibble形式になっている。 |
| trace | Fcurrentに対するmultiplierを様々に変えた場合の平衡状態（GT*20年後）における各種統計量 | | input.list | 各種管理基準値を計算するときに使用したfuture.vpaへの引数。do.call(future.vpa,引数)で計算の再現が可能 | |
*1: summaryまたはsummaryARのFref/Fcurrentが現行のFからのFの削減率になります（（Fref/Fcurrent-1)×100が資源評価票の要約表の「現状のF値からの増減％」に相当します）。この値にさらにβ（Btargetを上回る確率が５０％かつBlimitを上回る確率が９０％以上になるように調整する係数）と(B-Bban)/(Blim-Bban)を乗じたFをもとにABCが算定されます
# MSY管理基準値の計算
MSY.HS <- est.MSY(res.pma, # VPAの計算結果
fres.HS$input, # 将来予測で使用した引数
# nyear=NULL, # 何年計算するかは、指定しなければ関数内部で世代時間の20倍の年数を計算し、それを平衡状態とする
resid.year=0, # ARありの場合、最近何年分の残差を平均するかをここで指定する。ARありの設定を反映させたい場合必ずここを１以上とすること（とりあえず１としておいてください）。
N=100, # 将来予測の年数，繰り返し回数
PGY=c(0.9,0.6,0.1), # 計算したいPGYレベル。上限と下限の両方が計算される
onlylower.pgy=FALSE, # TRUEにするとPGYレベルの上限は計算しない（計算時間の節約になる）
B0percent=c(0.3,0.4)) # 計算したいB0%レベルEstimating MSY
F multiplier= 0.523244
Estimating PGY 90 %
F multiplier= 0.2633456
F multiplier= 0.9540607
Estimating PGY 60 %
F multiplier= 0.1131719
F multiplier= 1.019431
Estimating PGY 10 %
F multiplier= 0.01319282
F multiplier= 1.079894
Estimating B0 30 %
F multiplier= 0.4342835
Estimating B0 40 %
F multiplier= 0.3071907 
図：est.MSYのis.plot=TRUEで計算完了時に表示される図．Fの強さに対する平衡状態の親魚資源量（左）と漁獲量（右）．推定された管理基準値も表示．
# 結果の表示(tibbleという形式で表示され、最初の10行以外は省略されます)
(refs.all <- MSY.HS$summary_tb)# A tibble: 20 x 15
RP_name AR SSB SSB2SSB0 B U Catch Catch.CV `Fref/Fcur`
<chr> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 MSY FALSE 1.22e5 0.250 2.16e5 0.324 70056. 0.136 0.523
2 B0 FALSE 4.88e5 1 5.94e5 0 0 NaN 0
3 PGY_0.… FALSE 2.18e5 0.446 3.17e5 0.199 63051. 0.119 0.263
4 PGY_0.… FALSE 5.74e4 0.118 1.41e5 0.447 63052. 0.181 0.954
5 PGY_0.… FALSE 3.31e5 0.680 4.34e5 0.0968 42034. 0.106 0.113
6 PGY_0.… FALSE 3.52e4 0.0721 9.18e4 0.458 42054. 0.470 1.02
7 PGY_0.… FALSE 4.64e5 0.952 5.70e5 0.0123 7001. 0.0963 0.0132
8 PGY_0.… FALSE 5.53e3 0.0113 1.49e4 0.471 7008. 1.42 1.08
9 B0-30% FALSE 1.46e5 0.300 2.42e5 0.287 69534. 0.131 0.434
10 B0-40% FALSE 1.95e5 0.400 2.94e5 0.224 65777. 0.122 0.307
11 MSY TRUE 1.23e5 0.250 2.16e5 0.325 70210. 0.141 0.523
12 B0 TRUE 4.93e5 1 5.97e5 0 0 NaN 0
13 PGY_0.… TRUE 2.20e5 0.446 3.18e5 0.199 63416. 0.122 0.263
14 PGY_0.… TRUE 5.80e4 0.118 1.41e5 0.448 63031. 0.198 0.954
15 PGY_0.… TRUE 3.35e5 0.680 4.37e5 0.0969 42338. 0.109 0.113
16 PGY_0.… TRUE 3.61e4 0.0721 9.30e4 0.459 42704. 0.481 1.02
17 PGY_0.… TRUE 4.69e5 0.952 5.73e5 0.0123 7055. 0.101 0.0132
18 PGY_0.… TRUE 5.64e3 0.0113 1.54e4 0.469 7235. 1.70 1.08
19 B0-30% TRUE 1.48e5 0.300 2.43e5 0.288 69779. 0.135 0.434
20 B0-40% TRUE 1.97e5 0.400 2.95e5 0.224 66123. 0.125 0.307
# ... with 6 more variables: Fref2Fcurrent <dbl>, F0 <dbl>, F1 <dbl>,
# F2 <dbl>, F3 <dbl>, RP.definition <chr># 全データをじっくり見たい場合
View(refs.all)
# のちの使用のために、Bmsy, Blimit, Bban, Fmsyを定義しておく
refs <- list(BmsyAR=as.numeric(MSY.HS$summaryAR$SSB[1]),
BlimAR=as.numeric(MSY.HS$summaryAR$SSB[6]),
BbanAR=as.numeric(MSY.HS$summaryAR$SSB[8]),
Bmsy=as.numeric(MSY.HS$summary$SSB[1]),
Blim=as.numeric(MSY.HS$summary$SSB[6]),
Bban=as.numeric(MSY.HS$summary$SSB[8]),
Fmsy=as.numeric(MSY.HS$summary$"Fref/Fcur"[1]),
MSY=as.numeric(MSY.HS$summary$Catch[1]),
Umsy=as.numeric(MSY.HS$summary$Catch[1])/as.numeric(MSY.HS$all.stat$biom.mean[1]))
# また、各管理基準値に対して以下のようなラベルをつける
# どの管理基準値をどのように定義するか、ここで指定します
refs.all$RP.definition <- NA
refs.all$RP.definition[refs.all$RP_name=="MSY" & refs.all$AR==FALSE] <- "Btarget0" # RP_nameがMSYでARがなしのものをBtargetとする
refs.all$RP.definition[refs.all$RP_name=="B0-20%" & refs.all$AR==FALSE] <- "Btarget1" # たとえばBtargetの代替値をいちおう示す場合
refs.all$RP.definition[refs.all$RP_name=="PGY_0.95_lower" & refs.all$AR==FALSE] <- "Btarget2"
refs.all$RP.definition[refs.all$RP_name=="PGY_0.9_lower" & refs.all$AR==FALSE] <- "Blow0"
refs.all$RP.definition[refs.all$RP_name=="PGY_0.6_lower" & refs.all$AR==FALSE] <- "Blimit0"
refs.all$RP.definition[refs.all$RP_name=="PGY_0.1_lower" & refs.all$AR==FALSE] <- "Bban0"
refs.all$RP.definition[refs.all$RP_name=="Ben-19431" & refs.all$AR==FALSE] <- "Bcurrent"
refs.all$RP.definition[refs.all$RP_name=="Ben-63967" & refs.all$AR==FALSE] <- "Bmax"
refs.all$RP.definition[refs.all$RP_name=="Ben-24000" & refs.all$AR==FALSE] <- "Blimit1"
refs.all$RP.definition[refs.all$RP_name=="Ben-51882" & refs.all$AR==FALSE] <- "B_HS"
# 定義した結果を見る
refs.all %>% select(RP_name,RP.definition)# A tibble: 20 x 2
RP_name RP.definition
<chr> <chr>
1 MSY Btarget0
2 B0 <NA>
3 PGY_0.9_upper <NA>
4 PGY_0.9_lower Blow0
5 PGY_0.6_upper <NA>
6 PGY_0.6_lower Blimit0
7 PGY_0.1_upper <NA>
8 PGY_0.1_lower Bban0
9 B0-30% <NA>
10 B0-40% <NA>
11 MSY <NA>
12 B0 <NA>
13 PGY_0.9_upper <NA>
14 PGY_0.9_lower <NA>
15 PGY_0.6_upper <NA>
16 PGY_0.6_lower <NA>
17 PGY_0.1_upper <NA>
18 PGY_0.1_lower <NA>
19 B0-30% <NA>
20 B0-40% <NA> # refs.allの中からRP.definitionで指定された行だけを抜き出す
(refs.base <- refs.all %>%
filter(!is.na(RP.definition)) %>% # RP.definitionがNAでないものを抽出
arrange(desc(SSB)) %>% # SSBを大きい順に並び替え
select(RP.definition,RP_name,SSB,Catch,U,Fref2Fcurrent)) # 列を並び替え# A tibble: 4 x 6
RP.definition RP_name SSB Catch U Fref2Fcurrent
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Btarget0 MSY 122092. 70056. 0.324 0.523
2 Blow0 PGY_0.9_lower 57416. 63052. 0.447 0.954
3 Blimit0 PGY_0.6_lower 35186. 42054. 0.458 1.02
4 Bban0 PGY_0.1_lower 5534. 7008. 0.471 1.08 7. HCRをもとにした将来予測とABC計算
- 決定された管理基準値のもとで将来予測計算をおこないます
- 資源評価最終年＋２年目の「平均」漁獲量をABCとします
- そのときに用いる親魚資源量は資源評価最終年＋２年目の親魚資源量です
input.abc <- MSY.HS$input$msy # MSY計算で使った引数を使う
input.abc$N <- 1000 # 実際に計算するときは10000以上を使ってください
input.abc$HCR <- list(Blim = refs$Blim, Bban = refs$Bban, beta = 0.8) # 算定ルールのデフォルトは0.8
input.abc$nyear <- 30 # ABC計算時には長期間計算する必要はない
input.abc$ABC.year <- 2013 # ここでABC.yearを設定しなおしてください
input.abc$is.plot <- TRUE
fres.abc1 <- do.call(future.vpa, input.abc)
par(mfrow = c(1, 1))
hist(fres.abc1$ABC, main = "distribution of ABC") # ABCの分布
ABC <- mean(fres.abc1$ABC) # 平均値をABCとする
## SSBの将来予測結果
par(mfrow = c(1, 1))
plot.future(fres.abc1, what = c(FALSE, TRUE, FALSE), is.legend = TRUE, lwd = 2,
col = "darkblue", N = 5, label = rep(NA, 3))
draw.refline(cbind(unlist(refs[c(1, 1, 2, 3) + 3]), unlist(refs[c(1, 1, 2, 3)])),
horiz = TRUE, lwd = 1, scale = 1)
## 漁獲量の将来予測結果
par(mfrow = c(1, 1))
plot.future(fres.abc1, what = c(FALSE, FALSE, TRUE), is.legend = TRUE, lwd = 2,
col = "darkblue", N = 5, label = rep(NA, 3))
points(fres.abc1$input$ABC.year, ABC, pch = 20, col = 2, cex = 3)
text(fres.abc1$input$ABC.year + 1, ABC, "ABC", col = 2)
## 実際に、どんなFが将来予測で使われているか
boxplot(t(fres.abc1$faa[1, , ]/fres.abc1$faa[1, 1, ]), ylab = "multiplier to current F")
# どんなHCRなのか書いてみる
ssb.abc <- mean(fres.abc1$vssb[rownames(fres.abc1$vssb) %in% fres.abc1$input$ABC.year,
]) # ABC計算年のssbをとる
plot.HCR(beta = beta, bban = refs$Bban, blimit = refs$Blim, btarget = refs$Bmsy,
lwd = 2, xlim = c(0, refs$Bmsy * 2), ssb.cur = ssb.abc, Fmsy = refs$Fmsy,
yscale = 0.7, scale = 1000)
# 将来の親魚資源量がBMSYやBlimitを上回る確率の表示
plot(apply(fres.abc1$vssb > refs$Bmsy, 1, mean) * 100, type = "b", ylab = "Probability",
ylim = c(0, 100))
points(apply(fres.abc1$vssb > refs$BmsyAR, 1, mean) * 100, pch = 2, type = "b")
points(apply(fres.abc1$vssb > refs$Blim, 1, mean) * 100, pch = 1, col = 2, type = "b")
points(apply(fres.abc1$vssb > refs$BlimAR, 1, mean) * 100, pch = 2, col = 2,
type = "b")
abline(h = c(50, 90), col = c(1, 2))
legend("bottomright", col = c(1, 1, 2, 2), title = "Probs", pch = c(1, 2, 1,
2), legend = c(">Btarget_Eq", ">Btarget_AR", ">Blimit_Eq", ">Blimit_AR"))
# Kobe chart
plot.kobe(res.pma, Bmsy = refs$Bmsy, Umsy = refs$Umsy, Blim = refs$Blim)
8. ABC計算までのまとめ（時間がない人はここからスタート）
MSY管理基準値を計算は以下の手順でおこないます．
- データの読み込み
# 関数の読み込み →
# warningまたは「警告」が出るかもしれませんが，その後動いていれば問題ありません
source("../../rvpa1.9.2.r")
source("../../future2.1.r")
source("../../utilities.r", encoding = "UTF-8") # ggplotを使ったグラフ作成用の関数
# データの読み込み
caa <- read.csv("caa_pma.csv", row.names = 1)
waa <- read.csv("waa_pma.csv", row.names = 1)
maa <- read.csv("maa_pma.csv", row.names = 1)
dat <- data.handler(caa = caa, waa = waa, maa = maa, M = 0.5)
names(dat)- VPAの実施(vpa) → res.pma(VPAの結果)を得る
- current Fとしてどのような値を使うか、ここで設定しておく（fc.yearオプションで、何年から何年のFを平均するか指定)
# VPAによる資源量推定
res.pma <- vpa(dat, fc.year = 2009:2011, rec = 585, rec.year = 2011, tf.year = 2008:2010,
term.F = "max", stat.tf = "mean", Pope = TRUE, tune = FALSE, p.init = 1)- 再生産関係パラメータのあてはめ (fit.SR) → HS.par0 (HSにあてはめたときのパラメータ推定結果)を得る
- 残差の自己相関がある・なしを決める。ある場合はAR=1としたときの結果を用います。
# VPA結果を使って再生産データを作る
SRdata <- get.SRdata(res.pma)
head(SRdata)HS.par0 <- fit.SR(SRdata, SR = "HS", method = "L2", AR = 0, hessian = FALSE)
HS.par1 <- fit.SR(SRdata, SR = "HS", method = "L2", AR = 1, hessian = FALSE)
BH.par0 <- fit.SR(SRdata, SR = "BH", method = "L2", AR = 0, hessian = FALSE)
BH.par1 <- fit.SR(SRdata, SR = "BH", method = "L2", AR = 1, hessian = FALSE)
RI.par0 <- fit.SR(SRdata, SR = "RI", method = "L2", AR = 0, hessian = FALSE)
RI.par1 <- fit.SR(SRdata, SR = "RI", method = "L2", AR = 1, hessian = FALSE)
c(HS.par0$AICc, HS.par1$AICc, BH.par0$AICc, BH.par1$AICc, RI.par0$AICc, RI.par1$AICc)- HS.par0をもとに将来予測を実施する(future.vpa) → fres.HS (HSを仮定したときの将来予測結果)を得る
- 生物パラメータを平均する年,ABC計算年などのオプションを設定
- 資源量と年齢別の体重に相関がある場合はそれを将来予測の設定に取り込む(waa.fun=TRUE)
fres.HS <- future.vpa(res.pma,
multi=1, # res.pma$Fc.at.ageに掛ける乗数
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
is.plot=TRUE, # 結果をプロットするかどうか
seed=1,
silent=TRUE,
recfunc=HS.recAR, # 再生産関係の関数
# recfuncに対する引数
rec.arg=list(a=HS.par0$pars$a,b=HS.par0$pars$b,
rho=HS.par0$pars$rho, # ここではrho=0なので指定しなくてもOK
sd=HS.par0$pars$sd,resid=HS.par0$resid))- res.pmaとfres.HSを使ってMSY管理基準値を計算する (est.MSY) → MSY.HS (管理基準値の推定結果)を得る
# MSY管理基準値の計算
MSY.HS <- est.MSY(res.pma, # VPAの計算結果
fres.HS$input, # 将来予測で使用した引数
# nyear=NULL, # 何年計算するかは、指定しなければ関数内部で世代時間の20倍の年数を計算し、それを平衡状態とする
resid.year=0, # ARありの場合、最近何年分の残差を平均するかをここで指定する。ARありの設定を反映させたい場合必ずここを１以上とすること（とりあえず１としておいてください）。
N=100, # 将来予測の年数，繰り返し回数
PGY=c(0.9,0.6,0.1), # 計算したいPGYレベル。上限と下限の両方が計算される
onlylower.pgy=FALSE, # TRUEにするとPGYレベルの上限は計算しない（計算時間の節約になる）
B0percent=c(0.3,0.4)) # 計算したいB0%レベル# 結果の表示(tibbleという形式で表示され、最初の10行以外は省略されます)
(refs.all <- MSY.HS$summary_tb)
# 全データをじっくり見たい場合
View(refs.all)
# のちの使用のために、Bmsy, Blimit, Bban, Fmsyを定義しておく
refs <- list(BmsyAR=as.numeric(MSY.HS$summaryAR$SSB[1]),
BlimAR=as.numeric(MSY.HS$summaryAR$SSB[6]),
BbanAR=as.numeric(MSY.HS$summaryAR$SSB[8]),
Bmsy=as.numeric(MSY.HS$summary$SSB[1]),
Blim=as.numeric(MSY.HS$summary$SSB[6]),
Bban=as.numeric(MSY.HS$summary$SSB[8]),
Fmsy=as.numeric(MSY.HS$summary$"Fref/Fcur"[1]),
MSY=as.numeric(MSY.HS$summary$Catch[1]),
Umsy=as.numeric(MSY.HS$summary$Catch[1])/as.numeric(MSY.HS$all.stat$biom.mean[1]))
# また、各管理基準値に対して以下のようなラベルをつける
# どの管理基準値をどのように定義するか、ここで指定します
refs.all$RP.definition <- NA
refs.all$RP.definition[refs.all$RP_name=="MSY" & refs.all$AR==FALSE] <- "Btarget0" # RP_nameがMSYでARがなしのものをBtargetとする
refs.all$RP.definition[refs.all$RP_name=="B0-20%" & refs.all$AR==FALSE] <- "Btarget1" # たとえばBtargetの代替値をいちおう示す場合
refs.all$RP.definition[refs.all$RP_name=="PGY_0.95_lower" & refs.all$AR==FALSE] <- "Btarget2"
refs.all$RP.definition[refs.all$RP_name=="PGY_0.9_lower" & refs.all$AR==FALSE] <- "Blow0"
refs.all$RP.definition[refs.all$RP_name=="PGY_0.6_lower" & refs.all$AR==FALSE] <- "Blimit0"
refs.all$RP.definition[refs.all$RP_name=="PGY_0.1_lower" & refs.all$AR==FALSE] <- "Bban0"
refs.all$RP.definition[refs.all$RP_name=="Ben-19431" & refs.all$AR==FALSE] <- "Bcurrent"
refs.all$RP.definition[refs.all$RP_name=="Ben-63967" & refs.all$AR==FALSE] <- "Bmax"
refs.all$RP.definition[refs.all$RP_name=="Ben-24000" & refs.all$AR==FALSE] <- "Blimit1"
refs.all$RP.definition[refs.all$RP_name=="Ben-51882" & refs.all$AR==FALSE] <- "B_HS"
# 定義した結果を見る
refs.all %>% select(RP_name,RP.definition)
# refs.allの中からRP.definitionで指定された行だけを抜き出す
(refs.base <- refs.all %>%
filter(!is.na(RP.definition)) %>% # RP.definitionがNAでないものを抽出
arrange(desc(SSB)) %>% # SSBを大きい順に並び替え
select(RP.definition,RP_name,SSB,Catch,U,Fref2Fcurrent)) # 列を並び替え- 決定されたHCRを用いて20年程度の将来予測を実施し、ABC算出年のABCを計算する
- 10年後にBtargetを上回る確率なども計算
input.abc <- MSY.HS$input$msy # MSY計算で使った引数を使う
input.abc$N <- 1000 # 実際に計算するときは10000以上を使ってください
input.abc$HCR <- list(Blim = refs$Blim, Bban = refs$Bban, beta = 0.8) # 算定ルールのデフォルトは0.8
input.abc$nyear <- 30 # ABC計算時には長期間計算する必要はない
input.abc$ABC.year <- 2013 # ここでABC.yearを設定しなおしてください
input.abc$is.plot <- TRUE
fres.abc1 <- do.call(future.vpa, input.abc)
par(mfrow = c(1, 1))
hist(fres.abc1$ABC, main = "distribution of ABC") # ABCの分布
ABC <- mean(fres.abc1$ABC) # 平均値をABCとする
## SSBの将来予測結果
par(mfrow = c(1, 1))
plot.future(fres.abc1, what = c(FALSE, TRUE, FALSE), is.legend = TRUE, lwd = 2,
col = "darkblue", N = 5, label = rep(NA, 3))
draw.refline(cbind(unlist(refs[c(1, 1, 2, 3) + 3]), unlist(refs[c(1, 1, 2, 3)])),
horiz = TRUE, lwd = 1, scale = 1)
## 漁獲量の将来予測結果
par(mfrow = c(1, 1))
plot.future(fres.abc1, what = c(FALSE, FALSE, TRUE), is.legend = TRUE, lwd = 2,
col = "darkblue", N = 5, label = rep(NA, 3))
points(fres.abc1$input$ABC.year, ABC, pch = 20, col = 2, cex = 3)
text(fres.abc1$input$ABC.year + 1, ABC, "ABC", col = 2)
## 実際に、どんなFが将来予測で使われているか
boxplot(t(fres.abc1$faa[1, , ]/fres.abc1$faa[1, 1, ]), ylab = "multiplier to current F")# どんなHCRなのか書いてみる
ssb.abc <- mean(fres.abc1$vssb[rownames(fres.abc1$vssb) %in% fres.abc1$input$ABC.year,
]) # ABC計算年のssbをとる
plot.HCR(beta = beta, bban = refs$Bban, blimit = refs$Blim, btarget = refs$Bmsy,
lwd = 2, xlim = c(0, refs$Bmsy * 2), ssb.cur = ssb.abc, Fmsy = refs$Fmsy,
yscale = 0.7, scale = 1000)# 将来の親魚資源量がBMSYやBlimitを上回る確率の表示
plot(apply(fres.abc1$vssb > refs$Bmsy, 1, mean) * 100, type = "b", ylab = "Probability",
ylim = c(0, 100))
points(apply(fres.abc1$vssb > refs$BmsyAR, 1, mean) * 100, pch = 2, type = "b")
points(apply(fres.abc1$vssb > refs$Blim, 1, mean) * 100, pch = 1, col = 2, type = "b")
points(apply(fres.abc1$vssb > refs$BlimAR, 1, mean) * 100, pch = 2, col = 2,
type = "b")
abline(h = c(50, 90), col = c(1, 2))
legend("bottomright", col = c(1, 1, 2, 2), title = "Probs", pch = c(1, 2, 1,
2), legend = c(">Btarget_Eq", ">Btarget_AR", ">Blimit_Eq", ">Blimit_AR"))