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

[Marc Taylor]

Thanks Laurie. I would greatly appreciate the example if you are able to 
find it.

Cheers,
Marc


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 9:49 AM, laurie wrote:
> see, which is in FLCore
>
>
> rickerCa <- function() {
>   logl <- function(a, b, c, rec, ssb, covar)
>     loglAR1(log(rec), log(a * (1 - c * covar) * ssb * exp(-b * ssb)))
>
>   initial <- structure(function(rec, ssb) {
>         # The function to provide initial values
>     res  <-coefficients(lm(c(log(rec/ssb))~c(ssb)))
>     return(FLPar(a=max(exp(res[1])), b=-max(res[2]), c=1))},
>
>   # lower and upper limits for optim()
>     lower=rep(-Inf, 3),
>     upper=rep( Inf, 3))
>
>     model  <- rec ~ a * (1 - c * covar) * ssb * exp(-b * ssb)
>     return(list(logl=logl, model=model, initial=initial))
> }
>
>
> There is an example using it somewhere, I will see if I can find it.
>
> Laurie
>
> On 07/04/16 09:37, Marc Taylor wrote:
>> Dear FLR-Listers,
>>
>> I was hoping someone might help me figure out a way to create a 
>> custom SRModel that uses a covariate in addition to ssb and rec.
>>
>> From the documentation of SRModels (i.e. ?SRModels), reference is 
>> made to several models that would be of interest to me (e.g. 
>> 'ricker.c.a', 'ricker.c.b', etc.), but these do not appear to be 
>> presently included in FLCore.  I made an attempt (below) to define a 
>> custom model, but am getting errors when fitting with fmle().
>>
>> If anyone has any experience with such models, I would greatly 
>> appreciate any advise.
>>
>> Cheers,
>> Marc
>>
>> *Example script*:
>>
>> # packages 
>> ----------------------------------------------------------------
>> library(FLCore) # FLCore_2.5.20160107
>>
>>
>> # load data 
>> ---------------------------------------------------------------
>>
>> data("ple4")
>>
>>
>> # make srr obbject and add covariate 
>> --------------------------------------
>>
>> 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)
>>
>> ple4Sr at covar <- FLQuants(list(env1=FLQuant(env1, dim=dim(ple4Sr at ssb), 
>> dimnames=dimnames(ple4Sr at ssb))))
>> ple4Sr at covar$env1
>>
>>
>> # attempt with existing SRModel fails
>> model(ple4Sr) <- "ricker.c.a"
>> # Error in do.call(value, list()) : could not find function "ricker.c.a"
>>
>>
>>
>> # make new srr function 
>> ---------------------------------------------------
>>
>> # adapt ricker
>> ricker <- function(){
>>     logl <- function(a, b, rec, ssb) loglAR1(log(rec), log(a *
>>         ssb * exp(-b * ssb)))
>>     initial <- structure(function(rec, ssb) {
>>         res <- coefficients(lm(log(c(rec)/c(ssb)) ~ c(ssb)))
>>         return(FLPar(a = max(exp(res[1])), b = -max(res[2])))
>>     }, lower = rep(-Inf, 2), upper = rep(Inf, 2))
>>     model <- rec ~ a * ssb * exp(-b * ssb)
>>     return(list(logl = logl, model = model, initial = initial))
>> }
>>
>> # to this (incl. env. control term - 'env1')
>> ri.ec <- function(){
>>     logl <- function(a, b, e1, rec, ssb, env1) loglAR1(log(rec), log(a *
>>         ssb * exp(-b * ssb) * exp(e1 * env1)))
>>     initial <- structure(function(rec, ssb, env1) {
>>         res <- coefficients(lm(log(c(rec)/c(ssb)) ~ c(ssb) + c(env1)))
>>         return(FLPar(a = max(exp(res[1])), b = -max(res[2], env1 = 
>> max(res[3]))))
>>     }, lower = rep(-Inf, 3), upper = rep(Inf, 3))
>>     model <- rec ~ a * ssb * exp(-b * ssb) * exp(e1 * env1)
>>     return(list(logl = logl, model = model, initial = initial))
>> }
>>
>>
>> # fit SRR 
>> -----------------------------------------------------------------
>> model(ple4Sr) <- "ri.ec"
>> ple4Sr <- fmle(ple4Sr)
>> # Error in loglAR1(log(rec), log(a * ssb * exp(-b * ssb) * exp(e1 * 
>> env1))) :
>> #   error in evaluating the argument 'hat' in selecting a method for 
>> function 'loglAR1': Error: argument "e1" is missing, with no default
>> summary(ple4Sr)
>> plot(ple4Sr)
>> -- 
>> 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
>>
>>
>> _______________________________________________
>> flr-list mailing list
>> flr-list at flr-project.org
>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/flr-list
>
>
>
> _______________________________________________
> flr-list mailing list
> flr-list at flr-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/flr-list

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/flr-list/attachments/20160407/c3b00167/attachment-0001.html>