How to specify non centered parameterization for sigma in a distributional model in brms

Modeling
,
1

I want to fit a distributional model with brms along the following lines:

df <- data.frame() # some columns include:
  # value (continuous, zero-mean),
  # afactor (3 levels), bfactor (7 levels), 
  # gfactor (240 levels), hfactor (2 levels)
frml <- bf(
  value ~ 0,
  sigma ~ 0 + Intercept + afactor * bfactor + (1 | gfactor) + (1 | hfactor)
)
priors <- c(
  set_prior("normal(0,0.5)", class = "sd", dpar = "sigma"),
  set_prior("normal(0,0.5)", class = "b", dpar = "sigma"),
  set_prior("normal(3,1)", class = "b", coef = "Intercept", dpar = "sigma")
)
bfit <- brm(
  frml,
  data = df,
  family = gaussian(),
  prior = priors
)

I’m getting a large number of divergent transitions (consistently more than 100, sometimes more than 200). I’ve already tried to increase adapt_delta, but that made convergence (number of divergent transitions, Rhat, ESS) even worse.

Since sigma is the value being estimated by predictors here, I realized that this may be “The Devil’s Funnel” (I’ve been following Statistical Rethinking, 13.4.1), so I should use non-centered parameterization when estimating sigma from my predictors.

How can I specify that in brms, or how can I edit the stan code generated by brms? Alternatively, is this something I should achieve by manipulating {value, afactor, bfactor, gfactor, hfactor} in my data frame, and if so, what should I be doing to those columns?

2

hfactor (2 levels) might be your problem. In my experience, brm() usually makes a lot of divergent transitions when you fit random effects for factors with fewer than, say, 10 levels. Consider making hfactor a simple fixed effect instead.

3

Thanks, @Solomon, I’ll give letting hfactor be a fixed effect a try.

As a broader question about understanding brms, stan, and Bayesian modeling though, if I really wanted to use non-centered parameterization on sigma in a distributional model, how would I go about that?

4

I’m under the impression brms always uses the non-centered parameterization, even for distributional parameters. Have you run into documentation that states otherwise?

5

Not from documentation, but from inspecting the stan code generated by brms, it looks like it’s using non-centered parameterization for the group-level random effects but not the population-level fixed effects.

Here is the full stan code for my model:

// generated with brms 2.23.0
functions {
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> K_sigma;  // number of population-level effects
  matrix[N, K_sigma] X_sigma;  // population-level design matrix
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  array[N] int<lower=1> J_1;  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_sigma_1;
  // data for group-level effects of ID 2
  int<lower=1> N_2;  // number of grouping levels
  int<lower=1> M_2;  // number of coefficients per level
  array[N] int<lower=1> J_2;  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_2_sigma_1;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  vector[K_sigma] b_sigma;  // regression coefficients
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  array[M_1] vector[N_1] z_1;  // standardized group-level effects
  vector<lower=0>[M_2] sd_2;  // group-level standard deviations
  array[M_2] vector[N_2] z_2;  // standardized group-level effects
}
transformed parameters {
  vector[N_1] r_1_sigma_1;  // actual group-level effects
  vector[N_2] r_2_sigma_1;  // actual group-level effects
  // prior contributions to the log posterior
  real lprior = 0;
  r_1_sigma_1 = (sd_1[1] * (z_1[1]));
  r_2_sigma_1 = (sd_2[1] * (z_2[1]));
  lprior += normal_lpdf(b_sigma[1] | 3,1);
  lprior += normal_lpdf(b_sigma[2] | 0,0.5);
  lprior += normal_lpdf(b_sigma[3] | 0,0.5);
  lprior += normal_lpdf(b_sigma[4] | 0,0.5);
  lprior += normal_lpdf(b_sigma[5] | 0,0.5);
  lprior += normal_lpdf(b_sigma[6] | 0,0.5);
  lprior += normal_lpdf(b_sigma[7] | 0,0.5);
  lprior += normal_lpdf(b_sigma[8] | 0,0.5);
  lprior += normal_lpdf(b_sigma[9] | 0,0.5);
  lprior += normal_lpdf(b_sigma[10] | 0,0.5);
  lprior += normal_lpdf(b_sigma[11] | 0,0.5);
  lprior += normal_lpdf(b_sigma[12] | 0,0.5);
  lprior += normal_lpdf(b_sigma[13] | 0,0.5);
  lprior += normal_lpdf(b_sigma[14] | 0,0.5);
  lprior += normal_lpdf(b_sigma[15] | 0,0.5);
  lprior += normal_lpdf(b_sigma[16] | 0,0.5);
  lprior += normal_lpdf(b_sigma[17] | 0,0.5);
  lprior += normal_lpdf(b_sigma[18] | 0,0.5);
  lprior += normal_lpdf(b_sigma[19] | 0,0.5);
  lprior += normal_lpdf(b_sigma[20] | 0,0.5);
  lprior += normal_lpdf(b_sigma[21] | 0,0.5);
  lprior += normal_lpdf(sd_1 | 0,0.5)
    - 1 * normal_lccdf(0 | 0,0.5);
  lprior += normal_lpdf(sd_2 | 0,0.5)
    - 1 * normal_lccdf(0 | 0,0.5);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    // initialize linear predictor term
    vector[N] sigma = rep_vector(0.0, N);
    sigma += X_sigma * b_sigma;
    for (n in 1:N) {
      // add more terms to the linear predictor
      sigma[n] += r_1_sigma_1[J_1[n]] * Z_1_sigma_1[n] + r_2_sigma_1[J_2[n]] * Z_2_sigma_1[n];
    }
    sigma = exp(sigma);
    target += normal_lpdf(Y | mu, sigma);
  }
  // priors including constants
  target += lprior;
  target += std_normal_lpdf(z_1[1]);
  target += std_normal_lpdf(z_2[1]);
}
generated quantities {
}

At the beginning of the transformed parameters block, we see the following two lines “undoing” the non-centered paramterization for the group-level effects but not for the other effects:

  vector[N_1] r_1_sigma_1;  // actual group-level effects
  vector[N_2] r_2_sigma_1;  // actual group-level effects

My concern is that because the value being predicted is a distributional parameter, the Devil’s funnel may be relevant for all predictors, including the population-level effects.

6

I believe the “non-centered parameterization” is only meant to refer to random effects in brms, and I believe that’s also generally how McElreath used the term.

7

You might be equating the centered parameterization for random effects with mean centering a variable. If so, those are two very different things.

8

I believe the “non-centered parameterization” is only meant to refer to random effects in brms, and I believe that’s also generally how McElreath used the term.

I see. I was hoping there would be a way to reparameterize the model to deal with a potential effect of the Devil’s Funnel in the distributional model setting (even without random effects), but maybe I’m misunderstanding where that problem comes from.

9

Yes, I think you might be misunderstanding things. brms is already using the non-centered parameterization to avoid the devil’s funnel issue McElreath discussed in his text. There’s a lot of crossover with brms and rethinking, but they have a few very important differences. Another related issue is I believe adapt_delta defaults to 0.95 in rethinking, but to 0.8 in brms.

Based on what you’ve said about your data, using random effects for a factor with 2 levels is the first place I’d look to for complications. This would also be a problem with a frequentist function like lmer(), by the way. But anyway, the solution isn’t the centered vs non-centered parameterization issue. The solution is to rethink whether it’s worth the headache to model a standard deviation with two data points. Technically you can, but I wouldn’t recommend it if you have a convenient alternative.