Nonlinear Regression in SAS | SAS Library

This page was adapted from a page created by Oliver Schabenberger.
We thank Professor Schabenberger for
permission to adapt and distribute this page via our web site.

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.

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 = b₀ + b₁x
+ b₂x² + 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 b₀, b₁, and b₂: dy/db₀ = 1, dy/db₁ = x,
dy/db₂ = x² 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 R²-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 R

²

• .

• In nonlinear regression, such a measure is unfortunately, not readily defined. One of the

• problems with the R

²

• definition is that it requires the presence of an

• intercept, which most nonlinear models do not have. A measure, relatively closely corresponding to R

²

• in the nonlinear case is Pseudo-R

²

• = 1 –

• SS(Residual)/SS(Total

_Corrected

).

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 H

₀

• :

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 H

₀

• . You would have

• rejected H

₀

• :

b

• = 3.0 or H

₀

• :

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 H

₀

• :

d

• = c at significance level

n

• against H1:

d

• > c, for example, compare the test statistic t

_obs

• = (

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:

t_n,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 H₀: There are no differences in
log-logistic response among the two treatments The full model is y_j = d_j + (a_j – d_j)/(1 + exp{b_j
log(x/g_j)}) + 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 H₀.
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 y_j = 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 F_obs = (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

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

To test whether there is any difference among the two treatments,
calculate the sum of square reduction test: F_obs = { (7334.5 – 4465.2 ) / (108
– 104 ) } / 42.9342  = 16.70754 and the p-value under a F_4,104 distribution is
1.32E-10, H₀ 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 y_j = d_j
+ (a_j – d_j)/(1
+ exp{b_j log(x/g_j)})
+ e and the reduced model y_j = d_j + (a_j – d_j)/(1 + exp{b_j 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
y_j = d_j + (a_j
– d_j)/(1 + exp{b_j
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 y_j = d_j + (a_j – d_j)/(1 + exp{b_j
log(x/ g)}) + e. For treatment 2 (Tx=2) is true and the model
becomes y_j = d_j + (a_j
– d_j)/(1 + exp{b_j
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 H₀: g₁ = g₂. Here is the code (gammaD denotes D):

proc nlin data=weeds method=marquardt;
  parameters alpha1=100 delta1=4 beta1=2.0 gamma1 =0.2
             alpha2=100
delta2=4 beta2=2.0 gammaD=0;
  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
gammaD          -0.0693       0.0164     -0.1018     -0.0368

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
             alphaD=0
deltaD=0.0 betaD=0.0 gammaD=0;
  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
alphaD           1.9712       2.6662     -3.3161      7.2585
deltaD           4.3530       5.7602     -7.0698     15.7757
betaD           -0.7784       0.3062     -1.3857     -0.1712
gammaD          -0.0693       0.0164     -0.1018     -0.0368

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