devtools::load_all()
data(res_vpa)
# 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), ]) # 結果
# HSのうちR0が低いケース(12番)とAIC最小ケースとの比較
plot(SR.list[[1]]$pred,type="l",ylim=c(0,2000))
points(SR.list[[12]]$pred,type="l",col=2)
SRmodel.base <- SR.list[[1]] # AIC最小モデルを今後使っていく
SRmodel.R1 <- SR.list[[12]] # 候補となる別の加入シナリオ
future.Fcurrent <- future.vpa(res_vpa,
multi=1,
nyear=50, # 将来予測の年数
start.year=2018, # 将来予測の開始年
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))
future.Fcurrent_R1 <- future.vpa(res_vpa,
multi=1,
nyear=50, # 将来予測の年数
start.year=2018, # 将来予測の開始年
N=1000, # 確率的計算の繰り返し回数=>実際の計算では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.R1$pars$a,b=SRmodel.R1$pars$b,
rho=SRmodel.R1$pars$rho, # ここではrho=0なので指定しなくてもOK
sd=SRmodel.R1$pars$sd,resid=SRmodel.R1$resid))
plot.futures(list(future.Fcurrent,future.Fcurrent_R1))
# MSY管理基準値の計算; base caseのシナリオをもとにした管理基準値
MSY.base <- est.MSY(res_vpa, # VPAの計算結果
future.Fcurrent$input, # 将来予測で使用した引数
resid.year=0, # ARありの場合、最近何年分の残差を平均するかをここで指定する。ARありの設定を反映させたい場合必ずここを1以上とすること(とりあえず1としておいてください)。
N=100, # 確率的計算の繰り返し回数=>実際の計算では1000~5000回くらいやってください
calc.yieldcurve=TRUE,
PGY=c(0.6,0.1), # 計算したいPGYレベル。上限と下限の両方が計算される
onlylower.pgy=TRUE, # TRUEにするとPGYレベルの上限は計算しない(計算時間の節約になる)
B0percent=NULL,
Bempirical=NULL
) # 計算したいB0%レベル
(refs.all <- MSY.base$summary_tb)
# refs.allの中からRP.definitionで指定された行だけを抜き出す
(refs.base <- refs.all %>%
dplyr::filter(!is.na(RP.definition)) %>% # RP.definitionがNAでないものを抽出
arrange(desc(SSB)) %>% # SSBを大きい順に並び替え
select(RP.definition,RP_name,SSB,SSB2SSB0,Catch,Catch.CV,U,Fref2Fcurrent)) # 列を並び替え
# HCRによる将来予測
input.abc <- future.Fcurrent$input # Fcurrentにおける将来予測の引数をベースに将来予測します
input.abc$multi <- derive_RP_value(refs.base,"Btarget0")$Fref2Fcurrent # currentFへの乗数を"Btarget0"で指定した値に
input.abc$silent <- TRUE
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) # BlimitはBlimit0, BbanはBban0の値
input.abc$N <- 1000
future.default <- do.call(future.vpa,input.abc) # デフォルトルールの結果→図示などに使う
# MSEによる将来予測
source("future-diff.r")
input.mse <- input.abc
input.mse$N <- 100
input.mse$use.MSE <- TRUE # MSE仕様での将来予測
input.mse$is.plot <- FALSE
future.mse <- do.call(future.vpa,input.mse)
# use.MSEオプションで、ARありバージョンには十分対応していない→今後の課題
input.mse_R1 <- input.mse
input.mse_R1$N <- 100
input.mse_R1$MSE.options$recfunc <- future.Fcurrent_R1$recfunc
future.Fcurrent_R1$rec.arg$sd2 <- future.Fcurrent_R1$rec.arg$sd2[1:input.mse_R1$N]
input.mse_R1$MSE.options$rec.arg <- future.Fcurrent_R1$rec.arg
future.mse_R1 <- do.call(future.vpa,input.mse_R1)
plot_futures(res_vpa,list(future.default,future.mse,future.mse_R1),
future.name=c("default","mse","mse_R1"))
# 直近の漁獲量の比較
all.table <- purrr::map_dfr(list(future.mse,future.default,future.mse_R1),convert_future_table,.id="scenario")
all.table %>% dplyr::filter(stat=="catch",year<2025,year>2018) %>%
ggplot() +
geom_boxplot(aes(x=factor(year),y=value,fill=scenario))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.