- Email: [email protected]
Response Surfaces and Models
255
Chapter 16 Response Surfaces and Models 1. Introduction
The concepts developed in this chapter and the next three chapters are general and can be applied not only to analytical chemical systems (e.g. electrochemistry, spectroscopy, chromatography) but also to the systems which are being measured by the analytical chemist (e.g. pharmaceutical manufacturing, polymer production, water quality). While these concepts offer valuable insights for improving and understanding the responses from systems, it is well to keep in mind that we might never be able to discover the true behavior of a system. Uncertainties associated with the experimental measurements and with the mathematical description of the system will almost always exist. Thus, we will usually have a more or less “fuzzy” picture of how the system truly behaves. It follows that if we are to have a sufficiently clear view of how the system behaves, we must maximize our fundamental understanding of the system and minimize the uncertainties associated with our experimental measurements. These dual goals are achieved by the simultaneous application of precise measurement techniques, good chemical models, and proper experimental design.
2. Response surfaces A response surface is the graph of a system response plotted as a function of one or more of the system factors (i.e. “independent variables”, see Chap. 13). Response surfaces offer the chemometrician a convenient means of visualizing how various factors affect his or her measurement system. Consider the familiar analytical example shown in Fig. 1, a photometric method for the classical determination of manganese as permanganate. The amount of light transmitted through the system depends upon the concentration of MnO; in the sample. 2.1 True response surfaces
According to Beer’s law (a mathematical model), if the intensity of the light source is absolutely stable, if the light is monochromatic, if there is no stray light in the system, if the MnOi does not undergo chemical or physico-chemical changes as a function of concentration, and if there are no other absorbing species present in solution, then the intensity of transmitted light is given by
z = lox Reference p. 269
256
Light source
Photometric system
tntensi ty
Fig. 1. Systems theory diagram of a photometric method for the classical determination of manganese as permanganate.
and the transmittance, T, is given by
where I is the intensity or power of the transmitted light, I, is the intensity or
Fig. 2. Scaled photometric intensity as a function of the concentration of permanganate in a sample.
251
power of the light source, c is the molar absorptivity in L mol-’ cm-’,b is the path length through the sample in cm, and C is the concentration of MnOc in the sample in mol I-’. In this system, the intensity of transmitted light, I, is the system response and the concentration of MnO; in the sample, C, is the system factor. The true response surface (i.e. the graph of the system response plotted as a function of the single system factor) is shown in Fig. 2. Note that, for this true response surface, the physical understanding is exact and there is no experimental uncertainty. 2.2 Measured response surfaces Although many system diagrams are drawn in a manner similar to Fig. 1, such diagrams are incomplete because they do not explicitly take into account a separate but necessary measurement system. This is a subtle but important point. In the photometric determination of MnO; , for example, transmitted intensity is the important physical property of interest, but we seldom measure it directly. Instead, with modern photometric instruments, transmitted intensity is converted to an electrical current, then to a voltage, and finally to an analog or digital meter readout from which we visually “observe” the value of the transmitted intensity. A correct view of the double-beam photometric measurement system is given in Fig. 3. It is clear that the measurement system itself can introduce a set of uncertainties and biases quite separate from those introduced by the photometric system. For example, if electrical line voltage fluctuations affect the measurement system in such a way that the readout fluctuates, the originally well-defined response surface of Fig. 2 becomes uncertain again. In general, measured response surfaces are ‘‘fuzzy” because of uncertainties in the primary system itself and because of uncertainties introduced by the associated measurement system. It is usually desirable that this latter source of uncertainty should not
I T
0 f
Light source
Photometric
Measuremeo t
Readout
Fig. 3. Systems theory diagram showing the relationship between the photometric system and the measurement system for a double-beam photometer.
Reference p. 269
258
be greater than the former; hence the desire for precise (i.e. low uncertainty) analytical chemical methods. 2.3 Estimated response surfaces and models
Consider the set of experimental points shown in Fig. 6. It is probably unrealistic to try to write a mechanistic model (see Chap. 12) to explain these data. Our knowledge of metallurgy would be quickly exceeded. It is possible, however, to propose an empirical model that might provide a good description of the data. One possible empirical model might be the equation of a parabola
+
y l i = b, + b,x,, + b,,xfi eli (3) where y I i represents the system response (tensile strength) for experiment i, x l i represents the system factor (percent manganese in this case), b, is an offset parameter, b, is a first-order parameter, b,, is a second-order parameter, and el, is a residual or deviation between what is observed and what is predicted by the model. This is a single-factor model of a single-factor response surface. Matrix least-squares techniques (see Chap. 13) can be used to fit this model to the data shown in Fig. 6. As an example, suppose the data in Fig. 6 have the following numerical values (D is an “experimental design matrix” that has a number of rows equal to the number of experiments and a number of columns equal to the number of factors; each element of the D matrix thus represents the value of a particular factor in a particular experiment). 3 58 4 87 93 4 4 91 105 5 D= Y= 98 5 6 88 6 90 6 89 761
-.
Iron or Mongonese
system
meosuremen t
Tensile -strength meosuremen t
% Mn
Tensile strength
Fig. 4. Systems theory diagram showing the relationships among a steel-making system and two measurement systems.
259
Fig. 5. Tensile strength as a function of the percent manganese in steel (hypothetical data).
We can fit eqn. (3) by defining an X matrix, the columns of which contain the coefficients of the model parameters for each experiment. 1
3
1
4
1 4
X-
1 4 1 5 1 5 1 1 1 1
6 6 6 7
9 16 16 16 25 25 36 36 36 49
b = (X’X)-’(X’y)
264
[
1460 32.155
8388 - 13.183 5 SO62 1.2826 - 0.54348
1.2826
- 0.53438
0.054348
1
( x i y ) = ]::4:
=
[
22534
- 158.41 103.40
- 10.326
]
The resulting equation is y1=
- 158.41 + 103.40x1-
10.326~:
(4)
This model is plotted as the solid line in Fig. 6. The agreement between the model Reference p. 269
260
I
I
8
Per cent
mongonese
Fig. 6. Experimental results and fitted model for data taken from the system described in Fig. 4.
and the data is good; any deviations seem to be accounted for by the uncertainties that exist in the measurements themselves and not by any lack of fit between the model and the data. Presumably, the results shown in Fig. 6 came from the same system represented by Figs. 4 and 5 . Although the model represented by eqn. (4) is adequate over the domain of acquired data shown in Fig. 6, serious discrepancies between what this model predicts and the actual behavior of the system (Fig. 5 ) would exist both at lower and at higher percent manganese. This discrepancy, or lack of fit, is a source of uncertainty when trying to use an inadequate model to describe the true behavior of systems. In general, mechanistic models should be used if systems are adequately understood and there is some guarantee that the systems will not deviate greatly from their expected behavior. Otherwise, mechanistic models might be seriously misleading for predictive purposes; unbiased empirical models might prove to be better choices. Empirical models must be used when systems are not well understood and mechanistic models are not possible. Full second-order polynomial models [e.g. eqn. (4)] are very useful empirical models. We will use these models in much of the next few chapters because of their versatility in describing a wide variety of naturally occurring response surfaces over limited domains of the factors.
261
f-P
n-f
Fig. 7. Sums of squares and degrees of freedom tree for least-squares treatment of linear models containing a b, term.
The statistical basis of judging the adequacy of models has been given in Chap. 5 and is summarized here as the sum of squares and degrees of freedom tree shown in Fig. 7. The coefficient of multiple determination (R2)and the F-test for the significance of regression, FREG, offer useful means of evaluating how good the model is. The F-test for lack of fit, FLOF, is useful as a means of deciding if a better model can reasonably be expected to be found. We can calculate the values of the coefficient of multiple determination (i.e. the sum of squares due to regression divided by the sum of squares corrected for the mean), the F-ratio for the significance of regression, and the F-ratio for the lack of fit for the results of the previous example and can calculate at what levels the F-values are significant.
where n is the number of experiments in the set, f is the number of distinctly different factor combinations, and p is the number of parameters in the model. The F-ratio for regression is significant at the 99.99965%level of confidence. The F-ratio for lack of fit is significant at the 41% level of confidence, i.e. there is no reason to be concerned about the lack of fit, the model is probably adequate.
3. Two-factor response surfaces and models Many systems exhibit responses that are functions of not one but two factors. Examples of such systems are absorbance as a function of both the determinand Reference p. 269
262
and an interferent, tensile strength of steel as a function of both manganese and cobalt concentrations, and liquid chromatographic retention time as a function of both pH and ion-pairing reagent concentration. In this section, we show the application of empirical models to the two-factor case. The concepts developed here can be easily expanded to cover the multifactor case. 3.I First -order models
Let us consider again the photometric measurement technique described in Sect. 2.1, but rewrite Beer’s Law in the form A j = cjbC,
where A, is the absorbance (= -lo&Z/Zo]) of compound j at a given wavelength A, E, is the molar absorptivity of j at A, and is the molar concentration of j . If an interfering substance, k, also absorbs radiation at this wavelength, then the absorbance caused by this interfering substance, A,, may be written A , = €&bC&
(6)
where € & is the molar absorptivity of the interferent and c&is the concentration of the interferent. It is a characteristic of photometric absorption that the total absorbance, A,, is the sum of the individual absorbances of all absorbing species
+
A , = AJ + A , = cJbC, €kbCk
(7)
This model may be written in the form y11= blXl1
+ b 2 X 2 I + el,
(8)
where yl, = A , , b, = cJb,b2 = ckb, x l l = C, and x2, = c&. A response surface for this model is shown in Fig. 8. Suppose we have obtained experimental data at xll = 1, x2, = 1 and xI2= 2, x22= 2.
If we try to fit the model represented by eqn. (8) to this data, then
The determinant of (X’X)is equal to zero and the (X’X)-’ matrix is undefined. Therefore, the model cannot be fit to this data. An explanation for this is that the data points (and the origin) all lie in a straight line (i.e. they are co-linear). But the model is that of a plane. Thus, an infinite number of planes could pass equally well through the line of data and no unique solution exists.
263
Fig. 8. Two-factor response surface for the model y , = ;xl
+ fx2. See eqn. (8).
Let us try a different set of data obtained at (1, 2) and (2, 1).
The two data points plus the origin are not co-linear and the fitted plane is uniquely defined. The equation of the plane is .Y1=
fx, + 4 x 2
Note that, in Fig. 8, the slope of the response surface with respect to the factor x, is always equal to b,, independent of the values of x, and x,. Similarly, the slope of the response surface with respect to the factor x2 is always equal to b,. This is a proper interpretation of eqn. (8): b, and b, are the partial derivatives of the response y1 with respect to x, and x,, respectively.
Reference p. 269
264
A more generally useful form of the two-factor first-order polynomial model includes a bo term. = bo
'1
+ bix, + b 2 ~ 2
(11)
With this model, the response is not constrained to go through the origin. Addition of the b, term may be thought of as giving the response surface the degree of freedom required to move up or down in a "vertical" direction. 3.2 First-order models with interaction
In many systems, the effect of one factor will depend upon the value (or leuel) of another factor. This phenomenon is called factor interaction and may be written mathematically as
-a'' - - b, + b12x2
ax' That is, when the value of x2 is equal to zero, the slope of the response surface with respect to x1 is equal to b , , as before. However, when x 2 is not equal to zero, the slope of the response surface with respect to x 1 depends upon the value of x 2 . Substituting this value (i3y,/ax1) for the effects of x , in eqn. (11) gives Y1i
= b , + ( 6 , + b12~2i)xli+ b 2 ~ 2+ i e1i
or
+
+
(1 3)
+
y l i = b, + b,xli 6 1 2 ~ 1 j ~ 2 jel, (14) This is the equation of a full two-factor first-order model with interaction. It is interesting to note that eqn. (14) can be rearranged to give
-a'' - - b2 + b,,x, 3x2
as expected; that is, the effect of the factor x 2 also depends upon the value of the factor xl. As an example, suppose we want to investigate the effect of two variables, chart speed x1 and recorder sensitivity x 2 , on the precision of measuring chromatographic peaks. Assume we have gathered the following set of factor levels and corresponding responses.
The factor levels have been "coded" so that a real chart speed of, say, 1 cm min-' is
265
given a coded value of - 1 and a chart speed of 5 cm min-’ is given a coded value of 1. Similarly, coded recorder sensitivities of - 1 and + 1 might correspond to 1 mV full scale and 5 mV full scale, respectively. The responses (uncoded) might represent the standard deviation of peak heights from 10 repetitive injections. Let us fit the model given by eqn. (14) to this data. + I -1 -1 + 1 4 0 0 0 0 4 0 0 0 0 4 0 0 0 0 4
+
(X’y)=
“a 4
b = (X’X)-’\R.’y) =
[/]
1
0 0 0
1/4 0 0
0 1/4 0
O 0 I 1/41
A graph of the response surface for the fitted model y, = 4 + lx, lx, Ix,x,
+
+
(17)
is given in Fig. 9 where the concept of interaction is clearly evident. An interpreta-
Fig. 9. Two-factor response surface for the model yI = 4 + lx, + 1x2+ Ixlxz. See eqn. (14).
Reference p. 269
266
tion of the figure is that, at low recorder sensitivity (- 1 level of x , ) , the chart speed does not have much of an effect on precision.
+ 1x1 + 1x2 + = 4 + 1(-1) + h2 + 1(
y’= 4
lX,X2
- 1 ) ~ 2
=3
+
However, at high recorder sensitivity ( 1level of x ~ ) high , chart speed ( xz) gives worse precision than does low chart speed (- 1 level of x2). y, = 4
(18)
+ 1 level of
+ 1x1 + 1x2 + l X l X 2
=4+1(+1)+ 1x2+l(+l)xz =5
+ 2x2
3.3 Full second-order polynomial models
The concept of interaction can also be applied to individual factors. For example, in Fig. 6, the effect of percent manganese (factor x,) depends upon the value of that same factor. That is
a’ - - 6 , + b , , x , ax,
If a second factor were involved, it might also be true that the effect of this second
Fig. 10. Two-factor response surface showing a possible relationship between square wave polarographic sensitivity, voltage increment, and drop time. Scaled equation: y , = 0+0.8x, +0.8x2 - 0 . 0 8 ~ :- 0 . 0 8 ~ : .
267
1 0
Fig. 11. Two-factor response surface showing a possible relationship between the relative effect of bromine interference in the determination of chlorine, pH. and temperature. Scaled equation: yI = 40 . 8 -~0 .~8 ~ 2+O.ORx? +O.ORx;.
Fig. 12. Two-factor response surface showing a possible relationship between sample throughput efficiency of a laboratory, pay rate, and number of personnel. Scaled equation: .vl = 1.75 + 0.25~~ +
0.50~2 -0.04~:-0.03xf+0.02~1~,.
Reference p. 269
268
Fig. 13. Two-factor response surface showing a possible relationship between turbidity, temperature, and pH. Scaled equation: y l = - 9 - 0 . 4 x 1 + 7 x 2 - 0 . 0 4 x ~ - I x ~ +0.2xIx2.
factor depends upon its value
Substitution of these self-interaction effects into the model expressed by eqn. (14) gives yli
= bo
+ (b1 + b l l x , i ) x , i + (b2 + b 2 2 ~ 2 i I x 2 + i b12~1ix2i+ e1i
(22)
Fig. 14. Two-factor response surface showing a possible relationship between analytical recovery and the concentrations of two dilferent reagents. Scaled equation: v, = 2+0.8x1 - 0 . 8 ~-~0.08~:+ 0.08~:.
269
Fig. 15. Two-factor response surface showing a possible relationship between reaction rate and two different reactants. Scaled equation: y , = 0.04x1x,.
or
+
y l i = 6, b,x,, + b2x2,+ b , , x ~+ , bZ2xii+ b 1 2 ~ l i x+2 eli i (23) Equation (23) represents a full two-factor second-order polynomial model. Such equations are exceptionally versatile for use as empirical models in many systems over a limited domain of the factors [I]. Figures 10-15 show examples of response surfaces that can be fit by the model of eqn. (23). These figures might represent, for example, the following analytical chemical response surfaces: Fig. 10, square wave polarographic sensitivity as a function of voltage increment and drop time; Fig. 11, relative effect of bromine interference in the determination of chlorine as a function of pH and temperature; Fig. 12, sample throughput efficiency of a laboratory as a function of pay rate and number of personnel; Fig. 13, turbidity as a function of temperature and pH; Fig. 14, analytical recovery as a function of the concentrations of two different analytical reagents; and Fig. 15, reaction rate as a function of two different reactants.
Reference 1 G.E.P. Box and K.B. Wilson. On the experimental attainment of optimum conditions, J. R. Stat. Soc. Ser. B, 13 (1951) 1.
Recommended reading T.B. Barker, Quality by Experimental Design, Dekker, New York, 1985. L. von Bertalanffy, General System Theory. Foundations, Development, Applications, Braziller, New York, 1968. G.E.P. Box, W.G. Hunter and J.S. Hunter, Statistics for Experimenters. An Introduction to Design, Data Analysis, and Model Building, Wiley, New York, 1978.
270
W.G. Cochran and G.M. Cox, Experimental Designs, Wiley, New York, 1957. O.L. Davies (Ed.), Design and Analysis of Industrial Experiments, Hafner, New York, 2nd edn., 1956. S.N. Deming and S.L.Morgan, Experimental Design: A Chemometric Approach, Elsevier, Amsterdam, 1987.
A.J. Duncan, Quality Control and Industrial Statistics, Irwin, Homewood, IL, 1959. R.A. Fisher. The Design of Experiments, Hafner, New York, 1971. P.M.W. John, Statistical Design and the Analysis of Experiments, MacMillan, New York, 1971. J. Mandel. The Statistical Analysis of Experimental Data, Wiley. New York, 1964. W. Mendenhall, Introduction to Linear Models and the Design and Analysis of Experiments, Duxbury, Belmont, CA, 1968. M.G. Natrella, Experimental Statistics, National Bureau of Standards Handbook 91, U.S. Government Printing Office, Washington, DC,1963. G.W. Snedecor and W.G. Cochran, Statistical Methods, Iowa State University Press, Ames, IA, 6th edn., 1967.
G.M. Weinberg, An Introduction to General Systems Thinking, Wiley, New York, 1975. E.B. Wilson, Jr., An Introduction to Scientific Research, McGraw-Hill, New York, 1952. W.J. Youden, Statistical Methods for Chemists, Wiley, New York, 1951.
Chapter 16 Response Surfaces and Models 1. Introduction
The concepts developed in this chapter and the next three chapters are general and can be applied not only to analytical chemical systems (e.g. electrochemistry, spectroscopy, chromatography) but also to the systems which are being measured by the analytical chemist (e.g. pharmaceutical manufacturing, polymer production, water quality). While these concepts offer valuable insights for improving and understanding the responses from systems, it is well to keep in mind that we might never be able to discover the true behavior of a system. Uncertainties associated with the experimental measurements and with the mathematical description of the system will almost always exist. Thus, we will usually have a more or less “fuzzy” picture of how the system truly behaves. It follows that if we are to have a sufficiently clear view of how the system behaves, we must maximize our fundamental understanding of the system and minimize the uncertainties associated with our experimental measurements. These dual goals are achieved by the simultaneous application of precise measurement techniques, good chemical models, and proper experimental design.
2. Response surfaces A response surface is the graph of a system response plotted as a function of one or more of the system factors (i.e. “independent variables”, see Chap. 13). Response surfaces offer the chemometrician a convenient means of visualizing how various factors affect his or her measurement system. Consider the familiar analytical example shown in Fig. 1, a photometric method for the classical determination of manganese as permanganate. The amount of light transmitted through the system depends upon the concentration of MnO; in the sample. 2.1 True response surfaces
According to Beer’s law (a mathematical model), if the intensity of the light source is absolutely stable, if the light is monochromatic, if there is no stray light in the system, if the MnOi does not undergo chemical or physico-chemical changes as a function of concentration, and if there are no other absorbing species present in solution, then the intensity of transmitted light is given by
z = lox Reference p. 269
256
Light source
Photometric system
tntensi ty
Fig. 1. Systems theory diagram of a photometric method for the classical determination of manganese as permanganate.
and the transmittance, T, is given by
where I is the intensity or power of the transmitted light, I, is the intensity or
Fig. 2. Scaled photometric intensity as a function of the concentration of permanganate in a sample.
251
power of the light source, c is the molar absorptivity in L mol-’ cm-’,b is the path length through the sample in cm, and C is the concentration of MnOc in the sample in mol I-’. In this system, the intensity of transmitted light, I, is the system response and the concentration of MnO; in the sample, C, is the system factor. The true response surface (i.e. the graph of the system response plotted as a function of the single system factor) is shown in Fig. 2. Note that, for this true response surface, the physical understanding is exact and there is no experimental uncertainty. 2.2 Measured response surfaces Although many system diagrams are drawn in a manner similar to Fig. 1, such diagrams are incomplete because they do not explicitly take into account a separate but necessary measurement system. This is a subtle but important point. In the photometric determination of MnO; , for example, transmitted intensity is the important physical property of interest, but we seldom measure it directly. Instead, with modern photometric instruments, transmitted intensity is converted to an electrical current, then to a voltage, and finally to an analog or digital meter readout from which we visually “observe” the value of the transmitted intensity. A correct view of the double-beam photometric measurement system is given in Fig. 3. It is clear that the measurement system itself can introduce a set of uncertainties and biases quite separate from those introduced by the photometric system. For example, if electrical line voltage fluctuations affect the measurement system in such a way that the readout fluctuates, the originally well-defined response surface of Fig. 2 becomes uncertain again. In general, measured response surfaces are ‘‘fuzzy” because of uncertainties in the primary system itself and because of uncertainties introduced by the associated measurement system. It is usually desirable that this latter source of uncertainty should not
I T
0 f
Light source
Photometric
Measuremeo t
Readout
Fig. 3. Systems theory diagram showing the relationship between the photometric system and the measurement system for a double-beam photometer.
Reference p. 269
258
be greater than the former; hence the desire for precise (i.e. low uncertainty) analytical chemical methods. 2.3 Estimated response surfaces and models
Consider the set of experimental points shown in Fig. 6. It is probably unrealistic to try to write a mechanistic model (see Chap. 12) to explain these data. Our knowledge of metallurgy would be quickly exceeded. It is possible, however, to propose an empirical model that might provide a good description of the data. One possible empirical model might be the equation of a parabola
+
y l i = b, + b,x,, + b,,xfi eli (3) where y I i represents the system response (tensile strength) for experiment i, x l i represents the system factor (percent manganese in this case), b, is an offset parameter, b, is a first-order parameter, b,, is a second-order parameter, and el, is a residual or deviation between what is observed and what is predicted by the model. This is a single-factor model of a single-factor response surface. Matrix least-squares techniques (see Chap. 13) can be used to fit this model to the data shown in Fig. 6. As an example, suppose the data in Fig. 6 have the following numerical values (D is an “experimental design matrix” that has a number of rows equal to the number of experiments and a number of columns equal to the number of factors; each element of the D matrix thus represents the value of a particular factor in a particular experiment). 3 58 4 87 93 4 4 91 105 5 D= Y= 98 5 6 88 6 90 6 89 761
-.
Iron or Mongonese
system
meosuremen t
Tensile -strength meosuremen t
% Mn
Tensile strength
Fig. 4. Systems theory diagram showing the relationships among a steel-making system and two measurement systems.
259
Fig. 5. Tensile strength as a function of the percent manganese in steel (hypothetical data).
We can fit eqn. (3) by defining an X matrix, the columns of which contain the coefficients of the model parameters for each experiment. 1
3
1
4
1 4
X-
1 4 1 5 1 5 1 1 1 1
6 6 6 7
9 16 16 16 25 25 36 36 36 49
b = (X’X)-’(X’y)
264
[
1460 32.155
8388 - 13.183 5 SO62 1.2826 - 0.54348
1.2826
- 0.53438
0.054348
1
( x i y ) = ]::4:
=
[
22534
- 158.41 103.40
- 10.326
]
The resulting equation is y1=
- 158.41 + 103.40x1-
10.326~:
(4)
This model is plotted as the solid line in Fig. 6. The agreement between the model Reference p. 269
260
I
I
8
Per cent
mongonese
Fig. 6. Experimental results and fitted model for data taken from the system described in Fig. 4.
and the data is good; any deviations seem to be accounted for by the uncertainties that exist in the measurements themselves and not by any lack of fit between the model and the data. Presumably, the results shown in Fig. 6 came from the same system represented by Figs. 4 and 5 . Although the model represented by eqn. (4) is adequate over the domain of acquired data shown in Fig. 6, serious discrepancies between what this model predicts and the actual behavior of the system (Fig. 5 ) would exist both at lower and at higher percent manganese. This discrepancy, or lack of fit, is a source of uncertainty when trying to use an inadequate model to describe the true behavior of systems. In general, mechanistic models should be used if systems are adequately understood and there is some guarantee that the systems will not deviate greatly from their expected behavior. Otherwise, mechanistic models might be seriously misleading for predictive purposes; unbiased empirical models might prove to be better choices. Empirical models must be used when systems are not well understood and mechanistic models are not possible. Full second-order polynomial models [e.g. eqn. (4)] are very useful empirical models. We will use these models in much of the next few chapters because of their versatility in describing a wide variety of naturally occurring response surfaces over limited domains of the factors.
261
f-P
n-f
Fig. 7. Sums of squares and degrees of freedom tree for least-squares treatment of linear models containing a b, term.
The statistical basis of judging the adequacy of models has been given in Chap. 5 and is summarized here as the sum of squares and degrees of freedom tree shown in Fig. 7. The coefficient of multiple determination (R2)and the F-test for the significance of regression, FREG, offer useful means of evaluating how good the model is. The F-test for lack of fit, FLOF, is useful as a means of deciding if a better model can reasonably be expected to be found. We can calculate the values of the coefficient of multiple determination (i.e. the sum of squares due to regression divided by the sum of squares corrected for the mean), the F-ratio for the significance of regression, and the F-ratio for the lack of fit for the results of the previous example and can calculate at what levels the F-values are significant.
where n is the number of experiments in the set, f is the number of distinctly different factor combinations, and p is the number of parameters in the model. The F-ratio for regression is significant at the 99.99965%level of confidence. The F-ratio for lack of fit is significant at the 41% level of confidence, i.e. there is no reason to be concerned about the lack of fit, the model is probably adequate.
3. Two-factor response surfaces and models Many systems exhibit responses that are functions of not one but two factors. Examples of such systems are absorbance as a function of both the determinand Reference p. 269
262
and an interferent, tensile strength of steel as a function of both manganese and cobalt concentrations, and liquid chromatographic retention time as a function of both pH and ion-pairing reagent concentration. In this section, we show the application of empirical models to the two-factor case. The concepts developed here can be easily expanded to cover the multifactor case. 3.I First -order models
Let us consider again the photometric measurement technique described in Sect. 2.1, but rewrite Beer’s Law in the form A j = cjbC,
where A, is the absorbance (= -lo&Z/Zo]) of compound j at a given wavelength A, E, is the molar absorptivity of j at A, and is the molar concentration of j . If an interfering substance, k, also absorbs radiation at this wavelength, then the absorbance caused by this interfering substance, A,, may be written A , = €&bC&
(6)
where € & is the molar absorptivity of the interferent and c&is the concentration of the interferent. It is a characteristic of photometric absorption that the total absorbance, A,, is the sum of the individual absorbances of all absorbing species
+
A , = AJ + A , = cJbC, €kbCk
(7)
This model may be written in the form y11= blXl1
+ b 2 X 2 I + el,
(8)
where yl, = A , , b, = cJb,b2 = ckb, x l l = C, and x2, = c&. A response surface for this model is shown in Fig. 8. Suppose we have obtained experimental data at xll = 1, x2, = 1 and xI2= 2, x22= 2.
If we try to fit the model represented by eqn. (8) to this data, then
The determinant of (X’X)is equal to zero and the (X’X)-’ matrix is undefined. Therefore, the model cannot be fit to this data. An explanation for this is that the data points (and the origin) all lie in a straight line (i.e. they are co-linear). But the model is that of a plane. Thus, an infinite number of planes could pass equally well through the line of data and no unique solution exists.
263
Fig. 8. Two-factor response surface for the model y , = ;xl
+ fx2. See eqn. (8).
Let us try a different set of data obtained at (1, 2) and (2, 1).
The two data points plus the origin are not co-linear and the fitted plane is uniquely defined. The equation of the plane is .Y1=
fx, + 4 x 2
Note that, in Fig. 8, the slope of the response surface with respect to the factor x, is always equal to b,, independent of the values of x, and x,. Similarly, the slope of the response surface with respect to the factor x2 is always equal to b,. This is a proper interpretation of eqn. (8): b, and b, are the partial derivatives of the response y1 with respect to x, and x,, respectively.
Reference p. 269
264
A more generally useful form of the two-factor first-order polynomial model includes a bo term. = bo
'1
+ bix, + b 2 ~ 2
(11)
With this model, the response is not constrained to go through the origin. Addition of the b, term may be thought of as giving the response surface the degree of freedom required to move up or down in a "vertical" direction. 3.2 First-order models with interaction
In many systems, the effect of one factor will depend upon the value (or leuel) of another factor. This phenomenon is called factor interaction and may be written mathematically as
-a'' - - b, + b12x2
ax' That is, when the value of x2 is equal to zero, the slope of the response surface with respect to x1 is equal to b , , as before. However, when x 2 is not equal to zero, the slope of the response surface with respect to x 1 depends upon the value of x 2 . Substituting this value (i3y,/ax1) for the effects of x , in eqn. (11) gives Y1i
= b , + ( 6 , + b12~2i)xli+ b 2 ~ 2+ i e1i
or
+
+
(1 3)
+
y l i = b, + b,xli 6 1 2 ~ 1 j ~ 2 jel, (14) This is the equation of a full two-factor first-order model with interaction. It is interesting to note that eqn. (14) can be rearranged to give
-a'' - - b2 + b,,x, 3x2
as expected; that is, the effect of the factor x 2 also depends upon the value of the factor xl. As an example, suppose we want to investigate the effect of two variables, chart speed x1 and recorder sensitivity x 2 , on the precision of measuring chromatographic peaks. Assume we have gathered the following set of factor levels and corresponding responses.
The factor levels have been "coded" so that a real chart speed of, say, 1 cm min-' is
265
given a coded value of - 1 and a chart speed of 5 cm min-’ is given a coded value of 1. Similarly, coded recorder sensitivities of - 1 and + 1 might correspond to 1 mV full scale and 5 mV full scale, respectively. The responses (uncoded) might represent the standard deviation of peak heights from 10 repetitive injections. Let us fit the model given by eqn. (14) to this data. + I -1 -1 + 1 4 0 0 0 0 4 0 0 0 0 4 0 0 0 0 4
+
(X’y)=
“a 4
b = (X’X)-’\R.’y) =
[/]
1
0 0 0
1/4 0 0
0 1/4 0
O 0 I 1/41
A graph of the response surface for the fitted model y, = 4 + lx, lx, Ix,x,
+
+
(17)
is given in Fig. 9 where the concept of interaction is clearly evident. An interpreta-
Fig. 9. Two-factor response surface for the model yI = 4 + lx, + 1x2+ Ixlxz. See eqn. (14).
Reference p. 269
266
tion of the figure is that, at low recorder sensitivity (- 1 level of x , ) , the chart speed does not have much of an effect on precision.
+ 1x1 + 1x2 + = 4 + 1(-1) + h2 + 1(
y’= 4
lX,X2
- 1 ) ~ 2
=3
+
However, at high recorder sensitivity ( 1level of x ~ ) high , chart speed ( xz) gives worse precision than does low chart speed (- 1 level of x2). y, = 4
(18)
+ 1 level of
+ 1x1 + 1x2 + l X l X 2
=4+1(+1)+ 1x2+l(+l)xz =5
+ 2x2
3.3 Full second-order polynomial models
The concept of interaction can also be applied to individual factors. For example, in Fig. 6, the effect of percent manganese (factor x,) depends upon the value of that same factor. That is
a’ - - 6 , + b , , x , ax,
If a second factor were involved, it might also be true that the effect of this second
Fig. 10. Two-factor response surface showing a possible relationship between square wave polarographic sensitivity, voltage increment, and drop time. Scaled equation: y , = 0+0.8x, +0.8x2 - 0 . 0 8 ~ :- 0 . 0 8 ~ : .
267
1 0
Fig. 11. Two-factor response surface showing a possible relationship between the relative effect of bromine interference in the determination of chlorine, pH. and temperature. Scaled equation: yI = 40 . 8 -~0 .~8 ~ 2+O.ORx? +O.ORx;.
Fig. 12. Two-factor response surface showing a possible relationship between sample throughput efficiency of a laboratory, pay rate, and number of personnel. Scaled equation: .vl = 1.75 + 0.25~~ +
0.50~2 -0.04~:-0.03xf+0.02~1~,.
Reference p. 269
268
Fig. 13. Two-factor response surface showing a possible relationship between turbidity, temperature, and pH. Scaled equation: y l = - 9 - 0 . 4 x 1 + 7 x 2 - 0 . 0 4 x ~ - I x ~ +0.2xIx2.
factor depends upon its value
Substitution of these self-interaction effects into the model expressed by eqn. (14) gives yli
= bo
+ (b1 + b l l x , i ) x , i + (b2 + b 2 2 ~ 2 i I x 2 + i b12~1ix2i+ e1i
(22)
Fig. 14. Two-factor response surface showing a possible relationship between analytical recovery and the concentrations of two dilferent reagents. Scaled equation: v, = 2+0.8x1 - 0 . 8 ~-~0.08~:+ 0.08~:.
269
Fig. 15. Two-factor response surface showing a possible relationship between reaction rate and two different reactants. Scaled equation: y , = 0.04x1x,.
or
+
y l i = 6, b,x,, + b2x2,+ b , , x ~+ , bZ2xii+ b 1 2 ~ l i x+2 eli i (23) Equation (23) represents a full two-factor second-order polynomial model. Such equations are exceptionally versatile for use as empirical models in many systems over a limited domain of the factors [I]. Figures 10-15 show examples of response surfaces that can be fit by the model of eqn. (23). These figures might represent, for example, the following analytical chemical response surfaces: Fig. 10, square wave polarographic sensitivity as a function of voltage increment and drop time; Fig. 11, relative effect of bromine interference in the determination of chlorine as a function of pH and temperature; Fig. 12, sample throughput efficiency of a laboratory as a function of pay rate and number of personnel; Fig. 13, turbidity as a function of temperature and pH; Fig. 14, analytical recovery as a function of the concentrations of two different analytical reagents; and Fig. 15, reaction rate as a function of two different reactants.
Reference 1 G.E.P. Box and K.B. Wilson. On the experimental attainment of optimum conditions, J. R. Stat. Soc. Ser. B, 13 (1951) 1.
Recommended reading T.B. Barker, Quality by Experimental Design, Dekker, New York, 1985. L. von Bertalanffy, General System Theory. Foundations, Development, Applications, Braziller, New York, 1968. G.E.P. Box, W.G. Hunter and J.S. Hunter, Statistics for Experimenters. An Introduction to Design, Data Analysis, and Model Building, Wiley, New York, 1978.
270
W.G. Cochran and G.M. Cox, Experimental Designs, Wiley, New York, 1957. O.L. Davies (Ed.), Design and Analysis of Industrial Experiments, Hafner, New York, 2nd edn., 1956. S.N. Deming and S.L.Morgan, Experimental Design: A Chemometric Approach, Elsevier, Amsterdam, 1987.
A.J. Duncan, Quality Control and Industrial Statistics, Irwin, Homewood, IL, 1959. R.A. Fisher. The Design of Experiments, Hafner, New York, 1971. P.M.W. John, Statistical Design and the Analysis of Experiments, MacMillan, New York, 1971. J. Mandel. The Statistical Analysis of Experimental Data, Wiley. New York, 1964. W. Mendenhall, Introduction to Linear Models and the Design and Analysis of Experiments, Duxbury, Belmont, CA, 1968. M.G. Natrella, Experimental Statistics, National Bureau of Standards Handbook 91, U.S. Government Printing Office, Washington, DC,1963. G.W. Snedecor and W.G. Cochran, Statistical Methods, Iowa State University Press, Ames, IA, 6th edn., 1967.
G.M. Weinberg, An Introduction to General Systems Thinking, Wiley, New York, 1975. E.B. Wilson, Jr., An Introduction to Scientific Research, McGraw-Hill, New York, 1952. W.J. Youden, Statistical Methods for Chemists, Wiley, New York, 1951.