SAS Library
Nonlinear Regression in SAS


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.
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

    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
    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.

    3.3. Using Expressions inside NLIN

    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;

    3.4.  Fixing a parameter 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 3.6.  Choosing the Fitting Algorithm proc nlin data=whatever method=newton;
      <statements>
    run;
    proc nlin data=whatever method=marquardt;
      <statements>
    run;

    3.7.  Calculating predicted values and their confidence intervals 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;
    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

    3.8.  Calculating predicted values for plotting purposes data filler;
      do rate = 0.05 to 0.8 by 0.05;
        predict=1;
        output;
      end;
    run;

    data fitthis;
      set weeds filler;
    run;

    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;

    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

    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

    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
                 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.

    How to cite this page

    Report an error on this page or leave a comment

    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.