109 The MCMC Method of Expectation in SAEM
NONMEM 7 Technical Guide
In Markov Monte Carlo sampling, used in the SAEM and BAYES methods, samples are generated from a larger number of proposal densities than in importance sampling. As implemented in NONMEM, for a given set of population parameters mu and omega, proposed parameters phi for each individual are generated by a three mode process. The following is based on references (Lavielle 2007; Gilks et al. 1996).
During mode 1, a vector of model parameters is generated from the following proposal density or kernel:
\[ \begin{aligned} \log(k_1(\boldsymbol{\phi})) &= \log(N(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \\ &= \log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \\ &= -\frac{1}{2}(\boldsymbol{\phi} - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi} - \boldsymbol{\mu}_i) - \frac{1}{2}\log|\boldsymbol{\Omega}| \end{aligned} \tag{1}\]
For the acceptance test, we need to evaluate the above density along with the following backward density, at the present \(\boldsymbol{\phi}_i\) for subject \(i\):
\[ \begin{aligned} \log(k_1(\boldsymbol{\phi}_i)) &= \log(N(\boldsymbol{\phi}_i \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \\ &= \log(h(\boldsymbol{\phi}_i \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \\ &= -\frac{1}{2}(\boldsymbol{\phi}_i - \boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi}_i - \boldsymbol{\mu}_i) - \frac{1}{2}\log|\boldsymbol{\Omega}| \end{aligned} \tag{2}\]
Also, the joint density is evaluated at the present \(\boldsymbol{\phi}_i\):
\[ \begin{aligned} \log(\pi(\boldsymbol{\phi}_i)) &= \log(p(\mathbf{y}_i, \boldsymbol{\phi}_i \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \\ &= \log(l(\mathbf{y}_i \mid \boldsymbol{\phi}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})) + \log(h(\boldsymbol{\phi}_i \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \end{aligned} \tag{3}\]
And at the proposed \(\boldsymbol{\phi}\)
\[ \begin{aligned} \log(\pi(\boldsymbol{\phi})) &= \log(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \\ &= \log(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})) + \log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})) \end{aligned} \tag{4}\]
Then the test statistic is created:
\[ \begin{aligned} t_1 &= \log(\pi(\boldsymbol{\phi})) - \log(\pi(\boldsymbol{\phi}_i)) + \log(k_1(\boldsymbol{\phi}_i)) - \log(k_1(\boldsymbol{\phi})) = \\ &\quad \log(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})) - \log(l(\mathbf{y}_i \mid \boldsymbol{\phi}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}})) \end{aligned} \tag{5}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_1 \tag{6}\]
then the proposed sample set \(\boldsymbol{\phi}\) of parameters is accepted and becomes the new \(\boldsymbol{\phi}_i\) for subject \(i\).
Following mode 1 sampling, proposal kernel mode 1B is performed, in which the individual conditional mean and individual conditional variance collected from previous iterations as parameters to a gaussian proposal density, is used, starting after the 10th iteration. So the proposal kernel is a type of importance sampling kernel for SAEM.
Following mode 1B sampling, proposal kernel mode 1A sampling and testing is performed, in which a sample from one of the other subjects is randomly selected. It is assumed that the set of parameters among subjects is non-parametric, with a proposal density that is proportional to the width of the nearest left and right neighbor values of the chosen parameter. This method has limited use to assist certain subjects to find good parameter values. This mode should generally not be used. Alternatively, use mode 1A sampling at the beginning of an SAEM analysis for a few burn in iterations, then continue with a complete SAEM analysis with mode 1A sampling turned off, with more burn in and accumulated sampling iterations.
Following mode 1A sampling, proposal kernel mode 2 sampling and testing is performed, using the proposal density:
\[ \begin{aligned} \log(k_2(\boldsymbol{\phi} \mid \boldsymbol{\phi}_i)) &= \log(N(\boldsymbol{\phi} \mid \boldsymbol{\phi}_i, \mathbf{Z})) \\ &= -\frac{1}{2}(\boldsymbol{\phi} - \boldsymbol{\phi}_i)^{\prime} \mathbf{Z}^{-1}(\boldsymbol{\phi} - \boldsymbol{\phi}_i)' - \frac{1}{2}\log|\mathbf{Z}| \end{aligned} \tag{7}\]
where \(\boldsymbol{\phi}_i\) is the present set of parameters for individual \(i\) (it could have been the one accepted from the just completed mode 1 sampling), and where
\[ \mathbf{Z} = \kappa \boldsymbol{\Omega} \]
which includes a scaling factor \(\kappa\). This scaling factor is adjusted for each subject such that samples are accepted at a fractional value \(\rho_M\) = IACCEPT. This scaling factor \(\kappa\) is similar to the scaling factor \(\gamma\) in importance sampling, and is also subject to the boundary values of ISCALE_MIN and ISCALE_MAX (available in NONMEM 7.2). The backward density is
\[ \begin{aligned} \log(k_2(\boldsymbol{\phi}_i \mid \boldsymbol{\phi})) &= \log(N(\boldsymbol{\phi}_i \mid \boldsymbol{\phi}, \mathbf{Z})) \\ &= -\frac{1}{2}(\boldsymbol{\phi}_i - \boldsymbol{\phi})^{\prime} \mathbf{Z}^{-1}(\boldsymbol{\phi}_i - \boldsymbol{\phi})' - \frac{1}{2}\log|\mathbf{Z}| \\ &= k_2(\boldsymbol{\phi} \mid \boldsymbol{\phi}_i) \end{aligned} \tag{8}\]
so the test statistic is calculated as:
\[ \begin{aligned} t_2 &= \log(\pi(\boldsymbol{\phi})) - \log(\pi(\boldsymbol{\phi}_i)) + \log(k_2(\boldsymbol{\phi}_i \mid \boldsymbol{\phi})) - \log(k_2(\boldsymbol{\phi} \mid \boldsymbol{\phi}_i)) \\ &= \log(\pi(\boldsymbol{\phi})) - \log(\pi(\boldsymbol{\phi}_i)) \end{aligned} \tag{9}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_2 \tag{10}\]
then the proposed sample set \(\boldsymbol{\phi}\) of parameters is accepted and becomes the new \(\boldsymbol{\phi}_i\) for subject \(i\).
For proposal kernel mode 3, each parameter of the vector \(\boldsymbol{\phi}\) is sequentially sampled using the univariate density:
\[ \begin{aligned} k_3(\boldsymbol{\phi}_l \mid \boldsymbol{\phi}_i) &= \log(N(\boldsymbol{\phi}_l \mid \boldsymbol{\phi}_{i,l}, z_{l,l}^{-1})) \\ &= -\frac{1}{2}(\boldsymbol{\phi}_l - \boldsymbol{\phi}_{i,l})(\mathbf{Z}_{l,l}^{-1})^{-1}(\boldsymbol{\phi}_l - \boldsymbol{\phi}_{i,l})' + \frac{1}{2}\log\left|z_{l,l}^{-1}\right| \end{aligned} \tag{11}\]
where subscript \(l\) refers to the \(l\)th parameter, and \(z_{l,l}^{-1}\) is the \(l\)th diagonal element of \(\mathbf{Z}^{-1}\). The backward density is
\[ k_3(\boldsymbol{\phi}_i \mid \boldsymbol{\phi}_l) = k_3(\boldsymbol{\phi}_l \mid \boldsymbol{\phi}_i) \tag{12}\]
so the test statistic is:
\[ t_1 = \log(\pi(\bar{\boldsymbol{\phi}})) - \log(\pi(\boldsymbol{\phi}_i)) \tag{13}\]
Where \(\bar{\boldsymbol{\phi}}\) equals \(\boldsymbol{\phi}_i\) but with the \(l\)th element replaced with \(\boldsymbol{\phi}_l\):
\(\bar{\boldsymbol{\phi}} = \{\boldsymbol{\phi}_{i,\sim l}, \boldsymbol{\phi}_l\}\)
Once a parameter is tested, the result contributes to the new \(\boldsymbol{\phi}_i\) for the next parameter in the vector to be sampled.
The mode 1B kernel obtains samples using the individual conditional mean and individual conditional variance collected from previous iterations as proposal density (a type of importance sampling kernel for SAEM).
During the MCMC sampling process, the IACCEPT sets \(\rho_M\), ISAMPLE_M1 determines the number of mode 1 samplings, followed by ISAMPLE_M1B samplings, followed by ISAMPLE_M1A samplings, followed by ISAMPLE_M2 mode 2 samplings, followed by ISAMPLE_M3 mode 3 samplings. The final parameter set \(\boldsymbol{\phi}_i\) after the cycle of ISAMPLE_M1 + ISAMPLE_M1B + ISAMPLE_M2 + \(p\) ISAMPLE_M3 samplings (where \(p\) = “number of elements in vector \(\boldsymbol{\phi}_i\)”) serves as the results of one chain for subject \(i\). Then, as described with importance sampling, \(r\) = ISAMPLE separate chains of vectors \(\boldsymbol{\phi}\) may be collected. Then, as described with importance sampling, the following may be calculated:
\[ \bar{\boldsymbol{\phi}}_i = \sum_{k=1}^{r} \boldsymbol{\phi}_{k,i} \tag{14}\]
\[ \bar{\mathbf{B}}_i = \sum_{k=1}^{r} (\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_i)(\boldsymbol{\phi}_{k,i} - \bar{\boldsymbol{\phi}}_i)^{\prime} \tag{15}\]
Note that the acceptance/rejection process assured that the collection of \(\boldsymbol{\phi}_{k,i}\) reflect the distribution of the desired conditional density, and weights \(z\) are not needed.
During the stochastic mode, the updates to the population parameters (both MU and non-MU-modeled) are then performed as described in importance sampling. During the accumulation mode, update results from previous \(k\)-1 iterations are averaged into the updates of the present \(k\)th iteration.
For MU-modeled theta, and Omegas, the conditional means and variances are accumulatively updated and saved as follows:
\[ \bar{\boldsymbol{\phi}}_{A,i} = \frac{k-1}{k} \bar{\boldsymbol{\phi}}_{A,i} + \frac{1}{k} \bar{\boldsymbol{\phi}}_i \]
\[ \bar{\mathbf{S}}_{A,i} = \frac{k-1}{k}\left(\bar{\mathbf{B}}_{A,i} + \bar{\boldsymbol{\phi}}_{A,i} \bar{\boldsymbol{\phi}}_{A,i}^{\prime}\right) + \frac{1}{k}\left(\bar{\mathbf{B}}_i + \bar{\boldsymbol{\phi}}_i \bar{\boldsymbol{\phi}}_i^{\prime}\right) \]
\[ \bar{\mathbf{B}}_{A,i} = \bar{\mathbf{S}}_{A,i} - \bar{\boldsymbol{\phi}}_{A,i} \bar{\boldsymbol{\phi}}_{A,i}^{\prime} \]
followed by update of the main population parameters in the usual manner:
\[ \begin{aligned} \frac{\partial L}{\partial \boldsymbol{\theta}_\mu} & = -\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_{A,i} - \boldsymbol{\mu}_i) \\ &= \mathbf{g}_{\boldsymbol{\theta}_\mu} \\ &= \mathbf{0} \end{aligned} \tag{16}\]
\[ \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{17}\]
\[ \hat{\boldsymbol{\theta}}_\mu = \boldsymbol{\theta}_\mu + \mathbf{H}_{\boldsymbol{\theta}_\mu}^{-1} \mathbf{g}_{\boldsymbol{\theta}_\mu} \tag{18}\]
\(\hat{\boldsymbol{\mu}}_i = \boldsymbol{\mu}_i(\hat{\boldsymbol{\theta}}_\mu)\)
\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m} \sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_{A,i} - \hat{\boldsymbol{\mu}}_i)(\bar{\boldsymbol{\phi}}_{A,i} - \hat{\boldsymbol{\mu}}_i)^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \bar{\mathbf{B}}_{A,i} \tag{19}\]
For non-mu modeled theta, the thetas \(\hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}}\) of the present \(k\)th iteration are updated using equations 22–27 using the present iteration’s sampling process, followed by:
\[ \hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}, k} = \bigl((k -1) \hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}, k-1} + \hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}, k}\bigr)/k \tag{20}\]
First derivative gradients of non-MU-modeled thetas are also accumulated (for use in first order approximations of standard errors, see Appendix C):
\[ \hat{\mathbf{g}}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, k} = \frac{k-1}{k} \hat{\mathbf{g}}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, k-1} + \frac{1}{k} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \tag{21}\]
In general, the order of accuracy for the various methods is:
Monte Carlo EM (IMP, SAEM, DIRECT) > Laplace > FOCEI > ITS.