From 553b626310550038e459a85f2718992cf8e4be49 Mon Sep 17 00:00:00 2001 From: "kathryn.doering" Date: Fri, 4 Sep 2020 17:13:20 -0700 Subject: [PATCH] docs: add performance metric example, table of contents --- README.Rmd | 55 +++++++++++++++++-- README.md | 75 +++++++++++++++++++++++--- man/figures/README-plot_rel_SSB-1.png | Bin 0 -> 4776 bytes 3 files changed, 119 insertions(+), 11 deletions(-) create mode 100644 man/figures/README-plot_rel_SSB-1.png diff --git a/README.Rmd b/README.Rmd index 9b7ca2d3..0f04d05e 100644 --- a/README.Rmd +++ b/README.Rmd @@ -1,6 +1,9 @@ --- title: "README" -output: github_document +output: + github_document: + toc: true + toc_depth: 2 --- @@ -27,11 +30,12 @@ master: ************** -## This is a repository for the Stock Assessment Tool: SSMSE +**This is a repository for the Stock Assessment Tool: SSMSE** + - Supported by the NOAA Fisheries Integrated Toolbox -## Disclaimer +**Disclaimer** “The United States Department of Commerce (DOC) GitHub project code is provided on an ‘as is’ basis and the user assumes responsibility for its use. DOC has relinquished control of the information and no longer has responsibility to protect the integrity, confidentiality, or availability of the information. Any claims against the Department of Commerce stemming from the use of its GitHub project will be governed by all applicable Federal law. Any reference to specific commercial products, processes, or services by service mark, trademark, manufacturer, or otherwise, does not constitute or imply their endorsement, recommendation or favoring by the Department of Commerce. The Department of Commerce seal and logo, or the seal and logo of a DOC bureau, shall not be used in any manner to imply endorsement of any commercial product or activity by DOC or the United States Government.” @@ -69,7 +73,7 @@ Still having trouble installing SSMSE? Please don't hesitate to open an [issue]( ## An SSMSE example -Suppose we want to look at 2 scenarios, one where Steepness (H) is specified correctly and one where it is specified incorrectly in an estimation model (EM): +Suppose we want to look at how well we are able to achieve a specified management procedure under uncertainty in the operating model. We will look 2 scenarios, one where Steepness (H) is specified correctly and one where it is specified incorrectly in an estimation model (EM): 1. **H-ctl**: Cod operating model (H = 0.65) with correctly specified cod model EM (fixed H = 0.65) 2. **H-1**: Cod operating model (OM; H = 1) with misspecified cod model EM (fixed H = 0.65) @@ -200,6 +204,10 @@ sample_struct_list <- list("H-ctl" = sample_struct, "H-1" = sample_struct) Process error can be added through the recruitment deviations. In this case, `rec_dev_pattern = "rand"` in the call to `run_SSMSE` is used to use random recruitment deviations with the same standard deviation as the historical recruitment deviation pattern. Set `scope = 2` so that the same recruitment deviation patterns are used across scenarios, but different patterns are use across iterations in the same scenario. +### Performance metrics + +Quantitative performance metrics should be specified before conducting an MSE. Typically, a suite of performance metrics will be examined; however, for simplicity in this example, we will only look at what the achieved relative biomass was for the last 3 years of projection in the MSE to determine how it compares to the intended management target of 40% of unfished Spawning Stock Biomass. Note that we are only running our MSE projectsion for 6 years, but longer projections are typical in MSE analyses. + ### Run SSMSE Now, use `run_SSMSE` to run the MSE analysis loop (note this will take some time to run, ~ 20 min): @@ -226,7 +234,7 @@ run_SSMSE(scen_name_vec = c("H-ctl", "H-1"),# name of the scenario ``` Here, only 5 iterations are run per scenario; However, in a real MSE analysis, running 100+ iterations would be preferred, so that the full range of uncertainty (given observation and process errors) is visible in the results. -### Summarize results and view output +### Summarize results The function `SSMSE_summary_all` can be used to summarize the model results in a list of dataframes. @@ -268,6 +276,43 @@ ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod-1_OM", ggplot2::theme_classic() ``` +Now, calculate and plot the performance metric. +```{r plot_rel_SSB} +get_rel_SSB_avg <- function(summary, min_yr, max_yr) { + OM_vals <- unique(summary$ts$model_run) + OM_vals <- grep("_OM$", OM_vals, value = TRUE) + B_unfished <- summary$scalar %>% + filter(model_run %in% OM_vals) %>% + select(iteration, scenario,SSB_Unfished) + SSB_yr <- summary$ts %>% + filter(year >= min_yr & year <= max_yr) %>% + select(iteration, scenario, year, SpawnBio) + SSB_yr <- left_join(SSB_yr, B_unfished) %>% + mutate(Rel_SSB = SpawnBio/SSB_Unfished) %>% + group_by(iteration, scenario) %>% + summarize(avg_SSB = mean(Rel_SSB), .groups = "keep") %>% + ungroup() + SSB_yr +} +rel_SSB <- get_rel_SSB_avg(summary, min_yr = 104, max_yr = 106) + +# function to summarize data in plot +data_summary <- function(x) { + m <- mean(x) + ymin <- m-sd(x) + ymax <- m+sd(x) + return(c(y=m,ymin=ymin,ymax=ymax)) +} +ggplot(data = rel_SSB, aes(x = scenario, y = avg_SSB)) + + geom_hline(yintercept = 0.4, color = "gray") + + stat_summary(fun.data = data_summary, + position = position_dodge(width = 0.9), color = "blue") + + scale_y_continuous(limits = c(0, 0.8)) + + labs(title = "Long-term average relative SSB\n(years 101-106)", + x = "Scenario", y = "SSB/SSB_unfished") + + theme_classic() +``` + If you wish to delete the files created from this example, you can use: ```{r, eval = FALSE} unlink(run_SSMSE_dir, recursive = TRUE) diff --git a/README.md b/README.md index e66544f7..83ffe26a 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,16 @@ README ================ + - [(SSMSE) Management Strategy Evaluation for Stock Synthesis + (SS)](#ssmse-management-strategy-evaluation-for-stock-synthesis-ss) + - [Installing the SSMSE R + package](#installing-the-ssmse-r-package) + - [Troubleshooting Installation](#troubleshooting-installation) + - [An SSMSE example](#an-ssmse-example) + - [How can I contribute to SSMSE?](#how-can-i-contribute-to-ssmse) + - [Roadmap: Where is SSMSE headed + next?](#roadmap-where-is-ssmse-headed-next) + # (SSMSE) Management Strategy Evaluation for Stock Synthesis (SS) @@ -17,11 +27,11 @@ status](https://ci.appveyor.com/api/projects/status/github/nmfs-fish-tools/SSMSE ----- -## This is a repository for the Stock Assessment Tool: SSMSE +**This is a repository for the Stock Assessment Tool: SSMSE** - Supported by the NOAA Fisheries Integrated Toolbox -## Disclaimer +**Disclaimer** “The United States Department of Commerce (DOC) GitHub project code is provided on an ‘as is’ basis and the user assumes responsibility for its @@ -90,9 +100,10 @@ Still having trouble installing SSMSE? Please don’t hesitate to open an ## An SSMSE example -Suppose we want to look at 2 scenarios, one where Steepness (H) is -specified correctly and one where it is specified incorrectly in an -estimation model (EM): +Suppose we want to look at how well we are able to achieve a specified +management procedure under uncertainty in the operating model. We will +look 2 scenarios, one where Steepness (H) is specified correctly and one +where it is specified incorrectly in an estimation model (EM): 1. **H-ctl**: Cod operating model (H = 0.65) with correctly specified cod model EM (fixed H = 0.65) @@ -338,6 +349,17 @@ the historical recruitment deviation pattern. Set `scope = 2` so that the same recruitment deviation patterns are used across scenarios, but different patterns are use across iterations in the same scenario. +### Performance metrics + +Quantitative performance metrics should be specified before conducting +an MSE. Typically, a suite of performance metrics will be examined; +however, for simplicity in this example, we will only look at what the +achieved relative biomass was for the last 3 years of projection in the +MSE to determine how it compares to the intended management target of +40% of unfished Spawning Stock Biomass. Note that we are only running +our MSE projectsion for 6 years, but longer projections are typical in +MSE analyses. + ### Run SSMSE Now, use `run_SSMSE` to run the MSE analysis loop (note this will take @@ -369,7 +391,7 @@ analysis, running 100+ iterations would be preferred, so that the full range of uncertainty (given observation and process errors) is visible in the results. -### Summarize results and view output +### Summarize results The function `SSMSE_summary_all` can be used to summarize the model results in a list of dataframes. @@ -446,6 +468,47 @@ ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod-1_OM", ![](man/figures/README-plot_SSB-1.png) +Now, calculate and plot the performance metric. + +``` r +get_rel_SSB_avg <- function(summary, min_yr, max_yr) { + OM_vals <- unique(summary$ts$model_run) + OM_vals <- grep("_OM$", OM_vals, value = TRUE) + B_unfished <- summary$scalar %>% + filter(model_run %in% OM_vals) %>% + select(iteration, scenario,SSB_Unfished) + SSB_yr <- summary$ts %>% + filter(year >= min_yr & year <= max_yr) %>% + select(iteration, scenario, year, SpawnBio) + SSB_yr <- left_join(SSB_yr, B_unfished) %>% + mutate(Rel_SSB = SpawnBio/SSB_Unfished) %>% + group_by(iteration, scenario) %>% + summarize(avg_SSB = mean(Rel_SSB), .groups = "keep") %>% + ungroup() + SSB_yr +} +rel_SSB <- get_rel_SSB_avg(summary, min_yr = 104, max_yr = 106) +## Joining, by = c("iteration", "scenario") + +# function to summarize data in plot +data_summary <- function(x) { + m <- mean(x) + ymin <- m-sd(x) + ymax <- m+sd(x) + return(c(y=m,ymin=ymin,ymax=ymax)) +} +ggplot(data = rel_SSB, aes(x = scenario, y = avg_SSB)) + + geom_hline(yintercept = 0.4, color = "gray") + + stat_summary(fun.data = data_summary, + position = position_dodge(width = 0.9), color = "blue") + + scale_y_continuous(limits = c(0, 0.8)) + + labs(title = "Long-term average relative SSB\n(years 101-106)", + x = "Scenario", y = "SSB/SSB_unfished") + + theme_classic() +``` + +![](man/figures/README-plot_rel_SSB-1.png) + If you wish to delete the files created from this example, you can use: ``` r diff --git a/man/figures/README-plot_rel_SSB-1.png b/man/figures/README-plot_rel_SSB-1.png new file mode 100644 index 0000000000000000000000000000000000000000..6596678a5cd3785b2a6a6376355cf846a6200e61 GIT binary patch literal 4776 zcmds5XIN9&){YRRC{?5v3vCoEgc<|_BO*uuMGyii2nG-dJxUQLFvtLs9tkp-2mak`O=$H35`_+!NfC>si@_YCN2}J_HLy<%%9smW-FDN)13rCf}u_Pc0;Ex0WMtmQ>HBj() zB%X-GQ;K0;AOsd(LMkbtlS=6H5;~8F1QC%?A|3!RCg4E(j`uchj7EYph7yyHM&+(3SLG7)%Z99RD&E$><)Q^rGo(Jd4sjB9SJ( zwHVn}y}mO-is&e(e!X6OLXMWMI>!dA&ZRLb?ZqNOCNunz$&Nx8=mN*xt~dhY@Cj{L zn@BGa^wTn7N+_Q;>UVpQz5DBD3@vtcuhNIv*QQVMpSC#!I?QkriY*eZ?EA!3m?(p| zdUU_)31s(Y8Qrd7)EexY8}Rr!YkBhCHWeKiF(B9EKJ6&oz5jHsg^xZNthjIeT56R4 z2WPMaVumcnQ$NMzk#Z5-ufcLC`dLggCaSi96cYIk;;E!5$UFm~zA#bZ71zz1IejUT z?UKh*Rv8AVFNnEz51@hYI}fK$H};Px_z2d)n!GO%_U9em!sHbFwiQtZ=0+GFah=h` zMBk%Hz;%@TsFyxHc4EsKKf!fHmn-JSggC7tImNC9-+v|+CoR3WhAIby!MuROl z-BS*ITh2e8m|l?Rdc&O}Maa)loW((2-=SoS#{K-WVO%ctX4wq#Ima5CS);=R^PChj0=P?sAD3Vw0_MJt(6KGB>X>g&bIO~PO;tQTjMV%gM?3kvo^7Ha_& z*MCbX0*_MARL=NB{V7~#RM!5p)^9zJk%w=q?t$*eBgZ}}>rJ4U+M@kNW5_Z3W7;|LMO7JWU?=fbah;`Kw~9j2n;kBmbX-Z{3`wcUZ+E{6 zX^Z20$C0B{Xr=CXjvUf8H}h<9G~Jtc!%dw$bv(#4)VOUAH0)`6+!cw-&)6)DIg*L8 zQZuP&Sf}^VzFCp_R4zIh)9%fl+O#uI?4c;jm#)$Jr*NH^uTN?&5KvlBC$&8{aX<4S zV&wgoF(v{xkhSX}ACFQk7`*ZZ@u1yzwQO4alJnGa+C+_1_X`$yGb&>jdi1`^Oc|Y= zojA5*vN;^Fa-^ZroBYJ1`;&bg^MLXx&-O=NU2nuIoVQD>9*ub}y${iF6cBM2QoOaQ zd)x0iU&wqtEujZu7-Yg~>@-_b(%G8HaV^Q9ll}Ua9(UivO z{FRDFn+trO`*K0=bgVl3F)*eEY$y36dgdg0sCgxQAFn$#SVH>w>zD?Ub)hc}qsJ2N z@_Z=N@2GY}8vc2ed_#dxjaZ&9ThviOa~s1Svi#b5P}?&hIq zB=SFf=%4Xia)s^%kM1j#A1x2uYQD0+vF^ZBuH5Z74z8JGEfC~#6?feDX6O)N#M6R2 znQ|yb@Zr$;#a~!cbq5-C;Oj1&e&asDATcH7y#rT|M zwbbZ`6B*m<33<$*NR zxH8q?RBhE^x4XD2GIFi*afQwtAoEx2X5k||En{(nd1{_**igA8DYJ0Rea^G?N?Sv8 z+1jyD6S-^VsaLyviBgyEvX%Cdr4g)Xr(lUW=F4jxI-U{;E6)oO0btSpz>vh_!g8zs zQvQwT9Ja7*Ks>&9Y4ZY&+OB-$hnLT0(UBF)c@sv?+~;jeQRJhT>(Bba7c5LeYj9oh zZ@1|m;tT@8_-fYJt#s-eNIiI*C8_PIZC2Is>Lhu1aBe^f#$dR|qb-clCzi`o+7vtA ziYa=eJE-jPMU7tNpvvTWq-6pj>&8XyKV3Ccyp5Zc-{qSIvBjSy{NKWL3>N%)d-uV; zf0K?scO{Ck3^B4+mV$LbueMz#8NxDsDGGB;iVQ*$FmBoR*d^8*pymVmPU4UQ*4~0M zu)C-JlB%Wr+v&`pL*$3;6p9r=<55*rBJ8|oaPVGCx;>RDE1;WMTc-;~=taeN2v_*7 zZ4hOQob1p1At$Z`*mg|XtG%X?kpu5yPPYlq9djUDw%K0#a@X`r*hJkvga$bWP;vdV zP@Ai^Ck=Xp>99}%rg?S9wj+)?kRR{D`85d^{>m7W^(0zyuHbI=!b>RPp$ATRQLg7`y?oea=1mALAuF-SD00~sD1 z`|-@`qh67UGr}^XYeQO?bWx0zXZM;qx$ltExg-Az9YKV;Agzdn7SAJin`*HE|4jWD zis?EIMx~pO4|W-zCkp^EI>HH4N9)ab!gP44U?c^p7VuS58^iFZ7rfFfEDVg(ahsqU zynABQxTaOu>xFc+fl>VkRYCdEtlBP}A(FRdcwusrJ=v0gPd_bw({)u@C=oBxX8C;6 zmf)r#4zX0kSRoEqIixLE{;ZKld3l6DJMx1@c1ZQO!CBDxbkBHc?Pm3`?9iX$J_vhN zzY7rWGwuZQe{#uVi=Pfd23No-E_PQX5C)rFmG4w#^v+@vvV7Q{C?OS}p@cyZ8LNjv zpBI=Pe|GQfF%@uo)0=8tjCz#KAS5)|^S8UK>8I-bc}!m|5k0qfD@2&c#~HNb_kez8 za5~I%ZYo{koR5}B-#e!Kw>rEhQ6FV`pI3+;p-Lb2yTZ25&=^B z9TneUHZ~PKI{y#VqJT9{+{$lJ39xeUdCW$B?z{0et&N4Ho_m%AwYzj`DF%Y1@vWD|M~=c1a1RaAEAT;tH>nqb@o%hyVg9>dC2odt3Y1l-I4CUqmke$ZouSex;mj$6V4UBsCw_2!$xUL(JZby`H+GWQwb zq>inpfzpqoHD9323#pVR!J*eq`#w-c>TNcMKNW#K#Ed3tVXSPXuS@arzf2Yc6OL>+ z9zx}ez{x}iH~->6rDoZi)D+rPFgtz>UaqbjvdnaC^`HHUVa79w*fY96P^b>&F}CXs~ay1!htBx0RsE zclxV)`{S)U{fR{Wd-V&NV{1og;7|u;p6_`DkT%o8%)ve=T3F`T<2=fY>(d+h26q7# zU(mID4kWwHP`eOLNu5wxYFU0GhhsmDQXbHOnuw=FPW{ic0`H&>w@SBn`+w}9ytzY& zDz9FWBYgK?#MUAvlL-4GnN`5Jl!Ain!h#oTigz73C&JdcWlobjTI)PH208~YqP9CX zYHPP=V*$x!&1`(J$n*jh-7_+yvaf825G^b5F(;$~>fWkskNJct1++SY0{ zA;N(`$H^<)dn5j1w?zid)oL*wO$tZcb!)+$uhr7A97=_?r3EFR)$Z960<8$N*3D_$ zYwU2Uklr8f@7u#jW`GIS*(4{N(H*$n?H&fsdNTM*u9SjC4|1WR_hAq{wu_{HZS(9^cg?T)w&m8haEXze0JL+=M50LP?cOF>xw7`)%(9gW7DJzyg3`&{;JmqXua zVSNZXyMK|?CIR^p7