# Modeling heteroscedasticity (sigma) in a nonlinear hierarchical phenomenon

Hello guys,
It’s a pleasure to be here again.

I am modeling the heteroscedasticity of a nonlinear hierarchical Gaussian (sigmoidal) growth model f(x). Here’s the model:

w_{ij} \sim N(f(x,\theta_{ij}),\sigma_{ij})
f(x,\theta_{ij}) = nonlinear function
\sigma_{ij} = \tau_{ij} \cdot x \cdot f'(x,\theta_{ij})

where i is the index of the first level of the group and j second level of the group of a nested hierarchical design, w_{ij} is the weight, N() truncated normal distribution (positive), \theta_{ij} set of parameters governing the nonlinear model f(x,\theta_{ij}), \sigma_{ij} scale parameter in which I am trying to model the variability of the dependent variable w_{ij}, \tau_{ij} a weighting parameter, x = time the dependent variable and f'(x,\theta_{ij}) the first derivative of the function f(x,\theta_{ij}) with respect to x.

I’m working with the rstan interface with the brms package. The pseudo-code and simulated data are below and attached to the topic, including a reproducible script. It’s not working. Very likely the error is in the way I am declaring sigma (scale parameter) and in declaring the a priori probability density.

My question is:

1. Is it possible to declare sigma with this form in stan (brms package)?
2. If yes, how can I do it?

Thanks in advance for reading this topic. I hope I was clear enough in the model declaration and in the problem doubt.


rm(list = ls())

# 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)", nplar="tau", lb=0, dpar = "sigma") ) # Formula formula. <- bf( # location parameter peso ~ alpha - (alpha/(1+(time*kappa)^delta)), # f(x) # scale parameter sigma ~ time * tau * (alpha*delta*kappa*(kappa*x)^(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) # fitting a model fit <- brm( formula., family=gaussian(), data = bio, prior = prior., control = control., init=init, chains = 4, cores = getOption("mc.cores", 1) )  bio.RData (4.9 KB) 2_Script.R (1.5 KB) 1 Like 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(
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.

I believe I solved it. Simply wrapping the nonlinear expression of sigma with the nlf() function solved the problem.
I need to review how to declare priors for brms and this is easy to get from the stan and brms package manual and guide. Thank you, guys. I’m available for whatever you need to.

1 Like