Modeling heteroscedasticity (sigma) in a nonlinear hierarchical phenomenon

Modeling
, ,
2

I have found out through the forum with @paul.buerkner that non-linear formula except for the main one need to be wrapped in nlf() to work. But I still having trouble declaring a prior for the tau parameter. Perhaps due to lack of experience. Something is missing.

I noticed that the script below works when I declare the priors through the get_prior() method. So the error is precisely in the way I am declaring it, I believe.

And then I improved the pseudocode to:

rm(list = ls())
# loading the data
load(file="bio.RData")

# Preamble
library("brms")
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)


# Initial values
init <- function(chain_id=1) {
  list ( alpha  =  rep(80,n_cultivo) ,
         kappa  =  rep(0.043,n_cultivo) ,
         delta  =  rep(1.8,n_cultivo),
         tau  =  rep(0.055,n_cultivo))
}
n_cultivo = length(unique(bio$cultivo))

# NUTS configuration
control. <- list(
  adapt_engaged = TRUE,
  adapt_delta = 0.90, 
  stepsize = 0.05, 
  max_treedepth = 10
)

#### modelos ####
# Piors
prior. <- c(
  set_prior("normal(80, 5)", nlpar = "alpha",lb=0),
  set_prior("normal(0.043,0.02)", nlpar = "kappa", lb=0),
  set_prior("normal(1.8,0.15)", nlpar = "delta", lb=0),
  set_prior("gamma(1,1/.1)", nlpar="tau", class = "b", lb=0, dpar = "sigma")
)

# Formula  bf() = brmsformula()
formula. <- bf(
  # location parameter
  peso ~ alpha - (alpha/(1+(time*kappa)^delta)), # f(x)
  # scale parameter
  nlf(sigma ~ time * tau * (alpha*delta*kappa*(kappa*time)^(delta-1) )/( ((kappa*time)^delta)^2+2*(kappa*time)^delta+1) ), # time*tau*f'(x)
  # Nonlinear variables
  alpha + kappa + delta + tau ~ (1 | tq)+(1 | cultivo:tq),
  # Nonlinear fit
  nl = TRUE)

prior. <- get_prior(
  formula = formula.,
  data = bio,
  family=gaussian()
)

# fitting a model
fit <- brm(
  formula.,
  family=gaussian(), 
  data = bio,
  prior = prior.,
  control = control.,
  init=init,
  chains = 4,
  cores = getOption("mc.cores", 1)
)

It must be something so simple that I didn’t find.