45  Methods

1 Estimation Methods

1.1 The Laplacian Method

Let \(\lambda_i(\eta)\) be \(-2 \log l_i(\theta, \eta)\), and let \(g_i(\eta)\) and \(H_i(\eta)\) be the gradient (column) vector and hessian matrix, respectively, of \(\lambda_i\) evaluated at \(\eta\). An approximation to \(-2 \log L_i\) is given by

\[ \begin{aligned} &\hat\lambda_i + \log |\Omega| + \log \left|\Omega^{-1} + \tfrac{1}{2} \hat H_i\right| + \hat\eta_i' \Omega^{-1} \hat\eta_i \\ &\quad - \left(\tfrac{1}{2} \hat g_i + \Omega^{-1} \hat\eta_i\right)' \left(\Omega^{-1} + \tfrac{1}{2} \hat H_i\right)^{-1} \left(\tfrac{1}{2} \hat g_i + \Omega^{-1} \hat\eta_i\right) \end{aligned} \]

where \(\hat\eta_i\) is some estimate of \(\eta_i\), and \(\hat\lambda_i\), \(\hat g_i\), and \(\hat H_i\) are \(\lambda_i\), \(g_i\), and \(H_i\) all evaluated at \(\hat\eta_i\). This results from applying a general approximation approach to integrals, attributable to the French mathematician Laplace, and described by De Bruijn (1961). With \(\hat\eta_i\) equal to the conditional estimate obtained by maximizing the posterior density of \(\eta_i\) (in an unconstrained manner) - call this the unconstrained conditional estimate, this particular approximation has been used by others (Lindley 1980; Mosteller and Wallace 1964), although not with a function \(l_i\) that is as complicated as that which often arises in population pharmacokinetic and pharmacodynamic analyses. See also: Tierney and Kadane (1986). In this particular case, the last term of the approximation is 0. In general, the approximation can produce reasonable results as long the posterior distribution of \(\eta_i\) is dominated by a single mode. On occasion, a randomly dispersed parameter seems to have a multimodal distribution. See the discussion in section B concerning mixture models for a way to address this issue.

Each of the estimation methods uses a different variant of this approximation. However, with whatever variant is used, when in particular, the \(\hat\eta_i\) are taken to be conditional estimates of the \(\eta_i\) at \(\theta\) and \(\Omega\), the general method described in chapter I becomes what we call a conditional estimation method. When the approximation is used just as it is stated above, and when the \(\hat\eta_i\) are taken to be the unconstrained conditional estimates, the method is called the Laplacian estimation method, to honor the individual whose approximation plays such an essential role. However, the method itself involves an idea which is peculiar to NONMEM implementation. Namely, the approximation to \(L\) (the likelihood function of \(\theta\) and \(\Omega\)), resulting from using the Laplacian approximation, is maximized.

When mean-variance models are used, the assumption can be made that each intraindividual variance-covariance matrix \(C_i(\eta_i; \theta)\) is actually given by \(C_i(0; \theta)\), the matrix for the mean individual. With this particular assumption, there is said to be no \(\eta,\varepsilon\)-interaction; see chapter I. The \(\lambda_i\) are computed differently, depending on whether an \(\eta,\varepsilon\)-interaction is assumed, as are the posterior modes. With mean-variance models, by default, NONMEM implements the Laplacian method assuming that there is no \(\eta,\varepsilon\)-interaction. With the currently distributed NONMEM code it is possible to apply the Laplacian method when there is an \(\eta,\varepsilon\)-interaction, but this code and its usage are not supported by the NONMEM Project Group.

1.2 The FOCE Method

The matrix \(H_i\) can be approximated by another matrix. Suppose given \(\eta_i\), \(y_i\) is comprised of statistically independent subvectors \(y_i^{(1)}\), \(y_i^{(2)}\), etc., so that \(\lambda_i\) can be written as a sum over terms \(\lambda_i^{(1)}\), \(\lambda_i^{(2)}\), etc. Then each of \(g_i\) and \(H_i\) can be written as a sum over terms \(g_i^{(1)}\), \(g_i^{(2)}\), etc. and \(H_i^{(1)}\), \(H_i^{(2)}\), etc., respectively. An approximation \(\tilde H_i\) to \(H_i\) is obtained by replacing each \(H_i^{(k)}\) in the sum for \(H_i\) by \(\tilde H_i^{(k)} = \tfrac{1}{2} g_i^{(k)} (g_i^{(k)})'\). This is a type of first-order approximation; terms involving second derivatives have been dropped. It is called the first-order approximation. With this approximation, and when all the \(\hat\eta_i\) are taken to be equal to the unconstrained conditional estimates of the \(\eta_i\), the method is called the first-order conditional estimation (FOCE) method.

Actually, NONMEM allows the implementation of several versions of this method.

• When a mean-variance intraindividual model is used, by default, \(\tilde H_i^{(k)}\) is replaced by \(\tfrac{1}{2} E(g_i^{(k)} (g_i^{(k)})')\), where \(E\) represents the expectation over \(y_i\) under the intraindividual model. With the currently distributed NONMEM code it is possible to use the FOCE method without doing this, but this code and its usage are not supported by the NONMEM Project Group.

• The first-order conditional estimation method without interaction is the FOCE method applied with intraindividual mean-variance models and assuming no \(\eta,\varepsilon\)-interaction. When the intraindividual variance is assumed to be homoscedastic, and moreover, to be the same across individuals, then there is no \(\eta,\varepsilon\)-interaction, and in this case it may be shown that the FOCE method (without interaction) often produces results similar to those obtained with a method described by Lindstrom and Bates (1990). The first-order conditional estimation method with interaction is the FOCE method applied with intraindividual mean-variance models, but without the no interaction assumption. FOCE with and without interaction are both supported. With the currently distributed NONMEM code it is possible to apply the FOCE method with intraindividual models that are not mean-variance models, but this code and its usage are not supported by the NONMEM Project Group.

1.3 The FO Method

When the first-order approximation is used (with \(H_i^{(k)}\) replaced by \(\tfrac{1}{2} E(g_i^{(k)} (g_i^{(k)})')\)), but when all \(\hat\eta_i\) are taken to be 0 (the population mean value of \(\eta\)), the method is called the first-order (FO) estimation method. With the first-order method, the terms \(\hat\eta_i' \Omega^{-1} \hat\eta_i\) and \(\Omega^{-1} \hat\eta_i\) in the Laplacian approximation are 0. Note that since conditional estimates are not used, the first-order method is not a conditional estimation method.

It can be shown that when intraindividual mean-variance models are used, the method is equivalent to the first-order method as described, for example, in NONMEM Users Guide - Part I (also see e.g., Beal and Sheiner (1985)). Such an earlier description is also given below in section A.6. These earlier descriptions of the method apply only to mean-variance models. With the currently distributed NONMEM code it is possible to apply the FO method as defined above with intraindividual models that are not mean-variance models, but this usage is not recommended, and the code is not supported by the NONMEM Project Group.

1.4 The Hybrid Method

Suppose certain (but not all) elements of \(\eta\) are chosen to be in a set \(\Gamma\), that the elements of \(\hat\eta_i\) corresponding to the elements of \(\eta\) in \(\Gamma\) are taken to be 0, and that the remaining elements of \(\hat\eta_i\) are taken to be those given by the Bayes posterior mode of \(\eta_i\) under the restriction that all elements of \(\eta\) in \(\Gamma\) are 0. The conditional estimate thus defined is an example of a constrained conditional estimate. Suppose also that the first-order approximation is made. Then the method is a hybrid between the first-order method and the FOCE method. Accordingly, this conditional estimation method is called the hybrid method. Note that with the definition of the \(\hat\eta_i\) used with this method, in contrast with the definition used with the FOCE and Laplacian methods, the last term in the Laplacian approximation is not 0.

A hybrid method can be considered that uses a weaker version of the first-order approximation. Consider using the first-order approximation, but only for the submatrix of \(H_i\) consisting of just those partial second derivatives such that the two variables with respect to which the differentiation occurs are in \(\Gamma\). This method is not supported with the currently distributed NONMEM code.

When the intraindividual models are statistical linear models (linear in the parameters \(\eta_i\)), the first-order, first-order conditional, hybrid, and Laplacian methods are all the same method, the classical maximum likelihood method.

1.5 The Centering Methods

The \(\eta_i\) are assumed to be distributed in the population with mean 0. When the population model fits the data well, this will be reflected by the average, \(\bar\eta\), of the conditional estimates of the \(\eta_i\) across the sampled individuals (at the values of the population parameters given by the model) being close to 0. (The converse does not necessarily hold.) When \(\bar\eta\) is close to 0, the fit will be called centered. There is nothing about the methods defined above that insures that the fit will be centered. There are infrequently arising situations where the average is “far” from 0, where the model does not fit well (as judged e.g. by the differences \(y_{ij} - f_{ij}(0, \hat\theta)\) with mean-variance intraindividual models) and where a method that is designed to better center the fit might be tried (do see chapter III for some guidance). With a centering estimation method, the \(\hat\eta_i\) are taken to be the unconstrained conditional estimates, and the approximation to \(-2 \log L_i\) is given by

\[ -2 \log l_i(\theta, \hat\eta_i - \bar\eta) + \log |\Omega| + \log \left|\Omega^{-1} + \tfrac{1}{2} \hat H_i\right| + (\hat\eta_i - \bar\eta)' \Omega^{-1} (\hat\eta_i - \bar\eta) \]

With NONMEM, there are centering FOCE and Laplacian estimation methods (with no \(\eta,\varepsilon\)-interaction). A centering hybrid method is not implemented in NONMEM.

1.6 The Centering FOCE Method with the First-Order Model

The first-order model is the population model which results when for all \(i\), the \(i\)th given intraindividual model is a mean-variance model with mean \(E_i(\eta_i; \theta)\) and variance-covariance matrix \(C_i(\eta_i; \theta)\), and this model is replaced by another such model with mean

\[E_i(0; \theta) + \frac{\partial E_i}{\partial \eta_i}(0; \theta) \eta_i\]

and variance-covariance matrix \(C_i(0; \theta)\).

The linearity of the \(\eta_i\) under this model implies that the population expectation of \(y_{ij}\) is \(f_{ij}(0; \theta)\), the prediction obtained by taking \(\eta_i\) to be 0, its population mean. With mean-variance models, the FO estimation method is sometimes described as the application of the maximum likelihood method to the first-order model that results from the given model, and when using this method, it is usual to judge goodness of fit by the differences \(y_{ij} - f_{ij}(0; \hat\theta)\). When a conditional estimation method is used instead of the FO method, a centered fit may result, confirming that the population mean of the \(\eta_i\) is 0. However, the given intraindividual models are used, and they may be nonlinear in the \(\eta_i\). Therefore, conceivably, \(f_{ij}(0; \theta)\) may be a poor approximation to the population expectation of \(y_{ij}\), and for this reason alone, an apparent bias in the fit may result. Experience suggests, though, that this should not be a major concern (perhaps because the nonlinear effect is small relative to the size of intraindividual variability in the residuals). If one is concerned, there are a couple of strategies one might use.

First, the NONMEM program allows the expectation of the \(y_{ij}\) to be estimated by means of a couple different types of actual integration (and not just when the intraindividual models are of mean-variance kind); see NONMEM Users Guide - Part VIII. Second, when the intraindividual models are mean-variance models, NONMEM allows the first-order model to be obtained automatically from the given model and used with the centering FOCE method. (If the first-order model is used with the noncentering FOCE method, the result is the same as that obtained with the FO method.) When a conditional estimation method is needed (see chapter III), application of the centering FOCE method to the first-order model that results from the given model may yield adequate results, and of course, the expectation of \(y_{ij}\) under the first-order model is simply given by \(f_{ij}(0; \theta)\). Moreover, due to the linearity of the intraindividual models (of the first-order model) in the \(\eta_i\), the computational requirement is substantially less than that incurred with application of the (centering or noncentering) FOCE method to the given model. The savings in CPU time is achieved at the expense of possibly using too simple a model (and, of course is still not as great a savings as is achieved with the FO method).

The first-order model may be used with the centering FOCE method, but not with the centering Laplacian method (because due to the linearity, the result would be the same as that obtained with the centering FOCE method). Be aware that when this model is used with the centering FOCE method, the conditional estimates produced by the method are based on the first-order intraindividual models (unlike whenever the noncentering FOCE method is used, where the conditional estimates are based on the given intraindividual models). It is possible nonetheless to obtain posthoc estimates based on the given intraindividual models, at the population estimates obtained from using the centering FOCE method with the first-order model. A centering hybrid method is not implemented in NONMEM.

2 Mixture Models

On occasion, a model may need to incorporate a randomly dispersed parameter that has a possibly multimodal distribution. In this case a mixture model may be useful. This is a model where for each \(i\), there are several possible intraindividual models, \(M_1, M_2, \ldots, M_r\) for \(y_i\), and it is assumed that the particular model that actually describes \(y_i\) is one of these, but it is not known which one. It is assumed that the probability that it is \(M_k\) is \(p_k\), where \(p_1 + p_2 + \cdots + p_r = 1\). Loosely put, the \(i\)th individual is chosen randomly from a population divided into \(r\) subpopulations, their relative sizes either being known or unknown. The subpopulation of which the individual is a given member is not observable, but for each subpopulation, a model for data from an individual from the subpopulation is available. The mixing probabilities \(p_k\) correspond to the sizes of the subpopulations and are usually treated as parameters whose values are unknown and are estimated. With NONMEM, these probabilities can be modeled, i.e. related to covariables, and therefore, can vary between individuals. The parameters of these relationships can be estimated; they are included in \(\theta\). To indicate this generality, the \(p_k\) may be written \(p_{ik}(\theta)\) (the \(k\)th mixing probability for the \(i\)th individual).

Suppose, for example, that a clearance parameter of a pharmacokinetic model may be bimodally distributed in the population. Here is how this may be expressed with a population model. One may consider a mixture model with two intraindividual models for each individual: for the \(i\)th individual, one where the individual’s clearance is given by

\[CL_i = \theta_1 \exp(\eta_{i1}) \tag{1}\]

and another where it is given by

\[CL_i = \theta_2 \exp(\eta_{i2}) \tag{2}\]

(The parameters \(\eta_{i1}\) and \(\eta_{i2}\) are the first two elements of \(\eta_i\).) For each \(i\), the value \(\eta_i\) arises randomly (see chapter I). For each \(i\), a choice between the two intraindividual models is also viewed as one being made in a random fashion, according to probabilities \(p_1\) and \(p_2\) (\(p_1 + p_2 = 1\)). As a result of this choice, a value \(\eta_i^*\), which is either \(\eta_{i1}\) or \(\eta_{i2}\), is also “chosen”. (Consequently, if after \(\eta_{i1}\), say, is chosen, the value of \(\eta_{i2}\) does not influence the data.) From the point of view of not knowing what choices between intraindividual models were actually made, the distribution of the \(\eta_i^*\) across individuals is a mixture of two normal distributions, and the distribution of the \(CL_i\) is a mixture of two lognormal distributions.

The first two elements of the random variable \(\eta\) may have the same or different variances, i.e. \(\omega_{11}\) may or may not equal \(\omega_{22}\). If these variances are sufficiently small, while the parameters \(\theta_1\) and \(\theta_2\) are sufficiently far apart, and if both probabilities \(p_1\) and \(p_2\) are sufficiently large (however in this regard, the variances, the \(\theta\)’s, and the probabilities must actually be considered altogether), the distribution of \(CL_i\) is bimodal. Often, the data may not allow all of the different variances between mixture components, such as \(\omega_{11}\) and \(\omega_{22}\), to be well estimated, in which case the assumption might be made that these variances are the same (a homoscedastic assumption). With NONMEM, this can be done explicitly, or alternatively, the “same \(\eta\)” can be used with both mixture components, e.g. \(\eta_{i1}\) can be used in (1) and also in (2), instead of \(\eta_{i2}\). NONMEM will understand that \(\eta_{i1}\) is symbolizing two “different \(\eta\)’s”, each having the same variance.¹

Other examples of mixture models may be given. See NONMEM Users Guide - Part VI, section III.L.2 for an example where the mixture model describes a mixture of two joint lognormal distributions for clearance and volume, but which is not a bimodal distribution. The differences between the models \(M_k\) need not be differences concerning parameters; they could be differences in model form. They can be any set of differences whatsoever.

The likelihood for \(y_i\) under a mixture model is

\[L_i(\theta, \Omega) = \sum_{k=1}^{r_i} p_{ik}(\theta) L_{ik}(\theta, \Omega)\]

where \(L_{ik}\) is the likelihood function for \(y_i\) under the \(k\)th possible intraindividual model for individual \(i\). With a mixture model, any of the estimation methods described in section A uses the defining approximation for the method with each of the \(L_{ik}\), \(k = 1, \ldots, r\).

With a set of values for the population parameters \(\theta\) and \(\Omega\), NONMEM classifies each individual into one of the \(r\) subpopulations. The classification gives the most probable subpopulation of which the individual is a member. For each \(k\), the empirical Bayes (marginal) posterior probability that \(y_i\) is described by \(M_k\), given \(y_i\), is computed by \([p_{ik}(\theta) L_{ik}(\theta, \Omega)] / L_i(\theta, \Omega)\). The individual is classified into the \(k\)th subpopulation if the \(k\)th probability is the largest among these \(r\) values.

Beal, Stuart L., and Lewis B. Sheiner. 1985. “Methodology of Population Pharmacokinetics.” Chap. 4 in Drug Fate and Metabolism: Methods and Techniques, Vol. 5, edited by Edward R. Garrett and Jean L. Hirtz. Marcel Dekker.

De Bruijn, N. G. 1961. Asymptotic Methods in Analysis. North-Holland.

Lindley, D. V. 1980. “Approximate Bayesian Methods.” In Bayesian Statistics, edited by J. M. Bernardo, M. H. DeGroot, D. V. Lindley, and A. F. M. Smith. University Press.

Lindstrom, Mary J., and Douglas M. Bates. 1990. “Nonlinear Mixed Effects Models for Repeated Measures Data.” Biometrics 46 (3): 673–87. https://doi.org/10.2307/2532087.

Mosteller, Frederick, and David L. Wallace. 1964. Inference and Disputed Authorship: The Federalist. Addison-Wesley.

Tierney, Luke, and Joseph B. Kadane. 1986. “Accurate Approximations for Posterior Moments and Marginal Densities.” Journal of the American Statistical Association 81 (393): 82–86. https://doi.org/10.1080/01621459.1986.10478240.

---

1. With NONMEM Version IV, the same \(\eta\) can also be used, and NONMEM will understand that it is symbolizing two different \(\eta\)’s with the same variance, provided the first-order estimation method is used.