Intuition about parameter estimation in mixed models (variance parameters vs. conditional modes)

Question

I have read many times that random effects (BLUPs/conditional modes for, say, subjects) are not parameters of a linear mixed effects model but instead can be derived from the estimated variance/covariance parameters. E.g. Reinhold Kliegl et al. (2011) state:

Random effects are subjects’ deviations from the grand mean RT and subjects’ deviations from the fixed-effect parameters. They are assumed to be independently and normally distributed with a mean of 0. It is important to recognize that these random effects are not parameters of the LMM – only their variances and covariances are. [...] LMM parameters in combination with subjects’ data can be used to generate “predictions” (conditional modes) of random effects for each subject.

Can someone give an intuitive explanation how the (co)variance parameters of the random effects can be estimated without actually using/estimating the random effects?

Answer 1

Consider a simple linear mixed model, e.g. a random intercept model where we estimate the dependency of $y$ on $x$ in different subjects, and assume that each subject has their own random intercept:$$y = a + bx + c_i + \epsilon.$$ Here intercepts $c_i$ are modeled as coming from a Gaussian distribution $$c_i\sim \mathcal N(0, \tau^2)$$ and random noise is also Gaussian $$\epsilon \sim \mathcal N(0, \sigma^2).$$ In the lme4 syntax this model would be written as y ~ x + (1|subject).

It is instructive to rewrite the above as follows:

\begin{gather} y \mid c \sim \mathcal N(a + bx + c, \sigma^2) \\ c \sim \mathcal N(0, \tau^2) \end{gather}

This is a more formal way to specify the same probabilistic model. From this formulation we can directly see that the random effects $c_i$ are not "parameters": they are unobserved random variables. So how can we estimate the variance parameters without knowing the values of $c$?

Note that the first equation above describes the conditional distribution of $y$ given $c$. If we know the distribution of $c$ and of $y\mid c$, then we can work out the unconditional distribution of $y$ by integrating over $c$. You might know it as the Law of total probability. If both distributions are Gaussian, then the resulting unconditional distribution is also Gaussian.

In this case the unconditional distribution is simply $\mathcal N(a + bx, \sigma^2+\tau^2)$, but our observations are not i.i.d. samples from it because there are multiple measurements per subject. In order to proceed, we need to consider the distribution of the whole $n$-dimensional vector $\mathbf y$ of all observations: $$\mathbf y \sim \mathcal N(a+b\mathbf x, \boldsymbol\Sigma)$$ where $\boldsymbol\Sigma=\sigma^2 \mathbf I_n + \tau^2 \mathbf I_N \otimes \mathbf 1_M$ is a block-diagonal matrix composed of $\sigma^2$ and $\tau^2$. You asked for intuition so I want to avoid the math. The important point is that this equation does not have $c$ anymore! This is what one actually fits to the observed data, and that is why one says that $c_i$ are not the parameters of the model.

When the parameters $a$, $b$, $\tau^2$, and $\sigma^2$ are fit, one can work out the conditional distribution of $c_i$ for each $i$. What you see in the mixed model output are the modes of these distributions, aka the conditional modes.

Answer 2

You can easily estimate variance and covariance parameters without relying on random-effects by using fixed-effects (see here for a discussion fixed-effects vs. random-effects; be aware of the fact that there are different definitions of these terms).

Fixed-effects can be easily derived by adding a (binary) indicator variable for each group (or each time period or whatever you are thinking to use as random-effects; this is equivalent to the within transformation). This allows you easily to estimate the fixed-effects (which can be viewed as a parameter).

Fixed-effects assumption does not require you to make an assumption of the distribution of the fixed-effects, you can easily estimate the variance of the fixed-effects (although this extremely noise if the number of observation within each group is small; they minimize the bias for the expense of much larger variance compared to the random-effects because you lose one degree of freedom for each group through adding these indicator variables). You can also estimate covariances between different sets of fixed-effects, or between fixed-effects and other covariates. We have done that for example in a paper called Competitive Balance and Assortative Matching in the German Bundesliga to estimate whether better football players increasingly play for better teams.

Random-effects need a prior assumption about the covariance. In classical random-effects models, you assume that the random-effects are like an error and they are independent of the other covariates (so that you can ignore them and use OLS and get still consistent albeit inefficient estimates for the other parameter if the assumptions of the random-effects model hold true).

Further more technical information is available here. Andrew Gelman has also a lot of more intuitive work about this in his nice book Data analysis using regression and multilevel/hierarchical models