In laurieKell/biodyn: Biomass Dynamic Model

library(knitr)

source("/home/laurie/Desktop/flr/book/paths.R")

## Global options
options(max.print="75")
opts_chunk$set(fig.path=pathFig,
               echo=FALSE,
               cache=FALSE,
               cache.path="cache/",
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE)
opts_knit$set(width=75)

install.packages("FLCore",repos="http://flr-project.org/R")
install.packages("FLash", repos="http://flr-project.org/R")
install.packages("FLBRP", repos="http://flr-project.org/R")
install.packages("kobe")
install.packages("diags", repos="http://flr-project.org/R")

library(ggplot2)
library(kobe)

library(biodyn)
library(aspic)
library(diags)
library(plyr)

theme_set(theme_bw(base_size=16))

swon=aspic("/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/swon/2013/aspic/run2/aspic.inp")

swon=fit(swon)

plot(swon)

Bootstraps

swonBoot=boot(swon,boot=50)
plot(swonBoot)+
  theme(legend.position="bottom")

library(doParallel)
library(foreach)

cl=makeCluster(4)
registerDoParallel(cl)

files=c("http://rscloud.iccat.int/kobe/albn/2013/aspic/run1/aspic.inp",
        "http://rscloud.iccat.int/kobe/albn/2013/aspic/run2/aspic.inp",
        "http://rscloud.iccat.int/kobe/albn/2013/aspic/run3/aspic.inp",
        "http://rscloud.iccat.int/kobe/albn/2013/aspic/run4/aspic.inp")
files=c("/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/albn/2013/aspic/run1/aspic.inp",
        "/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/albn/2013/aspic/run2/aspic.inp",
        "/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/albn/2013/aspic/run3/aspic.inp",
        "/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/albn/2013/aspic/run4/aspic.inp")

#### Bug fix, include missing methods
setMethod('boot',  signature(object='aspics'),
          function(object, dir=tempdir(), package=class(object), exeNm="aspic",boot=500,
                   .combine=NULL,
                   .multicombine=T,.maxcombine=10,.packages=c("aspic","plyr")){

            if (is.null(.combine)) ..combine=list else ..combine=.combine
            res=foreach(i=names(object), .combine=..combine,
                        .multicombine=.multicombine,
                        .maxcombine  =.maxcombine,
                        .packages    =.packages) %dopar% {  
                          wkdir=tempfile('file', dir)
                          dir.create(wkdir, showWarnings = FALSE)
                          boot(object[[i]],dir=wkdir,boot=boot)}

            if (is.null(.combine)) {
              res=aspics(res)
              names(res)=names(object)}

            res})

# plotComp.R - 
# ggplotFL/R/plotComp.R

# Copyright 2003-2007 FLR Team. Distributed under the GPL 2 or later
# Maintainer: Iago Mosqueira, JRC, Laurie Kell, ICCAT
# $Id:  $

# whooow {{{
whooow  =function(x,fn,probs)
  as.data.frame(FLQuants(lapply(fn,
                                function(fn,x)
                                  quantile(fn(x), probs=probs, na.rm=T), x=x))) # }}}

# plotComp {{{
plotComp = function(x, fn=NULL, probs=c(0.75,0.50,0.25), size=c(0.5,1.0,0.5),
                    lty=c(2,1,2), facet=facet_wrap(~qname, scale="free"),worm=NA) {

  if (dims(x)$iter>=length(probs)){  
    res = whooow(x,fn,probs)
    p1  = ggplot(res) + geom_path(aes(x=year,y=data,group=iter,size=iter,lty=iter)) +
      scale_size_manual(    values=size, name="Quantile") +
      scale_linetype_manual(values=lty , name="Quantile")
  }else{
    res = whooow(x,fn,0.5)
    p1  = ggplot(res) + geom_path(aes(x=year,y=data))    }

  p1 = p1 +  expand_limits(y = 0) +
    xlab("Year") + ylab("") +
    facet

  if (length(worm) > 0)
    if (length(worm)<=dims(x)$iter)
      if (!(length(worm)==1 & is.na(worm[1])))  
        p1=p1+geom_path(aes(year,data,group=iter,colour=iter),
                        data=transform(subset(as.data.frame(FLQuants(lapply(fn,function(f,x) f(x), x=x))),iter %in% worm),iter=factor(iter)))

  p1
} # }}}

# plotComps {{{
plotComps = function(x, fn=NULL, probs=c(0.75,0.50,0.25), size=c(0.5,1.0,0.5),
                     lty=c(2,1,2), facet=facet_wrap(~qname,scale="free")) {

  x=x[seq(length(x))]

  if (max(laply(x,function(x) dims(x)$iter))>=length(probs))
    res = ldply(x, whooow, fn=fn, probs=probs)
  else
    res = ldply(x, whooow, fn=fn, probs=0.5)

  if ("X1" %in% names(res)) 
    names(res)[names(res)=="X1"]=".id"

  if (".id" %in% names(res)) 
    res$.id  = factor(res$.id)

  res$iter = factor(res$iter)

  if (length(unique(res$iter))>=length(probs)){
    p1  = ggplot(res) + geom_path(aes(x=year,y=data,group=.id:iter,size=iter,col=.id)) +
      scale_size_manual(    values=size, name="Quantile") +
      scale_linetype_manual(values=lty , name="Quantile")
  }else{
    p1  = ggplot(res) + geom_path(aes(x=year,y=data,group=.id,col=.id)) 
  }  

  p1= p1 + expand_limits(y = 0) +
    xlab("Year") + ylab("") +
    facet

  p1} 
# }}}

iter=iterators:::iter
albnBoots=boot(aspics(files),boot=10)
rm(iter)

plot(albnBoots)

Biodyn

Coercion of one class into another

bd  =as(swon,"biodyn")
cpue=index(swon,FALSE)

bd@control[c("b0","p"),"phase"]=-1

bd  =fit(bd,cpue) 

ASPIC v biodyn

plot(aspics(aspic=as(swon,"biodyn"),biodyn=bd))+
  theme(legend.position="bottom")

Comparison of uncertainty

Create a list class with the individual biodyn objects in

sims=biodyns("Deterministic"=bd)

Jackknife

sims[["Jack Knife"]]=fit(bd,jackknife(cpue))

plot(sims[["Jack Knife"]])

##############################################################
#' plotJack
#' 
#' Create a \code{ggplot2} plot based on a jack knifed biodyn and plots 
#' time series of biomass and harvest rate. 
#' The basic object can then be modified by adding ggpot2 layers.
#'
#' @param  \code{x}, an object of class \code{biodyn} that has been jack knifed, i.e. by 
#' providing a jack knifed CPUE series to fit
#' @param \code{y} the original \code{biodyn} object
#'
#' @return an \code{ggplot2} object
#' 
#' @export
#' @docType methods
#' @rdname plotJack
#'
#' @examples
#' \dontrun{
#' bd=fit(bd,cpue)
#' jk=fit(bd,jacknife(cpue))
#' plotJack(jk,bd)
#' }  
plotJack=function(x,y,ncol=1){

  js=function(x,y, ...) {

    n <- dims(x)$iter 

    mnU <- apply(x, 1:5, mean)   

    SS <- apply(sweep(x, 1:5, mnU,"-")^2, 1:5, sum)

    bias <- (n - 1) * (mnU - y)
    se <- sqrt(((n-1)/n)*SS)

    res=FLQuants(list(jack.mean=y, jack.se=se, jack.bias=bias))

    attributes(res)$jackknife=TRUE

    return(res)}

  df=rbind(cbind(qname="stock",  
                 model.frame(js(stock(  x),stock(  y)),drop=TRUE)),
           cbind(qname="harvest",
                 model.frame(js(harvest(x),harvest(y)),drop=TRUE)))

  # basic plot data vs. year
  p=ggplot(data=df, aes(x=year, y=jack.mean))+
    facet_wrap(~qname,ncol=ncol,scale="free_y")+
    geom_ribbon(aes(x=year, ymin=jack.mean-2*jack.se, 
                    ymax=jack.mean+2*jack.se),
                fill="blue", alpha = .20)+
    # line + xlab + ylab + limits to include 0 +
    geom_line(colour="red") + xlab("Year") + ylab("") + expand_limits(y=0) +
    geom_line(aes(year,jack.mean+jack.bias),colour="black")   +
    # no legend
    theme(legend.title = element_blank())

  return(p)}

plotJack(sims[["Jack Knife"]],sims[["Deterministic"]])

MCMC

Conduct Monte Carlo Markov Chain simulations, using the algorithm in ADMB.

sims[["MCMC"]]=fit(sims[[1]],cpue,cmdOps=c("-mcmc 1000000, -mcsave 80000"))

Diagnostics

acf(c(params(sims[["MCMC"]])["r"]))

variance/covariance matrix

Use the covariance matrix to conduct a Monte Carlo simulation of the parameter estimates to derive uncertainty in derived quantitries.

sims[["Vcov"]]=mvn(sims[["Deterministic"]],500,nms=c("r","k"),fwd=TRUE)

Can use the covariance matrix to estimate uncertainty in other parameters, e.g. in management quantities

sims[["Deterministic"]]@mng

sims[["Deterministic"]]@mngVcov

variance/covariance matrix

currentState   =sims[["Deterministic"]]@mng[c("bbmsy","ffmsy"),"hat",drop=T]
currentStateVar=sims[["Deterministic"]]@mngVcov[c("bbmsy","ffmsy"),
                                                c("bbmsy","ffmsy"),drop=T]

mvrnorm(10,currentState,currentStateVar)

ggplot(data=as.data.frame(mvrnorm(100,currentState,currentStateVar)),
       aes(bbmsy,ffmsy))+geom_point()+
  geom_hline(aes(yintercept=1))+
  geom_vline(aes(xintercept=1))+
  xlab(expression(B:B[MSY]))+
  ylab(expression(F:F[MSY]))

Covariance matrix

sims[["Vcov"]]=mvn(bd,100,fwd=TRUE)

Bootstrap CPUE and use this to simulate a new CPUE series

CPUE series

cpueSim=sims[[1]]@diags[,c("year","hat")]
names(cpueSim)[2]="data"
cpueSim=as.FLQuant(cpueSim)

cpueSim=rlnorm(100,log(cpueSim),0.25)
plot(cpueSim,na.rm=TRUE)

cpueSim[cpueSim==0]=NA

sims[["CPUE"]]=fit(propagate(sims[[1]],100),cpueSim)

library(kobe)
 setGeneric('kobe',          function(object,method,...)    standardGeneric('kobe'))

source('~/Desktop/flr/pkgs/biodyn/R/biodyn-kobe.R')

Marginal Density for Stock/BMSY

df=kobe(sims[-(1:2)])
ggplot(subset(df,year==2011)) + 
  geom_density(aes(x=stock, y=..count..), position = "stack",fill="red") +
  geom_vline(aes(xintercept=1))          +
  facet_wrap(~.id,scale="free_y",ncol=1)+
  scale_x_continuous(limits=c(0.5,1.5))

Marginal Density for Harvest/FMSY

ggplot(subset(df,year==2011)) + 
  geom_density(aes(x=harvest, y=..count..), position = "stack",fill="red") +
  geom_vline(aes(xintercept=1))          +
  facet_wrap(~.id,scale="free_y",ncol=1)+
  scale_x_continuous(limits=c(0.5,1.5))

Kobe Phase Plot

library(kobe)
kobePhase()+
  geom_point(aes(stock,harvest),data=subset(df,year==2011))+
  facet_wrap(~.id,ncol=1)

Projections

Use simulated data

bd   =biodyn:::sim()

## simulate HCRs, annual, tri-annula, F bound, TAC bound
bd=window(bd,end=29)
for (i in seq(29,49,1))
  bd=fwd(bd,harvest=hcr(bd,refYrs=i,yrs=i+1)$hvt)
sims=biodyns("1"=bd)

plot(sims)

bd=window(bd,end=29)
for (i in seq(29,49,3))
  bd=fwd(bd,harvest=hcr(bd,refYrs=i,yrs=i+1:3)$hvt)
sims[["3"]]=bd

plot(sims)

bd=window(bd,end=29)
for (i in seq(29,49,1))
  bd=fwd(bd,harvest=hcr(bd,refYrs=i,yrs=i+1,bndF=c(0.9,1.1))$hvt)
sims[["bound F"]]=bd

plot(sims)

bd=window(bd,end=30)
for (i in seq(29,49,1))
  bd=fwd(bd,catch=hcr(bd,refYrs=i,yrs=i+1,tac=T,bndTac=c(0.9,1.1))$tac)
sims[["bound TAC"]]=bd

plot(sims)

p.=plot(sims[c(1,3)])+
  theme(legend.position="bottom")+
  scale_colour_manual(values=c("cyan4","grey10"),
                      labels=c("HCR","10% Constraint on Effort"))+
  guides(col=guide_legend(title=NULL))+
  scale_x_continuous(limits=c(27,50),breaks=c(30,40,50),labels=c(1,11,21))+
  scale_y_continuous(breaks=NULL)

p.$data=transform(p.$data,qname=factor(qname,levels=c("Harvest","Yield","Stock")))
p.+geom_hline(aes(yintercept=X1),data=
                cbind("qname"=c("Yield","Harvest","Stock"),data.frame(refpts(bd))),col="red")

kb=ldply(sims[c(1)],function(x,sd=.1)
  model.frame(mcf(FLQuants(stock  =rlnorm(100,log(stock(  x)%/%bmsy(x)),sd),
                           harvest=rlnorm(100,log(harvest(x)%/%fmsy(x)),sd)))))
kb=subset(kb,year%in% 29:50)

pt=ldply(sims[c(1)],function(x)
  model.frame(mcf(FLQuants(stock  =stock(  x)%/%bmsy(x),
                           harvest=harvest(x)%/%fmsy(x)))))
pt =subset(pt,year%in% 1:50)
pt.=ddply(pt,.(year,.id),with,data.frame(stock=median(stock),harvest=median(harvest)))

i=40
print(kobePhase(subset(kb,year%in%i))+
        geom_line(aes(stock,harvest),data=biodyn:::hcrPlot(bd),col="brown",size=1.5)+
        ggtitle(paste("Year",i-29)) +
        theme(legend.position="none",plot.background=element_rect(fill="transparent",colour=NA),
              plot.title = element_text(lineheight=.8, face="italic"))+
        xlim(0,2)+ylim(0,2)+
        geom_point(aes(stock,harvest,col=.id,fill=.id), shape=21, size = 3) + 
        geom_point(aes(stock,harvest,group=.id,fill=.id),col="black", data=subset(pt.,year==i),shape=21,size=5)+
        geom_path( aes(stock,harvest,col=.id,group=.id) ,data=subset(pt, year<=i),size=1)+
        scale_fill_manual(  values=c("cyan1","green","red","yellow")) + 
        scale_colour_manual(values=c("cyan4","green","red","yellow"))+
        coord_cartesian(xlim=c(0,2),ylim=c(0,2)))

kb=ldply(sims[c(1,3)],function(x,sd=.1)
  model.frame(mcf(FLQuants(stock  =rlnorm(100,log(stock(  x)%/%bmsy(x)),sd),
                           harvest=rlnorm(100,log(harvest(x)%/%fmsy(x)),sd)))))
kb=subset(kb,year%in% 29:50)

pt=ldply(sims[c(1,3)],function(x)
  model.frame(mcf(FLQuants(stock  =stock(  x)%/%bmsy(x),
                           harvest=harvest(x)%/%fmsy(x)))))
pt =subset(pt,year%in% 1:50)
pt.=ddply(pt,.(year,.id),with,data.frame(stock=median(stock),harvest=median(harvest)))

i=40
kobePhase(subset(kb,year%in%i))+
  # geom_line(aes(stock,harvest),data=hcrPlot(bd),col="brown",size=1.5)+
  ggtitle(paste("Year",i-29)) +
  theme(legend.position="none",plot.background=element_rect(fill="transparent",colour=NA),
        plot.title = element_text(lineheight=.8, face="italic"))+
  xlim(0,2)+ylim(0,2)+
  geom_point(aes(stock,harvest,col=.id,fill=.id,size=.id), shape=21) + 
  geom_point(aes(stock,harvest,fill=.id,group=.id),data=subset(pt.,year==i),shape=21,col="black",size=5)+
  geom_path( aes(stock,harvest,col=.id,group=.id) ,data=subset(pt, year<=i),size=1)+
  scale_fill_manual(  values=c("cyan1","grey90","green","red","yellow")) + 
  scale_size_manual(  values=c(1,3)) + 
  scale_colour_manual(values=c("cyan4","grey20","green","red","yellow"))+
  coord_cartesian(xlim=c(0,2),ylim=c(0,2))

pe=rlnorm(100,FLQuant(0,dimnames=list(year=1:50)),0.5)
sims[[1]]=fwd(sims[[1]],harvest=harvest(sims[[1]])[,ac(1:50)],pe=pe)
sims[[3]]=fwd(sims[[3]],harvest=harvest(sims[[3]])[,ac(1:50)],pe=pe)
p.=plot(sims[c(1,3)])+
  theme(legend.position="bottom")+
  scale_colour_manual(values=c("cyan4","grey10"),
                      labels=c("HCR","10% Constraint on Effort"))+
  guides(col=guide_legend(title=NULL))+
  scale_x_continuous(limits=c(27,50),breaks=c(30,40,50),labels=c(1,11,21))+
  scale_y_continuous(breaks=NULL)

p.$data=transform(p.$data,qname=factor(qname,levels=c("Harvest","Yield","Stock")))
p.+geom_hline(aes(yintercept=X1),data=
                cbind("qname"=c("Yield","Harvest","Stock"),data.frame(refpts(bd))),col="red")

laurieKell/biodyn documentation built on May 20, 2019, 7:58 p.m.