### SAS Library Nonlinear Regression in SAS

The SAS System offers a powerful procedure to fit nonlinear regression models, PROC NLIN. Since I get many questions in statistical consulting sessions on how to fit a nonlinear regression and how to compare treatments in an experiments with nonlinear response models, I decided to put together some of the essentials.
Nonlinear Regression vs. Linear Regression
2. Fitting Nonlinear Regressions is an Iterative Process
3. PROC NLIN 4. Testing Hypotheses

1. Nonlinear Regression vs. Linear Regression

A regression model is called nonlinear, if the derivatives of the model with respect to the model parameters depends on one or more parameters. This definition is essential to distinguish nonlinear from curvilinear regression. A regression model is not necessarily nonlinear if the graphed regression trend is curved. A polynomial model such as y = b0 + b1x + b2x2 + e appears curved when y is plotted against x. It is, however, not a nonlinear model. To see this, take derivatives of y with respect to the parameters b0, b1, and b2: dy/db0 = 1, dy/db1 = x, dy/db2 = x2 None of these derivatives depends on a model parameter, the model is linear. In contrast, consider the log-logistic model y = d + (a - d)/(1 + exp{b log(x/g)}) + e Take derivatives with respect to d, for example: dy/dd = 1 - 1/(1 + exp{b log(x/g)}). The derivative involves other parameters, hence the model is nonlinear.

Fitting a nonlinear regression model to data is slightly more involved than fitting a linear model, but they have specific advantages:
• Nonlinear models are often derived on the basis of physical and/or biological considerations, e.g., from differential equations, and have justification within a quantitative conceptualization of the process of interest.

• The parameters of a nonlinear model usually have direct interpretation in terms of the process under study. In the case of the log-logistic model above, for example, the response takes on a sigmoidal shape between d and a. g is the value for which the response achieves (a + d)/2.

• Constraints can be built into a nonlinear model easily and are harder to enforce for linear models. If, e.g., the response achieves an asymptotic value as x grows, many nonlinear models have such behavior built in automatically.
• 2.Fitting Nonlinear Regressions is an Iterative Process

One of the disadvantages of nonlinear models is that the process is iterative. To estimate the parameters of the model, you commence with a set of user-supplied starting values. The software then tries to improve on the quality of the model fit to the data by adjusting the values of the parameters successively. The adjustment of all parameters is considered one iteration. In the next iteration, the program again attempts to improve on the fit by modifying the parameters. Once an improvement is not possible, the fit is considered converged. Care must be exercised in choosing good starting values. The fact that the program can not improve on the model fit between successive iterations may not indicate that the best parameter estimates have been found, but indicate lack of progress of the iterative algorithm. It is possible to send the algorithm off into regions of the parameter space, from which it can not escape, but that do not provide the best estimates. It is thus sensible to start the iterative process with different sets of starting values and to observe whether the program arrives at the same parameter estimates. If it does, your fine.
3. PROC NLIN

3.1. A simple example

The SAS procedure to fit nonlinear regression is PROC NLIN. It was improved dramatically in release 6.12 of The SAS System with the addition of a differentiator. Prior to release 6.12, if you wanted to fit a nonlinear model you had to supply the model specification as well as the formulas for the derivatives of the model. The latter was a real hassle, especially if the model is complicated. A method to circumvent the specification of derivatives was to choose a fitting algorithm that approximates the derivatives by differences. This algorithm, known as DUD (Does not Use Derivatives) was hence very popular. However, the algorithm is also known to be quite poor in computing good estimates. A method using derivatives is to be preferred. With release 6.12 SAS will calculate derivatives for you if you wish. The user still has the option to supply derivatives. It is, however, recommended to let The SAS System calculate them for you.

The minimum specification to fit a nonlinear regression with PROC NLIN demands that the user specifies the model and the parameters in the model. All terms in the model not defined as parameters are looked for in the data set that PROC NLIN processes. As an example, consider the following data set:
data weeds;
input tx rate y;
rate = rate * 1.12; /* convert from lbs/acre to kg/ha */
if rate < 1E-6 then rate = 1E-6;
datalines;
1        0.000           99
1        0.020           84
1        0.040           95
1        0.080           84
1        0.160           53
1        0.320            6
1        0.641            6
1        0.000          103
1        0.020           84
1        0.040           94
1        0.080           79
1        0.160           75
1        0.320           27
1        0.641            7
1        0.000          113
1        0.020           91
1        0.040           80
1        0.080           76
1        0.160           52
1        0.320            6
1        0.641            6
1        0.000           86
1        0.020           78
1        0.040           85
1        0.080           80
1        0.160           53
1        0.320           30
1        0.641            8
1        0.000          110
1        0.020          104
1        0.040           89
1        0.080           84
1        0.160           44
1        0.320           17
1        0.641            9
1        0.000           94
1        0.020          103
1        0.040           97
1        0.080           85
1        0.160           58
1        0.320           17
1        0.641            7
1        0.000           95
1        0.020          113
1        0.040           85
1        0.080           79
1        0.160           33
1        0.320           19
1        0.641            4
1        0.000          101
1        0.020          107
1        0.040          105
1        0.080           87
1        0.160           75
1        0.320           20
1        0.641           11
;
run;

It consists of a rate variable in kg/ha and a response variable y. Let's ignore the treatment variable TX for the time being. We wish to fit the log-logistic model y = d + (a - d)/(1 + exp{b log(x/g)}) + e to the data. Here is the minimal code:
proc nlin data=weeds;
parameters  alpha=100 delta=4 beta=2.0 gamma=0.2;
model y   = delta + (alpha-delta)  / (1 + exp(beta*log(rate/gamma)));
run;
The PARAMETERS statement defines which elements of the model statement are parameters and to be estimated. You can give them any valid SAS names (up to eight characters long, not starting with a number). I usually use mnemonics for the Greek letters in the statistical model specification. The values following the parameters are their starting values. The MODEL statement contains the mathematical expression of the model, apart from the error term. The only variable on the left hand side not defined in the PARAMETERS statement is RATE, which is looked for in the data set WEEDS. Here is the default output:
Non-Linear Least Squares Iterative Phase
Dependent Variable Y   Method: Gauss-Newton

Iter       ALPHA           DELTA           BETA            GAMMA    Sum of Squares
0     100.000000        4.000000       2.000000        0.200000    5284.207076
1      97.265416        1.238334       2.202450        0.191326    4468.645787
2      97.155009        0.894100       2.222681        0.193661    4464.406826
3      97.109336        1.039029       2.236293        0.193424    4464.263887
4      97.097903        1.047721       2.238380        0.193463    4464.257070
5      97.095319        1.052604       2.238974        0.193461    4464.256734
6      97.094723        1.053402       2.239095        0.193462    4464.256717
NOTE: Convergence criterion met.

Non-Linear Least Squares Summary Statistics     Dependent Variable Y

Source                 DF Sum of Squares     Mean Square
Regression              4    300475.74328     75118.93582
Residual               52      4464.25672        85.85109
Uncorrected Total      56   304940.00000
(Corrected Total)      55    74538.85714

Parameter      Estimate    Asymptotic              Asymptotic 95 %
Std. Error         Confidence Interval
Lower          Upper
ALPHA       97.09472257  2.2526644734  92.574427306   101.61501784
DELTA        1.05340175  5.4120397147  -9.806634348    11.91343785
BETA         2.23909502  0.3642121689    1.508250919    2.96993912
GAMMA        0.19346184  0.0162266928   0.160900642     0.22602303

Asymptotic Correlation Matrix

Corr               ALPHA              DELTA               BETA              GAMMA
----------------------------------------------------------------------------------
ALPHA                  1       -0.301160511         -0.52588196       -0.108538167
DELTA       -0.301160511                  1        0.7586141693       -0.754366204
BETA        -0.52588196        0.7586141693                   1       -0.465204141
GAMMA       -0.108538167       -0.754366204        -0.465204141                  1
The first table of output refers to the iteration history of the model. The row ITER 0 contains the starting values and the residual sum of squares that was achieved with these. Then SAS started to update these values in successive iteration. Notice how the residual sum of squares dropped until no improvement was realized after the sixth iteration. The model has converged.

The next table contains a little ANOVA table partitioning the total sum of squares for the model and data into a sum of squares explained by the model (SS(Regression)) and a residual sum of squares. The mean square error of this fit is 85.85109, your estimate of variability in the data when adjusted for the nonlinear, log-logistic trend.

The ANOVA table is followed by a table of parameter estimates. For each name listed in the PARAMETERS statement SAS prints the estimate, its standard error, and a asymptotic 95% confidence interval. Notice that all results in nonlinear regression are asymptotic. That means the standard error, for example, is only correct if you have an infinitely large sample size. For any finite sample size, the reported standard error is only an approximation which improves with increasing sample size.

Typical for nonlinear regression are correlations among the parameter estimates. Basically what this says is that if you change the value of a specific parameter, other parameters will also change. This is nothing to particularly worry about (as is often stated). However, if the correlations are very high, the fit of the model may be negatively affected by this. You can then choose another fitting algorithm.
3.2. Using a Starting Value Grid
If you are not sure about the starting values, you can use a grid by offering SAS more than one starting value. It will calculate the initial residual sum of squares for all combinations of starting values and start the iterations with the best set. For example,
proc nlin data=weeds;
parameters  alpha=100 delta=4
beta=1 to 2 by 0.5
gamma=0.1 to 0.4 by 0.1;
model y   = delta + (alpha-delta)  / (1 + exp(beta*log(rate/gamma)));
run;
will create 12 different sets of starting values. The iteration table looks as follows:

Non-Linear Least Squares Grid Search      Dependent Variable Y
ALPHA           DELTA           BETA            GAMMA    Sum of Squares
100.000000        4.000000       1.000000        0.100000   18289.313156
100.000000        4.000000       1.500000        0.100000   15482.810583
100.000000        4.000000       2.000000        0.100000   16378.826233
100.000000        4.000000       1.000000        0.200000   12021.327337
100.000000        4.000000       1.500000        0.200000    6699.920592
100.000000        4.000000       2.000000        0.200000    5284.207076
100.000000        4.000000       1.000000        0.300000   17317.275381
100.000000        4.000000       1.500000        0.300000   14778.328159
100.000000        4.000000       2.000000        0.300000   14874.379385
100.000000        4.000000       1.000000        0.400000   24644.017200
100.000000        4.000000       1.500000        0.400000   25883.425514
100.000000        4.000000       2.000000        0.400000   28375.884864

Non-Linear Least Squares Iterative Phase
Dependent Variable Y   Method: Gauss-Newton
Iter       ALPHA           DELTA           BETA            GAMMA    Sum of Squares
0     100.000000        4.000000       2.000000        0.200000    5284.207076
1      97.265416        1.238334       2.202450        0.191326    4468.645787
2      97.155009        0.894100       2.222681        0.193661    4464.406826
3      97.109336        1.039029       2.236293        0.193424    4464.263887
4      97.097903        1.047721       2.238380        0.193463    4464.257070
5      97.095319        1.052604       2.238974        0.193461    4464.256734
6      97.094723        1.053402       2.239095        0.193462    4464.256717

NOTE: Convergence criterion met.

The sixth set was the best one, it produced the smallest residual sum of squares (5284.2). SAS uses this set to commence the iterations (see values of parameters in iteration 0).
3.3. Using Expressions inside NLIN
In contrast to most SAS procedures, you can calculate variables and expressions inside PROC NLIN. This is especially useful if your model consists of many complicated terms. Here is an example:

proc nlin data=weeds;
parameters  alpha=100 delta=4 beta=2.0 gamma=0.2;
term = 1 + exp(beta*log(rate/gamma));
model y = delta + (alpha-delta)  / term;
run;

The messy denominator of the log-logistic model is stored in variable TERM which is used in the MODEL specification.

3.4.  Fixing a parameter
Sometimes you want to fix a parameter at a certain value, rather than estimating it. This is useful if you do not have sufficient data to estimate all parameters precisely and/or you know the value a parameter should take on. For the particular data set, theory tells us that a, the upper asymptote of the response, should be 100, since the response is expressed relative to a control value from which it should monotonically decrease. To fix a at 100, simply remove it from the PARAMETERS statement and give alpha the intended value:
proc nlin data=weeds;
parameters  delta=4 beta=2.0 gamma=0.2;
alpha = 100;
term = 1 + exp(beta*log(rate/gamma));
model y = delta + (alpha-delta)  / term;
run;

3.5.  Calculating a R2-Type Measure
Users of linear regression models are accustomed to expressing the quality of fit of a model in terms of the coefficient of determination, also known as R2. In nonlinear regression, such a measure is unfortunately, not readily defined. One of the problems with the R2 definition is that it requires the presence of an intercept, which most nonlinear models do not have. A measure, relatively closely corresponding to R2 in the nonlinear case is Pseudo-R2 = 1 - SS(Residual)/SS(TotalCorrected).
3.6.  Choosing the Fitting Algorithm
If your data and model are well behaved, it should not make a difference how you fit the nonlinear model to data. Unfortunately, this can not be said for all nonlinear regression models. You may have to choose carefully, which algorithm to use. In PROC NLIN different fitting algorithms are invoked with the METHOD= option of the PROC NLIN statement. Here are a few guidelines:
• If possible, choose a method that uses derivatives, avoid DUD
• If the parameters are highly correlated, choose the Levenberg-Marquardt method (keyword METHOD=MARQUARDT)
• Among the derivative dependent methods, prefer the Newton-Raphson (METHOD=NEWTON) over the Gauss (METHOD=GAUSS) method.

Unfortunately, if you do not specify derivatives and a METHOD= option, SAS will default to the DUD method. I therefore always use either of the following two PROC NLIN statements:

proc nlin data=whatever method=newton;
<statements>
run;
or if the parameter estimates show high correlations (> 0.8, < -0.8)
proc nlin data=whatever method=marquardt;
<statements>
run;

3.7.  Calculating predicted values and their confidence intervals
Predicted values are not displayed on screen in PROC NLIN. However, you can request to save them to a data set for later use. Along with the predicted values, you can calculate confidence bounds for the mean predictions, prediction intervals for an individual predictions and so forth. In the log-logistic example, here is the code:
proc nlin data=weeds method=newton;
parameters  alpha=100 delta=4 beta=2.0 gamma=0.2;
term = 1 + exp(beta*log(rate/gamma));
model y   = delta + (alpha-delta)  / term;
output out=nlinout predicted=pred l95m=l95mean u95m=u95mean
l95=l95ind u95=u95ind;
run;
proc print data=nlinout; run;
This will produce a data set NLINOUT containing predicted values (in variable PRED), upper and lower 95% confidence bounds for the mean prediction (variables U95MEAN and L95MEAN), and upper and lower 95% prediction intervals for an individual observation (variables U95IND and L95IND). All variables from the data set being processed are also copied into the output data set. Not repeating the standard PROC NLIN output, here is a printout of the data set NLINOUT:
OBS    TX      RATE        Y      PRED       L95MEAN    U95MEAN       L95IND     U95IND

1     1    0.00000      99    97.0946     92.3001     101.889     77.8936    116.296
2     1    0.02240      84    96.3318     92.2998     100.364     77.3069    115.357
3     1    0.04480      95    93.5968     90.1265      97.067     74.6830    112.511
4     1    0.08960      84    82.5522     76.8186      88.286     63.0954    102.009
5     1    0.17920      53    53.1811     47.8220      58.540     33.8314      72.531
6     1    0.35840       6    20.3498     14.8513      25.848      0.9611      39.739
7     1    0.71792       6     5.8938      -0.2713     12.059    -13.6944      25.482
8     1    0.00000    103     97.0946     92.3001    101.889      77.8936    116.296
9     1    0.02240      84    96.3318     92.2998     100.364     77.3069    115.357
10     1    0.04480      94    93.5968     90.1265      97.067     74.6830    112.511
11     1    0.08960      79    82.5522     76.8186      88.286     63.0954    102.009
12     1    0.17920      75    53.1811     47.8220      58.540     33.8314      72.531
13     1    0.35840      27    20.3498     14.8513      25.848      0.9611      39.739
14     1    0.71792       7     5.8938      -0.2713     12.059    -13.6944      25.482
15     1    0.00000    113     97.0946     92.3001    101.889      77.8936    116.296
16     1    0.02240      91    96.3318     92.2998     100.364     77.3069    115.357
17     1    0.04480      80    93.5968     90.1265      97.067     74.6830    112.511

and so forth.

3.8.  Calculating predicted values for plotting purposes
If you want to plot the predicted regression trend in publication quality you need a fairly dense set of x values for which to calculate the predictions, otherwise your graphed trend will look kinky. There are two ways to do this:

1. Take the values of the parameter estimates from SAS into a graphics program and calculate the predicted values for a set of dense x values there.
2. Fool SAS by appending a second data set that contains the values for which you want predictions.

I will discuss how to do the latter. Assume in the example, we want predicted values for rates ranging from 0 to 0.8 g/ha in 0.05 intervals. We want SAS to calculate predicted values, but of course we can not use these data points in fitting the nonlinear model. We can not pretend to have more data than there are. Set up a second data set containing the rates for which you want predictions, but no response:

data filler;
do rate = 0.05 to 0.8 by 0.05;
predict=1;
output;
end;
run;
Merge this with the original data:

data fitthis;
set weeds filler;
run;

and run PROC NLIN on this new data set:
proc nlin data=fitthis method=newton;
parameters  alpha=100 delta=4 beta=2.0 gamma=0.2;
term = 1 + exp(beta*log(rate/gamma));
model y   = delta + (alpha-delta)  / term;
output out=nlinout predicted=pred;
run;
proc print data=nlinout(where=(predict=1));
run;
SAS will exclude the observations coming from data set FILLER in fitting the model, since there is no information about the response Y in this data set. When calculating predicted values, it is only looking for variables on the right hand side of the model statement and will calculate predictions for the observations from FILLER. By including the variable predict in FILLER with value 1 which is not in data set WEEDS which contains the original data, we can pull out the predicted values in the FILLER data. Here are the observations of the output data set:

OBS    TX    RATE    Y     PREDICT      PRED

57     .    0.05    .        1       92.6668
58     .    0.10    .        1       79.2515
59     .    0.15    .        1       62.3952
60     .    0.20    .        1       47.2881
61     .    0.25    .        1       35.6570
62     .    0.30    .        1       27.2184
63     .    0.35    .        1       21.1817
64     .    0.40    .        1       16.8346
65     .    0.45    .        1       13.6563
66     .    0.50    .        1       11.2901
67     .    0.55    .        1        9.4958
68     .    0.60    .        1        8.1112
69     .    0.65    .        1        7.0252
70     .    0.70    .        1        6.1607
71     .    0.75    .        1        5.4632
72     .    0.80    .        1        4.8936

By changing the step size in the DO loop when creating data set FILLER, you can make the grid of predicted values arbitrarily fine.
4. Testing Hypotheses

4.1.  Testing hypotheses about a single parameter
The default NLIN output includes asymptotic 95% confidence intervals for every parameter in the model. These can be used to test hypotheses about the parameter. To perform a two-sided test at the 5% level, simply check if the hypothesized value is inside or outside of the confidence bounds. If it is outside, reject the null hypothesis. If it is inside the interval you fail to reject.
Parameter   Estimate    Asymptotic              Asymptotic 95 %
Std. Error         Confidence Interval
Lower          Upper
ALPHA      97.09472257  2.2526644734  92.574427306   101.61501784
DELTA       1.05340175  5.4120397147   -9.806634348   11.91343785
BETA        2.23909502  0.3642121689    1.508250919    2.96993912
GAMMA       0.19346184  0.0162266928    0.160900642    0.22602303
In this example, to test H0: b = 2.5, compare the value 2.5 against the confidence bounds for b. Since it is between lower and upper bound, do not reject H0. You would have rejected H0: b = 3.0 or H0: b = 1.2, for example.
To test one-sided hypotheses or test at significance levels other than 5%, you can use the parameter estimate itself and its standard error to construct tests. To test H0: d = c at significance level n against H1: d > c, for example, compare the test statistic tobs = (d - c)/se(d) against the cutoff from a t distribution with degrees of freedom equal to the residual degrees of freedom in the ANOVA table produced by NLIN: tn,df(Res). d denotes the estimate of the parameter d.

4.2. Comparing treatments

To compare two or more treatments in a nonlinear regression problem, we proceed similarly as in the case of a standard analysis of variance. First, we assess whether there are any differences between the treatments. If there are, we try to find out where the differences occur. To answer the first question, whether there are any treatment differences, a sum of square reduction test is used. This is a very general procedure that can be used to test all kinds of hypotheses. The idea is to fit two versions of the model. One is considered the full model and has more parameters. The reduced model with fewer parameters is a constrained version of the full model. You arrive at the reduced model by constraining some of the parameters of the full model. The constraint itself is implied by the hypotheses you test.

Assume that there are two different treatments for which the log-logistic model is to be fit. The hypothesis to be tested is H0: There are no differences in log-logistic response among the two treatments The full model is yj = dj + (aj - dj)/(1 + exp{bj log(x/gj)}) + e where subscript j identifies the treatment. This model has 2*4 = 8 parameters (two sets of four parameters, one set for each treatment). To derive the reduced model, we must invoke the hypothesis H0. If there are no differences in log-logistic response, the two treatments will share the same d parameter, the same a parameter, and so forth. Hence, the reduced model is yj = d + (a - d)/(1 + exp{b log(x/g)}) + e. To perform the sum of squares reduction test, you need to fit both the full and reduced model. Then, calculate the test statistic Fobs = (SS(Residual)Reduced - SS(Residual)Full)/(df(Residual)Reduced - df(Residual)Full) }/MSError(Full) and compare to the cutoffs from a F distribution with df(Residual)Reduced - df(Residual)Full numerator and df(Residual)Full denominator degrees of freedom. In the log-logistic example, assume that there are two treatments. The variable TX with possible values 1 and 2 identifies the treatment an observation is associated with. To fit the full model we need to fit a eight parameter model.
data weeds2;
input tx rate y;
rate = rate * 1.12; /* convert from lbs/acre to kg/ha */
if rate < 1E-6 then rate = 1E-6;
datalines;
1          0         99
1        .02         84
1        .04         95
1        .08         84
1        .16         53
1        .32          6
1       .641          6
1          0        103
1        .02         84
1        .04         94
1        .08         79
1        .16         75
1        .32         27
1       .641          7
1          0        113
1        .02         91
1        .04         80
1        .08         76
1        .16         52
1        .32          6
1       .641          6
1          0         86
1        .02         78
1        .04         85
1        .08         80
1        .16         53
1        .32         30
1       .641          8
1          0        110
1        .02        104
1        .04         89
1        .08         84
1        .16         44
1        .32         17
1       .641          9
1          0         94
1        .02        103
1        .04         97
1        .08         85
1        .16         58
1        .32         17
1       .641          7
1          0         95
1        .02        113
1        .04         85
1        .08         79
1        .16         33
1        .32         19
1       .641          4
1          0        101
1        .02        107
1        .04        105
1        .08         87
1        .16         75
1        .32         20
1       .641         11
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
2          0         99
2        .02         92
2        .04         82
2        .08         63
2        .16         40
2        .32         22
2       .641         12
;
run;

proc nlin data=weeds2 method=marquardt;
parameters alpha1=100 delta1=4 beta1=2.0 gamma1=0.2
alpha2=100 delta2=4 beta2=2.0 gamma2=0.2;
term1  = 1 + exp(beta1*log(rate/gamma1));
term2  = 1 + exp(beta2*log(rate/gamma2));
model y   = (delta1 + (alpha1-delta1 ) / term1 ) * (Tx = 1) +
(delta2 + (alpha2-delta2 ) / term2 ) * (Tx = 2) ;
run;
Notice how the MODEL statement is written. The expression (Tx =1) is a logical expression which returns the value 1 if true, 0 otherwise. For observations from treatment 1, the model statement becomes

model y   = (delta1 + (alpha1-delta1 ) / term1 );

for observations from treatment 2 it becomes

model y   = (delta2 + (alpha2-delta2 ) / term2 );

This technique can easily be extended to more than two treatments and more than one treatment factor. Here is the ANOVA table of this PROC NLIN run:

Sum of        Mean               Approx
Source                    DF     Squares      Square    F Value    Pr > F

Regression                 8      549963     68745.4     426.95    <.0001
Residual                 104      4465.2     42.9342
Uncorrected Total        112      554428

Corrected Total          111      132782

The reduced model can be fit in two ways. Either write a model statement with only four parameters:
proc nlin data=weeds method=marquardt;
parameters alpha=100 delta=4 beta=2.0 gamma=0.2;
term  = 1 + exp(beta*log(rate/gamma));
model y   = delta + (alpha-delta)  / term;
run;

or use the model statement from the full model and impose the necessary constraints:

proc nlin data=weeds method=marquardt;
parameters alpha1=100 delta1=4 beta1=2.0 gamma1=0.2;
alpha2 = alpha1;
beta2  = beta1;
delta2 = delta1;
gamma2 = gamma1;
term1  = 1 + exp(beta1*log(rate/gamma1));
term2  = 1 + exp(beta2*log(rate/gamma2));
model y   = (delta1 + (alpha1-delta1 ) / term1 ) * (Tx = 1) +
(delta2 + (alpha2-delta2 ) / term2 ) * (Tx = 2) ;
run;

Either way, the ANOVA table for the reduced model is

Sum of        Mean               Approx
Source                    DF     Squares      Square    F Value    Pr > F

Regression                 4      547094      136773     615.74    <.0001
Residual                 108      7334.5     67.9117
Uncorrected Total        112      554428

Corrected Total          111      132782

To test whether there is any difference among the two treatments, calculate the sum of square reduction test: Fobs = { (7334.5 - 4465.2 ) / (108 - 104 ) } / 42.9342  = 16.70754 and the p-value under a F4,104 distribution is 1.32E-10, H0 is rejected.
To calculate the p-value use the PROBF function in SAS:

data one;
prob = 1 - probf(16.70754,4,104);
run;

proc print;
run;
4.3.  Comparing single parameters among treatments
If you have multiple treatments or experimental conditions, you may want to compare a given parameter among the two. For example, are the g parameters between the treatments significantly different. Since g measures the rate at which the response achieves its half-way point, this is a very important hypothesis in dose-response investigations where g represents a lethal dosage that kills 50% of the subjects, or a dosage that reduces (enhances) growth by 50%, etc. To test this hypothesis, you could again perform a sum of square reduction test comparing the full model yj = dj + (aj - dj)/(1 + exp{bj log(x/gj)}) + e and the reduced model yj = dj + (aj - dj)/(1 + exp{bj log(x/g)}) + e. There is a simpler technique which answers the hypothesis in a single run of the model and provides an estimate of the actual difference between the two treatments. Parameterize the model as yj = dj + (aj - dj)/(1 + exp{bj log(x/ [g+D(Tx=2)])}) + e. g is the parameter for treatment 1 and the term D(Tx=2) does the trick. For treatment 1, (Tx=2) will be false and the term D(Tx=2) is zero. The model becomes yj = dj + (aj - dj)/(1 + exp{bj log(x/ g)}) + e. For treatment 2 (Tx=2) is true and the model becomes yj = dj + (aj - dj)/(1 + exp{bj log(x/ [g+D])}) + e. Consequently, D measures the difference in g parameters between the two treatments. Since this difference has become a parameter of the model you get an estimate for the difference, its standard error and confidence interval os part of the default NLIN output. Check whether the confidence interval contains 0. If it does not, reject the hypothesis H0: g1 = g2. Here is the code (gammaD denotes D):

proc nlin data=weeds method=marquardt;
parameters alpha1=100 delta1=4 beta1=2.0 gamma1 =0.2
gamma  = gamma1 + gammaD*(Tx = 2);
term1  = 1 + exp(beta1*log(rate/gamma));
term2  = 1 + exp(beta2*log(rate/gamma));
model y = (delta1 + (alpha1-delta1 ) / term1 ) * (Tx = 1) +
(delta2 + (alpha2-delta2 ) / term2 ) * (Tx = 2) ;
run;

and the output:

Source                    DF     Squares      Square    F Value    Pr > F

Regression                 8      549963     68745.4     426.95    <.0001
Residual                 104      4465.2     42.9342
Uncorrected Total        112      554428

Corrected Total          111      132782

Approx
Parameter      Estimate    Std Error    Approximate 95% Confidence Limits

alpha1          97.0946       1.5930     93.9355       100.3
delta1           1.0536       3.8272     -6.5360      8.6432
beta1            2.2391       0.2576      1.7284      2.7499
gamma1           0.1935       0.0115      0.1707      0.2162
alpha2          99.0658       2.1380     94.8260       103.3
delta2           5.4066       4.3049     -3.1303     13.9434
beta2            1.4607       0.1656      1.1322      1.7892

The estimate of D, denoted GAMMAD on the output is -0.0693. The confidence interval does not contain 0, hence the difference in g parameters between the two treatments is significant at the 5% level. Notice that an obvious starting value for GAMMAD is zero, implying there is no difference.

This technique can used to fit the full model by the way. Simply express the parameter for one treatment in terms of the parameters of the other treatment plus some parameter specific difference:

proc nlin data=weeds method=marquardt;
parameters alpha1=100 delta1=4.0 beta1=2.0 gamma1 =0.2
alpha  = alpha1 + alphaD * (Tx = 2);
beta   = beta1  + betaD  * (Tx = 2);
gamma  = gamma1 + gammaD * (Tx = 2);
delta  = delta1 + deltaD * (Tx = 2);

term  = 1 + exp(beta*log(rate/gamma));
model y = (delta + (alpha-delta ) / term );
run;

with output

Source                    DF     Squares      Square    F Value    Pr > F

Regression                 8      549963     68745.4     426.95    <.0001
Residual                 104      4465.2     42.9342
Uncorrected Total        112      554428

Corrected Total          111      132782

Approx
Parameter      Estimate    Std Error    Approximate 95% Confidence Limits

alpha1          97.0946       1.5930     93.9355       100.3
delta1           1.0536       3.8272     -6.5360      8.6432
beta1            2.2391       0.2576      1.7284      2.7499
gamma1           0.1935       0.0115      0.1707      0.2162

Notice that the ANOVA table is exactly the same as that of the full model in 4.2. This is the full model.

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.