## Global options library(rmarkdown) library(knitr) options(max.print="75") opts_chunk$set(prompt=FALSE, tidy=TRUE, comment=NA, message=FALSE, warning=FALSE)
est.MSY以降の計算についてはtidyverseパッケージのdplyrやggplot2ライブラリを多用しています。今のRのコーディングはこれを使ってやるのが主流のようです。
ここで必要なオブジェクトは以下です。
(レポート記述内容例)
library(frasyr) options(scipen=100) # 桁数表示の調整(1E+9とかを抑制する) # 再生産関係のプロット (g1_SRplot <- SRplot_gg(SRmodel.base))
(レポート記述内容例) (あくまで例です。今の例で、代替基準値を最大限選ぶとしたらどうするか、というものです)
表:さまざまな管理基準値
|ラベル | 管理基準値 | 説明 | |:-----------|:-----------------|:---------------------------------| | Btarget0 |目標 | 最大の平均漁獲量を得る時の親魚量(Bmsy)。過去最大親魚量の2倍となり、SSB>SSB_maxの範囲における不確実性が大きい懸念がある。 | | Btarget1 |目標(代替値候補1) | 説明を書く | | Btarget2|目標(代替値候補2) | 説明を書く| | Blimit0 |限界 |MSYの60%の平均漁獲量を得るときの親魚資源量 | | Blimit1 |限界(代替値候補1) |説明を書く| | Bban0 | 禁漁 |MSYの10%の平衡漁獲量を得るときの親魚資源量| | Bmax | 経験値 |過去最大親魚量 | | B_HS | 経験値 |HS再生産関係の折れ点 | | B_current |経験値 | 最近年の親魚量 |
# 管理基準値表 make_RP_table(refs.base) # 漁獲量曲線 # 再生産関係をもとにしたyield curveと管理基準値のプロット。 # 計算した全管理基準値を示す場合にはrefs.allを、厳選したものだけを示す場合にはrefs.baseを引数に使ってください # AR==TRUEにするとARありの結果もプロットされます # 将来予測と過去の漁獲量を追記した図 (g2_yield_curve <- plot_yield(MSY.base$trace,refs.base, future=list(future.Fcurrent,future.default), past=res.pma,AR=FALSE,xlim.scale=0.4,ylim.scale=1.3,lining=FALSE)) # xlimやylimを変更する場合 (g2.2 <- plot_yield(MSY.base$trace,refs.base,AR_select=FALSE,xlim.scale=0.5,ylim.scale=1.3,lining=FALSE)) # yield curveの元データが欲しい場合 yield.table <- get.trace(MSY.base$trace) yield.table <- yield.table %>% mutate(age=as.character(age)) %>% spread(key=age,value=value) %>% arrange(ssb.mean) # 神戸チャート # Btarget0として選ばれた管理基準値をベースにした神戸チャート4区分 # roll_meanで各年の値を何年分移動平均するか指定します (g3_kobe4 <- plot_kobe_gg(res.pma,refs.base,roll_mean=3,category=4, Blow="Btarget0", # Btargeと同じ値を入れておいてください Btarget="Btarget0")) # <- どの管理基準値を軸に使うのか指定。指定しなければ"0"マークがついた管理基準値が使われます # Btarget0, Blow0, Blimit0として選ばれた管理基準値をベースにした神戸チャート6区分 # Blowを使うかどうかは不明。とりあえず6区分の一番上の境界(Blowのオプション)は"Btarget0"と、targetで使う管理基準値の名前を入れて下さい (g4_kobe6 <- plot_kobe_gg(res.pma,refs.base,roll_mean=3,category=6,Blow="Btarget0"))
(レポート記述内容例) - 現状のFは2015年から2017年の年齢別Fの単純平均を用いた(など、どのようにcurrent Fを定義したかを書く)
# 親魚資源量と漁獲量の時系列の図示 (g5_future <- plot_futures(res.pma, #vpaの結果 list(future.Fcurrent,future.default), # 将来予測結果 future.name=c("F current",str_c("HCR(beta=",future.default$input$HCR$beta,")")), CI_range=c(0.1,0.9), maxyear=2045,ncol=2, Btarget=derive_RP_value(refs.base,"Btarget0")$SSB, Blimit=derive_RP_value(refs.base,"Blimit0")$SSB, # Blow=derive_RP_value(refs.base,"Blow0")$SSB, Bban=derive_RP_value(refs.base,"Bban0")$SSB, biomass.unit=10000, # バイオマスの単位(100, 1000, or 10000トン) font.size=12)) # フォントサイズ #(g5 <- g5+ggtitle("図5. 現行のFとデフォルトのHCRを用いた時の将来予測\n(実線:平均値、範囲:80パーセント信頼区間)")+ylab("トン")) # Fcurrentの図 (g6_Fcurrent <- plot_Fcurrent(res.pma,year.range=2007:2017)) # HCRの図 (g7_hcr <- plot_HCR(SBtarget=derive_RP_value(refs.base,"Btarget0")$SSB, SBlim=derive_RP_value(refs.base,"Blimit0")$SSB, SBban=derive_RP_value(refs.base,"Bban0")$SSB, Ftarget=1,biomass.unit=1000, beta=0.8)) ggsave("g7_hcr.png",g7_hcr,width=8,height=4,dpi=600)
ggsave("g1_SRplot.png",g1_SRplot,width=8,height=5,dpi=400) ggsave("g2_yield_curve.png",g2_yield_curve,width=8,height=5,dpi=400) ggsave("g3_kobe4.png",g3_kobe4,width=8,height=5,dpi=400) ggsave("g4_kobe6.png",g4_kobe6,width=8,height=5,dpi=400) ggsave("g5_future.png",g5_future,width=12,height=6,dpi=400) ggsave("g6_Fcurrent.png",g6_Fcurrent,width=6,height=4,dpi=400)
calc_kobeII_matrix
で計算した結果を使いますall.table <- bind_rows(catch.table, ssbtarget.table, ssblow.table, ssblimit.table, ssbmin.table) write.csv(all.table,file="all.table.csv")
library(formattable) catch.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","steelblue"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
library(formattable) Fsakugen.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","steelblue"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
ssbtarget.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","olivedrab"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
ssblimit.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","olivedrab"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
ssblimit.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","olivedrab"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
ssbmin.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","olivedrab"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
catch.aav.table %>% select(-stat_name) %>% formattable::formattable(list(area(col=-1)~color_tile("white","olivedrab"), beta=color_tile("white","blue"), HCR_name=formatter("span", style = ~ style(color = ifelse(HCR_name == "Btarget0-Blimit0-Bban0" & beta==0.8, "red", "black")))))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.