knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=10, fig.height=10 )
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最小モデルを今後使っていく
waa.year
, maa.year
, M.year
)。ABC計算年(ABC.year
)などの設定もここで。rec.fun
に関数名を、rec.arg
にリスト形式で引数を与えます。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))
nyear
で指定します)資源量やそれに対応するF等を管理基準値として算出しますres_future_Fcurrent
をMSY計算では使っていきます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)
# どの管理基準値をどのように定義するか。デフォルトから外れる場合はここで定義する 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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.