[FLR-list] fitting a custom SRModel with a covariate

Mon Apr 11 05:44:52 CEST 2016

Hi,

If you end up using it with some covariates that improve the fit, please 
do send me the code, as I would like to have an example for testing.

This is what we have right now on the model in the reviewed help file, 
still to be updated in the compiled package, avialable at

https://github.com/flr/FLCore/blob/topic/roxygen/R/SRmodels.R

#' Ricker model with one covariate. The covariate can be used, for 
example, to
#' account for an enviromental factor that influences the recruitment 
dynamics.
#' In the equations, \emph{c} is the shape parameter and \emph{X} is the
#' covariate.
#'
#' \itemize{
#' \item rickerCa: Ricker stock-recruitment model with one
#' multiplicative covariate.
#' \deqn{R = a (1- c X) S e^{-b S}}{R = a*(1-c*X)*S*e^{-b*S}} }

If you think we should be more precise, please feel free to suggest changes.

Thanks,

Iago

On 08/04/16 18:28, Marc Taylor wrote:
> Dear Iago,
>
> Thank you for updating this functionality. The fit does indeed converge
> now.
>
> This covariate was indeed arbitrary, so I have now provided an example
> below where the covariate is constructed based on the residuals of the
> ricker() fit. So, the rickerCa() should now make a significant
> improvement over the fit of the ricker(), but this is not clear by just
> looking at AIC.
>
> I'm also wondering where this model is described? I have come across
> another version of an adapted Ricker where the covariate term is
> exp(c*covar) rather than (1-c*covar). As is, the fitted rickerCa
> produces NAs.
>
> I've tried this alternate version and the fit looks quite good (as well
> as a substantially lower AIC).
>
> Anyway, thanks for fixing this so promptly.
>
> Cheers,
> Marc
>
> p.s. Perhaps you are now updating the documentation for ?SRModels - At
> the moment I am not able to access its help page in the new FLCore version.
>
> ###
>
> library(FLCore) # FLCore_2.6.0.20160401
> data("ple4")
>
>
>
> # fit basic ricker --------------------------------------------------------
> ple4SR <- as.FLSR(ple4)
> model(ple4SR) <- "ricker"
> ple4SR <- fmle(ple4SR)
> summary(ple4SR)
> plot(ple4SR)
>
>
> # create environmental covariate with residuals ---------------------------
> set.seed(1)
> env <- c(2 * residuals(ple4SR) + rnorm(dims(ple4SR at ssb)$year, 0, 1))
> env
> plot(c(env), c(residuals(ple4SR)))
> plot(c(env), t="l")
>
>
> # make new SR obj with covariate ------------------------------------------
> ple4SR2 <- as.FLSR(ple4)
> covar(ple4SR2) <- FLQuants(covar=FLQuant(env,
> dimnames=dimnames(ssb(ple4SR2))))
>
>
> # fit rickerCa
> model(ple4SR2) <- rickerCa
> ple4SR2 <- fmle(ple4SR2)
> plot(ple4SR2)
> params(ple4SR2)
> plot(c(covar(ple4SR2)$covar), c(residuals(ple4SR2)))
>   c(residuals(ple4SR2)) # produces NAs
>
>
>
> # Compare AIC -------------------------------------------------------------
> AIC(ple4SR)
> AIC(ple4SR2)
>
>
>
> # alt model (covariate term is exp(c*covar), initial=0)
> ---------------------------------------------------------------
>
> # compare nls
> START <- list(a=params(ple4SR)["a"], b=params(ple4SR)["b"], c=0)
> R <- c(ple4SR2 at rec)
> S <- c(ple4SR2 at ssb)
> covar <- c(covar(ple4SR2)$covar)
> fitnls <- nls(log(R) ~ log(a) + log(S) - b*S + c*covar, start=START)
> summary(fitnls)
> plot(log(R), predict(fitnls)); abline(0,1)
> plot(S, resid(fitnls)); abline(h=0)
> plot(covar, resid(fitnls)); abline(h=0)
>
>
> rickerCa.alt <- function ()
> {
>      logl <- function(a, b, c, rec, ssb, covar) loglAR1(log(rec),
>          log(a * ssb * exp(-b * ssb) * exp(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 = 0))
>      }, lower = rep(-Inf, 3), upper = rep(Inf, 3))
>      model <- rec ~ a * ssb * exp(-b * ssb) * exp(c*covar)
>      return(list(logl = logl, model = model, initial = initial))
> }
>
> ple4SR3 <- as.FLSR(ple4)
> covar(ple4SR3) <- FLQuants(covar=FLQuant(env,
> dimnames=dimnames(ssb(ple4SR3))))
> model(ple4SR3) <- rickerCa.alt
> ple4SR3 <- fmle(ple4SR3, start=as.list(coefficients(fitnls)) )
> plot(ple4SR3)
> params(ple4SR3)
> coefficients(fitnls)
>
> AIC(ple4SR)
> AIC(ple4SR2)
> AIC(ple4SR3)
>
>
>
>
> Achtung: Das Thünen-Institut hat die Domain gewechselt. Bitte ändern Sie
> meine Mailadresse in Ihrem Adressbuch!
> Notice: Thünen Institute has changed its domain. Please change my email
> address in your address book!
>
> Dr. Marc Taylor
> Marine Lebende Resourcen / Marine Living Resources
> Thünen-Institut für Seefischerei / Thünen Institute of Sea Fisheries
> Palmaille 9
> 22767 Hamburg, Germany
>
> Tel:  +49 40 38905-143
> Mail: marc.taylor at thuenen.de
> Web:  www.ti.bund.de
>
> On 4/7/2016 4:46 PM, Iago MOSQUEIRA (JRC) wrote:
>>
>> library(FLCore)
>>
>> data("ple4")
>>
>> ple4SR <- as.FLSR(ple4)
>>
>> env1 <- c(-0.81, -1.73, -0.09, 0.42, -1.01, -1.06, -0.79, -0.42,
>> -0.74, -0.3, 0.25, 0.09, -0.18, -0.29, 0.13, 0.09, 0.48, 0.33, 0.41,
>> 0.4, -0.21, -0.64, -0.97, -0.02, -0.19, 0.28, -0.03, -0.19, -0.47,
>> -0.97, -0.88, 0.23, 0.82, 0.92, 0.17, 0.67, -0.07, 0.32, 0.9, -0.64,
>> 0.97, 0.54, 1.18, 0.67, 0.44, 1.29, 1.54, 1.22, 0.47, 1.04, 1.31, 1.06)
>>
>> covar(ple4SR) <- FLQuants(covar=FLQuant(env1,
>> dimnames=dimnames(ssb(ple4SR))))
>>
>> model(ple4SR) <- rickerCa
>>
>> ple4SR <- fmle(ple4SR)
>>
>> plot(ple4SR)
>
>
>
>
> _______________________________________________
> flr-list mailing list
> flr-list at flr-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/flr-list
>

-- 
Dr Iago Mosqueira

European Commission – Joint Research Center
Institute for the Protection and Security of the Citizen (IPSC)
Maritime Affairs Unit G03
TP 051, Via Enrico Fermi 2749
I-21027 Ispra (VA), Italy

Office : +39 0332 785413
Fax: +39 0332 789658

iago.mosqueira at jrc.ec.europa.eu

http://fishreg.jrc.ec.europa.eu/home