PROC NLMIXED: One-Compartment Model with Pharmacokinetic Data :: SAS/STAT(R) 9.2 User's Guide, Second Edition

The NLMIXED Procedure

Example 61.1 One-Compartment Model with Pharmacokinetic Data

A popular application of nonlinear mixed models is in the field of pharmacokinetics, which studies how a drug disperses through a living individual. This example considers the theophylline data from Pinheiro and Bates (1995). Serum concentrations of the drug theophylline are measured in 12 subjects over a 25-hour period after oral administration. The data are as follows.

   data theoph;
     input subject time conc dose wt;
     datalines; 
    1  0.00  0.74 4.02 79.6
    1  0.25  2.84 4.02 79.6
    1  0.57  6.57 4.02 79.6
    1  1.12 10.50 4.02 79.6
    1  2.02  9.66 4.02 79.6
    1  3.82  8.58 4.02 79.6
    1  5.10  8.36 4.02 79.6
    1  7.03  7.47 4.02 79.6
    1  9.05  6.89 4.02 79.6
    1 12.12  5.94 4.02 79.6
    1 24.37  3.28 4.02 79.6
    2  0.00  0.00 4.40 72.4
    2  0.27  1.72 4.40 72.4
    2  0.52  7.91 4.40 72.4
    2  1.00  8.31 4.40 72.4
    2  1.92  8.33 4.40 72.4
    2  3.50  6.85 4.40 72.4
    2  5.02  6.08 4.40 72.4
    2  7.03  5.40 4.40 72.4
    2  9.00  4.55 4.40 72.4
    2 12.00  3.01 4.40 72.4
    2 24.30  0.90 4.40 72.4
    3  0.00  0.00 4.53 70.5
    3  0.27  4.40 4.53 70.5
    3  0.58  6.90 4.53 70.5
    3  1.02  8.20 4.53 70.5
    3  2.02  7.80 4.53 70.5
    3  3.62  7.50 4.53 70.5
    3  5.08  6.20 4.53 70.5
    3  7.07  5.30 4.53 70.5
    3  9.00  4.90 4.53 70.5
    3 12.15  3.70 4.53 70.5
    3 24.17  1.05 4.53 70.5
    4  0.00  0.00 4.40 72.7
    4  0.35  1.89 4.40 72.7
    4  0.60  4.60 4.40 72.7
    4  1.07  8.60 4.40 72.7
    4  2.13  8.38 4.40 72.7
    4  3.50  7.54 4.40 72.7
    4  5.02  6.88 4.40 72.7
    4  7.02  5.78 4.40 72.7
    4  9.02  5.33 4.40 72.7
    4 11.98  4.19 4.40 72.7
    4 24.65  1.15 4.40 72.7
    5  0.00  0.00 5.86 54.6
    5  0.30  2.02 5.86 54.6
    5  0.52  5.63 5.86 54.6
    5  1.00 11.40 5.86 54.6
    5  2.02  9.33 5.86 54.6
    5  3.50  8.74 5.86 54.6
    5  5.02  7.56 5.86 54.6
    5  7.02  7.09 5.86 54.6
    5  9.10  5.90 5.86 54.6
    5 12.00  4.37 5.86 54.6
    5 24.35  1.57 5.86 54.6
    6  0.00  0.00 4.00 80.0
    6  0.27  1.29 4.00 80.0
    6  0.58  3.08 4.00 80.0
    6  1.15  6.44 4.00 80.0
    6  2.03  6.32 4.00 80.0
    6  3.57  5.53 4.00 80.0
    6  5.00  4.94 4.00 80.0
    6  7.00  4.02 4.00 80.0
    6  9.22  3.46 4.00 80.0
    6 12.10  2.78 4.00 80.0
    6 23.85  0.92 4.00 80.0
    7  0.00  0.15 4.95 64.6
    7  0.25  0.85 4.95 64.6
    7  0.50  2.35 4.95 64.6
    7  1.02  5.02 4.95 64.6
    7  2.02  6.58 4.95 64.6
    7  3.48  7.09 4.95 64.6
    7  5.00  6.66 4.95 64.6
    7  6.98  5.25 4.95 64.6
    7  9.00  4.39 4.95 64.6
    7 12.05  3.53 4.95 64.6
    7 24.22  1.15 4.95 64.6
    8  0.00  0.00 4.53 70.5
    8  0.25  3.05 4.53 70.5
    8  0.52  3.05 4.53 70.5
    8  0.98  7.31 4.53 70.5
    8  2.02  7.56 4.53 70.5
    8  3.53  6.59 4.53 70.5
    8  5.05  5.88 4.53 70.5
    8  7.15  4.73 4.53 70.5
    8  9.07  4.57 4.53 70.5
    8 12.10  3.00 4.53 70.5
    8 24.12  1.25 4.53 70.5
    9  0.00  0.00 3.10 86.4
    9  0.30  7.37 3.10 86.4
    9  0.63  9.03 3.10 86.4
    9  1.05  7.14 3.10 86.4
    9  2.02  6.33 3.10 86.4
    9  3.53  5.66 3.10 86.4
    9  5.02  5.67 3.10 86.4
    9  7.17  4.24 3.10 86.4
    9  8.80  4.11 3.10 86.4
    9 11.60  3.16 3.10 86.4
    9 24.43  1.12 3.10 86.4
   10  0.00  0.24 5.50 58.2
   10  0.37  2.89 5.50 58.2
   10  0.77  5.22 5.50 58.2
   10  1.02  6.41 5.50 58.2
   10  2.05  7.83 5.50 58.2
   10  3.55 10.21 5.50 58.2
   10  5.05  9.18 5.50 58.2
   10  7.08  8.02 5.50 58.2
   10  9.38  7.14 5.50 58.2
   10 12.10  5.68 5.50 58.2
   10 23.70  2.42 5.50 58.2
   11  0.00  0.00 4.92 65.0
   11  0.25  4.86 4.92 65.0
   11  0.50  7.24 4.92 65.0
   11  0.98  8.00 4.92 65.0
   11  1.98  6.81 4.92 65.0
   11  3.60  5.87 4.92 65.0
   11  5.02  5.22 4.92 65.0
   11  7.03  4.45 4.92 65.0
   11  9.03  3.62 4.92 65.0
   11 12.12  2.69 4.92 65.0
   11 24.08  0.86 4.92 65.0
   12  0.00  0.00 5.30 60.5
   12  0.25  1.25 5.30 60.5
   12  0.50  3.96 5.30 60.5
   12  1.00  7.82 5.30 60.5
   12  2.00  9.72 5.30 60.5
   12  3.52  9.75 5.30 60.5
   12  5.07  8.57 5.30 60.5
   12  7.07  6.59 5.30 60.5
   12  9.03  6.11 5.30 60.5
   12 12.05  4.57 5.30 60.5
   12 24.15  1.17 5.30 60.5
   ;

Pinheiro and Bates (1995) consider the following first-order compartment model for these data:

$\text{[math]}$

where $\text{[math]}$ is the observed concentration of the $\text{[math]}$ th subject at time $\text{[math]}$ , $\text{[math]}$ is the dose of theophylline, $\text{[math]}$ is the elimination rate constant for subject $\text{[math]}$ , $\text{[math]}$ is the absorption rate constant for subject $\text{[math]}$ , $\text{[math]}$ is the clearance for subject $\text{[math]}$ , and $\text{[math]}$ are normal errors. To allow for random variability between subjects, they assume

	$\text{[math]}$	$\text{[math]}$
	$\text{[math]}$	$\text{[math]}$
	$\text{[math]}$	$\text{[math]}$

where the $\text{[math]}$ s denote fixed-effects parameters and the $\text{[math]}$ s denote random-effects parameters with an unknown covariance matrix.

The PROC NLMIXED statements to fit this model are as follows:

   proc nlmixed data=theoph;
      parms beta1=-3.22 beta2=0.47 beta3=-2.45 
            s2b1 =0.03  cb12 =0    s2b2 =0.4 s2=0.5;
      cl   = exp(beta1 + b1);
      ka   = exp(beta2 + b2);
      ke   = exp(beta3);
      pred = dose*ke*ka*(exp(-ke*time)-exp(-ka*time))/cl/(ka-ke);
      model conc ~ normal(pred,s2);
      random b1 b2 ~ normal([0,0],[s2b1,cb12,s2b2]) subject=subject;
   run;

The PARMS statement specifies starting values for the three $\text{[math]}$ s and four variance-covariance parameters. The clearance and rate constants are defined using SAS programming statements, and the conditional model for the data is defined to be normal with mean pred and variance s2. The two random effects are b1 and b2, and their joint distribution is defined in the RANDOM statement. Brackets are used in defining their mean vector (two zeros) and the lower triangle of their variance-covariance matrix (a general $\text{[math]}$ matrix). The SUBJECT= variable is subject.

The results from this analysis are as follows.

Output 61.1.1 Model Specification for One-Compartment Model

The NLMIXED Procedure

Specifications
Data Set	WORK.THEOPH
Dependent Variable	conc
Distribution for Dependent Variable	Normal
Random Effects	b1 b2
Distribution for Random Effects	Normal
Subject Variable	subject
Optimization Technique	Dual Quasi-Newton
Integration Method	Adaptive Gaussian Quadrature

The "Specifications" table lists the setup of the model (Output 61.1.1). The "Dimensions" table indicates that there are 132 observations, 12 subjects, and 7 parameters. PROC NLMIXED selects 5 quadrature points for each random effect, producing a total grid of 25 points over which quadrature is performed (Output 61.1.2).

Output 61.1.2 Dimensions Table for One-Compartment Model

Dimensions
Observations Used	132
Observations Not Used	0
Total Observations	132
Subjects	12
Max Obs Per Subject	11
Parameters	7
Quadrature Points	5

The "Parameters" table lists the 7 parameters, their starting values, and the initial evaluation of the negative log likelihood using adaptive Gaussian quadrature (Output 61.1.3). The "Iteration History" table indicates that 10 steps are required for the dual quasi-Newton algorithm to achieve convergence.

Output 61.1.3 Starting Values and Iteration History

Parameters
beta1	beta2	beta3	s2b1	cb12	s2b2	s2	NegLogLike
-3.22	0.47	-2.45	0.03	0	0.4	0.5	177.789945

Iteration History
Iter	Calls	NegLogLike	Diff	MaxGrad	Slope
1	5	177.776248	0.013697	2.873367	-63.0744
2	8	177.7643	0.011948	1.698144	-4.75239
3	10	177.757264	0.007036	1.297439	-1.97311
4	12	177.755688	0.001576	1.441408	-0.49772
5	14	177.7467	0.008988	1.132279	-0.8223
6	17	177.746401	0.000299	0.831293	-0.00244
7	19	177.746318	0.000083	0.724198	-0.00789
8	21	177.74574	0.000578	0.180018	-0.00583
9	23	177.745736	3.88E-6	0.017958	-8.25E-6
10	25	177.745736	3.222E-8	0.000143	-6.51E-8

NOTE: GCONV convergence criterion satisfied.

Output 61.1.4 Fit Statistics for One-Compartment Model

Fit Statistics
-2 Log Likelihood	355.5
AIC (smaller is better)	369.5
AICC (smaller is better)	370.4
BIC (smaller is better)	372.9

The "Fit Statistics" table lists the final optimized values of the log-likelihood function and information criteria in the "smaller is better" form (Output 61.1.4).

Output 61.1.5 Parameter Estimates for One-Compartment Model

Parameter Estimates
Parameter	Estimate	Standard Error	DF	t Value	Pr > \|t\|	Alpha	Lower	Upper	Gradient
beta1	-3.2268	0.05950	10	-54.23	<.0001	0.05	-3.3594	-3.0942	-0.00009
beta2	0.4806	0.1989	10	2.42	0.0363	0.05	0.03745	0.9238	3.645E-7
beta3	-2.4592	0.05126	10	-47.97	<.0001	0.05	-2.5734	-2.3449	0.000039
s2b1	0.02803	0.01221	10	2.30	0.0445	0.05	0.000833	0.05523	-0.00014
cb12	-0.00127	0.03404	10	-0.04	0.9710	0.05	-0.07712	0.07458	-0.00007
s2b2	0.4331	0.2005	10	2.16	0.0560	0.05	-0.01353	0.8798	-6.98E-6
s2	0.5016	0.06837	10	7.34	<.0001	0.05	0.3493	0.6540	6.133E-6

The "Parameter Estimates" table contains the maximum likelihood estimates of the parameters (Output 61.1.5). Both s2b1 and s2b2 are marginally significant, indicating between-subject variability in the clearances and absorption rate constants, respectively. There does not appear to be a significant covariance between them, as seen by the estimate of cb12.

The estimates of $\text{[math]}$ , $\text{[math]}$ , and $\text{[math]}$ are close to the adaptive quadrature estimates listed in Table 3 of Pinheiro and Bates (1995). However, Pinheiro and Bates use a Cholesky-root parameterization for the random-effects variance matrix and a logarithmic parameterization for the residual variance. The PROC NLMIXED statements using their parameterization are as follows, and results are similar.

   proc nlmixed data=theoph;
      parms ll1=-1.5 l2=0 ll3=-0.1 beta1=-3 beta2=0.5 beta3=-2.5 
         ls2=-0.7;
      s2   = exp(ls2);
      l1   = exp(ll1);
      l3   = exp(ll3);
      s2b1 = l1*l1*s2;
      cb12 = l2*l1*s2;
      s2b2 = (l2*l2 + l3*l3)*s2;
      cl   = exp(beta1 + b1);
      ka   = exp(beta2 + b2);
      ke   = exp(beta3);
      pred = dose*ke*ka*(exp(-ke*time)-exp(-ka*time))/cl/(ka-ke);
      model conc   ~ normal(pred,s2);
      random b1 b2 ~ normal([0,0],[s2b1,cb12,s2b2]) subject=subject;
   run;

Top of Page