113  Conditional Weighted Residuals (CWRES)

NONMEM 7 Technical Guide

We consider the following linear-epsilon residual error model:

\[ y_i = f_i(\boldsymbol{\eta}) + \sum_{m=1}^{M} q'_{im}(\boldsymbol{\eta}) \varepsilon_m \tag{1}\]

for data points \(i = 1\) to \(N\) of a particular subject, which takes into account intra-subject error components within a subject (such as homoscedastic EPS(1) mixed with heteroscedastic EPS(2)) plus possible intra-individual error interaction with other data points, as well as the possibility that \(q_{im}\) depends on \(\boldsymbol{\eta}\). We define the \(N \times 1\) normally distributed random vector:

\[ \boldsymbol{\varepsilon} = \{\varepsilon_m,\ m = 1\ \text{to}\ N\} \tag{2}\]

with the properties:

\[ \begin{aligned} E(\boldsymbol{\varepsilon}) &= \mathbf{0} \\ E(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}') &= \mathbf{I} \end{aligned} \tag{3}\]

\[ \begin{aligned} E(\varepsilon_i \eta_k) &= E(\varepsilon_i) E(\eta_k) \\ &= 0 \end{aligned} \tag{4}\]

and the \(N \times N\) matrix

\[ \mathbf{Q}' = \{q'_{im},\ i = 1\ \text{to}\ N,\ m = 1\ \text{to}\ N\} \tag{5}\]

with \(1 \times N\) row vectors

\[ \mathbf{q}'_i = \{q'_{im},\ m = 1\ \text{to}\ N\} \tag{6}\]

Then, for a given \(\boldsymbol{\eta}\), the expected value over all epsilon is

\[ \begin{aligned} E_\varepsilon(y_i) &= E(f_i(\boldsymbol{\eta})) + E\left(\sum_{m=1}^{N} q'_{im}(\boldsymbol{\eta})\varepsilon_m\right) \\ &= f_i(\boldsymbol{\eta}) + \sum_{m=1}^{N} q'_{im}(\boldsymbol{\eta}) E(\varepsilon_m) \\ &= f_i(\boldsymbol{\eta}) \end{aligned} \tag{7}\]

And residual variance is

\[ \begin{aligned} \mathrm{var}(y_i, y_j) &= E_\varepsilon((y_i - f_i(\boldsymbol{\eta}))(y_j - f_j(\boldsymbol{\eta}))) \\ &= E_\varepsilon\left(\left(\sum_{m=1}^{N} q'_{im}(\boldsymbol{\eta})\varepsilon_m\right)\left(\sum_{k=1}^{N} q'_{jk}(\boldsymbol{\eta})\varepsilon_k\right)\right) \\ &= E_\varepsilon\left(\sum_{m=1}^{N}\sum_{k=1}^{N} q'_{im}(\boldsymbol{\eta})\varepsilon_m q'_{jk}(\boldsymbol{\eta})\varepsilon_k\right) \\ & = \sum_{m=1}^{N}\sum_{k=1}^{N} q'_{im}(\boldsymbol{\eta}) q'_{jk}(\boldsymbol{\eta}) E(\varepsilon_m \varepsilon_k) = \sum_{m=1}^{N} q'_{im} q'_{jm} = \sum_{m=1}^{N} q'_{im} q_{mj} \end{aligned} \tag{8}\]

or

\[ \begin{aligned} \mathrm{var}(\mathbf{y}) &= \mathbf{Q}'\mathbf{Q} \\ &= \mathbf{V} \end{aligned} \tag{9}\]

To integrate over all \(\boldsymbol{\eta}\) and \(\boldsymbol{\varepsilon}\) and have an analytical solution, first define

\[ w'_{i(mk)} = \frac{\partial q'_{mi}}{\partial \eta_k} \tag{10}\]

as the series of \(N \times n\) matrices

\[ \mathbf{W}_i = \{w'_{i(mk)},\ m = 1\ \text{to}\ N,\ k = 1\ \text{to}\ n\} \tag{11}\]

and the \(N \times n\) matrix

\[ \mathbf{G} = \left\{g_{ik} = \frac{\partial f_i(\boldsymbol{\eta})}{\partial \eta_k},\ i = 1\ \text{to}\ N,\ k = 1\ \text{to}\ n\right\} \tag{12}\]

with \(1 \times n\) row vectors

\[ \mathbf{g}'_i = \left\{g_{ik} = \frac{\partial f_i(\boldsymbol{\eta})}{\partial \eta_k},\ k = 1\ \text{to}\ n\right\} \tag{13}\]

We now linearize by Taylor series expansion about \(\hat{\boldsymbol{\eta}}\) as follows:

\[ y_i = f_i(\hat{\boldsymbol{\eta}}) + \mathbf{g}'_i \boldsymbol{\eta} - \mathbf{g}'_i \hat{\boldsymbol{\eta}} + (\mathbf{q}'_i - \hat{\boldsymbol{\eta}}' \mathbf{W}'_i)\boldsymbol{\varepsilon} + \boldsymbol{\eta}' \mathbf{W}'_i \boldsymbol{\varepsilon} \tag{14}\]

If we now integrate over all \(\boldsymbol{\eta}\) and \(\boldsymbol{\varepsilon}\), we have the marginal density of \(\mathbf{y}\) with mean

\[ E_{\eta,\varepsilon}(y_i) = f_i(\hat{\boldsymbol{\eta}}) - \mathbf{g}'_i \hat{\boldsymbol{\eta}} \tag{15}\]

and total variance

\[ \mathrm{Var}(y_i, y_j) = E((y_i - f_i(\hat{\boldsymbol{\eta}}) + \mathbf{g}'_i \hat{\boldsymbol{\eta}})(y_j - f_j(\hat{\boldsymbol{\eta}}) + \mathbf{g}'_j \hat{\boldsymbol{\eta}})) = \tag{16}\]

\[ \begin{aligned} &E_{\eta,\varepsilon}\left((\mathbf{g}'_i \boldsymbol{\eta} + (\mathbf{q}'_i - \hat{\boldsymbol{\eta}}' \mathbf{W}'_i)\boldsymbol{\varepsilon} + \boldsymbol{\eta}' \mathbf{W}'_i \boldsymbol{\varepsilon})(\mathbf{g}'_j \boldsymbol{\eta} + (\mathbf{q}'_j - \hat{\boldsymbol{\eta}}' \mathbf{W}'_j)\boldsymbol{\varepsilon} + \boldsymbol{\eta}' \mathbf{W}'_j \boldsymbol{\varepsilon})\right) = \\ &\mathbf{g}'_i E_\eta(\boldsymbol{\eta}\boldsymbol{\eta}')\mathbf{g}_j + (\mathbf{q}'_i - \hat{\boldsymbol{\eta}}' \mathbf{W}'_i) E_\varepsilon(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}')(\mathbf{q}_j - \mathbf{W}_j \hat{\boldsymbol{\eta}}) + E_\eta(\boldsymbol{\eta}' \mathbf{W}'_i E_\varepsilon(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}')\mathbf{W}_j \boldsymbol{\eta}) \end{aligned} \tag{17}\]

But

\[ \begin{aligned} E(\boldsymbol{\eta}' \mathbf{W}'_i E(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}')\mathbf{W}_j \boldsymbol{\eta}) & = E(\boldsymbol{\eta}' \mathbf{W}'_i \mathbf{W}_j \boldsymbol{\eta}) \\ &= E\left(\sum_{m=1}^{n}\sum_{l=1}^{M}\sum_{k=1}^{n} \eta_k w'_{i(kl)} w_{j(lm)} \eta_m\right) \\ &= E\left(\sum_{m=1}^{n}\sum_{l=1}^{M}\sum_{k=1}^{n} \eta_k w'_{i(lk)} w_{j(lm)} \eta_m\right) \\ &= \sum_{m=1}^{n}\sum_{l=1}^{M}\sum_{k=1}^{n} w'_{i(lk)} E(\eta_k \eta_m) w_{j(lm)} \\ &= \sum_{l=1}^{M}\sum_{k=1}^{n}\sum_{m=1}^{n} w'_{i(lk)} \omega_{km} w_{j(lm)} \\ &= \mathrm{tr}(\mathbf{W}'_i \boldsymbol{\Omega} \mathbf{W}_j) = \mathrm{tr}(\mathbf{W}_j \boldsymbol{\Omega} \mathbf{W}'_i) = \mathrm{tr}(\mathbf{W}'_i \boldsymbol{\Omega} \mathbf{W}_j) \end{aligned} \tag{18}\]

And similarly,

\[ \hat{\boldsymbol{\eta}}' \mathbf{W}'_i \mathbf{W}_j \hat{\boldsymbol{\eta}} = \mathrm{tr}(\mathbf{W}_i \hat{\boldsymbol{\eta}} \hat{\boldsymbol{\eta}}' \mathbf{W}'_j) \tag{19}\]

So (CWRESI)

\[ \mathrm{Var}(y_i, y_j) = \mathbf{g}'_i \boldsymbol{\Omega} \mathbf{g}_j + \mathbf{q}'_i \mathbf{q}_j - (\mathbf{q}'_i \mathbf{W}_j + \mathbf{q}'_j \mathbf{W}_i) \hat{\boldsymbol{\eta}} + \mathrm{tr}\left[\mathbf{W}_i (\hat{\boldsymbol{\eta}} \hat{\boldsymbol{\eta}}' + \boldsymbol{\Omega}) \mathbf{W}'_j\right] \tag{20}\]

If \(\hat{\boldsymbol{\eta}} = \mathbf{0}\) (WRESI) then

\[ \mathrm{Var}(y_i, y_j) = \mathbf{g}'_i \boldsymbol{\Omega} \mathbf{g}_j + \mathbf{q}'_i \mathbf{q}_j + \mathrm{tr}\left[\mathbf{W}_i \boldsymbol{\Omega} \mathbf{W}'_j\right] \tag{21}\]

If \(\mathbf{W}_i = 0\), that is, interactive component is not taken into account, then (CWRES, (Hooker et al. 2007))

\[ \mathrm{Var}(y_i, y_j) = \mathbf{g}'_i \boldsymbol{\Omega} \mathbf{g}_j + \mathbf{q}'_i \mathbf{q}_j \tag{22}\]

In NONMEM 7.2, if $EST INTERACTION was specified prior to requesting $TABLE CWRES, then \(\mathbf{g}\) and \(\mathbf{q}\) are evaluated at \(\boldsymbol{\eta} = \hat{\boldsymbol{\eta}}\) in equation 22. If INTERACTION was not specified prior to requesting $TABLE CWRES, then \(\mathbf{g}\) and \(\mathbf{q}\) are evaluated at \(\boldsymbol{\eta} = 0\) in equation 22. In NONMEM 7.1.0 and 7.1.2, regardless of INTERACTION setting in a previous $EST statement, \(\mathbf{g}\) and \(\mathbf{q}\) are evaluated at \(\boldsymbol{\eta} = \hat{\boldsymbol{\eta}}\) in equation 22.

In NONMEM, the residual error is modeled as follows:

\[ y_i = f_i + h'_{i1} e_1 + h'_{i2} e_2 \ldots \tag{23}\]

for data point \(i\) of a particular subject, where \(e_k\) refers to the \(k\)th residual error component, that is in turn modeled to be normally distributed with variance

\[ E(e_{k_1} e_{k_2}) = \sigma_{k_1 k_2} \]

And

\[ h'_{ik} = \frac{\partial y_i}{\partial e_k} \]

Consider a problem where PK data are modeled with mixed homoscedastic error and heteroscedastic error, as is PD data, and there is a correlation between certain PK and PD data that are sampled at the same time. Such a correlation is indicated by the L2 variable listed in the data set. For such a problem, we could have:

\[ y_i = f_i + (2 - C_i) f_i e_1 + (2 - C_i) e_2 + (C_i - 1) f_i e_1 + (C_i - 1) e_2 \tag{24}\]

Where \(C = 1\) if the datum is PK, and \(C = 2\) if datum is PD. The Sigma matrix would be modeled as:

\[ \begin{pmatrix} \sigma_{11} & 0 & \sigma_{13} & 0 \\ 0 & \sigma_{22} & 0 & \sigma_{24} \\ \sigma_{31} & 0 & \sigma_{33} & 0 \\ 0 & \sigma_{42} & 0 & \sigma_{44} \end{pmatrix} \tag{25}\]

With correlation between certain paired PK,PD data, between their homoscedastic errors and heteroscedastic errors.

A grand matrix \(\mathbf{H}'\) is produced among all data points for a subject. Suppose a particular subject has four data points:

1. PK datum at time 1 hour
2. PD datum at time 1 hour
3. PK datum at time 2 hours
4. PD datum at time 3 horus

Data points 1 and 2 are coupled, and the others are not. An expanded \(4 \times 12\) matrix \(\mathbf{H}'\) is produced as follows:

\[ \begin{pmatrix} f_1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & f_2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & f_3 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & f_4 & 1 \end{pmatrix} \tag{26}\]

along with an expanded \(12 \times 12\) matrix \(\boldsymbol{\Sigma}_e\) of block diagonal form, consisting of the \(4 \times 4\) matrix \(\boldsymbol{\Sigma}\) duplicated three times along the diagonal. The placement of \(h'_{ik}\) in \(\mathbf{H}'\) determines if two data points are correlated within a shared block diagonal portion of \(\boldsymbol{\Sigma}_e\), as is the case with data points 1 and 2, or have separate block diagonals, and therefore uncorrelated, as with data points 3 and 4. In this example, \(M = 12\), \(N = 4\).

Defining matrix \(\boldsymbol{\Lambda}'\) as the lower triangular cholesky matrix to \(\boldsymbol{\Sigma}_e\) (earlier we defined \(\boldsymbol{\Lambda}\) as the lower triangular cholesky: we are redefining the nomenclature for this section for convenience):

\[ \boldsymbol{\Sigma}_e = \boldsymbol{\Lambda}'\boldsymbol{\Lambda} \tag{27}\]

It follows that we can construct

\[ \mathbf{Q}' = \mathbf{H}'\boldsymbol{\Lambda}' \tag{28}\]

\[ \mathbf{Q} = \boldsymbol{\Lambda}\mathbf{H} \tag{29}\]

Or, considering column vector \(\mathbf{h}_i\) of the \(i\)th column of \(\mathbf{H}\), then

\[ \mathbf{q}_i = \boldsymbol{\Lambda}\mathbf{h}_i \tag{30}\]

And the full intra-subject variance would be

\[ \begin{aligned} \mathbf{Q}'\mathbf{Q} &= \mathbf{V} \\ &= \mathbf{H}'\boldsymbol{\Lambda}'\boldsymbol{\Lambda}\mathbf{H} \\ &= \mathbf{H}'\boldsymbol{\Sigma}_e\mathbf{H} \end{aligned} \tag{31}\]

Furthermore, defining

\[ \mathbf{X}_i = \frac{\partial (\mathbf{h}_i)}{\partial \boldsymbol{\eta}} \tag{32}\]

then

\[ \begin{aligned} \mathbf{W}_i &= \frac{ \partial (\boldsymbol{\Lambda}\mathbf{h}_i) }{ \partial \boldsymbol{\eta} } \\ &= \boldsymbol{\Lambda} \frac{\partial (\mathbf{h}_i)}{\partial \boldsymbol{\eta}} \\ &= \boldsymbol{\Lambda}\mathbf{X}_i \end{aligned} \tag{33}\]

So,

\[ \mathrm{Var}(y_i, y_j) = \mathbf{g}'_i \boldsymbol{\Omega} \mathbf{g}_j + \mathbf{h}'_i \boldsymbol{\Sigma}_e \mathbf{h}_j - (\mathbf{h}'_i \boldsymbol{\Sigma}_e \mathbf{X}_j + \mathbf{h}'_j \boldsymbol{\Sigma}_e \mathbf{X}_i) \hat{\boldsymbol{\eta}} + \mathrm{tr}\left[(\hat{\boldsymbol{\eta}} \hat{\boldsymbol{\eta}}' + \boldsymbol{\Omega}) \mathbf{X}'_i \boldsymbol{\Sigma}_e \mathbf{X}_j\right] \tag{34}\]

An empirical method for evaluating the population weighted residual is to perform a Monte Carlo integration over all possible \(\boldsymbol{\eta}\). For a given subject, the expected population predicted values is

\[ \begin{aligned} \bar{\mathbf{f}}_{ij} &= E_\eta(\mathbf{f}(\boldsymbol{\eta})) \\ &= \int_{-\infty}^{\infty} \mathbf{f}(\boldsymbol{\eta}) p(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) d\boldsymbol{\eta} \end{aligned} \tag{35}\]

where

\[ p(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) = \frac{1}{\sqrt{2\pi |\boldsymbol{\Omega}|}} \exp\left(\frac{1}{2} \boldsymbol{\eta}' \boldsymbol{\Omega} \boldsymbol{\eta}\right) \tag{36}\]

The expected residual for an observed value is

\[ \begin{aligned} \bar{\mathbf{r}}_{ij} &= E_\eta(\mathbf{y} - \mathbf{f}(\boldsymbol{\eta})) \\ &= \int_{-\infty}^{\infty} (\mathbf{y} - \mathbf{f}(\boldsymbol{\eta})) p(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) d\boldsymbol{\eta} \\ &= \mathbf{y} - \bar{\mathbf{f}}_{ij} \end{aligned} \tag{37}\]

without using linearization methods on \(\boldsymbol{\eta}\). Now, since

\[ E_\varepsilon(\mathbf{y}) = \mathbf{f}(\boldsymbol{\eta}) \tag{38}\]

\[ E_\varepsilon((\mathbf{y} - \mathbf{f}(\boldsymbol{\eta}))(\mathbf{y} - \mathbf{f}(\boldsymbol{\eta}))') = \mathbf{Q}'\mathbf{Q} \tag{39}\]

Then the expected population variance is (without using linearization methods on eta):

\[ \begin{aligned} \mathbf{C}_{ij} &= E_{\eta,\varepsilon}((\mathbf{y} - \bar{\mathbf{f}}_{ij})(\mathbf{y} - \bar{\mathbf{f}}_{ij})') \\ &= E_\eta(E_\varepsilon((\mathbf{y} - \mathbf{f}(\boldsymbol{\eta}))(\mathbf{y} - \mathbf{f}(\boldsymbol{\eta}))')) + E_\eta((\mathbf{f}(\boldsymbol{\eta}) - \bar{\mathbf{f}}_{ij})(\mathbf{f}(\boldsymbol{\eta}) - \bar{\mathbf{f}}_{ij})') \\ &= \bar{\mathbf{V}}_{ij} + \mathbf{V}_f \end{aligned} \]

where

\[ \begin{aligned} \bar{\mathbf{V}}_{ij} &= E_\eta(\mathbf{Q}'(\boldsymbol{\eta})\mathbf{Q}(\boldsymbol{\eta})) \\ &= \int_{-\infty}^{\infty} \mathbf{Q}'(\boldsymbol{\eta})\mathbf{Q}(\boldsymbol{\eta}) p(\boldsymbol{\eta} \mid \mathbf{0}, \boldsymbol{\Omega}) d\boldsymbol{\eta} \end{aligned} \tag{40}\]

\[ \begin{aligned} \mathbf{V}_f &= E_\eta((\mathbf{f}(\boldsymbol{\eta}) - \bar{\mathbf{f}}_{ij})(\mathbf{f}(\boldsymbol{\eta}) - \bar{\mathbf{f}}_{ij})') \\ &= E_\eta(\mathbf{f}(\boldsymbol{\eta})\mathbf{f}(\boldsymbol{\eta})') - \bar{\mathbf{f}}_{ij}\bar{\mathbf{f}}'_{ij} \end{aligned} \tag{41}\]

To evaluate the expected weighted residual (EWRES),

\[ \mathbf{w}_{ij} = \mathbf{C}_{ij}^{-1/2} \bar{\mathbf{r}}_{ij} \]

where \(\mathbf{C}_{ij}^{-1/2}\) is the inverse square root of the expected population variance matrix.

An expected conditional (without interaction) weighted residual (ECWRES) can also be evaluated if we evaluate the intra-subject residual error at the conditional mean, such that

\[ \hat{\bar{\mathbf{V}}}_{ij} = \mathbf{Q}'(\hat{\boldsymbol{\eta}})\mathbf{Q}(\hat{\boldsymbol{\eta}}) \tag{42}\]

if $EST INTERACTION is specified followed by $TABLE ECWRES. But all other components are Monte Carlo integrated:

\[ \mathbf{C}_{ij} = \hat{\bar{\mathbf{V}}}_{ij} + \mathbf{V}_f \tag{43}\]

As of NONMEM 7.2, if INTERACTION in $EST was not specified, followed by $TABLE ECWRES, then

\[ \begin{aligned} \mathbf{V}_0 &= \mathbf{Q}'(\boldsymbol{\eta} \\ &= 0)\mathbf{Q}(\boldsymbol{\eta} \\ &= 0) \end{aligned} \tag{44}\]

and all other components are Monte Carlo integrated:

\[ \mathbf{C}_{ij} = \mathbf{V}_0 + \mathbf{V}_f \tag{45}\]

In NONMEM 7.1.0 and 7.1.2, regardless of INTERACTION setting from the previous $EST command, \(\hat{\bar{\mathbf{V}}}_{ij}\) is used to evaluate ECWRES.

1 NPDE

The NPDE is the normalized prediction distribution error (reference (Comets et al. 2008): takes into account within-subject correlations), also a Monte Carlo assessed diagnostic item. For the \(k\)th simulated vector of data \(\mathbf{y}_{sk}\):

\[ \mathbf{s}_{ijk} = \mathbf{y}_{sk} - \bar{\mathbf{f}}_{ij} \tag{46}\]

its decorrelated residual vector is calculated:

\[ \mathbf{w}_{sijk} = \mathbf{C}_{ij}^{-1/2} \mathbf{s}_{ijk} \tag{47}\]

and compared against the decorrelated residual vector of observed values \(\mathbf{w}_{ij}\) such that

\[ \mathbf{u} = \frac{1}{K} \sum_{k=1}^{K} \delta(\mathbf{w}_{ij} - \mathbf{w}_{sijk}) \tag{48}\]

For \(K\) random samples, where

\[ \delta(x) = \begin{cases} 1 & \text{for } x \geq 0 \\ 0 & \text{for } x < 0 \end{cases} \]

For each element in the vector. Then, an inverse normal distribution transformation is performed:

\[ \mathbf{w}_{\mathrm{npde}} = \Phi^{-1}(\mathbf{u}) \tag{49}\]

2 NPD

The NPD is the correlated normalized prediction distribution error (reference (Brendel et al. 2006): does not take into account within-subject correlations), also a Monte Carlo assessed diagnostic item. For each vector of data \(\mathbf{y}\):

\[ \mathbf{r}_{\eta_k} = \mathbf{V}(\boldsymbol{\eta}_k)^{-1/2} (\mathbf{y} - \mathbf{f}(\boldsymbol{\eta}_k)) \tag{50}\]

These are then averaged over all the random samples:

\[ \mathbf{u}_c = \frac{1}{K} \sum_{k=1}^{K} \Phi(\mathbf{r}_{\eta_k}) \tag{51}\]

Then, an inverse normal distribution transformation is performed:

\[ \mathbf{w}_{\mathrm{npdec}} = \Phi^{-1}(\mathbf{u}_c) \tag{52}\]

Brendel, Karl, Emmanuelle Comets, Celine Laffont, Christian Laveille, and France Mentré. 2006. “Metrics for External Model Evaluation with an Application to the Population Pharmacokinetics of Gliclazide.” Pharmaceutical Research 23 (9): 2036–49. https://doi.org/10.1007/s11095-006-9067-5.

Comets, Emmanuelle, Karl Brendel, and France Mentré. 2008. “Computing Normalized Prediction Distribution Errors to Evaluate Nonlinear Mixed Effects Models: The npde Add-on Package for R.” Computer Methods and Programs in Biomedicine 90 (2): 154–66. https://doi.org/10.1016/j.cmpb.2007.12.002.

Hooker, Andrew C., Christine E. Staatz, and Mats O. Karlsson. 2007. “Conditional Weighted Residuals (CWRES): A Model Diagnostic for the FOCE Method.” Pharmaceutical Research 24: 2187–97. https://doi.org/10.1007/s11095-007-9361-x.