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.
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)
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)