MSE visualisation
PL
03/11/2018
Processing data
Calculating performance measures
The results of simulations are processed to calculate four management measures, that are stored in a data frame called ‘results’.
The limit reference point (LRP) is chosed to be 20% of unfished SSB (spawning stock biomass).
The first measure ‘Kobe Green’ is the probability that in the future the stock will be in the green Kobe quadrand. That is that both SSB is above SSB_MSY and F is below F_MSY.
The second measure ‘Catch’ is the probability that catch is above 80% of Catch_MSY.
The third measure ‘Safety’ is the probability that the stock is above the LRP (>20% SSB_Virgin).
The fourth measure is ‘Stability’ which is represented by proximity to a 100% or by [100% - coefficient of variation (CV)], this is to make all measures comparable - ideally all four measures should be close to a 100%.
###################################################################################
## Analyse results for all mps and all OMs
##set mp
mp=pmb
name_mp<-"PMB"
##adjust the dim of results accordingly
results<-data.frame(row.names = c(1:(length(unique(mp$scen))*length(unique(mp$btrig))*length(unique(mp$ftar)))))
results<-mutate(results, OM =NA, MP = NA, Kobe = NA, Catch = NA, LRP_Virgin = NA)
i<-0
for (om in unique(mp$scen))
{
for (tar in unique(mp$ftar))
{
for (trigger in unique(mp$btrig))
{
i<-i+1
MP<-mp %>% filter (ftar == tar) %>% filter (btrig == trigger) %>% filter (scen == om) %>% select (year,iter,ssb,msy_ssb,catch, msy_yield, fbar, msy_harvest, virgin_ssb) %>% mutate(rel_ssb= ssb/msy_ssb) %>% mutate (rel_f= fbar/msy_harvest)%>% mutate(Below_LRP_MSY=(ssb<0.5*msy_ssb)) %>% mutate(Kobe_Green=(rel_ssb > 1&rel_f>1) ) %>% mutate(Target_Catch= catch >0.8*msy_yield) %>% mutate(Above_LRP_Virgin=(ssb>0.2*virgin_ssb)) %>% filter (year %in% c(2019:2038))
results$OM[i]<-paste("OM",om)
results$MP[i]<-paste("M: F_tar=",tar," B_trig=", trigger, sep="")
#Probability of being in the green Kobe quadrant
results$Kobe[i]<-sum(MP$Kobe_Green)*100/length(MP$iter)
#Probability of being above 80% of MSY for catch
results$Catch[i]<-sum(MP$Target_Catch)*100/length(MP$iter)
#Probability of being above LRP
results$LRP_Virgin[i]<-sum(MP$Above_LRP_Virgin)*100/length(MP$iter)
#CV of catch
results$Catch_Var[i]<-sd(MP$catch)*100/mean(MP$catch)
}
}
}
results<-mutate(results, Catch_Var = 100-Catch_Var)
dat=melt(results,id=c("MP","OM"))
## rename the performance measures
dat$variable=factor(dat$variable, levels=c("Kobe","Catch","LRP_Virgin","Catch_Var"),
labels=c("Kobe Green","Catch","Safety","Stability"))
## save for later
dat1=dat
##############
## make a table
file_name<-paste("Table_",name_mp,".png",sep ="")
png(file=file_name, width = 800, height = 1000, units = "px")
p1<-ggplot(aes(variable, value),data=dat)+geom_col(aes(fill = value))+
facet_grid(OM~MP,scale="free_y")+scale_fill_gradient2(low ="darkred", high ="yellow", mid = "red", midpoint = 50, limit = c(0,100),name = "Performance")+
ylab("Acceptability")+xlab("Performance Measure")+
ggtitle(paste("Management Procedure:",name_mp, sep ="")) +
labs(subtitle = "MP with 2 options for F target and B trigger")+
theme_light()+theme(axis.text.x=element_text(angle=45, hjust=1),plot.title = element_text(lineheight=.8, face="bold", colour="darkblue", size=12))
p1
dev.off()## quartz_off_screen
## 2##############
##############
mp=pmd
name_mp<-"PMD"
##adjust the dim of results accordingly
results<-data.frame(row.names = c(1:(length(unique(mp$scen))*length(unique(mp$gamma))*length(unique(mp$k1))*length(unique(mp$k2)))))
results<-mutate(results, OM =NA, MP = NA, Kobe = NA, Catch = NA, LRP_Virgin = NA)
i<-0
for (om in unique(mp$scen))
{
for (x in unique(mp$gamma))
{
for (y in unique(mp$k1))
{
for (z in unique(mp$k2))
{
i<-i+1
MP<-mp %>% filter (gamma == x) %>% filter (k1 == y) %>% filter (k2 == z) %>% filter (scen == om) %>% select (year,iter,ssb,msy_ssb,catch, msy_yield, fbar, msy_harvest, virgin_ssb) %>% mutate(rel_ssb= ssb/msy_ssb) %>% mutate (rel_f= fbar/msy_harvest)%>% mutate(Below_LRP_MSY=(ssb<0.5*msy_ssb)) %>% mutate(Kobe_Green=(rel_ssb > 1 & rel_f<1)) %>% mutate(Target_Catch= catch >0.8*msy_yield) %>% mutate(Below_LRP_Virgin=(ssb<0.2*virgin_ssb)) %>% filter (year %in% c(2019:2038))
results$OM[i]<-paste("OM",om)
results$MP[i]<-paste("D: G=",x," K1=", y, " K2=", z, sep="")
#Probability of being in the green Kobe quadrant
results$Kobe[i]<-sum(MP$Kobe_Green)*100/length(MP$iter)
#Probability of being above 80% of MSY for catch
results$Catch[i]<-sum(MP$Target_Catch)*100/length(MP$iter)
#Probability of being above LRP
results$LRP_Virgin[i]<-sum(1-MP$Below_LRP_Virgin)*100/length(MP$iter)
#CV of catch
results$Catch_Var[i]<-sd(MP$catch)*100/mean(MP$catch)
}
}
}}
results<-mutate(results, Catch_Var = 100-Catch_Var)
dat=melt(results,id=c("MP","OM"))
## rename the performance measures
dat$variable=factor(dat$variable, levels=c("Kobe","Catch","LRP_Virgin","Catch_Var"),
labels=c("Kobe Green","Catch","Safety","Stability"))
## save for later
dat2=dat
##############
## make a table
file_name<-paste("Table_",name_mp,".png",sep ="")
png(file=file_name, width = 800, height = 1000, units = "px")
p2<-ggplot(aes(variable, value),data=dat)+geom_col(aes(fill = value))+
facet_grid(OM~MP,scale="free_y")+scale_fill_gradient2(low ="darkred", high ="yellow", mid = "red", midpoint = 50, limit = c(0,100),name = "Performance")+
ylab("Acceptability")+xlab("Performance Measure")+
ggtitle(paste("Management Procedure:",name_mp, sep ="")) +
labs(subtitle = "MP with 2 options for F target and B trigger")+
theme_light()+theme(axis.text.x=element_text(angle=45, hjust=1),plot.title = element_text(lineheight=.8, face="bold", colour="darkblue", size=12))
p2## Warning: Removed 160 rows containing missing values (position_stack).dev.off()## quartz_off_screen
## 2##############
##############
mp=pmp
name_mp<-"PMP"
##adjust the dim of results accordingly
results<-data.frame(row.names = c(1:(length(unique(mp$scen))*length(unique(mp$k1))*length(unique(mp$k2)))))
results<-mutate(results, OM =NA, MP = NA, Kobe = NA, Catch = NA, LRP_Virgin = NA)
i<-0
for (om in unique(mp$scen))
{
for (y in unique(mp$k1))
{
for (z in unique(mp$k2))
{
i<-i+1
MP<-mp %>% filter (k1 == y) %>% filter (k2 == z) %>% filter (scen == om) %>% select (year,iter,ssb,msy_ssb,catch, msy_yield, fbar, msy_harvest, virgin_ssb) %>% mutate(rel_ssb= ssb/msy_ssb) %>% mutate (rel_f= fbar/msy_harvest)%>% mutate(Below_LRP_MSY=(ssb<0.5*msy_ssb)) %>% mutate(Kobe_Green=(rel_ssb > 1 & rel_f<1)) %>% mutate(Target_Catch= catch >0.8*msy_yield) %>% mutate(Below_LRP_Virgin=(ssb<0.2*virgin_ssb)) %>% filter (year %in% c(2019:2038))
results$OM[i]<-paste("OM",om)
results$MP[i]<-paste("P: K1=", y, " K2=", z, sep="")
#Probability of being in the green Kobe quadrant
results$Kobe[i]<-sum(MP$Kobe_Green)*100/length(MP$iter)
#Probability of being above 80% of MSY for catch
results$Catch[i]<-sum(MP$Target_Catch)*100/length(MP$iter)
#Probability of being above LRP
results$LRP_Virgin[i]<-sum(1-MP$Below_LRP_Virgin)*100/length(MP$iter)
#CV of catch
results$Catch_Var[i]<-sd(MP$catch)*100/mean(MP$catch)
}
}
}
results<-mutate(results, Catch_Var = 100-Catch_Var)
dat=melt(results,id=c("MP","OM"))
## rename the performance measures
dat$variable=factor(dat$variable, levels=c("Kobe","Catch","LRP_Virgin","Catch_Var"),
labels=c("Kobe Green","Catch","Safety","Stability"))
##save for later
dat3=dat
##############
## make a table
file_name<-paste("Table_",name_mp,".png",sep ="")
png(file=file_name, width = 800, height = 1000, units = "px")
p3<-ggplot(aes(variable, value),data=dat)+geom_col(aes(fill = value))+
facet_grid(OM~MP,scale="free_y")+scale_fill_gradient2(low ="darkred", high ="yellow", mid = "red", midpoint = 50, limit = c(0,100),name = "Performance")+
ylab("Acceptability")+xlab("Performance Measure")+
ggtitle(paste("Management Procedure:",name_mp, sep ="")) +
labs(subtitle = "MP with 2 options for F target and B trigger")+
theme_light()+theme(axis.text.x=element_text(angle=45, hjust=1),plot.title = element_text(lineheight=.8, face="bold", colour="darkblue", size=12))
p3
dev.off()## quartz_off_screen
## 2save(dat1,dat2,dat3, file = "Swordfish_MSE_Vis/data/results.RData")Visualise
Create a visual matrix where the rows are operating models and the columns are management procedures:
## make a table and save the image
file_name<-"Swordfish_MSE_Vis/www/Table.png"
png(file=file_name, width = 1200, height = 1000, units = "px")
##get rid of empty columns
dat2NA<-dat2 %>% filter(!is.na(value))
##combine all of the results
dat=rbind(dat1,dat2NA,dat3)
p<-ggplot(aes(variable, value),data=dat)+geom_col(aes(fill = value))+
facet_grid(OM~MP,scale="free_y", labeller=label_wrap_gen(10))+
scale_fill_gradient2(low ="#ea5c0c", high = "#57b88f", mid = "#f4eecd", midpoint = 50, limit = c(0,100),name = "Performance")+
ylab("Acceptability")+xlab("Performance Measure")+
ggtitle(paste("Management Procedures:",sep ="")) +
labs(subtitle = "Four variations of three mps: D = trend based, M = model based, and P = historical period based)")+
theme_light()+theme(axis.text.x=element_text(angle=45, hjust=1),plot.title = element_text(lineheight=.8, face="bold", colour="#204a60", size=14))
p
dev.off()## quartz_off_screen
## 2p
Reference points
Reference points which shape our perception on how well MP perform depend on the OM. Here is one way to illustrate this.
## Reference points
##############
rp2045 =dbGetQuery(conLK,"SELECT * from swonmpb2045")
str(rp2045)## 'data.frame': 4000 obs. of 16 variables:
## $ scen : chr "base" "base" "base" "base" ...
## $ ftar : num 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
## $ btrig : num 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ...
## $ iter : chr "1" "2" "3" "4" ...
## $ catch : num 12689 12716 12367 12671 13168 ...
## $ stock : num 75054 113661 58011 121406 57072 ...
## $ harvest: num 0.169 0.112 0.213 0.104 0.231 ...
## $ msy : num 13306 11972 13047 14820 13218 ...
## $ fmsy : num 0.257 0.146 0.304 0.232 0.325 ...
## $ bmsy : num 51836 82222 42895 63886 40633 ...
## $ r : num 0.1027 0.0582 0.1217 0.0928 0.1301 ...
## $ k : num 238707 378637 197533 294196 187119 ...
## $ p : num -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 ...
## $ b0 : num 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
## $ q1 : num 1.273 0.712 1.791 0.802 1.938 ...
## $ sigma1 : num 0.279 0.343 0.268 0.318 0.306 ...unique(rp2045$scen)## [1] "base" "m/0.1" "m/0.3" "m/lorenzen" "h/low"
## [6] "h/high" "sel/flat" "sel/dome" "wts/0.1" "sigmar/0.8"rpm<-rp2045%>% select (scen,fmsy, bmsy) %>% mutate (scen = factor(scen))
## Reference points, should not depend on the MP
rp<-mp %>% select (scen,msy_ssb, year, msy_harvest, virgin_ssb) %>% filter (year %in% c(2019:2038)) %>% mutate (scen = factor(scen))
str(rp)## 'data.frame': 80000 obs. of 5 variables:
## $ scen : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ msy_ssb : num 29091 29091 29091 29091 29091 ...
## $ year : num 2019 2020 2021 2022 2023 ...
## $ msy_harvest: num 0.124 0.124 0.124 0.124 0.124 ...
## $ virgin_ssb : num 171103 171103 171103 171103 171103 ...file_name<-paste("Swordfish_MSE_Vis/www/RefPoints.png",sep ="")
png(file=file_name, width = 900, height = 900, units = "px")
a<-ggplot(
rp,
aes(x = `msy_ssb`, y = `scen`)
) +
geom_density_ridges_gradient(
aes(fill = ..x..), scale = 3, size = 0.3
) +
scale_fill_gradientn(
colours = c("#816ba1", "#816ba1", "#816ba1"),
name = "SSB MSY"
)+
labs(title = 'OM: SSB MSY') + theme(legend.position="none") + ylab("Operating model") + xlab("t")+
xlim(c(0,90000))
b<-ggplot(
rp,
aes(x = `msy_harvest`, y = `scen`)
) +
geom_density_ridges_gradient(
aes(fill = ..x..), scale = 3, size = 0.3
) +
scale_fill_gradientn(
colours = c("#b6c932", "#b6c932", "#b6c932"),
name = "F MSY"
)+
labs(title = 'OM: F MSY') + theme(legend.position="none") + ylab("Operating model") + xlab("F")+
xlim(c(0,0.5))
c<-ggplot(
rpm,
aes(x = `bmsy`, y = `scen`)
) +
geom_density_ridges_gradient(
aes(fill = ..x..), scale = 3, size = 0.3
) +
scale_fill_gradientn(
colours = c("#c281b1", "#c281b1", "#c281b1"),
name = "Model SSB MSY"
)+
labs(title = 'BDM: SSB MSY') + theme(legend.position="none") + ylab("Operating model") + xlab("t")+
xlim(c(0,90000))
d<-ggplot(
rpm,
aes(x = `fmsy`, y = `scen`)
) +
geom_density_ridges_gradient(
aes(fill = ..x..), scale = 3, size = 0.3
) +
scale_fill_gradientn(
colours = c("#539c7b", "#539c7b", "#539c7b"),
name = "Model F MSY"
)+
labs(title = 'BDM: F MSY') + theme(legend.position="none") + ylab("Operating model") + xlab("t")+
xlim(c(0,0.5))
multiplot(a, b, c, d, cols=2)## Picking joint bandwidth of 4910## Picking joint bandwidth of 0.0198## Picking joint bandwidth of 2110## Picking joint bandwidth of 0.0112dev.off()## quartz_off_screen
## 2multiplot(a, b, c, d, cols=2)## Picking joint bandwidth of 4910## Picking joint bandwidth of 0.0198## Picking joint bandwidth of 2110## Picking joint bandwidth of 0.0112
Visualise changes over time in a dynamic plot
Looking at individual trajectories is a good way to communicate volatility over time. The historic part is the same for all of the simulations, but it is included in the plot to give a sense of how projected volatility compares to annual changes in catches observed in the past.
accumulate_by <- function(dat, var) {
var <- lazyeval::f_eval(var, dat)
lvls <- plotly:::getLevels(var)
dats <- lapply(seq_along(lvls), function(x) {
cbind(dat[var %in% lvls[seq(1, x)], ], frame = lvls[[x]])
})
dplyr::bind_rows(dats)
}
pmb_1<-pmb%>% filter (ftar == 0.5) %>% filter (btrig == 0.6) %>% filter (scen == 1) %>% select (year,iter,ssb,catch) %>% filter (year %in% c(1980:2038)) %>% mutate (Year = factor(year))%>% accumulate_by(~year)
traj<-round(runif(5)*100)
mp<-filter(pmb_1, pmb_1$iter %in% traj)
p <- mp %>%
plot_ly(
x = ~year,
y = ~catch,
split = ~iter,
frame = ~frame,
type = 'scatter',
mode = 'lines',
line = list(simplyfy = F)
) %>%
layout(
xaxis = list(
title = "Date",
zeroline = F
),
yaxis = list(
title = "Catch",
zeroline = F
),
showlegend = FALSE
) %>%
animation_opts(
frame = 100,
transition = 0,
redraw = FALSE
) %>%
animation_slider(
hide = T
) %>%
animation_button(
x = 1, xanchor = "right", y = 0, yanchor = "bottom"
)
p#Save data
save(pmb_1, file = "Swordfish_MSE_Vis/data/mp.RData")