# Modeling heteroscedasticity (sigma) in a nonlinear hierarchical phenomenon
**URL:** https://discourse.mc-stan.org/t/modeling-heteroscedasticity-sigma-in-a-nonlinear-hierarchical-phenomenon/27877
**Category:** Modeling
**Tags:** specification, fitting-issues, brms
**Created:** [June 18, 2022, 2:32am UTC](https://discourse.mc-stan.org/t/modeling-heteroscedasticity-sigma-in-a-nonlinear-hierarchical-phenomenon/27877 "2022-06-18T02:32:48Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@Carlos\_Zarzar](https://discourse.mc-stan.org/u/Carlos_Zarzar)
#### Post date: [June 18, 2022, 2:32am UTC](https://discourse.mc-stan.org/t/modeling-heteroscedasticity-sigma-in-a-nonlinear-hierarchical-phenomenon/27877/1 "2022-06-18T02:32:48Z")
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.
```stan
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)", 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](https://discourse.mc-stan.org/uploads/short-url/2nCO1E776iF0HfxdWRErNiI0lNV.RData) (4.9 KB)
[2\_Script.R](https://discourse.mc-stan.org/uploads/short-url/kEui1fXbeSs0RQDsjAEh2cGDuGF.R) (1.5 KB)
---
### Author:  [@Carlos\_Zarzar](https://discourse.mc-stan.org/u/Carlos_Zarzar)
#### Post date: [June 19, 2022, 2:13pm UTC](https://discourse.mc-stan.org/t/modeling-heteroscedasticity-sigma-in-a-nonlinear-hierarchical-phenomenon/27877/2 "2022-06-19T14:13:13Z")
I have found out through [the forum](https://discourse.mc-stan.org/t/non-linear-specification-variable-is-not-a-valid-distributional-or-non-linear-parameter/23521/3) 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:
```nohighlight
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.
---
### Author:  [@Carlos\_Zarzar](https://discourse.mc-stan.org/u/Carlos_Zarzar)
#### Post date: [June 20, 2022, 10:35am UTC](https://discourse.mc-stan.org/t/modeling-heteroscedasticity-sigma-in-a-nonlinear-hierarchical-phenomenon/27877/3 "2022-06-20T10:35:06Z")
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.