# 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 \(\bm{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

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:

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

where \(\tilde{\bm{b}}\) is the solution to

and \(\tilde{\bm{H}}=\bm{H}(\tilde{\bm{b}})\) is the value of the negative Hessian in \(\bm{b}\)

for \(\bm{b}=\tilde{\bm{b}}\).

Since this quadratic expansion, let us call it \(\ell_{\text{Lapl}}(\bm{y},\bm{b})\), is a (multivariate) quadratic function of \(\bm{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:

### Penalized quasi-likelihood (PQL)¶

If one disregards the dependence of \(\tilde{\bm{H}}\) on \(\bm{\alpha}\) and \(\bm{b}\), then \(\tilde{\bm{b}}\) maximizes not only \(\ell_{\text{cpl}}(\bm{y},\bm{b})\) but also \(\ell^*_{\text{Lapl}}\). This motivates the following IWLS/Fisher scoring equations for \(\hat{\bm{\alpha}}\) and \(\tilde{\bm{b}}\) (see Breslow, and Clayton 1993 and this page):

where

is the IWLS “working dependend variable” with \(\bm{\alpha}\), \(\bm{b}\), \(\bm{W}\), and \(\bm{\pi}\) computed in an earlier iteration.

Substitutions lead to the equations:

and

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

Here

and

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

or the “REML” variant:

This motivates the following algorithm, which is strongly inspired by the
`glmmPQL()`

function in Brian
Ripley’s *R* package MASS:

- Create some suitable starting values for \(\bm{\pi}\), \(\bm{W}\), and \(\bm{y}^*\)
- Construct the “working dependent variable” \(\bm{y}^*\)
- Minimize \(q_1\) (quasi-ML) or \(q_2\) (quasi-REML) iteratively (inner loop), to obtain an estimate of \(\bm{\Sigma}\)
- Obtain \(\hat{\bm{\alpha}}\) and \(\tilde{\bm{b}}\) based on the current estimate of \(\bm{\Sigma}\)
- Compute updated \(\bm{\eta}=\bm{X}\bm{\alpha} + \bm{Z}\bm{b}\), \(\bm{\pi}\), \(\bm{W}\).
- If the change in \(\bm{\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{\bm{\alpha}}\) and \(\tilde{\bm{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 for
this, because the weighting matrix \(\bm{W}\) is non-diagonal. Instead, \(q_1\) or
\(q_2\) are minized using the function `nlminb`

from the standard *R* package
“stats”.

## The Solomon-Cox approximation and MQL¶

### The Solomon-Cox approximation¶

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

where \(\bm{g}_0=\bm{g}(\bm{0})\) is the gradient of \(\ell(\bm{y}|\bm{b};\bm{\alpha})\)

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

Like before, the integral of the exponential this quadratic expansion (which we refer to as \(\ell_{\text{SC}}(\bm{y},\bm{b})\)) has a closed-form solution, as does its logarithm, which is:

### 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 \(\bm{y}^*\). With PQL it is constructed as \(\bm{y}^* = \bm{X}\bm{\alpha} + \bm{Z}\bm{b} + \bm{W}^{-}(\bm{y}-\bm{\pi})\) while the MQL working dependent variable is just

so that the algorithm has the following steps:

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