Getting Started
Nonlinear Growth Curves with Gaussian Data
As an introductory example, consider the orange tree data of Draper and Smith (1981). These data consist of seven measurements of the trunk circumference (in millimeters) on each of five orange trees. You can input these data into a SAS data set as follows :
data tree; input tree day y; datalines; 1 118 30 1 484 58 1 664 87 1 1004 115 1 1231 120 1 1372 142 1 1582 145 2 118 33 2 484 69 2 664 111 2 1004 156 2 1231 172 2 1372 203 2 1582 203 3 118 30 3 484 51 3 664 75 3 1004 108 3 1231 115 3 1372 139 3 1582 140 4 118 32 4 484 62 4 664 112 4 1004 167 4 1231 179 4 1372 209 4 1582 214 5 118 30 5 484 49 5 664 81 5 1004 125 5 1231 142 5 1372 174 5 1582 177 ;
Lindstrom and Bates (1990) and Pinheiro and Bates (1995) propose the following logistic nonlinear mixed model for these data:
Here, y ij represents the j th measurement on the i th tree ( i = 1 , , 5; j = 1 , , 7), d ij is the corresponding day, b 1 ,b 2 ,b 3 are the fixed-effects parameters, u i 1 are the random-effect parameters assumed to be iid 
, and e ij are the residual errors assumed to be iid 
and independent of the u i 1 . This model has a logistic form, and the random-effect parameters u i 1 enter the model linearly.
The statements to fit this nonlinear mixed model are as follows:
proc nlmixed data=tree; parms b1=190 b2=700 b3=350 s2u=1000 s2e=60; num = b1+u1; ex = exp(-(day-b2)/b3); den = 1 + ex; model y ~ normal(num/den,s2e); random u1 ~ normal(0,s2u) subject=tree; run;
The PROC NLMIXED statement invokes the procedure and inputs the TREE data set. The PARMS statement identifies the unknown parameters and their starting values. Here there are three fixed-effects parameters (B1, B2, B3) and two variance components (S2U, S2E).
The next three statements are SAS programming statements specifying the logistic mixed model. A new variable U1 is included to identify the random effect. These statements are evaluated for every observation in the data set when PROC NLMIXED computes the log likelihood function and its derivatives.
The MODEL statement defines the dependent variable and its conditional distribution given the random effects. Here a normal (Gaussian) conditional distribution is specified with mean NUM/DEN and variance S2E.
The RANDOM statement defines the single random effect to be U1, and specifies that it follows a normal distribution with mean 0 and variance S2U. The SUBJECT= argument defines a variable indicating when the random effect obtains new realizations; in this case, it changes according to the values of the TREE variable. PROC NLMIXED assumes that the input data set is clustered according to the levels of the TREE variable; that is, all observations from the same tree occur sequentially in the input data set.
The output from this analysis is as follows.
The NLMIXED Procedure Specifications Data Set WORK.TREE Dependent Variable y Distribution for Dependent Variable Normal Random Effects u1 Distribution for Random Effects Normal Subject Variable tree Optimization Technique Dual Quasi-Newton Integration Method Adaptive Gaussian Quadrature
The Specifications table lists some basic information about the nonlinear mixed model you have specified. Included are the input data set, dependent and subject variables , random effects, relevant distributions, and type of optimization.
The NLMIXED Procedure Dimensions Observations Used 35 Observations Not Used 0 Total Observations 35 Subjects 5 Max Obs Per Subject 7 Parameters 5 Quadrature Points 1
The Dimensions table lists various counts related to the model, including the number of observations, subjects, and parameters. These quantities are useful for checking that you have specified your data set and model correctly. Also listed is the number of quadrature points that PROC NLMIXED has selected based on the evaluation of the log likelihood at the starting values of the parameters. Here, only one quadrature point is necessary because the random-effect parameters u i 1 enter the model linearly.
The NLMIXED Procedure Parameters b1 b2 b3 s2u s2e NegLogLike 190 700 350 1000 60 132.491787
The Parameters table lists the parameters to be estimated, their starting values, and the negative log likelihood evaluated at the starting values.
The NLMIXED Procedure Iteration History Iter Calls NegLogLike Diff MaxGrad Slope 1 4 131.686742 0.805045 0.010269 -0.633 2 6 131.64466 0.042082 0.014783 -0.0182 3 8 131.614077 0.030583 0.009809 -0.02796 4 10 131.572522 0.041555 0.001186 -0.01344 5 11 131.571895 0.000627 0.0002 -0.00121 6 13 131.571889 5.549E-6 0.000092 -7.68E-6 7 15 131.571888 1.096E-6 6.097E-6 -1.29E-6 NOTE: GCONV convergence criterion satisfied.
The Iterations table records the history of the minimization of the negative log likelihood. For each iteration of the quasi-Newton optimization, values are listed for the number of function calls, the value of the negative log likelihood, the difference from the previous iteration, the absolute value of the largest gradient, and the slope of the search direction. The note at the bottom of the table indicates that the algorithm has converged successfully according to the GCONV convergence criterion, a standard criterion computed using a quadratic form in the gradient and inverse Hessian.
The NLMIXED Procedure Fit Statistics -2 Log Likelihood 263.1 AIC (smaller is better) 273.1 AICC (smaller is better) 275.2 BIC (smaller is better) 271.2
The Fitting Information table lists the final maximized value of the log likelihood as well as the information criteria of Akaike and Schwarz in two different forms. These statistics can be used to compare different nonlinear mixed models.
The NLMIXED Procedure Parameter Estimates Standard Parameter Estimate Error DF t Value Pr > t Alpha Lower b1 192.05 15.6473 4 12.27 0.0003 0.05 148.61 b2 727.90 35.2472 4 20.65 <.0001 0.05 630.04 b3 348.07 27.0790 4 12.85 0.0002 0.05 272.88 s2u 999.88 647.44 4 1.54 0.1974 0.05 -797.70 s2e 61.5139 15.8831 4 3.87 0.0179 0.05 17.4153 Parameter Estimates Parameter Upper Gradient b1 235.50 1.154E-6 b2 825.76 5.289E-6 b3 423.25 -6.1E-6 s2u 2797.45 -3.84E-6 s2e 105.61 2.892E-6
The Parameter Estimates table lists the maximum likelihood estimates of the five parameters and their approximate standard errors computed using the final Hessian matrix. Approximate t -values and Wald-type confidence limits are also provided, with degrees of freedom equal to the number of subjects minus the number of random effects. You should interpret these statistics cautiously for variance parameters like S2U and S2E. The final column in the output is the gradient vector at the optimization solution. Each element appears to be sufficiently small to indicate a stationary point.
Since the random-effect parameters u i 1 enter the model linearly, you can obtain equivalent results by using the first-order method (specify METHOD=FIRO in the PROC NLMIXED statement).
Logistic-Normal Model with Binomial Data
This example analyzes the data from Beitler and Landis (1985), which represent results from a multi-center clinical trial investigating the effectiveness of two topical cream treatments (active drug, control) in curing an infection. For each of eight clin-ics, the number of trials and favorable cures are recorded for each treatment. The SAS data set is as follows.
data infection; input clinictxn; datalines; 1 1 11 36 1 0 10 37 2 1 16 20 2 0 22 32 3 1 14 19 3 0 7 19 4 1 2 16 4 0 1 17 5 1 6 17 5 0 0 12 6 1 1 11 6 0 0 10 7 1 1 5 7 0 1 9 8 1 4 6 8 0 6 7 run;
Suppose n ij denotes the number of trials for the i th clinic and the j th treatment ( i = 1 ,..., 8 j = 0 , 1), and x ij denotes the corresponding number of favorable cures. Then a reasonable model for the preceding data is the following logistic model with random effects:
and
The notation t j indicates the j th treatment, and the u i are assumed to be iid 
.
The PROC NLMIXED statements to fit this model are as follows:
proc nlmixed data=infection; parms beta0=-1 beta1=1 s2u=2; eta = beta0 + beta1*t + u; expeta = exp(eta); p = expeta/(1+expeta); model x ~ binomial(n,p); random u ~ normal(0,s2u) subject=clinic; predict eta out=eta; estimate 1/beta1 1/beta1; run;
The PROC NLMIXED statement invokes the procedure, and the PARMS statement defines the parameters and their starting values. The next three statements define p ij , and the MODEL statement defines the conditional distribution of x ij to be binomial. The RANDOM statement defines U to be the random effect with subjects defined by the CLINIC variable.
The PREDICT statement constructs predictions for each observation in the input data set. For this example, predictions of · ij and approximate standard errors of prediction are output to a SAS data set named ETA. These predictions include empirical Bayes estimates of the random effects u i .
The ESTIMATE statement requests an estimate of the reciprocal of ² 1 .
The output for this model is as follows.
The NLMIXED Procedure Specifications Data Set WORK.INFECTION Dependent Variable x Distribution for Dependent Variable Binomial Random Effects u Distribution for Random Effects Normal Subject Variable clinic Optimization Technique Dual Quasi-Newton Integration Method Adaptive Gaussian Quadrature
The Specifications table provides basic information about the nonlinear mixed model.
The NLMIXED Procedure Dimensions Observations Used 16 Observations Not Used 0 Total Observations 16 Subjects 8 Max Obs Per Subject 2 Parameters 3 Quadrature Points 5
The Dimensions table provides counts of various variables. You should check this table to make sure the data set and model have been entered properly. PROC NLMIXED selects five quadrature points to achieve the default accuracy in the likelihood calculations.
The NLMIXED Procedure Parameters beta0 beta1 s2u NegLogLike -1 1 2 37.5945925
The Parameters table lists the starting point of the optimization.
The NLMIXED Procedure Iteration History Iter Calls NegLogLike Diff MaxGrad Slope 1 2 37.3622692 0.232323 2.882077 -19.3762 2 3 37.1460375 0.216232 0.921926 -0.82852 3 5 37.0300936 0.115944 0.315897 -0.59175 4 6 37.0223017 0.007792 0.01906 -0.01615 5 7 37.0222472 0.000054 0.001743 -0.00011 6 9 37.0222466 6.57E-7 0.000091 -1.28E-6 7 11 37.0222466 5.38E-10 2.078E-6 -1.1E-9 NOTE: GCONV convergence criterion satisfied.
The Iterations table indicates successful convergence in seven iterations.
The NLMIXED Procedure Fit Statistics -2 Log Likelihood 74.0 AIC (smaller is better) 80.0 AICC (smaller is better) 82.0 BIC (smaller is better) 80.3
The Fitting Information table lists some useful statistics based on the maximized value of the log likelihood.
The NLMIXED Procedure Parameter Estimates Standard Parameter Estimate Error DF t Value Pr > t Alpha Lower beta0 -1.1974 0.5561 7 -2.15 0.0683 0.05 -2.5123 beta1 0.7385 0.3004 7 2.46 0.0436 0.05 0.02806 s2u 1.9591 1.1903 7 1.65 0.1438 0.05 -0.8554 Parameter Estimates Parameter Upper Gradient beta0 0.1175 -3.1E-7 beta1 1.4488 -2.08E-6 s2u 4.7736 -2.48E-7
The Parameter Estimates table indicates marginal significance of the two fixed-effects parameters. The positive value of the estimate of ² 1 indicates that the treatment significantly increases the chance of a favorable cure.
The NLMIXED Procedure Additional Estimates Standard Label Estimate Error DF t Value Pr > t Alpha Lower Upper 1/beta1 1.3542 0.5509 7 2.46 0.0436 0.05 0.05146 2.6569
The Additional Estimates table displays results from the ESTIMATE statement. The estimate of 1 / ² 1 equals 1 / . 7385 = 1 . 3541 and its standard error equals 0 . 3004 / . 7385 2 = 0 . 5509 by the delta method (Billingsley 1986, Cox 1998). Note this particular approximation produces a t -statistic identical to that for the estimate of ² 1 .
Not shown is the ETA data set, which contains the original 16 observations and predictions of the · ij .
