104 FOCE and Laplace Methods
NONMEM 7 Technical Guide
Generally the integral of the joint density Equation 10 is very difficult to evaluate deterministically, but it may be approximated for classical methods FOCE and Laplace using a method described by Beal (Part VII of the NONMEM Users Guide). The derivation is given in (Wang 2007), while we will merely report the results. Classical NONMEM (FO, FOCE, and Laplace) does not require the use of \(\boldsymbol{\mu}\) modeling, so for this section, we will use the parameterization of \(\boldsymbol{\theta}, \boldsymbol{\eta}\), for all parameters, rather than distinguishing between \(\boldsymbol{\mu}\) modeled and non-\(\boldsymbol{\mu}\) modeled \(\boldsymbol{\theta}\). For example, the individual’s joint density may be alternately expressed as
\[ \begin{aligned} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i(\boldsymbol{\theta}_\mu), \boldsymbol{\Omega}) &= l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega}) \\ &= l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) \\ &= p(\mathbf{y}_i, \boldsymbol{\eta} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) \end{aligned} \tag{1}\]
and integration over all values of \(\boldsymbol{\eta}\) is equivalent to integrating over all values of \(\boldsymbol{\phi}\):
\[ L_i = -\log\left(\int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\eta} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) d\boldsymbol{\eta}\right) \tag{2}\]
In order to integrate the individual’s joint density over all \(\boldsymbol{\eta}\) using the approximation suggested by Beal, we wish first to determine the set of \(\boldsymbol{\eta}\) at the maximum of this joint density, or equivalently, at the minimum of the negative logarithm of the joint density:
\[ -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega})\bigr) = -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \frac{1}{2} \boldsymbol{\eta}^{\prime} \boldsymbol{\Omega}^{-1} \boldsymbol{\eta} \tag{3}\]
We minimize with respect to eta by evaluating
\[ \begin{aligned} \frac{ \partial -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\eta} } &= \frac{ \partial -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) }{ \partial \boldsymbol{\eta} } + \boldsymbol{\Omega}^{-1} \boldsymbol{\eta} \\ &= \mathbf{0} \end{aligned} \tag{4}\]
using typical search strategies. The \(\boldsymbol{\eta}\) at which equation 4 is satisfied is called the mode of the joint density for subject \(i\), and shall be designated \(\hat{\boldsymbol{\eta}}_i\) (the hat over the parameter shall refer to a mode or point estimate, whereas the line over a parameter refers to a mean). Finding the \(\hat{\boldsymbol{\eta}}_i\) parameters that provide the minimum of the individual’s joint density is called mode a posteriori (MAP) estimation. This is used to then evaluate an approximation of the negative logarithm of the individual’s integral of his joint density as follows:
\[ \begin{aligned} L_i &= -\log\left(\int_{-\infty}^{\infty} l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) d\boldsymbol{\eta}\right) \\ &\approx -\log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta})\bigr) + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \frac{1}{2} \hat{\boldsymbol{\eta}}_i^{\prime} \boldsymbol{\Omega}^{-1} \hat{\boldsymbol{\eta}}_i \\ &\quad + \frac{1}{2}\log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta}))\bigr) = L_{Ni} \end{aligned} \tag{5}\]
where \(\mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta})\) is the hessian or information matrix to the data density \(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta})\). The total approximate objective function is in turn, for \(m\) subjects:
\[ L_N = \sum_{i=1}^{m} L_{Ni} \tag{6}\]
(the subscript \(N\) refers to classical NONMEM) where the first three terms of Equation 5 are simply the negative logarithm of the joint density evaluated at the mode \(\hat{\boldsymbol{\eta}}_i\), and the last term is ½ of the negative logarithm of the determinant of the variance of \(\boldsymbol{\eta}\) under the joint density, \(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega})\). One can therefore think of the joint density evaluated at the mode (that is, the first three terms of equation 5) as the “height” of the joint density:
\[ -\log(H_i) = -\log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta})\bigr) + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \frac{1}{2} \hat{\boldsymbol{\eta}}_i^{\prime} \boldsymbol{\Omega}^{-1} \hat{\boldsymbol{\eta}}_i \tag{7}\]
where \(H_i\) is the “height” of individual \(i\)’s joint density. Similarly, one may think of one-half the determinant of the variance under the joint density (the last term of the equation 5) as its “width”:
\[ -\log(W_i) = \frac{1}{2}\log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta}))\bigr) \tag{8}\]
Thus equation 5 represents the negative logarithm of the height multiplied by the width, resulting in the “area” of the joint density. The \(\mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta})\) may be evaluated several ways. One method is to evaluate the second derivative, usually by finite difference methods:
\[ \mathbf{S}^{-1}(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) = \left\{\frac{ \partial^2 -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) }{ \partial \boldsymbol{\eta}_{k_1} \partial \boldsymbol{\eta}_{k_2} }, k_1 = 1\ \text{to}\ n, k_2 = 1\ \text{to}\ n\right\} \tag{9}\]
where \(\{\}\) means “matrix containing elements”. This evaluation is used in the Laplace method in NONMEM. This evaluation guarantees positive definiteness (assuming no numerical difficulties arise) when evaluated at the mode (see appendix B). Another method is by the cross product of the first derivatives of the individual data point densities:
\[ \mathbf{S}^{-1}(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) = \left\{\sum_{j=1}^{m_i} \frac{ \partial -\log\bigl(l(y_{ij} \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) }{ \partial \boldsymbol{\eta}_{k_1} } \frac{ \partial -\log\bigl(l(y_{ij} \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) }{ \partial \boldsymbol{\eta}_{k_2} }\right\} \tag{10}\]
where \(m_i\) is the number of data points for patient \(i\) (assuming all data are independent). Based on its structure, positive definiteness is guaranteed even when evaluated not at the mode (see appendix B). A third method for evaluating \(\mathbf{S}^{-1}(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta})\) is by the expected value of the second derivative:
\[ \begin{aligned} & \mathbf{S}^{-1}(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) \\ &= \left\{E\left[\frac{ \partial^2 -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) }{ \partial \boldsymbol{\eta}_{k_1} \partial \boldsymbol{\eta}_{k_2} }\right]\right\} \\ &= \left\{\frac{ \partial \mathbf{f}(t_{ij}, \boldsymbol{\eta}, \boldsymbol{\theta})^{\prime} }{ \partial \boldsymbol{\eta}_{k_1} } \mathbf{V}_i^{-1} \frac{ \partial \mathbf{f}(t_{ij}, \boldsymbol{\eta}, \boldsymbol{\theta}) }{ \partial \boldsymbol{\eta}_{k_2} } + \frac{1}{2}\operatorname{tr}\left(\mathbf{V}_i^{-1} \frac{ \partial \mathbf{V}_i }{ \partial \boldsymbol{\eta}_{k_1} } \mathbf{V}_i^{-1} \frac{ \partial \mathbf{V}_i }{ \partial \boldsymbol{\eta}_{k_2} }\right)\right\} \end{aligned} \tag{11}\]
which is also positive definite even when not evaluated at the mode. Equation 11 is used as the non-Laplace (CONDITIONAL) method in NONMEM. The \(\boldsymbol{\Omega}^{-1} + \mathbf{S}_i^{-1}\) is in turn the Hessian (information) matrix of the joint density:
\[ \begin{aligned} & \left\{E\left[\frac{ \partial^2 -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\eta}_{k_1} \partial \boldsymbol{\eta}_{k_2} }\right]\right\} \\ &= \left\{E\left[\frac{ \partial^2 -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta})\bigr) }{ \partial \boldsymbol{\eta}_{k_1} \partial \boldsymbol{\eta}_{k_2} }\right] + E\left[\frac{ \partial^2 -\log\bigl(h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\eta}_{k_1} \partial \boldsymbol{\eta}_{k_2} }\right]\right\} \\ &= \mathbf{S}^{-1}(\mathbf{y}_i \mid \boldsymbol{\eta}, \boldsymbol{\theta}) + \boldsymbol{\Omega}^{-1} \end{aligned} \tag{12}\]
and hence its inverse is the variance matrix of \(\boldsymbol{\eta}\) under the joint density, as mentioned earlier.
Because the sum of two positive definite matrices is itself positive definite, the variance of the joint density as evaluated above is positive definite. For joint densities that are exactly multivariate normally distributed with respect to \(\boldsymbol{\eta}\), equation 5 evaluates the joint area exactly. We shall also refer to \(\mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta})\) as \(\hat{\mathbf{S}}_i^{-1}\). The \(\hat{\mathbf{S}}_i^{-1}\) must be evaluated at the individual \(i\)’s mode of his joint density, at \(\hat{\boldsymbol{\eta}}_i\), and not at the mean population position of \(\boldsymbol{\eta} = \mathbf{0}\), so the INTERACTION option in NONMEM must be used.
Keep in mind that while the \(\hat{\boldsymbol{\eta}}_i\) represents the \(i\)th individual’s “best fit” parameters for its data, based on its joint density, it is only needed here to evaluate the area under his joint density using the above approximation. In other words, we really don’t need an individual’s best fit parameter set theoretically, but we need it practically, in order to evaluate the “height” of the density, and thus approximate his joint density area. There are alternative methods of finding the area without needing to know the individual’s “best fit” parameters, which we will explore later.
Following the evaluation of each individual’s objective function in the manner described above, these are summed to form the total approximate objective function \(L_N\). NONMEM optimizes \(L_N\) with respect to THETAS, OMEGAS, and SIGMAS using a variable metric method, in which \(L_N\) is evaluated at a series of values of \(\boldsymbol{\theta}\) and \(\boldsymbol{\Omega}\), to provide a directional search to find the set of \(\boldsymbol{\theta}\) and \(\boldsymbol{\Omega}\) that optimizes \(L_N\). The description of the variable metric method is beyond the scope of this document, but a good reference is (Dennis and Schnabel 1996).