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?


2 Answers 2


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.

  • 1
    $\begingroup$ I like this answer. I also liked the question. Personally, I am still struggling on the mechanism (I actually never cared about it to study the algorithms that solve LMEM's). So I guess that the difference of the random effects is being made by changing from $$\mathbf{y} \sim \mathcal{N}(a + b\mathbf{x}, \sigma^2 I)$$ to $$\mathbf{y} \sim \mathcal{N}(a + b\mathbf{x}, \Sigma) $$ I imagine that a tiny example that works this out might be nice. I am considering to make this myself, but maybe there are resources that already show such examples (anyone?). $\endgroup$ Commented Feb 23, 2018 at 16:27
  • $\begingroup$ @statmerkur Tau is a parameter; the last formula in my answer still includes tau. The crucial point is that the last formula does NOT include $c$. We simply combine two equations together such that $c$ falls out of there (technically, we integrate over $c$). Then we fit the model, which means we fit tau and other parameters. $\endgroup$
    – amoeba
    Commented Feb 24, 2018 at 22:41
  • $\begingroup$ I think I just don't get the integration step. As @Martijn Weterings pointed out a little (R code) example or reference were one can find this would be great! $\endgroup$
    – statmerkur
    Commented Feb 25, 2018 at 2:41
  • $\begingroup$ Thanks for accepting my answer and awarding me the bounty @statmerkur, but it's too bad that it remains unclear. I will try to think of an example. I will ping you when I update the answer. $\endgroup$
    – amoeba
    Commented Feb 25, 2018 at 20:52
  • $\begingroup$ @statmerkur In an answer to this question I demonstrate the manual calculation of a mixed effects model (manual in the sense of writing the likelihood function, the optimization is still done by a standard optimization function in R) stats.stackexchange.com/a/337348/164061 $\endgroup$ Commented Jun 25, 2018 at 13:54

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

  • 1
    $\begingroup$ I am referring to the (co)variance parameters of the random effects (see my edit). $\endgroup$
    – statmerkur
    Commented Jan 28, 2018 at 22:54
  • 2
    $\begingroup$ I don't think this answers the question. $\endgroup$
    – amoeba
    Commented Feb 19, 2018 at 11:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.