## 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 laboratory^{4, 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
$\left({\chi}^{2}\right)$
analysis to calculate parameter uncertainties extracted from a single spectrum.^{15}
${\chi}^{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 ${\chi}^{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
$(\text{row}\times \text{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:

The vector $\mathbf{s}$ is the observed measurement $(M\times 1)$ , the matrix $\mathbf{P}$ contains the model constituent vectors $(M\times N)$ and is full rank, the vector $\mathbf{c}$ contains the underlying coefficients of the model constituents $(N\times 1)$ , and the vector $\mathbf{w}$ represents noise in the system $(M\times 1)$ . Here $M$ represents the number of wavelengths and $N$ represents the number of model components. In other words, the measured vector $\mathbf{s}$ is a linear combination of the model components in $\mathbf{P}$ , weighted by the coefficients in $\mathbf{c}$ , and with the addition of random measurement noise $\mathbf{w}$ . The $\mathbf{w}$ is assumed to be a zero-mean Gaussian random vector with a known or measurable covariance matrix ${\mathbf{C}}_{\mathbf{w}}$ . Our goal is to analytically determine the best estimate $\widehat{\mathbf{c}}$ (fit coefficients) of the true underlying coefficients $\mathbf{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
$\left(\mathbf{s}\right)$
have been shown experimentally to be linear combinations of the distinct individual spectra of the underlying chemical species that can be measured in isolation
$\left(\mathbf{P}\right)$
and proportional to the concentration of those species
$\left(\mathbf{c}\right)$
.^{24, 25} Measurement noise
$\left(\mathbf{w}\right)$
is typically associated with the measurement system and detector, and its statistics
$\left({\mathbf{C}}_{\mathbf{w}}\right)$
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 Kowalski^{19} that can be specialized for our model as

## Eq. 2

$$\mathbf{s}=(\mathbf{P}+\delta \mathbf{P})\mathbf{c}+\mathbf{w}=\mathbf{P}\mathbf{c}+{\mathbf{w}}^{\prime}.$$## 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
$\mathbf{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

$$\widehat{\mathbf{c}}={\left({\mathbf{P}}^{T}{\mathbf{C}}_{w}^{-1}\mathbf{P}\right)}^{-1}{\mathbf{P}}^{T}{\mathbf{C}}_{w}^{-1}\mathbf{s}.$$^{23}that its covariance matrix is:

## Eq. 4

$$\mathrm{cov}\left(\widehat{\mathbf{c}}\right)={\left({\mathbf{P}}^{T}{\mathbf{C}}_{w}^{-1}\mathbf{P}\right)}^{-1}.$$The result of Eq. 4 can be specialized by assuming that
$\mathbf{w}$
is white (i.e., uncorrelated and identically distributed) Gaussian noise, so that
${\mathbf{C}}_{\mathbf{w}}={\sigma}^{2}\mathbf{I}$
. With this assumption, the MVU estimator and the covariance are given by:^{23}

## Eq. 5

$$\widehat{\mathbf{c}}={\left({\mathbf{P}}^{T}\mathbf{P}\right)}^{-1}{\mathbf{P}}^{T}\mathbf{s}.$$## Eq. 6

$$\mathrm{cov}\left(\widehat{\mathbf{c}}\right)={\sigma}^{2}{\left({\mathbf{P}}^{T}\mathbf{P}\right)}^{-1}.$$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.

### $\Delta 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 $\mathbf{P}=\mathbf{Q}\mathbf{S}$ , with $\mathbf{S}$ being a diagonal matrix, where the $k$ ’th diagonal entry ${s}_{k}$ is the Euclidean norm of the $k$ ’th component in $\mathbf{P}$ , as follows

and the columns of the matrix $\mathbf{Q}$ are thus normalized to unit length. This leads to a simple expression for the standard deviation $\Delta c$ of the $k$ ’th estimated parameter ${\widehat{c}}_{k}$ :## Eq. 8

$$\Delta c\equiv \mathrm{std}\left({\widehat{c}}_{k}\right)=\sigma {\left[{\left({\mathbf{P}}^{T}\mathbf{P}\right)}_{(k,k)}^{-1}\right]}^{1\u22152}=\frac{\sigma}{{s}_{k}}{\left[{\left({\mathbf{Q}}^{T}\mathbf{Q}\right)}_{(k,k)}^{-1}\right]}^{1\u22152}=\frac{\sigma}{{s}_{k}}\bullet {\mathit{olf}}_{k}.$$The overlap factor may take on values between 1 and
$\infty $
. If all of the columns of
$\mathbf{P}$
(or equivalently,
$\mathbf{Q}$
) are orthogonal (no overlap), then
${\mathit{olf}}_{k}=1$
. In the other extreme, if the
$k$
’th column of
$\mathbf{P}$
(or
$\mathbf{Q}$
) is nearly linearly dependent with one or more columns (almost complete overlap), then
${\mathbf{Q}}^{T}\mathbf{Q}$
is close to singular and its inverse does not exist. In a generalized sense,
${\mathit{olf}}_{k}$
then approaches
$\infty $
. (In the case of two columns being linearly dependent, one of the two columns should be removed so that
$\mathbf{P}$
becomes full rank, as was specified earlier.) In other words, when the model
$\mathbf{P}$
contains orthogonal constituent spectra, the estimator uncertainty
$\Delta c$
is equal to the ratio
$(\sigma \u2215{s}_{k})$
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
$\Delta c$
approaches
$\infty $
. For the more commonly encountered case of partial spectral overlap of the linearly independent spectral components of
$\mathbf{P}$
, we have
$1<{\mathit{olf}}_{k}<\infty $
, and
$\Delta c$
is a function of the three physically understandable quantities:
$\sigma $
,
${s}_{k}$
, and
${\mathit{olf}}_{k}$
. 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, $\sigma $ varies from sample to sample because shot noise is dependent on the sample-specific raw signal, whereas ${s}_{k}$ and ${\mathit{olf}}_{k}$ are sample-independent for a given spectroscopic technique and model.

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

## Eq. 9

$${\sigma}_{\mathrm{lim}}={\left\{\frac{1}{M}\sum _{i=1}^{M}\left[\frac{1}{L-1}\sum _{j=1}^{L}{({s}_{i,j}-{\overline{s}}_{i})}^{2}\right]\right\}}^{1\u22152},$$## Eq. 10

$${\sigma}_{\mathrm{act}}={\left[\frac{1}{M}{(\mathbf{s}-\mathbf{P}\widehat{\mathbf{c}})}^{T}(\mathbf{s}-\mathbf{P}\widehat{\mathbf{c}})\right]}^{1\u22152},$$Under the assumption of Eq. 1 that the only source of uncertainty is measurement noise, both approaches should yield the same value for $\sigma $ . However, if there are also modeling uncertainties as in Eq. 2, then only the second approach yields the actual $\sigma $ , as the residual includes measurement and modeling noise. Therefore, ${\sigma}_{\mathrm{act}}\u2a7e{\sigma}_{\mathrm{lim}}$ , and the difference between the actual and limiting $\sigma $ can serve to evaluate the accuracy of the modeling. In most of the following, we use Eq. 10 to calculate $\sigma $ . By extension, we can define $\Delta {c}_{\mathrm{act}}$ and $\Delta {c}_{\mathrm{lim}}$ using the values of ${\sigma}_{\mathrm{act}}$ and ${\sigma}_{\mathrm{lim}}$ , respectively.

## 2.5.

### Relation to Chi-Squared $\left({\chi}^{2}\right)$

A statistical method of calculating the uncertainty in extracted parameters can be implemented through
${\chi}^{2}$
analysis. Maximizing the likelihood of observing a particular measured spectrum (in random Gaussian noise) is equivalent to minimizing
${\chi}^{2}$
. The value of a parameter that minimizes
${\chi}^{2}$
is the optimal maximum likelihood (ML) value. The
${\chi}^{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
${\chi}^{2}$
function.^{15} Equivalently, an increase of one standard deviation of the parameter from the value at the minimum increases
${\chi}^{2}$
by unity.^{15, 28} The
${\chi}^{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 $\Delta c$ , Eq. 8, is equivalent to the one obtained statistically through ${\chi}^{2}$ analysis. This is to be expected, since the MVU estimator (or the least squares estimator) is equivalent to the ML estimator. Hence the ${\chi}^{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
$0\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}50\phantom{\rule{0.3em}{0ex}}\mathrm{mM}$
. The solutions were contained in a
$1\text{-}\mathrm{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\text{-}\mathrm{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
$217\phantom{\rule{0.3em}{0ex}}\mathrm{mW}$
with a spot area of
$\sim 1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$
. Backscattered Raman light was collected by the parabolodial mirror and passed through a
$2.5\text{-}\mathrm{in.}$
notch filter (Kaiser) to reject the Rayleigh peak at
$830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. The filtered light was focused into an optical fiber bundle composed of 65 fibers, core diameter
$396\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
,
$\mathrm{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\times 1340\mathrm{b}$
, Princeton Instruments, Trenton, NJ). The integration time per spectrum,
$2\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
, constituted one “frame,” and 30 consecutive frames were collected. Spectra from
$280\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1700\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, 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\text{-}\mathrm{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.

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 $\sim 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
$\Delta 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 system^{31} and Raman spectral probe.^{32} The excitation wavelength was
$830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, laser power was
$100\phantom{\rule{0.3em}{0ex}}\mathrm{mW}$
, spot area
$\sim 1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$
, and the collection time was
$5\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
, typically acquired in 20 consecutive measurements of
$0.25\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
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 $\widehat{\mathbf{c}}$ are to the reference values $\mathbf{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.321\phantom{\rule{0.3em}{0ex}}\mathrm{mM}$ , respectively.

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
$\Delta c$
, Eq. 8. The parameter values calculated for our spectra, the noise
$\left(\sigma \right)$
, signal
$\left({s}_{k}\right)$
, and overlap factor
$\left({\mathit{olf}}_{k}\right)$
, 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 $\Delta c$ . The $\Delta c$ value for each mixture in Fig. 4 was evaluated using an effective $\sigma $ equal to the root mean squared value of the individual $\sigma $ ’s that were calculated from each residual in the set of repeated measurements. The figure also indicates a $45\text{-}\mathrm{deg}$ line (black) to reference where $\Delta 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 $\Delta c$ . Note that now, every estimated analyte concentration (Fig. 3) can be associated with an error bar using the $\Delta c$ uncertainty (Fig. 4).

## 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 | |
---|---|---|

Cuvette | 36.1 | 1.61 |

Water | 59.8 | 2.44 |

Glucose | 3.58 | 1.51 |

Creatinine | 4.67 | 1.42 |

Urea | 3.32 | 1.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\text{-}\mathrm{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\text{-}\mathrm{s}$ each. This measured uncertainty was compared to the average uncertainty calculated by the $\Delta 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 $\sigma $ was 0.041, the ${s}_{k}$ ranged from 2.78 to 8.89, and ${\mathit{olf}}_{k}$ ranged from 1.05 to 5.83 for the set of spectral components. The measured and $\Delta c$ uncertainties for this particular sample never deviated from each other by more than a factor of 2.

A previously developed diagnostic algorithm^{33} 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 ${f}_{\mathrm{norm}}={a}^{*}{f}_{\text{raw}}$ , then the uncertainty propagates as $\Delta {f}_{\mathrm{norm}}={a}^{*}\Delta {f}_{\text{raw}}$ . 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 ${s}_{k}$ and ${\mathit{olf}}_{k}$ presented in Table 1 make physical sense. Of the three analytes of interest, creatinine has the greatest value of ${s}_{\text{creatinine}}$ , 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 ${\mathit{olf}}_{\text{urea}}$ , 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 ${\mathbf{C}}_{\mathbf{w}}$ 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 $(\sim 25)$ was much smaller than the number of points in the spectra $(\sim 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 $\Delta 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 $(\sim 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 $\sigma $ used in Eq. 8 is also calculated and is thus an estimate of the true $\sigma $ . 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 $\sigma $ calculated by Eqs. 9, 10. Across the 60 samples, the average value of ${\sigma}_{\mathrm{act}}$ as calculated by Eq. 10 is about 4% higher than the value of ${\sigma}_{\mathrm{lim}}$ , as calculated by Eq. 9. Although ${\sigma}_{\mathrm{act}}$ is greater than ${\sigma}_{\mathrm{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 ${\sigma}_{\mathrm{act}}$ and ${\sigma}_{\mathrm{lim}}$ include the $\mathbf{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 ${\mathbf{C}}_{\mathbf{w}}$ 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 ${\mathbf{C}}_{\mathbf{w}}$ .

## 5.2.

### Measurement of Diagnostic Parameters

The $\Delta 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. $\Delta 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 ${\sigma}_{\mathrm{act}}$ is on average 15% higher than ${\sigma}_{\mathrm{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 ${\sigma}_{\mathrm{act}}$ and ${\sigma}_{\mathrm{lim}}$ before and after interpolation, but found only a small difference $(<5\%)$ between the two. We attribute the remaining differences between ${\sigma}_{\mathrm{act}}$ and ${\sigma}_{\mathrm{lim}}$ to minor structure in the residuals that result in artificially high calculated values of ${\sigma}_{\mathrm{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 $\Delta c$ for a particular sample, we can be reasonably certain of detecting parameters (concentrations or fit coefficients) on the order of $\u2a7e3\Delta 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, ${\sigma}_{\mathrm{act}}\approx {\sigma}_{\mathrm{lim}}$ and the residuals are featureless. This implies that $\Delta {c}_{\mathrm{act}}\approx \Delta {c}_{\mathrm{lim}}$ , indicating that the measurements are being made with minimal uncertainty. Given the small difference between ${\sigma}_{\mathrm{act}}$ and ${\sigma}_{\mathrm{lim}}$ for the concentration measurements presented, we can conclude that these measurements fall in this category. The second case is that ${\sigma}_{\mathrm{act}}>{\sigma}_{\mathrm{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 ${\sigma}_{\mathrm{act}}>{\sigma}_{\mathrm{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
${\chi}^{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 $\Delta 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
${\chi}^{2}$
analysis. The
${\chi}^{2}$
value is defined as:

## Eq. 11

$${\chi}^{2}=\sum _{i=1}^{N}{\left(\frac{{\text{data}}_{i}-{\text{fit}}_{i}}{{\sigma}_{i}}\right)}^{2},$$## Eq. 12

$$\widehat{\mathbf{c}}=\underset{\widehat{\mathbf{c}}}{\mathrm{arg}\phantom{\rule{0.2em}{0ex}}\mathrm{min}}\left({\chi}^{2}\right).$$## Eq. 13

$${\chi}^{2}=\frac{{(\mathbf{s}-\mathbf{P}\widehat{\mathbf{c}})}^{T}(\mathbf{s}-\mathbf{P}\widehat{\mathbf{c}})}{{\sigma}^{2}}.$$## Eq. 14

$$\frac{{d}^{2}}{d{\widehat{\mathbf{c}}}^{2}}{\chi}^{2}=\frac{{d}^{2}}{d{\widehat{\mathbf{c}}}^{2}}\frac{{\left(\mathbf{P}\widehat{\mathbf{c}}\right)}^{T}\left(\mathbf{P}\widehat{\mathbf{c}}\right)}{{\sigma}^{2}}=\frac{{d}^{2}}{d{\widehat{\mathbf{c}}}^{2}}\frac{{\widehat{\mathbf{c}}}^{T}\left({\mathbf{P}}^{T}\mathbf{P}\right)\widehat{\mathbf{c}}}{{\sigma}^{2}},$$## Eq. 15

$$\frac{{d}^{2}}{d{\widehat{\mathbf{c}}}^{2}}{\chi}^{2}=\frac{2\left({\mathbf{P}}^{T}\mathbf{P}\right)}{{\sigma}^{2}}=2\bullet \mathrm{cov}{\left(\widehat{\mathbf{c}}\right)}^{-1}.$$## Eq. 16

$${\left[2{\left(\frac{{d}^{2}}{d{\widehat{\mathbf{c}}}^{2}}{\chi}^{2}\right)}_{(k,k)}^{-1}\right]}^{1\u22152}=\mathrm{std}\left({\widehat{c}}_{k}\right)\equiv \Delta c.$$If we had not made the simplifying assumption that ${\mathbf{C}}_{\mathbf{w}}={\sigma}^{2}\mathbf{I}$ 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 ${\chi}^{2}$ specifies the covariance, as given by the more general formula of Eq. 4. In this case, Eq. 15 would become:

## Eq. 17

$$\frac{{d}^{2}}{2d{\widehat{\mathbf{c}}}^{2}}{\chi}^{2}={\mathbf{P}}^{T}{\mathbf{C}}^{-1}\mathbf{P}=\mathrm{cov}{\left(\widehat{\mathbf{c}}\right)}^{-1},$$## 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

*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

*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

*in vivo*,” J. Raman Spectrosc., 33 (7), 552 –563 (2002). https://doi.org/10.1002/jrs.877 0377-0486 Google Scholar

*In situ*bioprocess monitoring of Escherichia coli bioreactions using Raman spectroscopy,” Vib. Spectrosc., 35 (1–2), 131 –137 (2004). 0924-2031 Google Scholar

*in vivo*disease diagnosis,” J. Biomed. Opt., 10 (3), 031113 (2005). https://doi.org/10.1117/1.1920247 1083-3668 Google Scholar