Non-centered parameterization slows down the model

[scijens]

Hi everyone,

I am running a Hurdle model with a hierarchical structure and I was having issues with poor convergence. Based on a helpful suggestion in this forum, I applied the non-centered parameterization (NCP) on the individual-specific intercepts, which helped to make the model converge. However, the model with NCP runs 3 times as long. Is there a way to speed the model with NCP up?

The Hurdle model is integrated in a hidden Markov model, which already makes the model computationally demanding. Here are the relevant parts of my model without NCP.

data{
  int<lower = 1> N; // number of observations
  int<lower = 0> S; // number of states
  int<lower = 1> H; // number of individuals
  int<lower = 1> id[N]; // identifier of individuals
  int<lower = 0> K1; // number of covariates inside delta in Bernoulli part
  int<lower = 0> K2; // number of covariates inside delta in Lognormal part
  matrix[N, K1] C1; // matrix of covariates for Bernoulli part
  matrix[N, K2] C2; // matrix of covariates for Lognormal part
  int<lower = 0, upper = 1> y[N]; // binary dependent variable
  real q[N]; // continuous variable we are modeling conditional on y
}

parameters {
  ordered[S] mu; // state-dependent intercepts in Bernoulli part
  vector[S] nu; // state-dependent intercepts in Lognormal part

  real alphaj[H]; // individual-specific intercept in Bernoulli part
  real alphai[H]; // individual-specific intercept in Lognormal part
  real<lower = 0> sigma_alphaj; // sd of individual-specific intercept in Bernoulli part
  real<lower = 0> sigma_alphai; // sd of individual-specific intercept in Lognormal part

  real<lower = 0> sigma_q;
  vector[K1] delta1;
  vector[K2] delta2;
}
model {
  // priors
  mu[1] ~ normal(0, 1)T[, mu[2]];
  mu[2] ~ normal(0, 2);

  nu[1] ~ normal(0, 1);
  nu[2] ~ normal(0, 2);

  alphai ~ normal(0, sigma_alphai);
  alphaj ~ normal(0, sigma_alphaj);
...
  for (t in 1:N) {
      target += log_sum_exp(gamma);
      for (k in 1:S){
        gamma_prev[k] =  bernoulli_logit_lpmf(y[t] | alphaj[id[t]] + mu[k] + C1[t]*delta1);
        if(y[t] == 1){
	  gamma_prev[k] += lognormal_lpdf(q[t] | alphai[id[t]] + nu[k] + C2[t]*delta2, sigma_q);
...

I just copied parts of the output summary, showing the problems with convergence indicated by Rhat and low ESS. I ran the model with 1000 warmup and 2000 iterations, treedepth at 11 and adapt_delta at .90.

The updated model with NCP looks like this:

...
parameters {
...
  real alpha0i;
  real alpha0j;
  vector[H] alphai_raw;
  vector[H] alphaj_raw;
  real<lower=0> sigma_alphai;
  real<lower=0> sigma_alphaj;
  real<lower = 0> sigma_q;
}

transformed parameters {
  vector[H] alphai;
  vector[H] alphaj;
    alphai = alpha0i + sigma_alphai * alphai_raw;
    alphaj = alpha0j + sigma_alphaj * alphaj_raw;
}

model {
...
  alphai_raw ~ std_normal();
  alphaj_raw ~ std_normal();

  alphai ~ normal(0, sigma_alphai);
  alphaj ~ normal(0, sigma_alphai);

  alpha0i ~ normal(0, 1);
  alpha0j ~ normal(0, 1);

I ran this updated model with NCP on the same data with same number of iterations. Here is the output:

Is there anything I can do to speed up the model with NCP?

Thanks in advance!

[andrjohns]

Hi Jens, the two parameterisations that you’ve given are implying different distributions.

In the first, you have:

  alphai ~ normal(0, sigma_alphai);
  alphaj ~ normal(0, sigma_alphaj);

Which is:

\alpha_i \sim N(0,\sigma_{\alpha_i}) \\ \alpha_j \sim N(0,\sigma_{\alpha_j})

And in the second you have:

    alphai = alpha0i + sigma_alphai * alphai_raw;
    alphaj = alpha0j + sigma_alphaj * alphaj_raw;

Which is:

\alpha_i \sim N(\alpha_{0i},\sigma_{\alpha_i}) \\ \alpha_j \sim N(\alpha_{0j},\sigma_{\alpha_j})

Additionally, you still have the original prior on alpha_i and alpha_j in your non-centered model:

  alphai ~ normal(0, sigma_alphai);
  alphaj ~ normal(0, sigma_alphai);

Which means that you’re effectively giving each parameter two different priors.

The simplest way to implement the non-centered parameterisation in your model is to use the offset and multiplier syntax when you declare the parameters and leave the rest of your model unchanged from the original:

parameters {
  ordered[S] mu; // state-dependent intercepts in Bernoulli part
  vector[S] nu; // state-dependent intercepts in Lognormal part

  real<lower = 0> sigma_alphaj; // sd of individual-specific intercept in Bernoulli part
  real<lower = 0> sigma_alphai; // sd of individual-specific intercept in Lognormal part
  vector<offset=0, multiplier=sigma_alphaj>[H] alphaj; // individual-specific intercept in Bernoulli part
  vector<offset=0, multiplier=sigma_alphai>[H] alphai; // individual-specific intercept in Lognormal part

  real<lower = 0> sigma_q;
  vector[K1] delta1;
  vector[K2] delta2;
}

Note that the offset=0 is redundant and unnecessary, but I’ve included it for clarity.

[scijens]

Thank you so much, Andrew.

This makes sense and I hope it’s going to run faster with this update.

Just to be sure: Using the parameterization you are proposing, I keep the priors on alphai and alphaj in the mode block and I do not add any additional priors?

  alphai ~ normal(0, sigma_alphai);
  alphaj ~ normal(0, sigma_alphaj);

I have to get used to this way of implementing the non-centered parameterization, but it seems less error-prone to me.

[andrjohns]

Yep that’s right, you just need to change the variable declarations in the parameters block. Everything else stays the same

[scijens]

It’s throwing that Identifier 'sigma_alphai' not in scope.. Any ideas what I can do against it?

[andrjohns]

Ah I should have clarified, you need to make sure that sigma_alpha is declared before including it in multiplier. The parameters block I posted above has this already done

[scijens]

This was my mistake, sorry, it started sampling now.