Is there a good way of hierarchically modeling correlation matrices?

[James Savage]

Hi all, 

Could anyone recommend some reading on hierarchically modeling correlation 

matrices? 

My problem looks like this. I have a multi-dimensional unbounded outcome Y. 

You could think of each observation as being a vector of GDP growth, investment 

growth and consumption growth. I have repeated observations for many countries. 

If I had a single country, I would just build a vector autoregression, which Stan 

handles nicely. Seeing as I have many countries, I want to build a hierarchical 

VAR. Something like

Y[n] ~ multi_normal(mu[country[n]], quad_form_diag(Omega[country[n]], tau[country[n]])) 

but where I have a model for each country's mu. 

The question is, I can easily specify hierarchical priors for mu and tau, but I don't 

know what I should be using for Omega. At the moment I'm modeling

Omega[country] ~ lkj(nu) 

nu ~ half cauchy(0, 5)

But this gives away mountains of information; might as well not be doing partial pooling. 

Any thoughts on approaches to take? 

Jim

[Michael Betancourt]

The most natural thing to do is a hierarchical model

on the latent, unconstrained space.  For a correlation

matrix that reduces to the N(N-1)/2 lower diagonal

elements of the Cholesky factor.  In other words,

have N(N-1)/2 independent univariate hierarchical

models then fill in a the lower diagonal elements

with those values, ones on the diagonal, and then

construct Omega.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[James Savage]

Mike - clever. Thanks! 

[Ben Goodrich]

On Saturday, November 26, 2016 at 10:34:28 PM UTC-5, Michael Betancourt wrote:

The most natural thing to do is a hierarchical model

on the latent, unconstrained space.  For a correlation

matrix that reduces to the N(N-1)/2 lower diagonal

elements of the Cholesky factor.  In other words,

have N(N-1)/2 independent univariate hierarchical

models then fill in a the lower diagonal elements

with those values, ones on the diagonal, and then

construct Omega.

This would be easier if the transformation from

real numbers -> canonical partial correlations -> Cholesky factors -> correlation matrices

were exposed in the Stan language. You could look at the implementations (yes, there are two for reasons I don't understand) in Stan or you could access them yourself with a little C++ in version 2.13.

But perhaps a more fundamental issue is that I don't think anyone knows what hierarchical prior to put on the real numbers in the unconstrained space. The LKJ thing is derived from a symmetric beta priors over the (-1,1) interval on the canonical partial correlations, but those have mean zero and so their derivation wouldn't go through in a hierarchical setting. The onion method is more promising since it is derived from unit vectors to the left of the diagonal on the Cholesky factor of a KxK correlation matrix, but there is no easy way in the Stan language to get as many unit vectors as there are countries for each value of k between 2 and K.

An alternative would be to set in transformed parameters

for (c in 1:N_countries) Omega[c] = alpha * Omega_global + (1 - alpha) * Omega_local[c];

where Omega_global is a correlation matrix with a LKJ(1) (i.e. joint uniform) prior, each Omega_local[c] is a correlation matrix with a LKJ(eta) prior for eta >> 1, and alpha is a weight over the [0,1] interval with a beta or Kumaraswamy prior. Or maybe set eta = alpha / (1 - alpha).

Ben

[Andrew Gelman]

Stan-users remains the best applied statistics journal out there.

[Bob Carpenter]

> On Nov 26, 2016, at 11:53 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Saturday, November 26, 2016 at 10:34:28 PM UTC-5, Michael Betancourt wrote:
> The most natural thing to do is a hierarchical model
> on the latent, unconstrained space. For a correlation
> matrix that reduces to the N(N-1)/2 lower diagonal
> elements of the Cholesky factor. In other words,
> have N(N-1)/2 independent univariate hierarchical
> models then fill in a the lower diagonal elements
> with those values, ones on the diagonal, and then
> construct Omega.

> This would be easier if the transformation from
>
> real numbers -> canonical partial correlations -> Cholesky factors -> correlation matrices
>
> were exposed in the Stan language.

Absolutely. That's been on the to-do list for a while. This one wouldn't
even need us to expose the Jacobian calculation.

> You could look at the implementations (yes, there are two for reasons I don't understand) in Stan or you could access them yourself with a little C++ in version 2.13.

If you're talking about two functions to do the constraining
transforms, it's because one computes Jacobians and one doesn't.

If it's deeper than that, it may just be a missed opportunity
to consolidate code.

> But perhaps a more fundamental issue is that I don't think anyone knows what hierarchical prior to put on the real numbers in the unconstrained space.

:-) I was just about to ask about that.

- Bob

[James Savage]

Ben, this is very helpful. 

A question about the partial pooling weights. In this specification, alpha is common across countries. My panel is very imbalanced, so I don't know that is the best assumption (it might be, I'm just not sure). Are you aware of any methods that would allow me to set/estimate these weights, taking into account the relative noisiness of the local and global correlations? 

Or is that implicitly set by the relative degrees of freedom of the local correlation prior? 

[Ben Goodrich]

On Sunday, November 27, 2016 at 9:13:52 AM UTC-5, James Savage wrote:

Ben, this is very helpful. 

A question about the partial pooling weights. In this specification, alpha is common across countries. My panel is very imbalanced, so I don't know that is the best assumption (it might be, I'm just not sure). Are you aware of any methods that would allow me to set/estimate these weights, taking into account the relative noisiness of the local and global correlations? 

Or is that implicitly set by the relative degrees of freedom of the local correlation prior? 

You could allow alpha and / or eta to vary by country. Certainly, if you have fewer observations on some countries than others, then the posterior distribution of the parameters for those countries will tend to be less concentrated. But I don't really see how that would affect your prior beliefs about alpha and / or eta prior to observing the data (and thus which countries have more observations).

Ben