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

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

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

• My understanding is that ridge regression (and its cousins such as the Lasso and elastic net) shrinks coefficients for correlated variables shared by all the observations in the regression (eg, student's socioeconomic status and GPA) while a random effects model performs shrinkage on coefficients for mutually exclusive levels or groups of correlated observations (such as student's student's socioeconomic status grouped by school id). – RobertF Oct 30 '14 at 18:59
• I think the best place for getting a unifying answer is to look at the keyword BLUP (for Best Linear Unbiased Predictor) esp. in the animal breeding literature. See for instance Robinson's survey in Statistical Science. Or Marvin Gruber's book – Xi'an Nov 6 '14 at 14:58
• @Xi'an: Thanks a lot, I have already found Gruber's book myself, and even though he certainly discusses both James-Stein and ridge regression a lot, I did not immediately find any direct comparison of the two (reading the whole book is not an option for me right now...). Thanks for the link to Robinson's survey, I will take a look; animal breeding! who would have thought. By the way, I have seen your comments on related threads, and guess that you might be one of the people who could actually provide a satisfactory answer here! This would be great; so far no answer makes me content. – amoeba Nov 6 '14 at 15:30
• @Xi'an: Well, your helpful comments below do make me miss an answer from you here. Anyway, I started reading Robinson and realized that "Best Linear Unbiased Predictor" is a biased estimator (obviously, as it implements shrinkage)! What a nice terminology. – amoeba Nov 7 '14 at 21:00
• They are good at names in animal breeding: after Casella & George 1992 "Gibbs for kids" had to change its title to get published, Wang & Gianola wrote a "Gibbs for pigs" introduction in 1993 at an European Association for Animal Production meeting! – Xi'an Nov 7 '14 at 21:04

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

• (+1) Thanks a lot! This is very helpful, and I will definitely look in Bishop's textbook which I know well and often consult. I did not expect to find anything on mixed models there, but it looks like Section 3.3 "Bayesian Linear Regression" is actually about that, just uses different terminology. Very good to know! But what is your take on my bullet questions? – amoeba Nov 2 '14 at 15:28
• You have too many questions in a post. :) 1) As I answered above, James-Stein estimator and ridge regression are equivalent when there is no covariates $\mathbf X$, or just an identity matrix. 2,3,4) as @James mentioned, the number of predictors ($p$ above) are not necessarily equal to the response dimension $m$. – Randel Nov 2 '14 at 15:54
• BTW, I cannot see sample average/mean is used in James-Stein estimator, it actually takes the estimator $\mathbf y$ and then shrinks it to $\mathbf 0$. – Randel Nov 2 '14 at 16:05
• The JS estimator and ridge regression are distinct. A ridge regression estimate of the $p$-dimensional location of a p-vector corresponds to the design matrix $\mathbf{I}_p$, which would lead to the estimate $(1+\lambda)^{-1}\mathbf{I}_p \mathbf y$, which is missing the (non-linear!) $\Vert \mathbf y \Vert^2$ term in the denominator of the JS-estimator – Andrew M Nov 4 '14 at 19:00
• I think it all depends what you call ridge estimator. In the early Hoerl and Kennard (1970) sense, there is indeed no dependence of $\lambda$ on the data. In the later sense of Casella's PhD thesis (1978), the manual determination of $\lambda$ is replaced by a function of the residual sum of squares. – Xi'an Nov 7 '14 at 19:58

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

• Thanks a lot, very interesting (+1). I can only wish your answer were more detailed... Re your footnote (3): are you saying that all Bayes estimators are admissible/indomitable (I like the word), independent of the prior? But James-Stein estimator can be derived from empirical Bayes; why is it inadmissible then? Also, it would mean that e.g. in ridge regression I can take a prior concentrated not around zero, but around some other value: $\beta \sim \mathcal N(\beta_0, 1/\lambda^2)$, and it will still be a reasonable regularization strategy? – amoeba Nov 4 '14 at 12:19
• For the reason why James-Stein estimator is inadmissible, you can find the answer here. There is also a detailed and interesting discussion in Lehmann & Casella (1998), Theory of Point Estimation. – Randel Nov 4 '14 at 16:31
• @Randel: yes, I know that it is inadmissible and have seen that reasoning, I am just wondering how it fits to Andrew's statement (given I understood it correctly) that all Bayes estimators are admissible, since James-Stein can be understood via Empirical Bayes... – amoeba Nov 4 '14 at 16:53
• @Amoeba: yes, any Bayes estimator that is the posterior under any proper prior leads to an admissible estimator. As far as empirical Bayes goes, such procedures are not in fact bonafide Bayes, because having the prior depend on the data can lead to pathologies. Sometimes they can be shown to be admissible, sometimes they are not--typically you have to work case-by-case. I have edited my answer to be a little bit more cagey on this point, because in fact I don't know if classical linear mixed models are admissible! – Andrew M Nov 4 '14 at 18:49
• Just need to point out that genuine proper Bayes estimators rarely work as James-Stein estimators because they are not minimax. Bill Strawderman showed for instance (in 1975) that there exists no minimax proper Bayes estimator in dimensions less than 5 for the usual normal mean problem that set it all. – Xi'an Nov 7 '14 at 19:50
• 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.

• Thank you for the answer. Regarding your first bullet: but what is being shrunken in ridge regression is $\beta$, which has as many dimensions as predictors, doesn't it? – amoeba Oct 30 '14 at 15:44
• Ok, then in theory JS should work better, assuming that it has been extended to the case when MSE is estimated and the variance-covariance matrix of beta is arbitrary. In that case, JS will not just take the point estimate of beta and multiply it by a scaling factor. Similary to Ridge Regression, different components of beta will be shrunk differently. – James Oct 30 '14 at 15:53
• Very good point about covariance matrix of $\beta$! I guess this answers (at least intuitively) my first bullet. – amoeba Oct 30 '14 at 16:07
• @James: linear models can be thought of projecting the sample (which lives in $R^n$) onto a $p$-dimensional subspace (the columns spanned by the design matrix.) In particular, we could always trivially project it onto the identity, which is the same as using sample mean of an $n$-vector when you only have a single observation. – Andrew M Nov 4 '14 at 9:20

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.

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