Bootstrap Confidence Regions for Optimal Operating Conditions in Response Surface Methodology
aa r X i v : . [ s t a t . M E ] N ov Bootstrap Confidence Regions forOptimal Operating Conditions inResponse Surface Methodology R oger D. Gibb ∗ , I- Li Lu † , and W alter H. Carter, Jr ‡ Abstract
This article concerns the application of bootstrap methodology to constructa likelihood-based confidence region for operating conditions associated withthe maximum of a response surface constrained to a specified region. Unlikeclassical methods based on the stationary point, proper interpretation of thisconfidence region does not depend on unknown model parameters. In addition,the methodology does not require the assumption of normally distributed er-rors. The approach is demonstrated for concave-down and saddle system casesin two dimensions. Simulation studies were performed to assess the coverageprobability of these regions.AMS 2000 subj classification: 62F25, 62F40, 62F30, 62J05.
Key words : Stationary point; Kernel density estimator; Boundary kernel. ∗ [email protected], The Procter & Gamble Company, Mason, Ohio 45040-9452. † Correspondence: [email protected], Applied Statistics, Phantom Works, The Boeing Com-pany, P.O. Box 3707, MC 7L-22, Seattle, WA 98124-2207. ‡ [email protected], Department of Biostatistics, Virginia Commonwealth University, Richmond,VA 23298-0032. One of the reasons for using a response surface analysis is to determine the operatingconditions which, without loss of generality, maximize a response y . Developmentof confidence regions for these operating conditions has been considered by severalinvestigators. The purpose of this report is to demonstrate that likelihood-basedbootstrap confidence region methodology is a powerful alternative which does notsuffer from some of the limitations associated with existing approaches.Over the experimental region it is assumed that the relationship of y and theregressors x , x , . . . , x k can be expressed y = g ( x , θ ) + ε, (1)where g is an unknown continuous, differentiable function and ε is a random sourceof variability not accounted for in g . Often g ( x , θ ) can be adequately approximatedwith a second-order polynomial, i.e. g ( x ) ≈ β + xβ + x ′ Bx , (2)where β = ( β , β , . . . , β k ) and B = β β / . . . β k / β / β . . . β k / β k / β k / . . . β kk / . (3)Box and Hunter (1954) constructed a confidence region for the stationary pointof a second-order response surface. The stationary point is given by x sp = − B − β , (4)ootstrap Confidence Regions 3the solution to the system of equations ∂y/∂ x = . The interpretation and relevanceof the stationary point depends on the nature of the response surface which is deter-mined by B . Consider the following cases: 1) the eigenvalues of B are mixed in sign,2) the eigenvalues of B are all positive and 3) the eigenvalues of B are all negative.In cases 1 and 2 the stationary point is a saddle point and the location of minimumresponse, respectively, not the location of maximum response and, therefore, is not ofinterest for the purpose of this report. In case 3 the stationary point is the locationof maximum response. In practice, the model parameters, including the elementsof B , are unknown but can be estimated from the data. Since there is uncertaintyassociated with any estimator of B there is also uncertainty in assessing the natureof the true response surface and, therefore, the interpretation of Box and Hunter’sconfidence region.Another issue of practical importance is the data are observed over a treatmentspace of finite dimensions, i.e. the experimental region, and one is generally unwillingto extrapolate to areas outside this region. Peterson (1992) utilized a transforma-tion technique to develop a confidence region for the stationary point constrainedto a specified region. However, as the approach is founded on the stationary point,the methodology suffers from the same interpretation difficulty as Box and Hunter’sconfidence region.In cases 1 and 2 the maximum response is undefined unless interest is restrictedto a subset of the treatment space. Ridge analysis was developed by Hoerl (1959)and refined by Draper (1963) to estimate the stationary point subject to the con-straint x ′ x = r . Under their approach the investigator must specify a value of theLagrangian multiplier µ which, if chosen greater than the largest eigenvalue of ˆ B ,facilitates determination of the location of the predicted constrained maximum. Sta-blein, Carter and Wampler (1983) constructed a confidence region for the constrained Roger Gibb, I-Li Lu, and Walter Carterstationary point conditional on the investigator’s choice of µ . However, proper in-terpretation of their confidence region depends on whether the choice of µ is greaterthan the largest eigenvalue of B . Since the eigenvalues of B can only be estimated,there is uncertainty involved in the confidence region’s interpretation.Let x cm be the operating conditions associated with maximum g ( x , θ ) subject tothe constraint that x is within the experimental region. Denote ˆ x cm and ˆ θ as estima-tors of x cm and θ respectively, then in practice, ˆ x cm can be calculated from g ( x ; ˆ θ )using a numerical optimization algorithm, such as the Nelder-Mead (1965) simplex.A favorable property of x cm is its interpretation does not depend on unknown modelparameters, unlike the stationary point or constrained stationary point. A drawback,however, is that in some instances ˆ x cm is not a consistent estimator. Consider thecase where g ( x , θ ) is continuous and ε i iid ∼ N(0 , σ ). If the model is correct, ˆ θ is aleast-squares estimator and x cm is unique, then ˆ x cm can be shown to be consistent(Kendall and Stuart, 1979, Chapter 18). To demonstrate an instance where x cm isnot unique and, therefore, ˆ x cm is not consistent, consider the k = 2 second-orderresponse surface where β = β = 0 and β = β >
0. If the experimental region isthe two dimensional direct product of the interval [ − a, a ], i.e. [ − a, a ] , then the max-imum response within the experimental region occurs at ( − a, − a ), ( a, − a ), ( − a, a )and ( a, a ). A formal assessment of whether x cm is unique could be made with testsof the appropriate hypotheses. For example, the likelihood of the symmetric modeldescribed earlier could be investigated by testing H : β = β = 0, β = β . If thereis insufficient evidence to reject H further investigation is warranted.Construction of a confidence region for x cm using classical methods would requireknowledge of the sampling distribution of ˆ x cm . As this distribution is unknown itis reasonable to explore the use of bootstrap confidence region methods that do notrequire its mathematical derivation. Unlike the confidence region methods describedootstrap Confidence Regions 5earlier, the bootstrap approach does not require normality assumptions on the modelerrors. Hall (1987, 1992) describes three methods for constructing likelihood-based bootstrapconfidence regions for a k -variate parameter vector ξ given a sample of size n , namely,the percentile- t method, the ordinary percentile (hybrid) method, and the percentilemethod. In regards to which of these methods is appropriate when ξ ≡ x cm , thereare several factors that merit consideration. First, the percentile- t method requiresan accurate estimate for V , the asymptotic covariance matrix of n (ˆ ξ − ξ ), where V is assumed to be positive definite. The analytic expression for ˆ V , an estimate of V , is unknown and accurate estimates using resampling techniques are unlikely toobtain in the small sample designed experiment setting. Second, and perhaps moreimportantly, only the percentile method preserves the range of ξ in all cases. Inparticularly, for ξ ≡ x cm , only confidence regions under the percentile method areguaranteed to be bounded by the experimental region. Therefore, for the purposesof this report, attention is focused exclusively on the percentile method to constructlikelihood-based bootstrap confidence regions for x cm . An adaptation of Hall’s (1987) likelihood-based bootstrap confidence region algorithmwhen ξ = x cm is given below, followed by several notes regarding its implementation.1. For the response surface model y = Xθ + ε , determine ˆ θ using the method ofleast-squares. Roger Gibb, I-Li Lu, and Walter Carter2. Standardize the residual vector ˆ ε with the elementwise operationˆ ε s = ˆ ε / p diag( I − X ( X ′ X ) − X ′ ).3. Compute b bootstrap estimates of x cm by following steps 3(a) - 3(d) a total of b times (see implementation notes for convenient choices of b ).(a) Generate a bootstrap sample of ˆ ε s , denoted by ˆ ε ∗ s .(b) Calculate the corresponding bootstrap sample of the data with y ∗ = X ˆ θ +ˆ ε ∗ s .(c) For the response surface model y ∗ = Xθ ∗ + ε , determine ˆ θ ∗ using themethod of least-squares.(d) Based on ˆ θ ∗ , calculate ˆ x ∗ cm .4. Fit a nonparametric density ˆ f to the ˆ x ∗ cm,i , i = 1, 2,. . . , b .5. The contour on ˆ f of smallest content that captures 100(1 − α )% of the b es-timates ˆ x ∗ cm is a 100(1 − α )% likelihood-based confidence region for x cm usingthe percentile method. Steps 1 − x cm . Bootstrapping residuals was proposed by Efron (1979) and is consid-ered appropriate when X is fixed and the model is correct with exchangeable errors.Standardization of the residual vector in step 2 ensures that the variance of eachˆ ε s,i matches that of the unobservable ε i . A simple bootstrap sample is generatedby drawing a sample of size n with replacement from the elements of ˆ ε s . The bal-anced bootstrap, proposed by Davison et al. (1986), was used in the examples andootstrap Confidence Regions 7simulation studies of section 2.4 to ensure that ˆ ε ∗ = ˆ ε . In step 4 the conditional dis-tribution of ˆ x cm is approximated nonparametrically and used as the basis for selectingthe confidence region boundary.In practice, the level contour on ˆ f of smallest content that captures 100(1 − α )%of the b estimates ˆ x ∗ cm was determined as follows: Choose b such that (1 − α ) b is anon-negative integer and assume that ˆ f (ˆ x ∗ cm,i ), i = 1, 2,. . . , b , are unique. Let ˆ f α bethe value of ˆ f corresponding to the level contour of smallest content on ˆ f ( x ) thatcaptures exactly (1 − α ) b of the ˆ x ∗ cm,i , i = 1, 2,. . . , b . To calculate ˆ f α consider that,by definition, ˆ f (ˆ x ∗ cm,i ) ≥ ˆ f α for all ˆ x ∗ cm,i captured by the contour and ˆ f (ˆ x ∗ cm,i ) < ˆ f α for every other ˆ x ∗ cm,i . Therefore, ˆ f α is the (1 − α ) b th element of the set with elementsˆ f (ˆ x ∗ cm,i ), i = 1, 2,. . . , b , sorted in descending order. In the event that all ˆ f (ˆ x ∗ cm,i ), i = 1, 2,. . . , b , are not unique, the contour ˆ f α ( x ) may capture more than (1 − α ) b of the ˆ x ∗ cm,i , i = 1, 2,. . . , b , in which case the confidence region is conservative.After identifying ˆ f α in this manner the confidence region boundary was graphicallydisplayed using the GCONTOUR procedure of SAS r . Likelihood-based confidence regions were constructed for x cm , where g ( x , θ ) was asecond-order polynomial ( k = 2) for the following cases: 1) concave down with x cm located inside the experimental region and 2) saddle system with x cm located, neces-sarily, on a boundary of the experimental region. The experimental region was definedas the two dimensional direct product of the interval [ − . , . x cm for both models are summarized in table 1.As indicated in section 2.1, a nonparametric estimate of the sampling distributionof ˆ x cm conditional on b bootstrap estimates of x cm serves as the basis for constructing Roger Gibb, I-Li Lu, and Walter Carterthe confidence region boundary. A product kernel density estimator for this purposeis given by ˆ f ( x ) = 1 bh h b X i =1 ( Y j =1 K (cid:18) x i − ˆ x ∗ cm,i,j h j (cid:19)) , (5)where K is a univariate kernel and h i is the bandwidth for the i th coordinate (see,e.g., Silverman, 1986, Scott, 1992). Gasser and M¨uller (1979) showed that if K issymmetric the bias of ˆ f ( x ) can be inflated near the boundaries. M¨uller (1988) showedthat uniform bias over the entire support of f ( x ) can be achieved if an asymmetricboundary corrected kernel is used in (5). Gasser and M¨uller (1979) and Jones (1993),among others, describe univariate boundary kernels designed for estimating densitieswith a single boundary. Densities with support over an interval may be estimatedusing two such kernels, one at each boundary, but this is appropriate only if thebandwidth is small relative to the interval length. Hart and Wehrly’s (1992) linearboundary kernel was implemented in (5), as this kernel is specifically designed toestimate densities with support on a interval, where the bandwidth need not be smallrelative the interval length. The biweight density was used as a basis for constructingthis boundary kernel.Two methods of multivariate bandwidth selection were explored: Normal rule-of-thumb (see, e.g., Scott, 1992) and the plug-in approach of Wand and Jones (1994,pp 107-113). The sample standard deviation in each coordinate of ˆ x ∗ cm was used toestimate the scale of the conditional distribution of ˆ x cm , as some measure of scaleis required by both bandwidth selectors. The possibility exists for all bootstrapestimates for x cm to be located on the experimental region boundary, a situation thatoccured almost exclusively with the saddle system. To ensure nonzero scale estimates,the standard deviation for the i th coordinate was calculated from a modified datasetootstrap Confidence Regions 9given by z i = x ∗ i,cm − sign( x ∗ i,cm ) δ if abs( x ∗ i,cm ) = 1 . x ∗ i,cm otherwise , (6)where δ is a random U (0 , .
05) deviate. The effect of using the modified datasetto estimate standard deviation was negligible in the concave-down system, as mostbootstrap estimates did not fall on a boundary. For the saddle system, using themodified dataset prevented nonzero bandwidth estimates and, thus, permitted use ofthe product kernel density estimator.
It is of interest to compare the bootstrap confidence region performance under bothmethods of bandwidth selection, i.e. Normal rule-of-thumb and the plug-in approach.Data for 500 experiments were simulated by adding a N(0 , ) deviate to the ‘true’response at operating conditions associated with a k = 2 rotatable central-compositedesign (5 center runs, n = 13). Confidence regions using both bandwidth selectorswere constructed using 2000 bootstrap samples in each simulated experiment. Band-widths under the plug-in method were slightly larger, on average, than those of theNormal rule-of-thumb (see table 2). This result is not unexpected, as simulationstudies of Wand and Jones (1994) indicate a tendency for the multivariate plug-inmethod to oversmooth in some cases.A contour of the true concave-down response surface and a representative con-fidence region for x cm with b = 2000 are provided in figures 1 and 2, respectively.Bootstrap estimates for x cm are indicated with points in figure 2, with a solid lineidentifying the confidence region boundary. Note that the confidence region bound-ary is within the experimental region, as required. Had the hybrid or percentile- t methods been implemented this would not have been true in all cases.0 Roger Gibb, I-Li Lu, and Walter CarterThe simulated experiments were arranged into five groups of 100 and the coverageprobability calculated in each group for confidence coefficients ranging from 90 − x cm are provided in figures 4 and 5, respectively. Note that in this case x cm islocated on the boundary of the experimental region. Since all bootstrap estimates for x cm are also on the boundary, visual identification of the confidence region is difficult.In the figure, two small triangles identify the location of lower and upper boundariesin the x coordinate. In the x coordinate the confidence region extends slightly awayfrom the experimental region boundary.Figure 6 summarizes the coverage probability under both bandwidth selectorsfor the saddle system. As with the concave-down response surface, the coverageprobability is statistically indistinguishable under the two approaches. The coverageprobability is closer to the nominal level than for the concave-down response surface.Calculation of the Normal rule-of-thumb bandwidths requires less time and com-putational resources than bandwidths under Wand and Jones’ plug-in method. Sincethe coverage probability was essentially the same under both methods, the Normalrule-of-thumb method was implemented in the simulation studies that follow.A simulation study was conducted to assess whether increasing the number ofbootstrap samples would improve the coverage probability. Confidence regions wereconstructed from the same simulated data described earlier using b = 4000 and b =6000. In neither the saddle system nor concave-down case was the coverage probabilityootstrap Confidence Regions 11significantly effected by increasing b .To explore the effect of sample size on coverage probability, confidence regions for x cm were constructed under both second-order response surfaces for n = 13, n = 26and n = 208, corresponding to 1, 2 and 16 experimental replications of a rotatablecentral-composite design with 5 center runs. The results for the concave-down andsaddle system with b = 2000 are summarized in figures 7 and 8, respectively. In bothcases the coverage probability approaches the nominal level with increasing samplesize. Existing confidence regions for operating conditions associated with the maximumof a response surface, either unconstrained or constrained to a specified region, arebased on the stationary point. These approaches are only applicable to second-ordermodels under the assumption of normally distributed errors. Also, the interpretationof these confidence regions can be ambiguous, since this requires assessment of theunknown elements of the B matrix.In contrast, the interpretation of a likelihood-based bootstrap confidence regionfor x cm does not depend on the nature of the response surface, assuming x cm isunique. For example, if g ( x , θ ) is a second-order polynomial the confidence regioninterpretation is the same whether g ( x , θ ) is concave-down, concave-up or a saddlesystem. In addition, the approach is not restricted to second-order models nor doesit require assumptions on the model error distribution, except for exchangeability ofthe errors. In principle, the methodology is also applicable to models where g ( x , θ ) isnonlinear in θ , though Hjorth (1994, pp 190) indicates that in non-linear regressionapplications a direct analogue of standardized residuals is generally not available.2 Roger Gibb, I-Li Lu, and Walter CarterSimulation results from section 2.4 (see figure 7) indicate the bootstrap confidenceregion coverage probability was less than nominal under the concave-down system( n = 13). Coverage probability was higher for the saddle system ( n = 13), but stillbelow the nominal level (see figure 8). A possible reason that higher coverage proba-bility was observed for the saddle system is the estimated conditional distribution ofˆ x cm is essentially restricted to one dimension in this case, thereby reducing the ‘curseof dimensionality’. Evidence that the coverage probability in both systems convergesto the nominal level with increasing sample size, an indication of the asymptoticaccuracy of the methodology, is apparent in figures 7 and 8.Loh (1987, 1991) proposed a bootstrap calibration technique to improve the ac-curacy of confidence sets. In brief, the approach consists of estimating the confidencecoefficient associated with the desired coverage probability. For example, to achieve acoverage probability of 90% in the concave-down response surface case ( n = 13), it isapparent in figure 7 that a confidence coefficient of approximately 99.5% is required.It is also apparent that in this case realization of coverage probabilities greater than ≈
91% are not feasible with this approach, unless the sample size is increased. Forthe saddle system ( n = 13) a confidence coefficient of approximately 95% is neededto achieve a coverage probability of 90% (see figure 8).The accuracy of bootstrap confidence regions for ˆ x cm using the percentile methodis largely determined by how accurately the conditional distribution of ˆ x cm , f , isestimated. Under the kernel method implemented in section 2, choice of bandwidthhas a critical impact on this accuracy. Hall (1987) indicated that bandwidths whichare optimal in a global sense, such as minimization of the mean integrated squarederror of ˆ f , can produce confidence regions that are “unduly bumpy”. He noted thatthis can be attributed to the fact that the ratio of the variance to the squared biasof ˆ f is greater in the tails than in the central region of the distribution. Therefore,ootstrap Confidence Regions 13the Normal rule-of-thumb and plug-in bandwidth selectors were investigated for theirtendency to avoid undersmoothing. The k = 2 concave-down and saddle systemcases explored in section 2.4 are examples where these methods provide reasonablebandwidth estimates. Jones (1993) followed a similar approach when he used a plug-inbandwidth selector in conjunction with a univariate boundary kernel.There are situations, however, in which bandwidth selection is more challenging.For example, if f is bimodal, use of the standard deviation to estimate the scale of f can result in considerable oversmoothing. Janssen et al. (1995) developed morerobust scale estimators for such cases. Consider, also, a situation where all bootstrapestimates for x cm are equally distributed on the two boundaries A = { x : x =1 . , x ∈ ( − . , . } and B = { x : x ∈ ( − . , . , x = 1 . } . On boundary A bandwidth h should be much smaller than h . However, on boundary B the oppositeis true. In this case, a variable kernel density estimator would seem more appropriatethan the fixed bandwidth kernel estimator implemented in section 2.4, as this wouldallow for different levels of smoothing depending on the location of bootstrap estimatesfor x cm .Application of the bootstrap confidence region methodology was restricted to the k = 2 case where the experimental region was rectangular in shape. The bootstrapand kernel density estimation techniques described in section 2 are not limited to twodimensions. Wand and Jones’ plug-in bandwidth selector can also be extended tohigher dimensions, though with greater computational expense. In light of the closecomparison of coverage probability under the Normal rule-of-thumb and plug-in meth-ods observed in section 2.4, rule-of-thumb methods may be adequate in many higher-dimensional applications. Graphical presentation of confidence regions for k > D , a simultaneous measure of the desirability of all the responses. Anarguable shortcoming of the methodology is it does not provide for estimation of thevariability of these operating conditions. Bootstrap confidence region methodology isdefined in sufficiently general terms to allow for construction of a confidence regionfor the operating condition that maximizes D constrained to the experimental region,under the assumption that this parameter is unique.Of practical interest is the computational resources and time required to constructa likelihood-based bootstrap confidence region. The programs that generated theresults of section 2.4 were written in SAS r and run on a 300 megahertz Pentium r IIdesktop computer with 128 megabytes of RAM. Approximately 2 minutes is requiredto generate a single confidence region with b = 2000, of which about 25 seconds areused for steps 1 − r . It is expected that steps 4 and 5 would requireconsiderably less time under their optimized procedures.ootstrap Confidence Regions 15 References
Box, G.E.P., and Hunter, J.S. (1954), “A Confidence Region for the Solution ofa Set of Simultaneous Equations With An Application to Experimental Design,”
Biometrika , 41, 190-199.Davison, A.C., Hinkley, D.V., and Schechtman, E. (1986), “Efficient Bootstrap Sim-ulation,”
Biometrika , 73, 555-566.Derringer, G., and Suich, R. (1980), “Simultaneous Optimization of Several ResponseVariables,”
Journal of Quality Technology , 12, 214-219.Draper, N.R. (1963), “Ridge Analysis for Response Surfaces,”
Technometrics , 5, 469-479.Efron, B. (1979), “Bootstrap Methods: Another Look at the Jackknife,”
The Annalsof Statistics , 7, 1-26.Gasser, T., and M¨uller, H.-G (1979), “Kernel Estimation of Regression Functions,” In
Smoothing Techniques for Curve Estimation (eds. T. Gasser and M. Rosenblatt),Heidelberg: Springer-Verlag, pp 23-69.Hall, P. (1987), “On the Bootstrap and Likelihood-Based Confidence Regions,”
Biometrika ,74, 481-493.Hall, P. (1992),
The Bootstrap and Edgeworth Expansion , New York: Springer-Verlag.Harrington, E.C. (1965), “The Desirability Function,”
Industrial Quality Control , 21,494-498.Hart, J.D. and Wehrly, T.E. (1992), “Kernel Regression When the Boundary Regionis Large, With An Application to Testing the Adequacy of Polynomial Models,”
The Journal of the American Statistical Assocication , 87, 1018-1024.Hoerl, A.E. (1959), “Optimum Solution of Many Variables Equations,”
ChemicalEngineering Progress , 55, 67-78.6 Roger Gibb, I-Li Lu, and Walter CarterHjorth, J.S.U. (1994),
Computer Intensive Statistical Methods, Validation Model Se-lection and Bootstrap , London: Chapman & Hall.Janssen, P., Marron, J.S., Veraverbeke, N., and Sarle, W. (1995), “Scale Measuresfor Bandwidth Selection,”
The Journal of Nonparametric Statistics , 5, 359-380.Jones, M.C. (1993), “Simple Boundary Correction for Kernel Density Estimation,”
Statistics and Computing , 3, 135-146.Kendall, M. and Stuart, A. (1979),
The Advanced Theory of Statistics, Volume II:Inference and Relationship,
The Journal of the Amer-ican Statistical Assocication , 82, 155-162.Loh, W.-Y. (1991), “Bootstrap Calibration for Confidence Interval Construction andSelection,”
Statistica Sinica , 1, 479-495.M¨uller, H.G. (1988).
Lecture Notes in Mathematics: Nonparametric Regression Anal-ysis of Longitudinal Data.
New York: Springer-Verlag, 46.Nelder, J.A., and Mead, R. (1965), “A Simplex Method for Function Minimization,”
Computer Journal , 7, 308-313.Peterson, J.J. (1992), “Confidence Regions for Constrained Response Surface Op-tima,” Presented at the August, 1992 Joint Statistical Meetings, Atlanta, GA.Scott, D.W. (1992),
Multivariate Density Estimation: Theory, Practice and Visual-ization , New York: Wiley.Silverman, B.W. (1986),
Density Estimation for Statistics and Data Analysis , London:Chapman & Hall.Stablein, D.M., Carter, W.H., and Wampler, G.L. (1983), “Confidence Regions forConstrained Optima in Response-Surface Experiments,”
Biometrics , 39, 759-763.Staniswalis, J.G., Messer, K., and Finston, D.R. (1993), “Kernel Estimators for Mul-tivariate Regression,”
Journal of Nonparametric Statistics , 103-121.ootstrap Confidence Regions 17Wand, M.P. (1994), “Fast Computation of Multivariate Kernel Estimators,” Journalof Computational and Graphical Statistics , 3, 433-445.Wand, M.P., and Jones M.C. (1994), “Multivariate Plug-In Bandwidth Selection,”
Computational Statistics Quarterly , 9, 97-116.Table 1: True second-order model parameters and associated x cm .ResponseSurface β β β β β β x cm concave down 86.850 5.242 4.778 -0.775 -2.781 -2.524 (0 . , . − . , . h h h h Normal rule-of-thumb 0.196 0.213 0.011 0.261(0.0033) (0.0037) (6 . × − ) (0.0049)Plug-in method 0.214 0.233 0.013 0.313(0.0034) (0.0038) (13 × − ) (0.0060)8 Roger Gibb, I-Li Lu, and Walter CarterFigure 1: Contour plot of the true concave-down response surface, where × identifies x cm .Figure 2: A 90% confidence region for x cm where the true response surface is concave-down and × identifies x cm .ootstrap Confidence Regions 19Figure 3: Comparison of coverage probability under the Normal rule-of-thumb andplug-in bandwidth selectors. The true response surface is concave-down.Figure 4: Contour plot of the true saddle system response surface, where × identifies x cm .0 Roger Gibb, I-Li Lu, and Walter CarterFigure 5: A 90% confidence region for x cm where the true response surface is saddlesystem and × identifies x cm .Figure 6: Comparison of coverage probability under the Normal rule-of-thumb andplug-in bandwidth selectors. The true response surface is a saddle system.ootstrap Confidence Regions 21Figure 7: Comparison of coverage probability for three different sample sizes wherethe response surface is concave-down and b = 2000 with Normal-rule-of-thumb band-widths.Figure 8: Comparison of coverage probability for three different sample sizes wherethe response surface is a saddle system and bb