The NLIN Procedure

Example 69.4 Affecting Curvature through Parameterization

The work of Ratkowsky (1983, 1990) has brought into focus the importance of close-to-linear behavior of parameters in nonlinear regression models. The curvature in a nonlinear model consists of two components: the intrinsic curvature and the parameter-effects curvature. See the section Relative Curvature Measures of Nonlinearity for details. Intrinsic curvature expresses the degree to which the nonlinear model bends as values of the parameters change. This is not the same as the curviness of the model as a function of the covariates (the x variables). Intrinsic curvature is a function of the type of model you are fitting and the data. This curvature component cannot be affected by reparameterization of the model. According to Ratkowsky (1983), the intrinsic curvature component is typically smaller than the parameter-effects curvature, which can be affected by altering the parameterization of the model.

In models with low curvature, the nonlinear least squares parameter estimators behave similarly to least squares estimators in linear regression models, which have a number of desirable properties. If the model is correct, they are best linear unbiased estimators and are normally distributed if the model errors are normal (otherwise they are asymptotically normal). As you lower the curvature of a nonlinear model, you can expect that the parameter estimators approach the behavior of the linear regression model estimators: they behave "close to linear."

This example uses a simple data set and a commonly applied model for dose-response relationships to examine how the parameter-effects curvature can be reduced. The statistics by which an estimator’s behavior is judged are Box’s bias (Box, 1971) and Hougaard’s measure of skewness (Hougaard, 1982, 1985).

The log-logistic model

$\mr{E}[Y|x] = \delta + \frac{\alpha - \delta }{1+\gamma \exp \left\{ \beta \ln (x)\right\} }$

is a popular model to express the response Y as a function of dose x. The response is bounded between the asymptotes $\alpha$ and $\delta$ . The term in the denominator governs the transition between the asymptotes and depends on two parameters, $\gamma$ and $\beta$ . The log-logistic model can be viewed as a member of a broader class of dose-response functions, those relying on switch-on or switch-off mechanisms (see, for example, Schabenberger and Pierce 2002, sec. 5.8.6). A switch function is usually a monotonic function $S(x,\btheta )$ that takes values between 0 and 1. A switch-on function increases in x; a switch-off function decreases in x. In the log-logistic case, the function

$S(x,[\beta ,\gamma ]) = \frac{1}{1+\gamma \exp \left\{ \beta \ln (x)\right\} }$

is a switch-off function for $\beta > 0$ and a switch-on function for $\beta < 0$ . You can write general dose-response functions with asymptotes simply as

$\mr{E}[Y|x] = \mu _{\min } + (\mu _{\max }-\mu _{\min })S(x,\btheta )$

The following DATA step creates a small data set from a dose-response experiment with response y:

data logistic;
   input dose y;
   logdose = log(dose);
   datalines;
0.009  106.56
0.035   94.12
0.07    89.76
0.15    60.21
0.20    39.95
0.28    21.88
0.50     7.46
;

A graph of these data is produced with the following statements:

proc sgplot data=logistic;
   scatter y=y x=dose;
   xaxis type=log logstyle=linear;
run;

Output 69.4.1: Observed Data in Dose-Response Experiment

When dose is expressed on the log scale, the sigmoidal shape of the dose-response relationship is clearly visible (Output 69.4.1). The log-logistic switching model in the preceding parameterization is fit with the following statements in the NLIN procedure:

proc nlin data=logistic bias hougaard nlinmeasures;
   parameters alpha=100 beta=3 gamma=300;
   delta = 0;
   Switch = 1/(1+gamma*exp(beta*log(dose)));
   model y = delta + (alpha - delta)*Switch;
run;

The lower asymptote $\delta$ is assumed to be 0 in this case. Since $\delta$ is not listed in the PARAMETERS statement and is assigned a value in the program, it is assumed to be constant. Note that the term Switch is the switch-off function in the log-logistic model. The BIAS and HOUGAARD options in the PROC NLIN statement request that Box’s bias, percentage bias, and Hougaard’s skewness measure be added to the table of parameter estimates, and the NLINMEASURES option requests that the global nonlinearity measures be produced.

The NLIN procedure converges after 10 iterations and achieves a residual mean squared error of 15.1869 (Output 69.4.2). This value is not that important by itself, but it is worth noting since this model fit is compared to the fit with other parameterizations later on.

Output 69.4.2: Iteration History and Analysis of Variance

The NLIN Procedure

Dependent Variable y

Method: Gauss-Newton

Iterative Phase
Iter	alpha	beta	gamma	Sum of Squares
0	100.0	3.0000	300.0	386.4
1	100.4	2.8011	162.8	129.1
2	100.8	2.6184	101.4	69.2710
3	101.3	2.4266	69.7579	68.2167
4	101.7	2.3790	69.0358	60.8223
5	101.8	2.3621	67.3709	60.7516
6	101.8	2.3582	67.0044	60.7477
7	101.8	2.3573	66.9150	60.7475
8	101.8	2.3571	66.8948	60.7475
9	101.8	2.3570	66.8902	60.7475
10	101.8	2.3570	66.8892	60.7475

NOTE: Convergence criterion met.

Note:

An intercept was not specified for this model.

Source	DF	Sum of Squares	Mean Square	F Value	Approx Pr > F
Model	3	33965.4	11321.8	745.50	<.0001
Error	4	60.7475	15.1869
Uncorrected Total	7	34026.1

The table of parameter estimates displays the estimates of the three model parameters, their approximate standard errors, 95% confidence limits, Hougaard’s skewness measure, Box’s bias, and percentage bias (Output 69.4.3). Parameters for which the skewness measure is less than 0.1 in absolute value and with percentage bias less than 1% exhibit very close-to-linear behavior, and skewness values less than 0.25 in absolute value indicate reasonably close-to-linear behavior (Ratkowsky, 1990). According to these rules, the estimators $\widehat{\beta }$ and $\widehat{\gamma }$ suffer from substantial curvature. The estimator $\widehat{\gamma }$ is especially "far-from-linear." Inferences that involve $\widehat{\gamma }$ and rely on the reported standard errors or confidence limits (or both) for this parameter might be questionable.

Output 69.4.3: Parameter Estimates, Hougaard’s Skewness, and Box’s Bias

Parameter	Estimate	Approx Std Error	Approximate 95% Confidence Limits		Skewness	Bias	Percent Bias
alpha	101.8	3.0034	93.4751	110.2	0.1415	0.1512	0.15
beta	2.3570	0.2928	1.5440	3.1699	0.4987	0.0303	1.29
gamma	66.8892	31.6146	-20.8870	154.7	1.9200	10.9230	16.3

The related global nonlinearity measures output table (Output 69.4.4) shows that both the maximum and RMS parameter-effects curvature are substantially larger than the critical curvature value recommended by Bates and Watts (1980). In contrast, the intrinsic curvatures of the model are less than the critical value. This implies that most of the nonlinearity can be removed by reparameterization.

Output 69.4.4: Global Nonlinearity Measures

Global Nonlinearity Measures
Max Intrinsic Curvature	0.2397
RMS Intrinsic Curvature	0.1154
Max Parameter-Effects Curvature	4.0842
RMS Parameter-Effects Curvature	1.8198
Curvature Critical Value	0.3895
Raw Residual Variance	15.187
Projected Residual Variance	5.922

One method of reducing the parameter-effects curvature, and thereby reduce the bias and skewness of the parameter estimators, is to replace a parameter with its expected-value parameterization. Schabenberger et al. (1999) and Schabenberger and Pierce (2002, sec. 5.7.2) refer to this method as reparameterization through defining relationships. A defining relationship is obtained by equating the mean response at a chosen value of x (say, $x^*$ ) to the model:

$\mr{E}[Y|x^*] = \delta + \frac{\alpha - \delta }{1+\gamma \exp \left\{ \beta \ln (x^*)\right\} }$

This equation is then solved for a parameter that is subsequently replaced in the original equation. This method is particularly useful if $x^*$ has an interesting interpretation. For example, let $\lambda _ K$ denote the value that reduces the response by $K\times 100$ %,

$\mr{E}[Y|\lambda _ K] = \delta + \left(\frac{100-K}{100}\right)(\alpha -\delta )$

Because $\gamma$ exhibits large bias and skewness, it is the target in the first round of reparameterization. Setting the expression for the conditional mean at $\lambda _ K$ equal to the mean function when $x=\lambda _ K$ yields the following expression:

$\delta + \left(\frac{100-K}{100}\right)(\alpha -\delta ) = \delta + \frac{\alpha - \delta }{1+\gamma \exp \left\{ \beta \ln (\lambda _ K)\right\} }$

This expression is solved for $\gamma$ , and the result is substituted back into the model equation. This leads to a log-logistic model in which $\gamma$ is replaced by the parameter $\lambda _ K$ , the dose at which the response was reduced by $K\times 100$ %. The new model equation is

$\mr{E}[Y|x] = \delta + \frac{\alpha -\delta }{1+K/(100-K)\exp \{ \beta \ln (x/\lambda _ K)\} }$

A particularly interesting choice is K = 50, since $\lambda _{50}$ is the dose at which the response is halved. In studies of mortality, this concentration is also known as the LD50. For the special case of $\lambda _{50}$ the model equation becomes

$\mr{E}[Y|x] = \delta + \frac{\alpha -\delta }{1+\exp \{ \beta \ln (x/\lambda _{50})\} }$

You can fit the model in the LD50 parameterization with the following statements:

proc nlin data=logistic bias hougaard;
   parameters alpha=100 beta=3 LD50=0.15;
   delta = 0;
   Switch = 1/(1+exp(beta*log(dose/LD50)));
   model y = delta + (alpha - delta)*Switch;
   output out=nlinout pred=p lcl=lcl ucl=ucl;
run;

Partial results from this NLIN run are shown in Output 69.4.5. The analysis of variance tables in Output 69.4.2 and Output 69.4.5 are identical. Changing the parameterization of a model does not affect the model fit. It does, however, affect the interpretation of the parameters and the statistical properties (close-to-linear behavior) of the parameter estimators. The skewness and bias measures of the parameter LD50 is considerably reduced compared to those values for the parameter $\gamma$ in the previous parameterization. Also, $\gamma$ has been replaced by a parameter with a useful interpretation, the dose that yields a 50% reduction in mean response. Also notice that the bias and skewness measures of $\alpha$ and $\beta$ are not affected by the $\gamma \rightarrow$ LD50 reparameterization.

Output 69.4.5: ANOVA Table and Parameter Estimates in LD50 Parameterization

The NLIN Procedure

Dependent Variable y

Method: Gauss-Newton

NOTE: Convergence criterion met.

Note:

An intercept was not specified for this model.

Source	DF	Sum of Squares	Mean Square	F Value	Approx Pr > F
Model	3	33965.4	11321.8	745.50	<.0001
Error	4	60.7475	15.1869
Uncorrected Total	7	34026.1

Parameter	Estimate	Approx Std Error	Approximate 95% Confidence Limits		Skewness	Bias	Percent Bias
alpha	101.8	3.0034	93.4752	110.2	0.1415	0.1512	0.15
beta	2.3570	0.2928	1.5440	3.1699	0.4987	0.0303	1.29
LD50	0.1681	0.00915	0.1427	0.1935	-0.0605	-0.00013	-0.08

To reduce the parameter-effects curvature of the $\beta$ parameter, you can use the technique of defining relationships again. This can be done generically, by solving

$\mu ^* = \delta + \frac{\alpha -\delta }{1+\exp \left\{ \beta \ln (x/\lambda )\right\} }$

for $\beta$ , treating $\mu ^*$ as the new parameter (in lieu of $\beta$ ), and choosing a value for $x^*$ that leads to low skewness. This results in the expected-value parameterization of $\beta$ . Solving for $\beta$ yields

$\beta = \frac{\log \left( \frac{\alpha -\mu ^*}{\mu ^*-\delta }\right)}{\log \left( x^* / \lambda \right)}$

The interpretation of the parameter $\mu ^*$ that replaces $\beta$ in the model equation is simple: it is the mean dose response when the dose is $x^*$ . Fixing $x^*=0.3$ , the following PROC NLIN statements fit this model:

proc nlin data=logistic bias hougaard nlinmeasures;
   parameters alpha=100 mustar=20 LD50=0.15;
   delta   = 0;
   xstar   = 0.3;
   beta    = log((alpha - mustar)/(mustar - delta)) / log(xstar/LD50);
   Switch  = 1/(1+exp(beta*log(dose/LD50)));
   model y = delta + (alpha - delta)*Switch;
   output out=nlinout pred=p lcl=lcl ucl=ucl;
run;

Note that the switch-off function continues to be written in terms of $\beta$ and the LD50. The only difference from the previous model is that $\beta$ is now expressed as a function of the parameter $\mu ^*$ . Using expected-value parameterizations is a simple mechanism to lower the curvature in a model and to arrive at starting values. The starting value for $\mu ^*$ can be gleaned from Output 69.4.1 at x = 0.3.

Output 69.4.6 shows selected results from this NLIN run. The ANOVA table is again unaffected by the change in parameterization. The skewness for $\mu ^*$ is significantly reduced in comparison to those of the $\beta$ parameter in the previous model (compare Output 69.4.6 and Output 69.4.5), while its bias remains on the same scale from Output 69.4.5 to Output 69.4.6. Also note the substantial reduction in the parameter-effects curvature values. As expected, the intrinsic curvature values remain intact.

Output 69.4.6: ANOVA Table and Parameter Estimates in Expected-Value Parameterization

The NLIN Procedure

Dependent Variable y

Method: Gauss-Newton

NOTE: Convergence criterion met.

Note:

An intercept was not specified for this model.

Source	DF	Sum of Squares	Mean Square	F Value	Approx Pr > F
Model	3	33965.4	11321.8	745.50	<.0001
Error	4	60.7475	15.1869
Uncorrected Total	7	34026.1

Parameter	Estimate	Approx Std Error	Approximate 95% Confidence Limits		Skewness	Bias	Percent Bias
alpha	101.8	3.0034	93.4752	110.2	0.1415	0.1512	0.15
mustar	20.7073	2.6430	13.3693	28.0454	-0.0572	-0.0983	-0.47
LD50	0.1681	0.00915	0.1427	0.1935	-0.0605	-0.00013	-0.08

Global Nonlinearity Measures
Max Intrinsic Curvature	0.2397
RMS Intrinsic Curvature	0.1154
Max Parameter-Effects Curvature	0.2925
RMS Parameter-Effects Curvature	0.1500
Curvature Critical Value	0.3895
Raw Residual Variance	15.187
Projected Residual Variance	5.9219