Is there a good way of hierarchically modeling correlation matrices?

191 views
Skip to first unread message

James Savage

unread,
Nov 26, 2016, 10:23:01 PM11/26/16
to Stan users mailing list
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

unread,
Nov 26, 2016, 10:34:28 PM11/26/16
to stan-...@googlegroups.com
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

unread,
Nov 26, 2016, 11:08:26 PM11/26/16
to Stan users mailing list
Mike - clever. Thanks! 

Ben Goodrich

unread,
Nov 26, 2016, 11:53:48 PM11/26/16
to Stan users mailing list
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

unread,
Nov 27, 2016, 12:30:24 AM11/27/16
to stan-...@googlegroups.com
Stan-users remains the best applied statistics journal out there.

Bob Carpenter

unread,
Nov 27, 2016, 2:39:45 AM11/27/16
to stan-...@googlegroups.com

> 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

unread,
Nov 27, 2016, 9:13:52 AM11/27/16
to Stan users mailing list
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

unread,
Nov 27, 2016, 4:50:08 PM11/27/16
to Stan users mailing list
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

 
Reply all
Reply to author
Forward
0 new messages