In ichimomo/frasyr: Fisheries Research Agency (FRA) provides the method for calculating sustainable yield (SY) with R

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=10,
  fig.height=10
)

res_vpaの読み込み

• 設定ポイント: vpa関数の引数fc.yearで指定した年数が今後current FのFとして扱われます。

library(frasyr)
library(tidyverse)

data(res_vpa)
res_vpa$Fc.at.age # 将来予測やMSY計算で使うcurrent Fを確認してプロットする
plot(res_vpa$Fc.at.age,type="b",xlab="Age",ylab="F",ylim=c(0,max(res_vpa$Fc.at.age)))

# 独自のFc.at.ageを使いたい場合は以下のようにここで指定する
# res_vpa$Fc.at.age[] <- c(1,1,2,2)

再生産関係の推定

# VPA結果を使って再生産データを作る
SRdata <- get.SRdata(res_vpa, years=1988:2016) 
head(SRdata)

## モデルのフィット(網羅的に試しています)
# 網羅的なパラメータ設定
SRmodel.list <- expand.grid(SR.rel = c("HS","BH","RI"), AR.type = c(0, 1), L.type = c("L1", "L2"))
SR.list <- list()
for (i in 1:nrow(SRmodel.list)) {
    SR.list[[i]] <- fit.SR(SRdata, SR = SRmodel.list$SR.rel[i], method = SRmodel.list$L.type[i], 
        AR = SRmodel.list$AR.type[i], hessian = FALSE)
}

SRmodel.list$AICc <- sapply(SR.list, function(x) x$AICc)
SRmodel.list$delta.AIC <- SRmodel.list$AICc - min(SRmodel.list$AICc)
SR.list <- SR.list[order(SRmodel.list$AICc)]  # AICの小さい順に並べたもの
(SRmodel.list <- SRmodel.list[order(SRmodel.list$AICc), ]) # 結果

SRmodel.base <- SR.list[[1]] # AIC最小モデルを今後使っていく

currentFによる将来予測

• 細かい設定の解説はこちら
  • 自己相関を考慮する場合
  • Frecオプション（目標の年に指定した確率で漁獲する）
  • 年齢別体重が資源尾数に影響される場合、などのオプションがあります
• 設定ポイント:　将来予測やMSY推定で使う生物パラメータをここで指定します（waa.year, maa.year, M.year）。ABC計算年（ABC.year）などの設定もここで。
• 設定ポイント:　再生産関係の関数型とパラメータも与えます。rec.funに関数名を、rec.argにリスト形式で引数を与えます。
• これはFcurrentでの将来予測を実施しますが、今後の管理基準値計算でもここで指定したオプションを引き継いで使っていきます
• 近年の加入の仮定(rec.new)や近年の漁獲量(pre.catch)を設定する場合にはここで設定してください
• 引数 silent == TRUE とすると、設定した引数のリストがすべて表示されます。意図しない設定などがないかどうか確認してください。

res_future_Fcurrent <- future.vpa(res_vpa,
                      multi=1,
                      nyear=50, # 将来予測の年数
                      start.year=2011, # 将来予測の開始年
                      N=100, # 確率的計算の繰り返し回数=>実際の計算では1000~5000回くらいやってください
                      ABC.year=2019, # ABCを計算する年
                      waa.year=2015:2017, # 生物パラメータの参照年
                      maa.year=2015:2017,
                      M.year=2015:2017,
                      is.plot=TRUE, # 結果をプロットするかどうか
                      seed=1,
                      silent=FALSE,
                      recfunc=HS.recAR, # 再生産関係の関数
                      # recfuncに対する引数
                      rec.arg=list(a=SRmodel.base$pars$a,b=SRmodel.base$pars$b,
                                   rho=SRmodel.base$pars$rho, # ここではrho=0なので指定しなくてもOK
                                   sd=SRmodel.base$pars$sd,resid=SRmodel.base$resid))

MSY管理基準値の計算

• MSY管理基準値計算では，上記の将来予測において，Fcurrentの値に様々な乗数を乗じたF一定方策における平衡状態時の（世代時間×20年をnyearで指定します）資源量やそれに対応するF等を管理基準値として算出します
• なので、ここまでのプロセスで、ABC計算のためにきちんとしたオプションを設定したfuture.vpaを実行しておいてください。その返り値res_future_FcurrentをMSY計算では使っていきます
• MSY.est関数の引数の詳細な解説はこちら
• オプションPGY(MSYに対する比率を指定) やB0percent(B0に対する比率を指定)、Bempirical(親魚資源量の絶対値で指定)で、別の管理基準値も同時に計算できます。
• 最近年の親魚量で維持した場合の管理基準値も、比較のためにあとで見るためBempiricalで指定しておいてください。また、B_HS(HSの折れ点)や最大親魚量などもここで計算しておいても良いかと。。。

# MSY管理基準値の計算
res_MSY <- est.MSY(res_vpa, # VPAの計算結果
                 res_future_Fcurrent$input, # 将来予測で使用した引数
                 N=100, # 確率的計算の繰り返し回数=>実際の計算では1000~5000回くらいやってください
                 calc.yieldcurve=TRUE,
                 PGY=c(0.95,0.9,0.6,0.1), # 計算したいPGYレベル。上限と下限の両方が計算される
                 onlylower.pgy=FALSE, # TRUEにするとPGYレベルの上限は計算しない（計算時間の節約になる）
                 B0percent=c(0.2,0.3,0.4),
                 Bempirical=c(round(tail(colSums(res_vpa$ssb),n=1)),
                              round(max(colSums(res_vpa$ssb))),
                              24000, # 現行Blimit
                              SRmodel.base$pars$b) # HSの折れ点
                 ) # 計算したいB0%レベル

結果の表示

• res_MSY$summaryにすべての結果が入っています。

# 結果の表示(tibbleという形式で表示され、最初の10行以外は省略されます)
options(tibble.width = Inf)
(refs.all <- res_MSY$summary)

# 全データをじっくり見たい場合
# View(refs.all)

管理基準値の選択

• 設定ポイント est.MSYで計算された管理基準値から、何をBtarget, Blimit, Bbanとして用いるかをチョイスします。
• 具体的には、refs.allにRP.definitionという新しい列をひとつ作って、その列にそれぞれの管理基準値をどのように使うかを指定します
• 「管理基準値名 + 0」はデフォルト規則による管理基準値
• デフォルトでは、ARなし、MSY="Btarget0", 0.9MSY="Blow0",0.6MSY="Blimit0", 0.1MSY="Bban0"になるようになっています
• 代替候補がある場合は「管理基準値名 + 数字」として指定
• たとえば目標管理基準値の第一候補はBmsyなのでRP_nameがMSYでARなしの行のRP.definitionには"Btarget0"と入力します
• Rコードがちょっと汚いですがご容赦ください。いい方法あったら教えてください。

# どの管理基準値をどのように定義するか。デフォルトから外れる場合はここで定義する
refs.all$RP.definition[refs.all$RP_name=="B0-20%" & refs.all$AR==FALSE] <- "Btarget1"  # たとえばBtargetの代替値としてB020%も候補に残しておきたい場合
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=="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 %>% dplyr::select(RP_name,RP.definition)

# refs.allの中からRP.definitionで指定された行だけを抜き出す
(refs.base <- refs.all %>%
    dplyr::filter(!is.na(RP.definition)) %>% # RP.definitionがNAでないものを抽出
    arrange(desc(SSB)) %>% # SSBを大きい順に並び替え
    dplyr::select(RP.definition,RP_name,SSB,SSB2SSB0,Catch,Catch.CV,U,Fref2Fcurrent)) #　列を並び替え

デフォルトルールを使った将来予測

# デフォルトのHCRはBtarget0,Blimit0,Bban0のセットになるので、それを使って将来予測する
input.abc <- res_future_Fcurrent$input # Fcurrentにおける将来予測の引数をベースに将来予測します
input.abc$multi <- derive_RP_value(refs.base,"Btarget0")$Fref2Fcurrent # currentFへの乗数を"Btarget0"で指定した値に
input.abc$HCR <- list(Blim=derive_RP_value(refs.base,"Blimit0")$SSB,
                      Bban=derive_RP_value(refs.base,"Bban0")$SSB,
                      beta=0.8,
                      year.lag=0)
future.default <- do.call(future.vpa,input.abc) # デフォルトルールの結果→図示などに使う

## 網羅的将来予測の実施
# default
kobeII.data <- beta.simulation(input.abc,
                               beta_vector=seq(from=0.1,to=1,by=0.1),
                               year.lag=0)

kobeII.table <- make_kobeII_table(kobeII.data,res_vpa,
                                 Btarget=derive_RP_value(res_MSY$summary,"Btarget0")$SSB,
                                 Blimit=derive_RP_value(res_MSY$summary,"Btarget0")$SSB,
                                 Bban=derive_RP_value(res_MSY$summary,"Btarget0")$SSB,
                 year.catch=2019:2030,
                 year.ssb  =2019:2030,
                 year.Fsakugen  =2019:2030,
                 year.ssbtarget  =2019:2030,
                 year.ssblimit  =2019:2030,
                 year.ssbban  =2019:2030,
                 year.ssbmax  =2019:2030,
                 year.ssbmin  =2019:2030,
                 year.aav  =2019:2030)

## csvファイルに一括して出力する場合
out.vpa(res=res_vpa,
        msyres=res_MSY,
        fres_current=res_future_current,
        fres_HCR=res_future_0.8HCR,
        kobeII=kobeII.table,
        filename="vpa")

ichimomo/frasyr documentation built on May 3, 2024, 1:30 a.m.