Fwd: universal prior for reg coefs; first step

[Ben Goodrich]

---------- Forwarded message ----------

From: Andrew Gelman <gel...@stat.columbia.edu>
Date: Wed, Dec 16, 2015 at 2:33 PM
Subject: universal prior for reg coefs; first step

The first first step was undertaken in 2010 when we twice discussed the implications of this prior for mi (see slides).

• For an LKJ distribution on the correlation matrix between X and Y, the distribution of the standardized coefficients is multivariate Pearson type II, which is a special case of an elliptical distribution with a beta generating function
  
• Each marginal standardized coefficient has an unnamed density that can be written in "closed" form with a Gaussian hypergeometric function
• The mean exists always, the variance exists if eta > 0.5, the skewness exists if eta > 1, the kurtosis exists if eta > 1.5, etc. Aside from possibly having heavy tails, it does not have a similar shape as a Cauchy density

Many, many steps have been taken since then. I've taught this prior for the last two years in the QMSS class. For a long time, stan_lm in rstanarm was doing this prior. The correlation matrix between X and Y implies a distribution --- that is beta(K/2, eta) --- for the R^2 so prior information about the location (mode, mean, or median) of the R^2 implies a value for eta. There are between one and two problems with this approach:

1. If you have a distribution of standardized coefficients, you have to unstandardize them to use in the likelihood. For the X variables, it is straightforward to divide the standardized coefficients by their standard deviation but you need to multiply by the marginal standard deviation of Y. I prefer making the marginal standard deviation of Y a primitive that has a Jeffreys prior and you prefer anything but that.
2. To get the standardized coefficients, you have to pre-multiply a vector by the inverse of a KxK matrix, which is slow in Stan for large K. So, stan_lm() does a QR-based reparameterization that avoids the need to do O(K^3) matrix operations. But now the prior is only over the covariance / correlation between the columns of Q and Y, rather than a whole correlation matrix. This is explained in the vignette: https://github.com/stan-dev/rstanarm/blob/master/vignettes/lm.Rmd
   

The main obstacles to generalizing these ideas to other distributional forms are:

1. There are currently no other distributions for a correlation matrix with known properties.
2. All known alternative ways of constructing a correlation matrix from unconstrained values besides the Stan variant of LKJ require inversion of submatrices. So whatever their properties are, they would need to be worth an O(K^4) transformation from the unconstrained space.
   

We do have some distributions for covariance matrices that can be rescaled into correlation matrices. I tend to not like the shapes of Wishart or inverse Wishart, but in this context you have the additional problem that there is no information in the data that informs the variances. This makes trouble for HMC but also it is hard to think about what the prior scatter matrix should be. 

 

To: Ben Goodrich <goodri...@gmail.com>, Michael Betancourt <betan...@gmail.com>


Hi, see attached.  I simulated corr matrices from LKJ distribution, Wishart, and inv-Wishart.  (for the latter two options, we get cov matrices so I rescaled to get corr matrices).  For each, I computed all the regression coefficients (i.e., for a KxK matrix, there are K sets of K-1 coefs, so that’s K*(K-1) coefs).  Just to be clear, I didn’t simulate any data or run any regressions; rather, I simulated corr matrices and from there used the formula to get the regression coefs.

Anyway, for each combination of model (LKJ, WIshart, or inv-Wishart), phi parameter (1, 2, 5, or 10), and dimension (2, 5, 10, 30, or 100), I did 100 simulations and plotted the histogram of the resulting 100*K*(K-1) coefs.  That’s what’s in the attached document.

Interestingly, some of these distributions look Cauchy-like.

See you
A

P.S.  For LKJ, the “nu” parameter is the same as “phi.”  For Wishart or inv-Wishart, I set nu <- phi + K - 1 to respect the constraint that nu must be at least K-1.

[Andrew Gelman]

Hi Ben.  Oh, I didn’t realize you’d done all this.  I did some more simulations and it looks like the distribution of these coefs may approach some fixed dist as K approaches infinity.  See attached figure, where I did LKJ(1), Wishart(K+1) and inv-Wishart(K+1) (in the latter case normalizing so that they gave corr matrices rather than cov matrices) and then I looked at the histogram of all the regression coefs.  For the LKJ, at least, it looks like the dist of the betas might be approaching a fixed limit as K increases.

Those histograms are hard to read, so also here are some quantiles of the betas for LKJ(1) for each of several K’s:

    K = 5, 1/5/25/50/75/95/99 percentiles are -2.1, -1.1, -0.5,  0.0,  0.5,  1.2,  2.1

  K = 10, 1/5/25/50/75/95/99 percentiles are -2.3, -1.2, -0.4,  0.0,  0.4,  1.2,  2.2

  K = 20, 1/5/25/50/75/95/99 percentiles are -2.5, -1.3, -0.5,  0.0,  0.5,  1.3,  2.5

  K = 50, 1/5/25/50/75/95/99 percentiles are -2.5, -1.3, -0.4,  0.0,  0.4,  1.3,  2.5

K = 100, 1/5/25/50/75/95/99 percentiles are -2.6, -1.3, -0.4,  0.0,  0.4,  1.3,  2.6

Strikingly stable, I’d say.  And reasonably fit by a t distribution with 2.5 df and scale .5:

> pfround(quantile(.5*rt(100000, 2.5), c(.01, .05, .25, .5, .75, .95, .99)), 1)
  1%   5%  25%  50%  75%  95%  99% 
-2.6 -1.3 -0.4  0.0  0.4  1.3  2.6 

See you

A

[Ben Goodrich]

On Wednesday, December 16, 2015 at 6:43:09 PM UTC-5, Andrew Gelman wrote:

Strikingly stable, I’d say.  And reasonably fit by a t distribution with 2.5 df and scale .5:

> pfround(quantile(.5*rt(100000, 2.5), c(.01, .05, .25, .5, .75, .95, .99)), 1)
  1%   5%  25%  50%  75%  95%  99% 
-2.6 -1.3 -0.4  0.0  0.4  1.3  2.6

I agree that it seems to have a limit but

• I don't see why we need to match a t distribution to it since we already know its marginal density and could work on deriving its limit as K goes to infinity
  
• The marginal density of a standardized regression coefficient is a bit of a red herring anyway
  
• I don't think setting eta = 1 is an appropriate universal hyperparameter

The critical property of the LKJ-based distribution you are simulating is that

beta' * Cor(X) * beta = R^2 ~ Beta(K/2, eta)

If you just put a bunch of independent t priors on the elements of beta, then you still need a distribution for the correlation matrix of X and even then the quadratic form doesn't equal anything in particular. So, in my opinion, the correct way to interpret this prior is as an informative prior on the unknown R^2 but conditional on the R^2, there is virtually no additional information about the coefficients imparted by the prior.

Thus, I don't think eta = 1 is a great default in light of the Beta(K/2, eta) distribution for the R^2. Clearly, any eta < 1 is usually inappropriate because of the infinite density at 1. Setting eta = 1 puts finite density at 1 but the mode is still 1, which doesn't really make sense for an R^2. For regressions in genetics, I don't think they would presume that the R^2 converges to 1 as K diverges to infinity. Rather, they seem to think that there are a finite number of predictors that actually matter and an infinite number that do not, so the R^2 is bounded away from 1.

The stan_lm() function forces the user to say something about their prior beliefs over the location of the R^2 and derives the value of eta that corresponds to that, which is not universal but I think it is appropriate and reasonably easy for users to deal with.

Ben

[Andrew Gelman]

Hi, yes, I agree there are some missing steps here.  One idea is that this LKJ thing is that it provides rhetorical support for a certain weak but not completely empty prior.  And it does give us some baseline.  I could try some simulations with other values of eta?

A

[Ben Goodrich]

On Thursday, December 17, 2015 at 6:45:44 PM UTC-5, Andrew Gelman wrote:

Hi, yes, I agree there are some missing steps here.  One idea is that this LKJ thing is that it provides rhetorical support for a certain weak but not completely empty prior.  And it does give us some baseline.  I could try some simulations with other values of eta?

It is a bit misleading to say that the prior is "weak". A Beta(K/2, eta) prior on the R^2 is often quite strong in the sense that you can make it quite concentrated around 0.6 if you thought in advance that the R^2 was 0.6. But all combinations of the prior correlation matrix among the predictors and standardized coefficients that yield a given R^2 are pretty much equally likely under the prior. Using the stochastic representation of the elliptical distribution, we have

beta = sqrt(R^2) * L^{-1} * u

where L * L' is the correlation matrix among the X variables and u is uniformally distributed on the surface of a K-hypersphere. Treating L as an unknown hyperparameter makes this a mixture of elliptical distributions but if you ignore that, you are basically left with the spherical normal R * u and there is no prior information about u except for the restriction that u'u = 1. So, almost all the information in the joint prior is piped through a Beta(K/2, eta) distribution for the R^2.

Ben

[Andrew Gelman]

Hmmm, interesting, I’ll have to process that!