Translator Disclaimer
1 November 2007 Determination of uncertainty in parameters extracted from single spectroscopic measurements
Author Affiliations +
Abstract
The ability to quantify uncertainty in information extracted from spectroscopic measurements is important in numerous fields. The traditional approach of repetitive measurements may be impractical or impossible in some measurements scenarios, while chi-squared analysis does not provide insight into the sources of uncertainty. As such, a need exists for analytical expressions for estimating uncertainty and, by extension, minimum detectable concentrations or diagnostic parameters, that can be applied to a single noisy measurement. This work builds on established concepts from estimation theory, such as the Cramér-Rao lower bound on estimator covariance, to present an analytical formula for estimating uncertainty expressed as a simple function of measurement noise, signal strength, and spectral overlap. This formalism can be used to evaluate and improve instrument performance, particularly important for rapid-acquisition biomedical spectroscopy systems. We demonstrate the experimental utility of this expression in assessing concentration uncertainties from spectral measurements of aqueous solutions and diagnostic parameter uncertainties extracted from spectral measurements of human artery tissue. The measured uncertainty, calculated from many independent measurements, is found to be in good agreement with the analytical formula applied to a single spectrum. These results are intended to encourage the widespread use of uncertainty analysis in the biomedical optics community.

1.

Introduction

Real-time analysis of spectroscopic measurements is essential in applications such as pharmacokinetics,1 bioreactor monitoring,2 and medical diagnosis.3 In our laboratory4, 5 and others,6, 7 real-time analysis of spectroscopic measurements acquired in vivo is under study to provide clinicians with immediate diagnoses, in lieu of histopathology. In medical applications, the confidence in the measurement of a particular diagnostic parameter can affect the course of disease management, with ramifications to the health of the patient. The uncertainties and associated confidence intervals of the parameters extracted from spectroscopic measurements serve to assess the accuracy, stability, and diagnostic value of the data. The importance of uncertainty is related to other figures of merit commonly mentioned in the chemometrics field: signal-to-noise ratio, precision, limit of detection, sensitivity, error propagation, and selectivity.8, 9 Note that measurement uncertainty (precision) is independent of measurement accuracy.

The most effective way to extract quantitative information from spectral data in a linear system is by utilizing the full spectrum (multivariate analysis).10 Consider, for example, measurement of the concentration of a particular species or analyte. This requires a model that, when applied to a measured spectrum, yields the concentration of interest. In most cases, the model can be conveniently expressed in terms of a regression spectrum or “b-vector” for a particular analyte; the analyte concentration of a prediction sample can then be expressed as the inner product of the measured spectrum and the b-vector. When all of the chemical components are known, the model can be based on the constituent spectra, measured directly, and ordinary least squares (OLS) analysis can be applied, yielding a b-vector for every component of interest. If the spectra cannot be measured directly, or if one or more components are not known, a calibration step is required to generate the b-vectors, and a direct calibration scheme such as classical least squares (CLS), or indirect calibration schemes such as partial least squares (PLS) or principal components regression (PCR) can be used. The calibration step requires a set of spectra with reference concentrations of the analyte(s) of interest. Similar approaches can be used to measure parameters extracted from biological tissue spectra that are used in disease diagnosis.11, 12 A concept closely related to the b-vector is the net analyte signal (NAS), introduced by Lorber,13 and extended by Lorber, Faber, and Kowalski,14 which is the portion of the signal for each analyte that is orthogonal to the other analyte spectra. The NAS is also useful in evaluating the figures of merit mentioned earlier.14

In principle, one can evaluate the parameter uncertainty by repeating the measurement many times and analyzing the standard deviation of parameters extracted from each of these multiple measurements. However, this is not practical for applications such as medical diagnosis, in which only one or a few measurements can be acquired. Alternatively, one can use chi-squared (χ2) analysis to calculate parameter uncertainties extracted from a single spectrum.15 χ2 analysis is a very useful technique, but it is statistical rather than analytical and provides little insight into the origins of uncertainty.

In this work, we adopt an alternate approach to uncertainty analysis and here present a method of analysis that can be used in conjunction with a single spectrum to provide physical insight into the sources of uncertainty. The analytical expression employed for this purpose describes concentration uncertainty as a function of measurement noise, signal strength, and spectral overlap—quantities easily extracted from spectroscopic measurements. As such, the method can guide improvements in data modeling, as well as the optimization of the instrument. This approach can be considered as an extension and a complement to previous work of our research group.16 In that study, we derived an analytical expression for the limiting uncertainty in analyte concentrations extracted from Raman spectra using PLS, and showed it to hold experimentally.16 Uncertainty was expressed as a function of measurement noise and the b-vector using PLS. Since then, we have worked on characterizing the lower bound on the uncertainty in extracted concentrations using a more generally applicable approach.17 Other groups have employed our framework in analyzing the uncertainties and sources of error in spectroscopic measurements.18

Our approach is a special case of the error analysis of Lorber and Kowalski,19 where calibration is very accurate and thus model uncertainty is negligible. Lorber and Kowalski have presented a complete and elegant treatment of error propagation associated with multivariate calibration. They derived a prediction error formula that depends on the noise in the spectrum of the prediction sample and the spectra and concentrations of the analyte of interest in the calibration set. The formula was successfully tested on near-infrared reflectance data analyzed by PCR.19 One practical shortcoming of the formula is its complexity and the difficulty in readily applying it to experimental data. Other groups have taken similar approaches to error analysis.20, 21, 22 We focus here on the important case in which the spectra and concentrations of the calibration dataset are measured more precisely than that of the prediction sample. This is often the case, because variables such as integration time may be increased or optimized for the calibration data. Therefore, in the limit in which calibration noise is small, uncertainty will be dominated by measurement noise in the prediction sample.

Although we are focusing on the case when measurement noise dominates, the analytical expression presented here can be used, in many cases, even when model uncertainties are not negligible. Therefore, we may calculate both actual uncertainty, which takes into account modeling and measurement noise, as well as the limiting uncertainty, where model noise is disregarded. The limiting uncertainty, the unavoidable uncertainty associated with the inherent spectral noise in the prediction sample and the spectral makeup of the model, also specifies the smallest concentration at which a constituent can be detected. This approach is applicable to concentration measurements and for determining diagnostic accuracy of parameters obtained for spectral diagnosis of disease, using linear spectroscopic techniques such as Raman scattering or fluorescence, and is very important in system design and evaluation.

In the following, we derive an analytical expression for the limiting uncertainty, and show its equivalence to χ2 analysis (see the Appendix in Sec. 7). We demonstrate that the concentration uncertainty calculated by the analytical formula is in good agreement with that measured experimentally from aqueous solutions of clinically relevant analytes. Furthermore, we demonstrate that in this case the calculated actual uncertainties are very close to the limiting uncertainties, which is indicative of the accuracy of our data acquisition and modeling. To illustrate the biomedical application of this analytical formula, we calculate uncertainties of parameters extracted from tissue spectra that are used in disease diagnosis. These results are intended to encourage the widespread use of uncertainty analysis in the biomedical optics community.

2.

Theory

We adopt a linear algebra approach and notation in this work. All vectors are column vectors and are denoted by bold lowercase letters. Similarly, matrices are bold uppercase letters comprising multiple column vectors, where size is indicated in parentheses (row×column) . Measurements and spectra are denoted as vectors, in that each element represents the response of a particular detector [e.g., charge-coupled device (CCD) pixel]. Many of the results presented below follow from Kay, a standard text in statistical signal processing.23

2.1.

Linear Model

We begin with the standard additive noise linear model:

Eq. 1

s=Pc+w.
The vector s is the observed measurement (M×1) , the matrix P contains the model constituent vectors (M×N) and is full rank, the vector c contains the underlying coefficients of the model constituents (N×1) , and the vector w represents noise in the system (M×1) . Here M represents the number of wavelengths and N represents the number of model components. In other words, the measured vector s is a linear combination of the model components in P , weighted by the coefficients in c , and with the addition of random measurement noise w . The w is assumed to be a zero-mean Gaussian random vector with a known or measurable covariance matrix Cw . Our goal is to analytically determine the best estimate ĉ (fit coefficients) of the true underlying coefficients c and the uncertainty specified by the standard deviation of that estimate.

The application of the model, Eq. 1, to characterize spectroscopic measurements is appropriate. Raman and fluorescence spectra (s) have been shown experimentally to be linear combinations of the distinct individual spectra of the underlying chemical species that can be measured in isolation (P) and proportional to the concentration of those species (c) .24, 25 Measurement noise (w) is typically associated with the measurement system and detector, and its statistics (Cw) can be readily measured; moreover, the Gaussian assumption also holds well in practice.

It should be noted that Eq. 1 is a special case of a more general formulation by Lorber and Kowalski19 that can be specialized for our model as

Eq. 2

s=(P+δP)c+w=Pc+w.
The formulation of Eq. 2 also includes modeling uncertainty δP , taking into account the uncertainties of the concentrations and measurements of the constituent spectra. If the modeling uncertainty is also assumed to be Gaussian, its effect can be lumped together with the measurement noise as w , where the covariance of w is greater than that of w in Eq. 1. As noted earlier, in most of this work we are interested in the limiting case in which the measurement noise of the prediction sample w is the dominant source of uncertainty. In this way, we calculate the limiting uncertainty.

2.2.

Estimator Performance

In determining the optimal estimator, we restrict our attention to unbiased estimators: those that, on average, accurately return the underlying parameters c . From the linear model with the assumptions described before and estimation theory, one can derive the minimum variance unbiased (MVU) estimator:23

Eq. 3

ĉ=(PTCw1P)1PTCw1s.
This estimator is desirable because out of all possible unbiased estimators, it is the one that achieves the minimum variance for all combinations of unknown underlying parameters c . For the MVU estimator of Eq. 3, it can be shown23 that its covariance matrix is:

Eq. 4

cov(ĉ)=(PTCw1P)1.
The diagonal entries of this matrix specify the variances of each ĉk fit coefficient. This is the most general result, as it specifies the covariance and, in turn, the uncertainty of our estimate for any particular noise covariance Cw and the model matrix P .

The result of Eq. 4 can be specialized by assuming that w is white (i.e., uncorrelated and identically distributed) Gaussian noise, so that Cw=σ2I . With this assumption, the MVU estimator and the covariance are given by:23

Eq. 5

ĉ=(PTP)1PTs.

Eq. 6

cov(ĉ)=σ2(PTP)1.
The estimator, Eq. 5, can also be recognized as the OLS solution for ĉ .

One final remark involves the concept of the Cramér-Rao lower bound (CRLB) from estimation theory.23 The CRLB is a lower bound on the covariance of any unbiased estimator. It can be shown that the covariance of the MVU estimator given in Eq. 5 is equal to the CRLB, and hence the estimator is deemed efficient.23 Moreover, for the linear model given before, the efficient MVU estimator implies that it is also the maximum likelihood (ML) estimator. We revisit this last point later.

2.3.

ΔC

The estimation framework and the CRLB described earlier are well-known, general concepts applicable to any type of linear system with the previously mentioned assumptions. However, physical insight can be provided by specializing Eq. 6 to elucidate variables relevant to spectroscopy. We can express P=QS , with S being a diagonal matrix, where the k ’th diagonal entry sk is the Euclidean norm of the k ’th component in P , as follows

Eq. 7

sk=[i=1M(Pi,k)2]12,
and the columns of the matrix Q are thus normalized to unit length. This leads to a simple expression for the standard deviation Δc of the k ’th estimated parameter ĉk :

Eq. 8

Δcstd(ĉk)=σ[(PTP)(k,k)1]12=σsk[(QTQ)(k,k)1]12=σskolfk.
The first factor on the right-hand side, σ , describes the measurement noise and sk quantifies the signal strength of the k ’th model component at unit concentration. The spectral overlap factor olfk indicates the amount of nonorthogonality (overlap) between the k ’th model component and the other (N-1) model components.

The overlap factor may take on values between 1 and . If all of the columns of P (or equivalently, Q ) are orthogonal (no overlap), then olfk=1 . In the other extreme, if the k ’th column of P (or Q ) is nearly linearly dependent with one or more columns (almost complete overlap), then QTQ is close to singular and its inverse does not exist. In a generalized sense, olfk then approaches . (In the case of two columns being linearly dependent, one of the two columns should be removed so that P becomes full rank, as was specified earlier.) In other words, when the model P contains orthogonal constituent spectra, the estimator uncertainty Δc is equal to the ratio (σsk) of the measurement noise to the signal strength for that particular component. In the extreme case of complete spectral overlap (two chemicals with very similar Raman or fluorescence bands across the wavelength range of interest), the estimate is unreliable, so Δc approaches . For the more commonly encountered case of partial spectral overlap of the linearly independent spectral components of P , we have 1<olfk< , and Δc is a function of the three physically understandable quantities: σ , sk , and olfk . The concept of overlap factors as defined here is directly related to variance inflation factors (VIF), first proposed by Marquardt.26 The relation of VIFs to the concept of condition number as well as other relevant methods for evaluation of spectral overlap are described by Kalivas.27 The spectral overlap is also related to the NAS of a particular analyte in that the former measures overlap (interference) while the latter is an indication of nonoverlap (orthogonality).

2.4.

Limiting Versus Actual Uncertainty

We note that of the three parameters, σ varies from sample to sample because shot noise is dependent on the sample-specific raw signal, whereas sk and olfk are sample-independent for a given spectroscopic technique and model.

The value of σ can be obtained in two different ways. One approach is to calculate σ from each pixel across many repeated measurements. The alternative approach, which requires only a single measurement, is to calculate σ from the residual between the observed spectrum and the best fit using Eq. 5. The former value of σ specifies the limiting uncertainty, while the latter value specifies the actual uncertainty through Eq. 8. Therefore, we have

Eq. 9

σlim={1Mi=1M[1L1j=1L(si,js¯i)2]}12,

Eq. 10

σact=[1M(sPĉ)T(sPĉ)]12,
where σlim is calculated as the root mean squared value (across M pixels) of the standard deviation of a representative pixel si calculated from L repeated measurements, while σact is calculated from the residual as the root mean squared between the data s and the fit Pĉ .

Under the assumption of Eq. 1 that the only source of uncertainty is measurement noise, both approaches should yield the same value for σ . However, if there are also modeling uncertainties as in Eq. 2, then only the second approach yields the actual σ , as the residual includes measurement and modeling noise. Therefore, σactσlim , and the difference between the actual and limiting σ can serve to evaluate the accuracy of the modeling. In most of the following, we use Eq. 10 to calculate σ . By extension, we can define Δcact and Δclim using the values of σact and σlim , respectively.

2.5.

Relation to Chi-Squared (χ2)

A statistical method of calculating the uncertainty in extracted parameters can be implemented through χ2 analysis. Maximizing the likelihood of observing a particular measured spectrum (in random Gaussian noise) is equivalent to minimizing χ2 . The value of a parameter that minimizes χ2 is the optimal maximum likelihood (ML) value. The χ2 function (of the underlying parameter) is parabolic in the vicinity of the minimum, and the curvature of the parabola is proportional to the uncertainty (standard deviation) in that parameter.15 Specifically, the variance of the parameter is equal to the reciprocal of the curvature of the χ2 function.15 Equivalently, an increase of one standard deviation of the parameter from the value at the minimum increases χ2 by unity.15, 28 The χ2 approach can be used in conjunction with many fitting procedures, regardless of whether the underlying fitting model is known (or directly measurable) or determinable through calibration.

The analytical expression for Δc , Eq. 8, is equivalent to the one obtained statistically through χ2 analysis. This is to be expected, since the MVU estimator (or the least squares estimator) is equivalent to the ML estimator. Hence the χ2 criterion, used in determining the ML estimate, should yield the same value for the parameter uncertainty as the analytical formula obtained for the MVU estimator. This is demonstrated in the Appendix in Sec. 7.

3.

Methods

In this section, we demonstrate the application of the previous error analysis formalism to the estimation of experimental uncertainty in a set of spectral measurements. Two experiments are performed, both using near-infrared Raman spectroscopy. In the first experiment, we prepare aqueous mixtures of known concentrations of clinically relevant analytes by dilution from stock solutions. Our goal is to extract the concentration measurements from the spectral data using OLS fitting of component spectra. The reference analyte concentrations are accurately known, and the spectral noise of the component spectra is minimal. Therefore, we can demonstrate both accurate extractions of concentrations in the prediction set and accurate assessment of their uncertainties via our formalism. In the second experiment, we record spectra of human artery tissue from which we obtain diagnostic parameters. Although the spectral components of the artery tissue spectra are known with minimal uncertainty, accurate reference concentrations are unavailable. Therefore, we again utilize OLS, but we obtain the relative (normalized) fit coefficients for each model parameter. From this, we demonstrate uncertainty assessment of measurements from biological tissue and the resulting confidence in a particular diagnosis.

3.1.

Concentration Measurements

Raman spectra were acquired from 60 aqueous solutions of glucose, creatinine, and urea with randomized analyte concentrations ranging from 0to50mM . The solutions were contained in a 1-cm fused silica cuvette that had been photobleached for one hour to deplete fluorescent impurities prior to the start of the experiment. The Raman system consisted of an 830-nm diode laser that was directed through a holographic bandpass filter (Kaiser Optical Systems, Incorporated, Ann Arbor, MI) and aperture to reduce emission outside the center wavelength. An external photodiode monitored the intensity of the laser beam and was used to correct for intensity variations. The laser beam was then passed through beam shaping optics and focused into the cuvette through a small hole in a gold-coated paraboloidal mirror (Perkin Elmer, Waltham, MA). The power at the sample was 217mW with a spot area of 1mm2 . Backscattered Raman light was collected by the parabolodial mirror and passed through a 2.5-in. notch filter (Kaiser) to reject the Rayleigh peak at 830nm . The filtered light was focused into an optical fiber bundle composed of 65 fibers, core diameter 396μm , NA=0.37 (Romack Fiber Optics, Williamsburg, VA). The input end was in the form of a circle, and the output end was a single row of 65 fibers, serving as the entrance slit of an f/1.4 spectrometer (Kaiser). The light was dispersed with a holographic grating onto a liquid nitrogen-cooled CCD detector ( 1300×1340b , Princeton Instruments, Trenton, NJ). The integration time per spectrum, 2s , constituted one “frame,” and 30 consecutive frames were collected. Spectra from 280to1700cm1 , occupying 1000 CCD pixels, were used in all data analysis. Owing to the large CCD size and the high-NA imaging system, the entrance slit image appeared curved on the CCD. Direct binning of vertical pixels would result in highly degraded spectral resolution.29 To correct the image curvature, a processing routine was developed that utilizes multiple spectral lines of a strong Raman-active material, such as acetaminophen for curvature calibration. The algorithm preserves instrumental diffraction-limited spectral resolution and improves wavelength accuracy of the measured spectra. The constituent spectra of the three chemically active species (glucose, creatinine, and urea) acquired at 53-mM concentration, as well as those of water and cuvette, are shown in Fig. 1 . By applying the OLS fitting specified by Eq. 5, the experimentally measured, offset-corrected total spectrum from each frame can be decomposed into the concentrations of the underlying constituents. [Alternatively, if the spectral noise is not white, Eq. 3 can be used in place of Eq. 5.] All of the spectral fitting was performed in the wavelength (CCD pixel) domain.

Fig. 1

Constituent Raman spectra—glucose, creatinine, urea, water and cuvette—plotted as functions of wavelength (CCD pixel). The corresponding wavenumber scale is indicated at the botttom.

064012_1_024706jbo1.jpg

Some of the raw experimental spectra contained a distortion in the middle of the spectral range. This artifact, which varied in size from frame to frame, is attributed to variations in the opening and closing of the mechanical shutter that gates the CCD camera, which allows relatively more (or less) light to be collected in the middle of the spectral range. Although the amplitude of this artifact was not very large in absolute terms (50 to 100 counts out of 7000 ), the changes in the predicted concentrations from these faulty frames were significant, creating statistical outliers from the mean of the 30 repeated measurements. The faulty frames were easily identified by looking at the shape of the residuals between the data and the fit, and were thus excluded from the ensuing analysis by setting a threshold on the amplitude of the residual. This resulted in the removal of 440 frames out of the original 1800, so that each sample contained a set of approximately 25 measurements with minimal experimental artifacts.

3.2.

Measurement of Diagnostic Parameters

We have applied the Δc analysis to experimental Raman spectra obtained from human artery tissue to illustrate the application of this analysis method to disease diagnosis. The experiment with excised human carotid artery tissue was part of a separate study and is described in detail elsewhere.30 The spectra were acquired using a clinical Raman system31 and Raman spectral probe.32 The excitation wavelength was 830nm , laser power was 100mW , spot area 1mm2 , and the collection time was 5s , typically acquired in 20 consecutive measurements of 0.25s each. The details of the system are described in Ref. 31. Raman spectra were extracted from the raw spectra by performing a white light correction, removing probe-related background, and subtracting tissue fluorescence.31, 32 A model composed of Raman-active tissue constituents, obtained from confocal Raman microscopy spectra of eight artery morphological structures,12 was used to fit the data using OLS. Prior to fitting, the Raman tissue spectra were interpolated and binned onto the same wavenumber scale as the spectral model constituent spectra. Only relative intensities of the Raman spectral components were obtained, and those relative fit coefficients from the eight spectral components were normalized to sum to unity.5

4.

Results

4.1.

Concentration Measurements

Figure 2 shows the data, the least squares fit using the spectral components, and the difference between the data and the fit (the residual) for one representative mixture. We can analyze how close the predicted parameters ĉ are to the reference values c by means of the plot of Fig. 3 . The root mean squared errors of prediction (RMSEP) across the 1360 total repeated measurements for glucose, creatinine, and urea are 0.488, 0.270, and 0.321mM , respectively.

Fig. 2

Representative data spectrum (blue), the least squares fit (red), and the residual between the data and the fit (black) obtained from a mixture solution. Spectral fitting is performed in the wavelength (CCD pixel) domain; the corresponding wavenumber scale is indicated at the bottom.

064012_1_024706jbo2.jpg

Fig. 3

Predicted concentrations using Eq. 5 versus the reference concentrations for the three analytes. The predicted concentrations closely follow the reference concentrations.

064012_1_024706jbo3.jpg

We next turn to uncertainty analysis. An empirical method to calculate the uncertainty associated with each of the extracted fit coefficients is to repeat the measurement many times, extract the parameters from each individual measurement, and then calculate the standard deviation across the entire set. We refer to this as the measured uncertainty. A faster and more broadly applicable approach to estimate the underlying parameter uncertainty is to employ the analytical formula for Δc , Eq. 8. The parameter values calculated for our spectra, the noise (σ) , signal (sk) , and overlap factor (olfk) , are given in Table 1 .

To analyze how accurately Eq. 8 characterizes the true measured uncertainty, Fig. 4 plots the measured uncertainty calculated from the set of repeated measurements for all 60 mixtures versus Δc . The Δc value for each mixture in Fig. 4 was evaluated using an effective σ equal to the root mean squared value of the individual σ ’s that were calculated from each residual in the set of repeated measurements. The figure also indicates a 45-deg line (black) to reference where Δc equals the measured uncertainty, as well as two additional lines (dotted red) to indicate the region where the measured uncertainty is within a factor of 1.5 of Δc . Note that now, every estimated analyte concentration (Fig. 3) can be associated with an error bar using the Δc uncertainty (Fig. 4).

Fig. 4

Measured uncertainty (standard deviation from repeated measurements) versus uncertainty calculated by the analytic formula Δc for the three analytes. The dotted lines indicate the region for which the measured uncertainty is within a factor of 1.5 of Δc .

064012_1_024706jbo4.jpg

Fig. 5

(a) Representative experimental Raman spectrum (blue), the least squares fit (red), and the residual between the data and the fit (black), obtained from a calcified carotid artery plaque in 0.25s . (b) Diagnostic algorithm, showing several representative Raman artery spectra including the spectrum from (a) (see text for details). The error bars in the two dimensions are calculated using the Δc equation. ( CP=calcified plaque, NCP=noncalcified plaque, IF=intimal fibroplasia.)

064012_1_024706jbo5.jpg

Table 1

Necessary parameters for calculating uncertainty using Eq. 8. The constituent glucose, creatinine, and urea spectra were measured at 53-mM concentration. The only value that varies from sample to sample is σ . The value of σ in the representative spectrum of Fig. 2 is 14.9.

sk (×103) olfk
Cuvette36.11.61
Water59.82.44
Glucose3.581.51
Creatinine4.671.42
Urea3.321.17

4.2.

Measurement of Diagnostic Parameters

Figure 5a shows a representative experimental Raman spectrum, the least squares fit using the artery morphological model, and the residual difference between the data and the fit obtained from a specimen of calcified carotid artery plaque with 0.25-s integration time. The measured uncertainty in the model fit coefficients was calculated by taking the standard deviation from the set of fit coefficients extracted from 20 consecutive measurements of 0.25-s each. This measured uncertainty was compared to the average uncertainty calculated by the Δc formula from Eq. 8 applied to any single one of the 20 independent measurements. For the representative spectrum of Fig. 5a, the value of σ was 0.041, the sk ranged from 2.78 to 8.89, and olfk ranged from 1.05 to 5.83 for the set of spectral components. The measured and Δc uncertainties for this particular sample never deviated from each other by more than a factor of 2.

A previously developed diagnostic algorithm33 was applied to the fit coefficients extracted from the spectrum plotted in Fig. 5a. The diagnostic algorithm uses the fit coefficients from three morphological components (calcium mineralization, cholesterol crystals, and foam cells/necrotic core) to classify the artery sample as being nonatherosclerotic (intimal fibroplasia), noncalcified plaque, or calcified plaque. Because only the relative intensities of the constituent Raman spectra were employed, the raw fit coefficients were normalized so that they sum to unity. In this way, the normalized fit coefficients represent the relative contributions of each morphological feature in the observed spectrum.33

The calculated uncertainties of each raw fit coefficient were similarly scaled to provide uncertainties of the normalized fit coefficients: if fnorm=a*fraw , then the uncertainty propagates as Δfnorm=a*Δfraw . Figure 5b shows the diagnosis for the artery specimen with the spectrum given in Fig. 5a along with diagnoses based on several other spectra from the artery dataset, using the diagnostic space described earlier. The error bars in the two directions indicate the uncertainty (one standard deviation) of the normalized diagnostic fit coefficients.

5.

Discussion

5.1.

Concentration Measurements

We first note that the values of sk and olfk presented in Table 1 make physical sense. Of the three analytes of interest, creatinine has the greatest value of screatinine , which indicates that it has a relatively larger Raman scattering cross section compared to the other analytes. Considering the overlap factors, we note that urea has the smallest olfurea , which can be understood qualitatively by the fact that its constituent spectrum does not overlap strongly with the other constituent spectra (Fig. 1).

Turning to the data presented in Fig. 3, we note that the relatively small values of RMSEP given the large range of concentrations means that the estimator is unbiased, as expected. In addition, we find that the RMSEP values are not very different (<4%) if the data are fit using Eq. 3, which takes into account wavelength-dependent noise variations, rather than Eq. 5. In making this comparison, Eq. 3 is evaluated using a noise covariance Cw matrix, whose diagonal elements are the calculated wavelength-dependent variances of each pixel across the set of repeated measurements, while the off-diagonal covariance terms are set to zero. [Because the number of repeated measurements (25) was much smaller than the number of points in the spectra (1000) , a direct calculation of the covariance matrix results in a rank-deficient matrix. Instead of attempting to correct this by artificially boosting the diagonal elements, we found it more sensible to just use the individual variances on the diagonals and constrain the off-diagonal covariance terms to be zero, since we know that our sensors are independent in any case.]

Figure 4 indicates that Eq. 8 provides an excellent estimate of the measured uncertainty. The measured uncertainties all lie within a factor of 1.5 of the Δc values calculated by Eq. 8, as indicated by the dotted red lines in Fig. 4, with an average deviation of only 11%. The largest contribution to the vertical spread is due to the fact that the measured uncertainty is calculated across a limited number of measurements (25) , and is thus subject to its own uncertainty. We calculated the standard deviation of the estimate and found it to be from 10 to 15% of the recorded measured uncertainty. Much of the horizontal spread is due to the fact that the value of σ used in Eq. 8 is also calculated and is thus an estimate of the true σ . This estimate is related to the discussion of actual versus limiting uncertainty. Any remaining deviations are likely due to subtle uncontrollable experimental factors.

We can quantify how close we are to the limiting uncertainty by comparing values of σ calculated by Eqs. 9, 10. Across the 60 samples, the average value of σact as calculated by Eq. 10 is about 4% higher than the value of σlim , as calculated by Eq. 9. Although σact is greater than σlim , the difference is very small, indicating that this measurement is very close to the limiting uncertainty. This good agreement is evidence that the linear model is valid and that measurement noise is the dominant source of uncertainty.

These results indicate that the analytical uncertainty analysis framework is an accurate and useful way of characterizing the experimental uncertainty obtainable from a single measurement. We note that great care was taken in accurately measuring the reference concentrations and minimizing spectral noise in the model, as well as in collecting the prediction spectra, thus fulfilling the necessary conditions of Eq. 1. Because the data preprocessing steps such as curvature correction are linear and deterministic, the assumption of uncorrelated noise holds as well for the corrected spectra as for the raw spectra. However, we observed that the distributions are only approximately Gaussian, and this may explain some of the small deviation. Other sources that account for imperfect agreement between σact and σlim include the P matrix being imperfect in modeling mixtures due to chemical interactions in the solution and perhaps other minor inaccuracies in the model.

The agreement presented in Fig. 4 is similar to that which would result had Eq. 4 been used to calculate the uncertainties while using Eq. 3 to predict the concentrations (data not shown), using Cw as described earlier. This result, together with consistent prediction accuracy described before, underscores the validity of the initial assumption that, when noise is relatively constant across pixels, Eqs. 5, 8 are valid practical approximations to the more analytically accurate Eqs. 3, 4 for the purposes of estimating parameters and their uncertainties, respectively. In this regime, use of Eqs. 5, 8 is advantageous, as it can be applied to a single spectrum, rather than requiring multiple repeated measurements to obtain Cw .

5.2.

Measurement of Diagnostic Parameters

The Δc analysis is particularly useful for calculating uncertainty in the parameters extracted from artery tissue. This uncertainty translates into diagnostic error bars [Fig. 5b] that indicate the confidence of the overall diagnosis. Note that in Fig. 5b, one of the diagnostic dimensions is the sum of two normalized fit coefficients: cholesterol crystals and foam cells/necrotic core. The uncertainty of the sum involves the uncertainties (variances) of each individual fit coefficients as well as the covariance of the two, which is specified by the off-diagonal terms of the matrix in Eq. 6.

Consider two specific specimens located on opposite sides of the decision line between calcified and noncalcified plaque, represented by open and solid squares in Fig. 5b. Without knowledge of the uncertainty in these assignments, one cannot be more or less confident of either classification assignment. Δc analysis allows for both qualitative and quantitative assessment of the confidence of the diagnosis by assignment of error bars, which effectively specify a probability distribution. If we assume a bivariate Gaussian distribution specified by mean corresponding to the fit coefficients and covariance matrix calculated by Eq. 6, we calculate the probability of the solid square specimen being calcified as 80%. Similarly, we calculate the probability of the open square specimen being noncalcified as 60%.

For simplicity, here we have considered that the classification algorithm is perfect, meaning there is absolute certainty about the decision, regardless of proximity to the decision line. However, in practice there is an additional probability associated with classification that arises from an imperfect decision line. Therefore, a more rigorous approach would be to use this classifier probability as a weighting factor on the data point probability in calculating diagnostic confidence.

We note that the value of σact is on average 15% higher than σlim for the 17 artery specimens examined, a somewhat larger discrepancy than that observed for the concentration measurements. Although data preprocessing steps are linear and deterministic, thus preserving the assumption about uncorrelated noise, the tissue Raman spectra are fit after interpolation and binning, which could undermine that assumption. Therefore, we calculated the changes in σact and σlim before and after interpolation, but found only a small difference (<5%) between the two. We attribute the remaining differences between σact and σlim to minor structure in the residuals that result in artificially high calculated values of σact . This finding is not surprising when considering the complex nature of tissue modeling, and indicates that there is room for improvement in the modeling.

5.3.

General Comments

Knowledge of the limiting uncertainty also provides the limit of detection. For example, under the stated assumptions, if we calculate a value of Δc for a particular sample, we can be reasonably certain of detecting parameters (concentrations or fit coefficients) on the order of 3Δc . This quantity specifies the lowest concentration of an analyte such as glucose that can be detected in a mixture solution or, equivalently, the smallest contribution of a morphological pure component from a tissue sample.

The differences between the actual and limiting uncertainties can be broken down into three cases. In the first case, σactσlim and the residuals are featureless. This implies that ΔcactΔclim , indicating that the measurements are being made with minimal uncertainty. Given the small difference between σact and σlim for the concentration measurements presented, we can conclude that these measurements fall in this category. The second case is that σact> σlim and the residuals are near featureless. This implies that there is noise in the model components that can be further reduced. This case would hold true in applications where the constituent spectra were not measured directly but rather obtained through PCR, for example; thus, the spectral components (principal components) may be noisy and add to measurement noise. Even for direct measurements of the spectral components, a particular component spectrum may contain more noise than others and may need to be collected again. Most of the artery tissue measurements fall in this category, and the uncertainty analysis should guide improvements in modeling until the limiting uncertainty is reached. Lastly, the third case, where σact> σlim and the residuals have structure, means that there are model components missing or there is some other error in the preprocessing of the data. Some infrequent tissue spectra fall in this category as tissue is very heterogeneous, especially when analyzing disease progression. In this case, careful understanding of the sample properties and variation, as well as accurate modeling, is needed to bring the uncertainty down to the limiting level.

As mentioned earlier, the limiting uncertainty can also be expressed in terms of parameters extracted by indirect calibration (such as PLS), that depend on measurement noise and the b-vector.16 The present work provides a natural extension to those results by demonstrating the applicability of a more general formula, Eq. 8, that arises naturally from the CRLB concept and that effectively breaks up the b-vector from Ref. 16 into the signal strength and spectral overlap contributions. When only indirect calibration is possible, such as for concentration measurements in solutions where individual spectral contributions of the constituents cannot be measured directly, the formula from Ref. 16 should be used. When both direct and indirect calibration are possible, both the formula from Ref. 16 and Eq. 8 can be used; this way, the formulas can be utilized to test and compare the robustness of indirect (PLS) and direct (OLS) prediction methods. Lastly, when indirect calibration is not possible but the spectral models are measurable, such as for extracting model parameters from samples including human tissue where reference contributions of particular morphological features are almost impossible to obtain, either formula can be used to calculate parameter uncertainty. In fact, the two equations provide exactly the same ultimate mathematical result. However, the advantage of using Eq. 8 is that it provides physical insight into the signal strength and spectral overlap effects on the b-vector. As described earlier, the χ2 approach can always be used to estimate the uncertainty by doing several constrained fits; however, the method is statistical and does not have a functional dependence, hence it cannot provide insight to the nature of the uncertainty.

The demonstration of experimental uncertainty being very close to Eq. 8 for the concentration measurements indicates the precision of our experimental apparatus, and can be used to guide instrument improvements. For example, increasing the slit width of a spectrograph increases the overlap factor by blurring the Raman peaks, but also increases the signal strength. This tradeoff should guide instrument optimization to yield the minimum extracted parameter uncertainty. Such improvements are crucial in developing multimodal spectroscopy systems that require real-time error assessments of parameters extracted from multiple spectral modalities.30

6.

Conclusion

We describe a simple and direct method for calculating the uncertainty from a single spectroscopic measurement and demonstrate its experimental usefulness, both for solution mixtures and human tissue. Not only does the analytic Δc expression, Eq. 8, provide a means of calculating parameter uncertainties, but it also assesses the calibration and consistency of the experimental apparatus. Because the expression from Eq. 8 is the CRLB, it represents the ultimate lower bound on the uncertainty of parameters extracted from a linear system by an unbiased estimator. Analytical expressions for characterizing uncertainty for nonlinear fitting, such as modeling diffuse reflectance spectroscopy measurements, are also presently under investigation.

Appendices

Appendix

The MVU estimator is equivalent to the ML estimator for the linear model described earlier.23 In this section, we demonstrate that the standard deviation of the estimator obtained analytically is the same as that obtained through χ2 analysis. The χ2 value is defined as:

Eq. 11

χ2=i=1N(dataifitiσi)2,
and the associated ML estimator is determined by minimizing this value:

Eq. 12

ĉ=argminĉ(χ2).
Specializing to our linear model and recognizing the fit as Pĉ , we can express χ2 in the form of an inner product:

Eq. 13

χ2=(sPĉ)T(sPĉ)σ2.
The second derivative of χ2 with respect to our estimate ĉ has a particularly useful form. Noting that:

Eq. 14

d2dĉ2χ2=d2dĉ2(Pĉ)T(Pĉ)σ2=d2dĉ2ĉT(PTP)ĉσ2,
and using the general property (d2dx2)xTAx=(A+AT) and Eq. 6, we can simplify the derivative quantity to:

Eq. 15

d2dĉ2χ2=2(PTP)σ2=2cov(ĉ)1.
Note that (d2dĉ2)χ2 is a Hessian matrix whose (i,j) entry specifies (d2dĉidĉj)χ2 . This matrix is used by standard optimization techniques and, as we see later, is particularly useful when evaluated at the minimum χ2 value. Dividing by 2, inverting, specializing to the k ’th diagonal component, and taking the square root, we obtain:

Eq. 16

[2(d2dĉ2χ2)(k,k)1]12=std(ĉk)Δc.
This expression demonstrates the connection between the analytical and standard χ2 analysis of error. That is, the curvature of χ2 (as function of ĉk ) is inversely proportional to the variance of the k ’th estimator ĉk .

If we had not made the simplifying assumption that Cw=σ2I for the noise vector w in Eq. 1, we could still follow an analogous mathematical approach, as shown before, to demonstrate that the curvature of χ2 specifies the covariance, as given by the more general formula of Eq. 4. In this case, Eq. 15 would become:

Eq. 17

d22dĉ2χ2=PTC1P=cov(ĉ)1,
and the final result of Eq. 16 still follows.

Acknowledgments

This research was conducted at the MIT Laser Biomedical Research Center under NIH grant number P41-RR-02594. The authors thank the reviewers for insightful comments. We also thank Sourav R. Dey for fruitful discussions on topics of estimation theory.

References

1. 

I. Gueorguieva, L. Aarons, K. Ogungbenro, K. M. Jorga, T. Rodgers, and M. Rowland, “Optimal design for multivariate response pharmacokinetic models,” J. Pharmacokinetics Pharmacodynam., 33 (2), 97 –124 (2006). Google Scholar

2. 

M. R. Riley, M. A. Arnold, D. W. Murhammer, E. L. Walls, and N. DelaCruz, “Adaptive calibration scheme for quantification of nutrients and byproducts in insect cell bioreactors by near-infrared spectroscopy,” Biotechnol. Prog., 14 (3), 527 –533 (1998). 8756-7938 Google Scholar

3. 

D. A. Benaron, “The future of cancer imaging,” Cancer Metastasis Rev., 21 (1), 45 –78 (2002). 0891-9992 Google Scholar

4. 

A. S. Haka, Z. Volynskaya, J. A. Gardecki, J. Nazemi, J. Lyons, D. Hicks, M. Fitzmaurice, R. R. Dasari, J. P. Crowe, and M. S. Feld, “In vivo margin assessment during partial mastectomy breast surgery using Raman spectroscopy,” Cancer Res., 66 (6), 3317 –3322 (2006). https://doi.org/10.1158/0008-5472.CAN-05-2815 0008-5472 Google Scholar

5. 

J. T. Motz, M. Fitzmaurice, A. Miller, S. J. Gandhi, A. S. Haka, L. H. Galindo, R. R. Dasari, J. R. Kramer, and M. S. Feld, “In vivo Raman spectral pathology of human atherosclerosis and vulnerable plaque,” J. Biomed. Opt., 11 (2), 021003 (2006). https://doi.org/10.1117/1.2190967 1083-3668 Google Scholar

6. 

I. J. Bigio, S. G. Bown, G. Briggs, C. Kelley, S. Lakhani, D. Pickard, P. M. Ripley, I. G. Rose, and C. Saunders, “Diagnosis of breast cancer using elastic-scattering spectroscopy: preliminary clinical results,” J. Biomed. Opt., 5 (2), 221 –228 (2000). https://doi.org/10.1117/1.429990 1083-3668 Google Scholar

7. 

G. M. Palmer, C. F. Zhu, T. M. Breslin, F. S. Xu, K. W. Gilchrist, and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part II: Application to breast cancer diagnosis,” Appl. Opt., 45 (5), 1072 –1078 (2006). https://doi.org/10.1364/AO.45.001072 0003-6935 Google Scholar

8. 

H. Kaiser, “Foundations for the critical discussion of analytical methods,” Spectrochim. Acta, Part B, 33 (9), 551 –576 (1978). 0584-8547 Google Scholar

9. 

A. C. Olivieri, N. K. M. Faber, J. Ferre, R. Boque, J. H. Kalivas, and H. Mark, “Uncertainty estimation and figures of merit for multivariate calibration,” Pure Appl. Chem., 78 (3), 633 –661 (2006). 0033-4545 Google Scholar

10. 

H. Martens, Multivariate Calibration, John Wiley and Sons, New York (1989). Google Scholar

11. 

K. E. Shafer-Peltier, A. S. Haka, M. Fitzmaurice, J. Crowe, J. Myles, R. R. Dasari, and M. S. Feld, “Raman microspectroscopic model of human breast tissue: implications for breast cancer diagnosis in vivo,” J. Raman Spectrosc., 33 (7), 552 –563 (2002). https://doi.org/10.1002/jrs.877 0377-0486 Google Scholar

12. 

H. P. Buschman, G. Deinum, J. T. Motz, M. Fitzmaurice, J. R. Kramer, A. van der Laarse, A. V. Bruschke, and M. S. Feld, “Raman microspectroscopy of human coronary atherosclerosis: biochemical assessment of cellular and extracellular morphologic structures in situ,” Cardiovasc. Pathol., 10 (2), 69 –82 (2001). https://doi.org/10.1016/S1054-8807(01)00064-3 1054-8807 Google Scholar

13. 

A. Lorber, “Error propagation and figures of merit for quantification by solving matrix equations,” Anal. Chem., 58 (6), 1167 –1172 (1986). https://doi.org/10.1021/ac00297a042 0003-2700 Google Scholar

14. 

A. Lorber, K. Faber, and B. R. Kowalski, “Net analyte signal calculation in multivariate calibration,” Anal. Chem., 69 (8), 1620 –1626 (1997). https://doi.org/10.1021/ac960862b 0003-2700 Google Scholar

15. 

P. R. Bevington and D. K. Robinson, Data Reduction and Error Analysis for the Physical Sciences, McGraw-Hill, Boston (2003). Google Scholar

16. 

A. J. Berger and M. S. Feld, “Analytical method of estimating chemometric prediction error,” Appl. Spectrosc., 51 (5), 725 –732 (1997). 0003-7028 Google Scholar

17. 

T.-W. Koo, “Measurement of blood analytes in turbid biological tissue using near-infrared Raman spectroscopy,” 261 Massachusetts Institute of Technology, (2001). Google Scholar

18. 

H. L. T. Lee, P. Boccazzi, N. Gorret, R. J. Ram, and A. J. Sinskey, “In situ bioprocess monitoring of Escherichia coli bioreactions using Raman spectroscopy,” Vib. Spectrosc., 35 (1–2), 131 –137 (2004). 0924-2031 Google Scholar

19. 

A. Lorber and B. R. Kowalski, “Estimation of prediction error for multivariate calibration,” J. Chemom., 2 (67), 93 –109 (1988). https://doi.org/10.1002/cem.1180020203 0886-9383 Google Scholar

20. 

K. Faber and B. R. Kowalski, “Propagation of measurement errors for the validation of predictions obtained by principal component regression and partial least squares,” J. Chemom., 11 (3), 181 –238 (1997). 0886-9383 Google Scholar

21. 

D. M. Haaland and D. K. Melgaard, “New prediction-augmented classical least-squares (PACLS) methods: Application to unmodeled interferents,” Appl. Spectrosc., 54 (9), 1303 –1312 (2000). https://doi.org/10.1366/0003702001951228 0003-7028 Google Scholar

22. 

B. Nadler and R. R. Coifman, “The prediction error in CLS and PLS: the importance of feature selection prior to multivariate calibration,” J. Chemom., 19 (2), 107 –118 (2005). https://doi.org/10.1002/cem.915 0886-9383 Google Scholar

23. 

S. M. Kay, Fundamentals of Statistical Signal Processing, PTR Prentice-Hall, Englewood Cliffs, NJ (1993). Google Scholar

24. 

R. Manoharan, J. J. Baraga, M. S. Feld, and R. P. Rava, “Quantitative histochemical analysis of human artery using Raman-spectroscopy,” J. Photochem. Photobiol., B, 16 (2), 211 –233 (1992). https://doi.org/10.1016/1011-1344(92)80009-K 1011-1344 Google Scholar

25. 

J. F. Brennan, T. J. Romer, R. S. Lees, A. M. Tercyak, J. R. Kramer, and M. S. Feld, “Determination of human coronary artery composition by Raman spectroscopy,” Circulation, 96 (1), 99 –105 (1997). 0009-7322 Google Scholar

26. 

D. W. Marquardt, “Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation,” Technometrics, 12 (3), 591 –612 (1970). https://doi.org/10.2307/1267205 0040-1706 Google Scholar

27. 

J. H. Kalivas, “Variance-decomposition of pure-component spectra as a measure of selectivity,” J. Chemom., 3 409 –418 (1989). 0886-9383 Google Scholar

28. 

R. A. Arndt and M. H. MacGregor, “Nucleon-nucleon phase shift analyses by chi-squared minimization,” Methods in Computational Physics, 253 –296 Academic Press, New York (1966). Google Scholar

29. 

J. Zhao, “Image curvature correction and cosmic removal for high-throughput dispersive Raman spectroscopy,” Appl. Spectrosc., 57 (11), 1368 –1375 (2003). https://doi.org/10.1366/000370203322554527 0003-7028 Google Scholar

30. 

O. R. Scepanovic, M. Fitzmaurice, J. A. Gardecki, G. O. Angheloiu, S. Awasthi, J. T. Motz, J. R. Kramer, R. R. Dasari, and M. S. Feld, “Detection of morphological markers of vulnerable atherosclerotic plaque using multimodal spectroscopy,” J. Biomed. Opt., 11 (2), 021007 (2006). https://doi.org/10.1117/1.2187943 1083-3668 Google Scholar

31. 

J. T. Motz, S. J. Gandhi, O. R. Scepanovic, A. S. Haka, J. R. Kramer, R. R. Dasari, and M. S. Feld, “Real-time Raman system for in vivo disease diagnosis,” J. Biomed. Opt., 10 (3), 031113 (2005). https://doi.org/10.1117/1.1920247 1083-3668 Google Scholar

32. 

J. T. Motz, M. Hunter, L. H. Galindo, J. A. Gardecki, J. R. Kramer, R. R. Dasari, and M. S. Feld, “Optical fiber probe for biomedical Raman spectroscopy,” Appl. Opt., 43 (3), 542 –554 (2004). https://doi.org/10.1364/AO.43.000542 0003-6935 Google Scholar

33. 

H. P. Buschman, J. T. Motz, G. Deinum, T. J. Romer, M. Fitzmaurice, J. R. Kramer, A. van der Laarse, A. V. Bruschke, and M. S. Feld, “Diagnosis of human coronary atherosclerosis by morphology-based Raman spectroscopy,” Cardiovasc. Pathol., 10 (2), 59 –68 (2001). https://doi.org/10.1016/S1054-8807(01)00063-1 1054-8807 Google Scholar
©(2007) Society of Photo-Optical Instrumentation Engineers (SPIE)
Obrad R. Scepanovic, Kate L. Bechtel, Abigail S. Haka, Wei-Chuan Shih, Tae-Woong Koo, Andrew J. Berger, and Michael S. Feld "Determination of uncertainty in parameters extracted from single spectroscopic measurements," Journal of Biomedical Optics 12(6), 064012 (1 November 2007). https://doi.org/10.1117/1.2815692
Published: 1 November 2007
JOURNAL ARTICLE
10 PAGES


SHARE
Advertisement
Advertisement
Back to Top