I have a random effects-only linear model as follows:

$$\mathbf{y_i = Z_ib_i + \epsilon_i,\quad i=1,...,N}$$ $$\mathbf{\epsilon_i \sim N(0,\sigma^2I)},\quad \mathbf{b_i \sim N(0,\sigma^2D)}$$

where $\mathbf{I}$ is an identity matrix and $\mathbf{D}$ a covariance matrix of size $n_i \times n_i$.

I am trying to find the MLE for $\mathbf{D}$ using the profile likelihood. For a sample $\mathbf{y_{i= 1,...,N}}$, the full likelihood function up to a constant term should be:

$$\ell(\theta) = -\frac{1}{2}\{(\sum_{i=1}^N n_i)\ln\sigma^2+\sum_{i=1}^N\left(\ln\det\left[\mathbf{I+Z_iDZ_i^T}\right] + \sigma^{-2}\mathbf{y_i^T\left[I+Z_iDZ_i^T\right]^{-1}y_i}\right)\}$$

Taking the profile with respect to $\sigma^2$ we end up with

$$\ell_p(\theta) = -\frac{1}{2}\{(\sum_{i=1}^N n_i)\ln\left[\sum_{i=1}^N y_i^T\left[\mathbf{I+Z_iDZ_i^T}\right]^{-1}y_i\right]+\sum_{i=1}^N\ln\det\left[\mathbf{I+Z_iDZ_i^T}\right]\}$$

Then the task at hand is to find the $\hat\theta=\mathbf{\hat D}$ which maximizes the expression.

My knowledge of matrix calculus is really basic. So far I've been able to find these potentially useful formulas:

  • $d\mathbf{_Ax^TAx} = \mathbf{x^Tx}$

  • $d_p\mathbf{A^{-1}} = \mathbf{-A^{-1}(}d_p\mathbf{A)A^{-1}}$

I assume some form of chain rule can then be applied using the above but I'm having a hard time parsing the whole thing.

Just in case someone's wondering, we were told to find the expression for $\ell_p(\theta)$ as homework but this part is an extra step I wanted to do.

  • 1
    $\begingroup$ Technically, you're not doing Maximum Likelihood, but Minimum Log-Likelihood; so you're not trying to maximize the given expression, but minimizing it. I'd like to understand what you're doing, but there are some missing pieces. What is meant by your model -- please state explicitly what the parameters, what the dependent and what the independent variable is. Also regarding the initial description, is the variance $\sigma^2$ really the same in both terms? And for the likelihood expression, what is $\theta$? I don't recognize any parameter. $\endgroup$
    – cherub
    May 7, 2018 at 12:08
  • $\begingroup$ Is it a typo that $D$ is $n_i \times n_i$? And is it indeed common to all $b_i$? $\endgroup$
    – jld
    May 7, 2018 at 21:29
  • $\begingroup$ @Chaconne you're correct, $D$ is actually $k_i \times k_i$ and may vary over $i$ $\endgroup$ May 8, 2018 at 19:19
  • $\begingroup$ Thanks for the update. Can we assume any structure at all on the $Z_i$ and $D_i$, or are they really completely arbitrary (aside from $D_i$ being positive definite and etc)? $\endgroup$
    – jld
    May 8, 2018 at 19:22
  • $\begingroup$ @cherub $\sigma^2$ is the variance of the $\epsilon_i$ but it's a scaling factor for the variance of the $b_i$. $y$ is the dependent variable, while all $Z_i$ are fixed independent variables. The vector $\theta$ can be seen as a column vector holding the minimum amount of parameters needed to specify the equation, so essentially the diagonals of all the $D_i$ and terms above or below it. $\endgroup$ May 8, 2018 at 19:23

1 Answer 1


So you really need to know three identities in matrix differentiation in order to get first order conditions for this problem. These are:

  1. For vectors $a$ and $b$ that do not depend on $t$ and $X$ a matrix that depends on $t$ we have $\frac{d}{dt}a'X(t)b=a'\frac{dX}{t} b$

  2. For a matrix $X$ that depends on $t$ we have $\frac{d}{dt}X^{-1}=-X^{-1}\frac{dX}{dt} X^{-1}$

  3. For a matrix $X$ that depends on $t$ we have $\frac{d}{dt}ln(det(X))=tr\bigg(X^{-1}\frac{dX}{dt}\bigg)$ Where 'tr' gives the tracee of the matrix.

Using these formulas to take derivatives we see that the score of your log likelihood (that is, the derivative of the profile likelihood w.r.t. $D$) is given by:

$-\bigg(\sum_{i=1}^N n_i\bigg)\bigg(\sum_{i=1}^N y_i ' (I+Z_i D Z_i ')^{-1} [Z_i Z_i'] (I+Z_i D Z_i ')^{-1} y_i \bigg) \times \bigg(\sum_{i=1}^N y_i ' (I+Z_i D Z_i ')^{-1} y_i \bigg)^{-1} - \bigg( \sum_{i=1}^N [I\times trace\big( (I+Z_i D Z_i ')^{-1} Z_i Z_i' \big)]\bigg) $

I would just find the $D$ that sets this equal to zero numerically by coding up the formula in MATLAB say. I don't think there will be a neat analytical solution. Anyway I thought I'd just work that out anyway in case you were curious. You can try applying those rules to take the derivatives yourself and you should get the same thing.


This site is temporarily in read-only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .