34 Model Building

1 What This Chapter is About

In this chapter, the simple phenobarbital example begun in Chapter 2 will be continued to illustrate how NONMEM is used to build a model for population data. The topic of model building, diagnosis and verification is a large one. This chapter can only give a very abbreviated example.

2 The Stages of Model Building

To analyze a population data set and build a model for it, one must proceed in logical stages. There are five stages, and their relationship to one another is presented diagrammatically in figure 11.1. One begins by checking the data. One then tries to find an adequate model incorporating the fixed effects; then an adequate model incorporating the random effects and describing random inter- and intra-individual variability. After a reasonably complete model is found, attempts are made to refine it, and finally, if desired, the various parts of the models (which often, in effect, simply assert the existence of certain relationships between independent variables and the dependent variable) can be subject to formal hypothesis tests, as described in Chapter 5. (However, it is well known by statisticians that formal hypothesis testing undertaken after model building is just an approximation for the type of hypothesis testing described in textbooks, which assumes that the model is the correct model.)

Figure 11.1. Stages in model building.

The five stages in order are:

Check-out — Fix data problems
Build Structural Model — Add structural features as indicated
Build Statistical Model — Add statistical model features as indicated
Refine Model — Remove/change suspicious features
Formally Test Model — Test major features by deleting one by one

Each stage uses NONMEM tables, scatterplots, estimation, and the covariance step as appropriate.

3 Check-out — Index Plots

The goal of this stage is to assure that the data are correct. There is no point to trying to model the data when gross errors are present. Most gross errors are encoding errors that cause certain values to be different from the intended value by a considerable amount (for example, a misplaced decimal point changes a value by a factor of 10), so that graphical display of the data is usually adequate to detect these. No numerical or statistical approaches are needed. Indeed, they are not usually useful, even for more subtle errors, as such errors cannot easily be detected by any means (how is a 10% error to be distinguished from inter- or intra-individual variability?).

To detect gross errors, then, one makes scatterplots of different data item types vs individuals’ identification numbers (i.e. the ID data item, or, if the values of this data item are arbitrary, another data item that identifies patients using sequential integer values; call this the sequence data item: SEQ). Such plots (of one data item versus ID or SEQ) are called here index plots, and are quite useful for revealing the structure of the data, as will be noted below, as well as for finding gross errors.

If NONMEM is used to make index plots, it will also be useful to implement the Tables Step, so that if a problem is noted in a scatterplot, one can refer to the table to try to find the datum that might account for the problem. To run NONMEM some model must be specified, even if all that is desired is an index plot. In such case, it makes little difference what model is used. It is easiest and useful to:

start with a simple ADVAN that is likely to provide at least a roughly satisfactory fit,
set each PK parameter to a (different) element of $\theta$,
use only one $\eta$ variable, modifying the scale parameter only, and one $\varepsilon$ variable, and
use roughly reasonable fixed initial estimates.

For the phenobarbital example, one might use ADVAN1 with $K = \theta_1$, $V = \theta_2 + \eta_1$, and $y = f + f\varepsilon_1$. Initial estimates might be: $K = .0057 \text{ hr}^{-1}$ (half-life = 5 days, a typical value for adults); $V = 1.44$ L. (The first patient has a concentration of 17.3 mg/L some few hours after an initial loading dose of 25 mg; $1.44 \times 17.3 = 25$); $\omega = .25$ (50% variability); $\sigma^2 = .04$ (20% variability).

Figures 11.2 and 11.3 show index plots that might be seen in a check-out run (gross errors have been added). In figure 11.2, DV is plotted vs ID, and a gross error occurring at about patient #13 is seen (an observation of about 24 mcg/ml was erroneously recorded as 240 mcg/ml). In figure 11.3, AMT is plotted vs ID, and patient #3 appears to have a grossly erroneous value (again, a decimal point error; a dose of 18 was misrecorded as 180).

Figure 11.2. Index plot of DV vs ID, showing a gross error at patient #13.

Figure 11.3. Index plot of AMT vs ID, showing a gross error at patient #3.

When such errors are found, they should be corrected (or, if confirmation from original data sources is needed, the affected records temporarily excluded using the IGNORE option of $DATA) and the check-out run re-done.

After correcting the errors, index plots of covariates (e.g. weight, APGAR score) against ID can help reveal the structure of the data. Figures 11.4 and 11.5 show corrected index plots of WT vs ID and APGR vs ID.

Figure 11.4. Corrected index plot of weight (WT) vs ID.

Figure 11.5. Corrected index plot of Apgar score (APGR) vs ID.

Note from figures 11.4 and 11.5 that Apgar scores covary with weight. In particular, a plot of WT vs APGR (figure 11.6) shows that the babies with the lowest weights tend to have the lowest Apgar scores.

Figure 11.6. Scatterplot of weight (WT) vs Apgar score (APGR).

4 Building the Structural Part of the Model

4.1 A General Approach

Once the data is clean, model building can begin. It is best to build incrementally. The primary purpose of the data items in $x$ is to help predict $y$. Start simple. Assume that PK parameters do not depend on concomitant variables (the $z$’s), and that $\varepsilon$ models the residual error for the typical individual. This first model is often called the minimal model. As problems with the fit of the minimal model arise, they should be addressed by adding additional structural features (more $z$’s). As problems with the statistical model arise, they should be addressed by adding additional variability terms ($\eta$’s, $\varepsilon$’s). This is the subject of the present section.

4.2 The Minimal Model

A minimal structural model for the phenobarbital data is:

\[Cl = \theta_1 + \eta_1 \tag{1}\]

\[V = \theta_2 + \eta_2 \tag{2}\]

\[y = f (1 + \varepsilon) \tag{3}\]

The NM-TRAN input for this model is shown in figure 11.7 (the CCV error model is chosen here).

$PROBLEM PHENOBARB
$INPUT ID TIME AMT WT APGR DV
$DATA PHENO
$SUBROUTINE ADVAN1
$PK
  TVCL=THETA(1)
  CL=TVCL+ETA(1)
  TVVD=THETA(2)
  V=TVVD+ETA(2)
  K=CL/V
  S1=V
$ERROR
  Y=F*(1+ERR(1))
$THETA (0,.0105) (0,1.05)
$OMEGA .25 .25
$SIGMA .04
$ESTIMATION
$TABLE       ID TIME AMT WT APGR
$SCATTERPLOT PRED VS DV UNIT
$SCATTERPLOT RES VS (PRED,WT,APGR)
$SCATTERPLOT WRES VS (PRED,WT,APGR)

4.3 Use of Constraints

It is important to realize that constraints on elements of $\theta$ or $\Omega$ may be part of a model. For example, constraining clearance to be positive is a modelling choice. One might implement this constraint in NONMEM using a lower bound on the $THETA record, and this would assure that the estimate of clearance will be positive. It may not be necessary to do this; even without the lower bound, the data might clearly force the estimate to be positive.

Often, however, analysts will constrain the range of a parameter in the belief that doing so will shorten computing time or stabilize the search for the minimum of the objective function. While this benefit may be gained, the data may force the parameter estimate to the constraint boundary even though this boundary may not, in fact, represent a true modeling choice. In this case the proper action is to relax the constraint and rerun the problem. To do otherwise, and leave the parameter estimate to be the boundary value, implies that at the outset the user assumes that the parameter must be within the boundary and elevates the constraint to the status of a modeling choice. If an estimate lies on a boundary, NONMEM will print a warning message (along with the standard message regarding the status of the termination of the Estimation Step). The reader is cautioned to look for such a message, and in general, it is a good idea to check the values of the final estimates against the boundary values. Alternatively, the implementation of constraints that are not intended to represent modeling choices might be used cautiously and only if they really seem necessary to stabilize a search.

4.4 Diagnostic Tools

4.4.1 Plot of DV vs `PRED`

Most useful diagnostic tools are graphical. For an overall sense of the fit, a useful diagnostic plot is DV vs PRED. When there are substantial and systematic deviations from the line of identity, this plot suggests that there are problems with the fit, but it does not suggest what exactly these problems might be or what to do about them. This plot for the fit of the phenobarbital data to model (1) is seen in figure 11.7.

Figure 11.7. Predictions from fit of model (1) to phenobarbital data vs observations themselves. The line of identity shows where the points should, ideally, fall.

Figure 11.7 reveals that there is a group of points where the observation is much lower than the prediction. To begin to determine why this is so, it will be useful to look at residual plots. Such plots are the basis of the most important of the diagnostic tools.

4.4.2 Residual Plots

As mentioned in Chapter 2, a residual is the difference between an observation and its prediction. The prediction in this case (the same prediction as denoted by PRED) is the population prediction, i.e. the prediction for the typical individual having the given values for all the concomitant variables.

With population data, weighted residuals are often more informative than (plain) residuals. The weighted residuals for an individual are formed by transforming the individual’s residuals so that under the population model, and assuming the true values of the population parameters are given by the estimates of those parameters, all weighted residuals have unit variance and are uncorrelated. Weighted residuals are more informative for several reasons. First, since they have unit variance, or what is the same, unit standard deviation, “large” weighted residuals are those with absolute values greater than 3 or so. Second, loosely speaking, although plain residuals remove the structural model from the data, allowing one to see what part of the data is not (yet) modeled, they do not remove the statistical model (formally, they are still correlated). Weighted residuals have both models ‘removed’ so any pattern in these is definitely not accounted for by the current model. This provides a more secure basis for future model building choices.

4.4.2.1 Index Plots of Residuals

Figure 11.8 is an index plot of residuals, which is a useful plot when combined with index plots of other data items. One can look for an association between unusual residuals and values of another data item. E.g. Are the largest discrepancies between model and data associated with certain (possibly extreme) values of the data item?

Figure 11.8. Index plot of residuals from fit of model (1) to phenobarbital data.

4.4.3 Plot of WRES vs Independent Variable

A plot of WRES vs an independent variable (e.g. weight) can suggest what independent variable might contribute to the structural model. Figure 11.9 shows WRES vs WT.

Figure 11.9. Weighted residuals vs weight from fit of model (1) to phenobarbital data.

The pattern in figure 11.9 is similar to that in figure 2.11 of Chapter 2: observations in low-weight infants tend to be underpredicted, while those in high-weight infants tend to be overpredicted. This suggests modeling CL and V as linear functions of WT.

4.5 Judging Goodness of Fit

4.5.1 A Global Measure — Change in the Objective Function

A particular model’s fit to the data results in a minimum value of the objective function, which is denoted by $l$. When a model is refined, e.g. by adding additional parameters, the minimum value of the objective function decreases. The magnitude of this decrease can be used to judge whether the refinement is worthwhile. Statistical theory (see Chapter 5, Section 5.2) says that a decrease greater than about 3.84 (equal to the 95th percentile of a chi-squared distribution with 1 degree of freedom) represents a significant improvement at the 5% level, when one parameter is added to the reduced model.

4.5.2 Decrease in Unexplained Variability

Unexplained variability is quantified by the estimates of $\omega^2$ and $\sigma^2$ (or elements of $\Omega$ and $\Sigma$). When refining a model by adding a new parameter (e.g., adding an effect of weight on CL), we would expect these variance estimates to decrease (because the new parameter is “explaining” more of the variability). Significant decreases are evidence that the new model is an improvement.

4.5.3 Improvement in Plots

Improvement in the fit can also be judged by examining diagnostic plots after refitting the new model. If systematic deviations in the plot of DV vs PRED disappear, or if patterns in residual plots are reduced, the refinement is worthwhile.

4.6 Using the Tools: Further Improvement

4.6.1 An Additional Effect of WT

Based on the pattern in figure 11.9, a new model incorporating WT is tried:

\[Cl = (\theta_1 + \theta_3 WT) + \eta_1 \tag{4}\]

\[V = (\theta_2 + \theta_4 WT) + \eta_2 \tag{5}\]

\[y = f(1 + \varepsilon) \tag{6}\]

This is the model of Chapter 2, figure 2.12. When this model is fit to the data, the objective function drops by over 100 points. (See Chapter 2, figure 2.13 for the results.)

However, when examining the plot of WRES vs WT for this model, the pattern seen in figure 11.9 has mostly disappeared. This is an indication that the new model adequately accounts for the influence of weight on the pharmacokinetic parameters.

As noted in Chapter 2, the final estimate of $\theta_1$ is very small (essentially 0). This suggests that $\theta_1$ should be set to 0, i.e., that clearance is simply proportional to weight rather than being a linear function with a non-zero intercept.

4.6.2 The Effect of APGR

After accounting for weight, a plot of WRES vs APGR can be made to see if Apgar score has additional effect. Figure 11.15 shows this plot.

Figure 11.15. Plot of weighted residuals vs Apgar score for fit of model (11.3) to phenobarbital data.

There is a weak suggestion from figure 11.15 that for APGR less than 3, the weighted residuals tend to be negative. Accordingly, a new model (7) can be proposed. One might define an indicator variable LA (Low Apgar) equal to 1 when APGR < 3, and 0 otherwise:

\[V = \theta_2 + \theta_4 WT + \theta_5 LA \cdot WT + \eta_2 \tag{7}\]

When this model is fit to the data, $\theta_5$ is estimated to be 1.18, implying that indeed, the volumes of distribution for infants with Apgar scores less than 3 are typically 18% higher than those of infants (of the same weight) with higher Apgar scores. The measures of goodness of fit do not change much: the objective function drops only slightly, while the variances of $\eta$ actually increase a few percent. Inspection of figure 11.6 suggests a reason for this: note that only 5 distinct individuals (separate symbols) have Apgar scores less than 3. There is simply not very much information about babies with low Apgar scores in this data set.

Figure 11.16. Plot of weighted residuals vs Apgar score for fit of model (7) to phenobarbital data.

5 Building the Statistical Model

5.1 Judging Among Alternatives

NONMEM can provide estimates of the $\eta$ variables for each individual (see Chapter 12, Sections 4.11-4.13). Plots of the estimated interindividual differences, which can be regarded as interindividual residuals, can be obtained. Plots of these residuals (associated with a particular parameter) can provide further help in building the structural part of the parameter model. Moreover, a plot of the estimated $\eta$ for one parameter vs. the estimated $\eta$ for another can reveal any correlation between them.

Lastly, help can be provided in the selection of a model for intraindividual error. Predictions for each individual based on conditional estimates of $\eta$ can be computed, and the residuals can be plotted versus the predictions. (This requires advanced techniques.) Again, if intraindividual errors are modeled with the additive model, and the resulting residuals suggest a CCV pattern, one would refine the model by choosing the CCV model for intraindividual errors.

When the interindividual error model is varied, the effect is primarily to influence the estimates of the variances of the corresponding $\eta$ variables. Estimates of population parameters (the $\theta$’s) are particularly precise. Sometimes, however, the variability in the random effects is an important part of the model. In these cases, care should be taken to choose the appropriate random effects model. This, however, may not be easy because it is an unfortunate, but generally true fact that statistical models are much harder to select on the basis of data than are structural models.

5.1.1 Unexplained Variability

When the statistical model is developed, a new $\eta$ variable may be added, or an old $\eta$ variable used differently. Then differences in the $\omega$s between models cannot really be used to judge improvement. (As discussed above, when fixed effects are added to an existing model, decreases in estimates of $\omega^2$ do indicate improvement.)

However, there is one sure sign that too many $\eta$s are in the model: NONMEM may estimate one or more of the $\omega$s to be zero, or very nearly zero. This can be disconcerting, particularly if the $\eta$ variable is the only such variable affecting Volume, for instance, as then this estimate seems to suggest that with respect to Volume, there is no interindividual variability, however. Rather, assuming the $\eta$ affecting Volume is the one most recently added, it indicates that the data do not have enough information to estimate this variance. The $\eta$ should be removed from the model.

5.1.2 Residual Plots

If the intraindividual error model is the correct one, there should be no pattern to the residuals plotted vs the predictions (a pattern with the predictions would indicate a defect in the structural model). This can be illustrated using the phenobarbital data. Figure 11.17 shows the plot of weighted residuals vs predictions for fit of model (7) to phenobarbital data, using a proportional (CCV) error model.

Figure 11.17. Plot of weighted residuals vs predictions for fit of model (7) to phenobarbital data (proportional error model).

Only the $ERROR statements of the control stream change to use an additive error model:

$ERROR
  Y=F+ERR(1)

Figure 11.18. Plot of weighted residuals vs predictions for fit of model (11.5) to phenobarbital data (additive error model).

Although the plots do not differ greatly, there is a small suggestion in figure 11.18 that a CCV error model is more appropriate, because the spread of weighted residuals appears to increase somewhat with increasing predictions. That this impression is valid is suggested also by the increase in objective function value.

6 Refine Model

Refinement is an iterative process of testing and adjusting. It differs from the Building stages in that refinement tests whether the model can be simplified or whether parameters can be reduced, while still retaining as good a fit to the data as possible.

6.1 Use of Standard Errors and Confidence Intervals

The basic technique at this stage is to run the Covariance Step and examine the standard errors of the parameter estimates. If a parameter has a very large standard error relative to its estimate, it may be removing the parameter (or a related one) without significantly degrading the fit.

Figures 11.19 and 11.20 show two pages of NONMEM output from a run fitting model (7) to the phenobarbital data, and implementing the Covariance Step. Figure 11.19 shows the standard errors. Figure 11.20 shows the correlation matrix.

7 Testing the Model

The final stage is formal testing of the model. This involves performing likelihood ratio tests (or other tests) to verify that each part of the model is actually supported by the data. This is done by, for each model parameter whose significance is in question, fitting a reduced model in which that parameter is removed (or fixed to 0). The difference in the objective function values between the full model and each reduced model can then be compared to the appropriate chi-squared percentile.

The importance of formal testing is that it provides an objective basis for the model, not just a subjective one. However, as noted at the start of this chapter, formal hypothesis testing after model building is problematic. The selected model may have been chosen (by inspection of plots) precisely because it fits the data well. Testing parts of this model against alternative, reduced models then tends to yield confirmatory results. A more rigorous approach would be to use the current data to build the model, and then to test the model’s significance using a different (independent) data set. This is often referred to as validation.

Model building and refinement is a difficult and iterative process. The tools NONMEM provides are valuable, but should be complemented by physiological knowledge of the drug being studied and reasonable prior expectations about the model.