TESTING LACK OF FIT IN A MIXTURE MODEL
BY
JOHN THOMAS SHELTON
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1982
To Nydra
and
My Parents
ACKNOWLEDGEMENTS
I would like to express my deepest appreciation to Drs.
Andre' Khuri and John Cornell for suggesting this topic to
me and for providing constant guidance and assistance. They
have made this project not only a rewarding educational
experience but an enjoyable one as well. A special word of
thanks goes to Mrs. Carol Rozear for her diligence in
transforming my handwritten draft into an expertly typed
manuscript.
iii
TABLE OF CONTENTS
page
ACKNOWLEDGEMENTS............ ............. .. .............iii
ABSTRACT .............................................. vii
CHAPTER
ONE INTRODUCTION................................... 1
1.1 The Response Surface Problem............... 1
1.2 The Mixture Problem....................... 5
1.2.1 Mixture Models..................... 6
1.2.2 Experimental Designs for Mixtures.. 12
1.3 The Purpose of this Research:
Investigation of Procedures for Testing
a Model Fitted in a Mixture System for
Lack of Fit ............................... 17
TWO LITERATURE REVIEWTESTING FOR LACK OF FIT..... 19
2.1 Introduction............................... 19
2.2 Partitioning the Residual Sum of Squares.. 21
2.3 Testing for Lack of Fit Without
Replicated ObservationsNear Neighbor
Procedures ................................ 26
2.4 Testing for Lack of Fit with Check Points. 33
THREE AN OPTIMAL CHECK POINT METHOD FOR TESTING
LACK OF FIT IN A MIXTURE MODEL................. 40
3.1 Introduction............................... 40
3.2 Testing for Lack of Fit in the Presence
of an External Estimate of Experimental
Error Variation............................ 41
3.2.1 The Test Statistic................. 41
3.2.2 The Testing Procedure and an
Expression for the Power of
The Test........................... 45
3.2.3 A Method for Locating Optimal
Check Points ....................... 47
3.3 Testing for Lack of Fit When MSE Is
Used to Estimate Experimental Error
Variation..................................... 51
3.3.1 The Test Statistic................. 51
3.3.2 The Rejection Region and its
Relation to the Power of the Test.. 53
iv
3.3.3 A Method for Locating Optimal
Check Points ....................... 56
3.3.4 Determining Whether the Test Is
Upper Tailed or Lower Tailed....... 58
3.4 Examples. ................................ 67
3.4.1 Theoretical Examples............... 67
3.4.2 Numerical Examples................. 83
3.5 Discussion. ............................... 95
FOUR USE OF NEAR NEIGHBOR OBSERVATIONS FOR
TESTING LACK OF FIT .............................. 99
4.1 Introduction ............................. 99
4.2 Notation......... ..........................101
4.3 Shillington's Procedure ...................106
4.3.1 Development of MSEB ................109
4.3.2 Development of MSEW................110
4.4 Development of SSE(weighted) .............112
4.5 Equivalence of SSEW and SSEW(weighted)....116
4.6 The Test Statistic .........................118
4.7 The Testing Procedure and its Power........122
4.8 Selection of Near Neighbor Groupings....... 125
4.8.1 Example 1Stack Loss Data.........129
4.8.2 Example 2Glass Leaching Data.....134
4.9 Discussion................................ 142
FIVE CONCLUSIONS AND RECOMMENDATIONS ................145
APPENDICES
1 INFLUENCE OF X1 ON
P{F" > F } ...................156
v 1'2 ;11 '12 a;v ,v2
2 A CONTROLLED RANDOM SEARCH PROCEDURE FOR
GLOBAL OPTIMIZATION ............................159
3 STATISTICAL INDEPENDENCE OF d'V d/o2
2 0
AND SSE/o .................................... 164
4 THEOREM 3.1................................... 168
5 THEOREM 3.2...................................... 169
6 AN APPROXIMATION TO THE DOUBLY NONCENTRAL
F DISTRIBUTION.................................. 171
7 EQUIVALENCE OF SSEB AND SSLOF WHEN
REPLICATES REPLACE NEAR NEIGHBOR
OBSERVATIONS...................................... 172
8 LEMMA 4.1...................................... 175
9 PROOF OF THEOREM 4.1(i) ............ ......... 178
10 PROOF OF THEOREM 4.1(ii)......................182
11 PROOF OF THEOREM 4.1(iii).............. ........185
12 PROOF OF THEOREM 4.2.............................191
1
13 PROOF OF THE EQUALITY SSE = d'V d + SSE .....193
A 0 
REFERENCES .... ....................................... 198
BIOGRAPHICAL SKETCH.. .......... .. ..... ............202
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in
Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy
TESTING LACK OF FIT IN A MIXTURE MODEL
By
John Thomas Shelton
May 1982
Chairman: Andre' I. Khuri
Cochairman: John A. Cornell
Major Department: Statistics
A common problem in modeling the response surface in
most systems, and in particular in a mixture system, is that
of detecting lack of fit, or inadequancy, of a fitted model
of the form E(Y) = X 1 in comparison to a model of the form
E(Y) = X1 + X 22 postulated as the true model. One method
for detecting lack of fit involves comparing the value of
the response observed at certain locations in the factor
space, called "check points," with the value of the response
that the fitted model predicts at these same check points.
The observations at the check points are used only for
testing lack of fit and are not used in fitting the model.
It is shown that under the usual assumptions of
independent and normally distributed errors, the lack of fit
test statistic which uses the data at the check points is an
vii
F statistic. When no lack of fit is present the statistic
possesses a central F distribution, but in general, in the
presence of lack of fit, the statistic possesses a doubly
noncentral F distribution. The power of this F test depends
on the location of the check points in the factor space
through its noncentrality parameters. A method of selecting
check points that maximize the power of the test for lack of
fit through their influence on the numerator noncentrality
parameter is developed.
A second method for detecting lack of fit relies on
replicated response observations. The residual sum of
squares from the fitted model is partitioned into a pure
error variation component and into a lack of fit variation
component. Lack of fit is detected if the lack of fit
variation is large in comparison to the pure error
variation. This method can be generalized when "near
neighbor" observations must be substituted for replicates.
In this case, the test statistic (assuming independent and
normally distributed errors) has a central F distribution
when the fitted model is adequate and a doubly noncentral F
distribution under lack of fit. The arrangement of near
neighbors is seen to affect the testing procedure and its
power.
viii
CHAPTER ONE
INTRODUCTION
1.1 The Response Surface Problem
A mixture problem is a special type of a response
surface problem. First we shall define the general response
surface problem and indicate the basic objectives sought in
its analysis, and follow this development with a discussion
of the mixture problem.
In the general response surface problem, we are inter
ested in studying the relationship between an observable
response, Y, and a set of q independent variables or
factors, xl, x2, ..., Xq, whose levels are assumed con
trolled by the experimenter. The independent variables are
quantitative and continuous. We express this relationship
in terms of a continuous response function, p, as
Y = (ul, x2 uq ) + u
where Yu is the uth of N observations of the response col
lected in an experiment, and xui represents the uth level of
the ith independent variable, u = 1, 2, ..., N, i = 1, 2,
..., q. The exact functional relationship, p, is unknown.
The term Eu is the experimental error of the uth
observation. It is assumed that E(cu) = 0, E(Eueu,) = 0,
for u u', and E(e ) = 2, for u = 1, 2, ..., N.
As the form of ( is unknown and may be quite complex, a
low order polynomial (usually first or second order) in the
independent variables xl, x2, ..., Xq is generally used to
approximate p. This may be justified by noting that such
polynomials constitute low order terms of a Taylor series
expansion of < about the point xl = x2 = ... = xq = 0,
(Myers, 1971, p. 62). Cochran and Cox (1957, p. 336) point
out that these low order polynomials may give a poor approx
imation to q when extrapolated beyond the experimental
region, and thus should not be used for this purpose.
A linear response surface model may be written in
matrix notation as
Y = X4 + E (1.1)
where Y is an Nxl vector of observable response values, X is
an Nxp matrix of known constants, a is a pxl vector of
unknown parameters (regression coefficients), and a is the
Nxl vector of random errors. When the model is a first or a
second degree polynomial, the columns of X correspond to the
first or second degree powers of the independent variables
xl, x2, ..., Xq, or their cross products. If the model
contains a constant term, 80, the first column of X will
correspond to this term, and will consist of N ones. Since
E(c) = 0, an alternative representation for the response
surface model of (1.1) is
E(Y) = X .
Once the form of the model that will be used to approx
imate ((xl, x2, ..., Xq) is chosen, the next step is to
estimate the regression coefficients, a, and then use the
estimated model to make inferences about the true response
function, (. The estimation of the elements of a is usually
accomplished by ordinary least squares techniques. For the
purpose of testing hypotheses concerning the regression
coefficients, ., it is assumed that L has a normal distribu
tion, that is, e ~ NN(O, a 2N).
Perhaps the most common objective in the exploration of
a response system is the determination of its optimum
operating conditions. By this we mean that it is desired to
find the settings of xl, x2, ..., xq that optimize (, which
in some applications may be interpreted as maximizing p,
while in other applications a minimum value of ( may be of
interest. It is also often desirable to determine the be
havior of the response function in the neighborhood of the
optimum. For second order response models, such an investi
gation can be carried out by performing a canonical analysis
of the second order surface as discussed in Myers (1971).
For simple systems having only one or two independent
variables, the response surface may be explored by just
plotting the fitted response values against values taken by
the independent variables. If q = 1, implying only one
independent variable, say xl, then a twodimensional plot of
the fitted response against xl may be used to locate the
optimum, as well as to investigate the response behavior in
other parts of the experimental range of xl. If q = 2, and
the two independent variables are xl and x2, then a plot of
the contours of constant response over the region specified
by the ranges of the values for xI and x2 can be used to
describe the response surface.
The properties that the fitted model possesses in terms
of its ability to represent the true surface, p, depend on
the settings of xl, x2, ..., Xq at which values of Y are
observed. Thus the experimental design is of great impor
tance. Much work has been done on the construction of
designs that are optimal with respect to one criterion or
another involving the fitted response and/or the true unfit
ted model. Box and Draper (1975) list fourteen criteria to
consider when choosing a design for investigating response
surfaces. Myers (1971) gives several designs for fitting
first and second order polynomial models. A discussion of
specific design considerations will not be attempted here,
as such a discussion is not the focus of this dissertation,
and would necessarily be lengthy.
The initial steps in the analysis of a response system
may be described as follows: First an attempt is made to
approximate the true response function, p(xl, x2, ..., Xq),
usually with a low order polynomial in xl, x2, ..., xq.
After the form of the model has been chosen, then comes the
selection of an appropriate experimental design, which
specifies the settings of the independent variables at which
observed values of the response will be collected. The
observed values of the response are used in estimating the
regression coefficients in the model, using, in general,
ordinary least squares. After a test for "goodness of fit"
of the model verifies the fitted model is adequate, the
fitted model is used in determining optimum operating condi
tions for the response system.
1.2 The Mixture Problem
A mixture system is a response system in which the
response depends only on the relative proportions of the
components or ingredients present in a mixture, and not on
the total amount of the mixture. For example, the response
might be the octane rating of a blend of gasolines where the
rating is a function only of the relative percentages of the
gasoline types present in the blend. The proportion of each
ingredient in the mixture, denoted by xi, must lie between
zero and unity, i = 1, 2, ..., q. The sum of the propor
tions of all the components will equal unity, that is,
q
0 < x. < 1, i = 1,2,...,q, E x. = 1. (1.2)
i=l
The factor space containing the q components is represented
by a (q 1)dimensional simplex. For q = 2 components, the
factor space is a straight line, whereas for q = 3
components, the factor space is an equilateral triangle, and
for q = 4 components, the factor space is represented by a
regular tetrahedron.
The objectives in the analysis of a mixture response
system are, in general, the same as in any response surface
exploration. That is, one seeks to approximate the surface
with a model equation by fitting an equation to observations
taken at preselected combinations of the mixture com
ponents. Another objective is to determine the roles played
by the individual components. We shall not concern our
selves with this but rather concentrate on the empirical
model fit. Once the model equation is deemed adequate an
attempt is made to determine which of the component combina
tions yield the optimal response. The models used to repre
sent the response in a mixture system are in most cases
different in form from the standard polynomial models. The
first type of model form that we discuss is the canonical
polynomial suggested by Scheffe.
1.2.1 Mixture Models
Scheffe (1958) introduced a canonical form of the poly
nomial model for representing the response in a mixture
system. These canonical polynomial models are derived from
the standard polynomials using the restrictions on the xi
shown in (1.2). With q = 2 mixture components, for example,
the standard second order polynomial model is of the form
2 2
E(Y) = a0 + alx1 + a2x2 + a12 x2 + all x + a22x2 (1.3)
Restrictions (1.2) imply that a0 = a0(xl + x2)'
x2 = l( and x2 = 2( thus (1.3) can be
1 xl(l x2), x2 x2(l xl)
written in the canonical form
E(Y) = 81x1 + 2x2 + 12XlX2
where 1= a0 + al + all' 82 = a0 + a2 + a22' and 812= al2
 a11 a22. There is no constant term in the above canoni
cal form and the pure quadratic terms in equation (1.3) have
been absorbed in the xixj terms.
The general form of the canonical polynomial of degree
d in q mixture components can be written as
q
E(Y) = E 8ixi, for d = 1,
i=l
q q
E(Y) = Z i.x. + Z E .x.x. for d = 2, and
i=l 1. i
q q q
E(Y) = a .x. + EZ .x.x. + E E 6..x.x.(x. x.)
i=l 1 i
q
+ E a ijkx.x.xk for d = 3. (1.4)
1Ki
The fourth degree canonical polynomial in q components is
given in Cornell (1981, p. 64). The general canonical poly
nomial of degree d > 4 in q components does not explicitly
appear in the literature, but is mentioned in Scheffe
(1958). If terms of the form 6ijxixj(xi xj) are removed
from the full cubic model (1.4), then the remaining terms
make up what is referred to as the special cubic model. For
example, for q = 3 components, the special cubic model is
E(Y) = 1x1 + 82x2 + 3x3 + a12 2 + 813X1X3
+ 823x2x3 + 123XlX2X3
Scheffe's canonical polynomial models are used for
approximating the response surface in many mixture systems.
Their popularity stems from the ease in interpreting the
coefficient estimates, especially when the models are fitted
to data collected at the points of the associated designs
(see Section 1.2.2). However, other models have been intro
duced which seem to better represent the response when the
components have strictly additive blending effects. We
present some of them now.
Becker (1968) introduced three forms of homogeneous
models of degree one which he recommends instead of the
polynomial models when one or more of the mixture components
have an additive effect or when one or more components are
inert. A function f(x, y, ..., z) is said to be homogeneous
of degree n when f(tx, ty, ..., tz) = tnf(x, y, ..., z), for
every positive value of t and (x, y, ..., z) (0, 0, ...,
0). These models, which Becker refers to as models HI, H2,
and H3, are of the form
q q
HI: E(Y) = Z B.x. + Z Bijmin(xi, x ) + ...
i=1 1 i
+ 12...qmin(xl, x 2, ..., x )
q q
q q 21
H2: E(Y) = E Bx. + Z Bix. .x ./(x. + x .) +
i= 13 1 3 1
i=1 14i
+ 812...qxlx2..x /(x + x2 + ... + xq)q
q q 1
H3: E(Y) = 6 .x. + Z . (x.x 1/2
i=1 l i
+ Bl2...q(X l2...xq)/q
Each term in the H2 model is defined to be zero when the
denominator of the term is zero.
Draper and St. John (1977) suggest a model which in
cludes inverse terms, 1/xi, in addition to terms in the
Scheffe polynomials. Such a term is used to model an
extreme change in the response as xi approaches zero. The
experimental region of interest is assumed to include the
region near the zero boundary (xi = 0), but does not include
the boundary itself. One example of this type of model is
the Scheffe linear polynomial model with inverse terms
q q 1
E(Y) = Z Bixi + Z B .x .
i=l i=l 
Another model form that is useful in the study of the
response in a mixture system is the model containing ratios
of the component proportions. A term such as xi/xj measures
the relationship of xi to xj rather than the percentage of
each in the blends. Snee (1973) points out that the ratio
model presents a useful alternative to the Scheffe and
Becker models in that the ratio model describes a different
type of curvature. He notes that the curvilinear terms for
the Scheffe and Becker models, when plotted as a function of
xi, are symmetric functions about xi = 1/2, whereas the
ratio term xi/xj is a monotone function when plotted against
xi.
The terms in the ratio models may also contain sums of
the components. For example, with q = 3 components, we
might express the second order model
q1 q1 q1
E(Y) = B + E izi + EE ijziz. + E 68iz.
i=l l
(note the constant term) where zl and z2 are defined as
z1 = xl/(x2 + x3) and z2 = x2/x3. Some terms will be unde
fined if points from the boundary of the experimental sim
plex are included in the design, for example, if x3 = 0,
then z2 = x2/x3 is not defined. Snee (1973) suggests adding
a small positive quantity, c, to each xi in this case.
This, of course, will not be of concern if the experimental
region is entirely inside of the simplex.
When one or more of the components is inactive, Becker
(1978) suggests that a ratio model that is homogeneous of
degree zero in the remaining components is appropriate. In
three components, such a model is of the form
E(Y) = 80 + 1X/(x1 + x2) + 2x2/(x2 + x3)
3
+ 83x /(x + x ) + Z BZ ..h..(x., x )
s3 3 1 3 1i
li
+ 8123h123(xl, x2, x3), (1.5)
where hij and h123 are specified functions that are homoge
neous of degree zero. The function h123 is intended to
represent the joint effect of all three components simulta
neously. If in fitting a model of the form (1.5) we deter
mine the model should be
E(Y) = 80 + a1X1/(X1 + x2) + 812h12(xl, x2)
then component three is said to be inactive and is removed
from further consideration. The model of equation (1.5) may
produce an extreme value near the vertices of the simplex
factor space when there are no inactive components. In this
case it is suggested that a model of the form (1.5) be used
only when the proportions are restricted so that no two of
the xi are simultaneously very close to zero. Becker notes
that other authors who have suggested ratio models have also
used them primarily over a subregion inside the simplex
factor space. Apparently this is where they are most appro
priate.
1.2.2 Experimental Designs for Mixtures
As in the general response surface problem, one of the
major concerns in exploring a mixture system is that of
choosing the experimental design for collecting observed
values of the response that will be used in fitting the
model. Scheffe (1958) proposed the {q,m} simplex lattice
designs for exploring the entire qcomponent simplex factor
space. In these designs, the proportions used for each
component have the m + 1 values spaced equally from zero to
one, xi = 0, 1/m, 2/m, ..., (m l)/m, 1, and all possible
mixtures with these proportions for each component are
used. The number of design points in the {q,m} simplex
lattice design is (m + q 1). The main appeal of these
m
designs is that they provide a uniform coverage of the fac
tor space. Another feature, which Scheffe (1958) demon
strates, is that the parameters of the canonical polynomial
of degree m in q components are expressible as simple linear
combinations of the true response values at the design
points of the {q,m} simplex lattice. The {3,2} simplex
lattice, which consists of six design points, is represented
on the two dimensional simplex in Figure 1 along with the
triangular coordinates (xl, x2, x3).
Scheffe (1963) also developed the simplex centroid
designs consisting of 2q 1 points, where the only mixtures
considered are the ones in which the components present
appear in equal proportions. That is, in a qcomponent
simplex centroid design, the design points correspond to the
q
q permutations of (1, 0, 0, ..., 0), the () permutations of
(1/2, 1/2, 0, ..., 0), the (3) permutations of (1/3, 1/3,
1/3, 0, ..., 0), ..., and the point (1/q, l/q, ..., l/q).
This design alleviates the problem inherent in the {q,m}
simplex lattice designs of observing responses at mixtures
containing at most m components. To give an example, the
q = 3 simplex centroid design is made up of 23 1 = 7
design points, and is equivalent to the {3,2} simplex
lattice design augmented by the center point (xl, x2, x3)
(1/3, 1/3, 1/3). This design is represented in Figure 2.
Scheffe (1963) mentions that the number of parameters
in the polynomial model of the form
q q q
E(Y) = .ix. + EE B. .x.x. + EE 8ikx. xxk
i=l 1 i
+ ... + B12...qxlx2 ... q (1.6)
is 2q 1 and therefore these models have a special rela
tionship with the simplex centroid design in q components.
This relationship is that the number of terms in the model
equals the number of points in the design and as a result
the parameters in model (1.6) are expressible as simple
functions of the responses at the 2q 1 points of the sim
plex centroid design. Polynomial models of the form (1.6)
(2 2
(0,1,0)
Figure 1.
(
(0,1 ,0)
(0,0,1 )
x =I ( x, I )
2 2'2 3
The {3,2}
simplex lattice design.
x =I
(0,0,1)
X3
The q = 3 simplex centroid design.
x I)
~~2 2'2
Figure 2.
therefore are natural models to fit using the simplex cen
troid design.
Ratio models may be desirable when the interest in one
or more of the mixture components is in terms of their rela
tionship to one another, rather than in terms of their per
centages in blends. Kenworthy (1963) proposed factorial
arrangements for ratio variables. An example of the use of
ratios is the following three component system in which the
mixture components are constrained by the upper and lower
bounds:
.2 < x1 < .4, .2 < x2 < .4, .3 < x3 < .5. (1.7)
The ratio variables of interest are zl = x2/xl and
z2 = x2/x3, and it is desired to fit either a first or a
second order polynomial model in zI and z2. For such a
problem, we can define a 22 and a 32 factorial design that
can be used for fitting the first and second order poly
nomial models, respectively, by taking as design points the
intersection of rays passing from two of the three vertices
of the twodimensional simplex through the region of
interest defined by the constraints (1.7). Kenworthy's 22
factorial design is shown in Figure 3.
Becker (1978) uses rays extending from one or more
vertices of the simplex factor space to the opposite bound
aries in developing "radial designs." These designs are
suggested for detecting the presence of an inactive
x :
0 Design Points
x :1 x :I
2 3
Figure 3. Kenworthy's 22 factorial design.
component or in another case a component which has an addi
tive effect, when models containing ratio terms that are
homogeneous of degree zero are fitted.
McLean and Anderson (1966) suggest an algorithm for
locating the vertices of a restricted region of the simplex
factor space which is defined by the placing of upper and
lower bounds on the mixture component proportions. The
vertices of the factor space and convex combinations of the
vertices are the candidates for design points for fitting a
first or second degree polynomial model in the mixture com
ponents. One drawback of the "extreme vertices" design is
that the design points are not uniformly distributed over
the factor space resulting in an imbalance in the variances
of Y(x), see Cornell (1973).
Another method that has been suggested for studying the
response over a subregion of the simplex mixture space is
to transform the q mixture components into q 1 independent
variables. Transforming to an independent variable system
was first suggested by Claringbold (1955) and later proposed
by Draper and Lawrence (1965a, 1965b) and Thompson and Myers
(1968). Standard response surface polynomial models in the
transformed variables can be fitted to data values collected
on standard designs and a design criterion such as the aver
age mean square error of the response can be employed when
distinguishing between designs. Thompson and Myers (1968)
suggest the use of rotatable designs (see also Cornell and
Good, 1970).
Designs other than rotatable designs, such as multiple
lattices and symmetricsimplex designs, to name a few, have
been suggested in the literature for fitting models to a
mixture system which may be appropriate depending on par
ticular experimental situations. However, as the intent
here is not to give an exhaustive list but only a sampling
of available designs, we shall not discuss designs further
but instead state the purpose of this work.
1.3 The Purpose of this Research:
Investigation of Procedures tor Testing a Model
Fitted in A Mixture System for Lack of Fit
A common problem in modeling the response in a mixture
system is that of detecting lack of fit, or inadequacy, of a
fitted model of the form E(Y) = XI when the true model is
of the form E(Y) = X2 + X2 2. The statistical literature
suggests several procedures for testing lack of fit, which
will be described in Chapter Two. In general, the authors
of these procedures are not specific in stating hypotheses
to be tested and do not adequately discuss the power of
their procedures.
The major purpose of this research is to investigate
the power of two of the testing procedures appearing in the
literature in detecting the inadequacy of a fitted model
when the general form of the true model is specified. Our
findings for a "check points" lack of fit testing procedure
are presented in Chapter Three while Chapter Four contains
findings for a "near neighbor" lack of fit testing proce
dure. For both procedures, explicit formulas for the power
of the test are given in terms of cumulative probabilities
of either the noncentral F or doubly noncentral F distribu
tion, which are derived by assuming that the response obser
vations are independent and normally distributed. Addition
ally, we propose methods for maximizing the power of the
testing procedures. In the final chapter, we make some
concluding comments concerning both of these procedures.
CHAPTER TWO
LITERATURE REVIEWTESTING FOR LACK OF FIT
2.1 Introduction
Let us return to the general response surface problem
and assume the true response is to be approximated by
fitting a model of the form
E(Y) = XI1 (2.1)
where Y is an Nxl vector of observable response values, X is
an Nxp matrix of known constants, and El is a pxl vector of
unknown regression coefficients. We wish to consider the
situation in which the true model contains terms in addition
to those in the fitted model. Then the true model has the
form
E(Y) = X1 + X2a2 (2.2)
where X2 is an Nxp2 matrix of known constants, and 12 is a
P2xl vector of unknown regression coefficients. We assume
that the vector Y has the normal distribution with
var(Y) = a 2N
It is desirable to determine the suitability of the
fitted model given by Eq. (2.1) when in reality the true
model is of the form given by Eq. (2.2). The process of
19
making this determination is referred to as testing for lack
of fit of the fitted model.
There are three general approaches to testing for lack
of fit. The first approach requires that there be replicate
observations of the response at one or more design points,
and involves partitioning the residual sum of squares from
the fitted model into a sum of squares due to lack of fit
and a sum of squares due to pure error. A large value for
the ratio of the mean square due to lack of fit to the mean
square due to pure error provides evidence for lack of fit.
If replicate observations are not available then the
above approach to testing for lack of fit cannot be used.
Green (1971), Daniel and Wood (1971), and Shillington (1979)
have proposed alternative methods that are applicable in
this case. Their approach is to group values of the
response which are observed at similar settings of the
independent variables and to call these grouped values
"pseudoreplicates" or "near neighbor observations." They
then treat these pseudoreplicates as they would treat true
replicates to form statistics for lack of fit testing,
although arriving at their respective statistics through
different approaches.
A third approach to testing for lack of fit involves
the use of "check points." In this method a model of the
form (2.1) is fitted to data at the design points and
additional observations are collected at other points in the
experimental region. The additional points other than the
design points are called check points, and the data at these
check points are not used in fitting the model. Lack of fit
is tested by comparing the values of the response observed
at the check points to the values of the response which the
fitted model predicts at these same check points.
We now discuss the first method mentioned above of
testing for lack of fit which involves partitioning the
residual sum of squares.
2.2 Partitioning the Residual Sum of Squares
The method for testing lack of fit which makes use of a
partitioning of the residual sum of squares from the fitted
model requires there be replicate observations of the
response at some of the design points (Draper and Smith,
1981, p. 120). When a model of the form (2.1) is fitted,
the residual sum of squares is defined as
n.
n 1 2
SSE = EZ (Y Y )
i=1 j=l 1
1
=Y'(I' X(X'X) X')Y
where n is the number of distinct design points, ni > 1 is
the number of replicate observations at the ith design
point, Yij is the jth observed value of the response at the
ith design point, Yi is the value which the model of the
form in Eq. (2.1), fitted by ordinary least squares
techniques, predicts for the response at the ith design
n
point, and N = Z n. Using the replicated observations
i=l 1
only, a pure error sum of squares can be calculated as
n "i 2
SSE = E E (Y. Y. )
pure i=l 13 '
i=1 j=1
where Yi. is the average of the values of the response
observed at the ith design point. The sum of squares due to
lack of fit can be obtained by taking the difference
LOF
=SSE SSE
pure
This partitioning of the residual sum of squares is
displayed in the analysis of variance table in Table 1.
Table 1. Analysis of Variance
Partitioning the Residual Sum of Squares.
Source
of Variation
Regression
Residual
Pure Error
Lack of Fit
Total(corrected)
Sum
of Squares
b{X'Y (1'Y)2/N
SSE
SSEpure
SSLOF
Y'Y (l'Y)2/N
Degrees
of Freedom
p 1
N p
N n
n p
N 1
bl represents the ordinary least squares estimator of B in
1
model (2.1), b = (X'X) X'Y, and 1 is an Nxl vector of
ones.
ones.
Mean
Square
MSE
MSEpure
MSLOF
To test the hypothesis of zero lack of fit, that is
HO: lack of fit = 0 or E(X) = XEI, an F statistic is formed
MS
LOF
F = F (2.3)
MSE
pure
which possesses a central F distribution if the true model
is of the form (2.1), but has a noncentral F distribution if
the true model is of the form (2.2). In other words
F ~ F
np,Nn
under H : E(Y) = X8_ and
F ~ F'
np,Nn;X2
under H : E(Y) = XB + X2_ where X2 is the noncentrality
a 1+ 22 2
parameter 2 = B(X2XA)'(X2XA)B2/22 and A = (X'X) X'X2
Under H E(MSLOF) = 2 + (X2 XA)'(X2 XA) 2/(np) and
E(MSEpure) = o2 (Draper and Smith, 1981, p. 120), hence HO
is rejected in favor of Ha if the value of F in (2.3)
exceeds the upper 100a percentage point of the central F
distribution, Fa;np,Nn. When HO is rejected, we conclude
that a significant lack of fit is present.
Draper and Herzberg (1971) demonstrated that the lack
of fit sum of squares can be partitioned into two
statistically independent sums of squares, SSL1 and SSL2'
when there are replicate observations at the center of the
response surface design and when the basic design without
center points is nonsingular. If the true model and the
fitted model are of the same form as in equation (2.1) then
the two F ratios FL1 = [SSLl/(n p 1)]/MSEpure and
FL2 = SSL2/MSEpure are both distributed as central F random
variates, with respective numerator and denominator degrees
of freedom (n p 1), (N n) for FL1 and 1, (N n) for
FL2. If the true model is of the form shown in equation
(2.2), then FL1 and FL2 are both distributed as noncentral F
random variates. The expected values of SSL1 and SSL2 are
used to show what functions of E2 are testable with FL1 and
FL2*
Two examples are presented by Draper and Herzberg to
illustrate this testing for lack of fit. The first example
makes use of a first order orthogonal design in k factors
augmented with center point replicates for fitting a first
order polynomial model. If the true model is of the second
order, then FL2 can be used to test a hypothesis concerning
the parameters associated with the pure quadratic terms in
the model. If all such parameters are zero, then FL1
provides a check on interaction terms. The second example
illustrates the fitting of a second order polynomial model
to a second order design with all odd design moments of
order six or less zero. If the true model is third degree,
then FL1 can be used to test the significance of the third
order terms, while FL2 tests terms of order greater than
three. The partitioning of SSLOF into SSL1 and SSL2 and the
corresponding tests of hypotheses are also given in Myers
(1971, p. 114119), for the special case of fitting a first
order polynomial model to a 2q factorial or a fraction of a
2q factorial design augmented with center point replicates
and the true model is of the second degree.
A more complete partitioning of the lack of fit sum of
squares in an attempt to obtain a more detailed diagnosis of
the lack of fit of the fitted model is given in a technical
report written by Khuri and Cornell (1981). The lack of fit
sum of squares, which has n p degrees of freedom, is
partitioned into n p independent sums of squares, each
having one degree of freedom. The expected values of these
single degreeoffreedom sums of squares are used to
identify at most n p linearly independent causes for the
lack of fit variation. Tests of significance are performed
on the assumed contributing causes. This method enables the
screening of all subsets of E2 in order to identify those
subsets which are most responsible for lack of fit of the
fitted model.
We shall now discuss the second general approach used
in lack of fit testing, which is to test for lack of fit by
making use of response values observed at points which are
near neighbors in the factor space when true replicate
observations are not available.
2.3 Testing for Lack of Fit Without
Replicated ObservationsNear Neighbor Procedures
Green (1971) suggests the following approach when
testing for lack of fit if there are no design points at
which replicate observations of the response are
available. The N observed values of a response, Y,
considered a function of only one variable, x, are divided
into g groups, by grouping observations which have similar
values of x. Green hypothesizes a model of the form Y= Ha +
e, where Y is an Nxl vector of observable responses, H is an
Nxm matrix whose columns correspond to known functions of
the variable, x, is an mxl vector of unknown regression
coefficients, and e is the Nxl vector of random errors,
2
N ~ N (0, o N).
Green's method assumes that the vector of differences
(EY Hg) can be well approximated by a dth order polynomial
in x within each of the g groups, d > 1. An alternative
model of the form
Y = H v + n +
is given, where is distributed as NN(Q, o21N), H1 is an
Nx [g(d + 1) + ml] matrix of known constants, u is a
[g(d + 1) + ml]xl vector of regression coefficients, and .,
as Green states is "a small vector." The first g(d + 1)
columns of H1 correspond to the polynomial terms for the g
groups (with (d + 1) terms for each group), the rightmost
mI < m columns in H1 correspond to terms that are in the
fitted model, but are not represented among the g(d + 1)
polynomial terms in the alternative model.
Under the assumption that a = Q, the presence of lack
of fit is tested by using the test statistic:
Y'[H (HHl) H H(H'H) 1H']Y/[g(d + 1) + mi m]
F = *
1
Y'[I H (HH1I) H{]Y/[N g(d + 1) mi]
(2.4)
This statistic is of the same form as the F statistic used
in the standard multiple regression test of a postulated
model against a more general one which includes the
postulated model as a special case. Lack of fit is
suspected if the calculated F ratio in (2.4) is greater
than Fa;g(d+l)+mlm, Ng(d+l)ml where this latter quantity
is the upper 100a percentage point of the central F
distribution.
Green notes that when there is no lack of fit, the
quadratic forms Y'[H1(HIH ) H H(H'H) H']Y and
Y'[I H (HIH,) Hi]Y are distributed independently as
o2X2 with g(d + 1) + mi m and N g(d + 1) mi degrees of
freedom, respectively. In this case the F ratio in (2.4)
possesses a central F distribution. If there is lack of fit
on the other hand, then these two quadratic forms are
distributed as noncentral chisquares, multiplied by 02,
with respective noncentrality parameters
I = [H + '[H (HiH1) H' H(H'H) H'][H + n]
and ; = n'[I H1(HH ) HI]n Thus the assumption that
n = 0 can affect the power of the test, since if n # 0 the
expected value of MSE is greater than a2, where MSE is the
quadratic form in the denominator of the F ratio. Hence if
n 0 the probability of calculating a large F value is
reduced, and we are less likely to detect lack of fit using
an upper tailed rejection region.
Daniel and Wood (1971) suggest another method for lack
of fit testing when replicated observations of the response
are not available. They make use of "near replicates" to
obtain an estimate of a, which is the standard deviation of
the observable responses in the true model. The value of
the estimate a is compared to the square root of the
residual mean square from the analysis of the fitted
model. Lack of fit is indicated if the square root of the
residual mean square is large compared to the estimate o.
To determine when observations are near replicates so that
an estimate of a can be found, they define the squared
distance between any two data points, j and j', to be
measured by
K
D2 = x. ij.)/5y]2
3J' i=1
where xij and xij, are the values of the ith independent
variable corresponding to the observations yj and yj,,
respectively, i = 1, 2, ..., K, and bi is the ordinary least
squares estimate of the ith regression coefficient. In the
denominator, s is the square root of the residual mean
square for the fitted model.
To obtain an estimate of a from near replicates, let
And = dj d ,j, n = 1, 2, ..., (2), where dj and d are
the residuals at points j and j', respectively, and where
there are N data observations in the experiment. Since the
expected value of the range for pairs of independent
observations from a normal distribution is 1.1280, a running
average of the And's is calculated and their average is
multiplied by .886 = (1/1.128) to get a running estimate,
sn, of a. That is, s = .886 n And/n The closest pair
of observations as judged by Dj, is used to begin the
running estimate, the next closest pair (next "nearest
neighbors") is used for A2d, and the procedure continues
until sn "stabilizes." The stabilized value of sn is used
to estimate a.
A third method for testing for lack of fit without
replication is given by Shillington (1979). The fitted
model is of the form
Y = XB + E (2.5)
where Y (Nxl), X (Nxp), and (pxl) are defined as in (1.2)
and E ~ N (0, a IN) The test for lack of fit of the
fitted model is a test for whether the true model has the
form
y = X1 + 6 + E 1
where 6 (Nxl) is a fixed effect quantifying the departure of
(2.5) from the true model.
Shillington assumes that the data can be grouped into g
cells, with nj observations in the jth cell, determined in
advance. Letting Cj refer to the jth cell, j = 1, 2, ...,
g, a vector of cell averages is written YC (gxl), where the
jth element of YC is the average of the observed responses
in Cj. The matrix XC of independent variables associated
with Y is the gxp matrix where the elements in the jth row
n.
are x'. = i x!./n. that is, row j of X is the row
i. =1 
vector x'. The matrix XC is assumed to be of full rank
p < g. Also within each cell are defined the differences
W.. = Y.. Y. i C. j = 1, 2, ..., g, where Y is
the jth element of Yc.
The two independent data sets, YC and {Wij} with g and
N g degrees of freedom, respectively, are used to find two
independent estimates of o2. The first estimate is written
as
9 2
MSEB = E n (Y x' /(g )
j=l B
where g4 is the weighted least squares estimate of g using
the regression of cell means, YC, on XC. The second
estimate of o2 uses the within cell deviations on cell
means, {Wij}, and is
n.
g n 2
MSEJ = Z (W. W. ) /(N g r),
j=1 i=l
where r is the rank of an Nxp matrix with rows equal to
x! x' i C C., j = 1, 2, .., g. If the matrix of
13j .j 3
independent variables, corrected for cell means, is of full
rank, then r = p. Here 7ij is the estimate of Wij from the
regression of cell residuals {Wij} on the associated vectors
of independent variates, x'. x'.
If the fitted model is the correct model, then MSEB and
MSEW are independent estimates of a2 and the ratio MSEB/MSEW
is an F statistic with g p and N g r degrees of
freedom. When all observations in a cell have the same
settings of the independent variables, that is, the
observations are truly replicates for all cells, then this F
statistic is identical to the F statistic in the usual lack
of fit test in which the residual sum of squares is
partitioned into lack of fit and pure error sums of squares,
as given in Draper and Smith (1981, p. 120).
If the true model is Y = X8 + 6 + E however, and if
we let X'6 = 0 and 2 = 6/N, then
E(MSEB) = o2 + [I XC(X'XC)1'] B/(gP)
n.
where B (gxl) has jth component equal to E 6../n .
B i=
Furthermore, with this latter true model form
E(MSE ) = 02 + 67(I XW(XJ XW) X)6W/(N g r)
where _W has the components ij .j, i e Cj, j = 1, 2,
..., g. The matrix XW (Nxp) has the rows x x'
i e C., j = 1, 2, ..., g. The power of the F test,
F = MSEB/MSEW, depends on the relative bias of the estimates
of 2, that is, the biases in MSEB and MSEW.
Shillington states that the power of the F test which
makes use of F = MSEB/MSEW is maximized by forming cells so
that the bias of E(MSEW) is minimized. This is the same as
forming cells so that the within cell variation in 6 is
minimized. Shillington (1979, p. 141) also states,
"Observations with near covariate (independent variable)
values might be expected to have similar 6 values, since we
assume that 6 varies in some continuous but unknown fashion
with X. This justifies the usual procedure of forming
groups by collapsing observations with adjacent covariate
values. Indeed, if covariates do not vary within cells we
have the usual lack of fit test and maximum power."
By imposing a further structure on the form of 6, it is
shown that if the F test has an upper tailed rejection
region, the power is maximized by selecting the group sizes
as nj = 2, j = 1, 2, ..., g. Finally, Shillington suggests
that in the presence of more than one independent variable
problems in grouping may arise, and in this case it may be
wise to perform a different lack of fit test for each
parameter. Following this approach, an example is given
which suggests testing lack of fit for each of the p
independent variables separately may be more powerful than
trying to form groups based on all independent variables at
once.
In summary, all the approaches we have discussed for
testing for lack of fit when replicate observations of the
response are not available at any of the settings of the
independent variables make use of grouping the observed
response values according to similar values of the
independent variables. The observations falling in such
groups are referred to as "pseudoreplicates" or "near
neighbor observations." These pseudoreplicates are used to
estimate the true variance of the observations, a2, but a
completely unbiased estimate of 02 cannot be attained unless
true replicate observations are available. In each case,
the power of the lack of fit testing procedure is reduced
because an unbiased estimate of a2 is not attainable. We
now turn to the use of check points for lack of fit testing.
2.4 Testing for Lack of Fit with Check Points
An alternative to the two approaches to lack of fit
testing already discussed is the method which makes use of
check points. We assume a model of the form E(Y) = Xl as
given in (2.1), is fitted in a response surface system, but
that the true model is of the form E(Y) = X1I + X22 as
given in (2.2). The parameters, a8, in the fitted model are
estimated by ordinary least squares techniques, making use
of the values of the response observed at the design
points. After the model is fitted, values of the response
are observed at additional points in the experimental region
called "check points." The observed response values at the
check points are compared to the values which the fitted
model predicts at these same check points. It is important
to note that the observed values of the response at the
check points are not used in fitting the model initially.
Snee (1977) gives four methods of validating regression
models, one of which is the collection of new data to check
predictions from a previously fitted model. In a designed
experiment these new data take the form of check points.
Snee suggests that the inclusion of a small number of check
points in any designed experiment is a "worthwhile"
procedure.
Scheffe (1958) proposed a test for lack of fit when the
{3,2} simplex lattice design is used for fitting a second
order canonical polynomial model in three mixture
components. It is desired to use the observed value of the
response at (1/3, 1/3, 1/3) as a check point blend. The
test statistic proposed is the t statistic of the form
Y Y
t = (2.6)
[var(Y )]
where Y is the observed value of the response at the check
point, and Y is the value of the response predicted at the
same point by the second order model which is fitted by
ordinary least squares techniques to the observed response
values at the six design points of the {3,2} simplex
lattice. The response value observed at the point
(1/3, 1/3, 1/3) is not used in fitting the model. Lack of
fit is inferred if the absolute value of the calculated t
value in equation (2.6) is larger than the corresponding
tabled t value.
In the denominator of the t test of equation (2.6), the
variance of the difference Y Y is shown to be
var(Y Y) = var(Y) + var(Y)
= (44/27r)o
when r replicates are taken at each design point. The
estimate of the variance of Y Y is (44/27r)o2, where o2 is
calculated from the replicated response values at the design
points.
Scheffe (1958) also alludes to a test for lack of fit
when several check points are used simultaneously. When
there are k check points, the test for lack of fit is an F
statistic of the form
1
F = 2 (2.7)
ka
where d' = (Y Y 2 2' ... Yk Yk) and V = o2V
var(d). Formulas are given for the elements of VO in the
special case when the check points are the design points of
the {3,2} simplex lattice. Lack of fit is suspected if the
calculated value of the F statistic given in (2.7) is larger
than the corresponding tabled F value.
Gorman and Hinman (1962) suggest the same t test in
equation (2.6) that Scheffe (1958) suggested for a check
point taken at (1/3, 1/3, 1/3) to test for lack of fit in a
second order polynomial model fitted from a {3,2} simplex
lattice design. They suggest using (1/3, 1/3, 1/3) as the
location of the check point because the observation at this
point may later be used to fit the next more complex model,
the special cubic, if the second order model is found to be
inadequate. They state that in general for the second order
polynomial model as well as higher order models, check
points should be taken in regions of particular interest, of
which there are usually many in any blending study.
Further, they suggest that the number of check points
depends on individual experimental situationstechnical
background, precision required, cost of materials and
analyses, and probability of requiring a more complex
model. However, no specific criterion is given by Gorman
and Hinman for selecting the location of the check points.
Gorman and Hinman (1962) indicate that a t test at a
check point other than at (1/3, 1/3, 1/3) takes the same
form as the statistic of equation (2.6),
Y Y
t 1/2
[var(Y) + var(Y)]/2
with the additional condition that if several check points
are taken, say for example k points, the method of checking
the fit is to compute the t value at each location and refer
these calculated t values to the 100(a/2k) percentage point
of the central t distribution rather than the 100(a/2)
percentage point.
Kurotori (1966) gives an example of a mixture
experiment where the response is the modulus of elasticity
of a rocket fuel, which is a mixture of three components,
binder (xl), oxidizer (x2), and fuel (x3). The factor space
of feasible mixtures is a subspace inside the two
dimensional simplex or triangle where all three components
are present simultaneously. "Pseudocomponents" are defined
and in the pseudocomponent system a special cubic model is
fitted to data collected at the points of the q = 3 simplex
centroid design (Figure 4). A check for adequacy of fit is
made by using three check points and the response values at
the check points are used only for testing the fit of the
model and not for fitting the model initially.
The reason for the choice of the particular check point
locations by Kurotori is that, as he states, "They are the
most remote mixtures from the seven design points." The
lack of fit test is an F statistic of the form
2
F = (2.8)
2 3
where s (Y Y ,) for the i = 1, 2, 3 check points
2 i=l
and 02 is an estimate of measurement error from a previous
analysis. Kurotori admits that the use of the F statistic
x :1
I
(1,0,0)
9 Design Points
O 0 Check Points
SE 3 ') o
S 2 I 112
e(o,,o) a the ser)
2 '2 '2)
Figure 4. Kurotori's rocket fuel example,
xI', x2', and x3' represent pseudocomponents.
in Eq. (2.8) for lack of fit testing may be risky because
the predicted values at the check points are correlated
(correlation of .5), although the observed values are not
correlated. Kurotori suggests individual t tests as
proposed by Scheffe (1958) might be the preferred procedure.
Snee (1971) repeats Kurotori's rocket fuel example
using the same F test for lack of fit as Kurotori and makes
the comment that the Yi's at the check points are
correlated. In stating that the F test is not an exact
test, he nevertheless offers no solution in the form of an
exact test.
39
In summary, only Scheffe refers to an exact F test when
several check points are considered simultaneously for
testing for possible lack of fit of a model fitted in a
mixture space, and his development is limited to the special
case where the check points are the design points used to
fit the model initially. No criterion is proposed by
Scheffe for selecting other locations for the check points.
CHAPTER THREE
AN OPTIMAL CHECK POINT METHOD FOR TESTING
LACK OF FIT IN A MIXTURE MODEL
3.1 Introduction
In Chapter Three we investigate the problem of testing
for lack of fit of a linear model fitted in a mixture
space. The testing is to be accomplished with the use of
check points. We assume that an experimental design is
specified, and that the fitted model is of the form
E(Y) = X1 (3.1)
where Y is an Nxl vector of observable response values, X is
an Nxp matrix of known constants and rank p, and 81 is a
vector of p unknown regression coefficients. The true model
is assumed to be of the form
E(Y) = X1 + X 2 (3.2)
1 22
where X2 is an Nxp2 matrix of known constants and g2 is a
vector of p2 unknown regression coefficients. Throughout
our development, we will assume that the random vector Y has
the normal distribution with variancecovariance matrix
equal to a IN.
In our investigation we wish to determine the proper
testing procedure to follow in deciding whether the fitted
model exhibits lack of fit. In order to optimize the lack
of fit testing procedure, we will determine the location of
the check points so that the power of the test is maximized.
3.2 Testing for Lack of Fit in the Presence of
an External Estimate of Experimental Error Variation
3.2.1 The Test Statistic
We wish to test the performance or fit of a fitted
model in a mixture space when the true model possibly
contains terms in addition to those in the fitted model.
The fit of the model is to be tested by a test which makes
use of the response values observed at certain locations
called "check points" in the experimental region, by
comparing them to the values which the fitted model predicts
at the same check points. The observed values at the check
points are not used for estimating the coefficients in the
fitted model and are assumed to represent the values of the
true surface at the check points.
Let us define the vector of differences
d = (Y* Y*)
(Y* Y Y Y* Y* y*)i
1 1' 2 21 k k
where Y* = 1, 2, ..., k are observed response values at
k check points and YV, i = 1, 2, ..., k are response values
predicted at the k check points by the fitted model,
VY = xt'b where b is the ordinary least squares estimator
of al and where x*' is the ith row of X*, the kxp matrix
1 1
whose columns are of the same form as the columns of X but
with its rows evaluated at the k check points. Note that
if 3 = 0, then E(d) = 0 and if a 2 0, then
2 2 2
E(d) = (X2 X*(X'X)X'X2). Let V represent the
variancecovariance matrix of the random vector d.
Then V = o V where
1
V = Ik + X*(X'X) X*'
and where Ik is the identity matrix of order kxk.
We assume that an unbiased estimate of a2 is available
^2
and we denote this estimate by a where the subscript ext
Sext'
^2
stands for external, and a is independent of the model
being fitted. The test statistic for the hypothesis of zero
lack of fit H : E(d) = 0 is
1
d'V d/k
0
F = (3.3)
^2
ext
(see Scheffe, 1958, p.358). It will be shown later in this
section that the F ratio in Eq. (3.3) possesses either a
central F distribution or a noncentral F distribution,
depending upon whether the true model is represented by Eq.
(3.1) or Eq. (3.2).
^2
The variance estimate a that appears in equation
ext
(3.3) is ordinarily generated from replicated observations
at some of the design points in the experiment. We assume
^2
that aext is a constant multiple of a central chisquare
random variable with v degrees of freedom. This is written
as
^2
a e = SSE /v
ext pure
2 2
= (a 2/v)(SSEpe/
pure
2 2
where SSEpure/ ~ X. Note that SSEpure denotes the
portion of the residual sum of squares due to replication
variation from the fitted model. The residual sum of
squares from the fitted model may be partitioned into
SSEpure and SSLOF only if replicated observations are
collected at one or more design points. For the case where
replicate observations are collected at all of the design
points
n n.
SSE = Z (Yij Yi. )
i=l j=1
where n is the number of distinct design points, n. > 2 is
1
the number of replicates at the ith design point, Yij is the
jth observation at the ith design point, and Y. is the
1.
average of the ni observations at the ith design point.
n
Here SSEpure has v = Z (ni 1) degrees of freedom.
i=l1
When the fitted model and the true model are of the
same form as defined by Eq. (3.1), the quantity d'Vld/a2
(3.lr th quntit d'
possesses a central chisquare distribution (Searle, 1971,
p.57, Theorem 2). However, when the true model is of the
1 2
form specified by Eq. (3.2), d'V d/a possesses a
noncentral chisquare distribution. Thus when the true
model is of the form in Eq. (3.1),
1 2 2
d'V d/ao 2 X
0 Xk'
but when the true model is of the form in Eq. (3.2),
1 2 2
d'V d/o2 ~ x ,
where in the second case the noncentrality parameter X has
the form
S= E(d)'V E(d)/2a2
1 0
1 2
= X* X*A)'Vo (X X*A) 2/2
2 2 0 2 2
1
The matrix A = (X'X) X'X is called the alias matrix and is
2
of order pxp2. In X1, the matrix X* is of order kxp2 and
has the same relationship to X2 as X* has to X.
2
Since SSE /a is statistically independent of
pure
1 2
d'V d/a then under model (3.1) the test statistic
0 0
1 2
d'V d/ko
F = 2
SSE pure/v
pure
1
d'V d/k
= ^2
ext
will have a central F distribution. When the true model
contains terms in addition to those in the fitted model then
F will have a noncentral F distribution. We write these two
cases as
F ~ Fk,
k,v
under model (3.1), and
F ~ F'
k,v;Xi
under model (3.2), where the noncentrality parameter is
x = (X X*A)'V0 (X X*A)_2/2 2.
3.2.2 The Testing Procedure and an Expression for the Power
of the Test
Given that the form of the fitted model is defined as
Eq. (3.1), the expected value of the numerator of the F
statistic in Eq. (3.3) will depend on the form of the true
model. For the case where the true model is expressed as
Eq. (3.2),
1
E(numerator) = E(d'V d/k)
0
2 2
= ( /k)EXx2
= (a /k)(k + 211)
2 2
= 2 + 20 X1/k
= o2 + BA1 /k, (3.4)
1
where A1 = (X* X*A)'V0 (X X*A). However, when the true
model is Eq. (3.1), 8 = 0 and in this case A = 0 so that
2 1
2 ^2
E(numerator) = 02. Also a is an unbiased estimator of
ext
o2 and
^2 2
E(a ) = a (3.5)
ext
Therefore the ratio E(numerator)/E(denominator) where
^2
the denominator is ext will equal unity under model (3.1),
that is, when there is no lack of fit. Under model (3.2),
the ratio will be greater than or equal to unity so lack of
fit should be suspected if the calculated F ratio in
equation (3.3) is large. We can thus use an upper tailed
rejection region to reject the hypothesis of zero lack of
fit. The power of the test is
PI k ,v ;X > a; k, v 1
P{Fv > F,; };k,I
where F is the upper 100a percentage point of the
a ;k,v
central F distribution with k numerator degrees of freedom
and v denominator degrees of freedom.
It is worth noting that from Eq. (3.4) and Eq. (3.5)
testing the hypothesis that 2 = 0 is equivalent to testing
the hypothesis that X1 = 0, assuming A1 is positive
definite. Thus testing a null hypothesis of zero lack of
fit using the proposed testing procedure involving the F
ratio in (3.3) may be expressed as a test of the hypotheses:
H: 1 = 0
H: a1 > 0.
a 1
3.2.3 A Method for Locating Optimal Check Points
Once a design for fitting model (3.1) in a mixture
space is chosen and the number of simultaneous check points
is decided on, say k > 1, the next step is to determine
where in the mixture space we should place the k check
points so as to maximize the power of the test for lack of
fit. The location of the check points is to be made
independently of the value of .
2
The power of the upper tailed F test for lack of fit is
an increasing function of X1 (see Appendix 1 for proof, with
X2 = 0). Therefore, to maximize the power of the test we
maximize the value of X1 defined as
1
Xl = S2AI82/2a2
where A = (X X*A)'V0 (X* X*A), by properly selecting
the k check points whose coordinates are defined in X*. To
maximize the value of X1, we shall concentrate on the matrix
A1.
The matrix A1 is a square matrix of order p2xp2 and is
a scalar quantity when p2 = 1. By maximizing the scalar
quantity A1 with respect to the k check points, the power is
maximized no matter what the value of .2 Maximizing the
2
scalar A1 can be accomplished by using The Controlled Random
Search Procedure given by Price (1977). This procedure is
described in Appendix 2. As a computational aid, A1 can be
expressed as
V + (X* X*A)(X* X*A)'
A = V2 1 (3.6)
1 V
when p2 = 1, where the symbol IBI denotes the determinant of
the square matrix B. Thus the computations reduce to
evaluating two determinants rather than inverting VO (see
Scheffe, 1959, Appendix V, p.417).
When p2 > 1 and A1 is no longer a scalar, maximizing X1
(and thus maximizing the power of the test) cannot be
accomplished without specifying 02 In this case we make
use of a lower bound for 1l (Graybill, 1969, p.330, Theorem
12.2.14(9)) defined as
lmin/22 < X1
(where min is the smallest eigenvalue of A1) to be used in
place of \l. Hence an approximate solution to the
maximization of X1 will be achieved by finding the k
simultaneous check points (using Price's procedure) that
maximize min' the smallest eigenvalue of A1. In other
words when p2 > 1, and in order to avoid specifying 2' we
seek to maximize a lower bound value for X1. This
maximization does not depend on the value of 02.
There are cases where the matrix A1 is of less than
full rank (less than rank p2) or equivalently where the
matrix A1 is positive semidefinite so that umin will be
equal to zero no matter which check points are selected.
One such case occurs when k < p2 (when the number of check
points is less than the number of parameters in the true
model which are not in the fitted model) since when k < p2
rank(Al) = rank[V I X* X*A)]
rank(X* X*A),
2
and so rank(Al) < min(k, p2) because the matrix (X* X*A)
is of order kxP2. Therefore when k < p2, the rank of A1 is
at most k so that A1 is of less than full rank. Since umin
must be equal to zero when Al is positive semidefinite, an
alternative method to that of maximizing umin to select
optimal check points must be found when A1 is positive semi
definite in order to produce a positive lower bound for X1.
In this pursuit, let us write X1 as
1 = 2AA2/202
= PAP'_2/202
= 8[P1:P2] diag[Al, A2=0][P1:P2] '2/2o2
= PlA 1P 12/202
where A is a diagonal matrix with elements equal to the
eigenvalues of Al, P is an orthogonal matrix whose columns
are orthonormal eigenvectors of AI, A1 and P1 correspond to
the positive eigenvalues of AI, while A2 = 0 and P2
correspond to zero eigenvalues of AI. Then by Theorem
12.2.14(9) in Graybill (1969) we can write
Sz'z/2o2 < 1 (3.7)
min p
where in is the smallest positive eigenvalue of A,, and
min
z = P'S Thus by Eq. (3.7), an approach to maximizing a
 12
positive lower bound for X1 when A1 is positive semi
definite is to select check points that maximize the
smallest positive eigenvalue of A1. It must be noted,
however, that this method can only be used when
a2 e n C(Pl), where C(P1) denotes the column space of P1
and n C(PI) denotes the intersection of all such spaces
which can be obtained at all possible check points
locations. This is because, in general, z'z in (3.7)
depends on the location of the check points through its
dependency on Pi. If, however, 2 e nC(P ), then
zz = P1 = P = 22 since 'P2 = 0.
It follows that when 2 e n C(P ), mn z'z/2o
2 1 min 
+ 2 +
= min 2+/2o and only mn depends on the location of the
check points.
3.3 Testing for Lack of Fit When MSE Is Used
to Estimate Experimental Error Variation
3.3.1 The Test Statistic
In this section we shall show that when an external
estimate of a2 is not available and the residual mean square
(MSE) from the fitted model of the form (3.1) must be used
as an estimate of a2, the test statistic
l
d'V d/k
F ME (3.8)
MSEpossesses a central F distribution when the true model is
possesses a central F distribution when the true model is
Eq. (3.1), but possesses a doubly noncentral F distribution
when the true model is Eq. (3.2).
In the initial section of this chapter, the quantity
d'V d/2 was said to possess a central chisquare
distribution or to possess a noncentral chisquare
distribution, depending on whether the true model was
specified by Eq. (3.1) or Eq. (3.2). Now, the residual sum
of squares from the fitted model is defined as
N
SSE = E (Y Y)
i=l
1
= Y'(I X(X'X) X')Y
and it is easy to show (Searle, 1971, p.57, Theorem 2) that
SSE/a2 possesses a central chisquare distribution if the
true model is Eq. (3.1), but under model (3.2), SSE/a2
possesses a noncentral chisquare distribution. This is
expressed as
2 2
SSE/a xp
Np
under model (3.1), and
2 2
SSE/a X.
Np,12
under model (3.2), where the noncentrality parameter \2 is
2
X2 = (X2 XA)'(X2 XA)2/2o0
The distributional form of the test statistic in Eq.
(3.8) is derived by knowing that the quantities
d'Vo d/o2 and SSE/o are statistically independent (see
Appendix 3), so that
d'vold/ko2
12
MSE/o
d'V d/k
MSE
is distributed as a central F when the true model is Eq.
(3.1), but when the true model is Eq. (3.2) the F ratio is a
doubly noncentral F, that is, under model (3.1),
F ~ F
k,Np
and under model (3.2),
F~ F"
k,Np;X ,X2
3.3.2 The Rejection Region and its Relation to the Power of
the Test
In Appendix 1 it is shown that if k, Np, and X2 are
fixed, then the power of the F test using the ratio (3.8) is
a function of the location of the rejection region (upper
tailed or lower tailed) of the test. The power increases
with increasing values of the numerator noncentrality
parameter, X1, when the test is an upper tailed test. The
power decreases with increasing values of X1 when the test
is a lower tailed test. This means that to study ways of
increasing the power of the test, we have to determine
whether the test is an upper tailed test or a lower tailed
test. Similarly, for fixed values of k, Np, and X1, the
power of the F test is a decreasing function of X2 for an
upper tailed test, and is an increasing function of X2 when
the F test is a lower tailed test (Scheffe, 1959, p. 136
137).
To decide if the test is an upper tailed test or a
lower tailed test, we recall from Section 3.2.2 that if the
true model is Eq. (3.1) then the expected value of the
numerator of the F statistic in (3.8) can be written as
E(numerator) = 02,
and if the true model is Eq. (3.2),
2 2
E(numerator) = a + 2a X1/k (3.9)
= 02 + QA 2/k
1
where the P2xp2 matrix A is A = (X* X*A)'V (X* X*A).
Similarly, it can be shown that if the true model is Eq.
(3.1), the expected value of the denominator of the F
statistic in (3.8), where the denominator equals MSE, is
E(denominator) = E(MSE)
= 2,
but if the true model is Eq.(3.2),
E(denominator) = E(MSE)
[o2/(N p)]Eyp,2
= [a2/(N p)][N p + 2X2]
= a2 + 22 2X/(N p) (3.10)
= 2 + a2A2 2/(N P)
where the P2xP2 matrix A2 is A2 = (X2 XA)'(X2 XA). Thus
the ratio E(numerator)/E(denominator) will equal unity if
the true model is Eq. (3.1), but if the true model is Eq.
(3.2), the ratio is greater than unity if B_~A12/k >
8_A 2 2/(N p). In this latter case we reject the null
hypothesis of zero lack of fit if the calculated value of
the F ratio in (3.8) is large. An upper tailed rejection
region seems reasonable for this test. When the true model
is Eq. (3.2), and if a2Al2/k < 2A2 2/(N p), then a lower
tailed rejection region is preferred.
3.3.3. A Method for Locating Optimal Check Points
Given a design for fitting a model of the form in Eq.
(3.1) in a mixture space (note that fixing the design fixes
A2 and (N p)), and given the number of simultaneous check
points desired, k > 1, we now wish to determine where in the
mixture space the k check points should be located so as to
maximize the power of the F test for lack of fit, where the
test statistic is given in Eq. (3.8). We also wish to
position the optimal check points in a manner that is
independent of the values of the elements in 2 .
2
The case of an upper tailed test. To help us find k
simultaneous check points that maximize the power of an
upper tailed test, we shall make use of the fact that the
power is an increasing function of Xi. Therefore to
maximize the power of the upper tailed F test, we shall seek
the locations of the k check points that maximize X1.
As in the case considered in Section 3.2.3, where the
test statistic had a noncentral F distribution, if the
number of extra terms in the true model is p2 = 1, then
maximizing X1 is equivalent to maximizing the scalar A1.
However, as before, if p2 > 1, then the P2xP2 matrix A1 is
not a scalar and we will have to approximate the
maximization of X1 by maximizing a lower bound for X1. This
is done by finding the maximum value of min' the smallest
eigenvalue of A1, since
minm2/22 X1
When the number of check points is less than the order
of the square matrix A1, that is, k < p2, then rank(A1) <
min(k, p2), and A1 will have Umin = 0. For this case, we
again try to maximize the smallest positive eigenvalue of A1
which we denote by in' while remembering from Section
min'
3.2.3 that this technique is limited to situations where
B2 e nC(P1)
21
The case of a lower tailed test. To find k check
points to maximize the power of a lower tailed test, we make
use of the fact that the power of the lower tailed F test
increases as X1 decreases. Then if P2 = 1 and A1 is a
scalar quantity, X1 can be minimized with respect to the k
check points by finding the check points that minimize A1.
If P2 > 1, then by Theorem 12.2.14(9) in Graybill (1969), we
see that an upper bound for X1 is
X1 max_22 /202 (3.11)
where umax is the largest eigenvalue of Al. An approximate
solution to minimizing X1 in (3.11) can be achieved by
minimizing max. It is not necessary to treat the case
'max
of k < p2 separately here, although X1 will equal zero if
_2 is in the column space of P2, where P2 is the matrix
whose columns are orthonormalized eigenvectors corresponding
to the zero eigenvalues of the matrix A .
3.3.4 Determining Whether the Test Is Upper Tailed or Lower
Tailed
The procedures outlined in Section 3.3.3 produce a set
of k check points that simultaneously maximize the power of
the upper tailed test as well as a second set of k check
points that simultaneously maximize the power of the lower
tailed test. The check points that are selected maximize
the power, given A2, k, and N p without specification
of a2' except that when A1 is positive semidefinite we
require that 82 e n C(P ).
It is now necessary to decide which of our two
candidates will be used for a lack of fit test. To choose
between the upper tailed test and the lower tailed test, let
us consider the quantity
R = [A /k] [A2/(N p)].
If R is positive definite when the true model is Eq. (3.2),
then no matter what the value of 82 is, the ratio
E(numerator)/E(denominator) will be greater than unity,
implying an upper tailed test is to be used. Similarly, if
R is negative definite, then a lower tailed test should be
used. Finally, if R is not definite, then neither an upper
nor a lower tailed test is implicated and further
investigation is necessary. The criterion of R = [A /k] 
[A2/(N p)] may yield any of the four following cases.
Case 1. If R = [A /k] [A2/(N p)] is positive
definite when A1 is generated by the k optimal upper tailed
test check points, and R is not negative definite when A1 is
generated by the k optimal lower tailed test check points,
then we recommend that the check points be used that yield
the optimal upper tailed test with an upper tailed rejection
region.
For Case 1 it is necessary for A1 to be positive
definite (see Appendix 4). Since A1 is a square matrix of
order p2xP2 with rank(A ) < min(k, p2), then A1 can be
positive definite only if k > p2. Thus, there must be at
least p2 check points for Case 1 to hold, where p2 is the
number of terms in the model of Eq. (3.2) that are not in
the model of Eq. (3.1).
From inspection of equations (3.9) and (3.10), it is
apparent that the testing for lack of fit in Case 1 is
equivalent to testing the hypothesis
1 2
H0 0 N p = (3.12)
against the alternative
1 2
H > 0
a k N p
60
since R = [A /k] [A2/(Np)] is positive definite when the
true model is Eq. (3.2). In Appendix 5(a) it is shown that
under Case 1, the hypothesis given by (3.12) is equivalent
to the hypothesis
H X = X = 0.
Case 2. In Case 2 we assume that R = [A /k] 
[A2/(N p)] is not positive definite for the k optimal
upper tailed test check points, but that R is negative
definite for the k optimal lower tailed test check points.
Here we recommend that the lower tailed test check points be
used with a lower tailed rejection region.
It is necessary for A2 to be positive definite for Case
2 to occur (see Appendix 4). However, A1 need not be
positive definite, and so k need not be greater than p2. In
Case 2 then, it is possible that lack of fit may be tested
with only one check point.
By inspection of equations (3.9) and (3.10), a
hypothesis of no lack of fit is equivalent to
X1 X
1 2
Hk = 0 (3.13)
0' k N p
while the alternative hypothesis that lack of fit is present
is equivalent to
X1 X
1 2
H N < 0
a k N p
since R = [A /k] [A2/(N p)] is negative definite. In
Appendix 5(b) it is shown that the hypothesis given by
(3.13) is equivalent to the hypothesis
H0: = 2 = 0.
Case 3. We assume R is positive definite for the k
optimal upper tailed test check points, and R is negative
definite for the k optimal lower tailed test check points.
Hence either an upper or lower tailed test may be considered
as a possible test for lack of fit. If the quantity
2
_'_2//o can be specified, then the minimum power for both
the optimal upper and optimal lower tailed tests can be
approximated, and the test with the greater minimum power is
recommended. In Appendix 4 it is shown that Case 3 can
occur only when A1 is positive definite for the upper tailed
test. Thus Case 3 can only occur when there are at least p2
check points.
The minimum power of the upper tailed test may be found
by calculating
P IF" > F ), (3.14)
k,Np;1 IL' 2U a;k,Np
where F;k,Np is the upper 100a percentage point of the
central F distribution,
2
IL = /min2/2a2
and
X2U = max22/2 02
where min is the smallest eigenvalue of A1 and max is the
mmn 1 max
largest eigenvalue of A2. Formula (3.14) yields a
conservative lower bound for the power of the optimal upper
tailed test. Note that A1 is generated using the optimal
upper tailed test check points. The cumulative distribution
function of F" can be approximated by multiplying the
cumulative probabilities of the central F distribution by a
constant (Johnson and Kotz, 1970, p.197). This
approximation is described in Appendix 6. Other
approximations for F" (such as the Edgeworth series
approximation suggested by Mudholkar, Chaubey, and Lin,
1976) exist which are generally more accurate, but we chose
to use the approximation given in Johnson and Kotz (1970,
p.197) due to its simplicity. Additionally, the
approximation of Mudholkar, Chaubey, and Lin (1976) produced
negative probabilties when only one degree of freedom was
available in either the numerator or denominator of F".
This problem was avoided by using the approximation given by
Johnson and Kotz (1970).
The minimum power of the optimal lower tailed test can
be approximated similarly (if B~22/o2 is specified) by
calculating
P IFF"
< (la);k,Npp
where
X lu = maxJ.2 2/2o 2
and
S2L = 6min22/20 2
with Umax equal to the largest eigenvalue of A1 and 6min
max n
equal to the smallest eigenvalue of A2. Note that A1 is
generated by using the optimal lower tailed test check
points. For the lower tailed test, A1 may be positive semi
definite, and if 82 is in the column space of P2 then A1 = 0.
In Case 3, the upper tailed test is a test of
HO: 1 = = 0
0 1 2
H 2 > 0
a k N p
while the lower tailed test is a test of
H: X =X2 0
H0 1 2 = 2
1 2
H P < 0.
a k N P < 0
Case 4. In Case 4 we assume that R = [Al/k] 
[A2/(N p)] is not positive definite for the k optimal
upper tailed test check points and R is not negative
definite for the k optimal lower tailed test check points.
Here it is useful to write the difference between the
expected value of the numerator and the expected value of
the denominator of the F ratio in (3.8) as
al[A /k A2/(N P)] = sns'S
2 1 2 2 = 2
= 8[S1:S2:S3] diag[l2'r2=0,03 [Sl:S2:S3] '2
= 8_2S"l I2 + 2S3 332
where 0 = diag(ll, 02' 23) is a diagonal matrix consisting
of the eigenvalues of R, 01 is a diagonal matrix of the
positive eigenvalues of R, 02 is a diagonal matrix of the
zero eigenvalues of R, and 03 is a diagonal matrix of the
negative eigenvalues of R. The orthogonal matrix S can be
expressed as S = [S1:S2:S3], where the matrices SI, S2, and
S3 have columns which are orthonormalized eigenvectors
corresponding to nl, 02, and 03, respectively.
In Case 4, neither the optimal upper tailed test nor
the optimal lower tailed test is applicable for all values
of _2 For completeness, we note that Case 4 actually
consists of nine subcases, where R may be positive semi
definite, negative semidefinite, or indefinite for either
the optimal upper tailed test or lower tailed test check
points. These subcases are listed in Table 2.
Table 2. Nine Subcases of Case 4.
RUpper RLower
Subcase Tailed Test Tailed Test
1 PSD PSD
2 PSD NSD
3 PSD I
4 NSD PSD
5 NSD NSD
6 NSD I
7 I PSD
8 I NSD
9 I I
PSD = positive semidefinite, NSD = negative semi
definite, I = indefinite.
If _2 lies in the column space of S2, then 8'[A /k 
A2/(N p)]s82 is zero, and therefore lack of fit is not
testable with either an upper or lower tailed test. A
sufficient condition for the test for lack of fit to be
upper tailed in Case 4 is that 0 be in the column space
2
of [S1:S2], but not entirely in the column space of S2. In
this case
;[A1/k A2/(N p)]_2 = 2S01812 + 2S303S32
= 2S1~lS 2 + 0
= S~Sl2lSI2,
and 8_[ Al/k 
indicating an
condition for
that _2 be in
in the column
A2/(N )]_2 will be greater than zero,
upper tailed test. Similarly, a sufficient
the test for lack of fit to be lower tailed is
the column space of [S2:S3], but not entirely
space of S2. Then
2[A1/k A2/(N P)]g2 = 0 + 2S33S3_2
= 2S3 3S 2
which makes _2[Al/k A2/(N p)] 2 less than zero,
indicating a lower tailed test.
To determine whether 2 is in the column space of
[S :S2], let us define the augmented matrix
Q1 = [f2:Sl:S2] If Q!Q1 has a zero eigenvalue, then 02 is
in the column space of [S :S2]. Similarly, if we define
Q2 = [82:S2] and Q3 = [82:S2:S3]' then a2 is in the column
space of S2 if Q'Q2 has a zero eigenvalue, and 2 is in the
column space of [S2:S3] if Q3Q3 has a zero eigenvalue.
Given that we are in a particular subcase of the nine
subcases described in Table 2, we recommend that lack of fit
be tested with the upper tailed test check points if it is
determined that 2 is such that '[A /k A2/(N p)]2 is
positive when A1 is generated from the upper tailed test
check points. Likewise, for the same given subcase, if the
value of 02 of interest is determined to produce a negative
value for i2[Al/k A2/(N P)]@2 when A1 is generated from
the lower tailed test check points, then we recommend that
lack of fit be tested with the lower tailed test.
We see then that Case 4 is an undesirable situation in
practice, since, in order to test for lack of fit, we must
assume a priori that any lack of fit is due to a nonzero
value of 82 that produces an upper tailed or lower tailed
rejection region. However, it would seem rare that such
knowledge would be available.
3.4 Examples
We now present several examples to illustrate the
technique for locating optimal check points to be used in
testing for lack of fit in a mixture model.
3.4.1 Theoretical Examples
Example 1. In this example a second order canonical
polynomial model is fitted in three mixture components using
the {3,2} simplex lattice design, which is presented in
Figure 1 of Chapter 1. The true model is assumed to be the
special cubic model containing the term 123x x 23 in
addition to the six terms of the fitted model. The expected
values of the response at the six design points are assumed
to be represented by the fitted model in the form
E(Y) = X01,
but with the true model the expectations are written as
E(Y) = X 1 + X22 ,
where X is a 6x6 matrix with rows that define the
coordinates of the six design points and columns that
correspond to the six terms in the fitted model (xi, x2, x3,
xlx2, xlx3' x2x3)' 8a is the 6x1 vector of regression
coefficients (81, 82, 83, 812, 813, 823), X2 is a 6x1
column vector containing the values of the term xlx2x3 at
the design points, and 82 is the single regression
coefficient 8123'
The {3,2} simplex lattice design consists of only six
design points, and since six parameters are estimated in the
second order fitted model, there are no degrees of freedom
remaining for obtaining an estimate of the experimental
error, 02. We assume therefore that an external estimate of
a is available, a2 which will be used in the denominator
of the lack of fit F statistic given in Eq. (3.3).
Since there is one term in the true model in addition
to those in the fitted model, that is p2 = 1, we know that
in order to locate k simultaneous check points that maximize
the power of the test for lack of fit it is necessary to
maximize the scalar quantity
1
A = (X* X*A)'V (X* X*A)
1 2 0 2
with respect to the coordinates of the k check points. Here
X* is a kelement column vector with ith element equal to
the value of x* x* x* at the ith check point, X* is a kx6
il i2 i3
matrix with ith row equal to the value of (x* x* x*
il 12' i3'
S x* x* x xt ) at the ith check point,
11 1i2' 11 Xi3' 12 i3
A = (X'X)X'X2 is the 6x1 alias vector, and
V = Ik + X*(X'X) X*'. This maximization is accomplished
by use of the Controlled Random Search Procedure (Price,
1977), which is described in Appendix 2.
When only a single (k = 1) optimal check point is
desired the Controlled Random Search Procedure locates a
point (x*, x*) which maximizes
1
A1 = (X* X*A)'V0 (X2 X*A),
where
S= x*x*x* = *x*( x* 
2 123 12 1 2
X* = (x*, x*, x* x* x*x*, x*x*)
= (x*, X* (1 x* x*), x*x*, x*(l *),
1 2 1 2 12' 1l 1 2
x*(l x* x*)),
2 1 2
l
and V0 = 1 + X*(X'X) *'. The value of A1 is calculated
using the formula of Eq. (3.6). Following this procedure,
we find that the single check point that maximizes A1, and
thus maximizes the power of the test, is the centroid of the
triangular factor space (1/3, 1/3, 1/3). The value of A1 at
this centroid point is A1 = 0.00084.
When the Controlled Random Search Procedure is used to
locate k = 2 simultaneous check points that maximize A,, the
centroid (1/3, 1/3, 1/3) is selected twice, and A1 =
0.00121. For three simultaneous optimal check points, the
centroid is selected three times, and A1 = 0.00142.
To test whether the second order model exhibits lack of
fit, when we suspect the special cubic model is the true
model, we form the F ratio
1
d'V d/k
0
F = ^ 2
^2
ext
with the single check point (1/3, 1/3, 1/3) where d =
Y* Y*, Y* is the observed response, Y* is the response
1 1 1 1
predicted by the second order fitted model at (1/3, 1/3,
1/3), and V0 = 1 + (1/3, 1/3, 1/3, 1/9, 1/9, 1/9)(X'X)1
(1/3, 1/3, 1/3, 1/9, 1/9, 1/9)'. If the calculated value of
the F ratio exceeds F where v equals the number of
degrees of freedom associated with a then we reject the
ext
null hypothesis that the second order model is the true
model in favor of the alternative hypothesis that the
special cubic model is the true model. Equivalently, we
reject H: I1 = 0 in favor of Ha: 1 > 0. For k = 2 or
k = 3 check points, the value of the F ratio is calculated
using the observed and predicted responses at the two or
three replicates at the centroid. The hypothesis
H0: X = 0 is rejected in favor of Ha: X > 0 if F
exceeds F;kv
a ;k,v
Example 2. In Example 2 we illustrate the second of
the four cases that could arise when MSE is used as an
estimate of 02 in the lack of fit test statistic (see
Section 3.3.4). We again fit a second order canonical
polynomial model in three mixture components, and assume the
true model is special cubic. The design to be used is the
q = 3 simplex centroid design, which consists of seven
design points, and is illustrated in Figure 2 of Chapter 1.
There are six parameters to be estimated and seven
design points hence one degree of freedom can be used to
calculate MSE. We shall use MSE to estimate a2. Optimal
upper and lower tailed test check points must be located,
and then a decision is made as to which test should
be used. The actual testing for lack of fit involves the F
statistic in (3.8).
As in Example 1, p2 = 1, since there is one term in the
true model in addition to those in the fitted model. Thus
Al is a scalar whose value we seek to optimize with respect
to the desired number of check points, k. When only a
single check point is sought for the purpose of testing lack
of fit, the Controlled Random Search Procedure has two
functions. First, the procedure is used to locate the
optimal candidate check point for an upper tailed test by
locating the check point that maximizes the scalar Al.
Secondly, the procedure is used to locate the optimal
candidate check point for a lower tailed test, which is
accomplished by locating the point that minimizes A1. The
quantity R = [Al/k] [A2/(N p)] is then calculated to
determine whether the upper or lower tailed test will be
used. If R is positive for the candidate check point for an
upper tailed test, then the test is upper tailed, and the
test is lower tailed if the candidate check point for a
lower tailed test produces a negative value for R. Note
that A2 = (X2 XA)'(X2 XA) is fixed once the design is
specified, since A2 does not depend on the check points.
Using the Controlled Random Search Procedure it is found
that the maximum value of A1 occurs at (xl, X*, x*) = (1/3,
1/3, 1/3), which will be the location for the check point
for the upper tailed test. Calculating A1 at this centroid
point, we find that R = [A /k] [A2/(N p)] = [(3.7258
x 104)/l] [(8.4175 x 104)/l] = 4.6917 x 104. Since R
is negative, the test is not upper tailed.
Using the Controlled Random Search Procedure to
minimize AI, we find that a subregion of the factor space
exists in which all points yield a near minimum value for
A1. We choose the point (0.0189, 0.9269, 0.0542) at random
from this subregion to be used as the optimal candidate for
a lower tailed test. Here R = 0 [(8.4175 x 104)/] =
8.4175 x 104.
Since R is negative for both the optimal upper tailed
test check point and for the optimal lower tailed test check
point, we have Case 2 of Section 3.3.4. The upper tailed
test check point is disregarded, and the lower tailed test
check point (0.0189, 0.9269, 0.0542) is used to test for
lack of fit. If the calculated F ratio,
1
d'V d
F =
MSE
is less than F ( );,then H: X = X = 0 is rejected in
favor of Ha: [X1/1] [X2/1] < 0, that is we conclude that
the second order model exhibits lack of fit, and the true
model is special cubic.
When two simultaneous check points are desired for
testing lack of fit, we can again use the Controlled Random
Search Procedure to locate the optimal settings. To
maximize the scalar A,, we find that both check points
should be selected at (1/3, 1/3, 1/3), for an upper tailed
test. With our calculations R = [(5.8275 x 104)/2] 
[(8.4175 x 104)/1] = 5.5038 x 104, but since R is
negative, the test is not upper tailed.
Minimizing A1 to locate two optimal lower tailed test
check points yields a subregion in the factor space of
optimal check points. The pair of check points (0.3749,
0.5752, 0.0499) and (0.5332, 0.4169, 0.0499) is selected at
random from this subregion, and these check points yield
R = 0 [(8.4175 x 104)/1] = 8.4175 x 104.
Since R is negative for the upper tailed test points
and the lower tailed test points, we have Case 2 of Section
3.3.4 again and the lower tailed test check points are used
to test for lack of fit. The hypothesis H0: X1 = 2 = 0 is
rejected in favor of Ha : [1/2] [x2/1] < 0 if the cal
1
culated value of F = (d'Vo d/2)/MSE is less than
F(1a);2,1, in which case we say lack of fit of the model is
present.
^2
If an external estimate aext had been available for
ext
this example, then the optimal upper tailed test check
points could have been used in the F ratio,
1 ^2
F = (d'V d/k)/o t and lack of fit would then be detected
0 k exta
if the calculated value of F exceeded F ,
a;k,v
Example 3. Example 3 illustrates the procedure for
locating optimal check points when there are two terms in
the true model in addition to those in the fitted model. A
second order canonical polynomial model in three mixture
components is fitted using a q = 3 simplex centroid design.
The true model is assumed to contain eight terms, six of
which are the same terms as in the fitted model, with the
additional two terms being the third order terms
612x x2(X1 x2) and B123x x2x3. As in Example 2, there is
one degree of freedom for MSE which is used to estimate
2 T
2. The test statistic, F = (d'V d/k)/MSE, is given in
equation (3.8).
Since p2 = 2 and A1 is a 2x2 matrix, locating the
optimal upper tailed test check points by the procedure of
maximizing X1 is assisted by the maximizing of a lower bound
for 1 namely maximizing u min /2a where u is the
I min22 mm
smallest eigenvalue of A1. Since 82 and a2 are unknown,
this is equivalent to maximizing min For min to exceed
zero, it is necessary that A1 be of full rank, and since
rank(A1) ( min(k, p2), it is necessary to select k > 2 check
points. If A1 is less than full rank, and thus is positive
semidefinite, only a subset of possible values of 2 could
be considered to make it possible to test for lack of fit
with an upper tailed test.
Using the Controlled Random Search Procedure, the
points that maximize min are found to be (0.418, 0.277,
mln
0.305) and (0.277, 0.418, 0.305). These points are thus
optimal candidates for upper tailed test check points. At
these check points we have umin = 5.1623 x 104, A1 =
diag[5.1623 x 104, 5.1916 x 104], A2 = diag[0, 8.4175 x
104], and R = [A1/2] [A2/I] = diag[2.5811 x 104, 5.8217
76
x 104]. Since the eigenvalues of R are 5.8217 x 104 and
2.5811 x 104, R is indefinite. Following the suggested
procedure for Case 4 of Section 3.3.4, we note that an upper
tailed test for lack of fit exists if the value of
'2 = [512' 8123]' is in the column space of [Sl:S2] but not
entirely in the column space of S2, where S1 is the matrix
whose columns are the orthonormalized eigenvectors of R
corresponding to the positive eigenvalues of R, and S2 is
the matrix whose columns are the orthonormalized
eigenvectors of R corresponding to the zero eigenvalues of
R. Since R has no zero eigenvalues in this example, S2 does
not exist, but S1 is the column vector, S1 = [1,0]'. Thus
if 2 is of the form 2 = [612 0]', where 12 0, then
a2 is in the column space of S1 and the test is upper
tailed.
The matrix A2 has rank one and therefore is positive
semidefinite. Hence it is impossible to locate two check
points that minimize umax and also make R = [A1/2] [A2/1]
negative definite (see Appendix 4), that is, it is
impossible to find a lower tailed test that is capable of
testing lack of fit for all values of 2. However, if we
use the Controlled Random Search Procedure to locate two
check points that minimize an upper bound for X1 which is
umax S2 /2 2, then by minimizing max, we find that any of
the check points in a particular subregion of the factor
space yield a near minimum for umax. One pair of points in
this subregion is selected as the points to be used as
optimal lower tailed test check points, namely the pair
consisting of the point (0.053, 0, 0.947) replicated twice.
Replicating this check point, we find max = 7.3900
x 1011, A = diag[0, 7.3900 x 1011], A2 = diag[0, 8.4175 x
104], and R = [A1/2] [A2/1] = diag[0, 8.4175 x 104].
The eigenvalues of R are 0 and 8.4175 x 104 implying that
R is negative semidefinite. The values of 8 that are in
2
the column space of [S2:S31 but not entirely in the column
space of S2 will provide a lower tailed test. Here, [S2:S3]
= diag[l,l] and S2 = [11,0]'. Thus, the test for lack of fit
is lower tailed if 8123 0.
For values of 8 that produce an upper tailed test we
2
use the check points (0.418, 0.277, 0.305) and (0.277,
0.418, 0.305) with the F ratio
1
d'Vo d/2
F =
MSE
and conclude there is lack of fit if the calculated value of
F exceeds Fa;21. For values of a2 that produce a lower
tailed test, we use two replicates of the check point
(0.053, 0, 0.947), and'conclude there is lack of fit if F is
less than F(1a);2,1, where again F is calculated by
1
F = (d'V d/2)/MSE.
0O
Example 4. Example 4 illustrates Case 3 of Section
3.3.4 in which MSE is used to estimate 02 in the lack of fit
test statistic. A second order canonical polynomial model
in three mixture components is fitted using the {3,3}
simplex lattice design, which appears in Figure 5. The true
model is assumed to be special cubic, thus p2 = 1 and A1 is
a scalar. The {3,3} design consists of ten design points
and since there are six parameters to be estimated in the
fitted model, 02 can be estimated by MSE with N p = 10 6
= 4 degrees of freedom.
We first suppose that a single check point is to be
used to test for lack of fit. Using the Controlled Random
Search Procedure we find the single check point that
maximizes the scalar
1
A= (X X*A)'V (X* X*A)
1 2 0 2
is located at the centroid of the simplex factor space.
Thus (x*, x*, x*) = (1/3, 1/3, 1/3) is the optimal candidate
for an upper tailed test check point. At this centroid
point, A1 = 4.9076 x 104. For the {3,3} design the scalar
quantity A2 = (X2 XA)'(X2 XA) is fixed and is equal to
A2 = 9.4062 x 104 and thus, R = [Al/k] [A2/(N p)] =
[(4.9076 x 104)/1] [(9.4062 x 104)/4] = 2.5560 x 104
The point that is the optimal candidate for a lower
tailed test check point is chosen randomly from a subregion
of points in the factor space, in which all points minimize
A1. The point selected has the value (x*, x*, x*) = (0.560,
0.410, 0.030). Here A1 = 9.6590 x 107 and R = [(9.6590 x
107)/1] [(9.4062 x 104)/4] = 2.3419 x 104.
(1,0,0)
3 3 a0o 0'x
( 0)/ 3 3
3' 3'3
(0o,,0)  (0,0,1)
x2 0, ) 3 3' ) x3= 1
Figure 5. The {3,3} simplex lattice design.
Since R is positive for the optimal upper tailed test
check point (1/3, 1/3, 1/3) and R is negative for the
optimal lower tailed test check point (0.560, 0.410, 0.030)
we are in Case 3 of Section 3.3.4. Either the upper or
lower tailed test could be used to test for lack of fit, but
if the quantity B B2/o2 can be specified, then we will
choose to use the test that has greater minimum power, since
greater power means that we are more likely to detect lack
of fit when in fact lack of fit exists. In this example
2 123*
For illustrative purposes, we arbitrarily choose
2
8'8 /C2 = 2000, so that an approximate conservative lower
bound for the power of the upper tailed test is found by
calculating
P {F" > F k
k,Np;XIL' 2U a;k,Np
where F;k,Np is the upper 100a percentage point of the
central F distribution, k is the number of check points, N
is the total number of response observations, p is the
number of parameters in the fitted model,
2 2
IlL = UminB2/2a and 2U 6 axI/20 The
quantity umin is the smallest eigenvalue of AI, where A1 is
evaluated at the optimal upper tailed test check point.
Since A is a scalar, mn = A. Likewise, 6 is the
1 mim 1 fmax
largest eigenvalue of A2, and since in this example A2 is a
scalar, 6max = A2. In this example we have k = 1, N p =
10 6 = 4, AlL = UminS2/2a2 = (4.9076 x 104)(2000/2)
1 2
= 4.9076 x 101, and 2U = 6 m /2a2
= (9.4062 x 10 )(2000/2) = 9.4062 x 101. Using the
approximation to the cumulative probabilities of the doubly
noncentral F distribution given by Johnson and Kotz (1970,
p.197) which is described in Appendix 6, and taking a = .05,
we find that a conservative lower bound for the power of the
optimal upper tailed test is approximately equal to .0649.
The minimum power for the optimal lower tailed test is
2
approximated (assuming 8_Y2/a = 2000) by calculating
P F" < FN
k,Np;lU' 2L (1a);k,Np
2
The quantities A1U and 2L are taken as A U = maxJ2 /2o
lJU 2L lU maxz22
and 2L = 6 mn /2o where max is the largest eigenvalue
2L min22 max
of A1 with A1 calculated using the optimal lower tailed test
check point, and where 6 in is the smallest eigenvalue of
A2. Since A1 and A2 are scalars, max = A and 6 min = A .
2 1 2 max 1 mm 2
In this example, k = 1, N p = 4,
AlU = (9.6590 x 10 )(2000/2) = 9.6590 x 104, and
4 i
x2L = (9.4062 x 10 )(2000/2) = 9.4062 x 101 Again if the
approximation to the doubly noncentral F distribution given
in Johnson and Kotz is used, an approximate conservative
lower bound for the power of the optimal lower tailed test
is .0555.
Having specified a 2/02 = 2000, the optimal upper
tailed test is chosen over the optimal lower tailed test,
because the approximate minimum power of the upper tailed
test is greater than the approximate minimum power of the
lower tailed test. Using the optimal upper tailed test
check point (1/3, 1/3, 1/3) in the test statistic
1
d'V d
F =
MSE
we conclude that lack of fit is significant if the
calculated value of F exceeds F a1,4 in which case we
a;1,4'
reject H0: XI = 2 = 0 in favor of Ha: A/l 12/4 > 0.
When two simultaneous check points are used for testing
lack of fit, the Controlled Random Search Procedure locates
the optimal upper tailed test and optimal lower tailed test
check points. It turns out that two replicates at (1/3,
1/3, 1/3) maximize Al, and are used as optimal check points
for an upper tailed test. The value of R = [A,/2] [A2/4]
is [(7.9210 x 104)/2] [(9.4062 x 104)/4] = 1.6090 x
104.
In searching for two optimal lower tailed test check
points, again a subregion of the factor space is found in
which any of the points nearly minimize A From this
subregion are chosen the points (0.6386, 0.3263, 0.0351) and
(0.7257, 0.2421, 0.0322) resulting in a value of R = [A,/2]
 [A2/4] of [(1.5216 x 109)/2] [(9.4062 x 104)/4] =
2.3516 x 104.
In conclusion, when two simultaneous check points are
used in the test for lack of fit in this example, R is
positive for the optimal upper tailed test and R is negative
for the optimal lower tailed test, and we have Case 3 of
2
Section 3.3.4. Selecting 002 /o = 2000 arbitrarily, we
found the approximate lower bound for the power of the upper
tailed test to be .0504, and the approximate lower bound for
the power of the lower tailed test to be .0612. Since the
power is higher with the lower tailed test it is our choice
for testing lack of fit when two check points are used
simultaneously. Lack of fit is detected and we reject
H0: 1 2 = 0 in favor of Ha: [1 /2] [ 2/4] < 0 if the F
1
ratio, F = (d'V0 d/2)/MSE, using the optimal lower tailed
test check points (0.6386, 0.3263, 0.0351) and (0.7257,
0.2421, 0.0322) is calculated to be less than F1a);2,4'
(ia);2 ,4"
3.4.2 Numerical Examples
Numerical Example 1. In this example we illustrate
numerically some of the findings in the first theoretical
example of Section 3.4.1. Data that were collected in a
rocket fuel experiment (Kurotori, 1966) will be used to
investigate the power of the lack of fit F test. The test
is set up to detect the inadequacy of a fitted second order
canonical polynomial model when the true model is special
cubic. Calculated values of the power of the test which
detects lack of fit through large values of
1
d'V0 d/k
F =
^2
ext
will be compared for several check point locations, includ
ing the location (1/3, 1/3, 1/3) at which the power was
found to be maximum in Example 1 of Section 3.4.1.
In Kurotori's experiment the modulus of elasticity (Y)
of a rocket fuel is expressed as a function of the
proportions of three componentsbinder (xl), oxidizer (x2),
and fuel (x3). Since lower bounds are placed on the
component proportions xl, x2, and x3, in the form of
0.20 < xl, 0.40 < x2, and 0.20 < x3, pseudocomponents (x!)
are defined in terms of the original components in the form
of x1 = (xl 0.20)/(1 .80), x' = (x2 0.40)/(1 .80),
and x' = (x3 0.20)/(1 .80). The true special cubic
model in the pseudocomponents, which is obtained by using
the data at the seven points of the simplex centroid design
in the pseudocomponent system, is
E(Y) = 2350x' + 2450x' + 2650x' + Ox'x'
1 2 3 1 2
+ lOOOx'x3 + 1600x'x' + 6150x'x'x'.
The second order canonical polynomial model that is fitted
to the six boundary points only, and which will be tested
for lack of fit, is given by
Y = 2350x' + 2450x' + 2650x'
3
+ 1000x'x' + 1600xNx'.
The configuration of the experimental design as well as the
check point locations are depicted in Figure 4 of Chapter 2
and the observed response values are given in Table 3 of
this chapter.
1 ^2
A value of the ratio F = [d'V d]/o is calculated at
0 ext
each of the four individual check points (1/3, 1/3, 1/3),
(2/3, 1/6, 1/6), (1/6, 2/3, 1/6), and (1/6, 1/6, 2/3)
"2 ^2
where ext is assumed to have the value et = 144 as
suggested by Kurotori (1966). We also assume without loss
of generality that the degrees of freedom associated
^2
with aet are v = 10. The power of the F test is calculated
at each of the four check points by using the approximation
to the cumulative probabilities of the noncentral F
Table 3. Observed Response Values at the
Pseudocomponent Settings for Kurotori's Rocket Fuel
ExperimentNumerical Example 1.
Observation
Number x'
______ _____
1
2
3
4
5
6
7*
8*
9*
10*
Binder Oxidizer Fuel
Modulus
x 3
0
1/2
1/2
0
1/3
2/3
1/6
1/6
0
1
0
1/2
0
1/2
1/3
1/6
2/3
1/6
0
0
1
0
1/2
1/2
1/3
1/6
1/6
2/3
of Elasticity
Y
2350
2450
2650
2400
2750
2950
3000
2690
2770
2980
*Check Points.
distribution given by Johnson and Kotz (1970, p. 197) to
evaluate
Power = P{FI' ,;,
11, 10;Xi
F.05;1,10}
2 ^2
where 1 = A 8 23/2 ext.
I 1 123 ex
The value of
A = (X* X*A)'V0 (X* X*A) is calculated for each check
point using the values of X*, X*, v and the value of
A = (X'X) X'X2 which is fixed by the {3,2} simplex lattice
design. Since the {3,2} simplex lattice consists of points
only on the boundaries of the triangle (and therefore at
each point at least one of the x! values is equal to zero),
1
then X2 = 0 and A = 0. From the true special cubic model,
123 = 6150.
The calculated value of F as well as the approximate
value for the power at each of the four check points is
given in Table 4. The check point (1/3, 1/3, 1/3) produced
the highest power of the four check points investigated,
supporting the previous results of Example 1 in Section
3.4.1 where (1/3, 1/3, 1/3) was selected as the check point
location with the maximum power when a second order
canonical polynomial was fitted using the {3,2} simplex
lattice design, but the true model was assumed to be special
cubic. Additional support for the point (1/3, 1/3, 1/3)
being optimal is given by the contour plot of values of A1
in Figure 6(d). The highest values of A1 appear near the
centroid (1/3, 1/3, 1/3) where high A1 values translate into
2 2
high X1 values, since Xi = A a123/22 which in turn implies
high power since we know the power is an increasing function
of X1.
As a second part of this example the power of the F
test that is obtained when three replicates are taken at
(1/3, 1/3, 1/3) is compared to the power of the F test that
is obtained when one replicate is taken at the three check
points (2/3, 1/6, 1/6), (1/6, 2/3, 1/6), and (1/6, 1/6,
2/3). These latter three point locations were suggested by
Kurotori for testing lack of fit of his fitted special
o 0 0 0
4
r1 4 r 4
, 4 , O
OC
11
mm
M N\N
i '
 '
0 *
0w
O
.JU
Q .H
3
Q)
4J)
U
4
0
U
0
m
4a
4J
(0
>1
cu
4Ja
CO
00
S0O
Cu
A 0
u r
0 0
c
o
4 r4
UO4
A U( 
S*
I 4 4 J
U K M
14
o
o
o
o
rI CN
, q 
CN .4 ,I
(m)
cMNrN
S S
0 
4J 4
Z 1
cubic model. The result of this comparison, see Table 4, is
that the three replicates at (1/3, 1/3, 1/3) produce the
test with greater power which again supports the findings of
Example 1 of Section 3.4.1.
All of the check point locations listed in Table 4
produce very high power values (> .999) which is due in part
to the high value of 123 (123 = 6150). If S1 were of
123 123 6123
lower magnitude, then the three replicates at (1/3, 1/3,
1/3) would show a still greater superiority in the power
value compared to the power using the other check points.
This superiority is demonstrated in Table 5 where values of
8123 are listed as 3000 and 1500 and the comparative power
values are listed as 0.998 compared to 0.795 and 0.662
compared to 0.249, respectively. Table 5 also demonstrates
the superior power value for the point (1/3, 1/3, 1/3) when
123 = 3000 or 123 = 1500 and each of the four check points
is used one at a time.
Finally, (1/3, 1/3, 1/3) being the optimal check point
location is seen in Figure 6(c), where contour plots of the
expected difference in the heights of the surfaces are
drawn. The differences in the heights are found by
subtracting the estimated height of the surface obtained
with the fitted second order model from the estimated height
of the surface obtained with the true special cubic model.
The expected difference between the true and fitted surfaces
approaches a maximum the closer one moves to the centroid of
the simplex factor space, so that the optimal check point
89
<4 0O CD CD C
a) m Co C MN C N c
C 0 *
U a4
I COD0 CO
U co co 1 o CD m
0
a4
41J 4
C a) m co o o co m a) 4)
0 (a D aN O Oo *10 10 m r En
C >1* *
0
U a)
0 a U1 m a) 0T co
m 2' 0 0
'0 Z O*o r * r* *< u
Z ; 0 O
Oi 0
*q 4J V
4 L 0 4 0 *
q m o .m .
0 )
En *' m
a) 0
0 Cr rn)Cr)Cr)
^ N. N. N.. N. 'NN.'N N
a) N N N N NNN N.N U
; 4 (N i 4 C '4 M 4 4
U .4'
x :I
I
X =1 x= I
2 3
(a) True special cubic surface.
X =I
x :1 x :1
2 3
(c) Expected difference between the
true special cubic surface and
the fitted second order surface.
x2
(b) Fitted second order surface.
x =1 X :1
2 3
(d) A( X X XA)' V'( X X*A)
Figure 6. Contour plots for Numerical Example 1.
location (1/3, 1/3, 1/3) coincides with the point where the
expected difference between the true special cubic surface
and the fitted second order surface is maximum.
Numerical Example 2. In this second numerical example,
we investigate the power of the F test for detecting lack of
fit when a second order canonical polynomial model is fitted
in a mixture system which is in truth represented by a
special cubic model. The true model is assumed to be
E(Y) = 2350x1 + 2450x2 + 2650x3
+ 1000x x3 + 1600x2x3 + 6150x 2X3
which is used to generate hypothetical response observations
at the seven points of the q = 3 simplex centroid design as
well as at three check points. The values of the response
are obtained by adding the value of a pseudorandom normal
variate with mean 0 and variance 144 to each true predicted
response value. The data are given in Table 6.
The response values at the seven points of the simplex
centroid design are used in the least squares normal
equations to obtain the fitted second order model
Y = 2341x1 + 2438x2 + 2630x3
+ 310x1x2 + 1304x x3 + 1970x2x3
Table 6. Generated Response ValuesNumerical Example 2.
xl x2 x3 Y
1 0 0 2357
0 1 0 2454
0 0 1 2646
1/2 1/2 0 2403
1/2 0 1/2 2747
0 1/2 1/2 2962
1/3 1/3 1/3 3013
1/3 1/3 1/3 2993
2/3 1/6 1/6 2693
.02 .93 .05 2550
*Check points.
which is to be tested for lack of fit using the test
1
statistic F = d'V0 d/MSE. The F statistic will be evaluated
at each of the three check points (1/3, 1/3, 1/3),
(2/3, 1/6, 1/6), and (.02, .93, .05), taken one at a time,
and the power of the test at the three check point locations
will be calculated and compared. The test is lower tailed
for all check point locations (since R = A1 A2 is negative
for all check point locations) and thus the power is defined
as
PI" ( F;, }.
1,1;XIX2 .95;1,1
