1.

Introduction

Fitting in vivo optical spectroscopic data with mathematical models that describe the underlying physical “principles” is common practice in biomedical optics. Examples of such analyses include the fitting of multicomponent fluorescence spectra in tissue (that can include tissue autofluorescence as well as exogenous fluorophores), the fitting of multicomponent Raman spectra to determine the tissue biochemistry, and the fitting of multicomponent extinction spectra of tissue chromophores using reflectance spectroscopic measurements. In vivo fluorescence, Raman, and reflectance spectroscopic measurements are all affected by the tissue optical properties. The effect of the optical properties on the measured spectra must first be accounted for using physical light-tissue interaction models to obtain the “intrinsic” fluorescence, Raman, and extinction spectra. These undistorted spectra can then be analyzed with linear component analysis techniques such as singular value decomposition or other least squares minimization techniques. However, if the physical models used in the “preprocessing step” are inaccurate, an error in the fitted component values will be introduced. Furthermore, even if the algorithms used to extract the intrinsic optical spectra from the measured spectra are completely accurate, a lack of knowledge of either the in vivo spectral shape or even the physical presence of some of the components may introduce another source of error in the analysis of in vivo optical spectroscopic measurements. Finally, even when the intrinsic spectra are determined correctly and full knowledge about the presence and spectral shapes of all the components is available, another source of error (uncertainty) in the analysis of optical spectroscopic measurements is introduced by a limited quality of the data (signal-to-noise ratio) that can prevent accurate determination of some or all of the components in the component analysis. The latter type of error is called a statistical error, while the first two types of errors are considered as systematic errors. These errors will introduce an uncertainty in the fitted parameter values extracted from the component analysis. Currently, it is not common practice in the field of biomedical optics to estimate the accuracy of the fitted parameter values from single spectroscopic measurements. However, we believe that it is useful to at least determine the statistical error in each of the fitted parameter values. For example, in the case of averaging multiple in vivo measurements, calculating the weighted average of the parameters using the statistical errors as weight factors would be more appropriate than simply calculating the unweighted average and ignoring the statistical uncertainties associated with the quality of the fits of the different spectra. Moreover, calculating the statistical errors in the fit parameters facilitates objective assessment of the quality of spectra and fits with the possibility to reject poor quality spectra or fit values from an in vivo dataset.

In principle, it is possible to calculate the parameter uncertainty by repeating the measurements and analyzing the standard deviation of parameters extracted from each of these multiple measurements. However, this is not practical for in vivo applications such as medical diagnosis, in which typically only one or a few measurements can be acquired. Scepanovic ¹ recently published an error analysis for single Raman spectroscopic measurements. They presented an analytical formula for estimating uncertainty expressed as a function of measurement noise, signal strength, and spectral overlap. Their analysis characterizes uncertainty for linear fitting functions such as Raman and fluorescence spectroscopy, but has not been tested for nonlinear fitting functions encountered in reflectance measurements. In this work we validate a general method for determining the statistical error or confidence interval (CI) on fitted parameters derived from modeling single optical reflectance spectroscopic measurements using synthetic datasets. The method estimates the parameter CIs as the square roots of the diagonal elements of the covariance matrix, which is obtained by multiplying the inverse of the second derivative matrix of ${\chi }^2$ with respect to its free parameters by ${\chi }^2∕v$ , with $v$ being the number of degrees of freedom. While the method is in principle valid for all spectroscopic techniques (fluorescence, Raman, and reflectance), we demonstrate the validity of the approach using reflectance spectroscopy as an example. Reflectance spectroscopy is sensitive to the absorption and scattering properties of tissue. The absorption coefficient is proportional to the concentration of the tissue chromophores, and the scattering coefficient reflects the size and density of scattering structures in the tissue. Several methods have been employed to extract physically relevant information (such as blood volume fraction and blood saturation) from reflectance spectra, e.g., using an analytical approximation of the transport equation such as the diffusion equation,^{2, 3, 4, 5, 6} Monte-Carlo-based models,^{7, 8} or empirical methods.^{9, 10, 11} Irrespective of the approach, a model function minimization or fitting routine must always be used to find the optimal values for the free parameters of the models. In general, a least squares $\left({\chi }^2\right)$ optimization routine is used where the (weighted) square of the difference between the data points and the fitted curve is minimized. Most commercial software packages that employ least squares fitting routines automatically generate the covariance (also called variance-covariance or error) matrix. In principle, the diagonal elements of this matrix represent the variances of the parameters if 1. the fitting function is linear in the parameters and 2. the measurement uncertainties of the individual data points are known exactly. However, neither of these points is fulfilled for a single reflectance spectroscopy measurement. Nevertheless, we demonstrate that our simple method to estimate the statistical error of the fit parameters from single reflectance spectroscopic measurements yields correct CIs as long as the model used to describe the data is correct. Furthermore, we analyze the effect of imperfections in the fitting model (systematic errors) and show that even minor imperfections in the model may introduce a bias in the fitted parameters that greatly exceeds the estimated statistical CIs. Finally, we investigate the use of methods to identify and subsequently minimize the bias in the fitted parameters associated with incorrect modeling.

2.

Theory

In this work we use the empirical model used to analyze a particular type of reflectance spectroscopy that we have been investigating over the past $5\text{years}$ , differential path-length spectroscopy (DPS),^{11, 12, 13, 14} as an example. The model to which the differential path-length spectra are fitted can be written as:

With our experimental setup,^{11, 12, 13, 14, 15, 16} the data acquired in each reflectance spectrum consist of $N=2048\text{pixel}$ values $R_{\text{data}}\left(j\right)$ corresponding to the differential reflectance at wavelengths ranging from $350\text{to}1000\mathrm{nm}$ . Since the transmission efficiency of the setup is different for different wavelengths, the signal-to-noise ratio (SNR) is wavelength dependent. Clearly it is desirable to put less weight on pixels with a low SNR in the fitting routine. To accomplish this, the spectra are first smoothed by binning $p$ pixel values and using the $N∕p=n$ average values $R_{\text{data}}\left(i\right)$ and corresponding standard deviations $s\left(i\right)$ in the fitting routine. In our simulations, we have typically used a bin width of $15\text{pixels}$ , which corresponds to an optical resolution of $6\mathrm{nm}$ , but we have also tested the effect of smaller and larger bin widths (10 and $30\text{pixels}$ , respectively). The data are then fitted to Eq. 1 using a least squares Levenberg-Marquardt fitting routine using the standard deviations of the binned data points as weight factors. Chi-square (the quantity that is minimized) is given by

The minimum value of ${\chi }^2$ corresponds to values for the seven fit parameters $(a_1$ to $a_7)$ that best describe the data. In principle, the CI of each of these parameters is given by the square root of the diagonal elements of the covariance matrix. The covariance matrix is the inverse of the second derivative matrix of ${\chi }^2$ with respect to its free parameters and is exact, if and only if; 1. the fitting function is linear in its parameters, and 2. the weight factors used in the minimization routine are exact.^{19, 20} In the case of our nonlinear model Eq. 1 with estimated weight factors $s\left(i\right)$ , we estimate the covariance matrix by multiplying the inverse of the second derivative matrix of ${\chi }^2$ with respect to its free parameters by ${\chi }^2∕v$ , with $v$ being the number of degrees of freedom.²⁰ The CI of each of the parameters is given by the square root of the diagonal elements of the estimated covariance matrix under the assumption that the model [Eq. 1] is correct. We validate that the estimated covariance matrix is reasonably accurate for our nonlinear fitting problem in Eq. 1.

The calculated CI represents the statistical error only and does not fully account for systematic errors due to inaccuracies in the physical model used to fit the data. For example, if an additional absorber was present in the sampling volume but not in the model fit function, or if the wavelength dependence of the scattering function is not modeled correctly, the fitted parameters will be biased with respect to their true values. We investigate the magnitude of the bias in the fitted parameters for the case of 1. a missing absorber in the model and 2. an error in the modeled scattering function, and investigate the use of methods to identify and subsequently minimize the bias in the fitted parameters associated with incorrect modeling.

3.

Simulations

We have generated 315 synthetic datasets, each consisting of $2048\text{pixel}$ reflectance values $R_{\mathrm{syn}}\left(j\right)$ using Eqs. 1, 2 with the parameters as described in Table 1 .

Table 1

Parameter values used for the synthetic datasets.

For each combination of parameter values, 100 different fits are performed by generating noise on the synthetic spectrum. The noise was generated by adding a random number between $-1$ and 1, multiplied by a specified noise amplitude to the pixels: $R_{\text{data}}\left(j\right)=R_{\mathrm{syn}}\left(j\right)+\mathrm{noiseamp}\cdot \text{random}[-1,1]$ . Note that the noise generated in this way is independent of wavelength. Since the data $R_{\mathrm{syn}}\left(j\right)$ are wavelength dependent, the relative noise amplitude is in fact wavelength dependent, which is typically the case for optical spectroscopic measurements. For our experimental setup, typical signal-to-noise ratios at ${\lambda }_0=800\mathrm{nm}$ are better than 100, and we have used $\mathrm{noiseamp}=0.01$ for most of the simulations, but also checked the validity of our approach for larger noise amplitudes (0.1 and 0.5). By definition, if the calculated confidence intervals are correct, the fitted parameter values should be within 2 CIs in 95 out of the 100 fits.

4.

Results

4.1.

Fitting with the Correct Model

We first analyze whether the CIs calculated by taking the square root of the diagonal elements of the estimated covariance matrix are correct, as a function of bin width and noise amplitude. The fits are performed using the correct model described in Eqs. 1, 2 with seven free parameters. Furthermore, we calculate the magnitude of the CIs as a function of blood volume fraction for this case.

4.1.1.

Influence of bin width

Table 2 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for the 315 synthetic datasets (standard deviations are indicated in parentheses) for different bin widths. If this number is close to 95, the CI was estimated correctly. If this number is much smaller, the CI was underestimated. The noise amplitude was set to 0.01.

Table 2

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values when the correct model is used, for bin widths of 10, 15, and 30pixels . Standard deviations are indicated in parentheses. Noise amplitude=0.01 .

	10	15	30
$a_1$	94 (4)	96 (4)	96 (7)
$a_2$	91 (4)	94 (4)	94 (5)
$a_3$	96 (3)	97 (4)	95 (8)
$a_4$	94 (3)	95 (5)	84 (25)
$a_5$	94 (3)	95 (3)	95 (7)
$a_6$	94 (3)	95 (5)	85 (22)
$a_7$	93 (4)	95 (3)	94 (9)

Table 2 shows that on average, close to 95% of the fitted parameter values is within 2 CIs, except for the largest bin width for parameters $a_4$ (betacarotene concentration) and $a_6$ (blood saturation). For this large bin width, the optical resolution of the binned data $\left(12\mathrm{nm}\right)$ becomes larger than the full width at half maximum (FWHM) of the narrowest betacarotene feature $\left(8\mathrm{nm}\right)$ , causing an underestimation in the CI of the betacarotene concentration $a_4$ . Since betacarotene absorbs in the same wavelength region as oxyhemoglobin and deoxyhemoglobin, the accuracy of the blood saturation $a_6$ is affected as well. From Table 2 it is clear that the best results were obtained for a bin width of 15, which is used in the remainder of this work.

4.1.2.

Influence of noise amplitude

Table 3 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for the 315 synthetic datasets (standard deviations are indicated in parentheses) for different noise amplitudes and a bin width of $15\text{pixels}$ . On average, close to 95% of the fitted parameter values is within 2 CIs, even for very high noise amplitudes. Table 4 shows the average CI for each of the parameters as a function of blood volume fraction $a_5$ for the three different noise amplitudes. The CI increases with increasing noise amplitude, as expected. Thus for large noise amplitudes (0.1 and 0.5, corresponding to a signal-to-noise ratio at $\lambda =800\mathrm{nm}$ of 10 and 2, respectively), the fitted parameters are also correctly within 2 CIs in not significantly less than 95% of the simulations, because the CIs are dramatically increased compared to the low noise amplitude (0.01, or signal to noise of 100). Note that for the low noise amplitude, which is typical for our in vivo data,^{11, 12, 13, 14} the CI of the fitted blood saturation $a_6$ is smaller than 1% for blood volume fractions $a_5>1\%$ . Even for a blood volume fraction as low as 0.2%, the CI of $a_6$ is only 5% when the correct model is used in the fit. These numbers indicate that our instrument can extract the in vivo microvascular saturation with high accuracy, under the assumption that the model used to fit the data is correct.

Table 3

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values when the correct model is used, for noise amplitudes of 0.01, 0.1, and 0.5. Bin width=15pixels .

	0.01	0.1	0.5
$a_1$	96 (4)	93 (3)	93 (3)
$a_2$	94 (4)	93 (3)	93 (3)
$a_3$	97 (4)	93 (3)	93 (3)
$a_4$	95 (5)	93 (3)	92 (3)
$a_5$	95 (3)	93 (4)	91 (6)
$a_6$	95 (5)	94 (2)	94 (3)
$a_7$	95 (3)	92 (3)	92 (3)

Table 4

Average CI for each of the parameters as a function of blood volume fraction a5 for the three different noise amplitudes.

	Noiseamplitudes	0%	0.2%	1%	2%	5%	10%	20%
$a_1$ [−]	0.010.100.50	0.30.52.5	0.00030.0022.4	0.00030.0020.01	0.00030.0020.01	0.00040.0020.01	0.00060.0030.02	0.0010.010.05
$a_2$ [−]	0.010.100.50	0.0020.010.07	0.0010.010.07	0.0010.010.06	0.0010.010.06	0.0020.010.07	0.0020.020.09	0.0040.030.2
$a_3$ [−]	0.010.100.50	0.00050.0020.08	0.00020.0010.02	0.00020.0010.006	0.00020.0010.006	0.00030.0020.008	0.00050.0030.01	0.0010.010.05
$a_4$ $\left[\mu \mathrm{M}\right]$	0.010.100.50	0.150.94.5	0.10.94.6	0.11.04.9	0.21.15.3	0.31.57.0	0.42.311	0.96.030
$a_5$ [%]	0.010.100.50	${10}^4$ ${10}^4$ ${10}^5$	0.0080.07 ${10}^4$	0.0090.080.4	0.010.080.4	0.010.10.5	0.020.140.7	0.040.31.4
$a_6$ [%]	0.010.100.50	${10}^6$ ${10}^8$ ${10}^7$	5.034 ${10}^7$	1.0737	0.6420	0.42.411	0.42.010	0.42.110
$a_7$ $\left[\mu \mathrm{m}\right]$	0.010.100.50	${10}^6$ ${10}^7$ ${10}^7$	0.86 ${10}^3$	0.21.37	0.10.74	0.050.42	0.040.31.7	0.080.73.4

Figures 1a and 1b show a typical fit and the residuals $(\equiv \text{data}\text{-}\text{fit})$ , respectively, for the synthetic spectrum with $a_3=0.01$ , $a_4=5\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ for a noise amplitude of 0.01 and a bin width of 15. The fitted values for the parameters and the CIs were $a_1=1.0003\pm 0.0002$ , $a_2=-1.000\pm 0.001$ , $a_3=0.0099\pm 0.0002$ , $a_4=5.0\pm 0.2\mu \mathrm{M}$ , $a_5=2.00\pm 0.01\%$ , $a_6=74.9\pm 0.7\%$ , and $a_7=10.04\pm 0.05\mu \mathrm{m}$ . All fit values are unbiased and within 2 CIs of the true values. The residual spectrum shows only noise and is featureless.

Fig. 1

(a) Typical fit and (b) the residuals for the synthetic spectrum with $a_3=0.01$ , $a_4=5\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ for a noise amplitude of 0.01 and a bin width of $15\text{pixels}$ .

In conclusion, when the correct model is used to fit the data, the fitted parameter values are unbiased, and the CIs are estimated correctly using the diagonal elements of estimated covariance matrix obtained by multiplying the inverse of the second derivative matrix of ${\chi }^2$ with respect to its free parameters by ${\chi }^2∕v$ , with $v$ being the number of degrees of freedom.

4.2.

Fitting with the Wrong Model: Missing Absorber

In this section, we investigate the magnitude of the bias in the fitted parameters introduced when the fitting model is missing an absorber. Two different situations are investigated: the absorption band of the missing absorber has 1. little or 2. significant spectral overlap with the other absorbers in the component fit. In Sec. 4.2.1 the missing absorber is an artificial Gaussian centered at $700\mathrm{nm}$ with a width of $25\mathrm{nm}$ . Such an absorber has only little overlap with the main absorption bands of oxyhemoglobin and deoxyhemoglobin. In Sec. 4.2.2 the missing absorber is betacarotene, which does overlap with the absorption bands of hemoglobin. Figure 2 shows the normalized absorption coefficients of the four absorbers used in the simulations. The noise amplitude is set to 0.01 and the bin width is 15 in this section.

Fig. 2

Normalized absorption coefficients of the four absorbers used in the simulations.

4.2.1.

Nonoverlapping Absorber: Gaussian Centered at $700nm$

An additional 210 synthetic datasets were generated using the Gaussian absorber instead of betacarotene in Eq. 2. Fitting was performed with only six free parameters by excluding $a_4\cdot {\mu }_a^{\mathrm{Bcar}}\left(\lambda \right)$ from Eq. 2. Figures 3a and 3b show a typical fit and the residuals, respectively, for the synthetic spectrum with $a_3=0.01$ , $a_4=5$ , $a_5=2\%$ , and $a_6=75\%$ . This “concentration” $(a_4=5)$ of Gaussian absorber results in a hardly visible maximum signal decrease of 2.5% at $\lambda =700\mathrm{nm}$ . The fitted values for the parameters and the CIs were $a_1=0.9977\pm 0.0007$ , $a_2=-0.992\pm 0.004$ , $a_3=0.0098\pm 0.0005$ , $a_5=1.88\pm 0.03\%$ , $a_6=74\pm 2\%$ , and $a_7=9.6\pm 0.2\mu \mathrm{m}$ .

Fig. 3

(a) and (c) typical fit and (b) and (d) the residuals for the synthetic spectrum with $a_3=0.01$ , $a_4=5$ and 30, $a_5=2\%$ , and $a_6=75\%$ , fitted without the Gaussian absorber in the fitting model.

Figures 3c and 3d show a typical fit and the residuals, respectively, for the synthetic spectrum with $a_3=0.01$ , $a_4=30$ , $a_5=2\%$ , and $a_6=75\%$ . This “concentration” $(a_4=30)$ of Gaussian absorber results in a clearly visible maximum signal decrease of 14% at $\lambda =700\mathrm{nm}$ . The fitted values for the parameters and the CIs were $a_1=0.982\pm 0.004$ , $a_2=-0.93\pm 0.02$ , $a_3=0.011\pm 0.003$ , $a_5=1.3\pm 0.2\%$ , $a_6=62\pm 11\%$ , and $a_7=7.4\pm 0.8\mu \mathrm{m}$ . Even though the statistical CIs of these fits are much larger than the statistical CIs of the parameters using the correct model (Sec. 4.1.2), only the Rayleigh amplitude $\left(a_3\right)$ and saturation $\left(a_6\right)$ fit values are within 2 CIs of the true values. The residual spectra very clearly show the features of hemoglobin absorption as well as the missing Gaussian absorber.

Table 5 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for the 105 synthetic datasets for different concentrations of Gaussian absorber.

Table 5

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values for different concentrations of missing Gaussian absorber.

	0	5	30
$a_1$	96 (4)	6 (21)	3 (16)
$a_2$	94 (3)	16 (36)	28 (44)
$a_3$	97 (4)	86 (31)	78 (40)
$a_4$
$a_5$	95 (3)	2 (7)	3 (15)
$a_6$	95 (5)	80 (36)	71 (44)
$a_7$	95 (3)	1 (3)	1 (3)

Table 6 shows the average bias ( $\equiv \text{absolute}$ value of absolute difference between true and fitted parameter values) and CI for each of the parameters for different concentrations of Gaussian absorber. Blood volume fractions $a_5<1\%$ were excluded from the calculation of the averages due to the extremely large bias and CI of the blood parameters ( $a_5$ , $a_6$ , and $a_7$ ) for small blood volume fractions (see for example Table 4). It is observed that, except for parameters Rayleigh amplitude $\left(a_3\right)$ and saturation $\left(a_6\right)$ , the average bias in the parameters is much larger than twice the average statistical CI of the parameters.

Table 6

Average bias and CI for each of the parameters for different concentrations of missing Gaussian absorber.

	Gauss absorberconcentration	Average bias	Average CI
$a_1$ [−]	530	0.0050.03	0.0020.008
$a_2$ [−]	530	0.020.13	0.0070.04
$a_3$ [−]	530	0.00160.009	0.00150.007
$a_4$
$a_5$ [%]	530	0.21.0	0.060.3
$a_6$ [%]	530	1.39	1.59
$a_7$ $\left[\mu \mathrm{m}\right]$	530	0.94.6	0.31.5

4.2.2.

Overlapping Absorber: Betacarotene

Similar to Sec. 4.2.1, fitting was performed with only six free parameters by excluding $a_4{\mu }_a^{\mathrm{Bcar}}\left(\lambda \right)$ from Eq. 2. Figures 4a and 4b show a typical fit and the residuals, respectively, for the synthetic spectrum with $a_3=0.01$ , $a_4=5\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ . This concentration $(a_4=5\mu \mathrm{M})$ of betacarotene results in a hardly visible maximum signal decrease of 2.5% at $\lambda =450\mathrm{nm}$ . The fitted values for the parameters and the CIs were $a_1=0.9987\pm 0.0006$ , $a_2=-0.984\pm 0.003$ , $a_3=0.0103\pm 0.0004$ , $a_5=2.02\pm 0.03\%$ , $a_6=91\pm 1\%$ , and $a_7=9.9\pm 0.1\mu \mathrm{m}$ . Even though the statistical CIs of this fit are larger than the statistical CIs of the parameters using the correct model (Sec. 4.1.2), only half of the fit values ( $a_3$ , $a_5$ , and $a_7$ ) are within 2 CIs of the true values. The fitted saturation $\left(a_6\right)$ is very biased with respect to the true value (91 versus 75%), while the CI is increased only marginally compared to the case where betacarotene was included in the fit. The residual spectrum clearly shows the features of hemoglobin absorption as well as betacarotene absorption.

Fig. 4

(a) and (c) typical fit and (b) and (d) the residuals for the synthetic spectrum with $a_3=0.01$ , $a_4=5$ and $30\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ , fitted without betacarotene in the fitting model.

Figures 4c and 4d show a typical fit and the residuals, respectively, for the synthetic spectrum with $a_3=0.01$ , $a_4=30\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ . This concentration $(a_4=30\mu \mathrm{M})$ of betacarotene results in a clearly visible maximum signal decrease of 14% at $\lambda =450\mathrm{nm}$ . The fitted values for the parameters and the CIs were $a_1=0.995\pm 0.003$ , $a_2=-0.91\pm 0.02$ , $a_3=0.012\pm 0.002$ , $a_5=2.4\pm 0.1\%$ , $a_6=155\pm 7\%$ , and $a_7=10.7\pm 0.6\mu \mathrm{m}$ . Even though the statistical CIs of this fit are much larger than the statistical CIs of the parameters using the correct model (Sec. 4.1.2), only half of the fit values ( $a_1$ , $a_3$ , and $a_7$ ) are within 2 CIs of the true values. The fitted saturation $\left(a_6\right)$ is very biased with respect to the true value (155 versus 75%), while the CI is increased only marginally compared to the case where betacarotene was included in the fit. Note that the fitted blood volume fraction $a_5$ has a smaller bias than expected from the figure; this is due to the fitted negative amount of deoxyhemoglobin (which also resulted in the physically impossible 155% saturation). The residual spectrum very clearly shows the features of hemoglobin absorption as well as betacarotene absorption.

Table 7 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for the 105 synthetic datasets for different concentrations of betacarotene.

Table 7

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values for different concentrations of missing betacarotene.

	0	5	30
$a_1$	96 (4)	24 (37)	23 (40)
$a_2$	94 (3)	48 (44)	42 (44)
$a_3$	97 (4)	43 (42)	38 (44)
$a_4$
$a_5$	95 (3)	58 (39)	46 (42)
$a_6$	95 (5)	10 (29)	22 (40)
$a_7$	95 (3)	45 (42)	39 (41)

Table 8 shows the average bias and CI for each of the parameters for different concentrations of betacarotene, for blood volume fractions $a_5>1\%$ .

Table 8

Average bias and CI for each of the parameters for different concentrations of missing betacarotene.

	Betacaroteneconcentration	Average bias	Average CI
$a_1$ [−]	530	0.0030.02	0.0010.004
$a_2$ [−]	530	0.0120.07	0.0040.02
$a_3$ [−]	530	0.0030.02	0.00070.003
$a_4$
$a_5$ [%]	530	0.060.38	0.030.14
$a_6$ [%]	530	844	16
$a_7$ $\left[\mu \mathrm{m}\right]$	530	0.52.1	0.21.2

It is observed that the average bias in the parameters is larger than twice the average statistical CI of the parameters, especially for the saturation $a_6$ .

4.3.

Fitting with the Wrong Model: Incorrect Scattering Model

In this section, we investigate the magnitude of the bias in the fitted parameters that is introduced when a wrong scattering function is used in the fitting routine. This is established by fitting with only six free parameters by excluding the Rayleigh scattering term $a_3{(\lambda ∕{\lambda }_0)}^{-4}$ in Eq. 1 from the fit function. Figures 5a and 5b show a typical fit and the residuals, respectively, for the synthetic spectrum with $a_3=0.01$ , $a_4=5\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ . The fitted values for the parameters and the CIs were $a_1=1.0092\pm 0.0009$ , $a_2=-1.052\pm 0.006$ , $a_4=6.0\pm 0.9\mu \mathrm{M}$ , $a_5=1.75\pm 0.05\%$ , $a_6=73\pm 4\%$ , and $a_7=9.6\pm 0.3\mu \mathrm{m}$ .

Fig. 5

(a) and (c) typical fit and (b) and (d) the residuals for the synthetic spectrum with $a_3=0.01$ and 0.1, $a_4=5\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ , fitted without Rayleigh scattering in the fitting model.

Figures 5c and 5d show a typical fit and the residuals, respectively, for the synthetic spectrum with $a_3=0.1$ , $a_4=5\mu \mathrm{M}$ , $a_5=2\%$ , and $a_6=75\%$ . The fitted values for the parameters and the CIs were $a_1=1.105\pm 0.003$ , $a_2=-1.35\pm 0.02$ , $a_4=-10\pm 3\mu \mathrm{M}$ , $a_5=1.2\pm 0.2\%$ , $a_6=65\pm 21\%$ , and $a_7=10\pm 2\mu \mathrm{m}$ . Even though the statistical CIs of these fits are much larger than the statistical CIs of the parameters using the correct model (Sec. 4.1.2), only half the fit values are within 2 CIs of the true values. The blood volume fraction is particularly underestimated with respect to the true value. The residual spectrum very clearly shows the features of hemoglobin absorption as well as a strong peak at the start of the spectrum associated with an incorrect scattering function.

Table 9 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for 70 synthetic datasets (only the extreme betacarotene concentrations 0 and $30\mu \mathrm{M}$ were considered) for different Rayleigh amplitudes.

Table 9

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values for different Rayleigh amplitudes with the Rayleigh term omitted from the fit function.

	0	0.01	0.1
$a_1$	91 (3)	0 (0)	0 (0)
$a_2$	91 (3)	0 (0)	0 (0)
$a_3$
$a_4$	94 (3)	62 (43)	9 (27)
$a_5$	94 (3)	11 (25)	1 (4)
$a_6$	95 (3)	86 (25)	80 (34)
$a_7$	94 (3)	62 (42)	69 (43)

Table 10 shows the average bias and CI for each of the parameters for different Rayleigh amplitudes. Again, blood volume fractions $a_5<1\%$ were excluded from the calculation of the averages due to the extremely large bias and CI of the blood parameters ( $a_5$ , $a_6$ , and $a_7$ ) for small blood volume fractions. For the largest Rayleigh amplitude, blood volume fractions $a_5<2\%$ were excluded from the calculation of the averages due to the large bias and CI of the blood parameters ( $a_5$ $a_6$ , and $a_7$ ) for these blood volume fractions as well.

Table 10

Average bias and CI for each of the parameters for different Rayleigh amplitudes.

	Rayleighamplitude	Average bias	Average CI
$a_1$ [−]	0.010.1	0.00980.103	0.00060.002
$a_2$ [−]	0.010.1	0.0420.30	0.0040.02
$a_3$ [−]
$a_4$ $\left[\mu \mathrm{M}\right]$	0.010.1	219	14
$a_5$ [%]	0.010.1	0.231.4	0.050.2
$a_6$ [%]	0.010.1	310	39
$a_7$ $\left[\mu \mathrm{m}\right]$	0.010.1	0.71.8	0.41.3

It is observed that the scattering parameters ( $a_1$ and $a_2$ ) as well as the absorber concentrations ( $a_4$ and $a_5$ ) are very biased compared to the calculated statistical CIs, whereas the blood saturation and vessel diameter ( $a_6$ and $a_7$ , respectively) are not affected as much by imperfections in the scattering function.

4.4.

Identification and Minimization of Bias

In the previous sections we have shown that when the fitting model is correct, the fitted parameter values are unbiased and the CIs can be estimated correctly using the reduced chi-square normalized covariance matrix. When the fitting model is incorrect, however, the fitted parameters are biased and the calculated statistical CIs cannot correctly account for this bias. This is expected, since the systematic errors associated with using the wrong model are not fully accounted for by our purely statistical analysis. It is practically impossible to account for systematic errors using statistical methods only, and hence it should be realized that the calculated CIs represent the lower limits of the true CIs of the parameters. In this section we present a method to identify and minimize systematic errors in the fitting routine.

Identification of systematic errors is easily performed based on the shape of the residual spectrum. When the fitting model is correct, the residual spectrum is featureless and shows random noise only [Fig. 1b]. When the fitting model is incorrect, the residual spectrum will always show features of the chromophore absorption coefficients, and in case the wrong scattering function is used, an additional (slow) background variation over the whole wavelength range can be observed [Figs. 3b, 3d, 5b, 5d]. The major difference in the residual spectra of fits with missing absorbers [Figs. 3b, 3d, 4b, 4d] and fits with an incorrect scattering function [Figs. 5b and 5d] is that in the first case, the residuals have a (strong) negative peak near the location of the absorption maximum of the missing absorber, while in the second case all peaks in the residuals are near the absorption maxima of the chromophores used in the fit. In Sec. 4.4.1, we discuss an approach to minimize the systematic errors associated with missing absorbers, and in Sec. 4.4.2, we discuss a method to minimize the systematic errors associated with an incorrect scattering function.

4.4.1.

Minimization of Bias: Missing Absorber

A missing absorber in the model can be identified by a strong negative peak in the residuals (if defined as data minus fit) at a location away from absorption maxima of the chromophores used in the fit. After identification of a missing absorber, one can reduce the systematic error by refitting the spectrum, but now including an additional absorber that peaks near the wavelength of the negative peak in the residuals. Since a Gaussian function is a good first approximation for the shape of an absorption band, we propose to use a Gaussian function for this unidentified absorber, following Eq. 4:

Eq. 4

$${\mu }_a^{\text{unknown}}\left(\lambda \right)=a_8\exp [-0.5{\left(\frac{\lambda -{\lambda }_{\text{peak}}-a_9}{a_{10}}\right)}^2].$$

Here ${\lambda }_{\text{peak}}$ is the location of the negative peak in the residual spectrum and $a_8$ , $a_9$ , and $a_{10}$ are additional free parameters representing the “concentration,” the peak shift, and the width of the unknown absorber, respectively. The peak shift parameter $a_9$ accounts for the fact that the location of the negative peak in the residual spectrum is not necessarily the exact location of the absorption maximum of the absorber, but will be close to it. Therefore, this parameter should be constrained to small values to force the location of the absorption maximum of the Gaussian to be near the location of the negative peak in the residual spectrum. We have used $a_9<10\mathrm{nm}$ in our analysis, but did not investigate in detail how the results depend on this constraint.

When this strategy is employed for the case of the missing absorbers in the fitting routine (Sec. 4.2), dramatic improvement in the fit quality is achieved. Table 11 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for the 105 synthetic datasets for different concentrations of Gaussian absorber. It is observed that on average, nearly 95% of the fitted parameter values is within 2 CIs of their true values, which indicates that the parameters and their CIs are estimated correctly.

Table 11

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values for different concentrations of Gaussian absorber with an additional absorber in the fit function.

	0	5	30
$a_1$	96 (4)	93 (4)	92 (6)
$a_2$	94 (3)	94 (6)	94 (7)
$a_3$	97 (4)	95 (8)	94 (10)
$a_4$
$a_5$	95 (3)	91 (8)	89 (12)
$a_6$	95 (5)	91 (10)	89 (11)
$a_7$	95 (3)	93 (5)	93 (6)

The location of the peak ${\lambda }_{\text{peak}}$ was determined by the mimimum of the residual spectra, and the absorption maximum of the unknown absorber is ${\lambda }_{\text{peak}}-a_9$ . For the small Gaussian absorber concentration $a_4=5$ , the average fitted values of the missing absorber were ${\lambda }_{\text{peak}}-a_9=701\mathrm{nm}$ and width $a_{10}=24\mathrm{nm}$ . For the larger Gaussian absorber concentration $a_4=30$ , the average fitted values of the missing absorber were ${\lambda }_{\text{peak}}-a_9=699\mathrm{nm}$ and width $a_{10}=27\mathrm{nm}$ . Thus in both cases the absorption maxima and widths of the unknown absorber reproduce the true absorption spectrum very accurately.

Table 12 shows the number of times each of the fitted parameter values is within 2 CIs of their true values, out of the 100 fits and averaged for the 105 synthetic datasets for different concentrations of betacarotene.

Table 12

Average number of times (out of 100) each of the fitted parameter values is within 2 CIs of their true values for different concentrations of betacarotene with an additional absorber in the fit function.

	0	5	30
$a_1$	96 (4)	72 (16)	39 (35)
$a_2$	94 (3)	53 (22)	11 (20)
$a_3$	97 (4)	38 (24)	9 (19)
$a_4$
$a_5$	95 (3)	42 (21)	17 (23)
$a_6$	95 (5)	60 (33)	55 (33)
$a_7$	95 (3)	44 (26)	25 (32)

The location of the peak ${\lambda }_{\text{peak}}$ was determined by the mimimum of the residual spectra, and the absorption maximum of the unknown absorber is ${\lambda }_{\text{peak}}-a_9$ . For the small betacarotene concentration $a_4=5\mu \mathrm{M}$ , the average fitted values of the missing absorber were ${\lambda }_{\text{peak}}-a_9=453\mathrm{nm}$ and width $a_{10}=40\mathrm{nm}$ . For the larger betacarotene concentration $a_4=30\mu \mathrm{M}$ , the average fitted values of the missing absorber were ${\lambda }_{\text{peak}}-a_9=455\mathrm{nm}$ and width $a_{10}=37\mathrm{nm}$ . These values were determined from the fits with blood volume fractions $a_5<10\%$ only; for larger blood volume fractions, the minimum of the residual spectra was colocated with the blood absorption peaks. As a result, no missing absorber could be identified for blood volume fractions $a_5⩾10\%$ . Figure 6 shows the normalized betacarotene absorption spectrum and the average extracted absorption spectra of the unknown Gaussian absorber determined from the fits. Good agreement in both the peak position and the width is found between the true betacarotene absorption coefficient and the absorption coefficients of the unknown Gaussians determined from the fits. Despite these promising results, Table 12 shows that on average, much less than 95% of the fitted parameter values is within 2 CIs of their true values. This indicates that there are still systematic errors and that the resulting fit values are biased. The disappointing results are also due to the fact that the method fails to identify a missing absorber for large blood volume fractions. Table 13 shows the average bias ( $\equiv \text{absolute}$ value of the absolute difference between true and fitted parameter values) and CI for each of the parameters for different concentrations of betacarotene, with an unknown Gaussian absorber in the fit that has a maximum absorption at a wavelength no more than $10\mathrm{nm}$ away from the minimum of the residuals, as explained before. Blood volume fractions $a_5<1\%$ were excluded from the calculation of the averages, as discussed previously. It is observed that the average bias in the parameters is mostly larger than twice the average statistical CI of the parameters, but much smaller than before (Table 8) when no additional Gaussian absorber was included in the fit. In particular, the bias on the saturation is decreased by a large amount.

Fig. 6

Normalized betacarotene absorption spectrum and the average extracted absorption spectra of the unknown Gaussian absorber determined from the fits. Solid line: true betacarotene spectrum. Dotted line: Gaussian approximation from fits with large betacarotene concentration. Dashed line: Gaussian approximation from fits with small betacarotene concentration.

Table 13

Average bias and CI for each of the parameters for different concentrations of betacarotene, with an unknown Gaussian absorber in the fit, which has a maximum absorption at a wavelength no more than 10nm away from the minimum of the residuals.

	Betacaroteneconcentration	Average bias	Average CI
$a_1$ [−]	530	0.0020.005	0.0080.001
$a_2$ [−]	530	0.0040.011	0.0040.005
$a_3$ [−]	530	0.0010.006	0.0010.001
$a_4$
$a_5$ [%]	530	0.040.24	0.020.05
$a_6$ [%]	530	1.53.5	0.61.5
$a_7$ $\left[\mu \mathrm{m}\right]$	530	0.31.6	0.10.3

4.4.2.

Minimization of Bias: Wrong Scattering Model

When a wrong scattering model is used, the residuals show features of the chromophore absorption coefficients on top of an additional background variation over the whole wavelength range. After identification of a wrong scattering model, one can reduce the systematic error by refitting the spectrum, but now including an additional scattering term. We propose to use a second-order polynomial function as an additional scattering term, since it provides a good first-order approximation to any nonperiodic continuous function, such as usually encountered in tissue scattering:

Eq. 5

$$R\left(\lambda \right)=[a_1{\left(\frac{\lambda }{{\lambda }_0}\right)}^{a_2}+a_3{\left(\frac{\lambda }{{\lambda }_0}\right)}^2+a_8\left(\frac{\lambda }{{\lambda }_0}\right)+a_9]\exp [-0.4{\mu }_a^{\text{total}}\left(\lambda \right)].$$

Thus the scattering function is now modeled by Mie scattering plus a second-order polynomial. Whereas the Mie scattering model was based on physical principles where the parameters $a_1$ and $a_2$ are associated with the density and size distribution of scatterers,²¹ respectively, the new scattering model is primarily designed to describe the scattering function mathematically without associating physical meaning to the parameters. In other words, the scattering function is modified with the aim to reduce the bias in the concentration estimates of the chromophores, at the cost of modeling with a physically uninterpretable scattering function. Table 14 shows the number of times the fitted microvascular parameter values are within 2 CIs of their true values, out of the 100 fits and averaged for 35 synthetic datasets (betacarotene was exluded in this subanalysis for convenience) for different Rayleigh amplitudes. It is observed that on average, the in the case of a small Rayleigh component, nearly 95% of the fitted parameter values is within 2 CIs of their true values, which indicates that the parameters and their CIs are estimated correctly. For the larger Rayleigh amplitude, the accuracy of the results decreases but is still much better than before (Table 9) without the additional second-order polynomial scattering term. As for the Mie and Rayleigh scattering parameters, it should not be expected that they are accurately reproduced, since the mathematical scattering function used in the fit is very different from the scattering function of the synthetic datasets. Finally, Fig. 7a shows the synthetic scattering function

$$S_{\mathrm{syn}}\left(\lambda \right)=1{\left(\frac{\lambda }{{\lambda }_0}\right)}^{-1}+0.1{\left(\frac{\lambda }{{\lambda }_0}\right)}^{-4},$$

and a fit of this function using the scattering function of Eq. 5. Figure 7b shows the residuals. Clearly the scattering function of Eq. 5 is a good but not perfect approximation for $S_{\mathrm{syn}}\left(\lambda \right)$ , leading to the results of Table 14.

Fig. 7

(a) The synthetic scattering function $S_{\mathrm{syn}}\left(\lambda \right)$ and a fit of this function using the scattering function of Eq. 5, and (b) the residuals.

Table 14

Average number of times (out of 100) the fitted microvascular parameter values are within 2 CIs of their true values for different Rayleigh amplitudes fitted with an additional second order polynomial scattering function.

	0	0.01	0.1
$a_5$	96 (3)	92 (5)	81 (20)
$a_6$	95 (3)	92 (3)	80 (10)
$a_7$	95 (3)	89 (6)	75 (24)

5.

Discussion and Conclusion

We have validated a simple method for determining the confidence intervals on fitted parameters derived from modeling optical reflectance spectroscopic measurements using synthetic datasets. The method estimates the parameter confidence intervals as the square roots of the diagonal elements of the covariance matrix, which is obtained by multiplying the inverse of the second derivative matrix of ${\chi }^2$ with respect to its free parameters by ${\chi }^2∕v$ , with $v$ being the degrees of freedom. Most commercial fitting software packages, employing least squares analysis, automatically generate the covariance matrix. CIs are easily extracted after multiplication of the diagonal elements with ${\chi }^2∕v$ to account for the uncertainty in the estimation of the statistical fluctuations and corresponding standard deviations $s\left(i\right)$ of the data points. We have shown in Sec. 4.1 that this method yields correct CIs for our nonlinear fitting function, as long as the model used to describe the data is correct. The statistical CIs do not depend on the number of pixels used to bin the data, as long as the optical resolution associated with this number does not exceed the optical resolution of the biological features. Furthermore, the CI of the fitted parameters correctly increases with decreasing signal-to-noise ratio, while in all cases the fitted parameter values are unbiased with respect to their true values. The CI calculated in this way represents the statistical error on the fit parameters and is a lower boundary of the true fit parameter uncertainty. We have also investigated the sensitivity of our parameters to systematic errors associated with incorrect models in results Secs. 4.2, 4.3. Even when there are only small amounts of missing absorber (Sec. 4.2), such that the fit does not deviate from the data by more than 2 to 3% over the entire wavelength range [Figs. 3a, 3b, 4a, 4b], the fitted parameters may be biased by large amounts compared to the calculated statistical CIs. Figures 4a and 4b are particularly instructive: the quality of the fit looks excellent, the residuals are smaller than 2.5% over the entire wavelength range, and yet the bias in the saturation is more than 15%. Thus imperfections in the fitting model introduce a bias in the fitted parameters that greatly exceeds the estimated statistical CIs. Since multicomponent spectral fitting of in vivo optical spectroscopic data is particularly sensitive to these problems, we have also introduced methods to identify and subsequently minimize the bias in the fitted parameters associated with incorrect modeling (Sec. 4.4). Identification of systematic errors is based on the shape of the residual spectrum. When the fitting model is correct, the residual spectrum is featureless and shows random noise only [Fig. 1b]. When the fitting model is incorrect, the residual spectrum always shows features of the included chromophore absorption coefficients, and in case the wrong scattering function is used, an additional background variation over the whole wavelength range can be observed [Figs. 3b, 3d, 5b, 5d]. In the case of a missing absorber, a strong negative peak in the residuals at a location away from absorption maxima of the chromophores used in the fit is observed. The systematic error associated with a missing absorber was reduced by refitting the spectrum, but now including an additional absorber that peaks near the wavelength of the minimum (maximum negative) in the residuals. Dramatic improvement in the fit quality was achieved when using a Gaussian absorber with three free parameters: amplitude, peak shift with respect to minimum of residuals, and peak width. In the case where the missing absorber truly was a Gaussian, the fits were perfect, and nearly 95% of the fitted parameter values was within 2 CIs of their true values, which indicates that the fitted parameters and their CIs were estimated correctly. Furthermore, the absorption maxima and widths of the unknown Gaussian absorber reproduced the true Gaussian absorption spectrum very accurately. In the case where the missing absorber was betacarotene, which has a non-Gaussian absorption spectrum overlapping the absorption spectra of oxy- and deoxyhemoglobin, the results were not as good. For large blood volume fractions, the minimum of the residual spectra was colocated with the blood absorption peaks, and no missing absorber could be identified for large $(a_5⩾10\%)$ blood volume fractions. For smaller blood volume fractions, the extracted absorption spectra of the unknown Gaussian absorber determined from the fits showed good agreement in both the peak position and the width with the true betacarotene absorption coefficient. Despite this good agreement, much less than 95% of the fitted parameter values was within 2 CIs of their true values for both high and low betacarotene concentrations. This indicates that there are still systematic errors and that the resulting fit values are biased. However, the bias in the parameters was shown to be much smaller than without an additional Gaussian absorber included in the fit; the bias on the saturation was decreased by a large amount (from 8 to 1.5% for small betacarotene concentration, and from 44 to 3.5% for large betacarotene concentration).

Palmer ²² have previously used a similar approach in a clinical dataset of malignant and nonmalignant breast tissue. They included a Gaussian absorber with fixed width and peak position in their fits based on the shape of their previously obtained residuals in a large clinical dataset, under the assumption that the spectral shape of the missing absorber is the same across the whole set of tissue samples measured. Interestingly, they found that although the quality of the fits was improved with the addition of the Gaussian absorber, the conclusions regarding the extracted absorber and scattering parameters from the fits were not significantly affected by its addition, and so all subsequent analysis was carried out without the inclusion of this Gaussian absorber. Based on our present research, this result is surprising. Possibly this statement was made with respect to the average of the fitted parameters over the entire dataset. In that case, the biases in the fitted parameters for the individual spectra may be canceled on averaging. Furthermore, the biological variation in the parameters for multiple samples is probably larger than the bias associated with the systematic error introduced by the incomplete fitting model, as we describe next.

The systematic error associated with a wrong scattering model was reduced by refitting the spectrum, but now including a second-order polynomial as an additional scattering term. The new scattering model is primarily designed to describe the scattering function mathematically without associating physical meaning to the scattering parameters, with the aim to reduce the bias in the concentration estimates of the chromophores. The implementation of this polynomial scattering function greatly reduced the bias in the microvascular parameters, although the approach has its limits for large deviations in the model scattering function from the true scattering function. The mathematical description of the “true” scattering function in case of biological (in vivo) tissue actually depends on the measurement geometry and on the tissue itself. In the case of diffuse measurements, the power law behavior of scattering seems to be an adequate description of the scattering function. In the case of nondiffuse measurements, more complicated functions (including periodic elements^{23, 24, 25}) may be necessary to fully describe the scattering function.

Other potential sources of systematic error that are not investigated in detail in this work are: 1. incorrect wavelength calibration of the spectrometer, and 2 incorrect or unknown specific absorption coefficients of the chromophores. Regarding the first point, we observed that a $1\text{-}\mathrm{nm}$ offset in the spectrometer calibration typically resulted in a 3% bias in the fitted saturation $a_6$ (data not shown). Since it should always be possible to calibrate the wavelength axis of a spectrometer to within $1\mathrm{nm}$ , we consider this source of error to be minor. Regarding the second point, the specific absorption coefficients of chromophores are usually measured when dissolved in water, ethanol, hexane, or chloroform. Apart from water, these environments are dramatically different from the in vivo environment, and it is not easy to predict how the environment changes the specific absorption coefficient of the chromophore. Spectral shifting as well as spectral broadening may occur, which would also result in biased fit parameters, the amount of bias depending on the amounts of shifting and broadening. Therefore, it is imperitive to use specific chromophore absorption coefficients in the fitting routine that were measured in an evironment that comes as close to the in vivo environment as possible. This is particularly important for component fits of in vivo fluorescence and Raman spectra, whose shapes can be altered even more dramatically than extinction spectra due changes in the microenvironment. Alternatively, any features in the residual spectra obtained from fitting in vivo spectra using ex vivo component spectra could be used to modify these component spectra and determine the true in vivo fluorescence, Raman, or extinction spectra. Moreover, environmental changes (such as pH variations) that are reflected in altered spectral shapes of the components will lead to features in the residual spectra. This provides an opportunity to characterize and monitor the microenvironment by analyzing the shape of the residual spectra.

As can be seen in Table 4, for a noise amplitude typical for our in vivo data, the CI of the fitted parameters is quite small (e.g., the CI for blood saturation $a_6$ is smaller than 1% for blood volume fractions $a_5>1\%$ ). This CI represents the statistical error associated with an individual spectrum. It is critical that this is interpreted differently from a standard deviation calculated from averaging multiple measurements. The latter standard deviation represents the biological variation in the parameters, which for most applications will be larger than the statistical uncertainty in the parameters related to the quality of the data. In fact, in our experience the biological variation in the parameters for multiple measurements is even larger than the bias associated with the systematic errors discussed in results Secs. 4.2, 4.3, and for most applications it will be sufficient to specify these biological standard deviations only. However, completely disregarding the statistical CIs can also cause problems. For example, if the blood volume fraction is very low, the statistical uncertainty of the blood related fit values becomes very large (Table 4). Therefore, in the case of averaging multiple measurements, we feel that calculating the weighted average of the parameters, with the statistical CIs as weight factors, is more appropriate than calculating the unweighted average and ignoring the statistical CIs. Furthermore, calculating the statistical errors in the fit parameters facilitates objective assessment of the quality of spectra, and fits with the possibility to reject poor quality spectra or fit values from an in vivo dataset, thereby avoiding unnecessary pollution of the dataset by poor quality spectra.

Table 8 shows that a missing absorber that overlaps with the spectra of oxy- and deoxyhemoglobin can cause a large bias in the extracted saturation. This is particularly true for small concentrations of blood, where the calculated statistical uncertainty of the saturation is much smaller than the bias in the saturation introduced by a missing spectrally overlapping absorber (data not shown). Therefore, to be safe, we have in previous publications restricted the calculations of average blood saturation to spectra with blood volume fractions larger than 1%.^{11, 12, 13, 14} It is important to note that the saturation CI was calculated for our specific DPS measurement geometry with $400\text{-}\mu \mathrm{m}$ fibers, featuring a path length of $0.4\mathrm{mm}$ . Since the saturation CI will depend on signal attenuation rather than on blood volume fraction, the saturation CI (and potential bias) will be smaller for similar blood volume fractions in the case of DPS with larger fiber diameters, or any other measurement geometry with a longer path length.

In this work, the parameters were not constrained, i.e., all free parameters were allowed to range from minus infinity to plus infinity. However, in reality the parameters are physically constrained to $a_1⩾0$ , $0⩽a_2⩽-4$ , $a_3⩾0$ , $a_4⩾0\mu \mathrm{M}$ , $a_5⩾0\%$ , $0\%⩽a_6⩽100\%$ , and $a_7⩾5\mu \mathrm{m}$ (the diameter of a red blood cell). We did not implement these boundary conditions here to fully explore the sources and magnitudes of the statistical and systematic errors. When fitting an in vivo spectrum, however, it would be appropriate to constrain the parameters by implementing these physical boundary conditions to exclude nonphysical parameter values. However, in that case, care should be taken that the final value of any parameter is not at one of these artificially imposed limits. Boundary conditions on the individual parameters may also be artificially imposed by transformations of the variables. For example, by making the transformation $a_6^{\prime }\to 1∕2+\left[\arctan \left(a_6\right)\right]∕\pi$ , the transformed parameter $a_6^{\prime }$ is restricted to $0<a_6^{\prime }<1$ for $-\infty <a_6<\infty$ . In this way the fit is performed without boundary conditions on the fit parameters while still imposing physical boundary conditions on the transformed parameters. If the transformation of any fit parameter $a$ is given by $a^{\prime }\to f\left(a\right)$ , the statistical error of the transformed parameter is given by ${\sigma }^{\prime }=\sigma [df\left(a\right)∕da]$ .²⁰ For our example, $a_6^{\prime }\to 1∕2+\left[\arctan \left(a_6\right)\right]∕\pi$ we find ${\sigma }_{a6}^{\prime }={\sigma }_{a6}∕\pi (1+a_6^2)$ , where ${\sigma }_{a6}$ is the square root of the sixth diagonal element of the reduced chi-square normalized covariance matrix obtained from the fit.

In conclusion, we show that we can correctly estimate the statistical error on parameters extracted from fits to reflectance spectroscopic data. We demonstrate the use of various methods to identify and subsequently minimize the bias in the fitted parameters associated with systematic errors. Detailed analysis of the residual spectrum is potentially useful not only for minimization of bias, but also for monitoring changes in the in vivo microenvironment.