Error when creating a method that uses RTMB to fit to data stored in and S4 object

[brandon...@gmail.com]

Hello,

I can't thank the developers enough for creating RTMB. Quite simply, it's everything TMB was, but a whole lot easier to debug. 

I didn't see anything in the issues associated with this. Simply put I'm creating a package that uses RTMB to estimate parameters of model. Using a couple of simple examples, I've found no issues with using RTMB in a method to estimate parameters of a model using the data stored an S4 object. However, when I use a loop to fill the values of a derived variable using parameters in the method, I get an error. In the example below, if I simply write chi <- rep(phi,nr), I do not get the error (albeit the model is incorrect).  The error is associated strictly with using a loop to fill the values in chi.

The error is: 

    # Error in advector(e2) :
    #   Invalid argument to 'advector' (lost class attribute?)

and I've paste in the short code below. 

I thought perhaps it was a scoping issue, so I replaced chi with this code,

    chi <- rep(1,nr)
    for(i in 1:nr){
      chi[i] <- phi
    }

Again, the model doesn't make sense, but it runs and produces estimates for the parameters, but only if I also change the likelihood to,

nll <-  nll - dbinom(1, 1, as.numeric(chi[i]), log = TRUE) * n[i]

and it gives the warning,

Warning messages: 1: In dbinom(1, 1, as.numeric(chi[i]), log = TRUE) : imaginary parts discarded in coercion 

Thank you in advance for an help on this, 

Brandon

#' Fit a Survival Model
#' Fits the survival model based on the provided mark-recapture data.
#'
#' @param object TaggingData object.
#' @return The fitted TagModel object.
#' @import RTMB
#' @export
setMethod("fit", "TagModel", function(object) {
  data <- object@data  # Capture history data, ensure it's a matrix
  

parameters <- list(f_p = 0, f_phi = 0)  # Initialize with an intercept value (log-odds)

  f <- function(parms) {
    getAll(data,parms)

    nll <- 0
    p <- plogis(f_p)
    phi <- plogis(f_phi)

    nc <- ncol(ch)
    nr <- nrow(ch)

    #This code bit here throws an error
    # Error in advector(e2) :
    #   Invalid argument to 'advector' (lost class attribute?)
    chi <- rep(1,nr)
    for(i in 1:nr){
      if(U[i]<nc){
        for(j in (nc-1):U[i]){
          chi[i] <- (1 - phi) + phi * (1 - p) * chi[i]
        }
      }
    }
   #end of problem
   
    for(i in 1:nr){
      if(U[i]>1){
        for(j in 2:U[i]){
          nll <-  nll - dbinom(1, 1, phi, log = TRUE) * n[i]
          nll <-  nll - dbinom(ch[i,j], 1, p, log = TRUE) * n[i]
        }
      }
      nll <-  nll - dbinom(1, 1, chi[i], log = TRUE) * n[i]
    }
    return(nll)  # Return the negative log-likelihood
  }

  obj <- MakeADFun(f, parameters)
  opt <- nlminb(obj$par, obj$fn, obj$gr)

  object@TMB$opt <- opt
  object@TMB$obj <- obj
  return(object)
})

 

[Kasper Kristensen]

Some useful (?) links.

Solve this issue using ADoverload:
https://groups.google.com/g/tmb-users/c/HlPqkfcCa1g

Mysterious behaviour is probably due to R's byte compiler:
https://groups.google.com/g/tmb-users/c/CkmgUr2WE40/m/x9zupcBGAwAJ

[brandon...@gmail.com]

Hello Kasper,

Thank you for the response. I found a very jinky solution. Within the method that creates the AD function and does the optimization for the S4 tmb object,  I've passed the tmb lists for the data and parameters, as well as the model function to the global environment. I will compare this solution to the links you've sent and try and come up with something a more elegant.

Thank you for your help and suggestions and I've attached a very simple solution below. There is a link to this simple package example here, https://github.com/bchasco/S4exampleTMB.

# Generic build method
setGeneric("build_tmb", function(object) standardGeneric("build_tmb"))

#' @param object A data object.
#' @return The tmb objective function output
#' @export
# Method building and runnning the tmb model
setMethod("build_tmb", "tmb_list", function(object) {
  tmb.data <<- list() #create in .GlobalEnv
  for(i in names(object@TMB$data)){
    tmb.data[[i]] <<-  object@TMB$data[[i]]
  }
  parameters <<- list()
  for(i in names(object@TMB$parameters)){
    parameters[[i]] <<-  object@TMB$parameters[[i]]
  }
  #pass function to global environment
  environment(f) <- .GlobalEnv
  obj <- RTMB::MakeADFun(f,

                         parameters)
  opt <- nlminb(obj$par, obj$fn, obj$gr)

  object@TMB$obj <- obj
  object@TMB$opt <- opt
  return(object)  # Return the updated object with results
})

#' @slot TMB
#' @export
setClass("tmb_list",
         slots = list(
           TMB = 'list'
         )
)

#Model function
f <- function(parms){
  RTMB::getAll(tmb.data,
               parms)
  nll <- -sum(RTMB::dnorm(y, beta, exp(sigma),log=TRUE))
  return(nll)
}

#create the class list

 tmb_list <- new("tmb_list",
                TMB = list(data = list(y = rnorm(100)),
                           parameters = list(beta = 0,
                                             sigma = 0))
)

#use the method to optimize
tmb_out <- build_tmb(tmb_list)
print(tmb_out@TMB$opt)

$opt $opt$par beta sigma 0.02425354 -0.01602493 $opt$objective [1] 140.2914 $opt$convergence [1] 0 $opt$iterations [1] 6 $opt$evaluations function gradient 11 7 $opt$message [1] "relative convergence (4)"