Guide to using the onion method

Modeling
1

Hello,

I am trying to implement the “onion method” for the construction of the var-cov matrix. (I think some of the problems I was running into earlier were due to initialising a 48x48 matrix).

I have found the code for the onion method here: Using the onion method in occupancy hierarchical model

And I would like some tips on using it when specifying a multi-level linear model (for example, following the code at Bayesian multilevel models using R and Stan (part 1) | Matteo Lisi)

The onion method generates a L[K, K] and Sigma[K, K]. Is the Sigma equivalent to the L_u parameter used in Lisi’s code?

Thanks

parameters {
  vector[2] beta;                   // fixed-effects parameters
  real<lower=0> sigma_e;            // residual std
  vector<lower=0>[2] sigma_u;       // random effects standard deviations

  // declare L_u to be the Choleski factor of a 2x2 correlation matrix
  cholesky_factor_corr[2] L_u;

  matrix[2,J] z_u;                  // random effect matrix
}

transformed parameters {
  // this transform random effects so that they have the correlation
  // matrix specified by the correlation matrix above
  matrix[2,J] u;
  u = diag_pre_multiply(sigma_u, L_u) * z_u;

}

model {
  real mu; // conditional mean of the dependent variable

  //priors
  L_u ~ lkj_corr_cholesky(1.5); // LKJ prior for the correlation matrix
  to_vector(z_u) ~ normal(0,2);
  sigma_e ~ normal(0, 5);       // prior for residual standard deviation
  beta[1] ~ normal(0.3, 0.5);   // prior for fixed-effect intercept
  beta[2] ~ normal(0.2, 2);     // prior for fixed-effect slope

  //likelihood
  for (i in 1:N){
    mu = beta[1] + u[1,Subject[i]] + (beta[2] + u[2,Subject[i]])*Days[i];
    RT[i] ~ normal(mu, sigma_e);
  }
}
2

I don’t know the answer to this.

As for general tips, here’s some questions:

  1. Is your code currently running or not running?
  2. What problems are you having with your code?

Lastly, here’s an example for the working onion method UMESC / quant-ecology / Occstanhm · GitLab (usgs.gov).

3

Hi,

my code is running, but not well. So I am backtracking to work up a simple multi-level linear model example that I can share and debug. (my main aim is a multi-level mixture model, which is getting quite complex!)

I was hoping that somebody would be able to help give some context as to what some of the parameters are, and how they related to the parameters in the “standard” (non-onion?) way of doing multi-level models in Stan.

I’ll try and get the simplified example code done later this week.

4

As a tip, please post a working example (ideally with simulated date) here. That will help others see what’s going on.

Some guesses include (none tested by me, only debugging tips):

  • You may not have enough data.
  • Your model may have identifiability issues.
  • You may have a coding problem in your model.

Simulating data will help you with all three.

5

My issue just now is none of the above… it is a more basic conceptual question about how I should replace the standard (or at least, more common?) method for specifying multi-level models with the methods used in the onion.stan example shared in the linked thread.

I know full well that there are identifiability and coding issues in my model, as I don’t understand what I am doing :)

So far, my research has turned up the simple onion.stan code that generates a var-cov matrix, and your far more complete and complex model. I was simply hoping for some pointers about how to use the onion method in a simpler example.

As it seems like this doesn’t yet exist, I will have a go at coding up myself and get back to you.

6

Sounds like a good plan.

And congratulations, you’ve reached the edge of Stan documentation and features! You’re becoming an expert.

7

Is there a reason you need to do this manually rather than using Stan’s built-in covariance parameterization, which is based on an onion method? @bgoodri is the Stan developer who knows the most about these things, I think—he’s the one who coded our constrained types and introduced the LKJ prior from the Lewandowski, Kurowicka, and Joe paper on the onion and vine methods.

1 Like
8

Maybe there is updated code/functions I don’t know about?

I was originally using something very similar to the code below, but was having real problems with initialisation when the number of dimensions increased (if I do a full maximal random effect structure, I need a 48x48 matrix!)

parameters {
  vector[2] beta;                   // fixed-effects parameters
  real<lower=0> sigma_e;            // residual std
  vector<lower=0>[2] sigma_u;       // random effects standard deviations

  // declare L_u to be the Choleski factor of a 2x2 correlation matrix
  cholesky_factor_corr[2] L_u;

  matrix[2,J] z_u;                  // random effect matrix
}

transformed parameters {
  // this transform random effects so that they have the correlation
  // matrix specified by the correlation matrix above
  matrix[2,J] u;
  u = diag_pre_multiply(sigma_u, L_u) * z_u;

}