Option to include fitted models (GLM, GAM) into a FLSR-object

[BernhardKuehn]

How to include general SR-models (besides BH, ricker, hockey...) like already fitted GLMs or GAMs into the FLSR-framework? For now everything seems to be strongly focused on providing a formula, which can be fitted with "fmle" or "nls".

If you provide anything other to the @model-slot in the FLSR-object than a "formula" it throws the following error:
"assignment of an object of class “...” is not valid for @‘model’ in an object of class “FLSR”; is(value, "formula") is not TRUE"

My idea is to get around this by not using the internal fitting routines to get my SR-model, but to just provide a fitted model as basis for further predictions. One way I thought about is to overwrite the slots the FLSR-object gets its predictions from, but so far I have no clue where to look for in the code.

Thanks for your help.

[iagomosqueira]

The model slot needs to be of class formula, yes, but that still allows you to use any function for prediction. If you can create a function that is able to predict recruitment based on any metric from the stock, predict(FLModel) will work on it.

For example, model could have rec ~ mygam(ssb) where mygam is a function that returns rec when called with an ssb input. If you want to store the model params there, rather than having them hard-coded in the function, you can also do so as long as you have them as function arguments rec ~ mygam(ssb, param1, param2).

Would this work for your models?

Note this will make it work in FLSR, but not for future projections using FLash or FLasher. The SRR models are hard coded in the C++ code, as we could not find any other way. Adding SRR models to FLasher is not that hard, but it is not trivial.

If you fit a GLM or GAM to SR data, what kind of formulation do you use to predict future recruitment?

[BernhardKuehn]

Thank you very much! I managed to make it work based on your suggestions and included some example code:

library(FLCore)
library(mgcv)

create some Stock-Recruitment-data with 2 environmental covariates

n = 30
set.seed(123)
ssb <- FLQuant(runif(n, 10000, 100000))
# environmental variable with positive autocorrelaion + trend
env1 <- FLQuant(seq(0.1,5,length.out = n) + arima.sim(n = n,model = list(order = c(1,0,0),ar = 0.7,sd = 0.1)),units = "env1")
# environmental variable with negative autocorrelation
env2 <- FLQuant(arima.sim(n = n,model = list(order = c(1,0,0),ar = -0.5,sd = 2)),units = "env2")
# recruitment model
rec <- exp(log(2) + 0.5*log(ssb) + 0.2*env1 -  0.4*env2^2 + rnorm(n,0,0.3))

create the FLSR-object

SR <- FLSR(ssb = ssb, rec = rec,covar = FLQuants(covar1 = env1,covar2 = env2)) 

Fit a gam to the stock-recruitment relationship

dat <- data.frame(rec = as.vector(rec),
                 ssb = as.vector(ssb),
                 covar1 = as.vector(env1),
                 covar2 = as.vector(env2))

gam.fit <- gam(rec ~ s(ssb,k = 4, bs = "tp") + 
                s(covar1,k = 4, bs = "tp")+ 
                s(covar2,k = 4, bs = "tp"),data = dat,family = gaussian(link = "log"))
par(mfrow = c(1,3))
plot(gam.fit)
summary(gam.fit)

Include the GAM in the FLSR-object

write the model as function

mygam = function(ssb,covar1,covar2){
  # convert to vector
  new.data = data.frame(ssb = as.vector(ssb),covar1 = as.vector(covar1),covar2 = as.vector(covar2))
  new.rec = predict(gam.fit,newdata = new.data,type = "response")
  # convert back to FLQuant-obj
  new.rec = FLQuant(as.vector(new.rec))
  
  return(new.rec)
}

include as formula

SR@model = as.formula(rec ~ mygam(ssb,covar))
predict(SR)

[iagomosqueira]

Glad it worked. A small suggestion would be to take care of the recruitment lag in the dimnames of the returned FLQuant. For example for a rec.age of 1

  new.rec = FLQuant(as.vector(new.rec),
    dimnames=list(year=as.numeric(dimnames(ssb)$year) + 1))