Hierarchical Non-linear Gompertz growth curve in BRMS

Modeling
,
1

I am trying to model growth of bird chicks over two seasons. The idea is to run the model on all data then to use post-hoc comparisons to contrast the growth rates between the two years. I am including a grouping variable (nestID) at the individual level because there are multiple samples from each individual. I have managed to get the model to converge but I am unsure if I have the formula structured correctly to answer the above question.

I interpret the Gompertz function as follows
y(t) = A * exp(-B * exp(-C * t))

where:
y(t) represents the value of the curve at time t (note: time is morphometric data in this case i.e. length of head and culmen which grows at a consistent rate)
A is the assymptote
B is a constant that determines displacement on the x axis (intercept)
C is a constant that determines the rate of change of the curve (growth rate)
y in this case is the weight of the chick

My two concerns here are:

  1. Do my linear formulas that are used to estimate each of the Gompertz function’s parameters accurately represent the question I am asking and the structure of the data? Particularly, does the grouping term (1|nestID) make sense for just the linear formula pertaining to the intercept parameter (B) or should it be calculated for all parameters? Similarly, do I need to include the categorical variable for year in just the parameter of interest (C, growth rate) or should it too be included in all linear predictors for the Gompertz function’s parameters?

So for the linear predictors, replace my current code with something like:
A ~ year + (1|nestID),
B ~ year + (1|nestID) ,
C ~ year + (1|nestID)

  1. My ultimate goal is to contrast the change in growth rate between years. Does the linear model convention work the same for non-linear models in BRMS?

E.g. something like hypothesis(model, “year2022 > year2021”)?


library(brms)

# Formula
form <- bf(finalBrdWt ~ A * exp( B * exp( C * head.culmen ) ),
             A ~ 1 ,
             B ~ 1 + (1 | nestID),
             C ~ year ,
             nl = TRUE)

# Priors
priors <- c(
  prior(normal(0.5, 2), nlpar = "A", lb=0.1),
  prior(normal(-5.6, 0.5), nlpar = "B", lb=-7, ub=-5.5),
  prior(normal(-2, 1), nlpar = "C", ub=0),
  prior(normal(0, 1), class = "sigma"))

brm.headcul<- brm(
  formula = form, 
  data = wts,
  prior = priors, 
  seed = 9987,
  family = lognormal(link = "identity", link_sigma = "identity"),
  chains = 3, 
  #warmup = 100,
  cores = 3, 
  iter = 3000,
  sample_prior = "yes")
2

Hello, if I’m understanding you correctly, here are some comments on your problem:

  1. I’ll presume that your implementation of the Gompertz function is correct and focus on the remainder of the model. The linear predictor expressions should be constructed based on both the question you’re asking, and your knowledge of the system at hand. For your current specification, A will be represented solely by a single estimated constant, C will be estimated by year (if it is binary, there will be two values of C), and B will have an estimated population value, with variation in around this value by nestID. The relevance of adding additional predictors for A and C as you have proposed is context-dependent. Consider for example if you want to propose, a priori, that the asymptote and growth rate are identical across both year and nestID, or conversely, if your design would support estimating these additional parameters with reasonable identifiability.

  2. I’m not sure how hypothesis() would work in this case as I don’t use this workflow. The way I would approach the effect of year would be either to summarize the effect of year on the relevant parameters on the linear predictor scale (this could be done via fixef(), or coef() to obtain summaries or posteriors), or to look at predictions of y(t) conditional on year, say via fitted() or predict().

  3. it’s probably not advisable to have so few warmups.

  4. that’s probably not enough priors. I imagine that you will want to set priors on all of the parameters of the linear predictors.

3

I’d just add that sometimes it’s helpful in these nonlinear formulas to wrap the parameters with exp() and inv_logit() if you know they need to either be positive or between 0 and 1. That lets your A ~ 1, B ~ 1, etc. estimates roam more freely without breaking down the computation. You won’t have to put upper or lower bounds on those priors.

For example, if you know A must positive since it’s an asymptote in a growth model, your formula might be:

finalBrdWt ~ exp(A) * exp( B * exp( C * head.culmen ) )

Additionally, the priors only make sense in the context of the entire model and the family. So I would definitely recommend looking at the output with the sample_priors = ‘only’ argument and adjust from there.