7.2.4. Nonlinear Regression
7.2.4.0. Overview
The Nonlinear Regression procedure provides a leastsquares method of fitting a userspecified function to a suitable data set. This regression function would usually be based on a theoretical model of the system under analysis, and can be written in terms of any number of independent variables and any number of parameters (subject to program restrictions). The user provides initial estimates for the parameter values, and the program then adjusts these over a number of iterations in order to obtain the best fit to the data (that is, minimise the residual sum of squares). The program ends by displaying a large number of output options regarding the parameters and fitted curve.
In addition to providing numerical values for the parameters and some related statistics, Nonlinear Regression can also be used to compare regression models. For example, with data describing an exponential decay curve, it is possible to assess (in a statistical sense) if the curve is monophasic:
or biphasic:
The same method can be used to compare two or more sets of data and determine if differences exist between them. This is explained in detail later in this section. First we shall look at the dialogues and facilities of this procedure.
7.2.4.1. Nonlinear Regression Variable Selection
Dependent: Selection of a dependent variable containing numeric data is compulsory. Only one dependent variable can be selected.
Variable: At least one independent variable must be selected. It is not necessary for all independent variable columns selected in the Variables Selected list to be used in the analysis. The regression function is written as an expression involving variables C1, C2, … or their labels, and parameters P(1), P(2), … or their labels. Therefore, which columns are used depends on how the regression function is constructed.
Weight: A column of spreadsheet containing numeric data can be selected as weights. Cases with the highest weight value have the greatest effect on the fitted curve and those with the lowest weight have the least effect. If a weight variable is specified, then all entries in that data column must be greater than zero and not greater than 1.0E+10.
Factor: If a factor column is selected, the program can fit two or more curves simultaneously on the subgroups defined by the levels of the factor. This is useful when making comparisons between sets of data. The program can deal with up to eight curves, using the first eight levels of a factor.
7.2.4.2. Nonlinear Regression Setup Dialogue
This dialogue allows you to enter all the information needed to specify how the regression is to be performed. All this information can be saved in a Nonlinear Regression setup file and loaded back when needed, which makes repeating the same or similar analyses on several data files much easier.
7.2.4.2.1. Buttons
The four buttons provided have the following tasks:
Param: Shortened for parameters. It will activate the Parameter Setup Dialogue.
Load: This will invoke the standard windows dialogue for opening a file. You can load a previously saved Nonlinear Regression setup file. The default extension is .USN.
The names and column numbers of all variables selected from the spreadsheet (dependent variable, all independent variables and if any, weight and factor variables) are recorded in the setup file. When the setup file is subsequently loaded, these will be checked against the actual selections made from the Variable Selection Dialogue so that there are no inconsistencies between the setup file and the variables in the spreadsheet. If there are inconsistencies, these will be reported, but it will still be possible to run an analysis. In this case, the selections made from the Variable Selection Dialogue will override the definitions in the setup file and it will be up to you to ensure that the variables used in various functions do refer to the correct columns in the spreadsheet.
Save: All current selections made in the Setup and Parameter dialogues, as well as the variable assignments made in the Variable Selection Dialogue will be saved in a binary file. The default extension is .USN and you are strongly recommended to keep that unchanged.
The [Save] option always saves the current values of parameters. Therefore, if this option is selected after the convergence has been achieved then the file may be saved with convergence parameters. Or more importantly, for particularly difficult problems, after breaking the run after a number of iterations, the current values of the parameters can be saved in the setup file making it possible to continue to run the model subsequently from where you left.
Clear: This will clear all controls in both Setup and Parameter dialogues for a fresh restart.
7.2.4.2.2. Controls
Title: This is a line of text you may enter to annotate a particular Nonlinear Regression setup and will also be printed as a sub title (after the procedure name) in the output.
Function: The regression function is entered as an expression involving the independent variable(s) and parameters. The independent variables are referred to with their labels or column numbers (C1, C2, etc.), and the parameters as P(1), P(2), etc.
As an example, to fit an exponential decay curve with the independent variable in C1, the function can be written as:
P(1)*Exp(C1*P(2))
If C1 has a label TIME and parameters 1 and 2 named ACTIVITY0 and RATE respectively, then the expression can be written as:
ACTIVITY0*Exp(TIME*RATE)
Weighting Function: Iterative reweighting can be used, which means that at each iteration the weights to be used in that iteration are calculated according to an expression (the weighting function) which is supplied by the user. This function is written in the same way as the regression function, except that a variable containing the fitted values V(0) can also be used.
The weight values generated in this way must be greater than zero and not greater than 1.0E+10, otherwise an error will be generated. The weighting function will only be used if the iterative reweighting option (see below) is enabled.
If a weights variable has been selected and iterative reweighting enabled, then an error will be reported. They cannot be used at the same time.
XAxis: A transform expression for the horizontal axis of the graphical display must be entered. In the simplest case this will be a single variable (e.g. TIME or C1 for the example described above), but any expression containing more than one independent variable can be used. This expression must not include parameters or the dependent variable.
YAxis: As above but applies to the vertical axis of the graphical display. The dependent variable must be included in the Yaxis expression.
Max Iterations: This is the maximum number of iterations to be performed (default value is 100). It is seldom necessary to alter this value, since the program will terminate when the best fit has been obtained and the residual sum of squares stops changing from one iteration to the next (see Converge below). Also, the program can be interrupted at any time. Setting the maximum number of iterations to zero has the same effect as selecting Run Mode = FINISH (see below). A negative value will cause an error to be reported.
Max Halving: At each iteration the program calculates a correction vector for each active parameter which is used to modify that parameter value. The new parameter values generated will usually result in a decrease in the residual sum of squares (RSS), but if the RSS increases then the correction vectors are halved and their signs reversed (subject to lower/upper limits on the parameters) until a decrease in the RSS is obtained or the maximum number of step halving has been performed. The default value is five and the minimum is zero. A negative value will cause an error to be reported. The setting for Max Halving is only relevant if Run Mode = STANDARD (see below).
Converge: The program will terminate when the absolute fractional change in the RSS over three successive iterations is less than the value specified for Converge. The default value is 0.0001, meaning that if the RSS changes by less than 0.01% over three iterations then it is assumed that the best fit (minimum RSS) has been obtained. The final calculations are made and the results are then displayed.
Tolerance: The pivoting tolerance value should not normally be changed from its default value of 1.0E11. However, when fitting a function which is more complex (has more parameters) than necessary the program may become unstable, evidenced by the RSS increasing rather than decreasing. This does not immediately indicate a program failure, but if it persists for 10 iterations or more, then some corrective action will probably be needed. This can usually be accomplished by specifying different initial estimates and / or supplying appropriate lower / upper limits for the parameters. If the instability is still present then increasing the tolerance value to 1.0E5 may be helpful.
Run Mode: This determines how the program operates. The four options available are:
1) FAST: The regression is performed using the MarquardtLevenburg algorithm, which in most cases is much faster than the conventional GaussNewton method. However, the MarquardtLevenburg method as implemented in this program is less suitable for difficult regression analyses and will fail if there is an excessively high correlation between two or more parameters. If this happens the program will automatically switch to STANDARD mode. Difficult regression analyses may include those in which a) for certain combinations of parameter values, a large change in RSS accompanies a very small change in the value(s) of one or more parameters; b) an inappropriate regression model has been specified, which might be done intentionally in order to compare it statistically with another model; and c) there is excessive correlation between the parameters.
2) STANDARD: The regression is performed using a slower but much more stable method based on the conventional GaussNewton algorithm, which can cope with difficult analyses as described above. Since the program will change from FAST to STANDARD mode itself when needed, it is not necessary to start in this mode. However, if an inappropriate regression model (particularly one that is more complex or has more parameters than necessary) is being used, then beginning the analysis in STANDARD mode may give somewhat different parameter values; the extra parameters tend not to change from their initial values in this mode.
3) FINISH: The final calculations are made on the current parameter values and the results are displayed. These statistics (the parameter standard errors, correlation matrix and predicted standard deviations of the fitted curve) are only valid if a good fit to the data has been obtained.
4) MANUAL: In this mode the user enters a set of parameter values and the program then shows the corresponding fitted curve and RSS. This process repeats until the user exits from this part of the program. The program does not alter the parameter values, except for evaluation of any parameter constraint expressions that have been defined. If weighting or iterative reweighting has been enabled then the weighted RSS is calculated. This mode can be used to obtain initial parameter estimates for the regression, or to examine the effect of changing the parameter values on the fitted curve.
Double Check: If the regression converges in FAST mode and the fitted curve satisfies a builtin goodnessoffit test (see 7.2.4.6. Technical Details) then the final calculations are made and the results displayed. If the fitted curve fails the goodnessoffit test then the STANDARD mode is entered automatically. However, if Double Check is enabled then the program will switch to STANDARD mode to see if any further improvement in fit can be obtained (which does sometimes happen). The default settings (Run Mode = FAST and Double Check = On) are recommended for most applications.
Iterative Reweighting: this enables iterative reweighting as explained above for weighting function.
7.2.4.2.3. Parameter Setup Dialogue
The program will determine the number of parameters used in the regression from the number of rows filled in by the user. The minimum is one and the maximum is fifty.
Label: These are optional, and can be used instead of the parameter numbers P(1), P(2), etc.
Value: In Nonlinear Regression the user must specify an initial value for each parameter which the program will then adjust in order to minimise the residual sum of squares. These initial estimates should be close to the true (or bestfit) values, but exactly how close they need to be depends on various factors such as the amount or error in the data and the complexity of the regression function. The initial values might be determined from the theory underlying the Regression Analysis, prior knowledge of the system, or graphical examination of the data.
The initial parameter values may be set to zero, but this will slow the program down quite markedly. A blank field does not have a zero value.
On completion of the Regression Analysis, the final parameter values are written back to the Parameter Setup Dialogue where they can be edited if necessary and used as starting values for further analysis.
Minimum: This is the minimum value that each parameter can take during the regression. A blank field means that no lower limit is imposed.
Maximum: This sets the maximum value for each parameter.
Although optional, it is recommended that lower and upper limits should be assigned for each parameter as a precaution against arithmetic errors. For example, if the regression function includes the square root of a parameter, then that parameter should have a lower limit of zero. Also, if the theory underlying Regression Analysis indicates that some parameter could not have a negative value, then that parameter would also have a lower limit of zero. Otherwise, quite broad limits can be set; parameters should not be restricted to a very narrow range unless there is special reason to do so.
If lower and / or upper limits have been set, then the initial value must be within the defined range otherwise an error will be reported when the Regression Analysis is started.
Constraint: This is the constraint expression referred to above. The expression is written in terms of other parameters only and must not include any of the input variables. The constraint expression will only be used if the corresponding parameter status is set to Constrained.
Status: This control determines how the parameter will be handled during the regression. There are three possible values:
1) Free: The value of the parameter will be adjusted in an attempt to obtain the best fit.
2) Fixed: The parameter value is fixed (constant) and it will not change.
3) Constrained: The value of the parameter is determined by other parameter values according to an expression (a parameter constraint expression) which you must supply. For example, if you wanted to perform the regression such that parameter 5 always had the same value as parameter 3, then parameter 5 would have status Constrained and the constraint expression should be P(3). If a parameter has status Constrained then the corresponding parameter constraint expression must be defined.
Parameters with status Free are termed active parameters. There must be at least one active parameter, otherwise an error will be reported.
7.2.4.3. Running the Program
After filling in the Nonlinear Regression Setup Dialogue and the Parameter Setup Dialogue, click [Next] to run the model. The program will first check all the entries and report any inconsistencies found. Otherwise a new window is opened showing a bivariate plot of data according to the Xaxis and Yaxis expressions defined in the Nonlinear Regression Setup Dialogue. At each iteration, the progress of fitting can be followed visually. The same window will also display the current values of number of iterations, residual sum of squares, as well as the current values of the parameters.
It is possible to interrupt or break the program while it is running by pressing <Escape>. However, the response may not be immediate, as the program will always finish the current iteration before stopping. When this happens, the Nonlinear Regression Setup Dialogue will be displayed. The current values can be saved to a regression setup file to continue later on, or change some of the parameters manually. When [Continue] is clicked, the program will resume iterations. Changing the regression function will result in a complete resetting of the regression run.
If the iterations continue with no sign of convergence, you may break the run, select the FINISH mode from the Nonlinear Regression Setup Dialogue and force the program to perform calculations necessary for output options.
7.2.4.4. Nonlinear Regression Output Options
After convergence has been achieved or the program is stopped before convergence and Run Mode = FINISH has been selected, the Output Options Dialogue will be displayed. For a sample printout of the full Nonlinear Regression output (see 7.2.4.5. Nonlinear Regression Examples).
The Actual and Fitted Values, Residuals, Confidence Intervals for Actual Y Values and Interpolation options that existed in earlier version of UNISTAT are now combined under the Case (Diagnostic) Statistics option.
Regression Results: The main regression results include a report on whether the convergence has been achieved and the result of a goodnessoffit test, which will be either OK or Failed test.
All values related to parameters are displayed in a table. These include minimum, maximum, coefficient estimate, standard error and status (Free, Fixed or Constrained). Confidence intervals for regression coefficients are also displayed for the userdefined confidence level. The entire table can be sent to the spreadsheet for further analysis.
Statistics reported include residual sum of squares, root mean square error, number of parameters, number of active parameters, number of cases, degrees of freedom (number of cases minus number of active parameters), Rsquared and DurbinWatson statistic.
Standard errors are only available for parameters with status Free. If a factor column has been selected, then most of the statistics will be displayed for each sub group separately.
ANOVA of Regression: The ANOVA table compares the regression model with the null hypothesis that “Y = constant”.
Correlation Matrix of Regression Coefficients: Only parameters with status Free are included in this table.
Covariance Matrix of Regression Coefficients: Only parameters with status Free are included in this table.
Case (Diagnostic) Statistics:
Predictions (Interpolations): As of this version of UNISTAT, the Interpolation output option is removed and the predicted Y values can be generated as in Linear Regression (see 7.2.1.2.2. Linear Regression Case Output). If, for a case, all independent variables are nonmissing, but only the dependent variable is missing, a value will be predicted for the fitted Y value from the estimated regression equation. When a case is predicted, its label will be prefixed by an asterisk (*).
Statistics available under Case (Diagnostic) Statistics option are as follows.
Actual Y: Observed values of the dependent variable.
Fitted Y: Estimated values of the dependent variable.
Confidence Intervals for Actual Yvalues: Confidence intervals for actual Y values are displayed. The default confidence level is 95%, but any value can be entered from the Variable Selection Dialogue.
Where is the critical value from the tdistribution for an α / 2 level of significance and 1 degree of freedom.
Standard Error of Fitted:
Residuals:
Plot of Actual and Fitted Values: Actual and fitted Y values, as well as their confidence intervals, are plotted. The Xaxis variable is as selected in the Nonlinear Regression Setup Dialogue.
As with plots in Linear Regression and Polynomial Regression, the data points here are also interactive, but pressing <Delete> when a point is highlighted will not omit this point and rerun the analysis.
Plot of Residuals: A scatter plot of residuals is drawn against the dependent variable.
7.2.4.5. Nonlinear Regression Examples
7.2.4.5.1. Example 1: Radioimmunoassay
This is a simple example based on a calibration curve for a type of assay. Open NONLINRG and select Statistics 1 → Regression Analysis → Nonlinear Regression. The data consists of 14 cases and two variables STANDARD and BOUND.
Select STANDARD (C1) as the independent variable by clicking [Variable], and BOUND (C2) as the [Dependent] variable. Then clicking [Next], the Nonlinear Regression Setup Dialogue will be opened. Click on the [Load] button and load the NONLIN1.USN setup file. The function involves three parameters:
which is entered as:
1/((P(1)*STANDARD+P(2))+P(3)
or could be entered as:
1/((P(1)*V(1)+P(2))+P(3)
The three parameters have initial values of 0.01, 0.1 and 0.1 respectively, all have status Free, and no lower or upper limits are set. With the default settings of Run Mode = FAST and Double Check = On, the program converges after 8 iterations, then switches to STANDARD mode and performs another 3 iterations, at which point it converges again. The following are the main results:
Nonlinear Regression
Regression Results
EXAMPLE 1: RADIOIMMUNOASSAY
BOUND = 1/(P(1)*STANDARD+P(2))+P(3)
Convergence achieved at iteration 11, No of halving = 0
Goodness of Fit: OK
Parameter 
Minimum 
Maximum 
Estimate 
Standard Error 
Status 
Lower 95% 
Upper 95% 
Par1 


0.0279 
0.0019 
Free 
0.0237 
0.0321 
Par2 


1.0914 
0.0171 
Free 
1.0538 
1.1291 
Par3 


0.0160 
0.0153 
Free 
0.0177 
0.0496 
Residual Sum of Squares = 
0.0018 
Root Mean Square Error = 
0.0126 
Number of Parameters = 
3 
No of Active Parameters = 
3 
Number of Cases = 
14 
Degrees of Freedom = 
11 
Rsquared = 
0.9983 
DurbinWatson Statistic = 
2.5336 
ANOVA of Regression
Due To 
Sum of Squares 
DoF 
Mean Square 
FStat 
Prob 
Regression 
1.005 
2 
0.503 
3141.719 
0.0000 
Error 
0.002 
11 
0.000 


Total 
1.007 
13 
0.077 


Correlation Matrix of Regression Coefficients

1 
2 
3 
1 
1.0000 
0.7328 
0.9336 
2 
0.7328 
1.0000 
0.8761 
3 
0.9336 
0.8761 
1.0000 
Covariance Matrix of Regression Coefficients

1 
2 
3 
1 
0.0000 
0.0000 
0.0000 
2 
0.0000 
0.0003 
0.0002 
3 
0.0000 
0.0002 
0.0002 
Case (Diagnostic) Statistics

Actual Y 
Fitted Y 
95% lb Actual Y 
95% ub Actual Y 
Standard Error of Fitted 
Residuals 
1 
0.9472 
0.9322 
0.8638 
1.0006 
0.0054 
0.0150 
2 
0.9255 
0.9322 
0.8638 
1.0006 
0.0054 
0.0067 
3 
0.8201 
0.8284 
0.7637 
0.8930 
0.0051 
0.0083 
4 
0.8312 
0.8284 
0.7637 
0.8930 
0.0051 
0.0028 
5 
0.7254 
0.7457 
0.6742 
0.8172 
0.0056 
0.0203 
6 
0.7530 
0.7457 
0.6742 
0.8172 
0.0056 
0.0073 
7 
0.5679 
0.5750 
0.4885 
0.6614 
0.0068 
0.0071 
8 
0.5910 
0.5750 
0.4885 
0.6614 
0.0068 
0.0160 
9 
0.4324 
0.4181 
0.3315 
0.5048 
0.0068 
0.0143 
10 
0.4105 
0.4181 
0.3315 
0.5048 
0.0068 
0.0076 
11 
0.2782 
0.2736 
0.2024 
0.3448 
0.0056 
0.0046 
12 
0.2634 
0.2736 
0.2024 
0.3448 
0.0056 
0.0102 
13 
0.1772 
0.1659 
0.1161 
0.2156 
0.0039 
0.0113 
14 
0.1546 
0.1659 
0.1161 
0.2156 
0.0039 
0.0113 
7.2.4.5.2. Example 2: Biphasic Exponential Decay
This example will illustrate the comparison of regression models. Open NONLINRG and select Statistics 1 → Regression Analysis → Nonlinear Regression. The data for this example consists of 21 cases and 3 variables.
Select TIME (C3) as the independent [Variable], ACTIVITY (C4) as the [Dependent] variable, and WEIGHT (C5) as the [Weight] variable. Click [OK] and after the Nonlinear Regression Setup Dialogue is displayed click on the [Load] button and load the NONLIN2.USN setup file.
The function involves 4 parameters (INIT1, RATE1, INIT2, RATE2) and describes the following biphasic exponential decay curve:
ACTIVITY = INIT1*Exp(RATE1*TIME)+INIT2*Exp(RATE2*TIME)
which could be entered, without variable or parameter names, as
P(1)*Exp(P(2)*C1)+P(3)*Exp(P(4)*C1)
The parameters have initial values of 10, 0.1, 5 and 0.01 respectively, all have status Free, and no lower or upper limits are set. The program will perform 8 iterations in FAST mode, then 3 iterations in STANDARD mode, and then the results are displayed.
Nonlinear Regression
Regression Results
EXAMPLE 2: BIPHASIC EXPONENTIAL DECAY
ACTIVITY = INIT1*EXP(RATE1*TIME)+INIT2*EXP(RATE2*TIME)
Convergence achieved at iteration 11, No of halving = 0
Goodness of Fit: OK
Parameter 
Minimum 
Maximum 
Estimate 
Standard Error 
Status 
Lower 95% 
Upper 95% 
INIT1 


14.0304 
1.5996 
Free 
10.6554 
17.4054 
RATE1 


0.2627 
0.0320 
Free 
0.3303 
0.1952 
INIT2 


7.4943 
0.1836 
Free 
7.1071 
7.8816 
RATE2 


0.0034 
0.0002 
Free 
0.0039 
0.0029 
With normalised weights: 

Residual Sum of Squares = 
1.4317 
Root Mean Square Error = 
0.2902 
With absolute weights: 

Residual Sum of Squares = 
40.0617 
Root Mean Square Error = 
1.5351 
Number of Parameters = 
4 
No of Active Parameters = 
4 
Number of Cases = 
21 
Degrees of Freedom = 
17 
Rsquared = 
0.9863 
DurbinWatson Statistic = 
0.8836 
The rest of the output options are not reproduced here.
The output is essentially the same as for Example 1, except for weights. When weights are used, two RSS values are shown; the first is obtained with a normalised set of weights used for internal calculations and the second with the absolute weight values from the data file (or those calculated using the weighting function). The normalised weights have a mean value of 1.0 and are proportional to the absolute values. Since the RMS error is obtained from the RSS, there are also two entries for RMSerror. If the weights only indicate reliability of each data point relative to the others, it would be more appropriate to use the normalised RSS for reporting purposes or further calculations.
Although the fitted curve passes the goodnessoffit test, the graphical display (plot of ACTIVITY against TIME) might suggest a slight systematic deviation of the data points about the curve. Interpretation of this would depend on the theoretical basis for the analysis, but suppose that a six parameter model was considered possible:
ACTIVITY = INIT1*Exp(RATE1*TIME)+INIT2
*Exp(RATE2*TIME)+INIT3*Exp(RATE3*TIME)
To establish if there is any statistical basis for accepting the six parameter model in preference to the fourparameter model, it is necessary to perform the regression again using the more complex function and then compare both sets of results on the basis of change in RSS with change in number of degrees of freedom for the regression using an extrasumofsquares test (sometimes referred to as an Fratio test).
Altering the regression setup is quite simple: assign labels and initial values to parameters 5 & 6 of INIT3 and RATE3 and 1.0 and 0.001 respectively, and include these parameters in the regression function as shown above. Also, assign a maximum value of zero for each RATE parameter. The regression can now be performed (taking about 60 iterations) using the six parameter model. The results will be as follows.
Nonlinear Regression
Regression results
EXAMPLE 2: BIPHASIC EXPONENTIAL DECAY
Convergence achieved at iteration 52, no of halving = 0
Goodness of Fit Test: OK
Parameter 
Minimum 
Maximum 
Estimate 
Std Error 
INIT1 


16.8864 
2.2845 
RATE1 

0 
0.4513 
0.0806 
INIT2 


3.1163 
0.5334 
RATE2 

0 
0.0438 
0.0168 
INIT3 


6.4984 
0.4217 
RATE3 

0 
0.0024 
0.0004 
With normalised weights: 

Residual Sum of Squares = 
0.4571 
Root Mean Square Error = 
0.1746 
With absolute weights: 

Residual Sum of Squares = 
12.7893 
Root Mean Square Error = 
0.9234 
Number of Parameters = 
6 
No of Active Parameters = 
6 
Number of Cases = 
21 
Degrees of Freedom = 
15 
Rsquared = 
0.9338 
DurbinWatson statistic = 
2.0314 
The null hypothesis is that “there is no difference in RSS between the two models”. The RSS has actually decreased by about 60%, but it is necessary to assign a significance level to this (the probability of rejecting the null hypothesis when it is true). For the following calculations it does not matter if the RSS values were obtained with normalised or absolute weights, except that you should use the same one throughout (normalised values are used below – this method is equally applicable to comparison of unweighted regression analyses). Defining;
Model 1 (4 parameters): RSS1 = 1.43174 with DF1 = 17
Model 2 (6 parameters): RSS2 = 0.45707 with DF2 = 15
an Fstatistic is calculated as:
with (DF1 – DF2) and DF2 degrees of freedom.
This gives F = 15.993 with 2 & 15 degrees of freedom. From statistical tables the associated pvalue is less than 0.001, and consequently the null hypothesis can be rejected. That is, the sixparameter model might be said to give a significantly better fit.
In fact, a fiveparameter model:
ACTIVITY = INIT1*Exp(RATE1*TIME)+INIT2
*Exp(RATE2*TIME)+CONSTANT
results in a significant decrease in RSS compared with the fourparameter model (P < 0.001), but there is only a slight difference in RSS between this and the sixparameter model (P > 0.05). This fiveparameter model can be specified either by rewriting the regression function in terms of five parameters, or using the sixparameter model with P(6) = 0 and with status Fixed (parameter 6 is held constant with a value of zero).
7.2.4.5.3. Example 3: DoseResponse Curves
This example will show how the program can be used to fit two (or more) curves simultaneously and to make comparisons between them using the Fratio method.
Open NONLINRG and select Statistics 1 → Regression Analysis → Nonlinear Regression. The data for this example consists of 84 cases and 4 variables. Select GROUP (C6) as the Factor, DOSE (C7) and CELLNO (C8) as the independent [Variable]s and RESPONSE (C9) as the [Dependent] variable. Click [OK] to proceed and when the Nonlinear Regression Setup Dialogue is displayed click on the [Load] button and load NONLIN3.USN regression setup file.
Suppose that in a pharmacology experiment, batches of cells are grown, subjected to some treatment or left untreated as the controls, then exposed to different concentrations of a drug and the response of the cells to that drug measured. The purpose of the experiment is to determine if the treatment has any effect on the drug response and if so, to describe that effect. The doseresponse curve is sigmoid and can be described by the logistic function:
where:
DS1=DOSE^SLOPE
ES1=ED50*SLOPE
DOSE represents the dose of the drug (an independent variable) and RESPONSE is the measured response (the dependent variable). The four parameters are BASAL (the minimum or resting response), ACTIVE (the activatible response or maximum rise above basal), ED50 (the concentration of drug giving halfmaximal activation) and SLOPE (which describes the steepness of the increase in response). Actually, ED50 is the logarithm of the dose giving halfmaximal activation; it is appropriate to estimate the logarithm of this parameter since it has a lognormal distribution (a plot of ACTIVITY against Log(DOSE) is symmetrical about the point DOSE = ED50).
Further, suppose that three experiments are performed, each involving a paired control and treated group. The cell preparation is homogeneous within each of the three experiments but can vary between experiments. The response is assumed to be proportional to the number of cells present, which is also measured (CELLNO, another independent variable). The function is now:
and therefore the parameters BASAL and ACTIVE represent the response per unit cell number.
Since the data file represents two curves (control and treated), it is necessary to indicate which cases belong to each curve. A data grouping variable (GROUP) is included in the data file with values of 1 (control) or 2 (treated). The order in which the data points appear in the file has no effect on the regression.
Since a fourparameter function is used to describe one doseresponse curve and the data describes two such curves, the regression will involve eight parameters. These are named BASAL:1, ACTIVE:1, ED50:1, SLOPE:1 (for group 1), and BASAL:2, ACTIVE:2, ED50:2 and SLOPE:2 (for group 2). It is now necessary to write the regression function in such a way that the value of GROUP is used to select which parameters are used in evaluating the function. The function is written:
(P((GROUP1)*4+1)+P((GROUP1)*4+2)*
(DOSE^P((GROUP1)*4+4)/10^(P((GROUP1)*4+3)*
P((GROUP1)*4+4)))/(1+DOSE^P((GROUP1)*4+4)/
10^(P((GROUP1)*4+3)*P((GROUP1)*4+4))))*CELLNO
Factor variable GROUP is used to calculate the parameter subscripts. When GROUP = 1, the function is equivalent to:
(P(1)+P(2)*(DOSE^P(4)/10^(P(3)*P(4)))/
(1+DOSE^P(4)/10^(P(3)*P(4))))*CELLNO
and when GROUP = 2 it is equivalent to:
(P(5)+P(6)*(DOSE^P(8)/10^(P(7)*P(8)))/
(1+DOSE^P(8)/10^(P(7)*P(8))))*CELLNO
In the setup file, the following parameter constraints are defined but not enabled (initially all parameters have status Free);
BASAL:2 = BASAL:1
ACTIVE:2 = ACTIVE:1
ED50:2 = ED50:1
SLOPE:2 = SLOPE:1
These will be used later in making comparisons between the two curves.
Finally, the transforms for the graphical display are defined:
XAxis = Log(DOSE)
YAxis = RESPONSE/CELLNO
The first regression analysis can now be performed. For the purposes of the following discussion, this model (8 active parameters) is model 1. The results are as follows:
Nonlinear Regression
Regression Results
EXAMPLE 3: DOSERESPONSE CURVES / MODEL 1
RESPONSE = (P((GROUP1)*4+1)+P((GROUP1)*4+2)*(DOSE^P((GROUP1)*4+4)/10^(P((GROUP1)*4+3)*P((GROUP1)*4+4)))/(1+DOSE^P((GROUP1)*4+4)/10^(P((GROUP1)*4+3)*P((GROUP1)*4+4))))*CELLNO
Convergence achieved at iteration 16, No of halving = 0
Goodness of Fit: OK
Parameter 
Minimum 
Maximum 
Estimate 
Standard Error 
Status 
Lower 95% 
Upper 95% 
BASAL:1 
0.0000 
10.0000 
0.8215 
0.1194 
Free 
0.5837 
1.0592 
ACTIVE:1 
0.0000 
10.0000 
3.4658 
0.2280 
Free 
3.0116 
3.9200 
ED50:1 
1.0000 
5.0000 
3.3174 
0.1232 
Free 
3.0720 
3.5629 
SLOPE:1 
0.1000 
5.0000 
0.8061 
0.1659 
Free 
0.4757 
1.1365 
BASAL:2 
0.0000 
10.0000 
0.8111 
0.1304 
Free 
0.5515 
1.0707 
ACTIVE:2 
0.0000 
10.0000 
3.3104 
0.1772 
Free 
2.9575 
3.6632 
ED50:2 
1.0000 
5.0000 
2.3012 
0.0923 
Free 
2.1173 
2.4851 
SLOPE:2 
0.1000 
5.0000 
1.1878 
0.2613 
Free 
0.6674 
1.7083 
Group 
No 
RSSQ(n/w) 
Rsquared 
DW stat 
Fit 
1 
42 
207905.2146 
0.9795 
1.9323 
OK 
2 
42 
226867.9214 
0.9826 
1.9267 
OK 
Residual Sum of Squares = 
434773.1359 
Root Mean Square Error = 
75.6353 
Number of Parameters = 
8 
No of Active Parameters = 
8 
Number of Cases = 
84 
Degrees of Freedom = 
76 
Rsquared = 
0.9381 
DurbinWatson Statistic = 
1.9277 
The rest of the output options are not reproduced here.
The output is largely as described previously except that a new list is included in the Regression Results for each data group; the number of cases, the RSS and the result of the goodnessoffit test. If weighting has been used then the RSS values are calculated using normalised weights as explained above.
From the graphical display there appears to be little difference in basal or maximal (basal plus activatible) responses, but there does seem to be a difference in ED50 and perhaps also in slope since, between minimum and maximum responses, the two curves are not completely parallel. To make a preliminary comparison between the two curves, we can determine if the difference between parameter values exceeds two standard errors of the difference, which would suggest a significant difference. That is, for parameters having values Pa and Pb with standard errors Ea and Eb, if:
then a genuine difference is likely. Using this method, the only apparent difference is between parameters 3 (ED50:1) and 7 (ED50:2). However, the Fratio test is more reliable, and should be used in preference as demonstrated below.
To test the null hypothesis that “there is no difference of any kind between control and treated groups” (model 2), all four constraints are enabled. This gives a model with four active parameters, and results in RSS = 801471.4 with DF = 80. The goodnessof fit test returns Fail both for the entire data set and for both groups, which certainly indicates that this model is not appropriate. To compare models 1 & 2, the Fstatistic is calculated as:
with (8076) and 76 degrees of freedom. This gives F = 16.025 with 4 & 76 degrees of freedom and a pvalue of less than 0.001. Therefore the null hypothesis is rejected and a genuine difference between the curves is established.
The preliminary analysis suggested no difference in BASAL or ACTIVE parameters between control and treated groups. To test this properly, the constraints for parameters 5 and 6 are enabled (thus BASAL:2 = BASAL:1 and ACTIVE:2 = ACTIVE:1. This is model 3, and results in RSS = 438789.1 with DF = 79, and all goodnessoffit tests are passed. Comparison of models 1 & 3 gives an Fstatistic of 0.351 with 2 & 76 degrees of freedom, and a corresponding pvalue of greater than 0.1. Thus the null hypothesis “no difference in BASAL or ACTIVE parameters” is not rejected.
With constraints BASAL:2 = BASAL:1, ACTIVE:2 = ACTIVE:1 and SLOPE:2 = SLOPE:1 (model 4), the RSS is 446337.8 with DF = 79. Comparison of this with model 1 gives F = 0.927 with 3 & 76 degrees of freedom and a pvalue of greater than 0.1. The null hypothesis “no difference in BASAL, ACTIVE or SLOPE parameters” is not rejected. It is sometimes necessary to carry out the comparisons in a stepwise manner. For example, having established no difference in BASAL or ACTIVE parameters, models 3 & 4 could be compared to further establish no difference in SLOPE parameters (which would give F = 1.341 with 1 & 78 degrees of freedom and p > 0.1, which is similar to the comparison of models 1 & 4).
Finally, with the single constraint ED50:2 = ED50:1 (model 5, null hypothesis “no difference in ED50 parameters”), we obtain RSS = 660274.5 with DF = 77. Comparison of this with model 1 gives F = 39.421 with 1 & 76 degrees of freedom and a pvalue of less than 0.001, and therefore the null hypothesis is rejected.
Thus the preliminary analysis was correct in that the only difference between control and treated groups is a change in ED50. The null hypothesis for the experiment would be that there is no difference between control and treated groups (model 2). Comparison of this with model 4 (difference in ED50 only) gives F = 62.857 with 1 & 79 degrees of freedom and a pvalue of less than 0.001. Using the results obtained with model 4, it can be stated that the treatment had no effect on BASAL, ACTIVE or SLOPE parameters (0.835 0.085, 3.319 0.129, and 1.044 0.156 respectively, p > 0.1) but reduced the ED50 from 3.253 0.091 to 2.322 0.091 (P < 0.001 using the Fvalue of 62.857).
7.2.4.6. Technical Details
The program seeks the set of parameter values that minimise the residual sum of squares (RSS):
or:
if weighting is used, where Y_{obs} and Y_{pred} are the observed and predicted values of the dependent variable, w is the weight, and the summation is over all cases. The predicted value Y_{pred} is obtained by evaluation of the regression function f(v, p), where v represents the set of independent variables and p the parameters.
The GaussNewton method is implemented in two forms. In both, the program obtains for each case the residual (z, = Y_{obs} – Y_{pred}), and for each case and each active parameter (p_{i}) the partial derivative
Using the approximation:
where the summation is over all active parameters and q_{i} represents a correction vector for each such parameter, solution of the set of linear equations so generated provides an estimate for each q_{i}. These correction vectors are used to generate a new set of parameter values for use in the next iteration, and this process is repeated until a termination condition is satisfied.
The partial derivatives are obtained arithmetically by evaluation of f(v, p) at parameter values p and p delta. The delta value is computed for each parameter as follows:
1) If p = 0, then delta starts at 1.0E30 and increments in steps of 10 until a suitable change in RSS is obtained or delta = p/10 or a boundary constraint is met.
2) If p0, then delta starts at p/1000 and increments in steps of Sqr(10) until a suitable change in RSS is obtained or delta = p/Sqr(10) or a boundary constraint is met.
In FAST mode, the delta values are computed on entry and subsequently recomputed only for any parameter that changes sign. The partial derivatives are obtained at parameter values p and (p + delta) or (p – delta) depending on any boundary constraints. In STANDARD mode, the delta values are computed at each iteration and the partial derivatives are obtained at parameter values (p + delta) and (p – delta) subject to any boundary constraints. If a suitable change in RSS is not obtained then the parameter is omitted from the current iteration (in FAST mode the delta value will be recomputed at the next iteration).
In FAST mode, solution of the set of linear equations is achieved by a matrix inversion method incorporating the MarquardtLevenburg procedure provided that excessive correlation between parameters is not detected. In the presence of such correlation or in STANDARD mode, a stepwise multivariate regression method is used with omission of parameters until any excessive correlation is removed.
The correction vectors obtained are used to update the parameter values; in STANDARD mode if a reduction in RSS is not achieved with the new set of parameter values then the correction vectors are halved and their signs reversed (subject to any boundary constraints) until a decrease in RSS is obtained or the maximum number of step halvings has been performed. At the end of each iteration, the RSS is compared with the recorded bestfit value and if necessary the bestfit value and associated parameter values are updated.
Parameters subject to an equality constraint are not included in the steps described above, but the constraints are evaluated at all stages where the parameter values are altered (e.g. in calculating the partial derivatives). The constraint expressions are evaluated in the order in which the parameters are numbered; therefore a constraint expression for parameter b should not use the result of constraint imposed on parameter a unless b > a.
Iteration continues until a) the maximum number of iterations has been performed, b) the convergence criterion is satisfied, c) the user interrupts the program, or an error occurs (termination conditions). Convergence is achieved when the absolute fractional change in RSS over 3 successive iterations is less than a userspecified value (Converge, default 0.0001). In a) or b), the parameters are reset to the recorded bestfit values. If convergence occurs in FAST mode and either the Double Check facility is enabled or the fitted curve fails the goodnessoffit test, further iterations are performed in STANDARD mode until a termination condition occurs again. If termination conditions a) or b) are met in FAST mode, then the parameter standard errors, correlations and standard deviations of the ypredicted values are calculated using the partial derivatives obtained in the last iteration providing that the RSS in that iteration is within 1 Â± Convergence of the bestfit value and the goodnessoffit test is passed. Otherwise (or in STANDARD mode) the partial derivatives are recomputed after resetting the parameters to the bestfit values and before the parameter standard errors etc. are calculated.
Standard errors and correlations are not computed for parameters that have zero or negligible effect on the fitted curve, for which an excessive correlation is present, or which are held constant or are the subject of an equality constraint (status Fixed or Constrained).
The goodnessoffit test consists of two parts; a) a Runs Test assuming a binomial distribution with p = 0.5 is applied to the signs of the residuals and fail registered at α ≤ 0.005; and b) the weighted mean and standard deviation of the residuals is calculated and Fail registered if zero is out with the interval mean 2 x standard deviations. The test always returns Fail if the number of cases is less than five.
For further reading see Ratkowski, D. A. (1983) and Wadsworth, H. M. (1990).