In laurieKell/mydas-pkg: Tools for conducting Management Strategy Evaluation for data limited stocks
Introduction
To provide advice on the status of data poor stocks ICES uses $MSY$ proxy reference points as part of a Precautionary Approach.
Data poor stocks include those for which only trends such as lpue, cpue, and mean length in the catch are available (Category 3), and stocks for which only reliable catch data are available (Category 4).
Methods currently approved by ICES for calculation of $MSY$ reference points for these stocks are
- Length based indicators
- Z derived from mean length
- Length based spawner per recruit; and
- Surplus Production models
Many approaches have emerged over the last few decades, for example Where length data are available methods include Length Based Spawning Potential Ratio (LBSPR), Length-Based Integrated Mixed Effects (LIME), and Length-Based Bayesian (LBB). While where only catch data are available methods include Catch-Maximum Sustainable Yield (Catch-MSY), State-Space Catch-Only Model (SSCOM), Depletion Based Stock Reduhction Analysis (DBSRA), and Simple Stock Synthesis (SSS) an extension of Catch-MSY (CMSY).
Empirical indicators and reference points can also be used to monitor stocks and these include
- $L_{max5\%}$ mean length of largest 5\%
- $L_{95\%}$ $95^{th}$ percentile
- $P_{mega}$ Proportion of individuals above $L_{opt} + 10\%$
- $L_{25\%}$ $25^{th}$ percentile of length distribution
- $L_{c}$ Length at $50\%$ of modal abundance
- $L_{mean}$ Mean length of individuals $> L_c$
- $L_{max_{y}}$ Length class with maximum biomass in catch
- $L_{mean}$ Meanlength of individuals $> L$
where potential reference points include
- $L_{opt} = L_{\infty}\frac{3}{3+\frac{M}{K}}$, assuming $M/K = 1.5$ gives $\frac{2}{3}L_{\infty}$
- $L_{F=M} = 0,75l_c+0.25l_{\infty}$
Methods
Simulation
Run scenarios with an increasing trend in F that leads to overfishing, then implement a recovery plan that brings fishing to the $F_{MSY}$ level then screen potential empirical MPs by
- Generating length, catch and CPUE indicators using an OEM
- Compare indicators to OM using Reciever Operating Characteristics (ROCs)
Receiver Operating Characteristics
Sort the observed outcomes by their predicted scores with the highest scores first, then calculate cumulative True Positive Rate (TPR) and True Negative Rate (TNR) for the ordered observed outcomes
knitr::opts_chunk$set(echo = FALSE)
library(knitr)
opts_chunk$set(comment =NA,
warning =FALSE,
message =FALSE,
error =FALSE,
echo =FALSE,
fig.width =6,
fig.height=6,
cache =TRUE,
fig.path ="tex/indicators-",
cache.path="cache/indicators/")
iFig=0
iTab=0
\newpage
library(FLCore)
library(FLBRP)
library(FLasher)
library(FLife)
library(mydas)
library(ggplotFL)
library(plyr)
library(dplyr)
library(reshape)
library(popbio)
library(spatstat)
par=lhPar(FLPar(c(linf= 59.1, k=0.28, t0=-0.4,
a=0.01111,b=3.15,a50=4.0, l50=43.25),units="NA"))
par["m1"]=.03
eq=lhEql(par)
fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(0.1,60),
seq(0.1,2,length.out=40)[-40],
seq(2,1,length.out=11),rep(1,20)))
set.seed(234)
srDev=rlnoise(100,fbar(eq)%=%0,0.3,0)
om=as(eq,"FLStock")
om=window(fwd(propagate(om,dim(srDev)[6]),fbar=fbar(eq)[,-1],sr=eq,residuals=srDev),
start=41)
plot(FLQuants(om,
"f" = function(x) fbar(x)%/%refpts(eq)["msy","harvest"],
"ssb" = function(x) ssb(x)%/%refpts( eq)["msy","ssb"],
"catch"=function(x) landings(x)%/%refpts(eq)["msy","yield"],
"rec" = function(x) rec(x)%/%refpts( eq)["msy","rec"])) +
geom_hline(aes(yintercept=1),col="red",linetype=2)+
theme_bw()
Figure r iFig=iFig+1; iFig Time series relative to MSY benchmarks.
\newpage
Mean length indicators
ind=mydas:::omSmry(om,par)
refs=mydas:::popdyn(par,eq)
rf1=FLQuants(
f =window(fbar(eq)/c(refs["fmsy"]),start=41),
clmsy=as.FLQuant(transmute(ind, year=year,iter=iter,data=cln/c(refs["clmsy"]))),
slmsy=as.FLQuant(transmute(ind, year=year,iter=iter,data=sln/c(refs["slmsy"]))),
clfm =as.FLQuant(transmute(ind, year=year,iter=iter,data=cln/c(refs["lfm"]))),
slfm =as.FLQuant(transmute(ind, year=year,iter=iter,data=sln/c(refs["lfm"]))))
quad=data.frame(x=c(40,80, 80, 40, 80,90,90,80, 90,105,105,90, 105,110,110,105, 110,130,130,110),
y=rep(c(-Inf,-Inf,Inf,Inf),5),
f=c(rep("green",4),rep("amber",4),rep("red",4),rep("amber2",4),rep("green2",4)))
p=plot(rf1)+
geom_polygon(aes(x,y,fill=f),data=quad,alpha=0.2)+
scale_fill_manual(values=c("orange","orange","green","green","red","blue"))+
facet_grid(qname~.,scale="free")+
geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 90,105)),col="red")+
geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 80,110)),col="orange")
p5=p$layers[[5]]
p$layers[[5]]=p$layers[[1]]
p$layers[[1]]=p5
p
Figure r iFig=iFig+1; iFig. Mean length indicators compared to $F:F_{MSY}$, vertical lines indicate 1 (green), 1.5 (orange) and 2 (red) times $F_{MSY}$.
\newpage
Length Frequencies
source('~/Desktop/sea++/mydas/pkg/R/oemLn.R')
## Matrix with P(len) by age for each set of params
agLs=alk(par,cv=0.1)
## Simulate length frequencies by year
lfd=lenSample(catch.n(om),iter(agLs,i),nsample=100)
lfd=subset(as.data.frame(lfd),data>0)
library(spatstat)
rln=ddply(subset(lfd,year%in%c(60,100,120)), .(year),
with,mydas:::lenInd(len,data,lopt=c(refs["lopt"]))[-5])
ggplot(subset(lfd,year%in%c(60,100,120)&iter==1))+
geom_histogram(aes(len,weight=data),binwidth=1)+
geom_vline(aes(xintercept=value,col=variable),data=melt(rln,id="year"))+
facet_grid(year~.)+
xlab("Length (cm)")+ylab("")+
scale_color_manual("Indicator",values=rainbow(4))
save(lfd,file="/home/laurence/tmp/t.RloadData")
Figure r iFig=iFig+1; iFig. Simulated length frequencies distributions with indicators.
save(lfd,refs,eq,par,file="/home/laurence/tmp/t.RData")
source('~/Desktop/sea++/mydas/pkg/R/oemLn.R')
inds=ddply(lfd, .(year,iter), with, lenInd(len,data,lopt=c(refs["lopt"])))
rf2=FLQuants(
f =window(fbar(eq)/c(refs["fmsy"]),start=41),
lbar =as.FLQuant(transmute(inds,year=year,iter=iter,data=lbar))/c(par["l50"]),
lfm =as.FLQuant(transmute(inds,year=year,iter=iter,data=lfm)),
pmega=as.FLQuant(transmute(inds,year=year,iter=iter,data=pmega)),
l95 =as.FLQuant(transmute(inds,year=year,iter=iter,data=l95/c(par["linf"]))),
l25 =as.FLQuant(transmute(inds,year=year,iter=iter,data=l25/c(par["linf"]))))
p=plot(rf2)+
facet_grid(qname~.,scale="free")+
geom_polygon(aes(x,y,fill=f),data=quad,alpha=0.2)+
geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 90,105)),col="red")+
geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 80,110)),col="orange")+
scale_fill_manual(values=c("orange","orange","green","green","red","blue"))
p5=p$layers[[5]]
p$layers[[5]]=p$layers[[1]]
p$layers[[1]]=p5
p
Figure r iFig=iFig+1; iFig. Time series of indicators compared to $F:F_{MSY}$, vertical lines indicate 1 (green), 1.5 (orange) and 2 (red) times $F_{MSY}$.
\newpage
Receiver Operating Characteristics
Detection of overfishing
simple_roc <- function(labels, scores){
labels <- labels[order(scores, decreasing=TRUE)]
data.frame(TPR=cumsum(labels)/sum(labels),
FPR=cumsum(!labels)/sum(!labels),
labels,
scores=sort(scores))}
ind1=ldply(rf1[-1],function(x) as.data.frame(window(x,start=70,end=100),drop=T))
ind1=ddply(ind1,.(.id),transform,overfished=as.numeric(!(year%in%90:105)))
ind1=ddply(ind1,.(.id),with, simple_roc(overfished,data))
tab=ddply(ind1,.(.id), with, {
best=max(TPR*(1-FPR))==TPR*(1-FPR);
data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])})
ggplot(ind1)+
geom_point(aes(FPR,TPR,col=.id))+
geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+
geom_point(aes(FPR,TPR),data=tab,size=2)+
xlab("False Negative Rate")+ylab("True Positive Rate")+
scale_color_manual("Indicator",values=rainbow(4))+
theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the four indicators of overfishing, points indicate the optimum value of the indicator.
tab
ind2=ldply(rf2[-1],function(x) as.data.frame(window(x,start=70,end=100),drop=T))
ind2=ddply(ind2,.(.id),transform,overfished=as.numeric(!(year%in%90:105)))
ind2=ddply(ind2,.(.id),with, simple_roc(overfished,data))
tab=ddply(ind2,.(.id), with, {
best=max(TPR*(1-FPR))==TPR*(1-FPR);
data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])})
ggplot(ind2)+
geom_point(aes(FPR,TPR,col=.id))+
geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+
geom_point(aes(FPR,TPR),data=tab,size=2)+
xlab("False Negative Rate")+ylab("True Positive Rate")+
scale_color_manual("Indicator",values=rainbow(5))+
theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the three indicators of overfishing, points indicate the optimum value of the indicator.
tab
Detection of recovery
ind1=ldply(rf1[-1],function(x) as.data.frame(window(x,start=105),drop=T))
ind1=ddply(ind1,.(.id),transform,overfished=as.numeric((year>110)))
ind1=ddply(ind1,.(.id),with, simple_roc(overfished,data))
tab=ddply(ind1,.(.id), with, {
best=max(TPR*(1-FPR))==TPR*(1-FPR);
data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])})
ggplot(ind1)+
geom_point(aes(FPR,TPR,col=.id))+
geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+
geom_point(aes(FPR,TPR),data=tab,size=2)+
xlab("False Negative Rate")+ylab("True Positive Rate")+
scale_color_manual("Indicator",values=rainbow(4))+
theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the four indicators of recovery, points indicate the optimum value of the indicator.
tab
ind2=ldply(rf2[-1],function(x) as.data.frame(window(x,start=105),drop=T))
ind2=ddply(ind2,.(.id),transform,overfished=as.numeric((year>105)))
ind2=ddply(ind2,.(.id),with, simple_roc(overfished,data))
tab=ddply(ind2,.(.id), with, {
best=max(TPR*(1-FPR))==TPR*(1-FPR);
data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])})
ggplot(ind2)+
geom_point(aes(FPR,TPR,col=.id))+
geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+
geom_point(aes(FPR,TPR),data=tab,size=2)+
xlab("False Negative Rate")+ylab("True Positive Rate")+
scale_color_manual("Indicator",values=rainbow(5))+
theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the three indicators of recovery, points indicate the optimum value of the indicator.
tab
\newpage
References {#References}
\newpage
Session Info
sessionInfo()
laurieKell/mydas-pkg documentation built on Nov. 8, 2019, 10:58 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Introduction
To provide advice on the status of data poor stocks ICES uses $MSY$ proxy reference points as part of a Precautionary Approach.
Data poor stocks include those for which only trends such as lpue, cpue, and mean length in the catch are available (Category 3), and stocks for which only reliable catch data are available (Category 4).
Methods currently approved by ICES for calculation of $MSY$ reference points for these stocks are
- Length based indicators
- Z derived from mean length
- Length based spawner per recruit; and
- Surplus Production models
Many approaches have emerged over the last few decades, for example Where length data are available methods include Length Based Spawning Potential Ratio (LBSPR), Length-Based Integrated Mixed Effects (LIME), and Length-Based Bayesian (LBB). While where only catch data are available methods include Catch-Maximum Sustainable Yield (Catch-MSY), State-Space Catch-Only Model (SSCOM), Depletion Based Stock Reduhction Analysis (DBSRA), and Simple Stock Synthesis (SSS) an extension of Catch-MSY (CMSY).
Empirical indicators and reference points can also be used to monitor stocks and these include
- $L_{max5\%}$ mean length of largest 5\%
- $L_{95\%}$ $95^{th}$ percentile
- $P_{mega}$ Proportion of individuals above $L_{opt} + 10\%$
- $L_{25\%}$ $25^{th}$ percentile of length distribution
- $L_{c}$ Length at $50\%$ of modal abundance
- $L_{mean}$ Mean length of individuals $> L_c$
- $L_{max_{y}}$ Length class with maximum biomass in catch
- $L_{mean}$ Meanlength of individuals $> L$
where potential reference points include
- $L_{opt} = L_{\infty}\frac{3}{3+\frac{M}{K}}$, assuming $M/K = 1.5$ gives $\frac{2}{3}L_{\infty}$
- $L_{F=M} = 0,75l_c+0.25l_{\infty}$
Methods
Simulation
Run scenarios with an increasing trend in F that leads to overfishing, then implement a recovery plan that brings fishing to the $F_{MSY}$ level then screen potential empirical MPs by
- Generating length, catch and CPUE indicators using an OEM
- Compare indicators to OM using Reciever Operating Characteristics (ROCs)
Receiver Operating Characteristics
Sort the observed outcomes by their predicted scores with the highest scores first, then calculate cumulative True Positive Rate (TPR) and True Negative Rate (TNR) for the ordered observed outcomes
knitr::opts_chunk$set(echo = FALSE) library(knitr) opts_chunk$set(comment =NA, warning =FALSE, message =FALSE, error =FALSE, echo =FALSE, fig.width =6, fig.height=6, cache =TRUE, fig.path ="tex/indicators-", cache.path="cache/indicators/") iFig=0 iTab=0
\newpage
library(FLCore) library(FLBRP) library(FLasher) library(FLife) library(mydas) library(ggplotFL) library(plyr) library(dplyr) library(reshape) library(popbio) library(spatstat)
par=lhPar(FLPar(c(linf= 59.1, k=0.28, t0=-0.4, a=0.01111,b=3.15,a50=4.0, l50=43.25),units="NA")) par["m1"]=.03 eq=lhEql(par) fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(0.1,60), seq(0.1,2,length.out=40)[-40], seq(2,1,length.out=11),rep(1,20))) set.seed(234) srDev=rlnoise(100,fbar(eq)%=%0,0.3,0) om=as(eq,"FLStock") om=window(fwd(propagate(om,dim(srDev)[6]),fbar=fbar(eq)[,-1],sr=eq,residuals=srDev), start=41)
plot(FLQuants(om, "f" = function(x) fbar(x)%/%refpts(eq)["msy","harvest"], "ssb" = function(x) ssb(x)%/%refpts( eq)["msy","ssb"], "catch"=function(x) landings(x)%/%refpts(eq)["msy","yield"], "rec" = function(x) rec(x)%/%refpts( eq)["msy","rec"])) + geom_hline(aes(yintercept=1),col="red",linetype=2)+ theme_bw()
Figure r iFig=iFig+1; iFig Time series relative to MSY benchmarks.
\newpage
Mean length indicators
ind=mydas:::omSmry(om,par)
refs=mydas:::popdyn(par,eq)
rf1=FLQuants( f =window(fbar(eq)/c(refs["fmsy"]),start=41), clmsy=as.FLQuant(transmute(ind, year=year,iter=iter,data=cln/c(refs["clmsy"]))), slmsy=as.FLQuant(transmute(ind, year=year,iter=iter,data=sln/c(refs["slmsy"]))), clfm =as.FLQuant(transmute(ind, year=year,iter=iter,data=cln/c(refs["lfm"]))), slfm =as.FLQuant(transmute(ind, year=year,iter=iter,data=sln/c(refs["lfm"]))))
quad=data.frame(x=c(40,80, 80, 40, 80,90,90,80, 90,105,105,90, 105,110,110,105, 110,130,130,110), y=rep(c(-Inf,-Inf,Inf,Inf),5), f=c(rep("green",4),rep("amber",4),rep("red",4),rep("amber2",4),rep("green2",4))) p=plot(rf1)+ geom_polygon(aes(x,y,fill=f),data=quad,alpha=0.2)+ scale_fill_manual(values=c("orange","orange","green","green","red","blue"))+ facet_grid(qname~.,scale="free")+ geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 90,105)),col="red")+ geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 80,110)),col="orange") p5=p$layers[[5]] p$layers[[5]]=p$layers[[1]] p$layers[[1]]=p5 p
Figure r iFig=iFig+1; iFig. Mean length indicators compared to $F:F_{MSY}$, vertical lines indicate 1 (green), 1.5 (orange) and 2 (red) times $F_{MSY}$.
\newpage
Length Frequencies
source('~/Desktop/sea++/mydas/pkg/R/oemLn.R') ## Matrix with P(len) by age for each set of params agLs=alk(par,cv=0.1) ## Simulate length frequencies by year lfd=lenSample(catch.n(om),iter(agLs,i),nsample=100) lfd=subset(as.data.frame(lfd),data>0)
library(spatstat) rln=ddply(subset(lfd,year%in%c(60,100,120)), .(year), with,mydas:::lenInd(len,data,lopt=c(refs["lopt"]))[-5])
ggplot(subset(lfd,year%in%c(60,100,120)&iter==1))+ geom_histogram(aes(len,weight=data),binwidth=1)+ geom_vline(aes(xintercept=value,col=variable),data=melt(rln,id="year"))+ facet_grid(year~.)+ xlab("Length (cm)")+ylab("")+ scale_color_manual("Indicator",values=rainbow(4)) save(lfd,file="/home/laurence/tmp/t.RloadData")
Figure r iFig=iFig+1; iFig. Simulated length frequencies distributions with indicators.
save(lfd,refs,eq,par,file="/home/laurence/tmp/t.RData")
source('~/Desktop/sea++/mydas/pkg/R/oemLn.R') inds=ddply(lfd, .(year,iter), with, lenInd(len,data,lopt=c(refs["lopt"]))) rf2=FLQuants( f =window(fbar(eq)/c(refs["fmsy"]),start=41), lbar =as.FLQuant(transmute(inds,year=year,iter=iter,data=lbar))/c(par["l50"]), lfm =as.FLQuant(transmute(inds,year=year,iter=iter,data=lfm)), pmega=as.FLQuant(transmute(inds,year=year,iter=iter,data=pmega)), l95 =as.FLQuant(transmute(inds,year=year,iter=iter,data=l95/c(par["linf"]))), l25 =as.FLQuant(transmute(inds,year=year,iter=iter,data=l25/c(par["linf"]))))
p=plot(rf2)+ facet_grid(qname~.,scale="free")+ geom_polygon(aes(x,y,fill=f),data=quad,alpha=0.2)+ geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 90,105)),col="red")+ geom_vline(aes(xintercept=x),data=data.frame(y="2FMSY",x=c( 80,110)),col="orange")+ scale_fill_manual(values=c("orange","orange","green","green","red","blue")) p5=p$layers[[5]] p$layers[[5]]=p$layers[[1]] p$layers[[1]]=p5 p
Figure r iFig=iFig+1; iFig. Time series of indicators compared to $F:F_{MSY}$, vertical lines indicate 1 (green), 1.5 (orange) and 2 (red) times $F_{MSY}$.
\newpage
Receiver Operating Characteristics
Detection of overfishing
simple_roc <- function(labels, scores){ labels <- labels[order(scores, decreasing=TRUE)] data.frame(TPR=cumsum(labels)/sum(labels), FPR=cumsum(!labels)/sum(!labels), labels, scores=sort(scores))} ind1=ldply(rf1[-1],function(x) as.data.frame(window(x,start=70,end=100),drop=T)) ind1=ddply(ind1,.(.id),transform,overfished=as.numeric(!(year%in%90:105))) ind1=ddply(ind1,.(.id),with, simple_roc(overfished,data)) tab=ddply(ind1,.(.id), with, { best=max(TPR*(1-FPR))==TPR*(1-FPR); data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])}) ggplot(ind1)+ geom_point(aes(FPR,TPR,col=.id))+ geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+ geom_point(aes(FPR,TPR),data=tab,size=2)+ xlab("False Negative Rate")+ylab("True Positive Rate")+ scale_color_manual("Indicator",values=rainbow(4))+ theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the four indicators of overfishing, points indicate the optimum value of the indicator.
tab
ind2=ldply(rf2[-1],function(x) as.data.frame(window(x,start=70,end=100),drop=T)) ind2=ddply(ind2,.(.id),transform,overfished=as.numeric(!(year%in%90:105))) ind2=ddply(ind2,.(.id),with, simple_roc(overfished,data)) tab=ddply(ind2,.(.id), with, { best=max(TPR*(1-FPR))==TPR*(1-FPR); data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])}) ggplot(ind2)+ geom_point(aes(FPR,TPR,col=.id))+ geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+ geom_point(aes(FPR,TPR),data=tab,size=2)+ xlab("False Negative Rate")+ylab("True Positive Rate")+ scale_color_manual("Indicator",values=rainbow(5))+ theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the three indicators of overfishing, points indicate the optimum value of the indicator.
tab
Detection of recovery
ind1=ldply(rf1[-1],function(x) as.data.frame(window(x,start=105),drop=T)) ind1=ddply(ind1,.(.id),transform,overfished=as.numeric((year>110))) ind1=ddply(ind1,.(.id),with, simple_roc(overfished,data)) tab=ddply(ind1,.(.id), with, { best=max(TPR*(1-FPR))==TPR*(1-FPR); data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])}) ggplot(ind1)+ geom_point(aes(FPR,TPR,col=.id))+ geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+ geom_point(aes(FPR,TPR),data=tab,size=2)+ xlab("False Negative Rate")+ylab("True Positive Rate")+ scale_color_manual("Indicator",values=rainbow(4))+ theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the four indicators of recovery, points indicate the optimum value of the indicator.
tab
ind2=ldply(rf2[-1],function(x) as.data.frame(window(x,start=105),drop=T)) ind2=ddply(ind2,.(.id),transform,overfished=as.numeric((year>105))) ind2=ddply(ind2,.(.id),with, simple_roc(overfished,data)) tab=ddply(ind2,.(.id), with, { best=max(TPR*(1-FPR))==TPR*(1-FPR); data.frame(TPR=TPR[best],FPR=FPR[best],scores=scores[best])}) ggplot(ind2)+ geom_point(aes(FPR,TPR,col=.id))+ geom_line(aes(x,y),data.frame(x=seq(0,1,0.01),y=seq(0,1,0.01)))+ geom_point(aes(FPR,TPR),data=tab,size=2)+ xlab("False Negative Rate")+ylab("True Positive Rate")+ scale_color_manual("Indicator",values=rainbow(5))+ theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig. ROC curve of the three indicators of recovery, points indicate the optimum value of the indicator.
tab
\newpage
References {#References}
\newpage
Session Info
sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.