112 MCMC Bayesian Analysis for Evaluating a Distribution of Population Parameters
NONMEM 7 Technical Guide
The Markov chain Monte Carlo (MCMC) Bayesian analysis can be used to obtain many thousands of population parameter and variance parameters that represent the distribution of population parameters
according to their ability to fit the data. This information is similar to what is obtained by boot strap methods, and MCMC Bayesian analysis can be used in their place. The Bayesian analysis may be performed with or without including prior information, but it is recommended that there at least be prior information for OMEGA.
There are two main types of Bayesian analysis available in NONMEM. The most efficient is the Gibbs sampling method, and is used to create samples of thetas that are linearly modeled with respect to their MU’s, and the inter-subject variances. This is performed in the manner of page 341 of (Bennett et al. 1996). Updating linearly modeled thetas (designated as \(\boldsymbol{\theta}_\mu\)) is done as follows.
Use the EM update method to obtain estimates \(\hat{\boldsymbol{\theta}}_\mu\):
\[ \mathbf{H}_{\boldsymbol{\theta}_\mu} = \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } + \boldsymbol{\Omega}_\theta^{-1} \tag{1}\]
Followed by
\[ \hat{\boldsymbol{\theta}}_\mu = \boldsymbol{\theta}_\mu + \alpha \mathbf{H}_{\boldsymbol{\theta}_\mu}^{-1} \mathbf{g}_{\boldsymbol{\theta}_\mu} \tag{2}\]
Next, sample from the following conditional density:
\[ [\boldsymbol{\theta}_\mu \mid \cdot] \sim N\bigl(\hat{\boldsymbol{\theta}}_\mu, \mathbf{H}_{\boldsymbol{\theta}_\mu}^{-1}\bigr) \tag{3}\]
For the Omegas:
\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m + \rho}\left[\sum_{i=1}^{m} (\bar{\boldsymbol{\phi}} - \boldsymbol{\mu})(\bar{\boldsymbol{\phi}} - \boldsymbol{\mu})^{\prime} + \sum_{i=1}^{m} \bar{\mathbf{B}}_i + \rho \boldsymbol{\Omega}_\Omega\right] \tag{4}\]
Followed by sampling from an inverse Wishart density:
\[ [\boldsymbol{\Omega}^{-1} \mid \cdot] \sim W\bigl(\bigl((m + \rho)\hat{\boldsymbol{\Omega}}\bigr)^{-1}, m + \rho\bigr) \tag{5}\]
A matrix with an inverse Wishart distribution of \(m + \rho\) degrees of freedom could be constructed as follows. Create \(k\) vectors of normally distributed random samples:
\[ \mathbf{x}_i \sim N(0, 1) \tag{6}\]
Then construct
\[ \mathbf{S}_{m+\rho} = \sum_{k=1}^{m+\rho} \mathbf{x}_i \mathbf{x}_i^{\prime} \tag{7}\]
\[ \boldsymbol{\Omega} = \mathbf{L}_\Omega \mathbf{S}^{-1} \mathbf{L}_\Omega^{\prime} \tag{8}\]
Where \(\mathbf{L}_\Omega\) is the cholesky of \(\hat{\boldsymbol{\Omega}}\). More efficient methods of creating an inverse Wishart matrix sample are available. Because these sample densities are also the conditional densities for the respective parameters, the samples are always accepted, and no acceptance/rejection analysis needs to be performed.
Sigma parameters (but not Sigma-like THETA parameters) that are isolated residual variance coefficients are updated as follows:
\[ \hat{\sigma}^2 = \sum_{i=1}^{m} \sum_{j=1}^{n_i} \frac{ (y_{ij} - f_j)^2 }{ (\partial y_{ij}/\partial \varepsilon)^2 } \tag{9}\]
Followed by sampling from an inverse chi-square:
\[ \sigma^{-2} = \chi^{-2}(N, \hat{\sigma}^{-2}) \]
Where \(N\) is the total number of data points involved in evaluation of that particular sigma.
Metropolis-Hastings sampling must be performed on all other types of theta parameters, as follows. For the first mode, for thetas not linearly modeled and cholesky decomposed sigma elements, designated collectively as \(\boldsymbol{\theta}\), proposed sample parameters for the \(k\)+1th iteration are created using
\[ \log(N(\boldsymbol{\theta} \mid \boldsymbol{\theta}_s, \mathbf{Z})) = -\frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}_s)^{\prime} (\mathbf{Z})^{-1} (\boldsymbol{\theta} - \boldsymbol{\theta}_s) - \frac{1}{2}\log|\mathbf{Z}| \tag{10}\]
\(\boldsymbol{\theta}_s\) and \(\mathbf{Z}\) vary according to how many samples have so far been created. During the first several hundred iterations of burn-in, \(\boldsymbol{\theta}_s\) are the initial thetas at the start of the analysis, and \(\mathbf{Z}\) is a diagonal matrix with diagonal elements that are equal to \((0.5 * \boldsymbol{\theta}_s)^2\). During the subsequent iterations of burn-in, \(\boldsymbol{\theta}_s\) and \(\mathbf{Z}\) are the sample means and variances of \(\boldsymbol{\theta}\) collected during the previous several hundred iterations. During the stationary phase, \(\boldsymbol{\theta}_s\) and \(\mathbf{Z}\) are the sample means and variances of all \(\boldsymbol{\theta}\) collected so far since the beginning of the stationary phase.
To reflect the probability of choosing these values, the following log density values are therefore calculated, based on the respective proposal densities, for mode 1
\[ \log(k_1(\boldsymbol{\theta} \mid \boldsymbol{\theta}_s)) = \log(N(\boldsymbol{\theta} \mid \boldsymbol{\theta}_s, \mathbf{Z}_s)) \tag{11}\]
The log likelihood of the \(k\)th set of population parameters with respect to the data, and with respect to positions of the \(k\)th set of individual parameters \(\boldsymbol{\phi}_{k,i}\) is evaluated also:
\[ \log(\pi(\boldsymbol{\theta}_k)) = \sum_{i=1}^{m} \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi}_{k,i} \mid \boldsymbol{\theta}_k, \boldsymbol{\mu}_{k,i}, \boldsymbol{\Omega}_k)\bigr) \tag{12}\]
The log likelihood of the proposed sample set of population parameters with respect to the data, and with respect to positions of the present \(k\)th set of individual thetas is evaluated also:
\[ \log(\pi(\boldsymbol{\theta})) = \sum_{i=1}^{m} \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi}_i \mid \boldsymbol{\theta}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) \tag{13}\]
The following test statistic is created:
\[ t_1 = -\log(k_1(\boldsymbol{\theta} \mid \boldsymbol{\theta}_s)) + \log(k_1(\boldsymbol{\theta}_s \mid \boldsymbol{\theta}_s)) - \log(\pi(\boldsymbol{\theta}_s)) + \log(\pi(\boldsymbol{\theta})) \tag{14}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_1 \tag{15}\]
then the proposed sample set of population parameters is accepted. If rejected, then the \(k\)th sample set is used as the \(k\)+1th sample set. This is done PSAMPLE_M1 (an option in NONMEM) times.
Next, during the second kernel density mode, the population parameters of the present position \(k\) may be used to create a sample for the next iteration:
\[ \log(N(\boldsymbol{\theta} \mid \boldsymbol{\theta}_k, \mathbf{Z})) = -\frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}_k)^{\prime} (w\mathbf{Z})^{-1}(\boldsymbol{\theta} - \boldsymbol{\theta}_k) - \frac{1}{2}\log|w\mathbf{Z}| \tag{16}\]
Where \(\boldsymbol{\theta}_k\) is the accepted theta of the \(k\)th iteration, \(w\) is a scaling parameter, which is adjusted throughout the analysis so that a fraction PACCEPT (option) of random sample sets are accepted. The PACCEPT (option) parameter is set by the user.
To reflect the probability of choosing these values, based on the respective proposal densities:
\[ \log(k_2(\boldsymbol{\theta} \mid \boldsymbol{\theta}_k)) = \log(N(\boldsymbol{\theta} \mid \boldsymbol{\theta}_k, w\mathbf{Z})) \tag{17}\]
as well as their backward density of mode 2:
\[ \begin{aligned} \log(k_2(\boldsymbol{\theta}_k \mid \boldsymbol{\theta})) &= \log(N(\boldsymbol{\theta}_k \mid \boldsymbol{\theta}, w\mathbf{Z})) \\ &= \log(k_2(\boldsymbol{\theta} \mid \boldsymbol{\theta}_k)) \end{aligned} \tag{18}\]
The test statistic is created:
\[ t_2 = -\log(k_2(\boldsymbol{\theta} \mid \boldsymbol{\theta}_k)) + \log(k_2(\boldsymbol{\theta}_k \mid \boldsymbol{\theta})) - \log(\pi(\boldsymbol{\theta}_k)) + \log(\pi(\boldsymbol{\theta})) \tag{19}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_2 \tag{20}\]
Then the sample is accepted. This is done PSAMPLE_M2 times.
As a third kernel sampling mode, samples on each parameter separately and sequentially may be made using the univariate distribution
\[ \log(N(\theta_l \mid \theta_{k,l}, z_{l,l}^{-1})) = -\frac{1}{2}(\theta_l - \theta_{k,l}) z_{l,l}^{-1} (\theta_l - \theta_{k,l}) + \frac{1}{2}\log\left|z_{l,l}^{-1}\right| \tag{21}\]
where \(z_{l,l}^{-1}\) is the \(l\)th diagonal element to \(\mathbf{Z}^{-1}\), for parameter \(l\). The other parameters are not moved when in this mode.
To reflect the probability of choosing these values, the following log density values are therefore calculated, based on the respective proposal densities:
\[ \log(k_3(\theta_l \mid \theta_{k,l})) = \log(N(\theta_l \mid \theta_{k,l}, w\mathbf{Z})) \tag{22}\]
and backward density of mode 3:
\[ \log(k_3(\theta_{k,l} \mid \theta_l)) = \log(N(\theta_{k,l} \mid \theta_l, w\mathbf{Z})) \tag{23}\]
The test statistic is created for each parameter \(l\):
\[ t_3 = -\log(k_3(\theta_l \mid \theta_{k,l})) + \log(k_3(\theta_{k,l} \mid \theta_l)) - \log(\pi(\boldsymbol{\theta}_k)) + \log(\pi(\theta_{k,l})) \tag{24}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_3 \tag{25}\]
then the proposed sample set of population parameters is accepted as the \(k\)+1th sample set. If rejected, then the \(k\)th sample set is used as the \(k\)+1th sample set.
The third mode is done for each parameter PSAMPLE_M3 times, for \(n\)*PSAMPLE_M3 times in a given iteration, where \(n\) is the number of population parameters in the vector \(\boldsymbol{\theta}\).
If the user has selected to perform Metropolis-Hastings samplings for Omega elements, then for each time that samples of population mean parameters and covariates are created, samples of population variances are also created using the inverse Wishart distribution.
For mode 1, using the starting position values (\(k\)=0) (OSAMPLE_M1 times):
\[ \begin{aligned} &\log\bigl(W^{-1}(\boldsymbol{\Omega}^{-1} \mid [(\rho + m)\boldsymbol{\Omega}_s]^{-1}, (\rho + m))\bigr) = \\ &\quad -\frac{1}{2}\left((\rho + m) \operatorname{tr}(\boldsymbol{\Omega}_s \boldsymbol{\Omega}^{-1}) + (\rho + m - n - 1)\ln(|\boldsymbol{\Omega}|) - (\rho + m)\bigl[\ln(|\boldsymbol{\Omega}_s|) + n\ln(\rho + m)\bigr]\right) \end{aligned} \tag{26}\]
To reflect the probability of choosing these values, the following log density values are therefore calculated, based on the respective proposal densities:
\[ \log(k_1(\boldsymbol{\Omega} \mid \boldsymbol{\Omega}_s)) = \log\bigl(W^{-1}(\boldsymbol{\Omega} \mid m \boldsymbol{\Omega}_s, \rho + m))\bigr) \tag{27}\]
The log likelihood of the \(k\) set of population parameters with respect to the data, and with respect to positions of the \(k\) set of individual parameters \(\boldsymbol{\phi}_i\) is evaluated also:
\[ \log(\pi(\boldsymbol{\Omega}_k)) = \sum_{i=1}^{m} \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi}_i \mid \boldsymbol{\theta}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}_k)\bigr) \tag{28}\]
The log likelihood of the proposed sample set of population variances with respect to the data, and with respect to positions of the present \(k\)th set of individual thetas is evaluated also:
\[ \log(\pi(\boldsymbol{\Omega})) = \sum_{i=1}^{m} \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi}_i \mid \boldsymbol{\theta}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) \tag{29}\]
During mode 1, the following test statistic is created:
\[ t_1 = -\log(k_1(\boldsymbol{\Omega} \mid \boldsymbol{\Omega}_s)) + \log(k_1(\boldsymbol{\Omega}_s \mid \boldsymbol{\Omega}_s)) - \log(\pi(\boldsymbol{\Omega}_s)) + \log(\pi(\boldsymbol{\Omega})) \tag{30}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_1 \tag{31}\]
then the proposed sample set of variances is accepted as the \(k\)+1th sample set. If rejected, then the \(k\)th sample set is used as the \(k\)+1th sample set.
For mode 2, the present position \(k\) is used (OSAMPLE_M2 times):
\[ \begin{aligned} &\log\bigl(W^{-1}(\boldsymbol{\Omega}^{-1} \mid [w(\rho + m)\boldsymbol{\Omega}_k]^{-1}, w(\rho + m))\bigr) = \\ &\quad -\frac{1}{2}\left(w(\rho + m) \operatorname{tr}(\boldsymbol{\Omega}_k \boldsymbol{\Omega}^{-1}) + (w(\rho + m) - n - 1)\ln(|\boldsymbol{\Omega}|) - (w(\rho + m))\bigl[\ln(|\boldsymbol{\Omega}_k|) + n\ln(w(\rho + m))\bigr]\right) \end{aligned} \tag{32}\]
where \(w\) is the scaling parameter (separate from that used for the normal distribution proposal density for the theta parameters) to allow OACCEPT acceptance rate.
To reflect the probability of choosing these values, the following log density values are therefore calculated, based on the respective proposal densities:
\[ \log(k_2(\boldsymbol{\Omega} \mid \boldsymbol{\Omega}_k)) = \log\bigl(W^{-1}(\boldsymbol{\Omega} \mid w(\rho + m)\boldsymbol{\Omega}_k, w(\rho + m))\bigr) \tag{33}\]
as well as their backward density of mode 2:
\[ \log(k_2(\boldsymbol{\Omega}_k \mid \boldsymbol{\Omega})) = \log\bigl(W^{-1}(\boldsymbol{\Omega}_k \mid w(\rho + m)\boldsymbol{\Omega}, w(\rho + m))\bigr) \tag{34}\]
The test statistic is created:
\[ t_2 = -\log(k_2(\boldsymbol{\Omega} \mid \boldsymbol{\Omega}_k)) + \log(k_2(\boldsymbol{\Omega}_k \mid \boldsymbol{\Omega})) - \log(\pi(\boldsymbol{\Omega}_k)) + \log(\pi(\boldsymbol{\Omega})) \tag{35}\]
A uniform random deviate \(u\) is then generated, log transformed, and if
\[ \log(u) < t_2 \tag{36}\]
then the proposed sample set of variances is accepted, and serves as the \(k\)+1th sample set. If rejected, then the \(k\)th sample set is used as the \(k\)+1th sample set. This is done OSAMPLE_M2 times.
A single iteration consists of: Gibbs sampling of THETAS, SIGMAS and OMEGAS, followed by Metropolis-Hastings sampling of other THETAS, PSAMPLE_M1 times for mode 1, then PSAMPLE_M2 times for mode 2, followed by OMEGAS sampled OSAMPLE_M1 times for mode 1, then OSAMPLE_M2 times for mode 2. The final sample set of THETAS, OMEGAS and SIGMAS after going through this process is then stored in the raw output file as the results for that particular iteration.