A graphical approach for evaluating and comparing designs for nonlinear models

COMPUTATIONAL STATISTICS & DATAANALYSIS ELSEVIER

Computational Statistics & Data Analysis 27 (1998) 433-443

A graphical approach for evaluating and comparing designs for nonlinear models Andr6 I. Khuri *, Juneyoung Lee Department of Statistics, University of Florida, P.O. Box 118545, 103 Griffin-Floyd Hall, Gainesville, FL 32611-8545, USA Received 1 October 1996; accepted 1 September 1997

Abstract

A graphical approach is proposed for the purpose of assessing the overall prediction capability of a nonlinear design inside a region of interest R. A visual description is given of the behavior of the meansquared error of prediction throughout the region R. More explicitly, for each of several concentric surfaces inside R, quantile plots of the so-called estimated scaled mean-squared error of prediction (ESMSEP) are obtained. In addition, when the number of control (input) variables in the associated model is equal to one, plots of the maximum and minimum of the scaled mean-squared error of prediction (SMSEP) over a subset of the model's parameter space are developed. The latter plots are primarily used to compare several nonlinear designs. Two examples are presented to demonstrate the usefulness of the proposed plots. (~) 1998 Elsevier Science B.V. All rights reserved.

Keywords." Nonlinear designs; Prediction bias; Prediction capability; Quantile plots; Scaled mean-squared error of prediction

I. Introduction The choice of design for a nonlinear model is based, for the most part, on criteria that pertain to the estimation of 0, the vector of unknown model parameters. Prominent among such criteria is local D-optimality, which was introduced by Box and Lucas (1959). If the mean response is approximately linear in the elements of 0 inside a neighborhood of a point 00 in the parameter space, then a locally D-optimal design minimizes the volume of a first-order approximation of the exact confidence * Corresponding author. Tel.: 352392 1941; fax: 3523925175; e-mail: [email protected]. 0167-9473/98/$19.00 (~) 1998 Elsevier Science B.V. All rights reserved

PII S 0 1 6 7 - 9 4 7 3 ( 9 8 ) 0 0 0 1 6 - 4

434

A.L Khuri, J. Lee~Computational Statistics & Data Analysis 27 (1998) 433-443

region for 0. Designs that minimize a second-order approximation of the volume of this region were introduced by Hamilton and Watts (1985). An undesirable feature of these designs is their dependency on the unknown value of 0. Furthermore, the designs consist of only p distinct points of support, where p is the number of elements of 0. It is therefore not feasible to check the adequacy of the hypothesized model. More recently, O'Brien (1992) used Hamilton and Watt's (1985) design criterion to construct designs with p + 1 support points. A review of optimal experimental designs for nonlinear models was given by Ford et al. (1989). See also Haines (1993). Due to the predominance of the local D-optimality criterion, it has become commonplace to evaluate and compare designs on the basis of their D-efficiency values. This criterion, however, produces a single-valued index, which may not be sufficient to adequately assess the quality of information that can be derived from a given design. The focus in D-optimality is on the variances of the estimates of the model's parameters. In some cases, however, it may be more meaningful to evaluate a design on the basis of its capability to provide good predictions throughout a region of interest R. Giovannitti-Jensen and Myers (1989) used the so-called variance dispersion graphs (VDGs) to assess the overall prediction capability of a design for fitting a linear model. These graphs consist of two-dimensional plots that display the maximum and minimum of the prediction variance function on each of several concentric spheres (hyperspheres, in general) inside a region of interest. More recently, Khuri et al. (1996) proposed the use of quantile plots of the prediction variance on concentric spheres. The latter plots provide more information concerning the quality of prediction than can be obtained from the VDGs. In this article, particular plots of the mean-squared error of prediction are developed throughout a region of interest for a nonlinear model. These plots can be used to evaluate the prediction capability of a given design, as will be seen in Section 3. They can also be used to compare designs. A demonstration of such plots is given in Section 4 for models with only one control (input) variable.

2. The mean-squared error of prediction In this section, an approximate expression is developed for the mean-squared error of a predicted response represented by a nonlinear model of the form y = f ( x , O ) + a,

(2.1)

where y is an observed response, f is a nonlinear function differentiable with respect to 0, x = (Xl,X2,... ,xk)' is a vector of k control (input) variables, 0 = (01,02,..., Op)' is a vector of p unknown parameters, and e is a random experimental error with mean zero and variance tr1. Suppose that model (2.1) is fitted using a design D consisting of n points, denoted by x~,x2 .... ,xn, inside a region of interest R. Let /) denote the least-squares estimate of 0. The predicted response 33 is given by .9 = f ( x , 0),

(2.2)

A.L Khuri, J. Lee~Computational Statistics & Data Analysis 27 (1998) 433-443

435

It is known that 0 is a biased estimator of 0. Several authors developed approximate expressions for the bias vector E(to), where tO----0- 0 (see Box, 1971; Cook et al., 1986). For example, from Cook et al. (1986, formula (3)) we approximately have

(2.3)

E(to) ~ (F'F)-IF'd,

where F is an n x p matrix having Of(xu, 0)/80i as its (u, i)th element, u = i, 2 , . . . , n; i = 1 , 2 , . . . , p , d is an n × 1 vector with elements --(0"2/2) tr[(F'F)-IH~], and Hu is a p × p matrix whose (i,j)th element is O2f(xu, O)/OOiOOj, u= 1,2,...,n; i,j= 1,2,..., ~o. Following Cook et al. (1986, formula (4)), a second-order approximation of f(x, O) in a small neighborhood around 0 is given by

f(x, O) ~--f(x, O) + g'(x, O)to + ~to 1 t Hto,

(2.4)

where

g'(x, O) = -Of (x,

, t3f (x,

" " ' ~f (x, O)

and H is a p × p matrix with the (i,j)th element O2f(x,O)/OOiOOj, i,j= 1 , 2 , . . . , p . Using (2.3) and (2.4), the prediction bias, E [ f ( x , O ) ] - f(x,O), is approximately given by

E[f(x, ~J)] - f(x, O) ~_g'(x, O)E(to) + ½E(to'Hto) =g'(x,O)E(to) +

½E(to')HE(to)+

i tr[H Var(0)].

(2.5)

Formula (2.5) follows from applying Theorem 1 in Searle (1971, p. 55), namely

E(to'Hto) = E(to')HE(d~)+tr[H Var(to)], and noting that Var(to) = Var(0). The latter variance-covariance matrix is approximately equal to (F'F)-10- 2. From (2.3) and (2.5) we then have

E[f(x, 0)] - f(x, O) ~ g'(x, O)(F'F)-~F'd + ½d'F(F'F)-IH(F'F)-~F'd +~1 tr[H(F,F)-10-2].

(2.6)

Now, the mean-squared error of f(x, O) is of the form M S E [ f ( x , ~J)] -- E[f(x, O) - f(x, 0)] 2 = Var[f(x, CJ)] + {E[f(x, ~J)] - f(x, 0)} 2, which can be approximately expressed as M S E [ f ( x , ~J)] _~ g'(x, O)(F'F)-~g(x, 0)0-2 + {g'(x, O)(F'F)-'F'd

+ ½d'F(F'F)-~H(F'F)-IF'd + l tr[H(F'F)-10-2]}2.

(2.7)

A.I. Khuri, 3. Lee I Computational Statistics & Data Analysis 27 (1998) 433-443

436

Henceforth, we shall consider the so-called scaled mean-squared error of prediction (SMSEP), zD(x, to), which we define as n

"cD(x,to) = ~

MSE[f(x, 0)1 n

t

~- ng'(x, O)(F'F)-lg(x, O) + -~{g (x, O)(F'F)-IF'd

+½a'r(r'V)-ln(F'V)-lF'd + l tr[H(F'F)-lrr2]} 2,

(2.8)

where to = (0', o'2)' and n is the number of observations. This results in a meansquared error free of scale and put on a "per observation" basis. Note that ZD(X,tO) depends on the elements of to, which are unknown. An estimate of zD(x, to) can be obtained by replacing to with tb=0J',~2) ', where ~ 2 = M S e is the residual mean square associated with the fitted model. We refer to zD(x, tb) as the estimated scaled mean-squared error of prediction (ESMSEP).

3. Evaluation of a design's prediction capability In this section, a graphical technique is presented to assess the prediction capability of a given design D over some region R. The technique is based on using plots of quantiles of the ESMSEP on each of several concentric surfaces within R. The surfaces are constructed by repeated reductions of the boundary of R using different values of a shrinkage factor 2. A demonstration of the proposed technique follows. Consider the nonlinear model y=

0103xl

+e,

(3.1)

1 + OlX 1 "~- 02X 2

which describes a certain chemical reaction, where y is the reaction rate, xl and x2 are partial pressures of a reactant and a product, respectively; 01 and 02 are absorption equilibrium constants for the reactant and the product, respectively, and 03 is the effective reaction rate constant. This model was cited by several authors, for example, Box (1971, p. 181), Draper and Smith (1981, p. 521), Fedorov (1972, p. 193), and Seber and Wild (1989, p. 259). Observations on y are obtained within the region R given by R = {(Xl,X2)10 _< Xl _< 3, 0 < x2 _< 2}.

(3.2)

The settings of a design D and the corresponding response values for model (3.1) are shown in Table 1. The same data were reported in Draper and Smith (1981, p. 521). The prediction capability of the design D is measured by the size of the scaled mean-squared error of prediction, zo(x, to), given in formula (2.8), where, in this example, to = (01, 02, 03, a2) '. Small values of zD(x, to) are desirable throughout the region R. Using the data in Table 1, the least-squares estimates o f 01,02, 03 are 0~ = 3.5690, 02 = 12.7957, 03 = 0.6295. Also, an estimate of a 2 is ~iven b y MSE = 0.000788 with 10 degrees of freedom. Replacing to with ~b=(OI,02,03,MSE)' in (2.8) produces the ESMSEP, rD(X,~b).

A.I. KhurL J. Lee/Computational Statistics & Data Analysis 27 (1998) 433-443

437

Table 1 Design settings and response values for model (3.1) xl

x2

y

1.0 2.0 1.0 2.0 0.1 3.0 0.2 3.0 0.3 3.0 3.0 0.2 3.0

1.0 1.0 2.0 2.0 0.0 0.0 0.0 0.0 0.0 0.8 0.0 0.0 0.8

0.126 0.219 0.076 0.126 0.186 0.606 0.268 0.614 0.318 0.298 0.509 0.247 0.319

S~

lambda=1.0 0.07

0.93 0.9

Q,7

0.6

0 i

0

1

2

3

Xl

Fig. 1. Concentric rectangles within the region R (e is location of maximum ESMSEP on S~ and o is location of minimum).

Since the region R in (3.2) is a rectangle, we consider values of ZD(X,tb) on each of several concentric rectangles centered at the point (1.5, 1.0). Each rectangle is obtained by a reduction of the boundary of R using a shrinkage factor 2, 0.50 < 2 < 1. The resulting rectangle is denoted by S~ (see Fig. 1). For each of several selected

438

A.I. Khuri, J. Lee~Computational Statistics & Data Analysis 27 (1998) 433~143

0

O

0.2

i

i

0.4

0.6

0.8

1.0

P

Fig. 2. Quantile plots of ESMSEP for model (3.1) on S~.

values of 2, points are generated at random on Sz using a uniform distribution (the number of points generated on each side of Sz is 2500 giving a total of 10000 points on S;,). Values of Zo(X,tb) are then obtained at the generated points on Sx for the given design. The quantiles of the resulting values can be computed using, for example, the S-PLUS language. Let Qo(p,2) denote the corresponding pth quantile. Plotting QD(p, 2) versus p for selected values of 2 gives rise to the so-called quantile plots of the ESMSEP. Fig. 2 shows a display of such plots for the selected values of 2. The actual plots were obtained by using a computer program written in the S-PLUS language. Cubic spline interpolation was used to connect points on the plots. An overall view of the distribution of Zo(X,tb) within the region R can be seen from Fig. 3, which gives the box plots of values of to(x, ¢b) on S;,, 2=0.6(0.1)1. For further details concerning the use of S-PLUS with regard to the proposed methodology, see Becker et al. [1988, p. 402 (boxplot), p. 588 (spline)]. From Fig. 2 we note that little variation exists in the quantile values for small to moderate values of 2. Relatively larger variations occur for values of 2 near 1. Similar information can be obtained from the box plots in Fig. 3. The maximum ESMSEP value over R is 4.25 and is attained on the boundary of R where 2--- 1. Thus, as we move away from the boundary closer to the center of R, the values of the ESMSEP become more stable. It is also interesting to note the stability in

A.1. KhurL J. Lee~Computational Statistics & Data Analysis 27 (1998) 433-443

439

, 4.280

3.254

3.176

3.082

o'J , 2.812 (a)

n IJJ O3 CO UJ

c~l

II

2.406 2.231 (b) 2.207 (c)

2.224

2,191 2,093

2.131

1.970 1.833

1.683

I 0.900

0.406 J 0.101

0

0.0

0.6

0.7

Fig. 3. Box plots of ESMSEP for model

0.8

( 3 . 1 ) (a = m a x i m u m , b =

0.9

median,

1

c = average, d =

minimum).

the values of the median. However, the average ESMSEP on S~ is somewhat more variable than the median (see Fig. 3). The advantage of using the quantile plots in Fig. 2 over the box plots is that the entire distribution of the ESMSEP on each S~ is displayed by the former while only specific values of that distribution are given by the latter. Fig. 1 gives locations of the maximum and minimum of the ESMSEP on S~ for 2 = 0.60, 0.70, 0.80, 0.90, 0.93, 0.97, 1.0. The values of the maxima increase with 2 while those for the minima decrease (see Fig. 3). Using these locations, it is possible to plot a path of increasing maximum values, and a path of decreasing minimum values of ESMSEP inside the region R. Note that the maximum value over R, namely 4.25, is attained at the lower right-hand comer of R. Furthermore, the minimum value over R, namely 0, occurs at all points on the left-hand side of the boundary of R where Xl = 0(2---- 1). This is true because at such points, the first-order and secondorder partial derivatives of f(x, O) for model (3.1) are equal to zero. Consequently, g(x,O),H, and hence MSE[f(x,O)] in formula (2.7), have the value zero.

4. Comparison of designs In this section, comparisons among nonlinear designs will be undertaken using plots of the scaled mean-squared error of prediction. No estimate of to will be assumed available. This is particularly needed here since, when comparing several designs, it may not be realistic to run a series of experiments for each design

440

A.L Khuri, J. Lee/Computational Statistics & Data Analysis 27 (1998) 433-443

considered. This is obviously a more difficult problem than the one in Section 3. We shall therefore restrict consideration of this problem to the special case of nonlinear models with only one control (input) variable. Consider, for example, the nonlinear model

y=Os(e - ° ' x - e - ° 2 x ) + e ,

x>_O,

(4.1)

where 0 < 0 1 < 0 2 , and 03>0. Model (4.1) represents a special case of a class of models that arise frequently in the study of pharmacokinetics, where x is time (in minutes) and y is concentration. This model is discussed more fully in Atkinson and Donev (1992, Example 18.4, p. 193). A slightly more general version of this model was used by Bates and Watts (1988, p. 101). Let us compare the prediction capabilities of the following two designs: DI: x = 0.166,0.333,0.5,0.666, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 24.0, 30.0, 48.0, D2: x = 0.2288, 0.2288, 0.2288, 0.2288, 0.2288, 0.2288, 1.3886, 1.3886, 1.3886, 1.3886, 1.3886, 1.3886, 18.417, 18.417, 18.417, 18.417, 18.417, 18.417. Design D1 was used by Atkinson and Donev (1992, Table 18.1, p. 202), and design D2 is locally D-optimal (see Atkinson and Donev, 1992, Table 18.2, p. 205). For both designs, the region of interest R is the closed interval [0, 50]. The comparison between D1 and O2 is based on the value of zD(x, to), the scaled mean-squared error of prediction (SMSEP) in formula (2.8). Since to = (01,02, 03, a2) ' is unknown, the following procedure is proposed: Let C denote a subset of the parameter space of to. In this example, we consider C as C = {tol0<01 ~ 0.10, 01 <02 ~ 6.0, 15 <03 < 25, 0
(4.2)

For each Di (i = 1,2) and each x E R, x = 0(5)50, the minimum and maximum values of rD,(X, to), with respect to to over C, are obtained. Consider the functions OD,(x) = micn{ZD,(x,to)},

i----- 1,2,

hDi(x) = max{zo~(x, to)},

i = 1,2.

EC

Plotting values of 9D,(X) and hD,(x) (i = 1,2) against x results in the so-called minimum and maximum SMSEP plots. A combination of such plots for designs DI and D2 is given in Fig. 4. We note that for x _< 25, the variability in the SMSEP values is small for both designs. Maximum dispersion in the values of SMSEP occurs on the boundary of R. From these plots we can clearly see that more stable and smaller maximum SMSEP values are achieved with design D1 than with the locally D-optimal design D2. This holds true especially for large values of x. Design D1 can therefore be considered superior to D2 in terms of predictive ability. Incidently, the D-efficiency of D1 is only 67.61 (see Atkinson and Donev, 1992, Table 18.3, p. 206).

441

A.I. Khuri, J. Lee/Computational Statistics & Data Analysis 27 (1998) 433-443

O CM

.,,"

.......... design D1

...," ...."

design D2 J

tO

/.,." ./," 13. LU

.,"

,,"

,j.'*

p;

.... .

f.O

tt~

C)

0

10

20

30

40

50

X

Fig. 4. Combined minimum and maximum SMSEP plots for designs D1 and D2 (based on the subset C in (4.2)).

In order to determine the sensitivity (or robustness) of the above results to changes in C, another subset, C*, of the parameter space of to is considered, namely, C* = {tolO < 01 _< 0.50, 01 < 02 _< 10.0, 25 < 03 _< 35, 0 < 0-2 ~ 5.0},

(4.3)

The corresponding combined plots for D 1 and D2 are given in Fig. 5. We note that the shapes of the plots have changed. However, the conclusion concerning the superiority of D1 to /92 remains valid. We also note that the SMSEP values corresponding to D2 are quite dispersed for x < 10. By contrast, D1 exhibits small values and a greater degree of stability throughout the region R.

5. Concluding remarks The graphical techniques introduced in Sections 3 and 4 provide an alternative procedure for assessing and comparing nonlinear designs on the basis of the mean-squared error of prediction criterion. The examples given in these sections demonstrate that the proposed plots can provide quite detailed information about the performance of designs over a region of interest R. By contrast, designs generated on the basis of a single-valued criterion, such as local D-optimality, may not perform as well as other designs in terms of the mean-squared error of prediction throughout the region R. This was clearly shown to be the case in the example given in

442

A.1. Khuri, J. Lee~Computational Statistics & Data Analysis 27 (1998) 433-443

0

""\

design D1 .......... design D2

0

tlJ

C~

IIIIT,

0 ,

,

0

10

20

T I . ..................:..................:..................: ,

,

1

30

40

50

Fig. 5. Combined minimum and maximum SMSEP plots for designs D1 and D2 (based on the subset C* in (4.3)).

Section 4 using plots of maximum and minimum values of the SMSEP. These plots also provide a clear depiction of the dependency of a nonlinear design's quality of prediction on the parameters of the associated model. It should be noted that the quantile plots used in Section 3 provide a description of the entire distribution of the ESMSEP on each of several concentric surfaces (S~) inside the region R. They are therefore very useful in assessing the overall prediction capability of a given design within R. These plots can be used with any number of control variables. However, they require the availability of an estimate of to [see formula (2.8)]. Such a requirement is not needed for the plots in Section 4, but their applicability is restricted to models with only one control variable. An extension of the latter plots to models with more than one control variable is possible, but is computationally very intensive as described below. The techniques used in Sections 3 and 4 can be combined in order to compare two designs, D1 and D2, without having to actually run experiments for the purpose of estimating to. This is accomplished by finding min~,ec{ZD,(x, to)} and maxo,ec{Zo,(X, to)} for each x on S~, as was done in Section 4, where Sa is one of several concentric surfaces inside R(i = 1,2). Quantile plots of values of each of these minimum and maximum functions on S;~ can then be obtained and compared for the two designs. This process, however, would require an extensive amount of computation, especially for models with several parameters and a large number of control variables.

A.1. Khuri, J. Lee~Computational Statistics & Data Analysis 27 (1998) 433-443

443

References Atkinson, A.C., Donev, A.N., 1992. Optimum Experimental Designs. Oxford University Press, Oxford. Bates, D.M., Watts, D.G., 1988. Nonlinear Regression Analysis and its Applications. Wiley, New York. Becker, R.A., Chambers, J.M., Wilks, A.R., 1988. The New S Language. Wadsworth & Brooks/Cole, Pacific Grove, CA. Box, G.E.P., Lucas, H.L., 1959. Design of experiments in non-linear situations. Biometrika 46, 77-90. Box, M.J., 1971. Bias in nonlinear estimation. J. Roy. Statist. Soc. B33, 171-190. Cook, R.D., Tsai, C.L., Wei, B.C., 1986. Bias in nonlinear regression. Biometrika 73, 615-623. Draper, N.R., Smith, H., 1981. Applied Regression Analysis. 2nd edn. Wiley, New York. Fedorov, V.V., 1972. Theory of Optimal Experiments. Academic Press, New York. Ford, I., Titterington, D.M., Kitsos, C.P., 1989. Recent advances in nonlinear experimental design. Technometrics 31, 49-60. Giovannitti-Jensen, A., Myers, R.H., 1989. Graphical assessment of the prediction capability of response surface designs. Technometrics 31, 159-171. Haines, L.M., 1993. Optimal design for nonlinear regression models. Comm. Statist. Theory Methods 22, 1613-1627. Hamilton, D.C., Watts, D.G., 1985. A quadratic design criterion for precise estimation in nonlinear regression models. Technometrics 27, 241-250. Khuri, A.I., Kim, H.J., Urn, Y., 1996. Quantile plots of the prediction variance for response surface designs. Comput. Statist. Data Anal. 22, 395-407. O'Brien, T.E., 1992. A note on quadratic designs for nonlinear regression models. Biometrika 79, 847-849. Searle, S.R., 1971. Linear Models. Wiley, New York. Seber, G.A.F., Wild, C.J., 1989. Nonlinear Regression. Wiley, New York.