-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathYFS_MCMC.Rmd
91 lines (65 loc) · 4.19 KB
/
YFS_MCMC.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
---
title: "YFS flatfish mcmc"
author: "SSMA Authors"
date: "10/9/2018"
output: word_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
source("../R/prelims.R")
mytheme = .THEME
#---------------------------------------------------------------
# Read in MCMC header and results
hdr <- read.table("/Users/ingridspies/admbmodels/flatfish/assessments/yfs/runs/data/hdr.dat",as.is=T,header=F)
mctmp2.df <-read_table2("/Users/ingridspies/admbmodels/flatfish/assessments/yfs/runs/m2/evalout.rep",col_names=FALSE)
mctmp1.df <-read_table2("/Users/ingridspies/admbmodels/flatfish/assessments/yfs/runs/m1/evalout.rep",col_names=FALSE)
names(mctmp1.df) <- hdr
names(mctmp2.df) <- hdr
mcdat=rbind(mctmp1.df,mctmp2.df)
mcdat$Model=c(rep("Model 18.1",nrow(mctmp1.df)),rep("Model 18.1a",nrow(mctmp1.df)))
p1=ggplot(data=mcdat, aes(x=Fmsyr,fill=Model)) + geom_density(aes(x=Fmsyr),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("Fmsy") + xlim(c(0,.3))+theme_bw()+ylab("Density")
p2=ggplot(data=mcdat, aes(x=1000*Bmsy,fill=Model)) + geom_density(aes(x=Bmsy),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("Bmsy (x1,000 t)") +theme_bw()+ylab("Density")+xlim(c(100,1300))
p3=ggplot(data=mcdat, aes(x=ln_mean_R,fill=Model)) + geom_density(aes(x=ln_mean_R),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("log(mean(Recruitment))x10^9") +theme_bw()+ylab("Density")
p4=ggplot(data=mcdat, aes(x=ABC_biom6,fill=Model)) + geom_density(aes(x=ABC_biom6),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("Age 6 Biomass (x1,000 t)") +theme_bw()+ylab("Density")+xlim(c(1000,7000))
p5=ggplot(data=mcdat, aes(x=SSB_2019,fill=Model)) + geom_density(aes(x=SSB_2019),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("Female Spawning Biomass (x1,000 t)") +theme_bw()+ylab("Density")+geom_vline(xintercept=FSB_this/1000)+xlim(c(300,1400))
p6=ggplot(data=mcdat, aes(x=1000*TotBiom_2019,fill=Model)) + geom_density(aes(x=1000*TotBiom_2019),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("Total Biomass (x1,000 t)") +theme_bw()+ylab("Density")
mctmp1.df$Model <- "Model 18.1a"
mctmp.df$ModNum <- 1
mctmp.df <- mctmp.df %>% mutate(Recruit_NLL=rec_like_1+rec_like_2+rec_like_3+rec_like_4, Selectivity_NLL=sel_like1+sel_like2+sel_like3)
mc.df <- rbind(mc.df,mctmp.df)
ggplot(data=mctmp1.df, aes(x=Fmsyr,fill=Model)) + geom_density(aes(x=Fmsyr),alpha=.4)+scale_y_continuous(breaks=NULL) + xlab("Fmsy") + xlim(c(0,.3))+theme_bw()+ylab("Density")
mc.df %>% select(Model,Fmsyr) %>%
ggplot(aes(x=Fmsyr,fill=Model)) + geom_density(alpha=.4) + mytheme+ scale_y_continuous(breaks=NULL) + xlab("Fmsy") + xlim(c(0,.3))
mc.df <-tibble()
mods <- c("Base","Model 18.1","FixQ", "Model 14.2")
i=2
for (i in c(2,4)){
mctmp.df <-read_table2(paste0("arc/mod",i,"_evalout.rep"),col_names=FALSE)
names(mctmp.df) <- hdr
mctmp.df$Model <- M18.1a
mctmp.df$ModNum <- 1
mctmp.df <- mctmp.df %>% mutate(Recruit_NLL=rec_like_1+rec_like_2+rec_like_3+rec_like_4, Selectivity_NLL=sel_like1+sel_like2+sel_like3)
mc.df <- rbind(mc.df,mctmp.df)
}
#---------------------------------------------------------------
```
## Marginals from MCMC runs
Model 2 (called Model 18_1 in SAFE) and Model 4 (called Model 14_2 in SAFE) need posterior of Fmsy, maybe on one figure?.
For Model 18.1 (new base) need posteriors of Fmsy Bmsy, ABC biomass 2018, ABC biomass 2019, FSB and mean log recruitment.
```{r fmsyr, echo=FALSE}
# All Fmsy
mc.df %>% select(Model,Fmsyr) %>%
ggplot(aes(x=Fmsyr,fill=Model)) + geom_density(alpha=.4) + mytheme+ scale_y_continuous(breaks=NULL) + xlab("Fmsy") + xlim(c(0,.3))
```
## Other plots
You can also embed plots, for example:
```{r Biomass, echo=FALSE}
# Biomass
mc.df %>% select(Model,Fmsyr,ABC_biom1) %>% ggplot(aes(x=ABC_biom1,fill=Model)) + geom_density(alpha=.4) + mytheme+ scale_y_continuous(breaks=NULL) + xlab("Biomass (kt)")
```
```{r catch, echo=FALSE}
# Catch
mc.df %>% select(Model,Fmsyr,ABC_biom1) %>% mutate(Catch = Fmsyr*ABC_biom1) %>% ggplot(aes(x=Catch,fill=Model)) + geom_density(alpha=.4) + mytheme+ scale_y_continuous(breaks=NULL) + xlab("2019 Fmsy x Biomass")
ABC.df <- mc.df %>% select(Model,Fmsyr,ABC_biom1) %>% summarise(FABC=1/mean(1/Fmsyr),Biomass=exp(mean(log(ABC_biom1))),ABC=FABC*Biomass)
kable(ABC.df,format="html")
```