Characterizing the Uncertainty in Near-Infrared Spectroscopic

Anal. Chem. 1998, 70, 2972-2982

Characterizing the Uncertainty in Near-Infrared Spectroscopic Prediction of Mixed-Oxygenate Concentrations in Gasoline: Sample-Specific Prediction Intervals Nicolaas (Klaas) M. Faber†

Center for Process Analytical Chemistry, University of Washington, Seattle, Washington 98195 David L. Duewer,* Steven J. Choquette, Terry L. Green, and Stephen N. Chesler‡

Analytical Chemistry Division, Chemical Science and Technology Laboratory, National Institute of Standards and Technology, Gaithersburg, Maryland 20899

Oxygenates are added to gasoline to reduce exhaust emission levels of carbon monoxide and to boost octane. The U.S. National Institute of Standards and Technology (NIST) provides 12 Standard Reference Materials (SRMs) for single oxygenates in reference gasoline. A previous study demonstrated the feasibility of nondestructively quantifying oxygenate concentration in SRM gasoline ampules using near-infrared (near-IR) spectroscopy combined with multivariate calibration techniques. A drawback of this approach has been that an average prediction uncertainty, rather than a sample-specific one, is obtained. Recent developments in multivariate calibration theory for prediction error variance cure this problem. This report characterizes the significant sources of uncertainties in multivariate calibration using principal component regression and partial least-squares, validating near-IR and other multivariate spectroscopic techniques for use in assigning certified values (expected value with specified uncertainty) to selected materials. This report interprets prediction results in terms of multivariate analytical figures of merit, enabling the visualization of complex multivariate models as univariate graphs. Recent U.S. clean air legislation mandates that gasoline sold during the winter months in many metropolitan areas must have an oxygen mass fraction of at least 2.7% O and that gasoline in some areas must be at least 2.0% O throughout the year.1 Ether oxygenates have been used as octane-number-boosting replacements for organometallic additives and may play a role in replacing aromatics in gasoline.2 Current analytical techniques for measur* Corresponding author. Tel.: 301-975-3935. Fax: 301-977-0587. Mail: Chemistry/222, Room A213. E-mail: [email protected]. † Current address: Netherlands Forensic Science Institute, Volmerlaan 17, 2288GD Rijswijk, The Netherlands. ‡ Current address: GEOMET Technologies, Inc., 20251 Century Blvd., Germantown, MD 20874-1192. (1) Kortum, D. J.; Miller, M. G. Prepr. Pap.sAm. Chem. Soc., Div. Fuel Chem. 1994, 39 (2), 267-272. (2) Piel, W. J. Prepr. Pap.sAm. Chem. Soc., Div. Fuel Chem. 1994, 39 (2), 273281.

2972 Analytical Chemistry, Vol. 70, No. 14, July 15, 1998

ing oxygenates in commercial gasoline include wet chemical methods, pyrolysis, atomic emission detection, gas chromatography (GC), mass spectrometry, nuclear magnetic resonance, and infrared spectrometry.3-5 The U.S. National Institute of Standards and Technology (NIST) provides 12 oxygenated gasoline Standard Reference Materials (SRMs) to help the petrochemical industry comply with legislative mandates. SRMs at nominal 2.0% O and 2.7% O are currently available for four different oxygenates: tert-amyl methyl ether, ethyl tert-butyl ether, methyl tert-butyl ether (MTBE), and ethanol (EtOH).6,7 These materials are gravimetrically prepared in large quantities, aliquots are sealed in 20-mL glass ampules, and the composition of subsamples of the ampule population is confirmed with GC techniques. Since the GC method requires opening each sampled ampule, the SRMs are batch-certified with a total %O uncertainty incorporating an estimate of batch heterogeneity. While the current SRMs all have an acceptable 0.04% O assigned batch uncertainty (95% confidence), a small number of outliers makes this uncertainty overly pessimistic for the large majority of individual oxygenate SRM ampules. Since oxygenates (3) Cooper, J. B.; Wise, K. L.; Welch, W. T.; Bledscoe, R. R.; Sumner, M. B. Appl. Spectrosc. 1996, 50, 917-921. (4) Manjarrez, A.; Capella, S.; Villaruel, E.; Pablo Garcia, F. Riv. Combust. 1992, 46, 337-340. (5) Hardman, J. S.; Hill, M. A.; Mills, G. A. Fuel 1993, 72, 1563-1566. (6) National Institute of Standards and Technology. Certificates of Analysis, Standard Reference Materials 2286sEthanol in Reference Gasoline, Nominal 2.0% Oxygen; 2287sEthanol in Reference Gasoline, Nominal 3.5% Oxygen; 2288stert-Amyl-methyl Ether in Reference Gasoline, Nominal 2.0% Oxygen; 2289stert-Amyl-methyl Ether in Reference Gasoline, Nominal 2.7% Oxygen; 2290sEthyl-tert-butyl Ether in Reference Gasoline, Nominal 2.0% Oxygen; 2291sEthyl-tert-butyl Ether in Reference Gasoline, Nominal 2.7% Oxygen; 2292sMethyl-tert-butyl Ether in Reference Gasoline, Nominal 2.0% Oxygen; 2293sMethyl-tert-butyl Ether in Reference Gasoline, Nominal 2.0% Oxygen; Standard Reference Materials Program, NIST: Gaithersburg, MD, 1995. (7) National Institute of Standards and Technology. Certificates of Analysis, Standard Reference Materials 2294sReformulated Gasoline, Nominal 11% MTBE, 35 mg/kg Sulfur; 2295sReformulated Gasoline, Nominal 15% MTBE, 300 mg/kg Sulfur; 2296sReformulated Gasoline, Nominal 13% ETBE, 35 mg/kg Sulfur; 2297sReformulated Gasoline, Nominal 10% Ethanol, 300 mg/ kg Sulfur; Standard Reference Materials Program, NIST: Gaithersburg, MD, 1998. S0003-2700(97)01270-5 CCC: $15.00

© 1998 American Chemical Society Published on Web 06/05/1998

are expensive and vast quantities of gasoline are produced, major cost savings can be obtained by making even a small reduction in oxygenate measurement uncertainty. Reducing the uncertainty in the values assigned to NIST oxygenate SRMs could contribute to improving the accuracy of oxygenate analysis. We have previously demonstrated the feasibility of nondestructively quantifying oxygenate concentration in ampules of SRM gasoline using near-infrared (near-IR) spectroscopy combined with multivariate calibration techniques.8 This technology has the potential for reducing the assigned batch homogeneity component through outlier screening or eliminating it altogether through individual ampule certification. Near-IR techniques can be precise and affordable, and can be employed remotely using fiber optics.9 However, the low-specificity near-IR signals typically require calibration sets of some 2-5 samples for each multivariate factor employed, a sometimes burdensome requirement.10 Since the oxygenate SRM ampules are produced in very large batches, moderately large calibration sets represent a relatively small resource commitment, particularly if models for different singleoxygenate SRMs use a single mixed-oxygenate calibration set. The predictive utility of a multivariate model must be validated by evaluating results for “test set” samples not included in the calibration set or through use of computationally intensive methods such as “cross-validation,” “jackknife,” or “bootstrap”.11-13 In this study, we use both a large test set and cross-validation, enabling separate evaluation of how well a given calibration model recognizes the samples from which it was constructed and how well it predicts similar but unknown samples. The commonly used terms for these performance summaries are root-mean-square error of calibration (RMSEC) and root-mean-square error of prediction (RMSEP). For cross-validation, where models based on different subsets of the calibration samples are evaluated with the remaining calibration samples, results from all the multiple models are combined and summarized as the root-mean-square error of cross-validation (RMSECV). All three of these summaries represent average model performance for a group of samples rather than providing sample-specific uncertainty estimates. There have been recent advances in the underlying principles of multivariate calibration methods such as partial least squares (PLS)14 and uncertainty propagation through complex systems.15 Sample-specific uncertainties for PLS predictions using known uncertainties in analyte concentrations can now be estimated.16-19 These concepts have been extended to principal component (8) Choquette, S. J.; Chesler, S. N.; Duewer, D. L.; Wang, S.; O’Haver, T. C. Anal. Chem. 1996, 68, 3525-3533. (9) Litani-Barzilai, I.; Sela, I.; Bulatov, V.; Zilberman, I.; Schechter, I. Anal. Chim. Acta 1997, 339, 193-199. (10) Skloss, T. W.; Kim, A. J.; Haw, J. F. Anal. Chem. 1994, 66, 536-542. (11) Martens, H.; Næs, T. Multivariate Calibration; Wiley: New York, 1989. (12) Hardy, A. J.; MacLaurin, P.; Haswell, S. J.; de Jong, S.; Vandeginste, B. G. M. Chemom. Intell. Lab. Syst. 1996, 34, 117-129. (13) Wehrens, R.; Van der Linden, W. E. J. Chemom. 1997, 11, 157-171. (14) Helland, I. S. Commun. Stat.sElem. Simul. Comput. 1988, 17, 581-607. (15) Magnus, J. R.; Neudecker, H. Matrix Differential Calculus with Applications in Statistics and Econometrics; Wiley: New York, 1988. (16) Denham, M. C. Ph.D. Dissertation, University of Liverpool, U.K., Liverpool, 1991. (17) Denham, M. C. J. Chemom. 1997, 11, 39-52. (18) Phatak, A. Ph.D. Dissertation, University of Waterloo, Waterloo, Ontario, Canada. 1993. (19) Phatak, A.; Reilly, P. M.; Penlidis, A. Anal. Chim. Acta 1993, 277, 495501.

regression (PCR) and to include spectral measurement uncertainties.20,21 The goals of this paper are to (1) establish which sources of analytical uncertainty significantly contribute to multivariate calibration prediction uncertainty, (2) construct sample-specific prediction intervals for oxygenate systems of interest, and (3) illustrate how the recent theoretical advances in multivariate calibration technology can be made practical. MTBE and EtOH appear to be the oxygenates of greatest current economic importance,2 with H2O as an unavoidable major contaminant of EtOH. We therefore use a triple-oxygenate system in MTBE, EtOH, and H2O as our test case. MATHEMATICAL BACKGROUND Model Assumptions. A linear relationship is assumed to hold between analyte concentrations (expressed in %O) and near-IR spectral data. Incomplete knowledge about possible interferents dictates the use of the inverse calibration model, with the analyte concentrations treated as the dependent or “response” variables and the spectral data as the independent or “predictor” variables.11 Spectral measurements are performed for n calibration samples at p wavelengths. In matrix-vector notation, the model for the calibration samples is

c ) Rb + e

(1)

where c (n × 1) is the vector of true analyte concentrations, R (n × p) is the matrix of true spectral responses, b (p × 1) is the unknown regression vector, and e (n × 1) is a vector of residuals. The true analyte concentrations and spectral responses are, in principle, unknowable because of uncertainties intrinsic to all measurement systems. We denote the total measurement error for the analyte concentrations as ∆c and that for the spectral responses as ∆R. We denote the knowable measured analyte concentrations as c˜ ) c + ∆c and the spectral responses as R ˜ ) R + ∆R, where the tilde indicates measured rather than true quantities. Equation 1 shows how the true but unknowable quantities are related. Taking the unavoidable measurement errors into account,

c ˜ ) c + ∆c ) Rb + e + ∆c )R ˜ b - ∆Rb + e + ∆c

(2)

distinguishes the relevant sources of uncertainty in the calibration data and provides a convenient starting point for further analysis. Estimation of the Regression Vector. Ordinary least squares (OLS) estimates the regression vector of eq 2 as

bˆ OLS ) (R ˜ TR ˜ )-1R ˜ Tc ˜

(3)

where the “hat” indicates estimated quantities, superscript -1 signifies matrix inversion, and superscript T signifies matrix or vector transposition. The form of eq 3 is well known, but the (20) Faber, N. M.; Kowalski, B. R. J. Chemom. 1997, 11, 181-238. (21) Faber, N. M.; Kowalski, B. R. Chemom. Intell. Lab. Syst. 1996, 34, 283292.

Analytical Chemistry, Vol. 70, No. 14, July 15, 1998

2973

tilde over the response and predictor variables is seldom explicit. Making this distinction is essential to the current work, since measurement errors lead to additional terms in the expression for prediction error (PE) variance. In contrast, it is widely recognized that all results based on the analysis of measured data are estimated or predicted. Since making this status explicit in the notation is often awkward, we use the “hat” notation only when it is necessary for clarity. Equation 3 can be solved with OLS only if the (p × p) crossproduct matrix R ˜ TR ˜ can be inverted, i.e., is of full rank p. This is seldom the case in spectroscopic applications, since the number of wavelengths p typically exceeds the number of calibration samples n. The cross-product matrix will have rank of, at most, n and may, if some samples are related, be less than n. Further, degrees of freedom can be removed with data preprocessing: subtracting the mean over all calibration spectra from each of the spectra (“mean centering”) ensures that the cross-product matrix has a rank of, at most, n - 1. Throughout this paper, it will be assumed that the data (spectra and concentrations) are mean centered. While a number of solutions have been proposed to overcome inversion problems, only two are widely used in analytical spectroscopy: wavelength selection and latent variables. Selection of wavelengths prior to estimating the regression vector by OLS presents the obvious problem of selecting a relatively small subset out of a possibly very large set of wavelengths. (For the data studied in this paper, n and p are 40 and 391, respectively.) Finding the optimum subset amounts to solving a huge combinatorial problem and is generally intractable. Suboptimal solutions, which often perform satisfactorily, can be found using approaches developed for stepwise regression procedures.22 Construction of a small number, e.g., A, where A ) rank(R ˜ TR ˜ ), of score vectors (“latent variables”) as linear combinations of the p measured spectral variables is often an easier approach. Replacing the spectral matrix by the score matrix leads to an A-dimensional model that contains nearly all the spectral information in the measured data. This is the rationale behind methods such as principal component regression (PCR) and partial least squares (PLS). PCR and PLS use different criteria for the construction of the score matrix: PCR uses only the spectral information, whereas PLS uses both the spectral and the concentration information. PLS often requires fewer score vectors than does PCR to yield equivalently predictive models. Construction of the A score vectors from the p measured spectral variables tends to reduce the random measurement error components (“noise averaging”). Having a large number of spectral variables is therefore not necessarily a problem and may even be advantageous.11 Avoiding the problems inherent in wavelength selection and obtaining good noise averaging properties make the use of PCR and PLS attractive from a theoretical as well as a practical point of view. For both PCR and PLS methods, the score vectors are given by

ˆta ) R ˜w â

a ) 1, ..., A

(4)

where ˆta (n × 1) is the ath score vector and w ˆ a (p × 1) is the (22) Miller, A. J. Subset Selection in Regression; Chapman and Hall: London, 1990.


associated weight vector. The components of the weight vector determine how much the p measured spectral variables contribute to a specific score vector. These weights can enable chemical interpretations of the latent variables and have been used for combining PLS and wavelength selection.23,24 In analogy with eq 4, the full score matrix for an A-dimensional PCR or PLS model is

˜W Â Tˆ A ) R

(5)

where Tˆ A (n × A) and W ˆ A (p × A) are the score and weight matrices. If the score vectors in eq 4 are normalized to unit length, the regression vector obtained for A-dimensional PCR or PLS is

bˆ A ) (W ˆ AW ˆ AT)R ˜ Tc ˜

(6)

This simple form for the estimated regression vector is possible for PLS only if the weight vectors are defined with respect to the measured spectral matrix.25 In the common presentation of PLS, the weight vectors are defined with respect to a residual matrix.11 However, working with the measured spectral matrix is the quintessence of de Jong’s “SIMPLS” formalism.25 SIMPLS and PLS are identical if only one concentration is modeled at a time; since we model only one oxygenate concentration at a time in our application, we follow the SIMPLS formalism. Derivation of Approximate Expression for Prediction Error Variance. Given the estimated regression vector, the predicted analyte concentration in an unknown sample is

cû ) h˜c + ˜r uTbˆ A

(7)

where ˜hc is the average analyte concentration for the calibration set, ˜ru ) ru + ∆ru is the measured (mean centered) spectral vector for the unknown sample, and the subscript u indicates an unknown sample. In analogy to eq 1, the true analyte concentration is

cu ) jc + ruTb + eu

(8)

where eu is the unknown sample residual. The PE is defined as the difference between the prediction and the true value:

PEu ≡ cû - cu

(9)

Inserting eqs 7 and 8 into eq 9 and neglecting products of error terms gives

PEu ) ˜hc - jc + ˜r uTbˆ A - ruTb - eu ≈ ∆cj + ruT∆bA + ∆ruTbA - eu

(10)

It has been shown elsewhere that, to a very good approximation, (23) Frank, I. Chemom. Intell. Lab. Syst. 1987, 1, 233-242. (24) Lindgren, F.; Geladi, P.; Ra¨nnar, S.; Wold, S. J. Chemom. 1994, 8, 349363. (25) de Jong, S. Chemom. Intell. Lab. Syst. 1993, 18, 251-263.

the variance of this PE is given by20,21 2 V(PEu) ) E[(PEu)2] ≈ (1/n + hu)[σ2e + σ∆c + 2 2 ] + σe2u + ||bA||2 σ∆r (11) ||bA||2 σ∆R u

where E[‚‚‚] denotes taking expectation, hu is the so-called leverage of the unknown sample, ||‚‚‚|| symbolizes the Euclidean norm, and V and σ denote the variance and standard deviation (the square root of variance), respectively, of the associated quantity. Taking expectation of the squared PE effectively converts an unknowable statistical error into a determinable measure of prediction uncertainty. The leverage indicates how close the unknown sample is to the calibration samples in multivariate space. Mathematically, leverage is the inner product of ru with itself, taking the covariance structure of the spectral matrix into account so that a correct measure of distance is obtained:

hu ) ruT(WAWAT)ru

(12)

Interpretation of Prediction Error Variance Formula. The approximate PE variance thus has three origins: the uncertainty in determining the model center and regression vector (“intercept” and “slope”), the residual of the unknown sample, and the measurement error in the unknown sample spectrum. Their respective contributions are given by the three terms on the far right-hand side of eq 11. Mean-centering of the data leads to the well-known (1/n) term, which is constant. The prediction sample term is also constant. The leverage term is the only samplespecific contribution. Equation 11 can be further simplified by assuming that the variance of the residuals and the variance of the spectral measurement errors are independent of the sample so that the subscript u can be dropped for the last two terms. This chemically reasonable assumption will be made in the remaining part of this paper. There are no cross-terms because all uncertainties are assumed to be mutually independent. Equation 10 essentially constitutes a first-order Taylor expansion around the true values (“local linearization”). This leads to an expression for PE variance that contains unknowable, errorless quantities. In practice, these unknown quantities have to be replaced by measured, estimated, or predicted values. The practical evaluation of eq 11 is discussed below. The derivation of eq 11 assumes that the spectral measurement errors are independently and identically distributed. Generally, this is not a good assumption in near-IR applications; however, the results presented later clearly demonstrate that the effect of the spectral errors on PE variance is negligible for the current application. More elaborate expressions have been derived that accommodate heteroscedastic as well as correlated noise in the spectra.20 Maximum likelihood-based calibration procedures have been developed that tend to outperform PCR and PLS in cases where the heteroscedasticity of the spectral noise is severe.26 (26) Wentzell, P. D.; Andrews, D. T.; Kowalski, B. R. Anal. Chem. 1997, 69, 2299-2311.

Equation 11 describes only PE variance. Since PCR and PLS may provide biased predictions, the magnitude of such bias must be evaluated for each application.20 If prediction bias is significant, overly optimistic values for the prediction uncertainty interval may be obtained. However, in “underdetermined” (p > n) applications, sufficient factors can always be included to model all variation in the predictor variables “within the noise”. Bias distortions observed in “overdetermined” (p < n) systems17 are thus not necessarily applicable to “well-modeled” spectroscopic applications. While the correct use of eq 11 “depends critically on the number of factors selected”,17 in underdetermined applications prediction bias is not an intrinsic problem but rather an issue of model sufficiency. Equation 11 results from simplifying a much more complicated expression.20 It is a “zeroth-order” approximation because firstand higher-order terms are neglected. This approach has been labeled “naive” in the context of classical regression models, distinguishing it from full linearization.16-19 Simulation results and theoretical arguments indicate that eq 11 is an excellent approximation if the percentage variance (%Var) of the spectral matrix that is explained by the model is high.20 We will show that the %Var is very high for all models used in the current application; thus, the approximate nature of eq 11 is of no practical importance here. However, %Var is an important applicationdependent criterion that must be monitored. For example, the currently advocated approach is likely to fail for a recently published data set where the optimum model explains just 60% of the spectral variance.27 Training samples in a given application range from being a pure random draw from a population (“natural calibration” using random variables) through following a strict design (“fixed calibration” using mathematical variables). Since eq 11 gives the prediction error variance conditional to the training data, how the training samples are defined is irrelevant: regardless of their status in the training stage, they are fixed (“frozen”) in the prediction stage. The interpretation of the contribution of the spectral measurement variance to eq 11 is greatly enhanced by the use of Lorber’s definition for multivariate sensitivity.28,29 This definition was originally formulated within the classical calibration framework, where spectral responses are modeled as a function of analyte concentration; however, the following relation holds between the size of the regression vector and the multivariate sensitivity (SEN):30

||b|| ) 1/SEN

(13)

The size of the regression vector is thus interpretable as the “inverse sensitivity” in the inverse model. If the multivariate net analyte signal (NAS) and (inverse) sensitivity are consistently defined, a multivariate model can be accurately represented as a univariate calibration graph (see Figure 6 below and Figure 17 of ref 31). (27) Ho ¨skuldsson, A. J. Chemom. 1995, 9, 91-123. (28) Lorber, A. Anal. Chem. 1986, 58, 1167-1172. (29) Lorber, A.; Faber, N. M.; Kowalski, B. R. Anal. Chem. 1997, 69, 16201626. (30) Lorber, A.; Kowalski, B. R. J. Chemom. 1988, 2, 67-79.

Analytical Chemistry, Vol. 70, No. 14, July 15, 1998

2975

The standard uncertainty in the estimated NAS of the unknown sample is given by31

σ(NA ˆ Su) ≈ σ∆ru

(14)

A region on the NAS axis is mapped onto a region on the analyte concentration axis, with the mapping amounting to a simple multiplication with the inverse sensitivity. It follows that the uncertainty in the NAS estimate of the unknown sample is converted to a contribution to the uncertainty in the predicted analyte concentration by multiplication with the inverse sensitivity:

contribution to V(PEu) ) (σ(NA ˆ Su)/SEN)2 ≈ (σ∆ru/SEN)2 (15)

which is consistent with eq 11. Interpreting abstract PE variance in terms of multivariate analytical figures of merit makes the mathematical form of eq 11 identical to that for univariate calibration with measurement errors in analyte concentration and a single spectral response.31 Practical Evaluation of Approximate Expression for Prediction Error Variance. Evaluation of eq 11 requires replacing unknown quantities with practical analogues. These unknown quantities can be divided into two classes: quantities that would (probably) be calculated anyway and those requiring additional experimental effort. Estimates for the regression vector are readily available, and the unknown sample leverage can be calculated using eq 12. The variances associated with measurements (concentrations as well as responses) have to be estimated from replication experiments that are independent of the data used to construct the model. The variance of the residuals is determined as follows. In the absence of prediction bias, the mean-squared error of calibration (MSEC) is given by21 2 2 MSEC ) σ2e + σ∆c + ||bA||2 σ∆R

(16)

In practice, MSEC is calculated from the squared fit errors for the calibration samples as

MSE ˆC )

1

Figure 1. Prediction uncertainty versus selected number of factors. The RMSEP (O) is the root-mean-square of the bias (B) and variance (V) components. The “optimum” model dimensionality is a compromise between bias (too few factors) and variance (too many factors).

While eq 16 contains only variances and, hence, assumes prediction bias to be absent, eq 17 includes squared bias in the case of underfactoring and thus introduces a bias term into eq 18. Inserting the result of eq 18 into eq 11 implicitly accounts for prediction bias if training and prediction samples are exchangeable. Formally, this approach provides a sample-specific meansquared error of prediction rather than a prediction error variance. Practically, simulations indicate that the prediction intervals obtained using this approach provide appropriate coverage (not shown). Selection of Optimum Model. In the preceding sections, the focus was on the variance in the predictions. However, if “biased regression techniques” such as PCR and PLS are used, attention must be given to prediction bias as well as prediction variance. Figure 1 shows how prediction variance and prediction bias are combined to a root-mean-square error of prediction (RMSEP). Constructing a model with an insufficient number of factors leads to underfitting bias, whereas constructing a model with too many factors gives overfitting variance. The “optimum” model dimensionality (number of factors) is thus a compromise between bias and variance. Analogous to eq 17, the total error of prediction for a particular model is calculated as

n

∑(cˆ - ˜c ) n-A-1

2

i

i

(17)

i)1

MSE ˆP )

1 nt

nt

∑(cˆ - ˜c )

2

i

i

(19)

i)1

The factor n - A - 1 accounts for the loss of degrees of freedom due to an A-dimensional model fit to the mean-centered data. While this adequately accounts for the loss of degrees of freedom with PCR, a further correction may be needed for PLS, since both spectral responses as well as analyte concentrations are used to construct the model.17,18 Rearranging eq 16 and replacing true values by their estimates gives the desired estimate for the variance of the residuals,

where nt denotes the number of test samples. Correction for loss of degrees of freedom is not required, since the summation is over samples that were not used for model building. A suitable compromise between bias and variance may be the first minimum in the RMSEP curve or a plateau.11 Equations 17 and 19 define error with respect to the measured rather than the true values, thus overestimating the true errors because the random uncertainty part of the measured values cannot be modeled.

2 2 σˆ 2e ) MSE ˆ C - σˆ ∆c - ||bˆ A||2 σˆ ∆R

EXPERIMENTAL SECTION Note: Certain commercial equipment, instruments, and materials are identified in this report to specify adequately the

(18)

(31) Faber, N. M.; Lorber, A.; Kowalski, B. R. J. Chemom. 1997, 11, 419-461.


experimental procedure. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose. Spectral Data. Full details on the spectrometer and spectra acquisition parameters are presented elsewhere.8 Briefly, nearIR absorbance spectra were acquired with a Bruker IFS66 Fourier transform infrared (FT-IR) spectrometer (Billerica, MA). All results presented here are based on 391 absorbance values evenly spaced in wavenumber space between 6000 and 9000 cm-1. All spectra were acquired within a 6-h period. All ampules were spun at 400 rpm to reduce baseline variation. In addition, several spectral preprocessing methods were applied to eliminate subtle substructures unrelated to chemical information (among others, first- and second-derivative). All results presented here are obtained using multiplicative scatter correction.32 Materials and Sample Preparation. A description of the industry-average gasoline (“RF-A”) and the MTBE and EtOH used to prepared the mixed-oxygenate samples is presented elsewhere.8 Three EtOH + H2O solutions were prepared to represent low, medium, and high levels of H2O contamination, using NIST housedistilled H2O. Preparation was randomized within each level. The H2O levels present in the RF-A, MTBE, and EtOH stock materials were not determined. All samples were prepared over a 3-day period, with a single EtOH + H2O level used throughout each day. Two different batches of RF-A obtained from the same supplier were used, one batch used throughout all of days 1 and 2 and through half of the samples of day 3. Samples were prepared gravimetrically: target volumes of the three components (RF-A, MTBE, and EtOH + H2O mixture) were delivered into an Ar-purged ampule placed on an analytical balance. The weight of each added component delivered was determined after a short equilibration interval. After delivery of all components, the ampule was cooled in liquid N2 and flamesealed under an Ar stream. Safety Considerations. The gasoline and organic oxygenates used are flammable and harmful to varying degrees with respect to inhalation, ingestion, and eye irritation. Standard laboratory precautions were taken to minimize exposure to liquids and vapors. Standard precautions were used in the sealing, transport, and storage of the ampuled materials. Sample Sets. The experimental design used to identify the target composition of the samples is discussed in the following section. Figure 2 displays the distribution of analyte levels among the 80 “normal” samples that were prepared. These 80 samples were sorted in order of increasing MTBE level and divided into two groups of 40 samples each by taking every second sample, one group designated “A” and used exclusively as the calibration set and the other group designated “B” and used exclusively as a test set. Table 1 summarizes the size, composition, and use of the sample sets. Two samples, designated “C1” and “C2”, that were expected from their preparation history to be unusual were included in the study to explore model performance. Sample C1 was the first ampule prepared and visibly evaporated while production “bugs” (32) Geladi, P.; McDougal, D.; Martens, H. Appl. Spectrosc. 1985, 39, 491-500.

Figure 2. Distribution of samples (calibration set A, and test sets B and C) over concentration range. Table 1. Description of Data Sets data set A B C

range of mass fraction %O delivered by analyte type

sizea

MTBE

EtOH

H 2O b

“normal” 40 × 391 0.1669-2.9389 0.0387-1.5190 0.0008-0.0655 “normal” 40 × 391 0.3501-3.0694 0.0386-1.4840 0.0013-0.0746 “suspect” 2 × 391 3.0165, 0.4744 0.0423, 0 0.0007, 0

a Samples × wavenumbers (unit, cm-1; ranging between 6001.4 and 9009.8, with spacing 7.7). b H2O present as contaminant of EtOH stock material is not included.

were worked out. Sample C2 is an “honest blunder”: immediately following preparation, it was noted that the EtOH + H2O mixture either had not been added or that the weight of the mixture had not been recorded. Prediction results presented below confirm the first explanation. Calculations. All calculations were performed in Matlab Version 4.2c (The Mathworks, Natick, MA). RESULTS AND DISCUSSION Experimental Design. The following design criteria were considered: (1) the number of the calibration samples must be sufficiently large and must be constructed with diverse enough analyte (and interference) concentrations in multivariate space to obtain a stable model with good predictive ability, (2) the number of test samples must be sufficiently large and must not just duplicate calibration samples to avoid obtaining overly optimistic prediction uncertainty estimates, (3) the target values of 2.0% O and 2.7% O should be well covered, (4) EtOH, which is typically added to a gasoline at the point of delivery, is likely to be the minor component of a mixed-oxygenate gasoline, and (5) the H2O concentration should not be unrealistically high. The influence of the number of calibration samples, n, on prediction uncertainty can be deduced from the first term on the right-hand side of eq 11. The variance in the model center is inversely proportional to n. The role of the unknown sample leverage is evaluated as follows. If the calibration samples cover the prediction range well, unknown samples will lie on interpolatAnalytical Chemistry, Vol. 70, No. 14, July 15, 1998

2977

Table 2. Summary of Prediction Results for Previously Constructed PLS Models single-oxygenate models

global model

dual-oxygenate model

analyte

A

RMSEP %O

RMSECV %O

A

RMSEP %O

RMSECV %O

A

RMSEP %O

RMSECV %O

MTBE EtOH

1 1

0.011 0.014

0.016 0.017

8 4

0.030 0.009

0.051 0.009

6 7

0.015 0.031

0.009 0.025

ing positions, and the leverage will be of the same magnitude as the average leverage of the calibration set. This average value, hh , follows by combining eqs 5 and 12:

1 1 Tr[RWAWATRT] ) Tr[TATAT] n n 1 A ) Tr[HA] ) n n

prediction uncertainty estimates is thus given by

σ(MSEP) x2nt ) ) MSEP nt

hh )

(20)

where Tr[‚‚‚] symbolizes the trace operator (sum of diagonal elements of a square matrix), and HA (n × n) denotes the “hat matrix” that converts the concentration vector of the calibration samples into the fitted values, cˆ A ) HAc. The hat matrix HA projects down from an n-dimensional space to an A-dimensional space. For a projection matrix, the trace is equal to its rank, which is given by the dimension of the vector space it projects down to, i.e., A. Thus, the prediction leverage will also be (approximately) inversely proportional to n. If the ratio A/n is small enough, on the order of 1/5 or less, adding samples to the calibration set will not lead to an appreciable reduction of the prediction uncertainty because the last two terms in eq 11, which are constant, start to overwhelm the model contribution. An idea of the model dimensionalities (A) that can be expected to be adequate for the current work is obtained from the prediction results obtained during our preliminary study.8 These are summarized in Table 2. The single-oxygenate models used were constructed from samples containing only MTBE or EtOH, the global model was constructed by merging the single-oxygenate sample sets, and the dual-oxygenate model was constructed using mixtures of MTBE and EtOH. (The RMSEP and RMSECV in Table 2 differ from those previously reported by the coverage factor used to account for batch heterogeneity.) As many as eight factors are necessary to obtain the optimum global model for MTBE; this compares well with the numbers reported for similar analyses by others.3,33 Thus, based on the information available at the start of the current work, we expected models of approximately eight factors. Consequently, 40 calibration samples were prepared. The number of test samples is derived as follows. Prediction uncertainty estimates such as RMSEP are highly variable themselves. If PEs are normally distributed, the MSEP will have a χ2 distribution35 leading to an expected value proportional to nt and a variance proportional to 2nt. The relative uncertainty in the (33) Swarin, S. J.; Drumm, C. A. Spectroscopy 1992, 7, 42-49. (34) Peng, C.; Lewis, K. C.; Stein, F. P. Fluid Phase Equilib. 1996, 116, 437444. (35) Stallard, B. R.; Garcia, M. J.; Kaushik, S. Appl. Spectrosc. 1996, 50, 334338.


x

2 nt

(21a)

The relative uncertainty in RMSEP follows from a linear approximation as half the relative uncertainty in MSEP:

σ(RMSEP)/RMSEP ≈ )

|∂RMSEP/∂MSEP|σ(MSEP) RMSEP (1/2RMSEP)σ(MSEP) ) RMSEP

x

1 2nt

(21b) Hence, 40 test samples enable a reasonably precise ([1/(2 × 40)]1/2 ) 0.112 ≈ 11%) estimate of prediction uncertainty. While perhaps more test samples than strictly necessary, the few large differences between RMSEP and RMSECV in Table 2 can be attributed to a too small test set. Given the 2.0% O and 2.7% O mandated threshold levels, the samples were designed to span the restricted range of 1.6-3.1% O. For the purposes of this demonstration study, it was assumed that a mixed MTBE-EtOH oxygenate system would arise from augmenting a 2.0% O MTBE gasoline with EtOH to reach the 2.7 total %O level. The EtOH contribution was hence constrained to contribute no more than 1.5% O. An approximate space-filling distribution of target MTBE and EtOH levels was generated under these constraints. The H2O levels were assigned by assuming 0.6, 1.2, or 1.9 mass fraction %H2O “contamination” of the added EtOH, giving a maximum 0.08% H2O level in the oxygenated gasoline. Higher levels of H2O approach the cloud-point for 2.0% O MTBE/EtOH mixtures at typical northern United States winter temperatures.34 Figure 3 is a principal component score plot11 that further illustrates that the design is “biased” in the sense that most samples lie in the region of medium to high MTBE content. The splitting of the “normal” samples in a calibration and test set is reasonably random. Sample C1 is found to be close to two calibration samples, while sample C2 is clearly identified as an outlier. Characterizing the Measurement Error in the Data. Two measurement errors have to be determined (independently from the data used to create or validate the model) to enable estimation of the variance associated with concentration residuals using eq 18: the standard deviation in analyte concentrations and the standard deviation in the near-IR absorbances. The analyte concentration standard deviation for the calibration samples is

Figure 3. Distribution of samples in plane spanned by PC scores 1 and 2: calibration set A (×), test set B (O), and test set C (C1 and C2).

derived from the sample preparation weighing errors:

σ(wt % O) ) σˆ (wt % analyte) ) 100 ×

(

MW(O) MW(analyte)

)

σˆ (weight analyte)/weight sample (22) MW(analyte)/MW(O)

where MW denotes the molar weight of the associated compound. Inserting the target weight of the samples (15 g), the molar weight of the analyte (88.15, 46.07, and 18.02 g for MTBE, EtOH, and H2O, respectively) and the molar weight of oxygen (16 g) in eq 22 leads to conversion factors for the three analytes of 1.21, 2.32, and 5.92, respectively. Replication experiments suggest that the weight-by-difference sample preparation procedure, in expert hands, has an estimated precision of σ(weight analyte) ) 2 × 10-5 g with low-volatility materials. Since all sample components (gasoline, MTBE, and EtOH + H2O mixture) are volatile, we use a pessimistic estimate of σ(weight analyte) ) 10-4 g for all components. The standard deviation in the background-corrected near-IR absorbances was estimated by replication to be σˆ ∆R ) 10-5 AU. This value coincides with a literature value for a Nicolet 800 Fourier transform spectrometer.35 We use σˆ ∆R ) 2 × 10-5 AU as a pessimistic estimate. Selection of Optimum Model. Only the results for MTBE are discussed in detail; similar observations are made for EtOH and H2O. Figure 4 displays RMSEC and RMSEP for MTBE for both PCR and PLS models. The PCR RMSEC reaches a stable plateau, suggesting that underfitting bias is negligible after extracting a certain number of factors, while the overfitting variance increases only slowly. However, there is a relatively large divergence observed between RMSEC and RMSEP for the fifth PCR component, indicating that the PE is overly optimistic for the five-dimensional model. We conclude that the “optimum” PCR model is six-dimensional. The PLS RMSEC steadily decreases, suggesting that the simple correction for degrees of freedom in eq 16 is insufficient for high-dimensional models (A > 12). Fit errors (RMSEC)

Figure 4. RMSEC (/) and RMSEP (O) obtained for MTBE versus model dimension for PCR and PLS. Table 3. Variance Components of Mean-Square Error of Calibration for Selected “Optimum” Models analyte method

A

MSEC (%O)2

2 σ∆c (%O)2

(σ∆R/SEN)2 (%O)2

σ2e (%O)2

MTBE

6 5 12 8 10 7

3.4 × 10-4 3.4 × 10-4 7.0 × 10-6 7.7 × 10-6 8.6 × 10-6 8.4 × 10-6

1.5 × 10-8 1.5 × 10-8 5.4 × 10-8 5.4 × 10-8 3.5 × 10-7 3.5 × 10-7

7.8 × 10-7 7.7 × 10-7 2.0 × 10-7 1.8 × 10-7 5.4 × 10-8 4.9 × 10-8

3.4 × 10-4 3.4 × 10-4 6.7 × 10-6 7.4 × 10-6 8.2 × 10-6 8.0 × 10-6

EtOH H2O

PCR PLS PCR PLS PCR PLS

becoming significantly smaller than the PEs (RMSEP) is a manifestation of overfitting: PLS is known to be more prone to overfitting than is PCR.25 Once all relevant structure is extracted from the concentration vector, only noise is built into subsequent score vectors. These additional score vectors do an exceptionally good job of “recognizing” the concentrations of the calibration set, but they cannot improve prediction for true unknowns. The PLS RMSEP gradually decreases for MTBE when going from a 5-dimensional to a 10-dimensional model; however, given the uncertainty in the RMSEP, this small decrease is probably not significant. Moreover, including these additional factors may give models that are not robust to small changes in RF-A composition. These considerations lead to the following optimum model dimensions for MTBE, EtOH and H2O: 6, 12, and 10 PCR dimensions and 5, 8, and 7 PLS dimensions. The results of crossvalidation (eight groups of 10 calibration samples each) are generally consistent with the prediction results for the independent test set. Characterization of Residuals. Estimating the standard deviation of the residuals amounts to breaking up MSEC into the individual contributing variance components of eq 16. Table 3 summarizes the results of this calculation, with the spectral variance component calculated from eq 15 and the pessimistic Analytical Chemistry, Vol. 70, No. 14, July 15, 1998

2979

Table 4. Summary Statistics for Selected “Optimum” Models analyte method A MTBE EtOH H2O


6 5 12 8 10 7

RMSECa RMSEPa %O %O 0.018(2) 0.018(2) 0.0026(4) 0.0028(4) 0.0029(4) 0.0029(4)

0.017(2) 0.015(2) 0.0035(4) 0.0036(4) 0.0023(3) 0.0024(3)

SEPa %O

SEN AU/%O

%Var

0.020(2) 0.020(2) 0.0031(4) 0.0031(4) 0.0033(4) 0.0031(4)

2.3 × 2.3 × 10-2 4.5 × 10-2 4.8 × 10-2 8.6 × 10-2 9.1 × 10-2

99.944 99.856 99.997 99.969 99.992 99.953

10-2

a The digits in parentheses denote standard deviation in units of the last reported digit of the preceding value.

Figure 5. RMSEC (/) and RMSEP (O) obtained for “optimum” MTBE PLS models (all of model dimension ) 5) as a function of the standard deviation of the random noise added to the spectral data.

values for measurement error standard deviations used to evaluate eq 18. The PCR and PLS results are highly consistent. For all models considered, the variance component associated with the residuals dominates the variance components associated with the measurement errors: the smallest ratio is found to be 8.0 × 10-6/ 3.5 × 10-7 ) 23. Since pessimistic values were used, we conclude that the measurement errors play a negligible role in this application: i.e., eq 1 is an adequate description for these particular experimental data sets and models. The adequacy of eq 1 is certainly not a general result for nearIR applications. The measurement error in the calibration set analyte concentration is usually much larger. Here, the analyte concentrations were determined gravimetrically, while in most applications reference values are obtained using less accurate techniques. Also, the spectral errors do not significantly contribute to PE variance, owing to the relatively small regression vector, which is equivalent to a relatively high multivariate sensitivity. This ensures a small amount of spectral error propagation in eq 11 for calibration and prediction samples alike. The norms for the regression vectors estimated by PCR and PLS are approximately 40, 20, and 10 (numerical values are sensitive to units used) for MTBE, EtOH, and H2O, respectively. We believe this to be exceptional for near-IR data, but the assertion is difficult to verify since norms of regression vectors or multivariate sensitivities are seldom reported in the literature. Measurement errors in the test spectra made a large contribution to PE variance in one recently reported application.36 Additional evidence for the conclusion that spectral measurement errors hardly contribute to PE variance is obtained by investigating the effect of adding random noise to the spectral data. Figure 5 displays PLS RMSEC and RMSEP for the selected “optimum” MTBE model as a function of the standard deviation of the noise artificially added in. The RMSEP remains approximately constant up to an added noise standard deviation of 6 × 10-5 and then rises slowly. Similar results are obtained for EtOH and H2O, with one additional factor needed for H2O when the added noise standard deviation exceeds 12 × 10-5. Since the pessimistic estimate for intrinsic spectral error is 2 × 10-5 (36) Berger, A. J.; Feld, M. S. Appl. Spectrosc. 1997, 51, 725-732.


AU, one-fourth of the noise level found to “make a difference” in the PCR and PLS noise addition studies, and considering that (independent) uncertainties are added in quadrature, the spectral measurement errors are, indeed, negligible for the current application. Since the effect of measurement errors is negligible, an approximate 100(1 - R)% coverage prediction interval (PI) for the true analyte concentration cu is

PIR(cu) ) cû ( tR/2,n-A-1V(PEu)1/2 ≈cû ( tR/2,n-A-1σˆ e(1 + 1/n + hu)

(23)

where tR/2,n-A-1 is the critical value of a t-distribution with n - A - 1 degrees of freedom. Because of the approximate character of eq 11, eq 23 is also approximate. With negligible measurement errors in the reference values, eqs 17 and 19 simplify to the form in which they are usually presented in the literature, with true values rather than measured values used for reference values. The terms “apparent” and “actual” prediction errors have been used to emphasize that such simplification overstates the PE if measurement error is appreciable.37 A procedure has recently been proposed to correct “apparent” PE estimates for the measurement error in reference values.38 Description of “Optimum” Models and Prediction Results. Table 4 summarizes results for the selected “optimum” models. The current prediction uncertainty results are comparable to the previous study’s (Table 2) for MTBE but are much better for EtOH. It was expected that the results for singleoxygenate samples would constitute a lower bound on prediction uncertainty for more complex situations. It was not expected that the predictions for three-component mixtures can be as good as (for MTBE) or better than (for EtOH) those for the simpler test cases. For H2O, a minor component in these mixtures, a reasonable estimate for the limit of detection is 3 × RMSEP ≈ 0.01% O. We attribute the excellent results currently obtained (especially for EtOH and H2O) to the adequate size and design of the calibration set. The standard error of prediction (SEP) values in Table 4 are calculated by averaging the PE variances from eq 11. If the input standard deviation of the error in the gravimetry, σˆ ∆c, is changed from 0 (ideal) to 10-4 (pessimistic estimate), the SEP values (37) DiFoggio, A. H. Appl. Spectrosc. 1995, 49, 67-75. (38) Faber, N. M.; Kowalski, B. R. Appl. Spectrosc. 1997, 51, 660-665.

Table 5. Prediction Results for 40 Samples of Test Data Set B Using Selected “Optimum” Models sample-specific prediction SD (mass fraction %O) analyte method A MTBE EtOH H2O


6 5 12 8 10 7

minimum

maximum

0.019 0.019 0.0028 0.0029 0.0031 0.0030

0.022 0.022 0.0037 0.0034 0.0039 0.0038

Table 6. Prediction Results for Test Data Set C Using Selected “Optimum” Models

leverage leverage

minimum maximum 0.065 0.022 0.110 0.067 0.106 0.068

0.413 0.345 0.921 0.513 0.760 0.691

change only in the fourth decimal digit. The standard deviations reported for RMSEC, RMSEP, and SEP are calculated using eq 21b. The number of degrees of freedom for RMSEP is nt ()40), but RMSEC and SEP have to be corrected for loss of degrees of freedom from estimating the model center and A factors, n - A - 1 ()39 - A). These standard deviations on the criteria help “explain” some of the irregularities in Table 4. It is usual to have some RMSEP values that exceed the corresponding RMSEC values, since the leverage of prediction samples tends to be larger than the average for the calibration set; none of the few such cases observed here are significant. The estimates of RMSEC and SEP are functionally related, unlike the estimates of RMSEC and RMSEP, with overestimation of RMSEC leading to a pessimistic estimate of SEP and vice versa. The variability in the uncertainty estimates is satisfactory for the current application (between 11 and 13%) because of the relatively high degrees of freedom. Table 4 also lists the multivariate sensitivities calculated as the inverse of the norm of the regression vector (eq 13) in units of AU/%O and the %Var obtained by reconstructing the spectral matrix from the model. The (multivariate) sensitivity of the instrument for the determination of these analytes is clearly ordered as H2O > EtOH > MTBE. Only slight differences between the PLS and PCR values are found; all models capture virtually 100% of the spectral variation, justifying the “naive” or “zeroth-order” approximation leading to eq 11. Results for full linearization have been calculated and typically differ with the results of this much more easily implemented, numerically robust, and interpretable approximation in the fourth decimal digit only. Table 5 summarizes prediction results for the test set samples. The predicted uncertainties and leverages are not homogeneous, confirming that calculating a constant prediction interval on the basis of average prediction uncertainty criteria such as RMSEP or RMSECV is not satisfactory. For example, the sample-specific prediction uncertainty for high-leverage samples is 15% larger than SEP; a SEP-based nominal 95% uncertainty interval would thus provide only about 90% coverage for these samples. Table 6 gives detailed prediction results for the two suspect samples, C1 and C2. Figure 3 shows that sample C1 is (probably) a normal sample, in the sense that it lies in the range spanned by the calibration set, whereas sample C2 is clearly an outlier. While of “normal” composition, C1 noticeably evaporated before the ampule was flame-sealed; the small PEs obtained for this sample suggest that this reduction in sample volume did not significantly change the relative composition of the mixture. The observed PEs are much larger for sample C2. They are, in many cases, larger than 3 × RMSEP, and they are more than 6 × RMSEP for

analyte method A MTBE EtOH H2O


6 5 12 8 10 7

C1 0.35 0.22 0.45 0.43 0.42 0.35

C2

observed prediction error (%O) C1

2.81 -0.007 2.59 -0.017 4.45 0.0028 2.85 0.0034 4.23 -0.0000 2.80 0.0000

sample-specific prediction SD (%O)

C2

C1

C2

-0.072 -0.056 -0.0229 -0.0177 0.0139 0.0088

0.022 0.021 0.0032 0.0033 0.0035 0.0034

0.036 0.035 0.0062 0.0055 0.0067 0.0056

the H2O PCR model. An average prediction uncertainty is not appropriate for this sample. A sample-specific prediction uncertainty, such as calculated from eq 11, is required. For example, the PE for the H2O PCR model is only slightly above 3 times the predicted standard deviation given in the last column. The same holds true for the EtOH PLS case. Given the extremely large leverages (more than 14 times the average for the calibration set) for these cases, all prediction resultssuncertainties as well as concentrationsshave to be considered with caution, since they were obtained in a region where the model might not be valid. Nevertheless, given that nonaddition of the EtOH + H2O mixture was, indeed, the production blunder, the proposed sample-specific prediction intervals also work satisfactorily for this difficult example. The differences in leverage for the same sample are due to the differences in model dimensionality needed for optimal prediction of the different analytes. Equation 20 shows that the average leverage for the calibration set is proportional to the model dimensionality A. It is reasonable to assume that such a dependency also approximately holds for prediction samples. Further insight can be obtained from a graphical representation of the model: plotting regression coefficients versus wavenumber explicitly provides information on all variables. Such plots are useful for diagnostic purposes, since problematic wavenumbers lead to irregularities (e.g., spikes) that are readily identified. No irregularities were observed in the estimated regression vectors (not shown); consequently, all wavenumbers are on an equal footing and can be lumped together using summary statistics such as the multivariate sensitivity.28 Figure 6 presents “univariate” calibration graphs for MTBE, EtOH, and H2O. The vertical projection of the experimental points onto the model line gives the exact observed PEs. Figure 6, therefore, provides the same information as the usual plots of predicted versus “known” values, as well as illustrating the effect of spectral measurement error on PE variance as quantified by eq 15. Visual inspection demonstrates that a horizontal shift of the experimental points over a distance of 2 × 10-5 (pessimistic spectral error estimate) has no palpable effect on the vertical projections. It is stressed that these plots describe models that are constructed using spectral data that have been pretreated by multiplicative scatter correction to remove variable background patterns. Other spectral pretreatment methods have been tested for the current application. For example, taking derivatives of the spectral data can remove small substructures unrelated to analyte concentration. Taking the first derivative will remove a Analytical Chemistry, Vol. 70, No. 14, July 15, 1998

2981

Figure 6. Mass fraction %O versus net analyte signal calculated using “optimum” PLS models for MTBE (model dimension ) 5), EtOH (dimension ) 8), and H2O (dimension ) 7): calibration set A (×) and test set B (O).

constant background, while taking the second derivative will eliminate a sloping background. It was found that taking derivatives leads to substantially lower sensitivities (higher inverse sensitivities) for the current application. Especially for EtOH, the sensitivity declined to the point that spectral error became the dominant PE variance component of eq 11. Plots similar to Figure 6 but displaying results for a single analyte using the same horizontal and vertical axis ranges across subplots readily illustrate the influence of spectral error on PE variance as a function of data pretreatment (not shown).


CONCLUSIONS A recently developed expression for PE variance has successfully characterized the uncertainty in near-IR spectroscopic predictions of mixed-oxygenate concentrations in gasoline. Samplespecific prediction intervals can be calculated, an immediate advantage over the use of traditional uncertainty estimates such as RMSEP that provide only a single (constant) prediction interval for all future unknowns. For the current application, nominal 95% coverage prediction intervals based on RMSEP provide less than 90% coverage for high-leverage “normal” samples. Calculation of a realistic prediction interval for each unknown sample enables the identification of outliers in a large batch of SRM for which the majority can be safely assumed to meet the requirements. Exclusion of these outliers from the batch will reduce the coverage factor needed to account for the batch heterogeneity component of the expanded uncertainties. The practical application of the proposed expression requires identification of all significant sources of measurement error (in concentration as well as spectra). While this is common practice in univariate calibration, until now equivalent data have seldom been collected and/or reported in applications involving multivariate calibration. Correct prediction uncertainty results can be expected to result only from stable models that are validated using a sufficiently large test set. Thus, practical guidelines for choosing the number of training and test samples were derived which have worked well in the current application. Since complicated expressions for PE variance are not appealing to the “practical analyst”, the interpretation of the spectral error contribution to the overall prediction uncertainty is enhanced by the use of an analytical figure of merit, i.e., the multivariate sensitivity. This interpretation enables a simple visualization of the multivariate model as a univariate calibration graph (Figure 6). The spectral error contribution is negligible for all analytes considered in the current application. ACKNOWLEDGMENT We thank Bruce Kowalski for initiating this collaboration between theory and practice, Tom Dean for proposing noise addition for evaluating the significance of spectral measurement error, and Marjan Faber-Meinders for her help in converting formulas into insight (and for tolerating a Maryland summer with grace and humor).

Received for review November 19, 1997. Accepted April 29, 1998. AC971270B

Characterizing the Uncertainty in Near-Infrared Spectroscopic

Recommend Documents