Model Fitting Algorithm

Linear model

For each subject i we observe a vector Y_i = (y_i1, …, y_imi)^⊤ ∈ ℝ^m_i and given a design matrix X_i ∈ ℝ^m_i × p and a corresponding coefficient vector β ∈ ℝ^p we assume that the observations are multivariate normal distributed: Y_i ∼ N(X_iβ, Σ_i) where the covariance matrix Σ_i ∈ ℝ^{m_i × m_i} is derived by subsetting the overall covariance matrix Σ ∈ ℝ^m × m appropriately by Σ_i = G_i^−1/2S_i^⊤ΣS_iG_i^−1/2 where the subsetting matrix S_i ∈ {0, 1}^m × m_i contains in each of its m_i columns contains a single 1 indicating which overall time point is matching t_ij. Each row contains at most a single 1 but can also contain only 0 if this time point was not observed. For example, assume a subject was observed on time points 1, 3, 4 out of total 5 then the subsetting matrix is $$ S_i = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix}. $$ G_i ∈ ℝ_> 0^{m_i × m_i} is the diagonal weight matrix, which is the identity matrix if no weights are specified. Note that this follows from the well known property of the multivariate normal distribution that linear combinations of the random vector again have a multivariate normal distribution with the correspondingly modified mean vector and covariance matrix.

Conditional on the design matrices X_i, the coefficient vector β and the covariance matrix Σ we assume that the observations are independent between the subjects.

We can write the linear model for all subjects together as Y = Xβ + ϵ where Y ∈ ℝ^N combines all subject specific observations vectors Y_i such that we have in total $N = \sum_{i = 1}^{n}{m_i}$ observations, X ∈ ℝ^N × p combines all subject specific design matrices and ϵ ∈ ℝ^N has a multivariate normal distribution ϵ ∼ N(0, Ω) where Ω ∈ ℝ^N × N is block-diagonal containing the subject specific Σ_i covariance matrices on the diagonal and 0 in the remaining entries.

Covariance matrix model

The symmetric and positive definite covariance matrix $$ \Sigma = \begin{pmatrix} \sigma_1^2 & \sigma_{12} & \dots & \dots & \sigma_{1m} \\ \sigma_{21} & \sigma_2^2 & \sigma_{23} & \dots & \sigma_{2m}\\ \vdots & & \ddots & & \vdots \\ \vdots & & & \ddots & \vdots \\ \sigma_{m1} & \dots & \dots & \sigma_{m,m-1} & \sigma_m^2 \end{pmatrix} $$ is parametrized by a vector of variance parameters θ = (θ₁, …, θ_k)^⊤. There are many different choices for how to model the covariance matrix and correspondingly θ has different interpretations. Since any covariance matrix has a unique Cholesky factorization Σ = LL^⊤ where L is the lower triangular Cholesky factor, we are going to use this below.

Weighted least squares estimator

Let’s have a look at the details of calculating the log likelihood above, including in particular the weighted least squares (WLS) estimator β̂.

Starting point is the linear equation above and the observation that both the left and right hand sides can be decomposed into subject-specific terms given the block-diagonal structure of Ω and therefore its inverse, W = Ω⁻¹: $$ X^\top \Omega^{-1} X = X^\top W X = \sum_{i=1}^{n} X_i^\top W_i X_i $$ and similarly $$ X^\top \Omega^{-1} Y = X^\top W Y = \sum_{i=1}^{n} X_i^\top W_i Y_i $$ where W_i = Σ_i⁻¹ is the weight matrix for subject i, the inverse of its covariance matrix.

Instead of calculating this inverse explicitly, it is always better numerically to work with the Cholesky factorization and solve linear equations instead. Here we calculate the factorization Σ_i = L_iL_i^⊤. Note that in the case where m_i = m, i.e. this subject has all time points observed, then Σ_i = Σ and we don’t need to calculate this again because we have already Σ = LL^⊤, i.e. L_i = L. Unfortunately, if m_i < m, then we need to calculate this explicitly, as there is no way to update the Cholesky factorization for a subset operation Σ_i = S_i^⊤ΣS_i as we have above. Given L_i, we solve L_iX̃_i = X_i for X̃_i with an efficient forward-solve, and similarly we solve L_iỸ_i = Y_i for Ỹ_i. Therefore we have X_i^⊤W_iX_i = X̃_i^⊤X̃_i and X_i^⊤W_iY_i = X̃_i^⊤Ỹ_i and we can thereby calculate the left and right hand sides for the WLS estimating equation. We solve this equation with a robust Cholesky decomposition with pivoting. The advantage is that we can reuse this decomposition for calculating the covariance matrix of β̂, i.e. K = (X^⊤WX)⁻¹, by supplying the identity matrix as alternative right hand side.

Determinant and quadratic form

For the objective function we also need the log determinant of Ω: where l_i, jj are the diagonal entries of the factor L_i and we have used that

• the determinant of a block diagonal matrix is the product of the determinants of the blocks,
• the determinant of the product of matrices is the product of the determinants,
• the determinant of the transposed matrix is the same as the original one,
• the determinant of a triangular matrix is the product of the diagonal.

And finally, for the quadratic form we can reuse the weighted response vector and design matrix: $$ (Y - X\hat{\beta})^\top \Omega^{-1} (Y - X\hat{\beta}) = \sum_{i=1}^{n} (Y_i - X_i\hat{\beta})^\top W_i (Y_i - X_i\hat{\beta}) = \sum_{i=1}^{n} (\tilde{Y}_i - \tilde{X}_i\hat{\beta})^\top (\tilde{Y}_i - \tilde{X}_i\hat{\beta}) $$

Completing the square

Let’s focus now on the quadratic form in the exponential function and complete the square with regards to β to obtain the kernel of a multivariate normal distribution: where we used K = (X^⊤WX)⁻¹ and could early on identify K as the covariance matrix of the kernel of the multivariate normal of β and then later β̂ as the mean vector.

With this, we know that the integral of the multivariate normal kernel is the inverse of the normalizing constants, and thus $$ \int_{\mathbb{R}^p} \exp\left\{ - \frac{1}{2}(Y - X\beta)^\top \Omega^{-1} (Y - X\beta) \right\} d\beta = \exp\left\{ -\frac{1}{2} Y^\top \Omega^{-1} Y + \frac{1}{2} \hat{\beta}^{-1} K^{-1} \hat{\beta} \right\} (2\pi)^{p/2} \det{K}^{1/2} $$ such that the integrated likelihood is $$ L(\theta; Y) = (2\pi)^{-(N-p)/2} \det(\Omega)^{-1/2} \det{K}^{1/2} \exp\left\{ -\frac{1}{2} Y^\top \Omega^{-1} Y + \frac{1}{2} \hat{\beta}^\top K^{-1} \hat{\beta} \right\}. $$

Objective function

As objective function which we want to minimize with regards to the variance parameters θ we again take the negative natural logarithm $$ f(\theta) = -\log L(\theta;Y) = \frac{N-p}{2} \log(2\pi) + \frac{1}{2}\log\det(\Omega) - \frac{1}{2}\log\det(K) + \frac{1}{2} \tilde{Y}^\top \tilde{Y} - \frac{1}{2} \hat{\beta}^\top \tilde{X}^\top \tilde{X} \hat{\beta} $$ It is interesting to see that computation of the REML objective function is only requiring a few additional calculations compared to the ML objective function. In particular, since we already have the matrix decomposition of K⁻¹, it is very easy to obtain the determinant of it.

Also here we use numeric optimization of f(θ) and the TMB library supports this efficiently through automatic differentiation.