Journal of Statistical North-Holland
Planning
and Inference
57
42 (1994) 57-77
Some alphabetic optimal designs for the logistic regression model William
R. Myers
Department of Biometrics OH 45241-2422, USA
Raymond Department
Walter
and Statistical
The Procter & Gamble Company,
Cincinnati,
H. Myers of Statistics,
Virginia Polytechnic
H. Carter
Department of Biostatistics, VA 23298-0032, USA
Received
Sciences,
24 August
Institute
and State University,
Blacksburg,
VA 24061, USA
Jr. Medical
College of Virginia,
1992; revised manuscript
received
Virginia Commonwealth
University,
Richmond,
8 April 1993
Abstract Alphabetic optimality has become an important component of experimental design in the case of the standard linear model. Many design criteria have been developed in order to produce optimal designs based on either parameter estimation or the prediction of a response. However, relatively little attention has been devoted to developing designs in the nonlinear model case (e.g. logistic regression model). D-optimality is the only criterion that has received much attention in the literature for the logistic regression model. This paper develops optimal designs for the logistic model based on the prediction of a response. Comparisons are made between different optimal designs by computing efficiencies based on certain optimality criteria. The robustness properties for some of the newly created designs are also investigated. Finally, the development of optimal designs for asymmetric regions of the logistic regression model is presented. AMS
Subject Classification:
Key words: Optimal-design
62H17. theory; Q-optimality;
average prediction
variance;
region of interest; dose levels.
1. Introduction The importance of optimal-design theory for the linear model has been well documented in the statistical literature. Kiefer and Wolfowitz (1959) were instrumental in the early development of this particular area of experimental design. These two
Correspondence to: Prof. W.H. Carter Jr., Professor and Chairman, Department of Biostatistics, 32, Medical College of Virginia, Virginia Commonwealth University, Richmond, VA 23298.
037%3758/94/$07X@ 0 1994-Elsevier SSDI 0378-3758(93)E0106-Q
Science B.V. All rights reserved.
P.O. Box
58
W.R. Myers et al. / Some alphabetic
optimal designsfor the logistic regression
model
authors as well as others developed numerous optimality criteria based on either parameter estimation or prediction of a response. From this evolved the term alphabetic optimality due to the use of alphabetic characters to categorize the various design criteria (e.g. A, D, G). Suppose we consider a general model of the type
(1.1) where f may in fact be nonlinear in /I and E has a known distribution, not necessarily normal. The vector p contains k unknown parameters and x contains regressor variables that can be controlled by the researcher. The precision of the maximum likelihood estimator of p is derived from the information matrix Z(B) which has (i,j) element
Here, of course, In L is the log likelihood. The design criterion that is most often discussed is D-optimality which is described by Max Det CIWNI , where the maximum is taken across designs and N is the design sample size. It is clear then that in the D-criterion the letter D stands for determinant. The motivation of D-optimality stems from the relationship between the size of the determinant of the information matrix and the quality of the maximum likelihood estimators of the model parameters. In fact, maximization of the determinant of the information matrix is equivalent to minimization of the volume of the lOO(1 -a)% asymptotic confidence ellipsoid on the model coefficients. Optimal-design theory has received some attention in applications to nonlinear modeling, with most of the focus confined to the D-criterion. But in the case of linear models with standard homogeneous variance assumptions on the error, the applications literature emphasizes not only D but several other criteria. Suppose in fact we consider Y=XB+&, (1.2) where the elements of E are i.i.d. N(0, a’). The X matrix represents a (N x k) general model matrix, where N is the number of experimental units and k is the number of model parameters. Y is a (N x 1) vector of observed responses and /3 is a (k x 1) vector of unknown parameters. The ordinary least-squares estimator of the parameters vector is /?=(X’X)-‘X’Y. Under normal theory, the maximum likelihood estimator (MLE) is equivalent to the ordinary least-squares estimator. It follows that where c2 is the unknown variance 1(/?)=X’X/a’ and indeed Varj?=o’(X’X)-‘, associated with each observation. Two points should be stressed in the standard linear model case: (1) 1(/I) is independent of model parameters and (2) the individual error
W.R. Myers et al. / Some alphabetic
variances problem
are homogeneous.
It is clear that these two conditions
in this case. For example, X’X ~ [ N
MaxDet
optimal designs for the logistic regression
the D-criterion
reduces
59
model
simplify
the design
to
1 ’
which is independent of model parameters. In the more general case depicited in (1.1) one is quite often confronted with an 1(p) that is not parameter free. Another criterion that has received attention in the linear models literature is Q-optimality, also known as IV optimality. IV refers to ‘integrated variance’. Here the average prediction variance over a specific design region, R, is minimized, i.e.,
sRVar Y(x)dx,
MinNF-l
-
2
an arbitrary where Var ?= a2x’(X’X)- ’ x, and F = &d x. The term x’ represents vector in the model space. Box and Draper (1959, 1963) first considered the idea of average prediction variance or integrated prediction variance over a defined region. Their work focused on an average weighted mean squared error, which is divided into the sum of the weighted variance and squared bias, averaged over the region R. As a result, a sharp distinction is drawn between D-optimality which focuses on parameter estimation and Q-optimality which emphasizes predicted values in a specific region of the factor space. A third criterion that will be discussed is G-optimality. The G-optimal design minimizes the maximum prediction variance over a specific design region R, i.e., Min
Var Chl
Max
XER Silvey (1980), Atkinson
ls2
.
(1982) and Box and Draper
discussion of these and other design criteria. It should be clear that, in addition to D-optimality,
(1987) provide
an informative
the Q- and G-criteria
are
independent of the model parameters in the case of the standard linear model. Further, as a result of the assumption of homogeneous variances, rr2 will have no effect on the development of the design. Consequently, X’X plays a predominant role in determining the optimal design. As a result, the design is solely dependent on the placement of the design points in the design space. Box and Lucas (1959) addressed the problem of experimental design in a general nonlinear model, f(x;/3). They were able to convey the complexity of obtaining a D-optimal design by pointing out that the information matrix is dependent on the model parameters. Therefore, if one is interested in obtaining a particular optimal design for fitting a nonlinear model, the parameters must be known. Only a small part of the literature involving optimal-design theory has been directed toward the logistic regression model. The logistic model is frequently used
W.R. Myers et al. / Some alphabetic
60
optimal designs for the logistic regression
model
when one works with binary data. It is expressed as 1
P=l+exp(-xx’P)
(1.2)
’
where p is the probability of a response at x. The vector x is a set of independent variables, and /? is a vector of unknown model parameters. Since the logistic model is nonlinear, the information matrix is a function of unknown model parameters. In the dose-response situation where binomial data are being observed, suppose for m groups, ri out of 12;subjects respond in the ith group. The log likelihood becomes the familiar binomial log likelihood given by lnL=
f In yi 0 I
+i$l~iln(Pi)t
:
i=l
@-ri)ln(1
-Pi).
i=l
The information matrix for a k-variable logistic model is expressed as
izlniPi(l
jlniPitl I(/4
-Pi)
-Pilxli
flniPi(l
-Pi)Xli
..’
2 i=l
-Pilxi
...
%Pi(l
i$lnipi(l
j14Pi(l
-Pijxki
-Pijxkixli
=
It is clear then that for the logistic model the determination of the optimal design for any of the criteria discussed requires the user to supply preliminary parameter estimates or guesses.
2. Overview of optimal designs for the logistic model Research in the area of experimental design for the logistic model has been primarily restricted to D-optimality. Simple interpretation and calculation are two reasons for its widespread use. White (1975) derived the D-optimal design for a quantal dose-response model with a symmetric tolerance distribution. Kalish and Rosenberger (1978) developed a two-level D-optimal design for the model parameterized as 1
‘=
1 +exp(-&+fi,x))’
W.R. Myers et al. 1 Some alphabetic optimal designs for the logistic regression model
These
authors
imposed
symmetry
value of x which produces tion matrix
a loop%
for the logistic
model
around response
the ED5,,.
The EDloop
rate. The determinant
in the case of a two point
represents
61
the
of the informa-
design
is expressed
as
Here p1 and p2 represent the probability of a response at the design points x1 and x2, respectively. As a result, the two-level D-optimal design is obtained by placing half the data at ED17.6 and the other half at EDs2.4. Abdelbasit and Plackett (1983) reviewed the two-level results and further pursued D-optimality for the logistic model by considering three-level designs. The three-level equal sample size D-optimal design, symmetric around the ED5,,, is formed by the and Plackett also investigated the doses ED13.6, ED50, and ED86.4. Abdelbasit robustness properties of these designs relative to poor ‘parameter guesses’. They also explored a multistage procedure for producing D-optimal designs in the one-variable logistic model case. Minkin (1987) then followed by showing that the previously created symmetric designs remain optimal over a general class of designs, eliminating the restriction of symmetry. Minkin further developed the multistage procedure of Abdelbasit and Plackett by considering conditional likelihood theory. In a recent book, Morgan (1992) gives a brief review of much of the work discussed here. Scientific investigators often try to characterize a pharmacological agent by determining its EDs0. Kalish (1990) developed optimal designs for estimating the EDs0, i.e., minimizing Var (ED,,), as well as the designs that minimize Var (jl) i.e., slope optimal. He performed an in depth comparison of these designs along with the D-optimal designs. Using the methods of Shenton and Bowman (1977), Kalish derived expressions for the variance of the MLEs of parameters in the logistic model, including terms to order 1/N2. He demonstrated that the ED,O-optimal design, using only the l/N term, requires all its data to be observed at the EDso. This degenerate design is impractical. Kalish further illustrates that for the ED,O-optimal design, a large discrepancy exists between the asymptotic results and the designs that are computed by including terms through order 1/N2. On the other hand, in the case of the slope-optimal and the D-optimal designs the asymptotic results are very good. This is especially true for the D-optimal design. Kalish demonstrated that the two-point ED5,-optimal design is not very efficient for global estimation (D-optimality). However, the D-optimal design is moderately efficient for estimation of the ED,o and highly efficient for estimating the slope. Kalish and Rosenberger (1978) developed a two-level G-optimal design. The criterion for their G-optimal design is Min Mey Var (ii),
62
W.R. Myers et aLlSome alphabetic
which places half of the observations Rosenberger interest R.
model.
probabilities
related works that encompass
One approach
model
at the EDT6,s and half at the ED23,2. Kalish and
(1978) define the entire range of response
There are other partially logistic
optimal designs for the logistic regression
is to use Bayesian
as the region
experimental
methodology
of
design for the
and employ
a prior
distribution on the unknown parameters. Chaloner and Larntz (1989) assumed that the ED50 and slope parameter have uniform and independent prior distributions on an interval. They indicate that this assumption may not be appropriate in all instances. The paper provides optimal designs given various prior distributions for both parameters. These distributions vary in the precision of information provided. Chaloner and Larntz observe that as the uncertainty in the prior distribution increases so do the number of design points. Zacks (1977) and Tsutakawa (1980) provide additional insight in the use of Bayesian techniques to arrive at experimental designs. Sequential designs have also been considered for estimating the percentiles of a quanta1 response curve. Wu (1985) developed an updating rule based on an efficient summary of all the data available via a parametric model, The procedure for the two parameter logistic model is analogous to the Robbins and Monro (1951) procedure in the case of binary data. Also presented is a nonparametric sequential design for estimating the median. Wu (1985) explains that the sequential nature of the design requires prompt responses so that the length of the experiment is not over extended. These sequential procedures will not apply to experimental situations where the time required for a response to occur is long. For further discussion on sequential designs see Wu and Wynn (1978), Wu (1985a) and Wu (1986). The final area of optimal-design theory that will be reviewed is the development of designs based on the sample size distribution assigned to preset values of the independent variable. In a clinical trial this would involve determining the number of patients to allocate to each treatment. Begg and Kalish (1984) provided D-optimal, Dsoptimal, which minimize the variance of the slope parameter, and G-optimal designs in the two-level
case. The three optimal
D-optimality The sample size distribution
designs
is independent
are described
in the following
way:
of p1 and p2, i.e., nl = n2.
D,-optimality nl CP2(1 -P2)11’2 ~=[p1(~-pI)]li2+[p2(~_p2)]li2~
Note that this design Bernoulli variance. G-optimality n1 77=~1(1-~1)+
allocates
more
P2(1 -P2) CPZ(~ --Pz)’
%=I-$. subjects
EL-$. N
to the treatment
with the smallest
63
W.R. Myers et al. 1 Some alphabetic optimal designs for the logistic regression model
Note the similarity, the D,-optimal the smallest
yet distinction,
between
the allocation
design. This design also allocates Bernoulli
3. Region-dependent
here and the allocation
more subjects
to the treatment
for with
variance.
optimal designs
The most common form of the logistic model is represented by equation (1.2). In the one-variable case, using the logit transformation, this model is expressed as logitp(xJ=bo+Brxi, Another
parameterization,
i=l,2
,..., N.
(3.1)
that will be used throughout
logitp(xi)=P,(xi-p)
this paper, is of the form
(i=l,2,...,N)
=Zi.
The terms b1 and p= -/IO/p1 correspond to the slope and the ED5,,, respectively. Simplicity is an attractive feature of this particular parameterization. In dealing with Q-optimal designs for logistic regression, it is important to use z, i.e., logit space to define the region of interest rather than the natural metric of x. This will be apparent later in this section. Myers (1991) provides further discussion of the different parameterizations of the logistic regression model. The major objective of this paper is the development of a new class of optimal designs based on the prediction of a response for the logistic regression model. This section will focus on creating Q-optimal designs for general regions and G-optimal designs for regions of interest centered at the ED5,, . The region of interest is the range of the independent variables over which the prediction variance will be averaged in the case of Q-optimality, or, over which the prediction variance will be maximized in the use of G-optimality. For example, in the case of a bioassay for a given compound, it could
be the range
of doses
that
produce
responses
in 10%
to 90%
of the
experimental subjects, or in logit space, from zr = - 2.150 to z2 = 2.150. The Q-optimality criterion addresses a design’s capability to have good average prediction in a specific region. The expression for the asymptotic average prediction variance
(APV) of logit ii across a region APV=NF-’
R is written
z’CZ(/?)*I-‘zdz,
as
(3.2)
where N is the total sample size, F is the volume of the region R, the region of interest in the logit space, and z’= [l z,], I(/?)* corresponds to the information matrix that is derivable from the two parameter fitted model logitlj(zi)=ft+fizi
(i=l,2,
. . ..N).
(3.3)
64
W.R. Myers et al. /Some alphabetic optimal designs
for the logistic regression
model
This conveniently allows the prediction variance to be a function of z. Here, of course, /?: and fi are merely estimating 0 and 1, respectively. For the symmetric design region R is [-a, a] and F = 2a. For an m-level design with ni runs at the ith level, equation (3.2) can be simplified to the form N
APV=Det[I(/?)*]
2 ‘(l-pi)Zf+~,~ [ i=ln’pC
nipi(lGpi)]. 1-t
(3.4)
It can be shown that the APV of logit $i is independent of the sample size and is represented by a single numeric value. The Q-optimal design is that which minimizes APV. One obvious advantage Q-optimality has over other design criteria, especially D-optimality, is its capability to consider various regions of interest. However, as we shall subsequently demonstrate, this also complicates its implementation. One should be reminded that D-optimal designs are invariant to region. A second design criterion that addresses prediction is G-optimality. Here, the G-optimal design is one which is associated with the minimum value of the maximum prediction variance (MPV) of logit $i over a predefined region of interest, i.e., MinMzr’[I(/?)*]-ls.
(3.5)
The only designs considered for G-optimality were symmetric around the ED5,,. For nonsymmetric designs no closed-form solution for the maximum prediction variance can be found. For the one-variable case, assuming a symmetric region of interest R, the maximum variance of logit fi design can be expressed as 1 MPV=Det[l(/?)*]
f ‘(l-pi)Zf+Ll’ i=l “‘I
i=l 1I f
nipi(l-Pi)
(3.6)
By observing (3.6) it is obvious that the MPV of logit $ is at the outer boundary of the design region. Therefore, the G-optimal design minimizes the prediction variance at the outer boundary of a defined region. It is interesting to note from (3.6) and (3.4) that G-optimality for the region (--a, a) is equivalent to Q-optimality for the region (-a&S). Two- and three-level Q- and G-optimal designs were constructed for a selected number of regions of interest. The three-level designs assume equal sample size at each design point and assume the EDso as the middle design level. There is no claim that the three-level designs chosen are optimal three-level designs. Indeed the optimal three-level design would assign no design runs to the center and reduce to the two-level optimal design. Fewer than N/3 runs in the center would result in the outer levels moving closer to the levels given in the two-level design. Figure 1 displays the results for the two- and three-level optimal equal sample size designs. The figure illustrates how the design points adjust depending on the region of interest. As the region begins to contract, the outer design points tend to gravitate toward the design center. For example, in the two-level case, for the region from ED01 to ED99 the two
65
W.R. Myers et al. / Some alphabetic optimal designs for the logistic regression model
0.05
0.01
I
I
I
0.05
0.10
0.15
1
Lower Perimeter
‘----_________-.
-.-._.-._.-.-~-
---------j-Levels
Fig. 1. Optimal
1
0.20
0.25
I
I
1
0.30
0.35
0.40
r 0.45
of Design Region 2-Levels
Q
Z-Levels
G
j-Levels
Q
G
designs
for symmetric
regions.
design points for the Q-optimal design are ED13.6 and ED86,4, while for the region from EDjo to EDTo they are ED27.9 and EDT2.r. It is interesting that for the wider regions the design points are within the region of interest while for narrow regions the the region is contained within the design points. The designs were obtained by using the Nelder and Mead (1965) simplex search procedure. Another interesting observation is the symmetric nature of each design. Although there is no restriction on the Q-optimal design to force symmetry, the designs are symmetric as long as the region is symmetric about the EDSo. In comparing two- and three-level Q-optimal designs note that the outside points for the three-level designs are beyond their two-level counterparts. The G-optimal designs illustrate similar qualities to those found in the Q-optimal designs. One obvious difference between the two optimality criteria is that the points for the G-optimal design are further from the design center. This is due to the requirement that they predict well at the perimeter of the experimental region. The same methodology used to produce the Q-optimal designs was used for the G-optimal designs.
W.R. Myers et al. / Some alphabetic optimal designs for the logistic regression model
66
Now that these particular optimal designs have been created it seems reasonable to make comparisons between designs. For example, the Q-efficiency of design A relative to design B within region R is defined as follows: Q-efficiency =
APV of design B for region R APV of design A for region R’
(3.7)
Through the Q-efficiency one can determine how well two- and/or three-level D-optimal designs are able to predict relative to the Q-optimal designs within various regions. One must keep in mind that the D-optimal design is invariant to the region of interest, thus intuitively it would seem that neither the two- nor three-level D-optimal designs would be effective in the prediction over all types of regions. Table 1 provides the Q-efficiency of the two- and three-level D-optimal designs for selected regions as well as a Q-efficiency for the three-level equal sample size Q-optimal designs. The results indicate that when the region is narrow, the D-optimal design does not compare favorably to the Q-optimal design. For example, when the region of interest is from ED30 to EDT,,, the Q-efficiencies of the two- and three-level D-optimal designs are only 82.9% and 89.7%, respectively. An efficiency can also be constructed for the G-optimality criterion, by replacing the APV with the MPV. Table 2 lists the Table 1 Q-efficiencies Region
relative to the two-level
of interest
Two-level
Q-optimal D-optimal
design Three-level
D-optima1
Three-level
(%)
(%)
(%)
EP3-ED99
95.5
ED,,-ED,, EDICED,, ED, s-ED*5 ED*,-ED,, ED,,-ED,, EDFED,,
99.8 99.3 96.6 92.8 88.2 82.9
80.2 89.7 94.0 95.3 94.7 92.4 88.1
83.3 90.4 94.0 96.1 97.4 98.4 99.0
Table 2 G-efficiencies Region
relative
of interest
ED,,-ED99 EDos-ED95 EDI,.-ED90 EDIs-ED~s ED,,-EDso EDx-ED,5 EDx-ED,,,
to the two-level
G-optima1
Two-level
D-optimal
Q-optimal
design
(%)
(%)
Three-level
D-optimal
Three-level symmetric G-optima1 (%)
87.6 94.0 97.9 99.8 99.8 97.9 93.9
70.1 77.9 84.4 89.3 93.0 95.1 94.9
74.9 81.5 86.4 90.1 93.1 95.4 97.1
61
W.R. Myers et al. / Some alphabetic optimal designs for the logistic regression model
G-efficiencies
of the two- and three-level
of the three-level
D-optimal
equal sample size Q-optimal
fare better in terms of G-efficiency the larger regions. A further study involves
designs as well as the G-efficiency
design. The D-optimal
for the more narrow
the relative
comparisons
regions
among
designs tend to
and Q-efficiency
the newly created
for two-
and three-level Q-optimal designs. It is apparent that the two-level designs are more effective in average prediction. For instance, the Q-efficiencies of the three-level designs relative to the two-level designs range from 83.3% to 99.0%. However, a three-level design will be considered more practical in the opinion of most investigators because of allowing a test for designs relative to G-efficiencies range
the importance of providing more than two design levels, thereby model adequacy. The G-efficiencies of the three-level G-optimal the two-level G-optimal designs are also shown in Table 2. The from 74.9% to 97.1%.
4. Robustness In order to construct
each of the optimal
designs that have been presented
in this
paper, the investigator must have prior knowledge of the true model parameters or the Because instances occur when limited knowledge of the true values of the ED,,,,,. parameters is available, it is necessary to determine the efficiency of these designs when one must guess the values of the independent variable to be used. In other words, how robust are the designs to poor guesses of the levels? Designs that are shown to be robust to inadequate choices of the optimal dose levels will obviously be more valuable to the investigator. This section will only consider the robustness of Q-optimal designs. Figure 2 provides contour plots to represent the Q-efficiency for all levels in the design region. The three-level
designs are not investigated
because
complexity.
of their dimensional
However,
by these graphical
they will be examined
methods via effici-
ency tables. The region of interest used for the purpose of illustration is EDlo to EDgo. Figure 2 shows, for example, that if an investigator over-estimates the optimal values (ED,,,,,ED,,,,) by choosing say ED,,,, and EDg2, the Q-efficiency is still about 90%. Many other poor guesses producing an efficiency over 90% are certainly apparent. By observing the contours which represent the Q-efficiencies, it becomes obvious that there is considerable room for error in the selection of the dose levels for producing a design with a high Q-efficiency. Other regions of interest were investigated and they all indicated good overall robustness with respect to poor choices of the dose levels. It is reassuring to note that the designs are quite robust in light of the fact that in practice one must always rely on preliminary data to arrive at guesses for the doses to be used in the experiment. Efficiency tables are presented in Table 3 for the region ED,,-ED,,. The objective of the Q-efficiency tables is twofold. First, they provide an understanding of how
W.R. Myers et al. 1 Some alphabetic optimal designs for the logistic regression model
The Probability of a Response at the Low Dose Level Fig. 2. Q-efficiency
Table 3 Comparing
the Q-efficiencies
relative
to design levels region of interest
of two- and three-level
Q-optimal
is EDlO-ED90.
designs from ED,,
to ED,,
B,(P-Po)/BI/BIo
0.85
0.90
0.95
1.00
1.05
1.10
1.15
Two-level 1.5 1.0 0.5 0.0
0.5168 0.7533 0.9151 0.9689
0.5579 0.7900 0.9389 0.9866
0.5965 0.8205 0.9553 0.9968
0.6316 0.8443 0.9644 1.00
0.6626 0.8614 0.9668 0.9969
0.6889 0.8719 0.9632 0.9884
0.7104 0.8762 0.9541 0.9749
0.5569 0.7464 0.8742 0.9182
0.5980 0.7773 0.8929 0.9311
0.6358 0.8025 0.9046 0.9382
0.6695 0.8220 0.9115 0.9404
0.6989 0.8364 0.9138 0.9384
0.7238 0.8458 0.9121 0.9329
0.7441 0.8509 0.907 1 0.9245
Three-level 1.5 1.0 0.5 0.0
design
design
W.R. Myers et al. / Some alphabetic
robust
the designs
are with respect
they permit a comparison
optimal designs for the logistic regression
to choosing
of the robustness
the parameters
of two- and three-level
model
69
fll and p. Secondly, designs. The tables
present the Q-efficiency of a design based on the parameters p1 and p. The terms fil and p are the true slope and ED,o, while PI,, and p. are the values that correspond to the preliminary parameter estimates. The Q-efficiency of a given design is relative to the optimal design (the two-level design) constructed in Section 3. The general nature of the results is the same for all regions. It is shown in the appendix how the APV of logit pi can be written in terms of the quantities /?1//?1o and /?l(p--o). This provides an ease in tabling Q-efficiencies. A summary of the table is as follows: (1) Both the two- and three-level designs appear to be fairly robust to poor initial guesses of parameters. (2) For both two- and three-level designs underestimating the true slope (pl/plo > 1) results in a higher Q-efficiency than does overestimating it. 1.0)) the effect of a poor (3) As the estimate of the EDso improves (fll(p-po)< estimate for the slope is not as severe. Some of the conclusions given above regarding design robustness are intuitively reasonable. For example, in the case of (2), an underestimate of the slope implies that the chosen design has a wider spread than a case in which the slope is overestimated, and it seems reasonable that a mistake of the latter variety would result in a greater reduction in efficiency. Conclusion (3) is certainly reasonable due to the importance of achieving symmetry around the ED50 in the design. The larger error in the preliminary estimate
of EDso the greater
the loss in symmetry
created
by an error in the
slope. Robustness studies were made for other regions but are not displayed here. For the three-level designs, the wider regions of interest seem to be more robust to poor guesses than the smaller ones. The same result is true for the two-level designs, only to a lesser degree. See Myers (1991). This makes good intuitive sense. Specific errors in the slope or the ED 5o for narrow regions results in more serious deviations regions.
from symmetry
and more serious deviations
from target ED’s than in wider
5. Designs for regions not centered at ED9 Each of the optimal designs that have been developed up to this point in the paper have involved a region of interest that is symmetric about p= EDso. However, instances will occur when an investigator is interested in a small, more restricted region. For example, a pharmacologist may be concerned with predicting well over the range of doses associated with the response rates from 10% to 30% or possibly from 20% to 40%. This section will focus on Q-optimal designs for regions that are not symmetric around the EDso. The expression for the APV of logit fii across a region not centered at p is similar to equation (3.4) with the only difference being the limits of integration. For the
70
W.R. Myers et al. / Some alphabetic
optimal designs for the logistic regression
model
one-variable case, the APV of logit $i for an m-level design is expressed as APV =
N
b-allDetCI(P)*I
where R= [aI, u2)]. Equation (5.1) is equivalent to (3.4) when a, = -u2, a region symmetric around z =O. While studying the “asymmetric” regions, one may focus, without loss of generality, on just one end of the logistic response curve. This implies that for our bioassay example, studying the design region from EDzo to ED4,, is identical to studying the design region from ED 60 to EDso. This is a result of the symmetric nature of the logistic function. 5.1. Two- and three-level designs with equal sample size
Two- and three-level Q-optimal designs with equal sample sizes are given for various design regions in Table 4. The particular regions chosen provide information that will be useful in understanding the type of optimal designs obtained. One should
Table 4 Equal sample size Q-optimal
designs for regions
Region of interest
Two-level
E&o-E&O
EDs.9 ED,,.9
EDrED,,
ED,,., ED,,.,
ED,,-ED,,
ED,., E&m
E&o-EDso
ED,., ED,,.,
ED,,-ED.,,
ED,,., ED,,.,
ED,,-ED.,,
E&.I ED 34.7
EDos-EDm
EDs.6 ED47.6
Q-optimal
not centered Three-level
ED,., E&s EDx.9 ED,,., Eb7.1 ED 51.4 ED,., EDs.4 ED 31.4 ED,,., ED,,., ED 38.5 ED,,., EDlg.3 ED52.1 ED,,., ED,,., EDw., ED,,., ED,,., ED 87.0
at EDsO Q-optimal
W.R. Myers et al. / Some alphabetic optimal designs for the logistic regression model
71
notice how the design points adjust depending upon the region of interest. Note how the two-level Q-optimal designs tend to concentrate observations in or around the region of interest. When examining the table one should notice the layout of the three-level Q-optimal design. It consists of a two-level design with two-thirds of the data at the low dose level and one-third at the high dose level. In other words, two of the design points are placed at identical dose levels. These particular designs produce a smaller APV of logit $i than the two-level Q-optimal designs. This suggests that a design will have a higher Q-efficiency if more observations are taken at the dose level with the smallest value of the Bernoulli variance [pi(l -Pi)]. 5.2. Optimal designs with unequal sample size The results obtained for the three-level designs with equal sample size suggest that some advantage may be gained by searching for the two- (and perhaps three-) level designs with unequal sample size. It is obvious that region coverage is important but the flexibility produced by unequal weighting is crucial when the region is not centered at the ED5,,. As a result, the Nelder-Mead procedure was used to construct first the optimal two-point design with the sample size allowed to be unequal. Table 5 shows these designs. Note that, as expected, the weighting is unequal at the two levels. Note also the general nature of the design. The left-hand point is placed within or close to the left boundary of the region while the right-hand point is placed, without exception, at the symmetric counterpart and with considerably smaller weighting. Thus, even though the region is not centered at z=O, the two design points do have
Table 5. Q-optimal designs regions not centered Region
of interest
(unequal at ED,, Dosage
sample
and weighting
W) EDICED EDO,CED,, ED,,-ED,, ED,,-ED,, ED*,-ED,, ED,,-ED,, ED,,-ED,,
size)
ED 17.4 (89) ED,,., (11) ED,, (88) ED,, (12) ED,,,, (82) ED,,., (18) ED,,., (79) ED,,., (21) ED,, (85) ED,, (15) ED,, (94) ED,, (6) ED,, (88) ED,, (12)
for
72
W.R. Myers et al. / Some alphabetic optimal designs for the logistic regression
model
the symmetry property, but with unequal weighting. The attempt to construct threelevel unequal sample size designs merely resulted in producing these same two-level designs. There certainly is a natural attraction to the two-level equal sample size designs in Table 4 since points are placed at or near the actual boundary of the region. It turns out that in some cases the equal sample size designs are reasonably efficient when compared to their unequal sample size counterpart. The efficiency of the equal sample size design is given and the Q-efficiencies of the two-level D-optimal designs are listed for the various regions in Table 6. The poor efficiency of the D-optimal designs is not surprising since the underlying strategy of D-optimality is to ‘cover’ the entire dose response curve rather than focusing on particular subregions. The equal sample size design is quite efficient when the region of interest is reasonably narrow and not centered close to either end of the dose-response curve. For example, the efficiency for a region ED 20-ED4,, is considerably better than the efficiency for the region ED ,,-EDZ5. Though the regions are of the same width, the 05-25 region contains a larger discrepancy in Bernoulli variance at the two end points, resulting in a lack of efficiency in the equal sample size design. This is quite reasonable intuitively since near equal Bernoulli variance at the end points emulates the ideal homogeneous variance structure. A rough rule is that an efficiency that exceeds 90% can be achieved when the ratio of the smallest to the largest Bernoulli variance at the end points exceeds about 0.6. With this type of region a reasonable prescription for the equal sample size design is to allot equal weighting to two design points that barely contain the region of interest. If the ratio of the smallest to the largest Bernoulli variance is smaller than 0.6, one should resort to a compromise given in the following section. 5.3. Compromise
designs
Even in ideal circumstances when the researcher has a reasonably good prior knowledge of the activity in the dose-response region, designs for logistic regression Table 6 Q-efficiencies
for regions
Region of interest
EDlo-E&o EDox-EDz, EDlo-ED20 ED,,-ED,, ED,,-ED,, ED,,-ED,, EDAo-ED40
not centered
at EDSo
D-optimal
Equal sample size Q-optimal
Design Cl
WI
WI
W)
62.4 59.4 56.79 70.1 74.9 65.3 61.3
85.9 71.6 89.1 71.9 74.7 95.6 99.1
95.6 91.9 94.4 80.1 59.7 88.1 88.5
W.R. Myers et al. / Some alphabetic optimal designs for the logistic regression
model
13
are difficult to implement due to the requirement of ‘guesses’ of parameters. Difficulties lie in the fact that the correct weighting depends on the region of interest. This suggests the consideration of ‘compromise designs’ that are not optimal but somewhat easier to implement. Let us first consider an attractive alternative to the optimal two-level designs listed in Table 5. In order to simplify implementation and not deal with weighting being a function of region, a study was conducted to produce a simpler but fairly efficient prescription. As a result, we offer a recommendation that places 85% of the weighting at the center of the region of interest and 15% of the weighting on the symmetric compliment. The efficiency of the design is quite good, particularly for narrow regions. Call this compromise design Cl. For some regions success will be achieved by placing the 85% weighting on a level a bit closer to the end point that is opposite the EDOS-ED4,, and EDo5-EDSo, efficiencies improve dramatically by using the levels that have slightly more spread, e.g., 85% at EDi and 15% at EDsS. A more complete description of recommendations appear in the following section.
6. Summary and recommendations This paper has several objectives. While optimal designs under the standard linear model can easily be implemented by the placement of the design points in the region, the implementation of optimal design theory for the logistic model is much more difficult. The information matrix for the logistic model is a function of unknown model parameters. We have given a brief review of the research that had been focused in this area. This involved the development of D-optimal designs as well as other related work. We also have developed a class of designs which focused on providing good prediction or estimation of logit values over predetermined regions of the logistic regression model. The capability to adjust depending on the region of interest illustrates a major benefit of these designs. The comparison between the D-optimal design and the Q-optimal design demonstrated that the D-optimal design does not perform well in terms of the average prediction in many design regions. Within this new class of designs a subgroup was developed to handle regions of interest that are not centered at the ED=,,,. This provides the user further flexibility when considering a particular region of interest. It must be emphasized that like D-optimal designs the actual Q-optimal design cannot be implemented because of the knowledge required of the researcher in choosing the level of x that produces the desired ED1OOp. For the case where the researcher is interested in regions that are centered at the ED5,,, the two-level D-optimal design is very Q-efficient for ‘wider’ regions (ED,,-EDg9 to EDlS-EDS5). On the other hand, for more ‘narrow’ regions of interest that are centered at ED5,,, a good rule of thumb is to use an equal sample size design with dosages placed at the end points of the region or slightly containing the region. There are situations in
W.R. Myers et al. / Some alphabetic
74
optimal designs for the logistic regression
model
which three-level designs are practical for the user. Again, for wider regions centered at EDs0, the three-level D-optimal design is quite reasonable. For a narrow region described by (EDiecp, EDNCW -J, with p roughly exceeding 0.15, a reasonable prescription for a three-level equal sample size design is to use (ED1,,0p--5, EDso, EDoooo -_p)+5) ). In other words, the three-level design performs best for narrow regions when the region is slightly contained by the design. Higher efficiency is obtained for the three-level design by using fewer runs at the EDso and levels that are slightly closer to the region end points. In the case of regions not centered at ED50, implementation may be more difficult but more important since the D-optimal design is very inefficient. A compromise design was given in Section 5.3. Design Cl can be used in cases where the researcher has sufficient prior knowledge to be able to guess levels to achieve EDloop values well outside the experimental region. A slight adjustment in implementation of Cl will improve results. If the lower end point of the region of interest exceeds the EDz5 or the upper end point exceeds the ED,,, the user should attempt to create somewhat more spread in the design than that suggested by Cl. For example, for the region ED30-ED50 an efficient strategy is to use 85% weighting at about ED33 and 15% at EDe7. If the researcher is uncomfortable with guessing levels which are far outside the region of interest he may be led to an equal sample size design as illustrated in Table 4. An efficient design can be obtained when the ratio of the smallest to the largest Bernoulli variance at the region end points exceeds roughly 0.6. The design prescription is to use equal sample size at two design points that slightly contain the region of interest.
Appendix
Illustrating how the APV of logitfii can be written as a function of pl/plo and Pr(p--cIo). From expression (3.4) the APV of logit$i in the two-level case is 2 PIU
-P1)P2(1 x
PIU
-Pz)CZ1
-%I2
-PI)z:+P*(l
-PJz:+;(PIu
--P1)+P2(l-P*))
1 (Al)
assuming n, = n2. Recall from Section 3 that for a two-level design z1 = /I1(x1 - p) and z2 = p1 (x2 - p), were zi = In [pi/(1 -pi)]. Preliminary parameter estimates p. and fir0 lead to the design points x1
1 =po--In P 10
Pzo
__ ( 1 -Pzo )
W.R. Myers et al. 1 Some alphabetic optimal designs for the logistic regression model
15
and x2=po+-In
1 B 10
The quantities
-
P20
( l-P20 ) .
pi0 and p20 are the probabilities
of a response
for the low and high dose
level, respectively, based on the Q-optimality criterion. It should also be pointed out that pi0 = 1 -pzo. By substituting the previous expressions for xi and x2 one can rewrite the various terms as follows: 1 ‘l=l
+exp(-Bl(~0-(1/P,o)ln(p201p,o)-~))
1 = 1 +exp(B1(~--o)(p201p,o)(81ia’0’)
’
1 P2=1+exp(-_~(~0+(1/B,0)ln(p201~~0)--))
1 =l +exp(P,(~++0)(p201p,0)'P"B'"')
=PTC-2(1/B~0)ln(p201Plo)12 = C-2(P1/Blo)ln(~20/~10)lz.
d=8~cx1-P12
By substituting
the above results into (Al) the APV of logit fii for the two-level
case
76
W.R. Myers et al. /Some alphabetic optimal designs for the logistic regression model
can be rewritten
as
642)
W.R. Myers et al. / Some alphabetic
Therefore,
given a region of interest
based on the quantities design relative
to another
optimal designs for the logistic regression R [-a,
a]
B1/p10 and pl(p-pO). can also be calculated
77
model
the APV of logit ji can be computed In addition,
the Q-efficiency
of one
based on these terms.
References Abdelbasit, K.M. and R.L. Plackett (1983). Experimental design for binary data. J. Am. Statist. Assoc. 78, 90-98. Begg, C.B. and L.A. Kalish (1984). Treatment allocation for nonlinear models in clinical trials: the logistic model. Biometrics 40, 409-420. Box, G.E.P. and N.R. Draper (1959). A basis for the selection of a response surface design. J. Am. Statist. Assoc. 54, 622-654. Box, G.E.P. and N.R. Draper (1963). The choice of a second order rotatable design. Biometrika 50,335-352. Box, G.E.P. and N.R. Draper (1987). Empirical Model-Building and Response Surfaces. Wiley, New York. Box, G.E.P. and H.L. Lucas (1959). Design of experiments in non-linear situations. Biometrika 46, 77-90. Chaloner, K. and K. Larntz (1989). Optimal Bayesian design applied to logistic regression experiments. J. Statist. Planning and Inference 21, 191-208. Kalish, L.A. (1990). Efficient design for estimation of median lethal dose and quanta1 dose-response curves. Biometrics 46, 737-748. Kalish, L.A. and J.L. Rosenberger (1978). Optimal designs for the estimations of the logistic function. Technical Report 33. Pennsylvania State University. Kiefer, J. and J. Wolfowitz (1959). Optimum designs in regression problems. The Ann. Math. Statist. 30, 271-294. Minkin, S. (1987). Optimal designs for binary data. J. Am. Statist. Assoc. 82, 1098-l 103. Morgan, B.J.T. (1992). Analysis of Quanta1 Response Data. Chapman & Hall, London. Myers, W.R. (1991). Optimal experimental designs for logistic regression. Ph.D. dissertation. Department of Biostatistics, Virginia Commonwealth University, Richmond, Va. Nelder, J.A. and R. Mead (1965). A simplex method for function minimization. Comput. J. 7, 308-313. Robbins, H. and S. Monro (1951). A stochastic approximation method. Ann. Math. Statist. 29, 400-407. Shenton, L.R. and K.O. Bowman (1977). Maximum Likelihood Estimation in Small Samples. Griffin, London. Silvey, S.D. (1980). Optimal Designs. Chapman & Hall, London. Tsutakawa, R.K. (1980). Selection of dose levels for estimating a percentage point of a logistic quanta1 response curve. Appl. Statist. 29, 25-33. White, L.V. (1975). The optimal design of experiments for estimation in nonlinear models. Unpublished Ph.D. thesis. Department of Mathematics, Imperial College, London. Wu, C.F.J. (1985a). Asymptotic inference from sequential design in a nonlinear situation. Biometrika 72, 553-558 Wu, C.F.J. (1985b). Efficient sequential designs with binary data. J. Am. Statist. Assoc. 80, 974-984. Zacks, S. (1977). Problems and approaches in design of experiments for estimation and testing in nonlinear models. In: P.R. Krishnaiah, Ed., Proc. 4th Int. Symp. on Multivariate Analysis. North-Holland, Amsterdam, 209-223.