44 Introduction
This document gives a brief description of the estimation methods for population type data that can be used with NONMEM Version V.1 These include, in particular, a few methods that are new with this version, the centered and hybrid methods. The more important changes from the earlier edition published in 1992, but not all changes, are highlighted with the use of vertical bars in the right margin. This document contains no information about how to communicate with the NONMEM program.
To read this document it may be helpful to have some familiarity with the notation used with the representation of statistical models for the NONMEM program. See discussions of models in NONMEM Users Guide - Part I, but if one’s interest is only in using NONMEM with PREDPP, see discussions of models in NONMEM Users Guides - Parts V and VI. Particular notation used in this Guide VII is given next.
The \(j\)th observation from the \(i\)th individual is denoted \(y_{ij}\). Each individual may have a different number of observations. Each observation may be measured on a different scale: continuous, categorical, ordered categorical, discrete-ordinal. An individual can have multivariate observations, each of different lengths. However, the multivariate nature of an observation is suppressed, as this is not relevant to the descriptions given in this document, and so the separate (scalar-valued) observations comprising the multivariate observations are all separately indexed by \(j\). Each multivariate observation may have a different length. The vector of all the observations from the \(i\)th individual is denoted \(y_i\).
It is assumed that there exists a separate statistical model for each \(y_i\). This model is called the intraindividual model, or the individual model for the \(i\)th individual. It is parameterized by \(\theta\), a (vector-valued) parameter common to all the separate intraindividual models, and \(\eta_i\), a (vector-valued) parameter specific to the intraindividual model for \(y_i\). Under this model, the likelihood of \(\eta_i\) for the data \(y_i\) (conditional on \(\theta\)) is denoted by \(l_i(\eta_i; \theta)\), the dependence on \(y_i\) being suppressed in the notation. This likelihood is called here the conditional likelihood of \(\eta_i\).
When all the elements of \(y_i\) are measured on a continuous scale, an often-used intraindividual model is given by the multivariate normal model with mean \(E_i(\eta_i; \theta)\) and variance-covariance matrix \(C_i(\eta_i; \theta)\) (usually, \(\theta\) is comprised of parameters \(\beta\) which are the only ones affecting \(E_i\), and other parameters which, along with \(\beta\), affect \(C_i\)). This type of model shall be referred to as the mean-variance model.2 It is usually expressed in terms of a multivariate normal vector \(\varepsilon\) with mean 0 and variance-covariance matrix \(\Sigma\). In the notation used here, the parameter \(\theta\) includes \(\Sigma\) (ignoring the matrix structure of \(\Sigma\)). For example,
\[y_{ij} = f_{ij}(\eta_i; \theta) + f_{ij}(\eta_i; \theta) \varepsilon_{ij}\]
where \(\varepsilon_{ij}\) is an instance of a univariate normal variable \(\varepsilon\) with variance \(\Sigma = \sigma^2\). (When \(\varepsilon\) is multivariate, the observation \(y_{ij}\) is modeled in terms of a single instance of this multivariate random vector. A few other observations as well may be modeled in terms of this same instance, and thus under the model, all such observations are correlated and comprise a multivariate observation.) In this example, \(E_{ij}(\eta_i; \theta)\) is \(f_{ij}(\eta_i; \theta)\) (the mean of \(y_{ij}\)), and \(C_{ijj}(\eta_i; \theta)\) is \(f_{ij}(\eta_i; \theta)^2 \sigma^2\) (the variance of \(y_{ij}\)). Since the ratio of the standard deviation of \(y_{ij}\) to the mean of \(y_{ij}\) is the constant \(\sigma\), this particular model is called the constant coefficient of variation model.
The dependence of \(C_i\) on \(\eta_i\) is often a consequence of the intraindividual variance depending on the mean function, as with the above example, which in turn depends on \(\eta_i\). This dependence represents an interaction between \(\eta_i\) and \(\varepsilon\). With the (homoscedastic) model expressed by
\[y_{ij} = f_{ij}(\eta_i; \theta) + \varepsilon_{ij}\]
there is no such interaction; \(C_{ijj}(\eta_i; \theta)\) is just \(\sigma^2\). There are two variants of the first-order conditional estimation method described in chapter II, one that takes this interaction into account and another that ignores it.
When an intraindividual model involving \(\varepsilon\) is presented to NM-TRAN (the “front-end” of the NONMEM system), the model is automatically transformed. A linearization of the right side of the equation is used: a first-order approximation in \(\varepsilon_{ij}\) about 0, the mean value of \(\varepsilon_{ij}\). Since the approximate model is linear in \(\varepsilon_{ij}\), it is a mean-variance model. Clearly, if the given model is itself a mean-variance model, the transformed model is identical to the given model. Consider, for example, an intraindividual model where the elements of \(y_i\) are regarded as lognormally distributed (because the normally distributed \(\varepsilon_{ij}\) appear as logarithms):
\[y_{ij} = f_{ij}(\eta_i; \theta) \exp(\varepsilon_{ij})\]
In this case the transformed model is the constant cv model given above. (Therefore, no matter whether the given intraindividual model or the constant cv model is presented to NM-TRAN, the results of the analysis will be the same.)
Alternatively, the user might be able to transform the data so that a mean-variance model applies to the transformed data, which can then be presented directly to NM-TRAN. With the above example, and using the log transformation on the data \(y_{ij}\), an appropriate mean-variance model to present to NM-TRAN would be
\[y_{ij} = \log f_{ij}(\eta_i; \theta) + \varepsilon_{ij}\]
(Actually, NM-TRAN allows one to explicitly accomplish the log transformation of both the data and the \(f_{ij}\).) The results of the analysis differ depending on whether or not the log transformation is used. Without the log transformation, the values of the \(f_{ij}\) are regarded as arithmetic means (under the approximate model obtained by linearizing), and with the log transformation, these values are regarded as geometric means. Use of the log transformation (when this can be done; when there are no \(y_{ij}\) or \(f_{ij}\) with value 0) can often lead to a better analysis.
It is also assumed that as individuals are sampled randomly from the population, the \(\eta_i\) are also being sampled randomly (and statistically independently), although these values are not observable. The value \(\eta_i\) is called the random interindividual effect for \(y_i\). It is assumed that the \(\eta_i\) are instances of the random vector \(\eta\), normally distributed with mean 0 and variance-covariance matrix \(\Omega\). The density function of this distribution (at \(\eta\)) is denoted by \(h(\eta; \Omega)\).
Often, some quantity \(P\) (viewed as a function of values of the covariates and the \(\eta_i\)) is common to different intraindividual models. For example, a clearance parameter may be common to different intraindividual models, but its value differs between different intraindividual models because the values of the covariates and the \(\eta_i\) differ. The randomness of the \(\eta_i\) in the population induces randomness in \(P\). The quantity \(P\) is said to be a randomly dispersed parameter. When speaking of its distribution, we are imagining that the values of the covariates are fixed, so that indeed, there is a unique distribution in question.
From the above assumptions, the (marginal) likelihood of \(\theta\) and \(\Omega\) for the data \(y_i\) is given by
\[L_i(\theta, \Omega) = \int l_i(\eta; \theta) h(\eta; \Omega) d\eta \tag{1}\]
In general, this integral is difficult to compute exactly. The likelihood for all the data is given by
\[L(\theta, \Omega) = \prod_i L_i(\theta, \Omega) \tag{2}\]
The first-order estimation method was the first population estimation method available with NONMEM. This method produces estimates of the population parameters \(\theta\) and \(\Omega\), but it does not produce estimates of the random interindividual effects. An estimate of \(\eta_i\) is nonetheless obtainable, conditional on the first-order estimates for \(\theta\) and \(\Omega\) (or on any other values for these parameters), by maximizing the empirical Bayes posterior density of \(\eta_i\), given \(y_i\): \([l_i(\eta; \theta) h(\eta; \Omega)] / L_i(\theta, \Omega)\), with respect to \(\eta\). In other words, the estimate is the mode of the posterior distribution. Since this estimate is obtained after values for \(\theta\) and \(\Omega\) are obtained, it is called the posthoc estimate. When a mean-variance model is used, and a request is put to NONMEM to compute a posthoc estimate, by default this estimate is computed using \(C_i(\theta, 0)\). In other words, the intraindividual variance-covariance is assumed to be the same as that for the mean individual, the hypothetical individual having the mean interindividual effect, 0, and sharing the same values of the covariates as has the \(i\)th individual. However, it is also possible to obtain the posterior mode without this assumption.
The posterior density can be maximized using any given values for \(\theta\) and \(\Omega\). Since the resulting estimate for \(\eta_i\) is obtained conditionally on these values, it is sometimes called a conditional estimate at these values, to emphasize its conditional nature.
In contrast with the first-order method, the conditional estimation methods to be described produce estimates of the population parameters and, simultaneously, estimates of the random interindividual effects. With each different method, a different approximation to the likelihood function (1) is used, and (2) is maximized with respect to \(\theta\) and \(\Omega\). The approximation to (1) at the values \(\theta\) and \(\Omega\) depends on an estimate \(\hat\eta_i\), and as this estimate itself depends on the values \(\theta\) and \(\Omega\), the approximation gives rise to a further dependence of \(L_i\) on the values of \(\theta\) and \(\Omega\), one not expressed in (1). Consequently, as different values \(\theta\) and \(\Omega\) are tried, different estimates \(\hat\eta_i\) are obtained as a part of the maximization of (2). The estimates \(\hat\eta_i\) at the values \(\theta\) and \(\Omega\) that maximize (2) constitute the estimates of the random interindividual effects produced by the method (except for the hybrid method).3 The estimate \(\hat\eta_i\) also depends on \(y_i\), and so, the approximation gives rise to a further dependence of \(L_i\) on \(y_i\), one also not expressed in (1).
In contrast with the first-order method, a conditional estimation method involves multiple maximizations within a maximization. The estimate \(\hat\eta_i\) is the value of \(\eta_i\) that maximizes the posterior distribution of \(\eta_i\) given \(y_i\) (except for the hybrid method).4 For each different value of \(\theta\) and \(\Omega\) that is tried by the maximization algorithm used to maximize (2), first a value \(\hat\eta_1\) is found that maximizes the posterior distribution given \(y_1\), then a value \(\hat\eta_2\) is found that maximizes the posterior distribution given \(y_2\), etc. Therefore, maximizing (2) is a very difficult and CPU intensive task. The numerical methods by which this is accomplished are not described in this document.
Fortunately, it often suffices to use the first-order method; a conditional estimation method is not needed, or if it is, sometimes it is needed only minimally during the course of a data analysis. Some guidance is given in chapter III. Briefly, the need for a conditional estimation method increases with the degree to which the intraindividual models are nonlinear in the \(\eta_i\). Population pharmacokinetic models are often actually rather linear in this respect, although the degree of nonlinearity increases with the degree of multiple dosing. Population pharmacodynamic models are more nonlinear. The potential for a conditional estimation method to produce results different from those obtained with the first-order estimation method decreases as the amount of data per individual decreases, since a conditional estimation method uses conditional estimates of the \(\eta_i\), which are all shrunken to 0, and the shrinkage is greater the less the amount of data per individual. Many population analyses involve little amounts of data per individual.
The conditional estimation methods that are available with NONMEM and which are described in chapter II are: the first-order conditional estimation method (with and without interaction when mean-variance models are used, and with or without centering), the Laplacian method (with and without centering), and the hybrid method (a hybrid between the first-order and first-order conditional estimation methods). For purposes of description here and in other NONMEM Users Guides, the term conditional estimation methods refers not only to these population estimation methods, but also to methods for obtaining conditional estimates themselves.
To summarize, each of the (population) conditional estimation methods involves maximizing (2), but each uses a different approximation to (1). Actually, \(-2 \log L\) is minimized with respect to \(\theta\) and \(\Omega\). This is called the objective function. Its minimum value serves as a useful statistic for comparing models. Standard errors for the estimates (indeed, an estimated asymptotic variance-covariance matrix for all the estimates) is obtained by computing derivatives of the objective function.
This document provides a description of estimation methods that can be used with observations of the same or different type. However, essentially, it neither contains any specific information about how to analyze observations of particular types, nor any information about how to communicate with NONMEM in order to do this.↩︎
Here and elsewhere in this section an explicit assumption concerning the normal probability distribution is made. This is done primarily to keep the discussion simple. To various degrees in different situations the normality assumption does not play as important a role as our formally making the assumption might indicate.↩︎
After obtaining the population parameter estimates with the hybrid method (see chapter II), NONMEM ignores the estimates of the \(\eta_i\) that have been produced simultaneously with the population parameter estimates, and as with the first-order method, the posthoc estimates (described above) are the ones reported as the estimates of the random interindividual effects.↩︎
With the hybrid method, a constrained maximum is computed.↩︎