111  Population Mixture Modeling

NONMEM 7 Technical Guide

Sometimes the data may be derived from two or more sub-populations, as evidenced by a distribution of a parameter among the subjects that appears to be bi-modal, or skewed. For example, suppose the data is first fit with a simple one compartment model, with volume \(V_c\) and rate constant of elimination \(k10\). A histogram analysis of the individual \(k10\)’s suggests a bimodal or skewed distribution. However, none of the known binary covariates (gender, for example) explains this bimodality. Under these circumstances, one can specify the probability of an individual belonging to a sub-group, without insisting on the certainty of belonging to that particular sub-group.

Consider that we have \(N_m\) sub-populations. Then for each subject \(i\), and for each sub-population \(j\) we have the probability

\[ p(\mathbf{y}_i \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) = \int_{-\infty}^{\infty} \sum_{j=1}^{N_m} a_j p_j(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{1}\]

where

\[ p_j(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) \tag{2}\]

is the density for sub-population model \(j\), for subject \(i\), and \(a_j\) is the probability of belonging to sub-population \(j\). Then define

\[ L_{ji} = -\log\left(\int_{-\infty}^{\infty} p_j(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) \tag{3}\]

so the negative log-likelihood of an individual is:

\[ L_i = -\log\left(\int_{-\infty}^{\infty} \sum_{j=1}^{N_m} a_j p_j(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) = \tag{4}\]

\[ -\log\left(\sum_{j=1}^{N_m} a_j \int_{-\infty}^{\infty} p_j(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) = \tag{5}\]

\[ -\log\left(\sum_{j=1}^{N_m} a_j \exp(-L_{ji})\right) \tag{6}\]

Consider that equations for updating the non-proportion (that is, non-\(a\)) population parameters \(\mathbf{q} = \{\boldsymbol{\theta}, \boldsymbol{\Omega}\}\) are derived from obtaining the partial derivatives of the objective function \(L\):

\[ \frac{\partial L_i}{\partial \mathbf{q}} = \sum_{j=1}^{N_m} \frac{ a_j \exp(-L_{ji}) \dfrac{\partial L_{ji}}{\partial \mathbf{q}} }{ \sum_{k=1}^{N_m} a_k \exp(-L_{ki}) } \tag{7}\]

or

\[ \begin{aligned} \frac{\partial L_i}{\partial \mathbf{q}} &= \sum_{j=1}^{N_m} \frac{ a_j r_{ji} \dfrac{\partial L_{ji}}{\partial \mathbf{q}} }{ \sum_{k=1}^{N_m} a_k r_{ki} } \\ &= \sum_{j=1}^{N_m} a_{ji} \frac{\partial L_{ji}}{\partial \mathbf{q}} \end{aligned} \tag{8}\]

where

\[ r_{ji} = \exp(-L_{ji}) \tag{9}\]

and

\[ a_{ji} = \frac{a_j r_{ji}}{\sum_{k=1}^{N_m} a_k r_{ki}} \tag{10}\]

is the probability or weight for individual \(i\), sub-population model \(j\). As an example,

\[ \begin{aligned} \frac{\partial L_{ji}}{\partial \boldsymbol{\mu}_i} &= -\boldsymbol{\Omega}^{-1} \sum_{i=1}^{m} \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \boldsymbol{\mu}_i) z_j(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \\ &= \mathbf{g}_{\boldsymbol{\mu}_i, j} \end{aligned} \tag{11}\]

where

\[ z_j(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) \tag{12}\]

is the conditional density for subject \(i\) modeled under sub-model \(j\), then the appropriate conditional mean for subject \(i\) would be

\[ \bar{\bar{\boldsymbol{\phi}}}_i = \sum_{j=1}^{N_m} a_{ji} \bar{\boldsymbol{\phi}}_{ji} \tag{13}\]

where

\[ \bar{\boldsymbol{\phi}}_{ji} = \int_{-\infty}^{\infty} \boldsymbol{\phi} z_j(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{14}\]

which are then used in the usual way to update the thetas. Similarly:

\[ \bar{\bar{\mathbf{B}}}_i = \sum_{j=1}^{N_m} a_{ji} (\bar{\mathbf{B}}_{ji} + \bar{\boldsymbol{\phi}}_{ji} \bar{\boldsymbol{\phi}}_{ji}^{\prime}) - \bar{\bar{\boldsymbol{\phi}}}_i \bar{\bar{\boldsymbol{\phi}}}_i^{\prime} \tag{15}\]

where

\[ \bar{\mathbf{B}}_{ji} = \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \bar{\boldsymbol{\phi}}_{ji})(\boldsymbol{\phi} - \bar{\boldsymbol{\phi}}_{ji})^{\prime} z_j(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{16}\]

is the conditional variance for individual \(i\), sub-model \(j\), whereupon the update is the usual:

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

The weighted average of the other expectation results are also performed, using the same weightings. The \(L_{ji}\), and therefore \(a_{ji}\) is readily obtained during the expectation step as the objective function to subject \(i\), under sub-population model \(j\). In practice therefore, the expectation step is done \(N_m\) times for each subject, collecting the resulting conditional means, variances, and objective function values to each sub-model, and then performing the weighted average, as shown above.

A method in keeping with minimizing the total objective function would be to construct partial derivatives and partial second derivatives, where for each subject \(i\):

\[ \begin{aligned} \frac{\partial L_i}{\partial a_j} &= -\frac{ \exp(-L_{ji}) - \exp(-L_{N_m, i}) }{ \sum_{k=1}^{N_m} a_k \exp(-L_{ki}) } \\ &= -\frac{r_{ji} - r_{N_m, i}}{\sum_{k=1}^{N_m} a_k r_{ki}} \end{aligned} \tag{18}\]

\[ \begin{aligned} -\frac{\partial^2 L_i}{\partial a_{j_1} \partial a_{j_2}} &= \frac{ (\exp(-L_{j_1, i}) - \exp(-L_{N_m, i}))(\exp(-L_{j_2, i}) - \exp(-L_{N_m, i})) }{ \left(\sum_{k=1}^{N_m} a_k \exp(-L_{ki})\right)^2 } \\ &= \frac{\partial L_i}{\partial a_{j_1}} \frac{\partial L_i}{\partial a_{j_2}} \end{aligned} \tag{19}\]

since

\[ a_{N_m} = 1 - \sum_{j=1}^{N_m - 1} a_j \tag{20}\]

Then, perform the usual Gauss-Newton update, where \(\boldsymbol{\theta}_a\) are all thetas that model the sub-population proportions in the $MIX module:

\[ \mathbf{g}_{\boldsymbol{\theta}_a} = \sum_{i=1}^{m} \frac{\partial L_i}{\partial \mathbf{a}} \frac{\partial \mathbf{a}}{\partial \boldsymbol{\theta}_a} \tag{21}\]

\[ \begin{aligned} \mathbf{H}_{\boldsymbol{\theta}_a} &= \sum_{i=1}^{m} \frac{\partial L_i}{\partial \boldsymbol{\theta}_a} \frac{\partial L_i}{\partial \boldsymbol{\theta}_a}^{\prime} \\ &= \sum_{i=1}^{m} \mathbf{g}_{\mathbf{a}, i} \mathbf{g}_{\mathbf{a}, i}^{\prime} \end{aligned} \tag{22}\]

\[ \boldsymbol{\theta}_{\text{new}} = \boldsymbol{\theta}_{\text{old}} - \bigl(\mathbf{H}_{\boldsymbol{\theta}_a}\bigr)^{-1} \mathbf{g}_{\boldsymbol{\theta}_a} \tag{23}\]