Setting simple binomial prior

Dear colleagues,

I am trying to fit a Baesyan von Bertalanffy curve using the code below:


#VON BERTALANFFY-3

#Set formula

formula <- brmsformula(LD ~ Linf * (1 - exp(-K * (TR - t0))),
                       Linf ~ 1, K ~ 1, t0 ~ 1, nl=TRUE)

#Set priors

priors <- prior(normal(0,1000), nlpar="Linf", lb=0) +
  prior(normal(0,1), nlpar="K", lb=0) +
  prior(normal( 0,10), nlpar="t0") +
  prior(student_t(3,0,40), class=sigma)

#Set starting values

initsLst <- function() list(
  Linf=runif(1, 200, 800),
  K=runif(1, 0.05, 3.00),
  t0=rnorm(1, 0, 0.5)
)

#Fit model

fit1 <- brm(formula,            # calls the formula object created above
            family=gaussian(),  # specifies the error distribution
            data=dat,         # the data object
            prior=priors,       # the prior probability object 
            init=initsLst,      # the initial values object
            chains=3,           # number of chains, typically 3 to 4
            cores=3,            # number of cores for multi-core processing. Typically set to
            # match number of chains. Adjust as needed to match the number
            # of cores on your computer and number of chains
            iter=3000,          # number of iterations
            warmup=1000,        # number of warm up steps to discard
            control=list(adapt_delta=0.80,    # Adjustments to algorithm to 
                         max_treedepth=15)) # improve convergence.

I need to set a binomial prior por the ā€œKā€ coefficient, but Iā€™m not getting it.

Thanks in advance for any help.

Regards,

Rodrigo

1 Like

Hi and welcome. I did a quick cleanup on your post using ``` to enclose your code. Itā€™s a bit easier to read that way.

I have a von B model setup in Stan but havenā€™t tried this in brms yet. I know there are a few Bayesian fish folks that hang out here and they might have some insight into how to do this is brms.

Thank you, Ara!

If you want some Stan code for this and simulated data let me know.

From your model formula, it looks like K is continuous. Why would you want to give it a bounded discrete prior like a binomial?

Not only are they discarded, like in burnin, the warmup iterations are used for adaptation in three phases (I: find bulk of probability mass, II: adapt mass matrix, III: fine-tune step size).

If you run default_prior(), giving it your model formula and data, it will show what the default priors are. With that information you should be able to see how to set a prior on K. Just look for the class associated with K and include that in the call to prior(). Note: I havenā€™t used the non-linear model specification in brms, so this might be off a little, but it should give you a decent starting point.

Thak you for your reply Bod, but ā€œKā€ is a coiefficient that I want the model to extract from other data in my dataset (namely, size, age and sex).

I am trying the binomial distribution because this coefficient is always <1 in other similar species of the animal I am working with (a stingray)

Thank you, I gonna try!

The original post was asking for a prior for K, and you had standard normal in the code, so I inferred that K was continuous. Iā€™m not sure where the binomial prior from the postā€™s title comes in.

Dear colleagues,

I tried to specify the binomial funcion in my code and I am getting the following error warning:

"Compiling Stan programā€¦
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
0
Semantic error in ā€˜stringā€™, line 29, column 12 to column 37:

27:    lprior += normal_lpdf(b_Linf | 0, 1000)
28:      - 1 * normal_lccdf(0 | 0, 1000);
29:    lprior += binomial_lpdf(b_K | 0, 1)
                 ^
30:      - 1 * binomial_lccdf(0 | 0, 1);
31:    lprior += normal_lpdf(b_t0 | 0, 10);

Function ā€˜binomial_lpdfā€™ is not implemented for distribution ā€˜binomialā€™, use ā€˜binomial_lpmfā€™ instead."

Iā€™m not sure exactly what is intended here, but there must be some confusion.

The binomial distribution places probability mass only over the integers greater than or equal to zero and less than or equal to N, one of the parameters of the binomial distribution.

There are numerous aspects of the your posts @RodrigoMartins that make me nearly certain that you do not actually want to put a binomial prior on K. These include:

  • In the brms model, K is a continuous quantity.
  • the inits function returns non-integer values for K
  • you say that K is always <1, but the only possible value for a binomially distributed K that is less than one is zero.
  • you write binomial_lpdf (which isnā€™t a function) instead of binomial_lpmf in your stan code.
  • you write binomial_lpdf(b_K | 0, 1). Even if you replace binomial_lpdf with the actual function binomial_lpmf, this refers to a binomial distribution with zero trials, which isnā€™t meaningful.

Is it possible that you donā€™t actually want a binomial distribution, and instead want some other distribution?

For K in the von B model we typically use something like normal(0.5,0.5). Iā€™d also recommend setting up a simulated data set with known parameters first. We use this code for our local minnow population

# Define parameters for the von Bertalanffy growth curve
set.seed(123)  # For reproducibility
N <- 500  # Number of samples
l_inf <- 90  # Asymptotic average length not the real max length
k <- 0.4  # Growth rate coefficient
t0 <- 0  # Hypothetical time at which length is zero
sd_error <- 2  # Standard deviation for additive measurement error

# Simulate time (age)
time <- runif(N, min = 0, max = 10)  # Uniform distribution for time (age) between 0 and 10

# Simulate length using von B growth function with some measurement noise (jitter)
length <- l_inf * (1 - exp(-k * (time - t0))) + rnorm(N, mean = 0, sd = sd_error)

# Combine into a data frame
df <- data.frame(time = time, len = length)

# Plot the simulated data
ggplot(df, aes(x = time, y = len)) +
  geom_point(alpha = 0.5) +
  labs(title = "Simulated Minnow Growth Data",
       x = "Time (Age)",
       y = "Length") +
  theme_minimal()

# Preview the simulated data
head(df)type or paste code here

then try running the model with the simulated data.

1 Like

So this is assuming you are using the simulated data above. The brms code should look like this

#VON BERTALANFFY

#Set formula

formula <- brmsformula(len ~ Linf * (1 - exp(-K * (time - t0))),
                       Linf ~ 1, K ~ 1, t0 ~ 1, nl=TRUE)

#Set priors

priors <- prior(normal(70, 5), nlpar="Linf", lb=0) +
  prior(normal(.6, .1), nlpar="K", lb=0) +
  prior(normal(-1, 0.5), nlpar="t0", lb=0)+
  prior(student_t(3,0,40), class=sigma)

#Fit model
fit1 <- brm(formula,            # calls the formula object created above
            family=gaussian(),  # specifies the error distribution
            data=df,         # the data object
            prior=priors,       # the prior probability object 
            chains=4,           # number of chains, typically 3 to 4
            cores=8,            # number of cores for multi-core processing. Typically set to
            # match number of chains. Adjust as needed to match the number
            # of cores on your computer and number of chains
            iter=2000,          # number of iterations
            warmup=1000,        # number of warm up steps to discard
            control=list(adapt_delta=0.85,    # Adjustments to algorithm to 
                         max_treedepth=15)) # improve convergence.

summary(fit1, prob=0.5)