# Issue with REML likelihood--logdet terms cancel

I'm trying to write an implementation of a linear mixed effects model using REML.

I'm working with a simple model:

$$y_{ij} = X_{ij}\beta + Z_{ij}b_i + \epsilon_{ij}$$ In my case, the covariate $$X$$ is a scalar, all groups are same size, and $$X=Z$$. Thus, we can assume that $$X_{ij} = X_j = Z_{ij}$$. Assume that $$\epsilon_{ij}\stackrel{iid}{\sim}N(0,\sigma^2)$$ where $$\sigma^2$$ is known, and $$b_i\sim N(0,\gamma)$$. There is no crossover, so the negative log-likelihood (up to a constant)

$$$$L(\gamma,\beta)=\sum_{i=1}^m\log\det(\sigma^2 I+\gamma XX^T) + \sum_{i=1}^m\log\det\left(X^T(\sigma^2 I+\gamma XX^T)^{-1}X)\right) + \\ +\sum_{i=1}^m (y_i - X\beta)^T(\sigma^2 I+\gamma XX^T)^{-1}(y_i - X\beta)$$$$ where the sums trivially simplify to $$$$L(\gamma,\beta)=m\log\det(\sigma^2 I+\gamma XX^T) + m\log\det\left(X^T(\sigma^2 I+\gamma XX^T)^{-1}X)\right) + \\ +\sum_{i=1}^m (y_i - X\beta)^T(\sigma^2 I+\gamma XX^T)^{-1}(y_i - X\beta)$$$$

However, the first two terms cancel to leave a constant, so that the only remaining non-constant term becomes the third, which is minimized at $$\gamma = \infty$$. This can be seen in the following calculation:

$$$$\det\left(I+\gamma XX^{T}\right)=1+\gamma\|X\|^{2} \\ \det\left(X^{T}\left(I+\gamma XX^{T}\right)^{-1}X\right) \\ =\det\left(X^{T}\left(\frac{X}{1+\gamma\|X\|^{2}}\right)\right)$$$$

For the second term, applying Sherman Morrison,

$$$$\left(I+\alpha XX^{T}\right)^{-1}=I-\frac{\alpha XX^{T}}{1+\alpha\|X\|^{2}} \\$$$$

Multiplying in $$X^T$$ and $$X$$ anda $$$$X^{T}\left(I+\gamma XX^{T}\right)^{-1}X=\|X\|^{2}-\frac{\gamma\|X\|^{4}}{1+\gamma\|X\|^{2}}=\|X\|^{2}\left(1-\frac{\gamma\|X\|^{2}}{1+\gamma\|X\|^{2}}\right) \\$$$$

$$$$=\|X\|^{2}\left(\frac{1+\gamma\|X\|^{2}-\gamma\|X\|^{2}}{1+\gamma\|X\|^{2}}\right)=\|X\|^{2}\left(\frac{1}{1+\gamma\|X\|^{2}}\right)$$$$

Upon taking logs, one finds that the terms which are nonconstant in $$\gamma$$ cancel.

I don't know how to reconcile this in my code--I originally happened upon this because my estimate for $$\gamma$$ was going to infinity. However, when I run the exact model in python statsmodels mixed effects, I get perfectly reasonable results, so it seems I am missing something about the REML objective.

• Is there no intercept in the fixed effect part of the model? Commented Jul 16 at 19:26
• And I don't see how the two first terms simplify to something not depending on $\gamma$. Commented Jul 16 at 19:31
• If you look at the it below, and take the log of each, the gamma terms will cancel. However I figured out the issue, I messed up a log-det
– Alex
Commented Jul 17 at 4:01

The issue is that the second term does not reduce to $$\sum_{i=1}^m \log\det(X^T(\sigma^2I + \gamma ZZ^T)X$$ as we can't separate the logdet because of the $$X^T$$ and $$X$$ sandwiching.
It should be $$\log\det (\sum_{i=1}^m(X^T(\sigma^2I + \gamma ZZ^T)X$$))--this makes sense as a degree of freedom correction, since the one parameter we're fitting should only kill one dof, and thus there should only be one log in correction (heuristically).