Most mixed effects models assume that the random effects $\gamma$ follow a $MVN(0, \Sigma_{\gamma})$ distribution. In some cases, specific structure is put on $\Sigma_{\gamma}$. For now, let's just assume it's an unstructured covariance matrix.

Now, if we don't apply any priors or penalties on $\Sigma_{\gamma}$, it seems to me that the likelihood is unbounded. In particular, if we set $\gamma = 0$, then as $\det(\Sigma_{\gamma}) \rightarrow 0$, the contribution of the log-density of $\gamma$ approaches infinity. As long as the other contributions to the likelihood are finite (which is typically the case), this implies the log-likelihood would be unbounded.

Clearly, for the MLE this creates an issue (although I know REML is a more popular alternative). For most MCMC algorithms, having unbounded log-densities can be problematic as well, even if the posterior is still proper.

How is this issue typically handled? Is there a canonical penalty/prior on $\Sigma_{\gamma}$?

  • 1
    $\begingroup$ I've wondered this too and would love to see an authoritative answer $\endgroup$
    – jld
    Aug 24, 2018 at 15:21
  • 1
    $\begingroup$ @jld: Much to my surprise (since mixed effects models are so common), I'm finding there doesn't seem to be an authoritative answer! I've been looking reading mostly in the Bayesian framework, and it seems for a long time people were using the conjugate Wishart prior. However, in recent years this has come in to question (hard to chose prior in informative way) and it seems now everyone has their own idea. I was surprised to see it seems the dust has still not settled! $\endgroup$
    – Cliff AB
    Aug 24, 2018 at 15:36
  • 2
    $\begingroup$ I am not sure I understand the reasoning here. What do you mean by "If we set gamma=0"? Consider the simplest mixed model y ~ (1 | id). It fits two parameters, the error variance and the intercept variance. Are you saying that if the intercept variance goes to 0, then the likelihood goes to infinity? $\endgroup$
    – amoeba
    Aug 24, 2018 at 16:43
  • 1
    $\begingroup$ I still don't understand why. In we marginalize out the random effects, the $y$ vector comes out from the Gaussian distribution with block-diagonal covariance matrix. It has $\sigma_\text{noise}^2+\sigma_\text{intercept}^2$ on the diagonal and $\sigma_\text{intercept}^2$ off-diagonal inside each block. If $\sigma_\text{intercept}=0$, then this reduces to the likelihood of the fixed-effects regression which is non-zero. I don't see what should go to infinity here. CC @jld. $\endgroup$
    – amoeba
    Aug 24, 2018 at 18:57
  • 1
    $\begingroup$ This problem occurs when one wants to cluster using multivariate normal mixtures. Here the likelihood is infinite when sampling units are clustered in their own singleton clusters - multivariate normals with mean equal to the observation and covariance matrix that has determinant zero. This paper "A smooth nonparametric estimate of a mixing distribution using mixtures of Gaussians" treats the problem by putting a lower limit to the determinants. In Bayesian stats, the prior doesn't put mass in matrices with zero determinants, and hence the posterior doesn't either. $\endgroup$
    – papgeo
    Aug 24, 2018 at 22:32

1 Answer 1


Taking the simple model from amoeba

$$y_{ij} \sim N(\mu_i,\sigma_f^2) \qquad \text{with} \qquad \mu_i \sim N(0,\sigma_r^2)$$

The probability density to observe a sample of $\mathbf{ y_{ij} }$ is:

$$f_{\mathbf{ Y_{ij} }}(\mathbf{ y_{ij} }) =det((2\pi)^k\Sigma)^{-\frac{1}{2}} e^{\mathbf{ y_{ij}^T\Sigma y_{ij}}}$$

With $\Sigma$ having a block structure like

$$\Sigma = \begin{bmatrix} J_1 & 0 & \dots &0 \\ 0 & J_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & J_n \\ \end{bmatrix}$$ and the blocks are like

$$J_i = \begin{bmatrix} \sigma_f^2+\sigma_r^2 & \sigma_r^2 & \dots & \sigma_r^2 & \sigma_r^2 \\ \sigma_r^2 & \sigma_f^2+\sigma_r^2 & \dots & \sigma_r^2 & \sigma_r^2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \sigma_r^2 & \sigma_r^2 & \dots & \sigma_f^2+\sigma_r^2 & \sigma_r^2 \\ \sigma_r^2 & \sigma_r^2 & \dots & \sigma_r^2 & \sigma_f^2+\sigma_r^2 \end{bmatrix} $$

Nothing goes wrong when $\sigma_r \to 0$

...except $f_{\mu_i}(0) \to \infty$ becomes a degenerate distribution, which is however not relevant for the calculation/expression of the distribution $f_{\mathbf{Y_{ij} }}$.

If you consider the space of points $Y_{ij},\mu_i$ then you can see all the probability concentrating on a hyper-surface with $\mu_i=0$ and the density $f_{\mu_i}$ (along with $f_{Y_{ij},\mu_i}$) goes to infinity on this surface. But instead of $f_{Y_{ij},\mu_i}$ you wish to calculate $f_{Y_{ij}}$ $$f_{Y_{ij}}(y_{ij}) = \int f_{Y_{ij},\mu_i}(y_{ij},\mu_i) d\mu_i = \int f_{Y_{ij}|\mu_i}(y_{ij},\mu_i) f_{\mu_i}(\mu_i) d\mu_i$$ This density distribution $f_{Y_{ij}|\mu_i}$ for the distribution of $Y_{ij}$ on the hyper-surfaces with coordinates $\mu_i$ does not go to infinity. Or from another viewpoint, you integrate $f_{\mu_i}$ over an infinitely thin surface.

The case might be that you use the following probability density / likelihood:

$$f_\mathbf{y_{ij}}(\mathbf{y_{ij}}\vert \mathbf{\mu_i}, \sigma_f, \sigma_r) = \frac{1}{\left( \sqrt{2 \pi \sigma_f^2} \right)^{n_j}} e^ { \frac{\sum_{j=1}^{n_{j}} (y_{ij}-\mu_i)^2}{2 \sigma_f^2} } \cdot \frac{1}{\left( \sqrt{2 \pi \sigma_r^2} \right)^{n_i}} e^{ \frac{\sum_{i=1}^{n_i} (\mu_i)^2}{2 \sigma_r^2} } $$

but I would say that this is badly defined (the math looks ok, but the interpretation is not). This is not a density that only needs to be integrated over $d y_{ij}$, but also over $d \mathbf{\mu_i}$. You should not turn this into a likelihood function like $\mathcal{L}(\mathbf{\mu_i}, \sigma_f, \sigma_r \vert \mathbf{y_{ij}})$ but instead $\mathcal{L}( \sigma_f, \sigma_r \vert \mathbf{y_{ij}}, \mathbf{\mu_i})$ (yet you do not observe $\mathbf{\mu_i}$).

It is incorrect to impose a relationship between the parameters in the likelihood function and add a corresponding density term to the expression of the likelihood function (that not surprisingly will blow up to infinity, in this way every unobserved variable may be added and becomes an infinite density somewhere, you could also add unobserved unicorns if you like).

  • $\begingroup$ Isn't this the REML solution rather than the ML solution, since we are integrating over the random effects rather than optimizing them? I know that REML is usually the default, but I don't understand how the ML solution exists. $\endgroup$
    – Cliff AB
    Aug 26, 2018 at 16:46
  • $\begingroup$ What do you mean by 'optimizing the random effects' and a difference between REML and ML? Could you spell out those terms. $\endgroup$ Aug 26, 2018 at 17:09
  • $\begingroup$ Note that by considering $f_{Y_{ij}}$, rather than $f_{Y_{ij} | u_i}$ and $f_{\mu_i}$, we are integrating over $u_i$ rather than picking the MLE of $u_i$. I believe this is the REML method, which has advantages over the ML method (such as an upper bounded on the likelihood!). But my question is in regards to the joint likelihood, not the marginal likelihood. $\endgroup$
    – Cliff AB
    Aug 26, 2018 at 18:09
  • $\begingroup$ Could you write the likelihoodfunction $f(y_{ij},\sigma_r,\sigma_f)$ for your case? $\endgroup$ Aug 26, 2018 at 18:25
  • $\begingroup$ @CliffAB No, this is not the reml. The reml log-likelihood objective function is the same as the ml one, but with an extra penalty term: $\log |X^' \Sigma X|$. What you are describing sounds like the hierarchical likelihood. $\endgroup$
    – papgeo
    Aug 27, 2018 at 0:12

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.