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
results<-data.frame(row.names = c(1:(length(unique(mpb$scen))*length(unique(mpb$btrig))*length(unique(mpb$ftar)))))
results<-mutate(results, OM =NA, MP = NA, Kobe = NA, Catch = NA, LRP = NA, Catch_Var = NA)
i<-0
for (om in unique(mpb$scen))
{
for (tar in unique(mpb$ftar))
{
for (trigger in unique(mpb$btrig))
{
i<-i+1
mp<-mpb %>% 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=(ssb<0.2*virgin_ssb)) %>% mutate(Kobe_Green=(rel_ssb > 1 & rel_f<1)) %>% mutate(Target_Catch= catch >0.8*msy_yield) %>% filter (year %in% c(2019:2038))
results$OM[i]<-paste("OM",om)
results$MP[i]<-paste("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[i]<-sum(1-mp$Below_LRP)*100/length(mp$iter)
#CV of catch
results$Catch_Var[i]<-sd(mp$catch)*100/mean(mp$catch)
}
}
}
str(results)## 'data.frame': 20 obs. of 6 variables:
## $ OM : chr "OM 1" "OM 1" "OM 1" "OM 1" ...
## $ MP : chr "F_tar=0.7 B_trig=0.7" "F_tar=0.7 B_trig=0.5" "F_tar=0.5 B_trig=0.7" "F_tar=0.5 B_trig=0.5" ...
## $ Kobe : num 10.8 10.8 64.8 64.8 0.8 ...
## $ Catch : num 99.5 99.5 95.7 95.7 99.5 ...
## $ LRP : num 22.7 22.7 70.8 70.8 96 ...
## $ Catch_Var: num 10.53 10.53 11.06 11.06 8.17 ...head(results)## OM MP Kobe Catch LRP Catch_Var
## 1 OM 1 F_tar=0.7 B_trig=0.7 10.8 99.45 22.70 10.532183
## 2 OM 1 F_tar=0.7 B_trig=0.5 10.8 99.45 22.70 10.532183
## 3 OM 1 F_tar=0.5 B_trig=0.7 64.8 95.65 70.75 11.058794
## 4 OM 1 F_tar=0.5 B_trig=0.5 64.8 95.65 70.75 11.058794
## 5 OM 2 F_tar=0.7 B_trig=0.7 0.8 99.55 96.05 8.171587
## 6 OM 2 F_tar=0.7 B_trig=0.5 0.8 99.55 96.05 8.171587head(mp)## year iter ssb msy_ssb catch msy_yield fbar msy_harvest
## 1 2019 1 37071.11 28118.77 13205.24 12845.98 0.1221912 0.1393908
## 2 2020 1 35965.48 28118.77 13205.24 12845.98 0.1216114 0.1393908
## 3 2021 1 35542.33 28118.77 13205.24 12845.98 0.1173420 0.1393908
## 4 2022 1 35474.85 28118.77 15349.39 12845.98 0.1463844 0.1393908
## 5 2023 1 34345.95 28118.77 15349.39 12845.98 0.1500020 0.1393908
## 6 2024 1 33425.92 28118.77 15349.39 12845.98 0.1516701 0.1393908
## virgin_ssb rel_ssb rel_f Below_LRP Kobe_Green Target_Catch
## 1 120642.1 1.318376 0.8766088 FALSE TRUE TRUE
## 2 120642.1 1.279056 0.8724496 FALSE TRUE TRUE
## 3 120642.1 1.264007 0.8418208 FALSE TRUE TRUE
## 4 120642.1 1.261607 1.0501730 FALSE FALSE TRUE
## 5 120642.1 1.221460 1.0761260 FALSE FALSE TRUE
## 6 120642.1 1.188740 1.0880927 FALSE FALSE TRUE#Save data
save(results, file = "Swordfish_MSE_Vis/data/results.RData")Visualise
Create a visual matrix where the rows are operating models and the columns are management procedures:
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","Catch_Var"),
labels=c("Kobe Green","Catch","Safety","Stability"))
## Example plot
png(file="Swordfish_MSE_Vis/www/Table.png", width = 700, height = 480, units = "px")
p<-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("Management Procedure:") +
labs(subtitle = "MPB 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))
p
dev.off()## quartz_off_screen
## 2p 
Visualise changes in distributions over time
One way to look at dynamics of the stock over time is to visualise the distributions.
mp<-mpb1000 %>% filter (ftar == 0.5) %>% filter (btrig == 0.5) %>% filter (scen == 1) %>% select (year,iter,ssb,msy_ssb,catch, msy_yield, fbar, msy_harvest) %>% mutate (rel_ssb= ssb/msy_ssb) %>% mutate (rel_f= fbar/msy_harvest)%>% mutate(Below_LRP=(ssb<0.5*msy_ssb)) %>% mutate(Kobe_Green=(rel_ssb > 1 & rel_f<1)) %>% mutate(Target_Catch= catch >0.8*msy_yield) %>% filter (year %in% c(1980, 1990, 2000, 2010, 2020, 2030)) %>% mutate (year = factor(year))
png(file="Swordfish_MSE_Vis/www/Distributions.png", width = 700, height = 480, units = "px")
p<-ggplot(
mp,
aes(x = `ssb`, y = `year`)
) +
geom_density_ridges_gradient(
aes(fill = ..x..), scale = 3, size = 0.3
) +
scale_fill_gradientn(
colours = c("#0D0887FF", "#CC4678FF", "#F0F921FF"),
name = "SSB"
)+
labs(title = 'Swordfish SSB, OM 1, MP with F_tar = 0.5 and B_trig = 0.5')
p## Picking joint bandwidth of 9090dev.off()## quartz_off_screen
## 2p ## Picking joint bandwidth of 9090
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
png(file="Swordfish_MSE_Vis/www/RefPoints.png", width = 700, height = 480, units = "px")
rp<-mpb %>% select (scen,msy_ssb, year, msy_harvest, virgin_ssb) %>% filter (year %in% c(2019:2038)) %>% mutate (scen = factor(scen))
str(rp)## 'data.frame': 40000 obs. of 5 variables:
## $ scen : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ msy_ssb : num 34172 34172 34172 34172 34172 ...
## $ year : num 2019 2020 2021 2022 2023 ...
## $ msy_harvest: num 0.125 0.125 0.125 0.125 0.125 ...
## $ virgin_ssb : num 189173 189173 189173 189173 189173 ...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("aquamarine2", "aquamarine3", "aquamarine4"),
name = "SSB MSY"
)+
labs(title = 'SSB MSY') + theme(legend.position="none") + ylab("Operating model") + xlab("t")
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("bisque2", "bisque3", "bisque4"),
name = "F MSY"
)+
labs(title = 'F MSY') + theme(legend.position="none") + ylab("Operating model") + xlab("F")
c<-ggplot(
rp,
aes(x = `virgin_ssb`, y = `scen`)
) +
geom_density_ridges_gradient(
aes(fill = ..x..), scale = 3, size = 0.3
) +
scale_fill_gradientn(
colours = c("cadetblue2", "cadetblue3", "cadetblue4"),
name = "F MSY"
)+
labs(title = 'Unfished SSB') + theme(legend.position="none") + ylab("Operating model") + xlab("t")
multiplot(a, b, c, cols=2)## Picking joint bandwidth of 7580## Picking joint bandwidth of 0.0157## Picking joint bandwidth of 36500dev.off()## quartz_off_screen
## 2multiplot(a, b, c, cols=2)## Picking joint bandwidth of 7580## Picking joint bandwidth of 0.0157## Picking joint bandwidth of 36500
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)
}
mp<-mpb1000 %>% filter (ftar == 0.5) %>% filter (btrig == 0.5) %>% filter (scen == 1) %>% select (year,iter,ssb,msy_ssb,catch, msy_yield, fbar, msy_harvest) %>% mutate (rel_ssb= ssb/msy_ssb) %>% mutate (rel_f= fbar/msy_harvest)%>% mutate(Below_LRP=(ssb<0.5*msy_ssb)) %>% mutate(Kobe_Green=(rel_ssb > 1 & rel_f<1)) %>% mutate(Target_Catch= catch >0.8*msy_yield) %>% filter (year %in% c(1980:2038)) %>% mutate (Year = factor(year))%>% filter (iter %in% c(36, 6, 11, 288, 642))%>%
accumulate_by(~year)
#Save data
save(mp, file = "Swordfish_MSE_Vis/data/mp.RData")
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