Is there a good way of hierarchically modeling correlation matrices?

191 views

James Savage

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

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

Nov 26, 2016, 10:34:28 PM11/26/16
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.

James Savage

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

Ben Goodrich

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

Nov 27, 2016, 12:30:24 AM11/27/16
Stan-users remains the best applied statistics journal out there.

Bob Carpenter

Nov 27, 2016, 2:39:45 AM11/27/16

> 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.

- Bob

James Savage

Nov 27, 2016, 9:13:52 AM11/27/16
to Stan users mailing list

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

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:

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