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 = 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.
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.
-
- 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 R
2
-
- .
-
- In nonlinear regression, such a measure is unfortunately, not readily defined. One of the
-
- problems with the R
2
-
- definition is that it requires the presence of an
-
- intercept, which most nonlinear models do not have. A measure, relatively closely corresponding to R
2
-
- in the nonlinear case is Pseudo-R
2
-
- = 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:
-
-
- 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.
- 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.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
0
-
- :
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
0
-
- . You would have
-
- rejected H
0
-
- :
b
-
- = 3.0 or H
0
-
- :
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
0
-
- :
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:
tn,df(Res)
-
- .
d
-
- denotes the
-
- estimate of the parameter
d
- .
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
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: 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.