From 267ae109a142b6f3437230070e8c0623a67826e6 Mon Sep 17 00:00:00 2001 From: Iago Mosqueira Date: Thu, 14 Apr 2016 16:51:31 +0200 Subject: [PATCH] Version 0.9.0 --- DESCRIPTION | 20 +- Read-and-delete-me | 9 - vignettes/FLife.Rmd | 866 ++++++++++++++++++++++++++++++++++++++++++++ vignettes/ref.bib | 28 ++ 4 files changed, 906 insertions(+), 17 deletions(-) delete mode 100644 Read-and-delete-me create mode 100644 vignettes/FLife.Rmd create mode 100644 vignettes/ref.bib diff --git a/DESCRIPTION b/DESCRIPTION index 37a23de..88c3b72 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,16 @@ Package: FLife Type: Package -Title: What the package does (short line) -Version: 1.0 -Date: 2015-09-09 +Title: Methods for modelling life history traits +Version: 0.9.0 +Date: 2016-04-13 Author: Laurence Kell -Maintainer: Laurence Kell -Description: More about what it does (maybe more than one line) +Maintainer: Laurence Kell +Description: Many studies have shown the relationships between life history traits for processes such as growth, maturity and natural mortality. Life history has been used to develop priors in stock assesments and to parameterise ecological models. Package has a variety of methods for modelling life history traits and processes. +LazyData: true License: GPL (>= 2) -Depends: methods -Imports: Rcpp (>= 0.12.0) -LinkingTo: Rcpp +Depends: + R(>= 3.0.1), + methods +Imports: + FLCore(>= 2.5), + FLBRP diff --git a/Read-and-delete-me b/Read-and-delete-me deleted file mode 100644 index 6d6236d..0000000 --- a/Read-and-delete-me +++ /dev/null @@ -1,9 +0,0 @@ -* Edit the help file skeletons in 'man', possibly combining help files for multiple - functions. -* Edit the exports in 'NAMESPACE', and add necessary imports. -* Put any C/C++/Fortran code in 'src'. -* If you have compiled code, add a useDynLib() directive to 'NAMESPACE'. -* Run R CMD build to build the package tarball. -* Run R CMD check to check the package tarball. - -Read "Writing R Extensions" for more information. diff --git a/vignettes/FLife.Rmd b/vignettes/FLife.Rmd new file mode 100644 index 0000000..fe1c03c --- /dev/null +++ b/vignettes/FLife.Rmd @@ -0,0 +1,866 @@ +--- +title: "Life History Relationship" +author: "Laurence Kell" +date: "August 13th, 2014" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{FLife} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +bibliography: ref.bib +--- + +```{r knitr_init, echo=FALSE, results="hide"} +library(knitr) +## Global options +opts_chunk$set(echo =!TRUE, + eval =TRUE, + cache =!FALSE, + cache.path="vignette/", + prompt =FALSE, + comment =NA, + message =FALSE, + tidy =FALSE, + warnings=FALSE, + fig.height=4.5, + fig.width =4.5, + fig.path ="tex/") +``` + +# Introduction + +Many studies have shown relationships between life history traits such as growth, maturity and natural mortality. This knowledge has been used to provide advice for data poor stocks, develop priors or fix values for data rich stock assesments and to parameterise ecological models. `FLife` package brings together a variety of methods for modelling life history traits and ecological processes and can be used to create `FLR` objects such as `FLBRP` and `FLStock` in order to model species, population or stock dynamics. + +```{r,echo=FALSE} +library(FLCore) +library(FLash) +library(FLBRP) +library(FLife) +library(popbio) +``` + +# Life history relationships + +## `lhPar` + +In data poor situtuations only the maximum size ($l_{max}$) may be known. Life history relationships can be used to derive the missing parameters [@gislason2008coexistence]. For example $k$ of the [@vonbert1957quantitative] growth equation + +\begin{equation} k=3.15l_{\infty}^{-0.64} \end{equation} + +and the length at which 50\% of the population mature + +\begin{equation} l_{50}=0.72l_{\infty}^{0.93} \end{equation} + +The `lhPar` method takes as its first argument an `FLPar` object with as a minimum a value for `linf` and uses these relationships to derive parameters such as `k` and $l50$ + +```{r lhPar, echo=TRUE} +par=lhPar(FLPar(linf=100)) +``` + +Natural mortality can be estimated from length + +\begin{equation} M=0.55(l-1.66l_{\infty}^{1.44}k\end{equation} + +or mass-at-age [@lorenzen2002density] + +\begin{equation} M=m_1*w^m_2\end{equation} + +where $m_1= 0.55(l_{\infty}^{1.44})k$ and $m_2=-1.61$ + +There are defaults for other values which can not be derived from life history theory. These include $a$ and $b$ from the length weight relationship $w=al^b$, `ato95` the age at which 95\% of fish are mature, offset to age at which 50\% are mature, `sl` selectivity-at-age parameter, standard deviation of lefthand limb of double normal, `sr` stock recruitment relationship, `s` steepness of stock recruitment relationship, `v` virgin biomass. + +Biological processes as growth, maturity and natural mortality can be modelled using functions such as `vonB`, `sigmoid`, and `lorezen` methods. These take as arguments an object for age, length or weight and an `FLPar` with the life history parameters. + +```{r age, echo=TRUE} +age=FLQuant(0:20,dimnames=list(age=0:20)) + +ln =vonB(age,par) +mat=sigmoid(age,par) +wt =par["a"]*ln^par["b"] +wt =len2wt(ln,par) +``` + + +```{r biol, fig.margin=TRUE, fig.cap="Biological age vectors",echo=FALSE} +ggplot(FLQuants(Length=ln,Mass=wt,Maturity=mat))+ + geom_line(aes(age,data))+ + facet_grid(qname~.,scale="free")+ + xlab("Age")+ylab("")+theme_bw() +``` + + +```{r age-m} +m =FLQuants(Lorenzen=lorenzen(wt,par[c("m1","m2")]), + Gislason=gislason(ln,par[-(8:9)])) +``` + +```{r, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +ggplot(m)+ + geom_line(aes(age,data,col=qname))+ + xlab("Age")+ylab("Natural Mortality")+theme_bw()+ + theme(legend.title=element_blank(),legend.position="bottom") +``` + +Selection pattern can be modelled as flat topped or dome shaped by using the double normal function + +```{r age-sel, fig.margin=TRUE, fig.height=4, fig.cap="Selection pattern",echo=FALSE} +sel=FLQuants("Flat"=dnormal(age,FLPar(a1="4",sl=3,sr=500)), + "Dome"=dnormal(age,FLPar(a1="4",sl=3,sr=5))) +``` + +```{r, fig.margin=TRUE, fig.height=4, fig.cap="Selection patterns",echo=FALSE} +ggplot(sel)+ + geom_line(aes(age,data,col=qname))+ + xlab("Age")+ylab("")+theme_bw() +``` + +## Creation of objects + +##FLBRP +```{r} +eql=lhEql(par) +``` + +```{r, fig.margin=TRUE, fig.height=4, fig.cap="Equilibrium",echo=FALSE} +plot(eql)+theme_bw()+theme(legend.position="bottom") +``` + +##FLStock +```{r om} +library(FLash) +om=fwd(eql) +``` + +```{r, fig.margin=TRUE, fig.height=4, fig.cap="Data",echo=FALSE,eval=FALSE} +library("GGally") +library(ggbiplot) + +data(teleost) + +ggpairs(teleost[c("l50","linf","k","b")], + lower=list(params=c(colour="blue")),#,continuous="smooth",se=FALSE), + diag=list(continuous="bar",params=c(colour="blue")), + upper=list(params=list(corSize=6)), axisLabels='none')+ + theme(legend.position = "none", + panel.grid.major = element_blank(), + axis.ticks = element_blank(), + panel.border = element_rect(linetype = 1, colour="black", fill=NA))+ + theme_bw(14) +``` + +```{r lhPar2, fig.margin=TRUE, fig.height=4, fig.cap="",echo=FALSE,eval=FALSE} +teleost=transform(teleost,l50linf=l50/linf) + +family=transform(teleost,what=ifelse(family%in%c("Clupeidae","Scombridae"),ac(family),"other")) + +pc=princomp(teleost[,c("linf","k","b","l50linf")], + cor=TRUE,use="pairwise.complete.obs") +gg=ggbiplot(pc, obs.scale=1, var.scale=1, + ellipse=TRUE, ellipse.prob=.5, circle=FALSE, + groups=factor(family[,"what"]))+ + kobe:::theme_ms(10,legend.position="bottom") + +#gg$layers[[2]]$mapping$colour=NULL +#gg$layers[[3]]$mapping$colour=NULL +gg+theme(legend.position="none") #guides(col=guide_legend(nrow = 3)) +``` + +```{r,eval=!FALSE, fig.margin=TRUE, fig.height=4, fig.cap="",echo=FALSE,eval=FALSE} +par=lhPar(teleost[,c("linf","k","t0","l50","a","b")]) + +attributes(par)$species=teleost[,"species"] +attributes(par)$family =teleost[,"family"] + +mGislason=function(length,params) + 0.55*(length^-1.66)%*%(params["linf"]^1.44)%*%params["k"] +ref=FLife:::lhRef(par,m=mGislason) +``` + +```{r lhPop, fig.margin=TRUE, fig.height=4, fig.cap="",echo=FALSE,eval=FALSE} +dat=ref[dimnames(ref[!is.na(ref[,"rc"]),])[[1]],] + +pc=princomp(dat[,c("r","rc","lopt","sk")],cor=TRUE,use="pairwise.complete.obs") +gg=ggbiplot(pc, obs.scale=1, var.scale=1, + ellipse=TRUE, ellipse.prob=.5, circle=FALSE, + groups=factor(family[,"what"]) )+ + kobe:::theme_ms(12,legend.position="bottom") + +#gg$layers[[2]]=NULL +#gg$layers[[2]]$mapping$colour=NULL +#gg$layers[[3]]$mapping$colour=NULL +gg+theme(legend.position="none") +``` + +# Stability + +An important factor determining a population’s response to perturbation is stability which can be measured in a variety of ways (e.g. Pimm 1984). In its simplest form, a population can be considered stable if it returns to equilibrium after a perturbation. Other definitions expand on this and involve the time taken to return to equilibrium after a perturbation, known as the characteristic return time or population resilience. The lower the characteristic return time, or higher the resilience, the more stable the population. The stability of a population is strongly influenced by the life history of the population and also the pattern of density dependence. For some population models the stability is a good indicator of a population’s response to noise (Taylor 1992), but generally stability is insufficient on its own to predict the response (Horwood 1993). Here we use it to indicate how quickly management can cause an effect in a population, e.g. to recover a stock to a level that would support MSY. In this way, stability can be used as a guide to how controllable the stock is. +For discrete, structured populations this can be calculated using the magnitude of the dominant eigenvalue of the Jacobian matrix evaluated at the equilibrium point (Beddington 1974; Caswell 2001). If the magnitude of this is less than 1 the population will return to equilibrium after a disturbance, with the stability decreasing as the magnitude approaches 1. When the magnitude of the dominant eigenvalue is 1 there is a bifurcation and past this point non-equilibrium dynamics, including extinction, are seen. + + +```{r,init} +library(ggplot2) +library(ggbiplot) +library("GGally") +library(plyr) +library(reshape) +library(FLCore) +library(FLBRP) +library(FLash) +library(FLife) +library(ggplotFL) +library(kobe) +library(popbio) + +data(teleost) +teleost=transform(teleost,l50linf=l50/linf) + +alb=FLPar(unlist(teleost[teleost$species=="Thunnus alalunga", + c("linf","k","t0","l50","a","b")])) +alb=lhPar(rbind(alb,FLPar(m1=0.15,m2=-0.288,s=0.75))) +``` + + +```{r fig1, fig.margin=TRUE, fig.height=4, fig.cap="Overfish",echo=FALSE} +ggpairs(log(100*teleost[,c("b","l50","linf","k","l50linf")]), + lower=list(continuous="smooth", alb=c(colour="blue",se=FALSE)), + diag=list(continuous="bar",alb=c(colour="blue")), + upper=list(alb=list(corSize=6)), axisLabels='none')+ + #scale_x_log10()+scale_y_log10() + theme(legend.position = "none", + panel.grid.major = element_blank(), + axis.ticks = element_blank(), + panel.border = element_rect(linetype = 1, colour="black", fill=NA)) + + theme_bw() +``` + + +```{r fig2, fig.margin=TRUE, fig.height=4, fig.cap="Rebuild",echo=FALSE} +teleost=transform(teleost,l50linf=l50/linf) + +family=transform(teleost, + what=ifelse(family%in%c("Scombridae"), + ac(family),"other")) + +pc=princomp(teleost[,c("linf","k","b","l50linf")], + cor=TRUE,use="pairwise.complete.obs") +gg=ggbiplot(pc, obs.scale=1, var.scale=1, + ellipse=TRUE, ellipse.prob=.5, circle=FALSE, + groups=factor(family[,"what"]))+ + kobe:::theme_ms(10,legend.position="bottom") +gg=gg+geom_point(aes(xvar,yvar),data=gg$data[130,],size=3) + +#gg$layers[[2]]$mapping$colour=NULL +#gg$layers[[3]]$mapping$colour=NULL +gg+theme(legend.position="none") + + scale_colour_manual(values=c("red","blue"))#guides(col=guide_legend(nrow = 3)) +``` + + +```{r fig3, fig.margin=TRUE, fig.height=4, fig.cap="SRR",echo=FALSE} +source('~/Desktop/flr/git/FLife/R/lhEql.R') +refpts=FLBRP:::refpts + +srr=FLBRPs("Beverton and Holt" =lhEql(alb,sr="bevholt"), + "Ricker" =lhEql(alb,sr="ricker"), + "Cushing" =lhEql(alb,sr="cushing"), + "Shepherd" =lhEql(rbind(alb,FLPar(c=1.5)),sr="shepherd"), + "Segmented \nRegression"=lhEql(alb,sr="segreg")) + +srr= + ldply(srr,function(x) { + refpts(x)=refpts(x)["msy"] + fbar(x)=seq(0,1,length.out=501) + res=brp(x) + subset(model.frame(FLQuants(res,"ssb","rec","catch"),drop=TRUE),ssb>=0)}) + +ggplot(melt(srr[,-5],id=c("year","ssb",".id")))+ + geom_vline(aes(xintercept=200))+ + geom_line(aes(ssb,value,col=.id))+ + theme_bw()+theme(legend.position="bottom")+ + scale_colour_manual("Stock Recruit \n Relationship", + values=c("red","green","yellow","blue","pink"))+ + xlab("Spawning Stock Biomass")+ylab("Recruits") +``` + + +```{r fig4, fig.margin=TRUE, fig.height=4, fig.cap="Production Functions",echo=FALSE} +ggplot(melt(srr[,-4],id=c("year","ssb",".id")))+ + geom_path(aes(ssb,value,col=.id))+ + theme_bw()+theme(legend.position="bottom")+ + scale_colour_manual("Stock Recruit \n Relationship", + values=c("red","green","yellow","blue","pink"))+ + xlab("Spawning Stock Biomass")+ylab("Yield") +``` + + + +```{r fig5} +par=lhPar(teleost[,c("linf","k","t0","l50","a","b")]) +attributes(par)$species=teleost[,"species"] +attributes(par)$family =teleost[,"family"] +family=transform(teleost, + what=ifelse(family%in%c("Scombridae"), + ac(family),"other")) + +scen=expand.grid(juve =!c(TRUE,FALSE), + dome =!c(TRUE,FALSE), + steep=c(0.75,.9), + m =c("Gislason","Constant")) + +species=attributes(par)$species + +mGislason=function(length,params) + 0.55*(length^-1.66)%*%(params["linf"]^1.44)%*%params["k"] +mConstant=function(length,params){ + length[]=params["l50"] + 0.55*(length^-1.66)%*%(params["linf"]^1.44)%*%params["k"]} + +ref=mdply(scen,function(dome,juve,steep,m){ + par["s"] =steep + par["a1",]=ifelse(juve,.75,1)*par["a1",] + par["sr",]=ifelse(dome,5,5000) + + if (m=="Gislason"){ + res=lhRef(par,m=mGislason)} + else{ + res=lhRef(par,m=mConstant)} + + data.frame(species,res)}) +``` + + +```{r, fig.margin=TRUE, fig.height=4, fig.cap="Population parameters",echo=FALSE,eval=FALSE} +dat=ref[dimnames(ref[!is.na(ref[,"rc"]),])[[1]],] + +pc=princomp(dat[,c("r","rc","lopt","sk","lro")], + cor=TRUE,use="pairwise.complete.obs") + +gg=ggbiplot(pc, obs.scale=1, var.scale=1, + ellipse=TRUE, ellipse.prob=.5, circle=FALSE, + groups=factor(rep(family[,"what"],16)))+ + kobe:::theme_ms(12,legend.position="bottom") + +gg=gg+geom_point(aes(xvar,yvar),data=gg$data[130,],size=3) + +gg$layers[[2]]=NULL +#gg$layers[[2]]$mapping$colour=NULL +#gg$layers[[3]]$mapping$colour=NULL +gg+theme(legend.position="none")+ + scale_colour_manual(values=c("red","blue")) +``` + + +```{r} +#alb["t0"]=-alb["t0"] +## Beverton and Holt recruitment +bh =lhEql(alb,m=lorenzen) +refpts(bh)=refpts(bh)["msy",] +source('~/Desktop/flr/git/FLBRP/R/plot-ggplot.R') +p=plot(bh)+theme_bw()+theme(legend.position="bottom")+ + scale_colour_manual("",values="red",label="MSY") +``` + +```{r DD,eval=FALSE} +library(FLife) + +source("/home/laurie/Desktop/flr/git/FLife/R/lorenzen.R") + +alb=FLPar(unlist(teleost[teleost$species=="Thunnus alalunga", + c("linf","k","t0","l50","a","b")])) +alb=lhPar(rbind(alb,FLPar(m1=0.15,m2=-0.288,s=0.75))) + +## Beverton and Holt recruitment +bh =lhEql(alb,m=lorenzen) + +fbar(bh)=FLQuant(c(seq(0,4,length.out=101)))* + refpts(bh)["msy","harvest"] +names(dimnames(fbar(bh)))[1]="age" +stk=FLStocks("SRR"=as(bh,"FLStock")) +stk[["SRR"]]=fwd(stk[["SRR"]],f=fbar(bh)[,-1],sr=bh) + +## Cushing recruitment +#sr=fmle(as.FLSR(stk[[1]][,1:10],model="geomean")) +cu =lhEql(alb,m=lorenzen,sr="cushing") +fbar(cu) =fbar(bh) +names(dimnames(fbar(cu)))[1]="age" +stk[["M"]] =as(cu,"FLStock") +stk[["M"]] =fwd(stk[["M"]],f=fbar(cu)[,-1],sr=cu) +stk[["Fecundity"]]=stk[["M"]] + +eql=FLBRPs("Beverton Holt"=bh,"Cushing"=cu) +rm(bh,cu) + +y =eql[["Cushing"]] +ref=stock.n(y)[,1] +n =stock.n(y)[,1] +m =m(y) + +mDD=mdply(data.frame(i=seq(dims(fbar(eql[["Cushing"]]))$year)), + function(i){ + + repeat{ + scale=(stock.n(y)[,i]%-%ref)%/%ref + m(y)[] =mdd(stock.wt(y),alb,scale,k=.9) + + if(sum((m-m(y))^2)<1e-6){break} + m=m(y)} + + subset(model.frame(FLQuants(y,"catch","ssb"),drop=TRUE), + year==i)}) + +p=ggplot(mDD)+geom_line(aes(ssb,catch)) + +y =eql[["Cushing"]] +ref=stock.n(y)[,1] +n =stock.n(y)[,1] +mat=mat(y)#[,max(i-1,dims(m(y))$year)] + +matDD=mdply(data.frame(i=seq(dims(fbar(eql[["Cushing"]]))$year)[-1]), + function(i){ + + repeat{ + scale=(stock.n(y)[,i]%-%ref)%/%ref + mat(y)[] =matdd(ages(scale),alb,scale,k=.9) + + if(sum((mat-mat(y))^2)<1e-6){break} + mat=mat(y)} + + subset(model.frame(FLQuants(y,"catch","ssb"),drop=TRUE), + year==i)}) +``` + + +```{r fig6, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE,eval=FALSE} +source("/home/laurie/Desktop/flr/git/FLife/R/m-dd.R") +data=rbind.fill( + cbind("Form"="M",mDD),cbind("Form"="Fecundity",matDD), + cbind("Form"="Stock-recruit",model.frame(FLQuants(eql[["Cushing"]],"catch","ssb"),drop=TRUE))) + +ggplot(data)+ + geom_line(aes(ssb,catch,colour=Form))+theme_bw()+theme(legend.position="bottom")+ + scale_colour_manual("",values=c("red","green","blue"))+ + xlab("Biomass")+ylab("Yield") +``` + +```{r fig7, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +source('~/Desktop/flr/git/FLife/R/m-dd.R') +source('~/Desktop/flr/git/FLife/R/mat-dd.R') + +k =1 +sr ="cushing" +s =0.7 +alb =FLPar(unlist(teleost[teleost$species=="Thunnus alalunga", + c("linf","k","t0","l50","a","b")])) +alb =lhPar(rbind(alb,FLPar(m1=0.3,m2=-0.288,s=s))) +eq =lhEql(alb,m=lorenzen,sr=sr) +fbar(eq)=FLQuant(c(seq(0,5,length.out=51)* + refpts(eq)["msy","harvest"])) +f =fbar(eq)[,-1] + +stkr=as(eq,"FLStock") +stkr=fwd(stkr,f=f,sr=eq) +ref =stock.n(eq)[,1] + +stkm=stkr +m( stkm)=mdd(stock.wt(stkm),alb,(stock.n(stkm)%-%ref)%/%ref,k) +stkm=fwd(stkm,f=f,sr=eq) + +stkf=stkr +mat(stkf)=matdd(ages(stock.wt(stkf)),alb, + (stock.n(stkf)%-%ref)%/%ref,k,TRUE) +stkf=fwd(stkf,f=f,sr=eq) + +for (i in seq(dims(fbar(eq))$year)[2:50]){ + scale =(stock.n(stkm)[,i]%-%ref)%/%ref + m(stkm)[,i]=mdd(stock.wt(stkm)[,i],alb,scale,k) + + stkm=fwd(stkm,f=f[,i],sr=eq)} + +for (i in seq(dims(fbar(eq))$year)[2:50]){ + scale =(stock.n(stkf)[,i]%-%ref)%/%ref + mat(stkf)[,i]=matdd(ages(scale),alb,scale,k,TRUE) + + stkf=fwd(stkf,f=f[,i],sr=eq)} + +plot(FLStocks("SRR"=stkr[,-(1:1)], + "M" =stkm[,-(1:1)], + "Fecundity"=stkf[,-(1:1)]))+ + theme_bw()+ + theme(legend.position="bottom") + +p=ggplot(as.data.frame(FLQuants("SRR"=ssb(stkr), + "M"=ssb(stkm), + "Fecundity"=ssb(stkf))))+ + geom_line(aes(year,data,col=qname))+ + theme_bw()+ + theme(legend.position="bottom") +``` + + +```{r fig8, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +k =1 +sr ="cushing" +s =0.7 +alb =FLPar(unlist(teleost[teleost$species=="Thunnus alalunga", + c("linf","k","t0","l50","a","b")])) +alb =lhPar(rbind(alb,FLPar(m1=0.3,m2=-0.288,s=s))) +eq =lhEql(alb,m=lorenzen,sr=sr) +fbar(eq)=FLQuant(c(seq(5,.5,length.out=51)* + refpts(eq)["msy","harvest"])) +f =fbar(eq)[,-1] + +stkr=as(eq,"FLStock") +stkr=fwd(stkr,f=f,sr=eq) +ref =stock.n(eq)[,1] + +stkm=stkr +m( stkm)=mdd(stock.wt(stkm),alb,(stock.n(stkm)%-%ref)%/%ref,k) +stkm=fwd(stkm,f=f,sr=eq) + +stkf=stkr +mat(stkf)=matdd(ages(stock.wt(stkf)),alb, + (stock.n(stkf)%-%ref)%/%ref,k,TRUE) +stkf=fwd(stkf,f=f,sr=eq) + +for (i in seq(dims(fbar(eq))$year)[2:50]){ + scale =(stock.n(stkm)[,i]%-%ref)%/%ref + m(stkm)[,i]=mdd(stock.wt(stkm)[,i],alb,scale,k) + + stkm=fwd(stkm,f=f[,i],sr=eq)} + +for (i in seq(dims(fbar(eq))$year)[2:50]){ + scale =(stock.n(stkf)[,i]%-%ref)%/%ref + mat(stkf)[,i]=matdd(ages(scale),alb,scale,k,TRUE) + + stkf=fwd(stkf,f=f[,i],sr=eq)} + +plot(FLStocks("SRR"=stkr[,-(1:1)], + "M" =stkm[,-(1:1)], + "Fecundity"=stkf[,-(1:1)]))+ + theme_bw()+ + theme(legend.position="bottom") + +p=ggplot(as.data.frame(FLQuants("SRR"=ssb(stkr), + "M"=ssb(stkm), + "Fecundity"=ssb(stkf))))+ + geom_line(aes(year,data,col=qname),size=2)+ + theme(legend.position="bottom") +``` + +```{r fig9, fig.margin=TRUE, fig.height=4, fig.cap="Year effects",echo=FALSE} +data(ple4) +res=noise(1,m(ple4)[1:8,ac(1980:2008)],burn=10,b=0.9,what="age") +ggplot()+ + geom_point(aes(year,age,size= data), + data=subset(as.data.frame(res),data>0))+ + geom_point(aes(year,age,size=-data), + data=subset(as.data.frame(res),data<=0),colour="red")+ + scale_size_area(max_size=4, guide="none")+ + theme_bw() +``` + +```{r fig10, fig.margin=TRUE, fig.height=4, fig.cap="Cohort Effects",echo=FALSE} +source('~/Desktop/flr/git/FLife/R/noise.R') +res=noise(1,m(ple4)[1:8,ac(1980:2008)],burn=10,b=0.9,what="cohort") +ggplot()+ + geom_point(aes(year,age,size= data), + data=subset(as.data.frame(res), + data>0))+ + geom_point(aes(year,age,size=-data), + data=subset(as.data.frame(res), + data<=0),colour="red")+ + scale_size_area(max_size=4, guide="none")+ + theme_bw() +``` + +```{r fig11, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +k =1 +sr ="cushing" +s =0.7 +alb =FLPar(unlist(teleost[teleost$species=="Thunnus alalunga", + c("linf","k","t0","l50","a","b")])) +alb =lhPar(rbind(alb,FLPar(m1=0.3,m2=-0.288,s=s))) +eq =lhEql(alb,m=lorenzen,sr=sr) +fbar(eq)=FLQuant(rep(1,1001)*refpts(eq)["msy","harvest"]) +f =fbar(eq)[,-1] +f =propagate(fbar(eq)[,-1],3) +f[,,,,,2]=0 +f[,,,,,3]=f[,,,,,1]*4 + +stkr=as(eq,"FLStock") +stkr=fwd(stkr,f=f,sr=eq) +ref =stock.n(eq)[,1] + +scale=noise(1,stock.wt(stkr),sd=0.3,b=0.9,what="cohort") +stkm=stkr + +m( stkm)=mdd(iter(stock.wt(stkm),1),alb,scale,k) +stkm=fwd(stkm,f=f,sr=eq) + +scale=noise(1,stock.wt(stkm),sd=0.3,b=0.9,what="cohort") +stkf=stkr +mat(stkf)=matdd(ages(stock.wt(stkf)),alb,scale,k,TRUE) +stkf=fwd(stkf,f=f,sr=eq) +stkr=fwd(stkr,f=f,sr=eq,sr.residuals=rlnorm(1,iter(f,1)*0,.3)) +dat =as.data.frame(FLQuants("SRR"=ssb(stkr), + "M"=ssb(stkm), + "Fecundity"=ssb(stkf))) +fishMat=ddply(subset(dat,year>50),.(iter,qname), transform, val=data/mean(data)) + +k =1 +sr ="cushing" +s =0.7 +alb =FLPar(unlist(teleost[teleost$species=="Thunnus alalunga", + c("linf","k","t0","l50","a","b")])) +alb =lhPar(rbind(alb,FLPar(m1=0.3,m2=-0.288,s=s))) +alb[c("a1","sr")]=c(0,5000) +eq =lhEql(alb,m=lorenzen,sr=sr) +fbar(eq)=FLQuant(rep(1,1001)*refpts(eq)["msy","harvest"]) +f =fbar(eq)[,-1] +f =propagate(fbar(eq)[,-1],3) +f[,,,,,2]=0 +f[,,,,,3]=f[,,,,,1]*4 + +stkr=as(eq,"FLStock") +stkr=fwd(stkr,f=f,sr=eq) +ref =stock.n(eq)[,1] + +scale=noise(1,stock.wt(stkr),sd=0.3,b=0.9,what="cohort") + +stkm=stkr +m( stkm)=mdd(iter(stock.wt(stkm),1),alb,scale,k) +stkm=fwd(stkm,f=f,sr=eq) +scale=noise(1,stock.wt(stkm),sd=0.3,b=0.9,what="cohort") +stkf=stkr +mat(stkf)=matdd(ages(stock.wt(stkf)),alb,scale,k,TRUE) +stkf=fwd(stkf,f=f,sr=eq) +stkr=fwd(stkr,f=f,sr=eq,sr.residuals=rlnorm(1,iter(f,1)*0,.3)) +dat =as.data.frame(FLQuants("SRR"=ssb(stkr), + "M"=ssb(stkm), + "Fecundity"=ssb(stkf))) +fishJuv=ddply(subset(dat,year>50),.(iter,qname), transform, val=data/mean(data)) + +dat=rbind(cbind("Selection"="Juvenile",fishJuv), + cbind("Selection"="Mature", fishMat)) +dat=transform(dat,F=factor(paste("F times",c(0,1,3))[iter])) + +ggplot(dat)+ + geom_line(aes(year,val,col=qname))+ + theme_bw()+ + theme(legend.position="bottom")+ + facet_grid(F~Selection,scale="free_y")+ + scale_x_continuous(limits=c(500,700))+ + xlab("Time (year)")+ylab("") +``` + + +```{r fig12, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +dat=rbind(cbind("Selection"="Juvenile",fishJuv), + cbind("Selection"="Mature", fishMat)) +dat=transform(dat,F=factor(paste("F times",c(0,1,3))[iter])) + +dat=ddply(dat,.(Selection,qname,F), with, +# as.data.frame(spectrum(data, spans = c(7,13), log = "dB", ci = 0.8,plot=FALSE)[c("freq","spec")])) + as.data.frame(spectrum(data, log = "dB", ci = 0.8,plot=FALSE)[c("freq","spec")])) + +dat=ddply(subset(dat,freq>0.05),.(F,Selection,qname),transform,val=spec/max(spec)) +ggplot(dat,aes(freq,val,col=qname))+ + geom_smooth(se=FALSE)+ + facet_grid(F~Selection,scale="free_y")+ + theme_bw()+ + theme(legend.position="bottom") +``` + + +```{r fig13, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality Mis-specification",echo=FALSE} +library(FLAssess) + +stk =iter(stkm[,101:250],1) +m(stk)=iter(m(stkr),1)[,101:250] +vpa =stk+VPA(stk) + +p=plot(FLStocks("Actual"=stk,"Model Estimate"=vpa)) +p$data=subset(p$data,qname%in%c("SSB","Rec")) +names(p$data)[8]="data" +p$data=ddply(p$data,.(qname), transform, val=data/mean(data)) +ggplot(subset(p$data,year>50))+ + geom_line(aes(year,val,col=qname))+ + facet_grid(qname~stock,scale="free")+ + theme_bw()+ + theme(legend.position="none")+ + xlab("Year")+ylab("") +``` + + +```{r fig14, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +library(FLife) +library(FLash) +library(FLBRP) + +par=lhPar(FLPar(linf=100)) +eql=lhEql(par) +mou=as(eql,"FLStock") +mou=FLash:::fwd(mou,f=fbar(eql)[,-1]/4,sr=eql) + +srDev=rlnorm(100,FLQuant(rep(0,136)),0.3) + +om=FLife:::mseSBT1(mou,eql,srDev, + start=dims(mou)$maxyear,end=dims(mou)$maxyear+20,lag=1,interval=3, + k1=1.5,k2=3.0,gamma=1,nyrs=5, + seed=7890,nits=100, + uCV =0.2) + +ggplot()+ + geom_line(aes(year,data,col=iter), + data=as.data.frame(FLQuants(iter(om[,ac(95:120)],5:7),c("Rec"=rec,"SSB"=ssb, + "Catch"=catch,"Harvest"=fbar)),drop=T))+ + facet_grid(qname~.,scale="free")+ + theme_bw()+xlab("")+ylab("")+ + theme(legend.position="bottom") +``` + + +```{r fig15, fig.margin=TRUE, fig.height=4, fig.cap="Natural Mortality",echo=FALSE} +cas=read.csv("/home/laurie/MEGAsync/mse/trade-offs/inputs/casSWOM8513_v1.csv", + colClasses=c("factor","numeric","factor","factor","factor","factor","numeric","numeric","numeric")) +names(cas)=c("species","year","flag","fleet","gear","stock","len","ln5","n") +powh=function(len,n){ + + fn=function(len,n){ + require(plyr) + + res=ddply(data.frame(n=n,len=len), .(len), function(x) sum(x$n)) + res=res[order(res$len),] + + csum =rev(cumsum(rev(res$V1))) + clsum=rev(cumsum(rev(res$len*res$V1))) + mn =clsum/csum + + data.frame(mn=mn,diff=mn-res$len,len=res$len,n=res$V1)} + + linf=function(x) -coefficients(x)[1]/coefficients(x)[2] + zk =function(x) (-1-coefficients(x)[2])/coefficients(x)[2] + + dat=fn(len,n) + + res=lm(diff~len,data=dat) + + params=c("linf"=linf(res),"zk"=zk( res)) + names(params)=c("linf","zk") + + return(list(params=params, + data =dat))} + +pw=ddply(subset(cas), .(year), + function(cas) powh(cas$len,cas$n)$data) + +pw=transform(pw, lustrum=(year%/%5)*5, + yr =year-(year%/%5)*5, + weight=ifelse(len>=100&len<=200,1,0)) + +ggplot(pw)+ + geom_line(aes(len,diff,colour=factor(yr),group=year))+ + scale_x_continuous(limits=c(0,300)) + + facet_wrap(~lustrum,ncol=2)+ + xlab("Length (cm)")+ + ylab("Difference between Length and Mean Size")+ + geom_smooth(aes(len,diff,weight=weight), + method="lm",col="red",size=1.25,alpha=.1)+ + theme_bw()+theme(legend.position="none") +``` + + +#Elasticity + +A measure of proportional effect, i.e., the effect that a change in a given matrix element has as a proportional to the change in that element + +```{r elasticity, fig.margin=TRUE, fig.height=4, fig.cap="",echo=FALSE,eval=FALSE} +library(ggplot2) +library(FLCore) +library(FLBRP) +library(FLife) +library(numDeriv) +library(popbio) + +source('~/Desktop/flr/git/FLife/R/lhEql.R') + loptAge=function(params, + m =function(length,params) params["m1"]%*%(exp(log(length)%*%params["m2"])), + growth=vonB, + ...){ + + loptFn=function(x,params,m){ + + age =0:ceiling(x) + dmns =list(age=age) + length=vonB(age=FLQuant(pmin(age+0.5,x),dimnames=dmns),params=params) + m. =FLQuant(m(length,params), dimnames=dmns) + mCum =FLQuant(aaply(m.,6,sum)) + n =exp(-mCum) + c(n*len2wt(length[ac(ceiling(x))],params))} + + dmns=dimnames(params) + dmns$params="lopt" + dm =dim(params) + + res=aaply(params,seq(length(dm))[-1],function(x){ + x.=FLPar(x) + rtn=try(optimise(loptFn,c(0,40),params=x.,maximum=TRUE,m=m)$maximum) + if ("character" %in% mode(rtn)) rtn=NA + rtn}) + + vonB(FLQuant(FLPar(array(res,dim=c(1, dm[-1]),dimnames=dmns))),params) + + } + +par=lhPar(FLPar(linf=100)) + +fn=function(par,all=FALSE){ + eql=lhEql(par) + + ref=refpts(eql)[,c(1:2,4:5)] + h =melt(ref[,2]/ref[,4])[,-3] + h$quantity="harvest" + + dimnames(ref)[[2]][1]="f" + ref=melt(ref[drop=T]) + + lop=loptAge(par) + + f =refpts(eql)["crash","harvest"] + r =log(lambda(leslie(eql,fbar=FLQuant(c(f)))[drop=T])) + + f =refpts(eql)["msy","harvest"] + rc =log(lambda(leslie(eql,fbar=FLQuant(c(f)))[drop=T])) + + res=rbind(ref,h, + cbind("refpt"="r", "quantity"="pop rate", "value"=melt(r[drop=T])), + cbind("refpt"="rc", "quantity"="pop rate", "value"=melt(rc[drop=T])), + cbind("refpt"="lopt","quantity"="length", "value"=melt(lop[drop=T]))) + + res=rbind(subset(res,refpt%in%c("r","rc","lopt")), + subset(res,refpt%in%c("msy")&quantity=="ssb")) + + res[,3]=log(res[,3]) + + if(all) res else c(res[,3])} + +val=fn(par,TRUE) +jcb=jacobian(fn,par) + +dimnames(jcb)=list(y=val$ref,x=dimnames(par)$params) +jcb=melt(jcb) +jcb$quantity=val$quantity + +#ggplot(subset(jcb,quantity!="pop rate"))+ +# geom_point(aes(x,y,size=log(abs(value)))) + +#ggplot(subset(jcb,quantity=="pop rate"))+ +# geom_point(aes(x,y,size=log(abs(value)))) + +ggplot(ddply(jcb,.(y),transform, value=value/max(abs(value))))+ + geom_point(aes(x,y,size=abs(value)))+ + scale_size_area()+ + theme_bw()+ + theme(legend.position="none") + +``` diff --git a/vignettes/ref.bib b/vignettes/ref.bib new file mode 100644 index 0000000..86a9605 --- /dev/null +++ b/vignettes/ref.bib @@ -0,0 +1,28 @@ +@article{gislason2008coexistence, + title = {Coexistence in North Sea fish communities: implications for growth and natural mortality}, + volume = {65}, + number = {4}, + journal = {{ICES} Journal of Marine Science: Journal du Conseil}, + author = {Gislason, H. and Pope, {J.G.} and Rice, {J.C.} and Daan, N.}, + year = {2008}, + pages = {514-530} +} + +@article{lorenzen2002density, + title={Density-dependent growth as a key mechanism in the regulation of fish populations: evidence from among-population comparisons}, + author={Lorenzen, Kai and Enberg, Katja}, + journal={Proceedings of the Royal Society of London. Series B: Biological Sciences}, + volume={269}, + number={1486}, + pages={49--54}, + year={2002}, + publisher={The Royal Society} +} + +@article{vonbert1957quantitative, + title = {Quantitative laws in metabolism and growth}, + journal = {Quarterly Review of Biology}, + author = {Von Bertalanffy, L.}, + year = {1957}, + pages = {217-231} +} \ No newline at end of file