Issue with REML likelihood--logdet terms cancel

Question

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)

\begin{equation} 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) \end{equation} where the sums trivially simplify to \begin{equation} 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) \end{equation}

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:

\begin{equation} \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) \end{equation}

For the second term, applying Sherman Morrison,

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

Multiplying in $X^T$ and $X$ anda \begin{equation} 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) \\ \end{equation}

\begin{equation} =\|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) \end{equation}

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.

Answer 1

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).