Example 60.5 Comparing Nonlinear Trends among Groups

When you model nonlinear trends in the presence of group (classification) variables, two questions often arise: whether the trends should be varied by group, and how to decide which parameters should be varied across groups. A large battery of tools is available on linear statistical models to test hypotheses involving the model parameters, especially to test linear hypotheses. To test similar hypotheses in nonlinear models, you can draw on analogous tools. Especially important in this regard are comparisons of nested models by contrasting their residual sums of squares.

In this example, a two-group model from a pharmacokinetic application is fit to data that are in part based on the theophylline data from Pinheiro and Bates (1995) and the first example in the documentation for the NLMIXED procedure. In a pharmacokinetic application you study how a drug is dispersed through a living organism. The following data represent concentrations of the drug theophylline over a 25-hour period following oral administration. The data are derived by collapsing and averaging the subject-specific data from Pinheiro and Bates (1995) in a particular, yet unimportant, way. The purpose of arranging the data in this way is purely to demonstrate the methodology.

   data theop;
     input time dose conc @@;
     if (dose = 4) then group=1; else group=2;
     datalines;
    0.00   4  0.1633  0.25   4   2.045
    0.27   4     4.4  0.30   4    7.37
    0.35   4    1.89  0.37   4    2.89
    0.50   4    3.96  0.57   4    6.57
    0.58   4     6.9  0.60   4     4.6
    0.63   4    9.03  0.77   4    5.22
    1.00   4    7.82  1.02   4   7.305
    1.05   4    7.14  1.07   4     8.6
    1.12   4    10.5  2.00   4    9.72
    2.02   4    7.93  2.05   4    7.83
    2.13   4    8.38  3.50   4    7.54
    3.52   4    9.75  3.53   4    5.66
    3.55   4   10.21  3.62   4     7.5
    3.82   4    8.58  5.02   4   6.275
    5.05   4    9.18  5.07   4    8.57
    5.08   4     6.2  5.10   4    8.36
    7.02   4    5.78  7.03   4    7.47
    7.07   4   5.945  7.08   4    8.02
    7.17   4    4.24  8.80   4    4.11
    9.00   4     4.9  9.02   4    5.33
    9.03   4    6.11  9.05   4    6.89
    9.38   4    7.14 11.60   4    3.16
   11.98   4    4.19 12.05   4    4.57
   12.10   4    5.68 12.12   4    5.94
   12.15   4     3.7 23.70   4    2.42
   24.15   4    1.17 24.17   4    1.05
   24.37   4    3.28 24.43   4    1.12
   24.65   4    1.15  0.00   5   0.025
    0.25   5    2.92  0.27   5   1.505
    0.30   5    2.02  0.50   5   4.795
    0.52   5    5.53  0.58   5    3.08
    0.98   5   7.655  1.00   5   9.855
    1.02   5    5.02  1.15   5    6.44
    1.92   5    8.33  1.98   5    6.81
    2.02   5  7.8233  2.03   5    6.32
    3.48   5    7.09  3.50   5   7.795
    3.53   5    6.59  3.57   5    5.53
    3.60   5    5.87  5.00   5     5.8
    5.02   5  6.2867  5.05   5    5.88
    6.98   5    5.25  7.00   5    4.02
    7.02   5    7.09  7.03   5   4.925
    7.15   5    4.73  9.00   5    4.47
    9.03   5    3.62  9.07   5    4.57
    9.10   5     5.9  9.22   5    3.46
   12.00   5    3.69 12.05   5    3.53
   12.10   5    2.89 12.12   5    2.69
   23.85   5    0.92 24.08   5    0.86
   24.12   5    1.25 24.22   5    1.15
   24.30   5     0.9 24.35   5    1.57
   ;

The following code plots the theophylline concentration data over time for the two groups (Output 60.5.1). In each group the concentration tends to rise sharply right after the drug is administered, followed by a prolonged tapering of the concentration.

   proc sgplot data=theop;
      scatter x=time y=conc / group=group;
      yaxis label='Concentration';
      xaxis label='Time';
   run;

Output 60.5.1 Observed Responses in Two Groups

In the context of nonlinear mixed models, Pinheiro and Bates (1995) consider a first-order compartment model for these data. In terms of two fixed treatment groups, the model can be written as

where is the observed concentration in group at time , is the dose of theophylline, is the elimination rate in group , is the absorption rate in group , is the clearance in group , and denotes the model error. Because the rates and the clearance must be positive, you can parameterize the model in terms of log rates and the log clearance:

In this parameterization the model contains six parameters, and the rates and clearance vary by group. This produces two separate response profiles, one for each group. On the other extreme, you could model the trends as if there were no differences among the groups:

In between these two extremes lie other models, such as a model where both groups have the same absorption and elimination rate, but different clearances. The question then becomes how to go about building a model in an organized manner.

To test hypotheses about nested nonlinear models, you can apply the idea of a "Sum of Squares Reduction Test." A reduced model is nested within a full model if you can impose constraints on the full model to obtain the reduced model. Then, if and denote the residual sum of squares in the reduced and the full model, respectively, the test statistic is

where are the number of observations used and are the number of parameters in the full model. The numerator of the statistic is the average reduction in residual sum of squares per constraint. The mean squared error of the full model is used to scale this average because it is less likely to be a biased estimator of the residual variance than the variance estimate from a constrained (reduced) model. The statistic is then compared against quantiles from an distribution with numerator and denominator degrees of freedom. Schabenberger and Pierce (2002) discuss the justification for this test and compare it to other tests in nonlinear models.

In the present application we might phrase the initial question akin to the overall F test for a factor in a linear model: Should any parameters be varied between the two groups? The corresponding null hypothesis is

where the first subscript identifies the type of the parameter and the second subscript identifies the group. Note that this hypothesis implies

If you fail to reject this hypothesis, there is no need to further examine individual parameter differences.

The reduced model—the model subject to the null hypothesis—is fit with the following PROC NLIN statements:

   ods select ConvergenceStatus Anova ParameterEstimates;
   proc nlin data=theop; 
      parms beta1=-3.22 beta2=0.47 beta3=-2.45;
      cl   = exp(beta1); 
      ka   = exp(beta2); 
      ke   = exp(beta3); 
      mean = dose*ke*ka*(exp(-ke*time)-exp(-ka*time))/cl/(ka-ke); 
      model conc = mean;
      ods output Anova=aovred(rename=(ss=ssred ms=msred df=dfred));
   run;

The clearance, the rates, and the mean function are formed independently of the group membership. The analysis of variance table is saved to the data set aovred, and some of its variables are renamed. This is done so that the data set can be merged easily with the analysis of variance table for the full model (see following).

The converged model has a residual sum of square of and a mean squared error of (Output 60.5.2). The table of parameter estimates gives the values for the estimated log clearance , the estimated log absorption rate , and the estimated log elimination rate .