Commit 103f45b3 authored by Berges, Benoit's avatar Berges, Benoit
Browse files

push changes to SRR estimation

parent 4afc2e8f
......@@ -3,9 +3,9 @@ library(TMB)
rm(list=ls())
path <- "C:/git/sanoba_nsas"
#path <- "C:/git/sanoba_nsas"
try(setwd(path),silent=TRUE)
#try(setwd(path),silent=TRUE)
################################################
# setup paths
......@@ -15,6 +15,7 @@ data.dir <- file.path(".","data")
script.dir <- file.path(".","side_scripts")
function.dir <- file.path(".","functions")
results.dir <- file.path(".","results")
model.dir <- file.path(".","model")
figures.dir <- file.path(".","figures")
runName <- '01b_SRR_fit'
......@@ -25,6 +26,11 @@ dir.create(file.path(figures.dir,runName),showWarnings = FALSE)
results.dir <- file.path(".","results",runName)
figures.dir <- file.path(".","figures",runName)
source(file.path(function.dir,'rickerCb.R'))
source(file.path(function.dir,'segRegBlim.R'))
load(file.path(model.dir,'NSAS_sanoba_assessment_results_sf.RData'))
################################################
# read data
################################################
......@@ -34,30 +40,27 @@ dat <- read.csv(file.path(results.dir,'..','00_NSAS_SST','SST_NSAS.csv'))
SSB_pred <- seq(from=min(dat$ssb),to=max(dat$ssb),by=1000)
SST_pred <- seq(from=min(dat$sst),to=max(dat$sst),by=0.5)#min(dat$sst)#
################################################
# compile TMB function
################################################
compile(file.path(function.dir,"C++","ricker_T.cpp"))
dyn.load(dynlib(file.path(function.dir,"C++","ricker_T")))
data <- list(SSB=dat$ssb,
logR=log(dat$rec),
SST=dat$sst,
SSB_pred=SSB_pred,
SST_pred=SST_pred)
parameters <- list(logA=0,
logB=0,
logC=0,
logSigma=0)
obj <- MakeADFun(data,parameters,DLL="ricker_T")
opt <- nlminb(obj$par,obj$fn,obj$gr)
rep <- sdreport(obj)
logR <- rep$par.fixed[1]+log(SSB_pred)-exp(rep$par.fixed[2])*SSB_pred+exp(rep$par.fixed[3])*SST_pred
plot(dat$ssb,log(dat$rec), main = "Ricker", cex.main = 2,cex.lab = 1.5,
xlab = "SSB",ylab = "R")
lines(SSB_pred, logR, lwd=3)
lines(SSB_pred, logR+2*rep$sd, lty="dotted")
lines(SSB_pred, logR-2*rep$sd, lty="dotted")
dat$sstResi <- (dat$sst-mean(dat$sst))/sd(dat$sst)
myFLSR <- as.FLSR(NSH)
myFLSR@covar <- FLQuants(covar=FLQuant(dat$sstResi,dimnames=list(year=dat$year)))
#model(myFLSR) <- rickerCb
#model(myFLSR) <- ricker
#model(myFLSR) <- segRegBlim
myFLSR <- fmle(myFLSR)
rickerParams <- myFLSR@params
model(myFLSR) <- rickerCa
myFLSR <- fmle(myFLSR,fixed=list(a=47.6))
rickerCaParams <- myFLSR@params
model(myFLSR) <- rickerCb
myFLSR <- fmle(myFLSR)
predict(myFLSR,ssb=FLQuant(c(250000,100000)),covar=FLQuant(c(0.2,0.1)))
library(FLCore)
library(TMB)
rm(list=ls())
#path <- "C:/git/sanoba_nsas"
#try(setwd(path),silent=TRUE)
################################################
# setup paths
################################################
data.dir <- file.path(".","data")
script.dir <- file.path(".","side_scripts")
function.dir <- file.path(".","functions")
results.dir <- file.path(".","results")
model.dir <- file.path(".","model")
figures.dir <- file.path(".","figures")
runName <- '01b_SRR_fit'
dir.create(file.path(results.dir,runName),showWarnings = FALSE)
dir.create(file.path(figures.dir,runName),showWarnings = FALSE)
results.dir <- file.path(".","results",runName)
figures.dir <- file.path(".","figures",runName)
source(file.path(function.dir,'rickerCb.R'))
################################################
# read data
################################################
dat <- read.csv(file.path(results.dir,'..','00_NSAS_SST','SST_NSAS.csv'))
SSB_pred <- seq(from=min(dat$ssb),to=max(dat$ssb),by=1000)
SST_pred <- seq(from=min(dat$sst),to=max(dat$sst),by=0.5)#min(dat$sst)#
dat$sstResi <- (dat$sst-mean(dat$sst))/sd(dat$sst)
load(file.path(model.dir,'NSAS_sanoba_assessment_results_sf.RData'))
myFLSR <- as.FLSR(NSH)
myFLSR@covar <- FLQuants(covar=FLQuant(dat$sstResi,dimnames=list(year=dat$year)))
#model(myFLSR) <- rickerCb
model(myFLSR) <- ricker
myFLSR <- fmle(myFLSR)
rickerParams <- myFLSR@params
model(myFLSR) <- rickerCa
myFLSR <- fmle(myFLSR,fixed=list(a=47.6))
rickerCaParams <- myFLSR@params
model(myFLSR) <- rickerCb
myFLSR <- fmle(myFLSR)
predict(myFLSR,ssb=FLQuant(c(250000,100000)),covar=FLQuant(c(0.2,0.1)))
################################################
# compile TMB function
################################################
compile(file.path(function.dir,"C++","ricker_T.cpp"))
dyn.load(dynlib(file.path(function.dir,"C++","ricker_T")))
data <- list(SSB=dat$ssb,
logR=log(dat$rec),
SST=dat$sstResi,
SSB_pred=SSB_pred,
SST_pred=SST_pred)
data <- list(SSB=dat$ssb,
logR=log(dat$rec),
SST=dat$sstResi)
parameters <- list(logA=0,
logB=0,
logC=0,
logSigma=0)
obj <- MakeADFun(data,parameters,DLL="ricker_T")
opt <- nlminb(obj$par,obj$fn,obj$gr)
rep <- sdreport(obj)
logR <- rep$par.fixed[1]+log(SSB_pred)-exp(rep$par.fixed[2])*SSB_pred+exp(rep$par.fixed[3])*SST_pred
plot(dat$ssb,log(dat$rec), main = "Ricker", cex.main = 2,cex.lab = 1.5,
xlab = "SSB",ylab = "R")
lines(SSB_pred, logR, lwd=3)
lines(SSB_pred, logR+2*rep$sd, lty="dotted")
lines(SSB_pred, logR-2*rep$sd, lty="dotted")
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_VECTOR(SSB)
DATA_VECTOR(SST)
DATA_VECTOR(logR);
PARAMETER(logA);
PARAMETER(logB);
PARAMETER(logC);
PARAMETER(logSigma);
vector<Type> pred_model=logA+log(1-exp(logC)*SST)+log(SSB)-exp(logB)*SSB;
Type nl=-sum(dnorm(logR,pred_model,exp(logSigma),true));
ADREPORT(pred_model)
/*ADREPORT(pred_model)*/
return nl;
}
rickerCb <- function (){
logl <- function(a, b, c, rec, ssb, covar) loglAR1(log(rec),
log(a * ssb * exp(-b * ssb+c*covar)))
initial <- structure(function(rec, ssb) {
res <- coefficients(lm(c(log(rec/ssb)) ~ c(ssb)))
return(FLPar(a = max(exp(res[1])), b = -max(res[2]),
c = 1))
}, lower = rep(-Inf, 3), upper = rep(Inf, 3))
model <- rec ~ a * ssb * exp(-b * ssb+c*covar)
return(list(logl = logl, model = model, initial = initial))
}
segRegBlim <- function ()
{
logl <- function(a, rec, ssb) {
loglAR1(log(rec), FLQuant(log(ifelse(c(ssb) <= 874198, a *
c(ssb), a * 874198)), dimnames = dimnames(ssb)))
}
model <- rec ~ FLQuant(ifelse(c(ssb) <= 874198, a * c(ssb), a *
874198), dimnames = dimnames(ssb))
initial <- structure(function(rec, ssb) {
return(FLPar(a = median(c(rec)/c(ssb), na.rm = TRUE)))
}, lower = rep(0, 0), upper = rep(Inf, 2))
return(list(logl = logl, model = model, initial = initial))
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment