106  Evaluating the Expectation Step: Importance Sampling (IMP)

NONMEM 7 Technical Guide

One can evaluate the area under the joint density and the other integrals by Monte Carlo techniques. The advantage to these methods is that the actual mathematical expression of the integral is not necessary for its computation, and the precision to which the integral is evaluated depends on the number of random samples generated to evaluate the integral. One Monte Carlo method is to use a sampling function that approximates the joint density, from which one obtains sample values of \(\boldsymbol{\eta}\) or \(\boldsymbol{\phi}\).

One possible sampling function is the multivariate normal density that has mean at \(\hat{\boldsymbol{\phi}}_i\) and variance of \((\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1}\) as described in the previous section. To get these values, therefore, one first maximizes the joint density with respect to \(\boldsymbol{\phi}\) (or \(\boldsymbol{\eta}\)) as one would for a MAP estimation. The negative logarithm of the area of this sampling function is exactly \(L_{Ni}\) of equation 5. Thus, the purpose of the randomization method is to modify \(L_{Ni}\) to the extent that the joint density deviates from this sampling density. In practice, one may start with a sampling function that is somewhat larger, by multiplying the variance by a value \(\gamma \ge 1\): \(\gamma(\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1}\). The area of this sampling function is then

\[ \begin{aligned} E_i &= -\log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i, \boldsymbol{\theta})\bigr) + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \\ &\quad \frac{1}{2}(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) + \frac{1}{2}\log\bigl(\det(\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})\bigr) - \frac{1}{2} n\log(\gamma) \end{aligned} \tag{1}\]

where \(n\) is the number of \(\boldsymbol{\phi}\) parameters to be integrated, since \(\boldsymbol{\phi}\) is integrated to form \(L_i\).

For the \(k\)th random sample selected from this sampling density, the parameter vectors \(\boldsymbol{\phi}_{k,i}\) are used to evaluate the logarithm of the joint density at that position:

\[ \log(\pi_i(\boldsymbol{\phi}_{k,i})) = \log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}_{k,i}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\boldsymbol{\phi}_{k,i} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) \tag{2}\]

To evaluate the normalized log of the joint density, we subtract

\[ \log(\pi_i(\hat{\boldsymbol{\phi}}_i)) = \log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\hat{\boldsymbol{\phi}}_i \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) \tag{3}\]

To obtain

\[ \log(\pi_i(\boldsymbol{\phi}_{k,i})) - \log(\pi_i(\hat{\boldsymbol{\phi}}_i)) \tag{4}\]

so that this normalized log joint density is 0 at the mode \(\hat{\boldsymbol{\phi}}_i\). We also evaluate the logarithm of the normalized sampling function (which is also equal to 0 at \(\hat{\boldsymbol{\phi}}_i\)):

\[ \log(e_i(\boldsymbol{\phi}_{k,i})) = -\frac{1}{2}(\boldsymbol{\phi}_{k,i} - \hat{\boldsymbol{\phi}}_i)^{\prime} (\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})(\boldsymbol{\phi}_{k,i} - \hat{\boldsymbol{\phi}}_i)/\gamma \tag{5}\]

The logarithm of the ratio between joint density and sampling density is then:

\[ q_{i,k,y} = \log(\pi_i(\boldsymbol{\phi}_{k,i})) - \log(\pi_i(\hat{\boldsymbol{\phi}}_i)) - \log(e_i(\boldsymbol{\phi}_{k,i})) \tag{6}\]

which evaluates to

\[ \begin{aligned} q_{i,k,y} & = \log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}_{k,i}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) - \log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) \\ &\quad + \frac{1}{2}(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) - \frac{1}{2}(\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_i) + \\ &\quad \frac{1}{2}(\boldsymbol{\phi}_{k,i} - \hat{\boldsymbol{\phi}}_i)^{\prime} (\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})(\boldsymbol{\phi}_{k,i} - \hat{\boldsymbol{\phi}}_i)/\gamma \end{aligned} \tag{7}\]

Its exponent is the probability of accepting this position by the joint density, relative to the sampling density, which we may consider as a weight:

\[ u_{i,k,y} = \exp(q_{i,k,y}) \tag{8}\]

Thus the following fraction,

\[ \psi_i = \frac{1}{r} \sum_{k=1}^{r} u_{i,k,y} \tag{9}\]

represents the ratio of the area of the conditional density to the area of the sampling density.

The \(r\)=ISAMPLE is the number of random samples selected for each individual. This fraction is now used to adjust the area of the sampling density \(E_i\), which is known, to obtain the true area of the conditional density, which is unknown:

\[ \begin{aligned} L_i &= E_i - \log(\psi_i) = \\ &\quad -\log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \frac{1}{2}(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) + \frac{1}{2}\log\bigl(\det(\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})\bigr) \\ &\quad -\log(\psi_i \gamma^{n/2}) \end{aligned} \tag{10}\]

so that \(-\log(\psi_i \gamma^{n/2})\) is the “correction factor” that the randomization method adds to our original area equation to improve its accuracy. In NONMEM, \(\gamma\) is continually adjusted so that \(\psi_i\) approximates IACCEPT, up to the limit of the boundaries of \(\gamma\) being between ISCALE_MIN and ISCALE_MAX (available in NONMEM 7.2).

The above derivation of sample weights and likelihood evaluation for importance sampling resulting in equations 7 and 10 was developed to demonstrate that they are based on general principles of obtaining integrals by Monte Carlo methods. These relationships can be simplified by moving all of the elements from \(E_i\) to \(q_{i,k,y}\), given that the components in \(E_i\) are constant for all random samples \(k\), so that the use of \(\exp(q_{i,k,y})\) as a weight factor will not be affected. Furthermore, we may generalize for all sampling densities \(\boldsymbol{\phi}_{k,i} \sim N(\bar{\boldsymbol{\mu}}_{p,i}, \gamma \boldsymbol{\Omega}_{p,i})\), by substituting a general sampling density mean \(\boldsymbol{\mu}_{p,i}\) in place of \(\hat{\boldsymbol{\phi}}_i\), and a general sampling density variance \(\gamma\boldsymbol{\Omega}_{p,i}\) in place of \(\gamma(\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1}\), so that we obtain:

\[ \begin{aligned} q_{i,k,y} & = \log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}_{k,i}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) - \frac{1}{2}(\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_i) - \frac{1}{2}\log(\det(\boldsymbol{\Omega})) \\ &\quad + \frac{1}{2}(\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_{p,i})^{\prime} (\gamma \boldsymbol{\Omega}_{p,i})^{-1}(\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_{p,i}) + \frac{1}{2}\log(\det(\gamma \boldsymbol{\Omega}_{p,i})) \end{aligned} \tag{11}\]

\[ u_{i,k,y} = \exp(q_{i,k,y}) \tag{12}\]

\[ \psi_i = \frac{1}{r} \sum_{k=1}^{r} u_{i,k,y} \tag{13}\]

\[ L_i = -\log(\psi_i) \tag{14}\]

The new and improved objective function is then

\[ L = \sum_{i=1}^{m} L_i \tag{15}\]

With this technique, we can also evaluate an improved mean \(\boldsymbol{\theta}\) and improved variance-covariance matrix. Letting

\[ z_{i,k,y} = \frac{u_{i,k,y}}{\sum_{k=1}^{r} u_{i,k,y}} \tag{16}\]

so that

\[ \sum_{k=1}^{r} z_{i,k,y} = 1 \tag{17}\]

then

\[ \bar{\boldsymbol{\phi}}_i = \sum_{k=1}^{r} z_{i,k,y} \boldsymbol{\phi}_{k,i} \tag{18}\]

\[ \bar{\mathbf{B}}_i = \sum_{k=1}^{r} z_{i,k,y} (\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_i)(\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_i)^{\prime} \tag{19}\]

(note also that these are means and variance about the means, as indicated by the line above the parameter). The update equations yield the following:

\[ \begin{aligned} \frac{\partial L}{\partial \boldsymbol{\theta}_\mu} &= \sum_{i=1}^{m} \frac{\partial L_i}{\partial \boldsymbol{\mu}_i} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \\ & = -\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \\ &= \mathbf{g}_{\boldsymbol{\theta}_\mu} \\ &= \mathbf{0} \end{aligned} \tag{20}\]

\[ \mathbf{H}_{\boldsymbol{\theta}_\mu} = \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta} } \boldsymbol{\Omega}^{-1} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta} } \tag{21}\]

The easiest way to maximize is to perform the following updates:

\[ \hat{\boldsymbol{\theta}}_\mu = \boldsymbol{\theta}_\mu + \alpha \mathbf{H}_{\boldsymbol{\theta}_\mu}^{-1} \mathbf{g}_{\boldsymbol{\theta}_\mu} \tag{22}\]

\(\hat{\boldsymbol{\mu}}_i = \boldsymbol{\mu}_i(\hat{\boldsymbol{\theta}})\)

And, according to equations 41, 39 and 32,

\[ \begin{aligned} \hat{\boldsymbol{\Omega}} &= \frac{1}{m} \sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}}_i)(\bar{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}}_i)^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \bar{\mathbf{B}}_i = \\ &\quad \frac{1}{m} \sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}}_i)(\bar{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}}_i)^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{r} z_{i,k,y} (\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_i)(\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_i)^{\prime} = \\ &\quad \frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{r} z_{i,k,y} (\boldsymbol{\phi}_{k,i} - \hat{\boldsymbol{\mu}}_i)(\boldsymbol{\phi}_{k,i} - \hat{\boldsymbol{\mu}}_i)^{\prime} \end{aligned} \tag{23}\]

This is equivalent to performing summary statistics on all of the random samples among all of the individuals. Note that the normalized weights \(z_{i,k,y}\) defined in equation 16 are obtained from sampled evaluations under the joint density \(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\), and are therefore empirical evaluations of the conditional density of equation 7. As the number of samples approaches infinity (\(r \to \infty\)), equations 22 and 23 approach the exact evaluation of updates that are required, as expressed in equations 14 and 41.

For subsequent iterations, the Monte Carlo evaluated conditional mean and variances of the previous iteration for that subject may be used as the parameters to the sampling density.

This multivariate density has mean as \(\bar{\boldsymbol{\phi}}_{p,i}\) and, and variance of \(\bar{\mathbf{B}}_{p,i}\), so we sample from \(\boldsymbol{\phi}_{k,i} \sim N(\bar{\boldsymbol{\phi}}_{p,i}, \gamma\bar{\mathbf{B}}_{p,i})\) where subscript \(p\) refers to previous iteration, so the pertinent weighting function is:

\[ \begin{aligned} q_{i,k,y} & = \log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}_{k,i}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) - \frac{1}{2}(\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi}_{k,i} - \boldsymbol{\mu}_i) - \frac{1}{2}\log(\det(\boldsymbol{\Omega})) \\ &\quad + \frac{1}{2}(\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_{p,i})^{\prime} (\gamma \bar{\mathbf{B}}_{p,i})^{-1}(\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_{p,i}) + \frac{1}{2}\log(\det(\gamma \bar{\mathbf{B}}_{p,i})) \end{aligned} \tag{24}\]

Followed by

\[ u_{i,k,y} = \exp(q_{i,k,y}) \tag{25}\]

\[ \psi_i = \frac{1}{r} \sum_{k=1}^{r} u_{i,k,y} \tag{26}\]

\[ L_i = -\log(\psi_i) \tag{27}\]

and the additional computations are carried out as before. Whether the parameters to the proposal density are obtained from a MAP estimation, or from conditional means and variances determined from a previous iteration, depends on whether METHOD=IMP or METHOD=IMPMAP is used, and the settings of MAPITER and MAPINTER (available in NONMEM 7.2).

The average Neffective is reported in the NONMEM iteration output (reported as Neff.), and is:

Neffective average =

\[ \frac{r}{m} \sum_{i=1}^{m} F_i \psi_i \tag{28}\]

\[ F_i = \pi_i(\boldsymbol{\mu}_i) / \max(\pi_i(\boldsymbol{\phi})) \tag{29}\]

where \(\max(\pi_i(\boldsymbol{\phi}))\) is the maximal joint density value found among the \(k=1\) to \(r\) samples for subject \(i\), and \(\pi_i(\boldsymbol{\mu}_i)\) is the density at the sample’s mean (\(\boldsymbol{\mu}_i\)) for subject \(i\). When the sampling density is centered at the mode of the joint density, as is the case with IMPMAP mode, \(F\) is set equal to 1. Otherwise, when the sampling density is centered at the conditional mean (\(\boldsymbol{\mu}_i\)) of

the previous iteration, then the \(F\) correction must be made for the fact that in general, the conditional mean will not be lined up with the conditional mode. Notice that if the joint density is identical at all sample positions to the sampling density, then Neffective=ISAMPLE=\(r\), since all summand terms are equal to 1.

The Smpl. is the total samples, usually ISAMPLE=\(r\).

The average relative fitness of the sampling density (fit.) to the joint density, reported in the iteration output of the NONMEM report file, is as follows:

Fitness average =

\[ 1 - \frac{0.5}{m} \sum_{i=1}^{m} F_i\left[\frac{1}{r} \sum_{k=1}^{r} \left|\frac{ \pi_i(\boldsymbol{\phi}_{k,i}) }{ \pi_i(\boldsymbol{\mu}_i) } - e_i(\boldsymbol{\phi}_{k,i})\right|\right] \tag{30}\]

where \(e_i(\boldsymbol{\phi}_{k,i})\) is the normalized sampling density at position \(\boldsymbol{\phi}_{k,i}\) for sample \(k\) of subject \(i\). Normalized here means that \(e_i(\boldsymbol{\mu}_{r,i}) = 1\). Notice that if the sampling and joint densities are identical, then the difference between them is always 0, and the Average fitness will be 1.

For those \(\boldsymbol{\theta}\) that are not expressed in the model through \(\boldsymbol{\mu}\), the \(\boldsymbol{\theta}\) may appear anywhere in the likelihood. To optimize the population objective function in these \(\boldsymbol{\theta}\) as well, we need to perform a finite difference on the entire likelihood for each non-mu modeled \(\theta_{{\enclose{horizontalstrike}{\mu}}, j} \in \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\)

\[ \frac{ \partial L_i }{ \partial \theta_{{\enclose{horizontalstrike}{\mu}}, j} } \approx \frac{ L_i(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j}) - L_i(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) }{ \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j} } = \tag{31}\]

\[ \int_{-\infty}^{\infty}\left[\frac{ -\log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) + \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j} }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \approx \tag{32}\]

\[ \sum_{k=1}^{r} z_{i,k,y}\left[\frac{ -\log\bigl(p(\mathbf{y}_{i,k} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) + \log\bigl(p(\mathbf{y}_{i,k} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j} }\right] \approx \tag{33}\]

\[ E_{\boldsymbol{\phi}, i}\left(\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \theta_{{\enclose{horizontalstrike}{\mu}}, j} }\right) = \mathbf{g}_{\theta_{{\enclose{horizontalstrike}{\mu}}, j}} \tag{34}\]

\[ \begin{aligned} \frac{ \partial L }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } &= \sum_{i=1}^{m} \frac{ \partial L_i }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } \\ &= \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \\ &= \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \\ &= \mathbf{0} \end{aligned} \tag{35}\]

where \(\mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}}\) is the vector of all \(\mathbf{g}_{\theta_{{\enclose{horizontalstrike}{\mu}}, j}}\) for which \(\theta_{{\enclose{horizontalstrike}{\mu}}, j} \in \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\). A Hessian matrix may be constructed as follows:

\[ \mathbf{H}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} = \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i}^{\prime} \tag{36}\]

\[ \hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}} = \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \alpha \mathbf{H}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}}^{-1} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \tag{37}\]

We now consider the computational expense for importance sampling required to update mu modeled theta parameters versus non-mu modeled parameters. For complex PK/PD problems that use the numerical integration ($DES), the greatest computational expense is in evaluating the predictive function \(\mathbf{f}_i(\mathbf{t}, \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\). The evaluation of the individual likelihood is evaluated \(r\) times in the evaluation of the conditional means and variances, per subject per iteration, regardless of how many mu modeled parameters are to be evaluated, according to equation 18. Once the conditional means and variances are determined, the individual objective function is no longer needed to evaluate the update for these thetas, according to equations 20 and 22.

For non-MU modeled parameters, however, equations 33, 34, 35, 36, and 37 suggest that \(n_{\enclose{horizontalstrike}{\mu}} r\) individual objective function calls are required, where \(n_{\enclose{horizontalstrike}{\mu}}\) is the number of non-mu modeled parameters, one for each \(\log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid (\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j}), \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr)\) evaluation.

There is a sub-class of non-mu modeled parameters for which some computation efficiency can be made, and these are the SIGMA parameters, or Sigma-like, theta parameters. Such parameters are not used in evaluating the predictive function \(\mathbf{f}_i(\mathbf{t}, \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\), the most computationally expensive component, but only in evaluating the residual variance

\(\mathbf{V}_i(\mathbf{f}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\), so NONMEM uses \(\mathbf{f}_i(\mathbf{t}, \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\) in evaluating \(\log(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}))\) as well as \(\log(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j}))\) during the finite difference step, and will not re-evaluate \(\mathbf{f}\):

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

\[ \begin{aligned} &-2\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j})\bigr) = \\ &\quad [\mathbf{y}_i - \mathbf{f}_i(\boldsymbol{\phi})]^{\prime} \mathbf{V}_i^{-1}(\mathbf{f}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j})[\mathbf{y}_i - \mathbf{f}_i(\boldsymbol{\phi})] + \log\bigl(\det(\mathbf{V}_i(\mathbf{f}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta\theta_{{\enclose{horizontalstrike}{\mu}}, j}))\bigr) \end{aligned} \tag{39}\]

Note that for these parameters only the \(\mathbf{V}_i(\mathbf{f}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\) has to be re-evaluated (as the Y value in the NONMEM control stream file), which is usually a simple algebraic relation. SIGMA parameters are automatically recognized by NONMEM as those for which it can make this short-cut. THETA parameters that are used only in evaluating the residual variance (in the evaluation of Y in the control stream file) but not, directly or indirectly, in evaluating the predictive function (in the evaluation of F in the control stream file), may be given an S designation in the GRD setting of $EST, and only then will NONMEM utilize the short cut for evaluating its partial derivative.

Sigma parameters (but not Sigma-like THETA parameters) can be additionally updated efficiently by evaluating their partial derivative gradient contributions analytically, as given in Appendix F (available in NONMEM 7.2). However, if the user specifies that Sigma’s GRD value with an N, then their partial derivatives are evaluated numerically by finite difference method.

In general therefore, it is best to mu model THETA parameters whenever possible, to take advantage of the efficiency available for EM methods, and to specify when thetas may be considered Sigma-like, or to take advantage of modeling residual variances via SIGMA parameters, as much as possible.