The NLIN Procedure

The degree to which parameter estimators exhibit close-to-linear behavior can be assessed with Box’s bias (Box 1971) and Hougaard’s measure of skewness (Hougaard 1982, 1985). The bias and percentage bias measures are available through the BIAS option in the PROC NLIN statement. Box’s bias measure is defined as

where if the SIGSQ option is not set. Otherwise, is the value you set with the SIGSQ option. is the diagonal weight matrix specified with the _WEIGHT_ variable (or the identity matrix if _WEIGHT_ is not defined) and is the Hessian matrix at the ith observation. In the case of unweighted least squares, the bias formula can be expressed in terms of the acceleration array ,

As the preceding formulas illustrate, the bias depends solely on the parameter-effects array, thereby permitting its reduction through reparameterization. Example 88.4 shows how changing the parameterization of a four-parameter logistic model can reduce the bias. Ratkowsky (1983, p. 21) recommends that you consider reparameterization if the percentage bias exceeds 1%.

In addition to Box’s bias, Hougaard’s measure of skewness, (Hougaard 1982, 1985), is also provided in PROC NLIN to assess the close-to-linear behavior of parameter estimators. This measure is available through the HOUGAARD option in the PROC NLIN statement. Hougaard’s skewness measure for the ith parameter is based on the third central moment, defined as

where the sum is a triple sum over the number of parameters and

The term denotes the value in row i, column j of the matrix . (Hougaard (1985) uses superscript notation to denote elements in this inverse.) The matrix is a three-dimensional array

The third central moment is then normalized using the standard error as

The previous expressions depend on the unknown values of the parameters and on the residual variance . In order to evaluate the Hougaard measure in a particular data set, the NLIN procedure computes

Following Ratkowsky (1990, p. 28), the parameter is considered to be very close to linear, reasonably close, skewed, or quite nonlinear according to the absolute value of the Hougaard measure being less than 0.1, between 0.1 and 0.25, between 0.25 and 1, or greater than 1, respectively.

Bates and Watts (1980) formulated the maximum parameter-effects and maximum intrinsic curvature measures of nonlinearity to assess the close-to-linear behavior of nonlinear models. Ratkowsky (1990) notes that of the two curvature components in a nonlinear model, the parameter-effects curvature is typically larger. It is this component that you can affect by changing the parameterization of a model. PROC NLIN provides these two measures of curvature both through the STATS plot-option and through the NLINMEASURES option in the PROC NLIN statement.

The maximum parameter-effects and intrinsic curvatures are defined, in a compact form, as

where and denote the maximum parameter-effects and intrinsic curvatures, while and stand for the parameter-effects and intrinsic curvature arrays. The maximization is carried out over a unit-vector of the parameter values (Bates and Watts 1980). In line with Bates and Watts (1980), PROC NLIN takes as the convergence tolerance for the maximum intrinsic and parameter-effects curvatures. Note that the preceding matrix products involve contraction of the faces of the three-dimensional acceleration arrays with the normalized parameter vector, . The corresponding expressions for the RMS (root mean square) parameter-effects and intrinsic curvatures can be found in Bates and Watts (1980).

The statistical significance of and and the corresponding RMS values can be assessed by comparing these values with , where F is the upper quantile of an F distribution with p and degrees of freedom (Bates and Watts 1980).

One motivation for fitting a nonlinear model in a different parameterization is to obtain a particular interpretation and to give parameter estimators more close-to-linear behavior. Example 88.4 shows how changing the parameterization of a four-parameter logistic model can reduce the parameter-effects curvature and can yield a useful parameter interpretation at the same time. In addition, Example 88.6 shows a nonlinear model with a high intrinsic curvature and the corresponding diagnostics.

In contrast to linear regression, there are several measures of leverage in nonlinear regression. Furthermore, in nonlinear regression, the effect of a change in the ith response on the ith predicted value might depend on both the size of the change and the ith response itself (St. Laurent and Cook 1992). As a result, some observations might show superleverage —namely, leverages in excess of one (St. Laurent and Cook 1992).

PROC NLIN provides two measures of leverages: tangential and Jacobian leverages through the PLOTS option in the PROC NLIN statement and the H= and J= options of OUTPUT statement. Tangential leverage, , is based on approximating the nonlinear model with a linear model that parameterizes the tangent plane at the least squares parameter estimators. In contrast, Jacobian leverage, , is simply defined as the instantaneous rate of change in the ith predicted value with respect to the ith response (St. Laurent and Cook 1992).

The mathematical formulas for tangential and Jacobian leverages are

where is the vector of residuals, is the diagonal weight matrix if you specify the special variable _WEIGHT_ and otherwise the identity matrix, i indexes the corresponding quantities for the ith observation, and the matrix is defined in the section Notation for Nonlinear Regression Models. The brackets indicate column multiplication as defined in Bates and Watts (1980). The preceding formula for tangential leverage holds if the gradient, Marquardt, or Gauss methods are used. For the Newton method, the tangential leverage is set equal to the Jacobian leverage.

In a model with a large intrinsic curvature, the Jacobian and tangential leverages can be very different. In fact, the two leverages are identical only if the model provides an exact fit to the data () or the model is intrinsically linear (St. Laurent and Cook 1993). This is also illustrated by the leverage plot and nonlinearity measures provided in Example 88.6.

St. Laurent and Cook (1993) suggest using , the direction that yields the maximum normal curvature, to assess the local influence of an additive perturbation to the response variable on the estimation of the parameters and variance of a nonlinear model. Defining the normal curvature components

where is the Jacobian leverage matrix and , you choose the that results in the maximum of the two curvature components (St. Laurent and Cook 1993). PROC NLIN provides through the PLOTS option in the PROC NLIN statement and the LMAX= option in the OUTPUT statement. Example 88.6 shows a plot of for a model with high intrinsic curvature.

If a nonlinear model is intrinsically nonlinear, using the residuals for diagnostics can be misleading (Cook and Tsai 1985). This is due to the fact that in correctly specified intrinsically nonlinear models, the residuals have nonzero means and different variances, even when the original error terms have identical variances. Furthermore, the covariance between the residuals and the predicted values tends to be negative semidefinite, complicating the interpretation of plots based on (Cook and Tsai 1985).

Projected residuals are proposed by Cook and Tsai (1985) to overcome these shortcomings of residuals, which are henceforth called raw (ordinary) residuals to differentiate them from their projected counterparts. Projected residuals have zero means and are uncorrelated with the predicted values. In fact, projected residuals are identical to the raw residuals in intrinsically linear models.

PROC NLIN provides raw and projected residuals, along with their standardized forms. In addition, the mean or expectation of the raw residuals is available. These can be accessed with the PLOTS option in the PROC NLIN statement and the OUTPUT statement options PROJRES=, PROJSTUDENT=, RESEXPEC=, RESIDUAL= and STUDENT=.

Denote the projected residuals by and the expectation of the raw residuals by . Then

where is the ith observation raw residual, is an n-dimensional identity matrix, is the projector onto the column space of , and . The preceding formulas are general with the projectors defined accordingly to take the weighting into consideration. In unweighted least squares, reduces to

with being the last columns of the matrix in the QR decomposition of and the dimensional vector being defined in terms of the intrinsic acceleration array

Standardization of the projected residuals requires the variance of the projected residuals. This is estimated using the formula (Cook and Tsai 1985)

The standardized raw and projected residuals, denoted by and respectively, are given by

The use of raw and projected residuals for diagnostics in nonlinear regression is illustrated in Example 88.6.

The global measures of nonlinearity, discussed in the preceding section, are very useful for assessing the overall nonlinearity of the model. However, the impact of global nonlinearity on inference regarding subsets of the parameter set cannot be easily determined (Cook and Tsai 1985). The impact of the nonlinearity on the uncertainty of individual parameters can be efficiently described by profile t plots and confidence curves (Bates and Watts 1988; Cook and Weisberg 1990).

A profile t plot for parameter is a plot of the likelihood ratio pivotal statistic, , and the corresponding Wald pivotal statistic, (Bates and Watts 1988). The likelihood ratio pivotal statistic is defined as

with

where is the profile parameter and refers to the remaining parameters. is the sum of square errors where the profile parameter is constrained at a given value and is the least squares estimate of conditional on a given value of . In contrast, is the sum of square errors for the full model. For linear models, matches , which is defined as

where is the constrained value, is the estimated value for the parameter, and is the standard error for the parameter. Usually a profile t plot is overlaid with a reference line that passes through the origin and has a slope of one. PROC NLIN follows this convention.

A confidence curve for a particular parameter is useful for validating Wald-based confidence intervals for the parameter against likelihood-based confidence intervals (Cook and Weisberg 1990). A confidence curve contains a scatter plot of the constrained parameter value versus . The Wald-based confidence intervals are overlaid as two straight lines that pass through with a slope of . Hence, for different levels of significance, you can easily compare the Wald-based confidence intervals against the corresponding confidence intervals that are based on the likelihood ratio. Cook and Weisberg (1990) recommend that you report a single Wald-based confidence interval only if there is a good agreement between the Wald and likelihood intervals up to at least the 99% confidence level.

Compared to local influence, the leave-one-out method performs a more complete analysis of the influence of observations on the values of parameter estimates. In this method, jackknife resampling removes each observation in turn and fits the model to the remaining data set. Hence, a data set with n observations will have n corresponding data sets with n–1 observations. The impact of each observation on a parameter is assessed by the absolute relative percentage change in the value of the parameter compared with the reference value from the full data.

Bootstrap resampling and estimation methods can be used to produce confidence intervals and covariance matrix estimates that have second-order, accuracy, where n is the number of observations (DiCiccio and Efron 1996). In contrast, the standard Wald-based confidence interval has first-order, accuracy. Bootstrap methods achieve this higher accuracy at the cost of an orders-of-magnitude increase in numerical computation compared to standard asymptotic approximations. However, dramatic increases in performance and decreases in the cost of numerical computation have made bootstrap methods very attractive for routine statistical inference (MacKinnon 2002).

PROC NLIN samples as many bootstrap sample data sets (replicates) as you specify in the NSAMPLES= option and performs least squares fit on each replicate. For each least squares fit, PROC NLIN uses the original parameter estimates as starting values for the model parameters. The statistics from the least squares fits that converge are collected and used to produce confidence intervals, covariance and correlation matrices, and histograms and scatter plots of the bootstrap parameter estimates.

Each replicate is obtained by sampling the residuals instead of the input data set. The sampled residuals are used to simulate the response by using a bootstrap data generating process (DGP). PROC NLIN’s bootstrap DGP produces replicates that contain the same number of observations as the number of usable observations in the input data set. In fact, the bootstrap DGP for a particular replicate starts by discarding the input data set observations that PROC NLIN deems unusable during the original least squares estimation. The next step of the bootstrap DGP for a replicate can be represented by the formula

where is the ith simulated response, is the original least squares parameter estimate, is the ith regressor vector, and is the ith simulated error.

PROC NLIN makes several bootstrap DGP types available in the DGP= option of the BOOTSTRAP statement. These bootstrap DGP types differ only in how they obtain the ith simulated error, . In general, can be represented as

where r is a uniform random integer between 1 and the number of usable observations, is a scale factor that depends on the chosen bootstrap DGP options, and is the rth residual obtained from the rth usable observation of the original least squares fit. The scale factor, , that captures the differences among the various bootstrap DGP types takes one of the following values:

In the preceding formula, is the number of usable observations, is the number of model parameters, and and are the rth tangential and Jacobian leverages, respectively. For the WILD bootstrap DGP (Wu 1986), which is the only bootstrap DGP type exclusively available for weighted least squares, is a random number given by

PROC NLIN makes three types of bootstrap confidence intervals available in the BOOTCI option in the BOOTSTRAP statement. These confidence intervals are the percentile, normal, and bias-corrected bootstrap confidence intervals. The computational details of these confidence intervals for , the ith model parameter, follow. For simplicity of notation, denote as . Also, let B represent the number of replicates for which the least squares fit converges

The option that computes the percentile bootstrap confidence interval, BOOTCI(PERC), does so by computing the th and th percentiles from the bootstrap parameter estimates. is from the ALPHA= option in the PROC NLIN statement. These percentiles are computed as follows. Let … represent the ordered values of the bootstrap estimates for . Let the kth weighted average percentile be y, set , and let

where j is the integer part of and g is the fractional part of . Then the kth percentile, y, is given by

which corresponds to the default percentile definition of the UNIVARIATE procedure.

In contrast, the BOOTCI(NORMAL) option in the BOOTSTRAP statement computes the normal bootstrap confidence interval by approximating the distribution of the bootstrap parameter estimates as a normal distribution. Consequently, the normal bootstrap confidence interval, for level, is given by

where is the original least squares parameter estimate, is the standard deviation of the bootstrap parameter estimates, ; and is the th percentile of the standard normal distribution.

The BOOTSTRAP statement option that computes the bias-corrected bootstrap confidence interval, BOOTCI(BC), does so by making use of the cumulative distribution function (CDF), , of the bootstrap parameter estimates to correct for the upper and lower endpoints of the level. The bias-corrected bootstrap confidence interval is given by

where is the standard normal CDF, , and is a bias correction given by

where is the estimate of from the original least squares fit and is the number of bootstrap estimates, , that are less than or equal to .

In addition, PROC NLIN produces bootstrap estimates of the covariance and correlation matrices. For the ith and jth model parameters, the covariance is estimated by

where the sum runs over the nonmissing bootstrap parameter estimates. The bootstrap correlation matrix is estimated by scaling the bootstrap covariance matrix

where and are the bootstrap standard deviations for the ith and jth parameters.