Error- cannot create vector with constraints made up of elements of a cholesky factor correlation matrix


To constrain hierarchical parameters arising from a bivariate normal distribution ( \beta_{1,i} ,\beta_{2,i}) \sim N( [\mu_{1}, \mu_{2}] , [\sigma_{1}, \sigma_{2}] ) such that \beta_{1,i} \ge \beta_{2,i} we can constrain the parameters such that:

z_{i,1} \ge (\mu_{2} - \mu_{1} + \sigma_{2}*z_{i,2}*\sqrt{1-\rho}) / (\sigma_{1} - \rho*\sigma_{2})

using off-centred parameterisations.

I am trying to do this in Stan by declaring:

parameters {
              ordered[2] a1_m_raw[num_refs]; 
              vector<lower=0>[2] sd[nt]; 
              cholesky_factor_corr[2] L_Omega[nt];
              vector<lower=(to_vector(a1_m_raw[ , 1]) - to_vector(a1_m_raw[ , 2]) + sd[1,2]*z2*L_Omega[1,2,2])/(sd[1,1] - L_Omega[1,1,2]*sd[1,2])>[NS] z1; //  Se + Fp > 1

However, I am getting complication errors when attempting to compile the code. It is very long error message but I will copy-paste a screenshot of part of the code below:

Would anyone have any idea why this is happening? When I declare a similar constraint for 2 univariate normals (i.e. no L_Omega’s needed in the constraint) it works fine


1 Like

I get similarly confusing errors on my computer. I made an issue:

Even if this isn’t possible (not clear why it would be impossible tho), we should be able to give a better error message.

Edit: Thanks for reporting it!

1 Like

Thanks for making the issue!

Regardless of whether it should work or not (in terms of programming), perhaps there is a more sensible way of declaring a bivariate normal with the constraints mentioned above in Stan?

You’ll probably have to code these constraints manually unfortunately.

Check the docs here and here.

Here’s an example of doing manual constraints on vectors: How to specify variable parameter constraints based on other model parameters

1 Like

Thanks for these links!

So for my example I think it would be this:

parameters {
              ordered[2] a1_m_raw[num_refs];
              vector<lower=0>[2] sd_ref[num_refs];
              cholesky_factor_corr[2] L_Omega_ref[num_refs];
              vector[n_studies] z2_unconstrained; 
              vector[n_studies] z1_unconstrained; 

transformed parameters {
  vector[n_studies] z2; 
  vector[n_studies] z1; 

             for (s in 1:n_studies) {  
                      z2[s] = exp(z2_unconstrained[s]);
                      z1[s] = exp(z1_unconstrained[s]) + (a1_m_raw[ref[s] , 1] - a1_m_raw[ref[s], 2] + sd_ref[1,2]*z2[s]*L_Omega_ref[1,2,2])/(sd_ref[1,1] - L_Omega_ref[1,1,2]*sd_ref[1,2]) ; //  Se + Fp > 1 

model {

       z1 ~ std_normal(); 
       z2 ~ std_normal(); 

   for (s in 1:n_studies) {
           target += z1_unconstrained[s];
           target += z2_unconstrained[s];


Is it OK to set the priors on these transformed parameters (z1 and z2) in this way?

1 Like

Yup, that’s what the doing the Jacobian transform allows you to do. You’ll get warnings about putting priors on transformed parameters, but you’ve added your Jacobian adjustment so assuming it’s bug-free you’re good to go.

For z2 I guess you can just use Stan’s built-in constrain stuff, but maybe that was just a test that you were coding it correctly.

1 Like

Thanks - z2 has no constraint so I guess I can just declare it normally. I don’t get any Jacobian warnings strangely

1 Like