Modeling heteroscedasticity (sigma) in a nonlinear hierarchical phenomenon

Modeling
, ,
1

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())
# 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 (4.9 KB)
2_Script.R (1.5 KB)

1 Like
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.

3

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