Unified view on shrinkage: what is the relation (if any) between Stein's paradox, ridge regression, and random effects in mixed models?

Question

Consider the following three phenomena.

1. Stein's paradox: given some data from multivariate normal distribution in $\mathbb R^n, \: n\ge 3$, sample mean is not a very good estimator of the true mean. One can obtain an estimation with lower mean squared error if one shrinks all the coordinates of the sample mean towards zero [or towards their mean, or actually towards any value, if I understand correctly].

   NB: usually Stein's paradox is formulated via considering only one single data point from $\mathbb R^n$; please correct me if this is crucial and my formulation above is not correct.

2. Ridge regression: given some dependent variable $\mathbf y$ and some independent variables $\mathbf X$, the standard regression $\beta = (\mathbf X^\top \mathbf X)^{-1} \mathbf X^\top \mathbf y$ tends to overfit the data and lead to poor out-of-sample performance. One can often reduce overfitting by shrinking $\beta$ towards zero: $\beta = (\mathbf X^\top \mathbf X + \lambda \mathbf I)^{-1} \mathbf X^\top \mathbf y$.

3. Random effects in multilevel/mixed models: given some dependent variable $y$ (e.g. student's height) that depends on some categorical predictors (e.g. school id and student's gender), one is often advised to treat some predictors as 'random', i.e. to suppose that the mean student's height in each school comes from some underlying normal distribution. This results in shrinking the estimations of mean height per school towards the global mean.

I have a feeling that all of this are various aspects of the same "shrinkage" phenomenon, but I am not sure and certainly lack a good intuition about it. So my main question is: is there indeed a deep similarity between these three things, or is it only a superficial semblance? What is the common theme here? What is the correct intuition about it?

In addition, here are some pieces of this puzzle that don't really fit together for me:

• In ridge regression, $\beta$ is not shrunk uniformly; ridge shrinkage is actually related to singular value decomposition of $\mathbf X$, with low-variance directions being shrunk more (see e.g. The Elements of Statistical Learning 3.4.1). But James-Stein estimator simply takes the sample mean and multiplies it by one scaling factor. How does that fit together?

  Update: see James-Stein Estimator with unequal variances and e.g. here regarding variances of $\beta$ coefficients.

• Sample mean is optimal in dimensions below 3. Does it mean that when there are only one or two predictors in the regression model, ridge regression will always be worse than ordinary least squares? Actually, come to think of it, I cannot imagine a situation in 1D (i.e. simple, non-multiple regression) where ridge shrinkage would be beneficial...

  Update: No. See Under exactly what conditions is ridge regression able to provide an improvement over ordinary least squares regression?

• On the other hand, sample mean is always suboptimal in dimensions above 3. Does it mean that with more than 3 predictors ridge regression is always better than OLS, even if all the predictors are uncorrelated (orthogonal)? Usually ridge regression is motivated by multicollinearity and the need to "stabilize" the $(\mathbf X^\top \mathbf X)^{-1}$ term.

  Update: Yes! See the same thread as above.

• There are often some heated discussion about whether various factors in ANOVA should be included as fixed or random effects. Shouldn't we, by the same logic, always treat a factor as random if it has more than two levels (or if there are more than two factors? now I am confused)?

  Update: ?

---

Update: I got some excellent answers, but none provides enough of a big picture, so I will let the question "open". I can promise to award a bounty of at least 100 points to a new answer that will surpass the existing ones. I am mostly looking for a unifying view that could explain how the general phenomenon of shrinkage manifests itself in these various contexts and point out the principal differences between them.

Answer 1

Connection between James–Stein estimator and ridge regression

Let $\mathbf y$ be a vector of observation of $\boldsymbol \theta$ of length $m$, ${\mathbf y} \sim N({\boldsymbol \theta}, \sigma^2 I)$, the James-Stein estimator is, $$\widehat{\boldsymbol \theta}_{JS} = \left( 1 - \frac{(m-2) \sigma^2}{\|{\mathbf y}\|^2} \right) {\mathbf y}.$$ In terms of ridge regression, we can estimate $\boldsymbol \theta$ via $\min_{\boldsymbol{\theta}} \|\mathbf{y}-\boldsymbol{\theta}\|^2 + \lambda\|\boldsymbol{\theta}\|^2 ,$ where the solution is $$\widehat{\boldsymbol \theta}_{\mathrm{ridge}} = \frac{1}{1+\lambda}\mathbf y.$$ It is easy to see that the two estimators are in the same form, but we need to estimate $\sigma^2$ in James-Stein estimator, and determine $\lambda$ in ridge regression via cross-validation.

Connection between James–Stein estimator and random effects models

Let us discuss the mixed/random effects models in genetics first. The model is $$\mathbf {y}=\mathbf {X}\boldsymbol{\beta} + \boldsymbol{Z\theta}+\mathbf {e}, \boldsymbol{\theta}\sim N(\mathbf{0},\sigma^2_{\theta} I), \textbf{e}\sim N(\mathbf{0},\sigma^2 I).$$ If there is no fixed effects and $\mathbf {Z}=I$, the model becomes $$\mathbf {y}=\boldsymbol{\theta}+\mathbf {e}, \boldsymbol{\theta}\sim N(\mathbf{0},\sigma^2_{\theta} I), \textbf{e}\sim N(\mathbf{0},\sigma^2 I),$$ which is equivalent to the setting of James-Stein estimator, with some Bayesian idea.

Connection between random effects models and ridge regression

If we focus on the random effects models above, $$\mathbf {y}=\mathbf {Z\theta}+\mathbf {e}, \boldsymbol{\theta}\sim N(\mathbf{0},\sigma^2_{\theta} I), \textbf{e}\sim N(\mathbf{0},\sigma^2 I).$$ The estimation is equivalent to solve the problem $$\min_{\boldsymbol{\theta}} \|\mathbf{y}-\mathbf {Z\theta}\|^2 + \lambda\|\boldsymbol{\theta}\|^2$$ when $\lambda=\sigma^2/\sigma_{\theta}^2$. The proof can be found in Chapter 3 of Pattern recognition and machine learning.

Connection between (multilevel) random effects models and that in genetics

In the random effects model above, the dimension of $\mathbf y$ is $m\times 1,$ and that of $\mathbf Z$ is $m \times p$. If we vectorize $\mathbf Z$ as $(mp)\times 1,$ and repeat $\mathbf y$ correspondingly, then we have the hierarchical/clustered structure, $p$ clusters and each with $m$ units. If we regress $\mathrm{vec}(\mathbf Z)$ on repeated $\mathbf y$, then we can obtain the random effect of $Z$ on $y$ for each cluster, though it is kind of like reverse regression.

Acknowledgement: the first three points are largely learned from these two Chinese articles, 1, 2.

Answer 2

I'm going to leave it as an exercise for the community to flesh this answer out, but in general the reason why shrinkage estimators will *dominate*$^1$ unbiased estimators in finite samples is because Bayes$^2$ estimators cannot be dominated$^3$, and many shrinkage estimators can be derived as being Bayes.$^4$

All of this falls under the aegis of Decision Theory. A exhaustive, but rather unfriendly reference is "Theory of point estimation" by Lehmann and Casella. Maybe others can chime in with friendlier references?

---

$^1$ An estimator $\delta_1(X)$ of parameter $\theta \in \Omega$ on data $X$ is dominated by another estimator $\delta_2(X)$ if for every $\theta \in \Omega$ the Risk (eg, Mean Square Error) of $\delta_1$ is equal or larger than $\delta_2$, and $\delta_2$ beats $\delta_1$ for at least one $\theta$. In other words, you get equal or better performance for $\delta_2$ everywhere in the parameter space.

$^2$ An estimator is Bayes (under squared-error loss anyways) if it is the the posterior expectation of $\theta$, given the data, under some prior $\pi$, eg, $\delta(X) = E(\theta | X)$, where the expectation is taken with the posterior. Naturally, different priors lead to different risks for different subsets of $\Omega$. An important toy example is the prior $$\pi_{\theta_0} = \begin{cases} 1 & \mbox{if } \theta = \theta_0 \\ 0 & \theta \neq \theta_0 \end{cases} $$ that puts all prior mass about the point $\theta_0$. Then you can show that the Bayes estimator is the constant function $\delta(X) = \theta_0$, which of course has extremely good performance at and near $\theta_0$, and very bad performance elsewhere. But nonetheless, it cannot be dominated, because only that estimator leads to zero risk at $\theta_0$.

$^3$ A natural question is if any estimator that cannot be dominated (called admissible, though wouldn't indomitable be snazzier?) need be Bayes? The answer is almost. See "complete class theorems."

$^4$ For example, ridge regression arises as a Bayesian procedure when you place a Normal(0, $1/\lambda^2$) prior on $\beta$, and random effect models arise as an empirical Bayesian procedure in a similar framework. These arguments are complicated by the fact that the vanilla version of the Bayesian admissibility theorems assume that every parameter has a proper prior placed on it. Even in ridge regression, that is not true, because the "prior" being placed on variance $\sigma^2$ of error term is the constant function (Lebesgue measure), which is not a proper (integrable) probability distribution. But nonetheless, many such "partially" Bayes estimators can be shown to be admissible by demonstrating that they are the "limit" of a sequence of estimators that are proper Bayes. But proofs here get rather convoluted and delicate. See "generalized bayes estimators".

Answer 3

• James-Stein assumes that the dimension of response is at least 3. In the standard ridge regression the response is one-dimensional. You are confusing the number of predictors with the response dimension.

• That being said, I see the similarity among those situations, but what exactly to do, e.g. whether a factor should be fixed or random, how much shrinkage to apply, if at all, is dependent on the particular dataset. E.g., the more orthogonal the predictors are, the less it makes sense to pick Ridge regression over standard regression. The larger the number of parameters, the more it makes sense to extract the prior from the dataset itself via Empirical Bayes and then use it for shrinking the parameter estimates. The higher the signal-to-noise ratio, the smaller the benefits of shrinkage, etc.

Answer 4

As others have said, the connection between the three is how you incorporate the prior information into the measurement.

1. In case of the Stein paradox, you know that the true correlation between the input variables should be zero (and all the possible correlation measures, since you want to imply independence, not just uncorrelatedness), hence you can construct a variable better than the simple sample mean and suppress the various correlation measures. In the Bayesian framework, you can construct a prior that literally down weighs the events that lead to correlation between the sample means and up weighs the others.
2. In case of ridge regression you want to find a good estimate for the conditional expectation value E(y|x). In principle this is a infinite-dimensional problem and ill-defined since we have only finite number of measurements. However, the prior knowledge is that we are looking for a continuos function that models the data. This is still ill-defined, since there are still infinitely many ways to model continuos functions, but the set is somewhat smaller. Ridge regression is just one simple way to sort the possible continuos functions, test them and stop at a final degree of freedom. An interpretation is the VC-dimension picture: during the ridge regression, you check that how well a f(x, p1, p2... ) model with a given degree of freedom describes the uncertainty inherent in the data. Practically, it measures how well can the f(x, p1, p2 ... ) and the empirical P(p1,p2...) can reconstruct the full P(y|x) distribution and not just E(y|x). This way the models with too many degree of freedom (which usually overfit) are weighed down, since more parameter mean after a certain degree of freedom will give larger correlations between the parameters and consequently much wider P(f(x, p1, p2... ) ) distributions. An other interpretation is that the original loss function is a measure value as well, and it the evaluation on a given sample comes with an uncertainty, so the real task is not minimizing the loss function but to find a minimum that is significantly lower than the others (practically changing from one degree of freedom to an other is a Bayesian decision, so one changes the number of parameters only if they give a significant decrease in the loss function). The ridge regression can be interpreted as an approximation to these two pictures (CV-dimension, expected loss). In some cases you want to prefer higher degrees of freedoms, for example in particle physics you study particle collision where you expect the produced number of particles to be a Poisson distribution, so you reconstruct the particle track from on an image (a photo for example) in a way that prefers a given number of tracks and suppresses models which has smaller or higher track-number-interpretation of the image.
3. The third case also tries to implement a prior information into the measurement, namely that it is known from previous measurements that the students' height can be modeled very well by Gaussian distributions and not by a Cauchy, for example.

So in short, the answer is that you can shrink the uncertainty of a measurement if you know what to expect and categorize the data with some previous data (the prior information). This previous data is what constrains your modeling function that you use to fit the measurements. In simple cases you can write down your model in the Bayesian framework, but sometimes it is impractical, like in integrating over the all the possible continuos functions to find the one that has the Bayesian Maximal A Posterior value.

Answer 5

James Stein estimator and Ridge regression

Consider

$\mathbf y=\mathbf{X}\beta+\mathbf{\epsilon}$

With $\mathbf{\epsilon}\sim N(0,\sigma^2I)$

Least square solution is of the form

$\hat \beta= \mathbf S^{-1}\mathbf{X}'\mathbf{y}$ where $\mathbf S= \mathbf X'\mathbf X$.

$\hat \beta $ is unbiased for $\beta$ and has covriance matrix $\sigma^2 \mathbf S^{-1}$. Therefore we can write

$\hat \beta \sim N(\beta, \sigma^2\mathbf S^{-1})$ Note that $\hat \beta $ are the the Maximum likelihood estimates, MLE.

James Stein

For simplicity for the Jame Stein we will assume $\mathbf S=\mathbf I$. James and Stein will then add a prior on the $\beta$, of the form

$\beta \sim N(0,a\mathbf I)$

And will get a posterior of the form $\frac{a}{a+\sigma^2}\hat \beta=(1-\frac{\sigma^2}{a+\sigma^2})\hat \beta$, they will then estimate $\frac{1}{a+\sigma^2}$ with $\frac{p-2}{\|\hat \beta\|^2}$ and get a James Stein estimator of the form

$\hat \beta=(1-\frac{p-2}{\|\hat \beta\|^2})\hat \beta$.

Ridge Regression

In ridge regression $\mathbf X$ is usually standadised (mean 0, vairance 1 for each column of $\mathbf X$ ) so that the regression parameters $\beta=(\beta_1,\beta_2,\ldots, \beta_p)$ are comparable. When this is $S_{ii}=1$ for $i=1,2,\ldots,p$.

A ridge regression estimate of $\beta$ is defined as, $\lambda\geq0$, to be

$\hat \beta (\lambda) =(\mathbf S+\lambda I)^{-1}\mathbf X'\mathbf y=(\mathbf S +\lambda\mathbf I)^{-1}\mathbf S \hat \beta$ note that $\hat \beta$ is the MLE.

How was $\hat \beta (\lambda)$ derived ?? Recall

$\hat \beta \sim N(\hat \beta, \sigma^2\mathbf S^{-1})$ and if we add a Bayesian prior

$\beta\sim N(0,\frac{\sigma^2}{\lambda}\mathbf I)$

Then we get

$\text{E}\left(\beta|\hat \beta\right)=(\mathbf S +\lambda\mathbf I)^{-1}\mathbf S \hat \beta$

Same as the ridge regression estimate $\hat \beta (\lambda)$. So the original form of the James Stein given here takes $\mathbf S=\mathbf I$ and $a=\frac{\sigma^2}{\lambda}$.