Non-linear Discrete von Bertalanffy Growth Function

Interfaces brms
, ,
1

I am wondering if anyone has any suggestions for how I could go about creating a non-linear model for discrete growth in brms. A colleague (much wiser than I) has simulated the data and it appears to function well, I just cannot find a way to translate it into brms.

#Create simulated data for 5 IDs
sim_data <- tibble(
  ID = rep(1:5, each = 5),
  Year = rep(2008:2012, times = 5),
  Age = rep(1:5, times = 5), #Ages 1:5, assuming age 0 has size 0
  Capture_date = rep(NA, 25), 
  L = rep(NA, 25), #Length
  L_prev = rep(NA, 25),
  delta_t = rep(NA, 25) #Some change in time
)

sim_data <- sim_data %>%
  group_by(Year) %>%
  mutate(Capture_date = as.Date(paste0(Year, "-08-01")) + 
           sample(0:90, size = n(), replace = TRUE)) %>% #Random select some dates they were captured
  group_by(ID) %>% 
  mutate(L = ifelse(Age == 1, sample(200:270),
                    ifelse(Age == 2, sample (370:430),
                           ifelse(Age == 3, sample (430:460), 
                                  ifelse(Age == 4, sample (460:500), 
                                         ifelse(Age == 5, sample (500:510), NA))))), #Radnom select some Lengths
         L_prev = lag(L), #previous leg measurement
         delta_t = difftime(Capture_date, lag(Capture_date), units = "days"), 
         delta_t = as.numeric(delta_t)/365.25) %>% #convert to years due to limitations of difftime() function
  ungroup()


vb_discrete <- function(L_prev, L_max, r, deltat){
  L_now <- L_prev*exp(-r * deltat) + L_max * (1 - exp(-r * deltat))
  return(L_now)
} #Basic discrete VB growth function 

loop_down_vb_discrete <- function(size_vec, Lmax, r, deltaT) {
  L <- size_vec
  
  for(j in 2:length(L)){
    if (!is.na(L[j-1]) & !is.na(deltaT[j])) {
      L[j] <- vb_discrete(L_prev = L[j-1], 
                          L_max = Lmax,
                          r = r, 
                          deltat = deltaT[j])
    }
  } 
  return(L)
} #second loop to take previous estimated Length as start point for the following measure


sim_data2<- sim_data%>% 
  mutate(L2 = loop_down_vb_discrete(L,
                                    Lmax = 570, #Assume individuals reach asymptote at  570
                                    r = .6,
                                    deltaT = delta_t))

#Though this is wrong, my model would look something like this, where I can loop in the previously estimate value at each time sequence, similar to the simulation:

Discrete_vb_equation <- bf(Lnow ~ exp(Lprev) * exp(-exp(loggrowthrate) * deltat) + exp(logLmax ) * (1 - exp(-exp(loggrowthrate) * deltat)),
                                          Lprev~ 1 + (1|r|ID), 
                                          loggrowthrate ~ 1 + (1|r|ID), 
                                          logLmax ~ 1 + (1|r|ID),
                                          nl = TRUE,
                                          family = gaussian())

Also: Is it possible to estimate missing values in non-linear brms models? I was having difficulties and was suggested to use raw stan, which is beyond my skill levels.

  • Operating System: Windows 10x64 build 22621
  • brms Version: ‘2.18.0’
2

If you are comfortable with R your skills should be enough to write it in raw Stan (a.k.a. Stan). I don’t normally use R, but it looks like it’s a discrete-time model that could be implemented as a deterministic model in any language as it is or as an ODE model. I’m not sure it fit can be done in brms (but maybe it can, since it now does so much beyond linear models), but it can almost certainly be implemented in Stan, the only caveat is it must be a deterministic model (although there are some rare exceptions).

The code helps, but if you have a specific model in mind it could be helpful to write out the mathematical formulation, so anyone could understand its structure, and then translate it to Stan (or brms) or any other language. There are plenty of brms wizards here on Discourse, and they’d be able to tell your write away if your model can be written as a brms model or not.