The NLIN Procedure

Measures of Nonlinearity, Diagnostics and Inference

Subsections:

A "close-to-linear" nonlinear regression model, in the sense of Ratkowsky (1983, 1990), is a model in which parameter estimators have properties similar to those in a linear regression model. That is, the least squares estimators of the parameters are close to being unbiased and normally distributed, and they have minimum variance.

A nonlinear regression model sometimes fails to be close to linear due to the properties of one or several parameters. When this occurs, bias in the parameter estimates can render inferences that use the reported standard errors and confidence limits invalid.

PROC NLIN provides various measures of nonlinearity. To assess the nonlinearity of a model-data combination, you can use both of the following complementary sets of measures:

  • Box’s bias (Box, 1971) and Hougaard’s skewness (Hougaard, 1982, 1985) of the least squares parameter estimates

  • curvature measures of nonlinearity (Bates and Watts, 1980).

Furthermore, PROC NLIN provides residual, leverage, and local-influence diagnostics (St. Laurent and Cook, 1993).

In the following several sections, these nonlinearity measures and diagnostics are discussed. For this material, several basic definitions are required. Let $\mb{X}$ be the Jacobian matrix for the model, $\mb{X} = \frac{\partial \mb{f}}{\partial \bbeta }$, and let $\mb{Q}$ and $\mb{R}$ be the components of the QR decomposition of $\mb{X} = \mb{Q}\mb{R}$ of $\mb{X}$, where $\mb{Q}$ is an $(n\times n)$ orthogonal matrix. Finally, let $\mb{B}$ be the inverse of the matrix constructed from the first p rows of the $(n \times p)$ dimensional matrix $\mb{R}$ (that is, $\mb{B} = \mb{R}^{-1}_{p}$). Next define

\begin{align*} \left[\mb{H}_ j\right]_{kl} & = \frac{\partial ^2 \mb{f}_ j}{\partial \bbeta _ k \partial \bbeta _ l}\\ \left[\mb{U}_ j\right]_{kl} & = \sum _{mn} \mb{B}’_{km} \left[\mb{H}_ j\right]_{mn} \mb{B}_{nl} \\ \left[\mb{A}_ j\right]_{kl} & = \sqrt {p \times \mr{mse} }\sum _{m} \mb{Q}’_{jm} \left[\mb{U}_ m\right]_{kl}\, , \end{align*}

where $\mb{H}$, $\mb{U}$ and the acceleration array $\mb{A}$ are three-dimensional $(n \times p \times p)$ matrices. The first p faces of the acceleration array constitute a $(p \times p \times p)$ parameter-effects array and the last $(n-p)$ faces constitute the $(n-p \times p \times p)$ intrinsic curvature array (Bates and Watts, 1980). The previous and subsequent quantities are computed at the least squares parameter estimators.

Box’s Measure of Bias

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

\begin{align*} \widehat{\mr{E}}\left[ \widehat{\beta } - \beta \right] & = -\frac{\sigma ^2}{2} \left(\mb{X}’ \bW \mb{X}\right)^{-1} \sum ^{n}_{i=1} w_ i \mb{x}’_ i \, \text {Tr}\left( \left(\mb{X}’ \bW \mb{X}\right)^{-1} \left[\mb{H}_{i}\right]\right) \end{align*}

where $\sigma ^2 = \hbox{mse}$ if the SIGSQ option is not set. Otherwise, $\sigma ^2$ is the value you set with the SIGSQ option. $\mb{W}$ is the diagonal weight matrix specified with the _WEIGHT_ variable (or the identity matrix if _WEIGHT_ is not defined) and $\left[\mb{H}_{i}\right]$ is the $(p \times p)$ 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 $\mb{A}$,

\begin{align*} \widehat{\mr{E}}\left[ \widehat{\beta }_ i - {\beta }_ i\right] & = - \frac{\sigma ^2}{2 p \times \hbox{mse}} \sum ^{p}_{j,\, k=1} \mb{B}_{ij} \left[\mb{A}_ j\right]_{kk} \end{align*}

As the preceding formulas illustrate, the bias depends solely on the parameter-effects array, thereby permitting its reduction through reparameterization. Example 69.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%.

Hougaard’s Measure of Skewness

In addition to Box’s bias, Hougaard’s measure of skewness, $g_{1i}$ (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

\[ \mr{E} \left[ \widehat{\beta }_ i - \mr{E}\left(\widehat{\beta }_ i\right)\right]^3 = - \left(\sigma ^2\right)^2 \sum _{jkl} [\mb{L}]_{ij}[\mb{L}]_{ik}[\mb{L}]_{il} \left( \left[\mb{V}_ j\right]_{kl}+[\mb{V}_ k]_{jl}+[\mb{V}_ l]_{jk}\right) \]

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

\[ \mb{L} = \left(\mb{X}’\mb{X}\right)^{-1} = \left( \frac{\partial \mb{f}}{\partial \bbeta ^\prime } \frac{\partial \mb{f}}{\partial \bbeta } \right)^{-1} \]

The term $[\mb{L}]_{ij}$ denotes the value in row i, column j of the matrix $\mb{L}$. (Hougaard (1985) uses superscript notation to denote elements in this inverse.) The matrix $\mb{V}$ is a three-dimensional $(p \times p \times p)$ array

\begin{align*} \left[\mb{V}_ j\right]_{kl} & = \sum _{m=1}^ n \frac{\partial F_ m}{\partial \beta _ j} \frac{\partial ^2 F_ m}{\partial \beta _ k \partial \beta _ l} \end{align*}

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

\[ G_{1i} = \mr{E}\left[ \widehat{\beta }_ i - \mr{E}(\widehat{\beta }_ i)\right]^3 / \left(\sigma ^2 \times \left[ \mb{L} \right]_{ii} \right)^{3/2} \]

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

\begin{align*} g_{1i} & = \widehat{\mr{E}}\left[ \widehat{\beta }_ i - \mr{E}(\widehat{\beta }_ i)\right]^3 /\left( \mr{mse} \times [\widehat{\mb{L}}]_{ii} \right) ^{3/2} \\ \widehat{\mr{E}}\left[ \widehat{\beta }_ i - \mr{E}(\widehat{\beta }_ i)\right]^3 & = - \mr{mse}^2 \sum _{jkl} [\widehat{\mb{L}}]_{ij} [\widehat{\mb{L}}]_{ik} [\widehat{\mb{L}}]_{il} \left( [\widehat{\mb{V}}_ j]_{kl} + [\widehat{\mb{V}}_ k]_{jl} + [\widehat{\mb{V}}_ l]_{jk} \right) \end{align*}

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

Relative Curvature Measures of Nonlinearity

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

\begin{align*} \mb{T}^{\tau } & = \text {max}_{\beta } \big {|}\big {|}\beta ^{'} \mb{A}^{\tau } \beta \big {|}\big {|}\\ \mb{T}^{\eta } & = \text {max}_{\beta } \big {|}\big {|}\beta ^{'} \mb{A}^{\eta } \beta \big {|}\big {|} \end{align*}

where $\mb{T}^{\tau }$ and $\mb{T}^{\eta }$ denote the maximum parameter-effects and intrinsic curvatures, while $\mb{A}^{\tau }$ and $\mb{A}^{\eta }$ 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 $10^{-4}$ 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, $\beta $. 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 $\mb{T}^{\tau }$ and $\mb{T}^{\eta }$ and the corresponding RMS values can be assessed by comparing these values with $1/\sqrt {F}$, where F is the upper $\alpha \times 100\% $ quantile of an F distribution with p and $n-p$ 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 69.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 69.6 shows a nonlinear model with a high intrinsic curvature and the corresponding diagnostics.

Leverage in Nonlinear Regression

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, $\mb{H}_{i}$, 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, $\mb{J}_{i}$, 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

\begin{align*} \mb{H}_ i & = w_ i \mb{x}_ i(\mb{X}^{\prime }\mb{WX})^{-1} \mb{x}_ i^{\prime }\\ \mb{J}_ i & = w_ i \mb{x}_ i (\mb{X}^{\prime }\mb{WX}-[\mb{W\, e}][\mb{H}])^{-1} \mb{x}_ i^{\prime }\, , \end{align*}

where $\mb{e}$ is the vector of residuals, $\mb{W}$ is the diagonal weight matrix if you specify the special variable _WEIGHT_ and otherwise the identity matrix, and i indexes the corresponding quantities for the ith observation. 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 ($\mb{e} = 0$) 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 69.6.

Local Influence in Nonlinear Regression

St. Laurent and Cook (1993) suggest using $l_{\hbox{max}}$, 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

\begin{align*} \mb{C}_\beta & = \hbox{max}_ l \, \frac{2}{\sigma ^2} l’ \mb{J} l\\ \mb{C}_\sigma & = \hbox{max}_ l\, \frac{4}{\sigma ^2} l’ \mb{P}_ e l \end{align*}

where $\mb{J}$ is the Jacobian leverage matrix and $\mb{P}_ e = \mb{e}\mb{e}’/(\mb{e}’\mb{e})$, you choose the $l_{\hbox{max}}$ that results in the maximum of the two curvature components (St. Laurent and Cook, 1993). PROC NLIN provides $l_{\hbox{max}}$ through the PLOTS option in the PROC NLIN statement and the LMAX= option in the OUTPUT statement. Example 69.6 shows a plot of $l_{\hbox{max}}$ for a model with high intrinsic curvature.

Residuals in Nonlinear Regression

If a nonlinear model is intrinsically nonlinear, using the residuals $\mb{e}=\mb{y} - \mb{\hat{y}}$ 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 $\mb{e}$ (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 $\mb{e}_{p}$ and the expectation of the raw residuals by $\mr{E}[\mb{e}]$. Then

\begin{align*} \mb{e}_{p} & = \left(\mb{I}_ n - \mb{P}_{xh}\right)\mb{e}\, \\ \mr{E} \left[\mb{e}_ i\right] & = -\frac{\sigma ^2}{2}\, \sum ^{n}_{j=1} \tilde{\mb{P}}_{x,\, ij} \text {Tr}\left(\, \left[\mb{H}_ j\right] \, \left(\mb{X}’\mb{X}\right)^{-1}\, \right) \end{align*}

where $\mb{e}_ i$ is the ith observation raw residual, $\mb{I}_ n$ is an n-dimensional identity matrix, $\mb{P}_{xh}$ is the projector onto the column space of $\left(\mb{X}|\mb{H}\right)$, and $\tilde{\mb{P}}_ x = \mb{I}_ n - \mb{P}_ x$. The preceding formulas are general with the projectors defined accordingly to take the weighting into consideration. In unweighted least squares, $\mr{E} \left[\mb{e}\right]$ reduces to

\begin{align*} \mr{E} \left[\mb{e}\right] & = -\frac{1}{2} \sigma ^2 \tilde{\mb{Q}}\, \mb{a}\, \end{align*}

with $\tilde{\mb{Q}}$ being the last $n-p$ columns of the $\mb{Q}$ matrix in the QR decomposition of $\mb{X}$ and the $(n-p)$ dimensional vector $\mb{a}$ being defined in terms of the intrinsic acceleration array

\begin{align*} \mb{a}_ i & = \sum ^{p}_{j=1}\left[\mb{A}_{i+p}\right]_{jj} \end{align*}

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

\begin{align*} \sigma ^2_{p} & = \frac{\mb{e}'_{p} \mb{e}_{p}}{\text {Tr}\left(\mb{I}_ n - \mb{P}_{xh}\right)} \end{align*}

The standardized raw and projected residuals, denoted by $\tilde{\mb{e}}$ and $\tilde{\mb{e}}_{p}$ respectively, are given by

\begin{align*} \tilde{\mb{e}} & = \frac{\sqrt {w_ i}\, \mb{e}}{ \sigma \sqrt {1 - \mb{P}_{x,ii}}}\\ \tilde{\mb{e}}_{p} & = \frac{\sqrt {w_ i} \, \mb{e}_ p}{ \sigma \sqrt {1 - \mb{P}_{xh,ii}}} \end{align*}

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

Profiling Parameters and Assessing the Influence of Observations on Parameter Estimates

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 $\beta $ is a plot of the likelihood ratio pivotal statistic, $\tau (\beta )$, and the corresponding Wald pivotal statistic, $\sigma (\beta )$ (Bates and Watts, 1988). The likelihood ratio pivotal statistic is defined as

\[ \tau (\beta ) = \mbox{sign}(\beta - \hat{\beta })\, \mbox{L}(\beta ) \]

with

\[ \mbox{L}(\beta ) = \left( \frac{\mbox{SSE}(\beta , \tilde{\Theta }) - \mbox{SSE}(\hat{\beta }, \hat{\Theta }) }{\mbox{mse}} \right)^{1/2} \]

where $\beta $ is the profile parameter and $\Theta $ refers to the remaining parameters. $\mr{SSE}(\beta , \tilde{\Theta })$ is the sum of square errors where the profile parameter $\beta $ is constrained at a given value and $\tilde{\Theta }$ is the least squares estimate of $\Theta $ conditional on a given value of $\beta $. In contrast, $\mr{SSE}(\hat{\beta }, \hat{\Theta })$ is the sum of square errors for the full model. For linear models, $\tau (\beta )$ matches $\sigma (\beta )$, which is defined as

\[ \sigma (\beta ) = \frac{\beta - \hat{\beta }}{\mbox{stderr}_\beta } \]

where $\beta $ is the constrained value, $\hat{\beta }$ is the estimated value for the parameter, and $\mbox{stderr}_\beta $ 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 $\mr{L}(\beta )$. The Wald-based confidence intervals are overlaid as two straight lines that pass through $(0, \hat{\beta })$ with a slope of $\pm \mbox{stderr}_\beta $. 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

Bootstrap resampling and estimation methods can be used to produce confidence intervals and covariance matrix estimates that have second-order, $\mathcal{O}(1/n)$ accuracy, where n is the number of observations (DiCiccio and Efron, 1996). In contrast, the standard Wald-based confidence interval has first-order, $\mathcal{O}(1/{\sqrt {n}})$ 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

\[ \tilde{Y}_ i = f(\hat{\bbeta }; \mb{z}_{i}^\prime ) + \tilde{\epsilon }_ i \]

where $\tilde{Y}_ i$ is the ith simulated response, $\hat{\bbeta }$ is the original least squares parameter estimate, $\mb{z}_ i^\prime $ is the ith regressor vector, and $\tilde{\epsilon }_ i$ 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, $\tilde{\epsilon }_ i$. In general, $\tilde{\epsilon }_ i$ can be represented as

\[ \tilde{\epsilon }_ i = {\mb{s}}_ r {\hat{\mb{e}}}_{r} \]

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

\[ \mb{s}_ r = \left\{ \begin{array}{ll} 1 & \quad \text {if DGP=RESIDUAL(RAW)} \cr \sqrt {\frac{\mb{n}}{\mb{n} - \mb{p}}} & \quad \text {if DGP=RESIDUAL(ADJSSE)} \cr \frac{1}{\sqrt {1 - \mb{H}_ r}} & \quad \text {if DGP=RESIDUAL(TAN)} \cr \frac{1}{\sqrt {1 - \mb{J}_ r}} & \quad \text {if DGP=RESIDUAL(JAC)} \cr \frac{\gamma }{\sqrt {1 - \mb{H}_ r}} & \quad \text {if DGP=WILD} \end{array} \right. \]

In the preceding formula, $\mb{n}$ is the number of usable observations, $\mb{p}$ is the number of model parameters, and $\mb{H}_ r$ and $\mb{J}_ r$ 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, $\gamma $ is a random number given by

\[ \gamma = \left\{ \begin{array}{ll} -\frac{\sqrt {5} - 1}{2} & \, \, \, \mr{with}\, \, \mr{probability} \, \, \, \, \, \frac{\sqrt {5} + 1}{2\sqrt {5}} \cr \frac{\sqrt {5} + 1}{2} & \, \, \, \mr{with}\, \, \, \mr{ probability} \, \, \, \, \, \frac{\sqrt {5} - 1}{2\sqrt {5}} \end{array} \right. \]

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 $\beta _ i$, the ith model parameter, follow. For simplicity of notation, denote $\beta _ i$ as $\beta $. 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 $100 \alpha /2$ and $100 (1-\alpha /2)$ percentiles from the bootstrap parameter estimates. $\alpha $ is from the ALPHA= option in the PROC NLIN statement. These percentiles are computed as follows. Let $\tilde{\beta }_{1}, \tilde{\beta }_{2},$ …$,\tilde{\beta }_{B}$ represent the ordered values of the bootstrap estimates for $\beta $. Let the kth weighted average percentile be y, set $p=\frac{k}{100}$, and let

\[ n p = j + g \]

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

\[ y = \left\{ \begin{array}{ll} \frac{1}{2} \big ( \tilde{\beta }_{j} + \tilde{\beta }_{j+1}\big ) & \, \, \, \mr{if} \, \, \, \, \, g = 0 \cr \tilde{\beta }_{j+1} & \, \, \, \mr{if} \, \, \, \, \, g > 0 \end{array} \right. \]

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 $\alpha $ level, is given by

\[ \bar{\tilde{\beta }} \pm \mbox{stdb}_{\tilde{\beta }} \times \mr{z}^{(1-\alpha /2)} \]

where $\bar{\tilde{\beta }}$ is the mean of the bootstrap estimates, $\tilde{\beta }$; $\mbox{stdb}_{\tilde\beta }$ is the standard deviation of the bootstrap parameter estimates; and $\mr{z}^{(1-\alpha /2)}$ is the $100 (1-\alpha /2)$ 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), $G(\tilde{\beta })$, of the bootstrap parameter estimates to correct for the upper and lower endpoints of the $\alpha $ level. The bias-corrected bootstrap confidence interval is given by

\[ G^{-1} \big (\Phi (2 z_0 {\pm } z_{\alpha })\big ) \]

where $\Phi $ is the standard normal CDF, $z_\alpha = \Phi ^{-1} (\alpha )$, and $z_0$ is a bias correction given by

\[ z_0 = \Phi ^{-1}\bigg({\frac{\mr{N}\big (\tilde{\beta } \leq \hat{\beta }\big )}{\mr{B}}}\bigg) \]

where $\hat{\beta }$ is the estimate of $\beta $ from the original least squares fit and $\mr{N}(\tilde{\beta } \leq \hat{\beta })$ is the number of bootstrap estimates, $\tilde{\beta }$, that are less than or equal to $\hat{\beta }$.

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

\[ \mbox{Cov}(\tilde{\beta }_ i, \tilde{\beta }_ j) = \sum ^{k=\mr{B}}_{k=1}\frac{(\tilde{\beta }_{ik} - \bar{\tilde{\beta }}_{i}) (\tilde{\beta }_{jk} - \bar{\tilde{\beta }}_{j})}{\mr{B} - 1} \]

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

\[ \mbox{Corr}(\tilde{\beta }_ i, \tilde{\beta }_ j) = \frac{\mbox{Cov}(\tilde{\beta }_ i, \tilde{\beta }_ j)}{\mr{stdb}_{\tilde{\beta }_ i} \mr{stdb}_{\tilde{\beta }_ j}} \]

where $\mbox{stdb}_{\tilde{\beta }_ i}$ and $\mbox{stdb}_{\tilde{\beta }_ j}$ are the bootstrap standard deviations for the ith and jth parameters.