Approximate Inference for Multinomial Logit Models with Random Effects

The problem

A crucial problem for inference about non-linear models with random effects is that the likelihood function for such models involves integrals for which no analytical solution exists.

For given values $\boldsymbol{b}$ of the random effects the likelihood function of a conditional logit model (and therefore also of a baseline-logit model) can be written in the form

$$\mathcal{L}_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b}) = \exp\left(\ell_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b})\right) =\exp \left( \ell(\boldsymbol{y}|\boldsymbol{b};\boldsymbol{\alpha}) -\frac12\ln\det(\boldsymbol{\Sigma}) -\frac12\boldsymbol{b}'\boldsymbol{\Sigma}^{-1}\boldsymbol{b} \right)$$

However, this “complete data” likelihood function cannot be used for inference, because it depends on the unobserved random effects. To arrive at a likelihood function that depends only on observed data, one needs to used the following integrated likelihood function:

$$\mathcal{L}_{\text{obs}}(\boldsymbol{y}) = \int \exp\left(\ell_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b})\right) \partial \boldsymbol{b} = \int \exp \left( \ell(\boldsymbol{y}|\boldsymbol{b};\boldsymbol{\alpha}) -\frac12\ln\det(\boldsymbol{\Sigma}) -\frac12\boldsymbol{b}'\boldsymbol{\Sigma}^{-1}\boldsymbol{b} \right) \partial \boldsymbol{b}$$

In general, this integral cannot be “solved”, i.e. eliminated from the formula by analytic means (it is “analytically untractable”). Instead, one will compute it either using numeric techniques (e.g. using numerical quadrature) or approximate it using analytical techniques. Unless there is only a single level of random effects numerical quadrature can become computationally be demanding, that is, the computation of the (log-)likelihood function and its derivatives can take a lot of time even on modern, state-of-the-art computer hardware. Yet approximations based on analytical techniques hand may lead to biased estimates in particular in samples where the number of observations relative to the number of random offects is small, but at least they are much easier to compute and sometimes making inference possible after all.

The package “mclogit” supports to kinds of analytical approximations, the Laplace approximation and what one may call the Solomon-Cox appoximation. Both approximations are based on a quadratic expansion of the integrand so that the thus modified integral does have a closed-form solution, i.e. is analytically tractable.

The Laplace approximation and PQL

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 $\boldsymbol{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 $\ell^*_{\text{Lapl}}(\boldsymbol{y},\boldsymbol{b})$—is a (multivariate) quadratic function of $\boldsymbol{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 $\boldsymbol{\alpha}$ and $\boldsymbol{b}$, then $\tilde{\boldsymbol{b}}$ maximizes not only $\ell_{\text{cpl}}(\boldsymbol{y},\boldsymbol{b})$ but also $\ell^*_{\text{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

$$\boldsymbol{y}^* = \boldsymbol{X}\boldsymbol{\alpha} + \boldsymbol{Z}\boldsymbol{b} + \boldsymbol{W}^{-}(\boldsymbol{y}-\boldsymbol{\pi})$$

is the IWLS “working dependend variable” with $\boldsymbol{\alpha}$, $\boldsymbol{b}$, $\boldsymbol{W}$, and $\boldsymbol{\pi}$ 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

$$(\boldsymbol{Z}'\boldsymbol{W}\boldsymbol{Z} + \boldsymbol{\Sigma}^{-1})\boldsymbol{b} = \boldsymbol{Z}'\boldsymbol{W}(\boldsymbol{y}^*-\boldsymbol{X}\boldsymbol{\alpha})$$

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

Here

$$\boldsymbol{V} = \boldsymbol{W}^-+\boldsymbol{Z}\boldsymbol{\Sigma}\boldsymbol{Z}'$$

and

$$\boldsymbol{V}^- = \boldsymbol{W}- \boldsymbol{W}\boldsymbol{Z}'\left(\boldsymbol{Z}'\boldsymbol{W}\boldsymbol{Z}+\boldsymbol{\Sigma}^{-1}\right)^{-1}\boldsymbol{Z}\boldsymbol{W}$$

Following Breslow and Clayton (1993) the variance parameters in $\boldsymbol{\Sigma}$ are estimated by minimizing

$$q_1 = \det(\boldsymbol{V})+(\boldsymbol{y}^*-\boldsymbol{X}\boldsymbol{\alpha})\boldsymbol{V}^-(\boldsymbol{y}^*-\boldsymbol{X}\boldsymbol{\alpha})$$

or the “REML” variant:

$$q_2 = \det(\boldsymbol{V})+(\boldsymbol{y}^*-\boldsymbol{X}\boldsymbol{\alpha})\boldsymbol{V}^-(\boldsymbol{y}^*-\boldsymbol{X}\boldsymbol{\alpha})+\det(\boldsymbol{X}'\boldsymbol{V}^{-}\boldsymbol{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 $\boldsymbol{\pi}$, $\boldsymbol{W}$, and $\boldsymbol{y}^*$
2. Construct the “working dependent variable” $\boldsymbol{y}^*$
3. Minimize $q_1$ (quasi-ML) or $q_2$ (quasi-REML) iteratively (inner loop), to obtain an estimate of $\boldsymbol{\Sigma}$
4. Obtain $hat{\boldsymbol{\alpha}}$ and $\tilde{\boldsymbol{b}}$ based on the current estimate of $\boldsymbol{\Sigma}$
5. Compute updated $\boldsymbol{\eta}=\boldsymbol{X}\boldsymbol{\alpha} + \boldsymbol{Z}\boldsymbol{b}$, $\boldsymbol{\pi}$, $\boldsymbol{W}$.
6. If the change in $\boldsymbol{\eta}$ 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 $\boldsymbol{W}$ is non-diagonal. Instead, $q_1$ or $q_2$ are minimized using the function nlminb from the standard R package “stats” or some other optimizer chosen by the user.

The Solomon-Cox approximation and MQL

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 $\boldsymbol{g}\_0=\boldsymbol{g}(\boldsymbol{0})$ is the gradient of $\ell(\boldsymbol{y}\|\boldsymbol{b};\boldsymbol{\alpha})$

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

at $\boldsymbol{b}=\boldsymbol{0}$, while $\boldsymbol{H}\_0=\boldsymbol{H}(\boldsymbol{0})$ is the negative Hessian at $\boldsymbol{b}=\boldsymbol{0}$.

Like before, the integral of the exponential this quadratic expansion (which we refer to as $\ell_{\text{SC}}(\boldsymbol{y},\boldsymbol{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 $\boldsymbol{y}^*$. With PQL it is constructed as $$\boldsymbol{y}^* = \boldsymbol{X}\boldsymbol{\alpha} + \boldsymbol{Z}\boldsymbol{b} + \boldsymbol{W}^{-}(\boldsymbol{y}-\boldsymbol{\pi})$$ while the MQL working dependent variable is just

$$\boldsymbol{y}^* = \boldsymbol{X}\boldsymbol{\alpha} + \boldsymbol{W}^{-}(\boldsymbol{y}-\boldsymbol{\pi})$$

so that the algorithm has the following steps:

1. Create some suitable starting values for $\boldsymbol{\pi}$, $\boldsymbol{W}$, and $\boldsymbol{y}^*$
2. Construct the “working dependent variable” $\boldsymbol{y}^*$
3. Minimize $q_1$ (quasi-ML) or $q_2$ (quasi-REML) iteratively (inner loop), to obtain an estimate of $\boldsymbol{\Sigma}$
4. Obtain $\hat{\boldsymbol{\alpha}}$ based on the current estimate of $\boldsymbol{\Sigma}$
5. Compute updated $\boldsymbol{\eta}=\boldsymbol{X}\boldsymbol{\alpha}$, $\boldsymbol{\pi}$, $\boldsymbol{W}$.
6. If the change in $\boldsymbol{\eta}$ 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}}$.

References

Breslow, Norman E., and David G. Clayton. 1993. “Approximate Inference in Generalized Linear Mixed Models.” Journal of the American Statistical Association 88 (421): 9–25.

Harville, David A. 1997. Matrix Algebra from a Statistician’s Perspective. New York: Springer.

Pinheiro, José C., and Douglas M. Bates. 2000. Mixed-Effects Models in s and s-PLUS. New York: Springer. https://doi.org/10.1007/b98882.

Solomon, P. J., and D. R. Cox. 1992. “Nonlinear Component of Variance Models.” Biometrika 79 (1): 1–11. https://doi.org/10.1093/biomet/79.1.1.

Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with s. Fourth. New York: Springer. https://www.stats.ox.ac.uk/pub/MASS4/.