ices-tools-prod / msy

A collection of methods to estimate equilibrium reference points for fish stocks
GNU General Public License v2.0
13 stars 11 forks source link

exctraction of uncertainty around stock-recruitment relationship #8

Closed Isabella84 closed 8 years ago

Isabella84 commented 8 years ago

Is it possible to extract quantitively the uncertainty around the stock-recruitment relationship? For example CI of the parameters? Thank you

einarhjorleifsson commented 8 years ago

check out the sr.sto dataframe returned as one of the list object when calling eqsr_fit. einar

Isabella84 commented 8 years ago

Thank you to your answer Einar, I had found it already; I saw that for each run there is a CV associated. But what I need are "the blue lines" in the eqsr_plot. prova dps 9-10 eqsim Isabella

einarhjorleifsson commented 8 years ago

check out e.g. this code that resides within the eqsr_plot function

x <- fit$sr.det
    ssb <- seq(1,round(max(maxSSB)),length=100)
    z <- sapply(1:nrow(x), function(i) rec <- exp(match.fun(as.character(x$model[i])) (x[i,], ssb)))
    modelLines <- as.data.frame(cbind(ssb,z))
    names(modelLines) <- c("ssb",paste(x$model,x$prop))
    modelLines <- reshape2::melt(modelLines,id.var="ssb",variable.name="Model",value.name="rec")

    out$ssb <- out$ssb/Scale
    out$rec <- out$rec/Scale
    out$mid.grp <- out$mid.grp/Scale
    Percentiles$ssb <- Percentiles$ssb/Scale
    Percentiles$p50 <- Percentiles$p50
    Percentiles$p05 <- Percentiles$p05
    Percentiles$p95 <- Percentiles$p95

    modelLines$ssb <- modelLines$ssb/Scale
    modelLines$rec <- modelLines$rec/Scale

    fit$rby$ssb <- fit$rby$ssb/Scale
    fit$rby$rec <- fit$rby$rec/Scale
    i <- sample(nrow(out),n)

ggplot2::ggplot(out[i,]) + 
      ggplot2::theme_bw() +
      ggplot2::geom_point(ggplot2::aes(x=ssb,y=rec,colour=Model),size=1) +
      ggplot2::geom_line(data=Percentiles,ggplot2::aes(x=ssb,y=p05),colour="yellow") +
      ggplot2::geom_line(data=Percentiles,ggplot2::aes(x=ssb,y=p95),colour="yellow") +
      ggplot2::geom_line(data=Percentiles,ggplot2::aes(ssb,p50),col="yellow",lwd=2) +
      ggplot2::geom_line(data=modelLines,ggplot2::aes(ssb,rec,colour=Model),lwd=1) +
      ggplot2::coord_cartesian(ylim=c(0,quantile(out$rec[i],0.99))) +
      ggplot2::geom_path(data=fit$rby,ggplot2::aes(ssb,rec),col="black",linetype=2) +
      ggplot2::geom_text(data=fit$rby,ggplot2::aes(ssb,rec,label=substr(year,3,4)),size=4,col="black",angle=45) +
      ggplot2::theme(legend.position = c(0.20,0.85)) +
      ggplot2::labs(x="Spawning stock biomass",y="Recruitment",colour="Model")

the last step is just the plot, but what precedes it this should have the elements you are interested in. einar

Isabella84 commented 8 years ago

Thank you very much Einar. For me, you can close the issue.