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

[Marc Taylor]

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) 

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