103  Basic Theory of Nonlinear Mixed Effects

NONMEM 7 Technical Guide

Individual parameters phi (\(\boldsymbol{\phi}\)) to a PK/PD model are assumed to have a random distribution in a population of subjects, typically a normal distribution with mean MU (\(\boldsymbol{\mu}\)) and variance OMEGA (\(\boldsymbol{\Omega}\)). The mean \(\boldsymbol{\mu}\) may in turn be modeled as a function of a set of unknown but to be estimated fixed effects parameters THETA (\(\boldsymbol{\theta}\)), and a set of covariates, or information about individual \(i\), \(\mathbf{x}_i\). The deviation of the individual parameter \(\boldsymbol{\phi}\) from its mean is designated \(\boldsymbol{\eta}\), so that the following relation holds:

\[ \boldsymbol{\phi} = \boldsymbol{\mu}_i(\boldsymbol{\theta}_\mu, \mathbf{x}_i) + \boldsymbol{\eta} \tag{1}\]

where \(\boldsymbol{\theta}_\mu\) are those thetas that are related to etas through a MU function, of the above format. Thus, the distribution of \(\boldsymbol{\phi}\) can be described as

\[ h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega}) \propto \frac{1}{\det(\boldsymbol{\Omega})^{1/2}} \exp\left(-\frac{1}{2}(\boldsymbol{\phi} - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi} - \boldsymbol{\mu}_i)\right) \tag{2}\]

The population parameter density \(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\) is the probability that the particular \(\boldsymbol{\phi}\) would occur for an individual, given mean population parameters \(\boldsymbol{\mu}_i\) and its inter-individual covariance \(\boldsymbol{\Omega}\). The distribution of \(\boldsymbol{\eta}\) is therefore centered about zero (\(\mathbf{0}\)), and can be described as

\[ h(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) \propto \frac{1}{\det(\boldsymbol{\Omega})^{1/2}} \exp\left(-\frac{1}{2} \boldsymbol{\eta}^{\prime} \boldsymbol{\Omega}^{-1} \boldsymbol{\eta}\right) \tag{3}\]

Not all fixed effects theta are involved in an eta (\(\boldsymbol{\eta}\)) relationship as shown above. For those theta that are not exclusively expressed in the PK/PD model or error model via MU (\(\boldsymbol{\mu}\)), these are considered non-MU-modeled. We shall designate these thetas as \(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\). The entire vector of thetas is then

\[ \boldsymbol{\theta} = \{\boldsymbol{\theta}_\mu, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\} \tag{4}\]

The parameters as designated in NONMEM are THETA for \(\boldsymbol{\theta}\), ETA for \(\boldsymbol{\eta}\), and OMEGA for \(\boldsymbol{\Omega}\). There are also a set of parameters designated as SIGMA in NONMEM, which are never MU-modeled, and because in our discussion they will be treated in exactly the same way as non MU-modeled theta, we shall include them in \(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\) to reduce the complexity of the nomenclature. There may also be some etas that are not related to a MU-model, in which case

\[ \boldsymbol{\phi} = \boldsymbol{\eta} \tag{5}\]

Thus the phi vector includes all etas, whether or not they are involved in a MU function.

For observed data that are modeled as normally distributed, a predictive function may be evaluated using the individual PK/PD parameters phi, and/or may be modeled directly from fixed effects parameters not be phi/MU-modeled, \(\mathbf{f}_i(\boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\). In addition, a residual variance matrix \(\mathbf{V}\) describes the uncertainty of the observed values, and may be directly a function of the predicted value \(\mathbf{f}_i\), sigma parameters and other non-MU-modeled thetas, and rarely, individual parameters \(\boldsymbol{\phi}\): \(\mathbf{V}_i(\mathbf{f}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\). The normal data density can be expressed as

\[ l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) \propto \frac{ \exp\left(-\dfrac{1}{2}(\mathbf{y}_i - \mathbf{f}_i)^{\prime} \mathbf{V}_i^{-1}(\mathbf{y}_i - \mathbf{f}_i)\right) }{ \sqrt{\det(\mathbf{V}_i)} } \tag{6}\]

where \(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\) is the individual data density, the probability of data \(\mathbf{y}_i\) occurring for individual \(i\), given individual PK/PD parameters \(\boldsymbol{\phi}\), and fixed effect parameters \(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\) that are not MU-modeled.

The joint density of data \(\mathbf{y}_i\) and \(\boldsymbol{\phi}\) for an individual is then

\[ 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{\theta}_\mu), \boldsymbol{\Omega}) \tag{7}\]

The \(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\) is the joint likelihood density of \(\boldsymbol{\phi}\) and \(\mathbf{y}_i\) for a given individual. It is integrated over all possible values of \(\boldsymbol{\phi}\) for each individual, so that the “best” population parameters \(\boldsymbol{\theta}\) and \(\boldsymbol{\Omega}\) are determined by taking into account the joint probability to an individual’s data over the entire parameter space of \(\boldsymbol{\phi}\), rather than at just one particular location, such as at the individual’s best fit. We are therefore interested in evaluating the marginal density of \(\mathbf{y}_i\) for any given \(\boldsymbol{\theta}\) and \(\boldsymbol{\Omega}\) (or \(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}\)):

\[ \begin{aligned} p(\mathbf{y}_i \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) &= \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \\ &= \int_{-\infty}^{\infty} l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \end{aligned} \tag{8}\]

for each subject \(i\). The total marginal density for all \(m\) subjects is then

\[ p(\mathbf{y} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) = \prod_{i=1}^{m} \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{9}\]

It is convenient at this stage to use the negative logarithm of the density, and refer to this as the objective function, for each individual:

\[ L_i = -\log\left(\int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) \tag{10}\]

and for the total data set:

\[ \begin{aligned} L &= -\log\bigl(p(\mathbf{y} \mid \boldsymbol{\theta}, \boldsymbol{\Omega})\bigr) \\ &= \sum_{i=1}^{m} L_i \end{aligned} \tag{11}\]

Thus, the negative logarithm of the parameter density is

\[ -\log\bigl(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) = \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \frac{1}{2}(\boldsymbol{\phi} - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi} - \boldsymbol{\mu}_i) \tag{12}\]

And the negative logarithm of the data density is:

\[ -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) = \frac{1}{2}(\mathbf{y}_i - \mathbf{f}_i)^{\prime} \mathbf{V}_i^{-1}(\mathbf{y}_i - \mathbf{f}_i) + \frac{1}{2}\log\bigl(\det(\mathbf{V}_i)\bigr) \tag{13}\]

To fit a model with mean population parameters \(\boldsymbol{\theta}\) and population variance \(\boldsymbol{\Omega}\) to data \(\mathbf{y}\), the marginal density Equation 9 is to be maximized with respect to \(\boldsymbol{\theta}\) and \(\boldsymbol{\Omega}\). These parameter values are then considered the most likely for the observed data \(\mathbf{y}\). Therefore, the maximization of the marginal density with respect to \(\boldsymbol{\theta}\) and \(\boldsymbol{\Omega}\) is called the maximum likelihood method. In practice, as an equivalent process, the negative logarithm of the marginal density Equation 11 is minimized. This is the goal of the first order (FO), first order conditional estimation (FOCE/FOCEI), Laplace, iterative two stage (ITS), and expectation maximization (EM) methods.