0

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’

0 Answers0