Modelling sigma

Modeling
6

You could model the mean and the standard deviation of the lognormal distribution with some reparametrization. To do that in brms you have to specify a custom family. I submitted an issue to the brms github repo a while back doing exactly that. https://github.com/paul-buerkner/brms/issues/1027

Here is the code, which is slightly rewritten to take into account that one might want to model the standard deviation as well:

custom_family(name = "lognormal7", dpars = c("mu", "sigma"), 
              links = c("log", "log"), 
              lb = c(0, 0), type = "real") ->
  lognormal7


stan_funs <- "
  real lognormal7_lpdf(real y, real mu, real sigma) {
    return lognormal_lpdf(y | log(mu/sqrt(1+sigma^2/mu^2)), sqrt(log(1+sigma^2/mu^2)));
  }
  real lognormal7_rng(real mu, real sigma) {
    return lognormal_rng(log(mu/sqrt(1+sigma^2/mu^2)), sqrt(log(1+sigma^2/mu^2)));
  }
"


stanvars <- stanvar(scode = stan_funs, block = "functions")

log_lik_lognormal7 <- function(i, prep) {
  mu <- prep$dpars$mu[, i]
  if(NCOL(prep$dpars$sigma)==1){sigma <- prep$dpars$sigma}else
  {sigma <- prep$dpars$sigma[, i]}   ## [, i] if sigma is modelled, without otherwise
  y <- prep$data$Y[i]
  lognormal7_lpdf(y, mu, sigma)
}


posterior_predict_lognormal7 <- function(i, prep, ...) {
  mu <- prep$dpars$mu[, i]
  if(NCOL(prep$dpars$sigma)==1){sigma <- prep$dpars$sigma}else
  {sigma <- prep$dpars$sigma[, i]}   ## [, i] if sigma is modelled, without otherwise
  lognormal7_rng(mu, sigma)
}

posterior_epred_lognormal7 <- function(prep) {
  mu <- prep$dpars$mu
  return(mu)
}

and how to use it:

brm(data = mtcars, formula = bf(mpg ~ wt), 
    family = lognormal7, stanvars = stanvars) -> model_test1

brm(data = mtcars, formula = bf(mpg ~ wt, sigma ~ wt), 
    family = lognormal7, stanvars = stanvars) -> model_test2

expose_functions(model_test1, vectorize = TRUE) # exposes the Stan functions to R
model_test1 %>% conditional_effects()
model_test2 %>% conditional_effects()

model_test1 %>% conditional_effects(method = "posterior_predict")
model_test2 %>% conditional_effects(method = "posterior_predict")

loo(model_test1, model_test2)
1 Like