Home

The MyDas project simulation tested a range of assessment models and methods in order to establish Maximum Sustainable Yield (MSY), or proxy MSY reference points across the spectrum of data-limited stocks. A main output of the project was the development of two R packages mydas and FLife, both are implemented as FLR packages with example vignettes are available.

MyDas was funded by the Irish exchequer and EMFF 2014-2020

Contact: Laurence Kell and Coilin Minto

---

Introduction

The mydas package is part of FLR family of R packages, a collection of tools for quantitative fisheries science, developed in the R language, that facilitates the construction of bio-economic simulation models of fisheries systems.

These can be installed from the website

install.packages(c("FLCore","FLFishery","FLasher","FLBRP","FLife","mydas"), 
             repos="http://flr-project.org/R")

A range of tutorials are also available.

The mydas package can be installed from this repository using devtools if you want the latest development version, e.g.

install.packages("devtools",dependencies=TRUE)

library(devtools)

Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="TRUE")

devtools::install_github("flr/mydas")

Libraries

library(plyr)
library(reshape)
library(ggplot2)

FLR libraries

library(FLCore)
library(FLasher)
library(FLBRP)
library(FLife)
library(mydas)

Plotting

library(ggplotFL)

Plotting is done using ggplot2 which provides a powerful alternative paradigm for creating both simple and complex plots in R using the ideas the Grammar of Graphics. The idea of the grammar is to specify the individual building blocks of a plot and then to combine them to create the desired graphic.

ggplotFL implements a number of plots for the main FLR objects.

The ggplot functions expects a data.frame for its first argument, data; then a geometric object geom that specifies the actual marks put on to a plot and an aesthetic that is "something you can see" have to be provided. Examples of geometic Objects (geom) include points (geom_point, for scatter plots, dot plots, etc), lines (geom_line, for time series, trend lines, etc) and boxplot (geom_boxplot, for, well, boxplots!). Aesthetic mappings are set with the aes() function and, examples include, position (i.e., on the x and y axes), color ("outside" color), fill ("inside" color), shape (of points), linetype and size.

Operating Model Conditioning

Life history parameters

Life history parameters can be loaded from the fishnets github repository, a library of multivariate priors for fish population dynamics parameters

load(url("https://github.com//fishnets//fishnets//blob//master//data//fishbase-web//fishbase-web.RData?raw=True"))

Select turbot as a case study

lh=subset(fb,species=="Psetta maxima")

names(lh)[c(14,17)] = c("l50","a50")
lh=lh[,c("linf","k","t0","a","b","a50","l50")]

head(lh)

Get the means and create an FLPar object

lh=apply(lh,2,mean,na.rm=T)
lh=FLPar(lh)

Values can be replace with any better estimates available.

Then lhPar fills in missing values using life history theory, while quantities like selection-at-age are set using defaults, in this case set to be the same as maturity-at-age.

?lhPar

par=lhPar(lh)
par

The parameters can then be used by lhEql to simulate equilibrium dynamics by combining the spawner/yield per recruit relationships with a stock recruitment relationship, by creating an FLBRP object

eq=lhEql(par)

The FLBRP object created has vectors for mass, maturity, natural mortality and selectivity-at-age

sel<-function(x) 
  catch.sel(x)%/%fapex(catch.sel(x))

dat=FLQuants(eq,"Mass"=stock.wt,"M"=m,"Maturity"=mat,"Selectivity"=sel)

res=ldply(dat,function(x) cast(as.data.frame(quantile(x,probs=c(0.025,0.25,0.5,0.75,0.975))),
                               age~iter,value="data"))

ggplot(res)+
  geom_ribbon(aes(age,ymin=`25%`,ymax=`75%`),alpha=0.5,fill="red")+
  geom_ribbon(aes(age,ymin=`2.5%`,ymax=`97.5%`),alpha=0.1,fill="red")+
  geom_line(aes(age,`50%`))+
  facet_wrap(~.id,scale="free")+
  scale_x_continuous(limits=c(0,10))+
  xlab("Age")+ylab("")

plot(eq,refpts="msy")

To model time series the FLBRP object created by lhEql is coerced to an FLStock object and then projected forward

For example a fishing time series is simulated that represents a stock that was originally lightly exploited, then effort is increased until the stock is overfished and then fishing pressure was reduced to recover the stock biomass to $B_{MSY}$.

fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(.1,19),
                                              seq(.1,2,length.out = 30),
                                              seq(2,1.0,length.out = 10)[-1],
                                              rep(1.0,61)))[,1:105]
plot(fbar(eq))

Coercing the FLBRP object into an FLStock

om=as(eq,"FLStock")

We can simply plot the forward F projection catch, recruitment and biomass estimates etc.. We call fwd() with the stock, the control object (fbar) and the stock recruitment relationship, and look at the results

om=fwd(om,fbar=fbar(om)[,-1], sr=eq)

can add units

units(catch(om)) = units(discards(om)) = units(landings(om)) = units(stock(om)) = 'tonnes'
units(catch.n(om)) = units(discards.n(om)) = units(landings.n(om)) = units(stock.n(om)) = '1000'
units(catch.wt(om)) = units(discards.wt(om)) = units(landings.wt(om)) = units(stock.wt(om)) = 'kg'

plot(om)

Stochasticity

To add recruitment variability into add the residuals argument, e.g. for 100 Monte Carlo simulations with a CV of 0.3

nits=100
srDev=rlnoise(nits, rec(om) %=% 0, 0.3)

om=propagate(om,nits)

om=fwd(om,fbar=fbar(om)[,-1], sr=eq, residuals=srDev)

plot(om)