108  Iterative Two Stage (ITS)

NONMEM 7 Technical Guide

Iterative two stage (ITS) approximates the expectation step by using the conditional modes and approximate conditional variances that are evaluated during the MAP estimation method that is also used in FOCE or LAPLACE methods, as described earlier. We can consider an approximate update for MU-modeled thetas that is applied in ITS, by evaluating:

\[ \begin{aligned} \frac{\partial L_i}{\partial \boldsymbol{\mu}_i} &= -\boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \\ &\approx -\boldsymbol{\Omega}^{-1}(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \\ &= -\boldsymbol{\Omega}^{-1} \hat{\boldsymbol{\eta}}_i \\ &= \mathbf{g}_{\boldsymbol{\mu}_i} \\ &= \mathbf{0} \\ &= \frac{\partial L_{Ai}}{\partial \boldsymbol{\mu}_i} \end{aligned} \tag{1}\]

Where subscript A refers to “approximate”. This is an approximation to that extent that the mode \(\hat{\boldsymbol{\phi}}_i = \boldsymbol{\mu}_i(\boldsymbol{\theta}_\mu) + \hat{\boldsymbol{\eta}}_i\) approximates the mean \(\bar{\boldsymbol{\phi}}_i\). Then as before,

\[ \begin{aligned} \frac{\partial L_A}{\partial \boldsymbol{\theta}_\mu} &= \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \frac{\partial L_{Ai}}{\partial \boldsymbol{\mu}_i} \\ &= -\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1} \hat{\boldsymbol{\eta}}_i \\ &= \mathbf{g}_{\boldsymbol{\theta}_\mu} \\ &= \mathbf{0} \end{aligned} \tag{2}\]

We then perform a Gauss-Newton update:

\[ \mathbf{H}_{\boldsymbol{\theta}_\mu} = \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\mu}_i} \mathbf{g}_{\boldsymbol{\mu}_i}^{\prime} \tag{3}\]

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

Then, updating the MUs: \(\hat{\boldsymbol{\mu}}_i = \boldsymbol{\mu}_i(\hat{\boldsymbol{\theta}})\)

Similarly, to update Omega, iterative two stage approximates update 41:

\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m} \sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \bar{\mathbf{B}}_i \tag{5}\]

\[ \approx \frac{1}{m} \sum_{i=1}^{m} (\hat{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}}_i)(\hat{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}}_i)^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \hat{\mathbf{B}}_i \tag{6}\]

where

\[ \hat{\mathbf{B}}_i = (\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1} \tag{7}\]

is the approximate conditional variance evaluated during the FOCE or LAPLACE integration step.

The approximate optimization of the ITS method is related to an approximate optimization of FOCE’s or Laplace’s \(L_N\). To consider optimizing \(L_N\) for MU-modeled thetas, we can conveniently rephrase equation 5 as

\[ \begin{aligned} L_{Ni} &= -\log\bigl(l(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) \\ &\quad + \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) \\ &\quad + \frac{1}{2}\log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i))\bigr) \end{aligned} \tag{8}\]

Since \(\hat{\boldsymbol{\phi}}_i\) is at the mode of the posterior density, then

\[ \left(\frac{ \partial -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) }{ \partial \boldsymbol{\phi} }\right)_{\hat{\boldsymbol{\phi}}_i} + \boldsymbol{\Omega}^{-1}(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) = \mathbf{0} \tag{9}\]

It follows that:

\[ \begin{aligned} \frac{\partial L_N}{\partial \boldsymbol{\theta}_\mu} &= \sum_{i=1}^{m} \frac{\partial L_{Ni}}{\partial \boldsymbol{\theta}_\mu} \\ &= \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\phi}}_i }{ \partial \boldsymbol{\theta}_\mu } \left[ \left(\frac{ \partial -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) }{ \partial \boldsymbol{\phi} }\right)_{\hat{\boldsymbol{\phi}}_i} + \frac{1}{2} \frac{ \partial (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) }{ \partial \hat{\boldsymbol{\phi}}_i } \right] \\ &\quad - \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i^{\prime} }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \\ &\quad + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\phi}}_i }{ \partial \boldsymbol{\theta}_\mu } \frac{ \partial \log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i))\bigr) }{ \partial \hat{\boldsymbol{\phi}}_i } \end{aligned} \tag{10}\]

\[ \begin{aligned} \frac{\partial L_N}{\partial \boldsymbol{\theta}_\mu} &= -\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i^{\prime} }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \\ &\quad + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\phi}}_i }{ \partial \boldsymbol{\theta}_\mu } \frac{ \partial \log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\phi}}_i))\bigr) }{ \partial \hat{\boldsymbol{\phi}}_i } \\ &= \mathbf{0} \end{aligned} \tag{11}\]

where the term in parentheses cancels because of equation 9. Comparing equation 11 with 2 shows that ITS only approximates the optimization of FOCE’s \(L_N\) because it does not include the contribution of change in the information matrix of the joint density with respect to theta (the \(\log(\det)\) term in equation 11).

Similarly, we consider differentiating FOCE’s objective function with respect to OMEGA:

\[ \begin{aligned} \frac{\partial L_N}{\partial \boldsymbol{\Omega}^{-1}} &= \sum_{i=1}^{m} \frac{\partial L_{Ni}}{\partial \boldsymbol{\Omega}^{-1}} \\ &= \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\phi}}_i }{ \partial \boldsymbol{\Omega}^{-1} } \left[ \left(\frac{ \partial -\log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})\bigr) }{ \partial \boldsymbol{\phi} }\right)_{\hat{\boldsymbol{\phi}}_i} + \frac{1}{2} \frac{ \partial (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) }{ \partial \hat{\boldsymbol{\phi}}_i } \right] \\ &\quad - \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \log(\det(\boldsymbol{\Omega}^{-1})) }{ \partial \boldsymbol{\Omega}^{-1} } \\ &\quad + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) }{ \partial \boldsymbol{\Omega}^{-1} } \\ &\quad + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta}))\bigr) }{ \partial \boldsymbol{\Omega}^{-1} } \\ &\quad + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\eta}}_i }{ \partial \boldsymbol{\Omega}^{-1} } \frac{ \partial \log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta}))\bigr) }{ \partial \hat{\boldsymbol{\eta}}_i } \end{aligned} \tag{12}\]

\[ \begin{aligned} \frac{\partial L_N}{\partial \boldsymbol{\Omega}^{-1}} &= \sum_{i=1}^{m} \left[ -\frac{1}{2} \boldsymbol{\Omega} + \frac{1}{2} \hat{\boldsymbol{\eta}}_i \hat{\boldsymbol{\eta}}_i^{\prime} + \frac{1}{2}\bigl(\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}\bigr)^{-1} \right] \\ &\quad + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\eta}}_i }{ \partial \boldsymbol{\Omega}^{-1} } \frac{ \partial \log(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1})) }{ \partial \hat{\boldsymbol{\eta}}_i } \\ &= \mathbf{0} \end{aligned} \tag{13}\]

Here we differentiate the objective function with respect to \(\boldsymbol{\Omega}^{-1}\) for convenience. When the gradient with respect to \(\boldsymbol{\Omega}^{-1}\) equals \(\mathbf{0}\), then the gradient with respect to \(\boldsymbol{\Omega}\) also equals \(\mathbf{0}\). The details of the linear algebra manipulations leading to the last part of 13 are given in appendix A. Then we can express the exact minimization of \(L_N\) with respect to Omega as:

\[ \begin{aligned} \frac{\partial L_N}{\partial \boldsymbol{\Omega}^{-1}} &= \sum_{i=1}^{m} \left[ -\frac{1}{2} \boldsymbol{\Omega} + \frac{1}{2} \hat{\boldsymbol{\eta}}_i \hat{\boldsymbol{\eta}}_i^{\prime} + \frac{1}{2} \hat{\mathbf{B}}_i \right] + \frac{1}{2} \sum_{i=1}^{m} \frac{ \partial \hat{\boldsymbol{\eta}}_i }{ \partial \boldsymbol{\Omega}^{-1} } \frac{ \partial \log(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1})) }{ \partial \hat{\boldsymbol{\eta}}_i } \\ &= \mathbf{0} \end{aligned} \tag{14}\]

Note that \(\hat{\mathbf{B}}_i\) represents a linearized approximation to the true conditional variance \(\bar{\mathbf{B}}_i\). We may consider an approximate partial gradient to \(L_N\) with respect to Omega as:

\[ \begin{aligned} \frac{\partial L_N}{\partial \boldsymbol{\Omega}^{-1}} &= \sum_{i=1}^{m} \frac{\partial L_{Ni}}{\partial \boldsymbol{\Omega}^{-1}} \\ &\approx \sum_{i=1}^{m}\left[-\frac{1}{2}\boldsymbol{\Omega} + \frac{1}{2}\hat{\boldsymbol{\eta}}_i \hat{\boldsymbol{\eta}}_i^{\prime} + \frac{1}{2}\hat{\mathbf{B}}_i\right] \\ &= \mathbf{0} \end{aligned} \tag{15}\]

Solving for the next estimate of Omega from equation 15:

\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m} \sum_{i=1}^{m} \hat{\boldsymbol{\eta}}_i \hat{\boldsymbol{\eta}}_i^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \hat{\mathbf{B}}_i \tag{16}\]

which is similar to the update of Omega in the ITS algorithm 6. Again, the ITS method updates \(L_N\) only approximately, as it does not take into account the \(\log(\det)\) term in equation 14.

To summarize, in ITS, \(\boldsymbol{\theta}\) is updated here using the average of the modes of the individual joint densities, which serves only as an approximation to the more precise update of the average of the means of the \(\boldsymbol{\phi}\) under the joint density, as dictated by the exact equation 14. If there is a skewness to each individual’s joint density, such that the means tend to differ systematically from the modes, then the iterative two stage update may yield biased results.

For non-MU-modeled theta, ITS in NONMEM falls back on a forward difference evaluation of the full likelihood:

\[ \frac{ \partial L_N }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } \approx \frac{ L_N(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \Delta \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) - L_N(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) }{ \Delta \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } = \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \tag{17}\]

Followed by a single step Gauss-Newton update:

\[ \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{18}\]

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

To summarize, the NONMEM FOCE method optimizes \(L_N\), which is an approximation to the true objective function \(L\), and ITS further approximates the optimization of \(L_N\).

All iterative update methods that rely on updating the population parameters \(\boldsymbol{\theta}\) using the average to the individual estimates guarantee “centeredness” of the population parameters about the individual parameters by definition. However, because the FOCE NONMEM

method uses a search algorithm on an approximate objective function, it does not guarantee centeredness. One can impose the CENTERING option to the estimation process in NONMEM, which then optimizes a modified objective function of equation 5:

\[ \begin{aligned} L'_{Ni} &\approx - \log\bigl(l(\mathbf{y}_i \mid (\hat{\boldsymbol{\eta}}_i - \bar{\boldsymbol{\eta}}), \boldsymbol{\theta})\bigr) \\ &\quad + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) \\ &\quad + \frac{1}{2}(\hat{\boldsymbol{\eta}}_i - \bar{\boldsymbol{\eta}})^{\prime} \boldsymbol{\Omega}^{-1} (\hat{\boldsymbol{\eta}}_i - \bar{\boldsymbol{\eta}}) \\ &\quad + \frac{1}{2}\log\bigl(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^{-1}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_i, \boldsymbol{\theta}))\bigr) \end{aligned} \tag{20}\]

where

\[ \bar{\boldsymbol{\eta}} = \frac{1}{m} \sum_{i=1}^{m} \hat{\boldsymbol{\eta}}_i \tag{21}\]

to ensure statistical centering, although not exact centering.