Finding the MLE for covariance matrix in random effects model

Question

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.

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.