Laplace approximation

The (first-order) Laplace approximation is based on the quadratic expansion the logarithm of the integrand, the complete-data log-likelihood

$$ \ell_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b})\approx \ell(\boldsymbol{y}|\tilde{\boldsymbol{b}};\boldsymbol{\alpha}) - \frac12 (\boldsymbol{b}-\tilde{\boldsymbol{b}})' \tilde{\boldsymbol{H}} (\boldsymbol{b}-\tilde{\boldsymbol{b}}) -\frac12\ln\det(\boldsymbol{\Sigma}) -\frac12(\boldsymbol{b}-\tilde{\boldsymbol{b}})'\boldsymbol{\Sigma}^{-1}(\boldsymbol{b}-\tilde{\boldsymbol{b}}) $$

where $\tilde{\boldsymbol{b}}$ is the solution to

$$ \frac{\partial\ell_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b})}{\partial\boldsymbol{b}} = 0 $$

and $\tilde{\boldsymbol{H}}=\boldsymbol{H}(\tilde{\boldsymbol{b}})$ is the value of the negative Hessian with respect to b

$$ \boldsymbol{H}(\boldsymbol{b})=-\frac{\partial^2\ell(\boldsymbol{y}|\boldsymbol{b};\boldsymbol{\alpha})}{\partial\boldsymbol{b}\partial\boldsymbol{b}'} $$

for $\boldsymbol{b}=\tilde{\boldsymbol{b}}$.

Since this quadratic expansion—let us call it ℓ_Lapl^*(y, b)—is a (multivariate) quadratic function of b, the integral of its exponential does have a closed-form solution (the relevant formula can be found in Harville (1997)).

For purposes of estimation, the resulting approximate log-likelihood is more useful:

$$ \ell^*_{\text{Lapl}} = \ln\int \exp(\ell_{\text{Lapl}}(\boldsymbol{y},\boldsymbol{b})) \partial\boldsymbol{b} = \ell(\boldsymbol{y}|\tilde{\boldsymbol{b}};\boldsymbol{\alpha}) - \frac12\tilde{\boldsymbol{b}}'\boldsymbol{\Sigma}^{-1}\tilde{\boldsymbol{b}} - \frac12\ln\det(\boldsymbol{\Sigma}) - \frac12\ln\det\left(\tilde{\boldsymbol{H}}+\boldsymbol{\Sigma}^{-1}\right). $$

Penalized quasi-likelihood (PQL)

If one disregards the dependence of $\tilde{\boldsymbol{H}}$ on α and b, then $\tilde{\boldsymbol{b}}$ maximizes not only ℓ_cpl(y, b) but also ℓ_Lapl^*. This motivates the following IWLS/Fisher scoring equations for $\hat{\boldsymbol{\alpha}}$ and $\tilde{\boldsymbol{b}}$ (see Breslow and Clayton (1993) and this page):

$$ \begin{aligned} \begin{bmatrix} \boldsymbol{X}'\boldsymbol{W}\boldsymbol{X} & \boldsymbol{X}'\boldsymbol{W}\boldsymbol{Z} \\ \boldsymbol{Z}'\boldsymbol{W}\boldsymbol{X} & \boldsymbol{Z}'\boldsymbol{W}\boldsymbol{Z} + \boldsymbol{\Sigma}^{-1}\\ \end{bmatrix} \begin{bmatrix} \hat{\boldsymbol{\alpha}}\\ \tilde{\boldsymbol{b}}\\ \end{bmatrix} = \begin{bmatrix} \boldsymbol{X}'\boldsymbol{W}\boldsymbol{y}^*\\ \boldsymbol{Z}'\boldsymbol{W}\boldsymbol{y}^* \end{bmatrix} \end{aligned} $$

where

y^* = Xα + Zb + W⁻(y − π)

is the IWLS “working dependend variable” with α, b, W, and π computed in an earlier iteration.

Substitutions lead to the equations:

$$ (\boldsymbol{X}\boldsymbol{V}^-\boldsymbol{X})\hat{\boldsymbol{\alpha}} = \boldsymbol{X}\boldsymbol{V}^-\boldsymbol{y}^* $$

and

(Z′WZ + Σ⁻¹)b = Z′W(y^* − Xα)

which can be solved to compute $\hat{\boldsymbol{\alpha}}$ and $\tilde{\boldsymbol{b}}$ (for given Σ)

Here

V = W⁻ + ZΣZ′

and

V⁻ = W − WZ′(Z′WZ + Σ⁻¹)⁻¹ZW

Following Breslow and Clayton (1993) the variance parameters in Σ are estimated by minimizing

q₁ = det (V) + (y^* − Xα)V⁻(y^* − Xα)

or the “REML” variant:

q₂ = det (V) + (y^* − Xα)V⁻(y^* − Xα) + det (X′V⁻X)

This motivates the following algorithm, which is strongly inspired by the glmmPQL() function in Brian Ripley’s R package MASS (Venables and Ripley 2002):

1. Create some suitable starting values for π, W, and y^*
2. Construct the “working dependent variable” y^*
3. Minimize q₁ (quasi-ML) or q₂ (quasi-REML) iteratively (inner loop), to obtain an estimate of Σ
4. Obtain hatα and $\tilde{\boldsymbol{b}}$ based on the current estimate of Σ
5. Compute updated η = Xα + Zb, π, W.
6. If the change in η is smaller than a given tolerance criterion stop the algorighm and declare it as converged. Otherwise go back to step 2 with the updated values of $\hat{\boldsymbol{\alpha}}$ and $\tilde{\boldsymbol{b}}$.

This algorithm is a modification of the IWLS algorithm used to fit conditional logit models without random effects. Instead of just solving a linear requatoin in step 3, it estimates a weighted linear mixed-effects model. In contrast to glmmPQL() it does not use the lme() function from package nlme (Pinheiro and Bates 2000) for this, because the weighting matrix W is non-diagonal. Instead, q₁ or q₂ are minimized using the function nlminb from the standard R package “stats” or some other optimizer chosen by the user.

The Solomon-Cox approximation

The (first-order) Solomon approximation (Solomon and Cox 1992) is based on the quadratic expansion the integrand

$$ \ell_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b})\approx \ell(\boldsymbol{y}|\boldsymbol{0};\boldsymbol{\alpha}) + \boldsymbol{g}_0' \boldsymbol{b} - \frac12 \boldsymbol{b}' \boldsymbol{H}_0 \boldsymbol{b} -\frac12\ln\det(\boldsymbol{\Sigma}) -\frac12\boldsymbol{b}'\boldsymbol{\Sigma}^{-1}\boldsymbol{b} $$

where g_0 = g(0) is the gradient of ℓ(y∥b; α)

$$ \boldsymbol{g}(\boldsymbol{b})=-\frac{\partial\ell(\boldsymbol{y}|\boldsymbol{b};\boldsymbol{\alpha})}{\partial\boldsymbol{b}} $$

at b = 0, while H_0 = H(0) is the negative Hessian at b = 0.

Like before, the integral of the exponential this quadratic expansion (which we refer to as ℓ_SC(y, b)) has a closed-form solution, as does its logarithm, which is:

$$ \ln\int \exp(\ell_{\text{SC}}(\boldsymbol{y},\boldsymbol{b})) \partial\boldsymbol{b} = \ell(\boldsymbol{y}|\boldsymbol{0};\boldsymbol{\alpha}) - \frac12\boldsymbol{g}_0'\left(\boldsymbol{H}_0+\boldsymbol{\Sigma}^{-1}\right)^{-1}\boldsymbol{g}_0 - \frac12\ln\det(\boldsymbol{\Sigma}) - \frac12\ln\det\left(\boldsymbol{H}_0+\boldsymbol{\Sigma}^{-1}\right). $$

Marginal quasi-likelhood (MQL)

The resulting estimation technique is very similar to PQL (again, see Breslow and Clayton 1993 for a discussion). The only difference is the construction of the “working dependent” variable y^*. With PQL it is constructed as y^* = Xα + Zb + W⁻(y − π) while the MQL working dependent variable is just

y^* = Xα + W⁻(y − π)

so that the algorithm has the following steps:

1. Create some suitable starting values for π, W, and y^*
2. Construct the “working dependent variable” y^*
3. Minimize q₁ (quasi-ML) or q₂ (quasi-REML) iteratively (inner loop), to obtain an estimate of Σ
4. Obtain $\hat{\boldsymbol{\alpha}}$ based on the current estimate of Σ
5. Compute updated η = Xα, π, W.
6. If the change in η is smaller than a given tolerance criterion stop the algorighm and declare it as converged. Otherwise go back to step 2 with the updated values of $\hat{\boldsymbol{\alpha}}$.