R2 prior in stan_lm

[Kevin Van Horn]

I just got done reading the vignette "Estimating Regularized Linear Models with rstanarm," and I'm having a hard time making sense of its description / derivation of the R2 prior. Is there a paper that goes into this prior in more detail?

[Ben Goodrich]

On Monday, February 15, 2016 at 8:31:37 PM UTC-5, Kevin Van Horn wrote:

I just got done reading the vignette "Estimating Regularized Linear Models with rstanarm," and I'm having a hard time making sense of its description / derivation of the R2 prior. Is there a paper that goes into this prior in more detail?

No, there isn't. I would say that the less you think about it, the easier it is to understand. What do you think the prior location (mode by default if K > 2) of the R^2 is? The answer to that question pins down the single hyperparameter to the joint prior on the coefficients. Marginally, each coefficient has a symmetric prior centered at zero with fairly heavy tails (although how heavy depends on the hyperparameter).

Ben

[Andrew Gelman]

We are planning to write a paper on this!

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

[Dustin Tran]

Is it possible to derive the equivalent prior on R2 given a prior on beta? (or conversely) The doc says a normal prior on beta corresponds to a non-uniform prior on R2, although it doesn’t say what it is more specifically.

Dustin

[Kevin Van Horn]

R^2 is defined for a particular design matrix X, coefficient vector beta, and response vector y... so it's unclear to me what is meant by a prior for R^2. Do you mean

1. The marginal prior for
    R^2(X, beta, \tilde{y}),
where
    \tilde{y} ~ Normal(X * beta, \sigma_{\eps})
and X and \sigma_{\eps} are taken as fixed?

2. The prior for
    ER2(beta) = E[R^2 | X, beta, \sigma_{\eps}],
where the expectation is over \tilde{y}, and X and \sigma_{\eps} are taken as fixed?

Or do you mean something else entirely?

[Ben Goodrich]

On Monday, February 15, 2016 at 10:22:12 PM UTC-5, Kevin Van Horn wrote:

R^2 is defined for a particular design matrix X, coefficient vector beta, and response vector y... so it's unclear to me what is meant by a prior for R^2.

That is the usual way of looking at the R^2 in a OLS context, yes. However, the R^2 is invariant to a rotation of the design matrix X, so we can utilize its QR decomposition to go from

y_hat = X * beta = (QR) * beta = Q * (R * beta) = Q * theta

where theta = R * beta. So, you can think of this as a regression of y on Q and then beta = R^{-1} * theta is a generated quantity.

Okay, so what is the prior distribution on theta, the vector of coefficients for Q? Since Q'Q = I, theta_k can be conceptualized as the correlation between y and the k-th column of Q multiplied by the standard deviation of y, divided by the standard deviation of the k-th column of Q which is the square root of N - 1.

That pushes back the question to what is your prior on the correlations between y and the columns of Q? Here we make the assumption that this vector of correlations has a spherical distribution, which can be represented as \sqrt{R^2} * u, where u is uniformally distributed on the surface of a hypersphere. Since u'u = 1, the sum-of-squared correlations is R^2, which is the proportion of variance in the outcome attributable to the predictors in a linear model.

And that brings us to the question, what is your prior on R^2? The LKJ paper implies that if you have a correlation matrix among both the predictors and the outcome, the distribution of the R^2 is beta with shape parameters K / 2 and eta. So, in principle, a user would need to choose eta to reflect their subjective beliefs about the R^2. However, the R2() function that is called by stan_lm() reverses this around and asks what do you think the prior mode / mean / median / log-expectation of the R^2 is? Given the answer to that question, it can solve for the implied value of eta.

All together, the primitive parameters are R^2, u, and sigma_y (the marginal standard deviation of Y). Since sigma_y is a scale parameter and not of much interest, it is given a Jeffreys' prior, which is uniform on the log scale. Given these three things (and the fixed R), we can calculate beta and sigma (the conditional standard deviation of Y) as generated or intermediate quantities, where sigma = sigma_y * sqrt(1 - R^2).

So, in summary, this is really just a way of eliciting an informative prior on the R^2 that happens to imply a pretty reasonable set of joint prior beliefs about beta (it's centered at the origin, has fat tails, is round at the mode, etc.). Consequently, relative to OLS the posterior medians will usually be shrunk toward zero a bit, which is usually a good thing for out-of-sample predictions.

Do you mean

1. The marginal prior for
    R^2(X, beta, \tilde{y}),
where
    \tilde{y} ~ Normal(X * beta, \sigma_{\eps})
and X and \sigma_{\eps} are taken as fixed?

You can choose the prior location of the R^2 given what you know about X, but it is not conditioning on beta or y or sigma. Rather, beta and sigma are functions of the R^2.
 

2. The prior for
    ER2(beta) = E[R^2 | X, beta, \sigma_{\eps}],
where the expectation is over \tilde{y}, and X and \sigma_{\eps} are taken as fixed?

If you specify what = "mean", then you are saying "Here is what I think the prior mean of R^2 is given X", but again it does not condition on beta or sigma. But the default is what = "mode", so it is asking what is the most likely value of R^2 given X.
 

Or do you mean something else entirely?

I think the key point is that beta and sigma are generated / intermediate quantities that are both functions of the primitive R^2. So, the prior for R^2 does not explicitly condition on beta or sigma in the mathematical sense. Nevertheless, saying "I think the most likely value for R^2 is 0.75" conveys some information about beta and sigma.

As with any stan_ function in rstanarm, you can get a sense for the prior distribution(s) by specifying prior_PD = TRUE, in which case it will run the model but not condition on the data so that you just get draws from the prior. In the case of stan_lm, the Jeffreys' prior on sigma_y is improper, so it just sets sigma_y = 1 when prior_PD = TRUE. So, the resulting distribution of beta is implicitly assuming a standardized outcome but it should be centered on zero with some extreme positive and negative values.

Ben

[Kevin Van Horn]

There is an inconsistency between your definition of theta and the vignette's definition of theta:

1. Your definition: theta = R * beta, or equivalently, theta = Q' * y_hat

2. Vignette's definition: theta = Q' * y

y and y_hat are not the same thing. Which definition is correct?

"theta_k can be conceptualized as the correlation between y and the k-th column of Q multiplied by..."

Here you seem to be using theta = Q' * y, instead of your earlier definition of theta = Q' * y_hat.

"You can choose the prior location of the R^2 given what you know about X, but it is not conditioning on beta or y or sigma. Rather, beta and sigma are functions of the R^2."

Borrowing some terminology from Interaction Design, you're talking about the implementation model, whereas I'm talking about the conceptual model. The conceptual definition of R^2 is in terms of X, beta, and y:

  R^2 = 1 - SS_r / SS_t

  SS_r = SUM(i:: (y_i - yhat_i)^2)

  SS_t = SUM(i:: (y_i - ybar)^2)

  ybar = empirical mean of y_i

or, an alternative and slightly different definition uses expected values according to the model:

  SS_r = E[SUM(i:: (y_i - yhat_i)^2)] = n * sigma_{eps}^2

  SS_t = E[SUM(i:: (y_i - ybar)^2)] = n * sigma_{eps}^2 + n * V_{i=1}^n yhat_i

  yhat_i = X_i * beta

I'm guessing that underlying the derivation of the R2 prior is one of these two conceptual definitions for R^2.

[Ben Goodrich]

On Tuesday, February 16, 2016 at 11:53:36 PM UTC-5, Kevin Van Horn wrote:

There is an inconsistency between your definition of theta and the vignette's definition of theta:

1. Your definition: theta = R * beta, or equivalently, theta = Q' * y_hat

2. Vignette's definition: theta = Q' * y

y and y_hat are not the same thing. Which definition is correct?

"theta_k can be conceptualized as the correlation between y and the k-th column of Q multiplied by..."

Here you seem to be using theta = Q' * y, instead of your earlier definition of theta = Q' * y_hat.

If you prefer, I could have said it with an error term rather than with the y_hat

y = X * beta + epsilon = (QR) * beta + epsilon = Q * (R * beta) + epsilon = Q * theta + epsilon

where theta = R * beta. But this is all prior to observing the y vector. We can ask ourselves "What are our prior beliefs about theta = (Q'Q)^{-1} Q'y = Q'y? For fixed R, whatever those prior beliefs about theta are they imply some prior beliefs about beta = R^{-1} theta. Specifically, if theta has a spherical prior distribution, then beta has an elliptical prior distribution. And the stochastic representation of a spherical distribution implies

theta = sigma_y * sqrt(N - 1) * sqrt(R^2) * u

where u is uniform on the surface of a hypersphere. So, you need some prior distribution for sigma_y and some prior distribution for R^2, which can be interpreted as the proportion of variance in y attributable to Q (or X) in a linear model.
 

"You can choose the prior location of the R^2 given what you know about X, but it is not conditioning on beta or y or sigma. Rather, beta and sigma are functions of the R^2."

Borrowing some terminology from Interaction Design, you're talking about the implementation model, whereas I'm talking about the conceptual model. The conceptual definition of R^2 is in terms of X, beta, and y:

  R^2 = 1 - SS_r / SS_t

  SS_r = SUM(i:: (y_i - yhat_i)^2)

  SS_t = SUM(i:: (y_i - ybar)^2)

  ybar = empirical mean of y_i

or, an alternative and slightly different definition uses expected values according to the model:

  SS_r = E[SUM(i:: (y_i - yhat_i)^2)] = n * sigma_{eps}^2

  SS_t = E[SUM(i:: (y_i - ybar)^2)] = n * sigma_{eps}^2 + n * V_{i=1}^n yhat_i

  yhat_i = X_i * beta

I'm guessing that underlying the derivation of the R2 prior is one of these two conceptual definitions for R^2.

Yes, in any textbook on regression from a frequentist perspective, you see those definitions for R^2. But in that setup, the R^2 is not a parameter in the model, it is a statistic that is calculated afterward as a function of the data and the estimates of beta. What I am saying is that in lm.stan, the R^2 is a primitive parameter and the posterior distribution of beta = R^{-1} theta is a function of the R^2 (via theta) and the data (which define the R matrix in QR = X). So, in my view, it does not make sense to conceptually define the R^2 in terms of an estimate of beta because beta is defined in the lm.stan code in terms of the primitive parameter R^2. But conceptually, the R^2 can still be interpreted as the proportion of variance in the outcome attributable to the predictors for the purpose of expressing your prior beliefs about its location. In other words, since sigma is an intermediate quantity that is defined as sigma_y * sqrt(1 - R^2), the model implies that the variance of the outcome is

(sigma_y)^2 * (R^2 + (1 - R^2)) = (sigma_y)^2

where sigma_y is also a primitive unknown that is estimated but is not especially interesting. So, you can ascribe some proportion of that to the predictors and the complementary proportion to the errors.

To conceptualize the R^2 as a primitive parameter, it helps to envision a random correlation matrix among both the predictors and the outcome. Let L be the Cholesky factor of that correlation matrix. The so-called onion method for constructing L implies that the bottom-right element of L is sqrt(1 - R^2) and the elements in the last row of L to the left of the diagonal can be written as sqrt(R^2) * u, so that the sum-of-squares along the last row of L is 1. If your beliefs about that correlation matrix are described by the LKJ distribution with shape parameter eta, then your beliefs about R^2 are a beta distribution with first shape parameter K / 2 and second shape parameter eta. That is why some prior information about R^2 is solicited in order to pin down the second shape parameter, eta.

Ben