- Email: [email protected]
A symmetrical subtraction combined with interpolated values for eliminating scattering from fluorescence EEM data
Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Contents lists available at ScienceDirect
Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy journal homepage: www.elsevier.com/locate/saa
A symmetrical subtraction combined with interpolated values for eliminating scattering from fluorescence EEM data Jing Xu a,⁎, Xiaofei Liu b, Yutian Wang a a b
Measurement Technology and Instrumentation Key Lab of Hebei Province, Yanshan University, Qinhuangdao, Hebei 066004, China Hebei Provincial Key Laboratory of Parallel Robot and Mechatronic System, Yanshan University, Qinhuangdao, Hebei 066004, China
a r t i c l e
i n f o
Article history: Received 19 January 2016 Received in revised form 29 March 2016 Accepted 3 April 2016 Available online 09 April 2016 Keywords: Symmetrical subtraction Interpolated values Combination Eliminating scattering Parallel factor analysis
a b s t r a c t Parallel factor analysis is a widely used method to extract qualitative and quantitative information of the analyte of interest from fluorescence emission-excitation matrix containing unknown components. Big amplitude of scattering will influence the results of parallel factor analysis. Many methods of eliminating scattering have been proposed. Each of these methods has its advantages and disadvantages. The combination of symmetrical subtraction and interpolated values has been discussed. The combination refers to both the combination of results and the combination of methods. Nine methods were used for comparison. The results show the combination of results can make a better concentration prediction for all the components. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Fluorescence spectroscopy is widely used for material analysis [1]. Parallel factor analysis (PARAFAC) method is widely used for many fields [2,3],it is a useful method with fully used tri-linear structure of fluorescence emission-excitation matrix (EEM). Scattering appearing in fluorescence EEM data will destroy the structure of fluorophore, bringing problematic spectra. Removing this kind of non-trilinear factor is necessary to improve the solution of PARAFAC. Existing methods include: interpolation on the band of scattering [4], fitting Gaussian curve to the scattering peak [5], modeling 1st Rayleigh scattering as a separate [6], insertion of zeros outside the data area [7], automatical identification of scattering using robust techniques [8], subtraction of a standard [9], MILES [10], using missing values on the band of scattering [11], using constraints [12], avoiding scattering area [13] and correction method [14]. Each method has its advantages, also with some disadvantages, but that is not a big deal, as special problem calling special method. For example, subtraction of a standard is a good method for eliminating Raman scattering, but sometimes a standard is not available especially for some real environment data sets. Avoiding scattering area can avoid non-trilinear effects but more chemical information will be lost. Insertion of interpolations will make spectra smoother, which is easy to fit, while ⁎ Corresponding author. E-mail address: [email protected] (J. Xu).
http://dx.doi.org/10.1016/j.saa.2016.04.006 1386-1425/© 2016 Elsevier B.V. All rights reserved.
interpolations are only the estimation of true spectra and sometimes they are not real, especially when a wrong interpolation area is chosen. Insertion of missing values would preserve the original shape of spectra, while too much missing values will bring problems to PARAFAC. A comparative study of two most common methods (missing values and interpolations) shows interpolations can improve converge speed, make the solution more meaningful and prevent local minima when a suitable blank is available [15].It seems combination of some of methods may bring some more advantages for the problem to be analyzed. A method of symmetrical subtraction combined with interpolated values for eliminating scattering from fluorescence EEM data is proposed in this study. Two kinds of combinations are applied in this study. One is to use each method individually and choose the best result to combine that is the combination of results. The other one is to combine symmetrical subtraction and interpolated values in one method that is symmetrical subtraction for suited scattering and interpolated values for other scattering. The method of symmetrical subtraction is based on the observation of the symmetry of scattering. When we can obtain half of the scattering band not overlapped with other signals, this part of scattering signal can be used for subtraction from the other half of its corresponding scattering band. It is often the case for first Rayleigh scattering. Other kind of scattering which cannot obtain nonoverlapping half scattering band can be removed using interpolations if necessary. Also the width of non-overlapping half scattering band can be used for the parameter for interpolations.
2
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
2. Theory Interpolation on the band of scattering has been introduced in detail in reference [4], so it will not introduce in detail here. Symmetrical subtraction is described in detail below. It is obvious to observe from the EEM spectra that the 1st order Rayleigh scattering is symmetric about emission wavelength which equals corresponding excitation wavelength at each emission spectrum. Due to the loss of energy, there is no emission will appear below its corresponding excitation. There the signals are only scattering and noise. While PARAFAC is a good method to extract tri-linear structure from the signal with noise, some large non-trilinear structure will bring some problematic solution. In most cases, the 1st order Rayleigh scattering is a bigger disturbance than other scattering because of its large amplitude. The appropriate diminution of 1st order Rayleigh scattering will improve the problem. Raman scattering also can influence the solution, but it is easy to remove by subtracting a solvent baseline from the emission spectra [5]. When a solvent baseline is unavailable, interpolations on the scattering area may be a good choice. When interpolations are used, a right width of scatter band is necessary. If the width is chosen too narrow, the scattering cannot be eliminated completely. Residual scattering will lead interpolations to fit scattering again and sometimes will bring a bigger fitting scattering. On the contrary, if the chosen width is too wide, a lot of chemical information will be lost. Also it has a chance to make wrong interpolations. The appropriate choice of the band width is the key point to apply the interpolation method. Unfortunately, the common used method to decide the width is though observational method. Because of the character of 1st order Rayleigh scattering, a symmetrical subtraction can be used. Even though it cannot eliminate 1st order Rayleigh scattering completely, it can reduce scattering to a great extent and the effect of residual scattering can be ignored to some extent. If other scattering band has nonoverlapping half scattering band, a symmetrical subtraction also can be used. If not, interpolations will be used. If the amplitude of Raman scattering is not so big to influence the results of PARAFAC, the symmetrical subtraction method can be used alone without combination method. There are some advantages of symmetrical subtraction method. First, it can make the spectra closer to the real shape. Second, compared the missing values, it has a faster speed and more stable. Third, if the stokes shift is short, some useful information will overlap with 1st order Rayleigh scattering. Interpolations will lose this information, while subtraction can retain them as much as possible. The method of combination of symmetrical subtraction and interpolated values in one method has advantages above three and another one. Forth it can reduce the number of input parameter, because at least the parameters for the width of 1st order Rayleigh scattering are needless. Also if a completely eliminating 1st order Rayleigh scattering is needed, the width of nonoverlapping half scattering band can be used as the input parameter for interpolations. The steps of the symmetrical subtraction method
interpolation function with the known values in the nonoverlapping half scattering band. Note because there are not always emission wavelengths which equal excitation wavelengths and not every symmetric emission wavelength is measured in measured EEM, interpolation is necessary. If we can obtain fluorescence values at every symmetric emission wavelength exactly, these values can be used for subtraction directly. Also, the best situation of applying cubic spline interpolation function is all the symmetric points within range of EEM, that is to say extrapolation should be avoided. When extrapolation must be used, please use accurate range of half scattering band for subtraction or insert zeros into the half band of interpolation without measured values for a better result. In practice, modest extrapolation won’t affect the results too much, which will show in experiment section. 3. The fluorescence values on the symmetric positions are subtracted from the values on the corresponding point Emo. If the new values on Emo less than zero, make them equal zeros. 4. Make the fluorescence values on non-overlapping half scattering band equal zeros. The steps of the combination of symmetrical subtraction and interpolated values in one method are shown below. Step1 to step4 are the same as above, and step5 is shown below 5. If it is necessary to remove scattering band without non-overlapping half scattering band, scattering band is replaced by interpolation as described in reference [4]. 3. Experimental 3.1. Data set The data set is made of fifteen mixtures of five fluorophores. They are catechol, hydroquinone, indole, tryptophane and tyrosine. Every sample in the data set contains two to four fluorophores. The fluorophores have highly overlapping spectra. The excitation range is 230–320 nm with 5 nm interval. The emission range is 230– 500 nm with 2 nm interval. The data set get from http://www. models.kvl.dk/source [6]. The information of the data set in more detail is explained by Rinnan and Andersen, which can be found in reference of reference [6]. A sample in the data set is shown in Fig. 1. 1st Rayleigh scattering, 1st Raman scattering and 2nd Rayleigh scattering can be observed obviously. Compared to Rayleigh scattering,
1. According to the corresponding excitation wavelength Exi, determine the length of emission wavelength greater than Exi. These points lie in the length of emission wavelength called the Emo. The length needn't very accurate. The only request is the length need longer than the half width of 1st order Rayleigh scattering because the signals below 1st order Rayleigh scattering are zeros. Note when other scattering with non-overlapping half scattering band uses symmetrical subtraction, the length of emission wavelength needed subtraction should be no longer than the distance of nonoverlapping area. 2. According to symmetry, calculate the symmetric positions of emission wavelength which equals Exi. Calculate the fluorescence values on the symmetric positions using the cubic spline
Fig. 1. The contour of No.16 sample in data set which contains four components.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
the intensity of Raman scattering is low. 1st Rayleigh scattering is overlapped with chemical components. This data set can be used to show the effect of removing overlapping scattering and nonoverlapping scattering and the effect of ignoring little residual scattering to PARAFAC modeling. First the data set is reduced as it was handled in demo.m in flumod_ver1.0_2005Feb (the demo file for reference [6]). The last 25 emission wavelengths were reduced to remove high-emission points because flumod.m does not currently handle those. The reason of the data set handled as the demo file is for the comparison of scores and spectra. Then the whole emission wavelengths were used to show the result for removing 1st Rayleigh scattering, 1st Raman scattering and 2nd Rayleigh scattering. 3.2. Software Matlab (The MathWorks, Natick, MA) R2014a is used. The algorithms used for interpolation get from www.models.kvl.dk [4]. The algorithms used for PARAFAC get from N-way Toolbox which also can be got from www.models.kvl.dk [16]. 4. Results and discussion Nine methods of data processing for comparison are shown below. Method1: handle reduced data set without eliminating scattering using PARAFAC with non-negative constrains applied on all three modes. Method2: handle reduced data set with modeling 1st order Rayleigh scattering as a separate component provided by reference [6]. Method3: handle reduced data set with eliminating 1st order Rayleigh scattering by symmetrical subtraction. The range of interpolation or subtraction is 10 points. Then use PARAFAC with nonnegative constrains applied on all three modes. Method4: handle reduced data set with eliminating 1st order Rayleigh scattering by symmetrical subtraction and eliminating 1st order Raman scattering by interpolation. The range of interpolation or subtraction is 10 points. The parameters of the width of 1st order Raman scattering are [10 10]. Then use PARAFAC with non-negative constrains applied on all three modes. Method5: handle reduced data set with eliminating 1st order Rayleigh scattering and 1st order Raman scattering by interpolation. The parameters of the widths of scatterings are [10 10], [10 10], respectively. Then use PARAFAC with non-negative constrains applied on all three modes. Method6: handle whole data set without eliminating scattering using PARAFAC with non-negative constrains applied on all three modes. Method7: handle whole data set with eliminating 1st and 2nd order Rayleigh scattering by symmetrical subtraction. The range of interpolation or subtraction is 4 points. Then use PARAFAC with nonnegative constrains applied on all three modes. Method8: handle whole data set with eliminating 1st and 2nd order Rayleigh scattering by symmetrical subtraction and eliminating 1st order Raman scattering by interpolation. The range of interpolation or subtraction is 4 points. The parameters of the width of 1st order Raman scattering are [10 10]. Then use PARAFAC with non-negative constrains applied on all three modes. Method9: handle whole data set with eliminating 1st and 2nd order Rayleigh scattering and 1st order Raman scattering by interpolation. The parameters of the widths of scatterings are [10 10], [10 10], [10
3
10], respectively. Then use PARAFAC with non-negative constrains applied on all three modes. Figs. 1 to 9 (in Appendix A) shows the plots of scores and spectra using above nine methods. The figures in Appendix A show all nine methods shows the similar results with minor differences in low emission wavelengths of reduced data set and both low and high wavelengths of whole data set. From the emission loadings it is easy to see all the methods of eliminating scattering can decrease scattering in different degrees. The method using symmetrical subtraction can decrease 1st order and 2nd Rayleigh scattering to the degree of 1st order Raman scattering. In this study, some extrapolation is applied when the excitation wavelength is low where there are not enough points for complete symmetrical subtraction, which produce some negative subtracted from the symmetrical peaks to make a larger amplitude, but that doesn’t influence the result too much which can be drawn from the figures in Fig. 3 and Fig. 4 in Appendix A. For more perfect results, a more exact interpolation range of half scattering band for subtraction should be used to avoid too much extrapolation (Fig. 7 and Fig. 8). Also, for this data set 1st Rayleigh scattering corresponding to the lowest excitation wavelength and 2nd Rayleigh scattering corresponding to the highest excitation are not eliminated, which makes the little residual in Fig. 7 (in Appendix A). For better results, these two kinds of residual can eliminating by other methods, such as subtraction by half symmetrical band of adjacent emission spectrum or setting these data to zero or missing data because most of time these values far away from useful information which can be seen from the Fig. 1. Method2, Method5 and Method9 can make corresponding scattering area smoother due to precise modeling scattering by times of iterations (for Method2) and the properties of the interpolation function (for Method5 and Method9). It is not the main purpose to get smooth spectra because it is easy to see from the spectra which parts are the scatterings. Most of the time, we use PARAFAC to obtain the information of qualitative and quantitative. Little residual scattering will not influence the qualitative analysis. We still easy to obtain the information of the analytes we need. The purpose we need to remove scatterings before PARAFAC because they will influence the algorithm of convergence to stable and accurate solutions. All the figures in Appendix A show the scattering will not influence the general shapes of spectra, that is, the qualitative analysis is feasible, so some quantitative indexes are compared. A method of removing scattering that can speed up PARAFAC has practical application significance, so the run time of functions and the number of iterations are compared. The results of run time and number of iterations of PARAFAC of nine methods are shown in Table 1. In Table 1 t1 means the run time of eliminating scattering. t2 means the run time of function parafac. t means the sum of t1 and t2 that is the total time of obtaining the final spectra of components. However, for method2, t1 means the run time of function fluwid (a function for estimating the width of Rayleigh peak) and t2 means the run
Table 1 The run time and number of iterations of PARAFAC of eight methods. t1 Method1 Method2 Method3 Method4 Method5 Method6 Method7 Method8 Method9
0.1230 s 0.1300 s 0.2730 s 0.1270 s 0.2060 s 0.2970 s 0.1400 s
t2
t
it
6.5260 s 32.9190 s 4.4400 s 4.1670 s 4.5230 s 6.2350 s 4.8590 s 5.6000 s 4.9200 s
6.5260 s 33.0420 s 4.5700 s 4.4400 s 4.6500 s 6.2350 s 5.0650 s 5.8970 s 5.0600 s
59 73 65 64 60 72 63 59
4
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Table 2 The sum of squared residuals of the model and concentrations.
Method1 Method2 Method3 Method4 Method5 Method6 Method7 Method8 Method9
Model
Component 1
Component 2
Component 3
Component 4
Component 5
608750 51060 109070 73805 64388 852930 93432 104940 69631
3.3478e−13 1.8576e−13 1.9886e−13 1.8912e−13 1.5970e−13 4.2127e−13 2.0166e-13 1.8898e-13 1.6117e−13
3.2485e−12 5.0681e−13 3.3389e−13 3.0836e−13 2.8415e−13 4.1706e−12 3.3954e-13 3.2234e-13 2.8693e−13
3.0442e−12 2.5005e−12 2.8955e-12 2.8348e−12 2.8630e−12 3.2065e−12 2.8841e-12 2.8106e-12 2.8679e−12
7.0051e−11 4.8974e−11 3.9426e−11 9.8311e−11 1.0076e−10 7.4465e−11 3.8785e-11 9.0075e-11 1.0439e−10
1.3247e−12 9.3308e−13 8.8540e−13 9.9235e−13 9.7026e−13 2.2970e−12 8.7337e-13 9.7550e-13 9.7586e−13
time of function flumod (a function for modeling 1st order Rayleigh scattering as a separate component and decomposing an EEM). Parameter ‘it’ means the number of iterations of PARAFAC to achieve convergence, which is an assistant analysis of t2. For method2, it is an iterative process of PARAFAC, so the number of iterations of PARAFAC is not given. The sum of squared residuals of the model and the data set together with the sum of squared residuals of the predicted concentrations and real concentrations are shown in Table 2. Table 1 shows Method2 has the longest time to obtain the solution because it is an iterative process of PARAFAC that is several times of function parafac needed to be done during the process. Even though it is a litter time-consuming, the advantage of this method is modeling the scattering. Comparing t1 of Method3, 4, 5, 7, 8, 9, the more scattering to be removed by symmetrical subtraction, the more time the algorithm needs, because more symmetrical points and subtraction need to be calculated. Fortunately, the time of symmetrical subtraction is short. Conversely, no matter how many scattering need to be removed, the time of interpolation are almost the same. This is mainly due to the programming of algorithms. Except Method2, t2 of other methods are in the same order of magnitude. Due to random initial values, there are minor differences in run time of PARAFAC. The methods with eliminating scattering take less time. Note the data set used in the sum of squared residuals of the model and data set is the data set with different pretreatment of removing scattering that is they have different standard data sets, so direct comparison of these values makes little sense to evaluate model, which can only show PARAFAC can model smooth data set better and will ignore non-trilinear structure. For above reasons, the sums of squared residuals of the predicted concentrations and real concentrations are used. Table 2 shows PARAFAC without eliminating scattering has the maximum sum of squared residuals of concentrations among component1 to component3 and component5, which suggests scatterings will influence the models and they should be removed. Method3 reduces the residual of all the components comparing to Method1 which shows symmetrical subtraction can remove scattering effectively. Method7 has the minimum sum of squared residuals of concentrations among component4 to component5 due to a more accurate range of half scattering band for interpolation and subtraction. The combination of symmetrical subtraction and interpolation in one method obtains a better result of component3 than either of the single method. Table 2 shows methods containing interpolation of Raman scattering make better predicted concentrations of component1 to component3. The reason might be Raman scattering go through the spectra of these spectra and influence the modeling of these components. Also methods containing interpolation of Raman scattering make worse predicted concentrations of component4 and component5. The reason might be interpolation make a poor estimation of these two components and more information
was lost using interpolation. Because the amplitude of 1st order Raman scattering is not big compared to the signal of components, keep Raman scattering in data set will not influence the result seriously. The combination of symmetrical subtraction and interpolated values in one method can improve some estimation of some components, but they are not all better than the best estimations in single method, because the combination bring advantage and disadvantage of both methods. From Table 2, we can draw a conclusion that different methods of removing scattering have different predictive accuracies of different components. No one method is best for all different kinds of components in all situations. The effect of different methods depends on the position of the scattering and components, the scattering amplitude and other factors. The most suitable method depends on the type of the components of interest and their environment. If all the components need higher predictive accuracies, the combination of the results of different method can get better results. For this data set, using symmetrical subtraction for component4 and component5 combined with interpolated values for component1 and component2 and Method2 for component3 can have a better result of all the components than using either of the methods. Also, Table 2 shows all the methods of eliminating scattering proposed here are useful. 5. Conclusion Nine methods are used to compare the results of PARAFAC after applied different methods of eliminating scattering. The results show no one method can have all best estimation of predicted concentrations. The combination of the results can have the best results of all components or only choose the best method for just some components of interest. The combination of symmetrical subtraction and interpolated values in one method can average the results of both methods. This comparison is only suit this kind of data set and other different data set need to compare again to find suitable method of eliminating scattering. When a single method of eliminating scattering cannot get the best results for all components of interest, a combination method may be a good choice. Both combination of results and combination of methods can be tried to obtain a good result. Acknowledgments This project is supported by the National Natural Science Foundation of China (Grant No. 61471312) and the National Natural Science Foundation of Hebei (Grant No. F2015203072, Grant No. C2014203212) and Graduate Student Innovation Fund Project of Hebei. Thanks the authors of the references for the data set and program and the platform of www.models.kvl.dk for providing much useful information.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Appendix A
Fig. 1. The results of Method1. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
5
6
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 2. The results of Method2. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 3. The results of Method3. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
7
8
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 4. The results of Method4. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 5. The results of Method5. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
9
10
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 6. The results of Method6. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
11
Fig. 7. The results of Method7. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
12
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 8. The results of Method8. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
13
Fig. 9. The results of Method9. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
14
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
References [1] K. Kumar, A.K. Mishra, Chemom. Intell. Lab. Syst. 116 (2012) 78–86. [2] C.H. Santos, G. Nicolodelli, R.A. Romano, et al., J. Braz. Chem. Soc. 26 (6) (2015) 1136–1142. [3] F. Alarcon, M.E. Baez, M. Bravoetc, Talanta 103 (2013) 361–370. [4] M. Bahram, et al., J. Chemom. 20 (2006) 99–105. [5] C.I. David, et al., Anal. Biochem. 139 (1984) 316–318. [6] A. Rinnan, K. Booksh, R. Bro, Anal. Chim. Acta 537 (2005) 349–358. [7] A. Rinnan, C.M. Andersen, Chemom. Intell. Lab. Syst. 76 (2005) 91–99. [8] S. Engelen, Chemom. Intell. Lab. Syst. 86 (2007) 35–51.
[9] D.M. McKnight, E.W. Boyer, P.K. Westerhoff, P.T. Doran, T. Kulbe, D.T. Andersen, Limnol. Oceanogr. 46 (2001) 38–48. [10] R. Bro, N.D. Sidiropoulos, A.K. Smilde, J. Chemom. 16 (2002) 387–400. [11] J. Christensen, V.T. Povlsen, J. Sørensen, J. Dairy Sci. 86 (2003) 1101–1107. [12] C.M. Andersen, R. Bro, J. Chemom. 17 (2003) 200–215. [13] R. Bro, Chemom. Intell. Lab. Syst. 46 (1999) 133–147. [14] R.G. Zeppa, W.M. Sheldonb, M.A. Moran, Mar. Chem. 89 (2004) 15–36. [15] S. Elcoroaristizabal, et al., Chemom. Intell. Lab. Syst. 142 (2015) 124–130. [16] C.A. Andersson, R. Bro Chemom, Intell. Lab. Syst. 52 (1) (2000) 1–4.
Contents lists available at ScienceDirect
Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy journal homepage: www.elsevier.com/locate/saa
A symmetrical subtraction combined with interpolated values for eliminating scattering from fluorescence EEM data Jing Xu a,⁎, Xiaofei Liu b, Yutian Wang a a b
Measurement Technology and Instrumentation Key Lab of Hebei Province, Yanshan University, Qinhuangdao, Hebei 066004, China Hebei Provincial Key Laboratory of Parallel Robot and Mechatronic System, Yanshan University, Qinhuangdao, Hebei 066004, China
a r t i c l e
i n f o
Article history: Received 19 January 2016 Received in revised form 29 March 2016 Accepted 3 April 2016 Available online 09 April 2016 Keywords: Symmetrical subtraction Interpolated values Combination Eliminating scattering Parallel factor analysis
a b s t r a c t Parallel factor analysis is a widely used method to extract qualitative and quantitative information of the analyte of interest from fluorescence emission-excitation matrix containing unknown components. Big amplitude of scattering will influence the results of parallel factor analysis. Many methods of eliminating scattering have been proposed. Each of these methods has its advantages and disadvantages. The combination of symmetrical subtraction and interpolated values has been discussed. The combination refers to both the combination of results and the combination of methods. Nine methods were used for comparison. The results show the combination of results can make a better concentration prediction for all the components. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Fluorescence spectroscopy is widely used for material analysis [1]. Parallel factor analysis (PARAFAC) method is widely used for many fields [2,3],it is a useful method with fully used tri-linear structure of fluorescence emission-excitation matrix (EEM). Scattering appearing in fluorescence EEM data will destroy the structure of fluorophore, bringing problematic spectra. Removing this kind of non-trilinear factor is necessary to improve the solution of PARAFAC. Existing methods include: interpolation on the band of scattering [4], fitting Gaussian curve to the scattering peak [5], modeling 1st Rayleigh scattering as a separate [6], insertion of zeros outside the data area [7], automatical identification of scattering using robust techniques [8], subtraction of a standard [9], MILES [10], using missing values on the band of scattering [11], using constraints [12], avoiding scattering area [13] and correction method [14]. Each method has its advantages, also with some disadvantages, but that is not a big deal, as special problem calling special method. For example, subtraction of a standard is a good method for eliminating Raman scattering, but sometimes a standard is not available especially for some real environment data sets. Avoiding scattering area can avoid non-trilinear effects but more chemical information will be lost. Insertion of interpolations will make spectra smoother, which is easy to fit, while ⁎ Corresponding author. E-mail address: [email protected] (J. Xu).
http://dx.doi.org/10.1016/j.saa.2016.04.006 1386-1425/© 2016 Elsevier B.V. All rights reserved.
interpolations are only the estimation of true spectra and sometimes they are not real, especially when a wrong interpolation area is chosen. Insertion of missing values would preserve the original shape of spectra, while too much missing values will bring problems to PARAFAC. A comparative study of two most common methods (missing values and interpolations) shows interpolations can improve converge speed, make the solution more meaningful and prevent local minima when a suitable blank is available [15].It seems combination of some of methods may bring some more advantages for the problem to be analyzed. A method of symmetrical subtraction combined with interpolated values for eliminating scattering from fluorescence EEM data is proposed in this study. Two kinds of combinations are applied in this study. One is to use each method individually and choose the best result to combine that is the combination of results. The other one is to combine symmetrical subtraction and interpolated values in one method that is symmetrical subtraction for suited scattering and interpolated values for other scattering. The method of symmetrical subtraction is based on the observation of the symmetry of scattering. When we can obtain half of the scattering band not overlapped with other signals, this part of scattering signal can be used for subtraction from the other half of its corresponding scattering band. It is often the case for first Rayleigh scattering. Other kind of scattering which cannot obtain nonoverlapping half scattering band can be removed using interpolations if necessary. Also the width of non-overlapping half scattering band can be used for the parameter for interpolations.
2
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
2. Theory Interpolation on the band of scattering has been introduced in detail in reference [4], so it will not introduce in detail here. Symmetrical subtraction is described in detail below. It is obvious to observe from the EEM spectra that the 1st order Rayleigh scattering is symmetric about emission wavelength which equals corresponding excitation wavelength at each emission spectrum. Due to the loss of energy, there is no emission will appear below its corresponding excitation. There the signals are only scattering and noise. While PARAFAC is a good method to extract tri-linear structure from the signal with noise, some large non-trilinear structure will bring some problematic solution. In most cases, the 1st order Rayleigh scattering is a bigger disturbance than other scattering because of its large amplitude. The appropriate diminution of 1st order Rayleigh scattering will improve the problem. Raman scattering also can influence the solution, but it is easy to remove by subtracting a solvent baseline from the emission spectra [5]. When a solvent baseline is unavailable, interpolations on the scattering area may be a good choice. When interpolations are used, a right width of scatter band is necessary. If the width is chosen too narrow, the scattering cannot be eliminated completely. Residual scattering will lead interpolations to fit scattering again and sometimes will bring a bigger fitting scattering. On the contrary, if the chosen width is too wide, a lot of chemical information will be lost. Also it has a chance to make wrong interpolations. The appropriate choice of the band width is the key point to apply the interpolation method. Unfortunately, the common used method to decide the width is though observational method. Because of the character of 1st order Rayleigh scattering, a symmetrical subtraction can be used. Even though it cannot eliminate 1st order Rayleigh scattering completely, it can reduce scattering to a great extent and the effect of residual scattering can be ignored to some extent. If other scattering band has nonoverlapping half scattering band, a symmetrical subtraction also can be used. If not, interpolations will be used. If the amplitude of Raman scattering is not so big to influence the results of PARAFAC, the symmetrical subtraction method can be used alone without combination method. There are some advantages of symmetrical subtraction method. First, it can make the spectra closer to the real shape. Second, compared the missing values, it has a faster speed and more stable. Third, if the stokes shift is short, some useful information will overlap with 1st order Rayleigh scattering. Interpolations will lose this information, while subtraction can retain them as much as possible. The method of combination of symmetrical subtraction and interpolated values in one method has advantages above three and another one. Forth it can reduce the number of input parameter, because at least the parameters for the width of 1st order Rayleigh scattering are needless. Also if a completely eliminating 1st order Rayleigh scattering is needed, the width of nonoverlapping half scattering band can be used as the input parameter for interpolations. The steps of the symmetrical subtraction method
interpolation function with the known values in the nonoverlapping half scattering band. Note because there are not always emission wavelengths which equal excitation wavelengths and not every symmetric emission wavelength is measured in measured EEM, interpolation is necessary. If we can obtain fluorescence values at every symmetric emission wavelength exactly, these values can be used for subtraction directly. Also, the best situation of applying cubic spline interpolation function is all the symmetric points within range of EEM, that is to say extrapolation should be avoided. When extrapolation must be used, please use accurate range of half scattering band for subtraction or insert zeros into the half band of interpolation without measured values for a better result. In practice, modest extrapolation won’t affect the results too much, which will show in experiment section. 3. The fluorescence values on the symmetric positions are subtracted from the values on the corresponding point Emo. If the new values on Emo less than zero, make them equal zeros. 4. Make the fluorescence values on non-overlapping half scattering band equal zeros. The steps of the combination of symmetrical subtraction and interpolated values in one method are shown below. Step1 to step4 are the same as above, and step5 is shown below 5. If it is necessary to remove scattering band without non-overlapping half scattering band, scattering band is replaced by interpolation as described in reference [4]. 3. Experimental 3.1. Data set The data set is made of fifteen mixtures of five fluorophores. They are catechol, hydroquinone, indole, tryptophane and tyrosine. Every sample in the data set contains two to four fluorophores. The fluorophores have highly overlapping spectra. The excitation range is 230–320 nm with 5 nm interval. The emission range is 230– 500 nm with 2 nm interval. The data set get from http://www. models.kvl.dk/source [6]. The information of the data set in more detail is explained by Rinnan and Andersen, which can be found in reference of reference [6]. A sample in the data set is shown in Fig. 1. 1st Rayleigh scattering, 1st Raman scattering and 2nd Rayleigh scattering can be observed obviously. Compared to Rayleigh scattering,
1. According to the corresponding excitation wavelength Exi, determine the length of emission wavelength greater than Exi. These points lie in the length of emission wavelength called the Emo. The length needn't very accurate. The only request is the length need longer than the half width of 1st order Rayleigh scattering because the signals below 1st order Rayleigh scattering are zeros. Note when other scattering with non-overlapping half scattering band uses symmetrical subtraction, the length of emission wavelength needed subtraction should be no longer than the distance of nonoverlapping area. 2. According to symmetry, calculate the symmetric positions of emission wavelength which equals Exi. Calculate the fluorescence values on the symmetric positions using the cubic spline
Fig. 1. The contour of No.16 sample in data set which contains four components.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
the intensity of Raman scattering is low. 1st Rayleigh scattering is overlapped with chemical components. This data set can be used to show the effect of removing overlapping scattering and nonoverlapping scattering and the effect of ignoring little residual scattering to PARAFAC modeling. First the data set is reduced as it was handled in demo.m in flumod_ver1.0_2005Feb (the demo file for reference [6]). The last 25 emission wavelengths were reduced to remove high-emission points because flumod.m does not currently handle those. The reason of the data set handled as the demo file is for the comparison of scores and spectra. Then the whole emission wavelengths were used to show the result for removing 1st Rayleigh scattering, 1st Raman scattering and 2nd Rayleigh scattering. 3.2. Software Matlab (The MathWorks, Natick, MA) R2014a is used. The algorithms used for interpolation get from www.models.kvl.dk [4]. The algorithms used for PARAFAC get from N-way Toolbox which also can be got from www.models.kvl.dk [16]. 4. Results and discussion Nine methods of data processing for comparison are shown below. Method1: handle reduced data set without eliminating scattering using PARAFAC with non-negative constrains applied on all three modes. Method2: handle reduced data set with modeling 1st order Rayleigh scattering as a separate component provided by reference [6]. Method3: handle reduced data set with eliminating 1st order Rayleigh scattering by symmetrical subtraction. The range of interpolation or subtraction is 10 points. Then use PARAFAC with nonnegative constrains applied on all three modes. Method4: handle reduced data set with eliminating 1st order Rayleigh scattering by symmetrical subtraction and eliminating 1st order Raman scattering by interpolation. The range of interpolation or subtraction is 10 points. The parameters of the width of 1st order Raman scattering are [10 10]. Then use PARAFAC with non-negative constrains applied on all three modes. Method5: handle reduced data set with eliminating 1st order Rayleigh scattering and 1st order Raman scattering by interpolation. The parameters of the widths of scatterings are [10 10], [10 10], respectively. Then use PARAFAC with non-negative constrains applied on all three modes. Method6: handle whole data set without eliminating scattering using PARAFAC with non-negative constrains applied on all three modes. Method7: handle whole data set with eliminating 1st and 2nd order Rayleigh scattering by symmetrical subtraction. The range of interpolation or subtraction is 4 points. Then use PARAFAC with nonnegative constrains applied on all three modes. Method8: handle whole data set with eliminating 1st and 2nd order Rayleigh scattering by symmetrical subtraction and eliminating 1st order Raman scattering by interpolation. The range of interpolation or subtraction is 4 points. The parameters of the width of 1st order Raman scattering are [10 10]. Then use PARAFAC with non-negative constrains applied on all three modes. Method9: handle whole data set with eliminating 1st and 2nd order Rayleigh scattering and 1st order Raman scattering by interpolation. The parameters of the widths of scatterings are [10 10], [10 10], [10
3
10], respectively. Then use PARAFAC with non-negative constrains applied on all three modes. Figs. 1 to 9 (in Appendix A) shows the plots of scores and spectra using above nine methods. The figures in Appendix A show all nine methods shows the similar results with minor differences in low emission wavelengths of reduced data set and both low and high wavelengths of whole data set. From the emission loadings it is easy to see all the methods of eliminating scattering can decrease scattering in different degrees. The method using symmetrical subtraction can decrease 1st order and 2nd Rayleigh scattering to the degree of 1st order Raman scattering. In this study, some extrapolation is applied when the excitation wavelength is low where there are not enough points for complete symmetrical subtraction, which produce some negative subtracted from the symmetrical peaks to make a larger amplitude, but that doesn’t influence the result too much which can be drawn from the figures in Fig. 3 and Fig. 4 in Appendix A. For more perfect results, a more exact interpolation range of half scattering band for subtraction should be used to avoid too much extrapolation (Fig. 7 and Fig. 8). Also, for this data set 1st Rayleigh scattering corresponding to the lowest excitation wavelength and 2nd Rayleigh scattering corresponding to the highest excitation are not eliminated, which makes the little residual in Fig. 7 (in Appendix A). For better results, these two kinds of residual can eliminating by other methods, such as subtraction by half symmetrical band of adjacent emission spectrum or setting these data to zero or missing data because most of time these values far away from useful information which can be seen from the Fig. 1. Method2, Method5 and Method9 can make corresponding scattering area smoother due to precise modeling scattering by times of iterations (for Method2) and the properties of the interpolation function (for Method5 and Method9). It is not the main purpose to get smooth spectra because it is easy to see from the spectra which parts are the scatterings. Most of the time, we use PARAFAC to obtain the information of qualitative and quantitative. Little residual scattering will not influence the qualitative analysis. We still easy to obtain the information of the analytes we need. The purpose we need to remove scatterings before PARAFAC because they will influence the algorithm of convergence to stable and accurate solutions. All the figures in Appendix A show the scattering will not influence the general shapes of spectra, that is, the qualitative analysis is feasible, so some quantitative indexes are compared. A method of removing scattering that can speed up PARAFAC has practical application significance, so the run time of functions and the number of iterations are compared. The results of run time and number of iterations of PARAFAC of nine methods are shown in Table 1. In Table 1 t1 means the run time of eliminating scattering. t2 means the run time of function parafac. t means the sum of t1 and t2 that is the total time of obtaining the final spectra of components. However, for method2, t1 means the run time of function fluwid (a function for estimating the width of Rayleigh peak) and t2 means the run
Table 1 The run time and number of iterations of PARAFAC of eight methods. t1 Method1 Method2 Method3 Method4 Method5 Method6 Method7 Method8 Method9
0.1230 s 0.1300 s 0.2730 s 0.1270 s 0.2060 s 0.2970 s 0.1400 s
t2
t
it
6.5260 s 32.9190 s 4.4400 s 4.1670 s 4.5230 s 6.2350 s 4.8590 s 5.6000 s 4.9200 s
6.5260 s 33.0420 s 4.5700 s 4.4400 s 4.6500 s 6.2350 s 5.0650 s 5.8970 s 5.0600 s
59 73 65 64 60 72 63 59
4
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Table 2 The sum of squared residuals of the model and concentrations.
Method1 Method2 Method3 Method4 Method5 Method6 Method7 Method8 Method9
Model
Component 1
Component 2
Component 3
Component 4
Component 5
608750 51060 109070 73805 64388 852930 93432 104940 69631
3.3478e−13 1.8576e−13 1.9886e−13 1.8912e−13 1.5970e−13 4.2127e−13 2.0166e-13 1.8898e-13 1.6117e−13
3.2485e−12 5.0681e−13 3.3389e−13 3.0836e−13 2.8415e−13 4.1706e−12 3.3954e-13 3.2234e-13 2.8693e−13
3.0442e−12 2.5005e−12 2.8955e-12 2.8348e−12 2.8630e−12 3.2065e−12 2.8841e-12 2.8106e-12 2.8679e−12
7.0051e−11 4.8974e−11 3.9426e−11 9.8311e−11 1.0076e−10 7.4465e−11 3.8785e-11 9.0075e-11 1.0439e−10
1.3247e−12 9.3308e−13 8.8540e−13 9.9235e−13 9.7026e−13 2.2970e−12 8.7337e-13 9.7550e-13 9.7586e−13
time of function flumod (a function for modeling 1st order Rayleigh scattering as a separate component and decomposing an EEM). Parameter ‘it’ means the number of iterations of PARAFAC to achieve convergence, which is an assistant analysis of t2. For method2, it is an iterative process of PARAFAC, so the number of iterations of PARAFAC is not given. The sum of squared residuals of the model and the data set together with the sum of squared residuals of the predicted concentrations and real concentrations are shown in Table 2. Table 1 shows Method2 has the longest time to obtain the solution because it is an iterative process of PARAFAC that is several times of function parafac needed to be done during the process. Even though it is a litter time-consuming, the advantage of this method is modeling the scattering. Comparing t1 of Method3, 4, 5, 7, 8, 9, the more scattering to be removed by symmetrical subtraction, the more time the algorithm needs, because more symmetrical points and subtraction need to be calculated. Fortunately, the time of symmetrical subtraction is short. Conversely, no matter how many scattering need to be removed, the time of interpolation are almost the same. This is mainly due to the programming of algorithms. Except Method2, t2 of other methods are in the same order of magnitude. Due to random initial values, there are minor differences in run time of PARAFAC. The methods with eliminating scattering take less time. Note the data set used in the sum of squared residuals of the model and data set is the data set with different pretreatment of removing scattering that is they have different standard data sets, so direct comparison of these values makes little sense to evaluate model, which can only show PARAFAC can model smooth data set better and will ignore non-trilinear structure. For above reasons, the sums of squared residuals of the predicted concentrations and real concentrations are used. Table 2 shows PARAFAC without eliminating scattering has the maximum sum of squared residuals of concentrations among component1 to component3 and component5, which suggests scatterings will influence the models and they should be removed. Method3 reduces the residual of all the components comparing to Method1 which shows symmetrical subtraction can remove scattering effectively. Method7 has the minimum sum of squared residuals of concentrations among component4 to component5 due to a more accurate range of half scattering band for interpolation and subtraction. The combination of symmetrical subtraction and interpolation in one method obtains a better result of component3 than either of the single method. Table 2 shows methods containing interpolation of Raman scattering make better predicted concentrations of component1 to component3. The reason might be Raman scattering go through the spectra of these spectra and influence the modeling of these components. Also methods containing interpolation of Raman scattering make worse predicted concentrations of component4 and component5. The reason might be interpolation make a poor estimation of these two components and more information
was lost using interpolation. Because the amplitude of 1st order Raman scattering is not big compared to the signal of components, keep Raman scattering in data set will not influence the result seriously. The combination of symmetrical subtraction and interpolated values in one method can improve some estimation of some components, but they are not all better than the best estimations in single method, because the combination bring advantage and disadvantage of both methods. From Table 2, we can draw a conclusion that different methods of removing scattering have different predictive accuracies of different components. No one method is best for all different kinds of components in all situations. The effect of different methods depends on the position of the scattering and components, the scattering amplitude and other factors. The most suitable method depends on the type of the components of interest and their environment. If all the components need higher predictive accuracies, the combination of the results of different method can get better results. For this data set, using symmetrical subtraction for component4 and component5 combined with interpolated values for component1 and component2 and Method2 for component3 can have a better result of all the components than using either of the methods. Also, Table 2 shows all the methods of eliminating scattering proposed here are useful. 5. Conclusion Nine methods are used to compare the results of PARAFAC after applied different methods of eliminating scattering. The results show no one method can have all best estimation of predicted concentrations. The combination of the results can have the best results of all components or only choose the best method for just some components of interest. The combination of symmetrical subtraction and interpolated values in one method can average the results of both methods. This comparison is only suit this kind of data set and other different data set need to compare again to find suitable method of eliminating scattering. When a single method of eliminating scattering cannot get the best results for all components of interest, a combination method may be a good choice. Both combination of results and combination of methods can be tried to obtain a good result. Acknowledgments This project is supported by the National Natural Science Foundation of China (Grant No. 61471312) and the National Natural Science Foundation of Hebei (Grant No. F2015203072, Grant No. C2014203212) and Graduate Student Innovation Fund Project of Hebei. Thanks the authors of the references for the data set and program and the platform of www.models.kvl.dk for providing much useful information.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Appendix A
Fig. 1. The results of Method1. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
5
6
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 2. The results of Method2. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 3. The results of Method3. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
7
8
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 4. The results of Method4. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 5. The results of Method5. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of emission loadings, and (d) excitation loadings.
9
10
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 6. The results of Method6. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
11
Fig. 7. The results of Method7. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
12
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
Fig. 8. The results of Method8. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
13
Fig. 9. The results of Method9. (a) Plot of sample scores, (b) emission loadings, (c) partial enlarged plot of low emission loadings, (d) partial enlarged plot of high emission loadings, and (e) excitation loadings.
14
J. Xu et al. / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 165 (2016) 1–14
References [1] K. Kumar, A.K. Mishra, Chemom. Intell. Lab. Syst. 116 (2012) 78–86. [2] C.H. Santos, G. Nicolodelli, R.A. Romano, et al., J. Braz. Chem. Soc. 26 (6) (2015) 1136–1142. [3] F. Alarcon, M.E. Baez, M. Bravoetc, Talanta 103 (2013) 361–370. [4] M. Bahram, et al., J. Chemom. 20 (2006) 99–105. [5] C.I. David, et al., Anal. Biochem. 139 (1984) 316–318. [6] A. Rinnan, K. Booksh, R. Bro, Anal. Chim. Acta 537 (2005) 349–358. [7] A. Rinnan, C.M. Andersen, Chemom. Intell. Lab. Syst. 76 (2005) 91–99. [8] S. Engelen, Chemom. Intell. Lab. Syst. 86 (2007) 35–51.
[9] D.M. McKnight, E.W. Boyer, P.K. Westerhoff, P.T. Doran, T. Kulbe, D.T. Andersen, Limnol. Oceanogr. 46 (2001) 38–48. [10] R. Bro, N.D. Sidiropoulos, A.K. Smilde, J. Chemom. 16 (2002) 387–400. [11] J. Christensen, V.T. Povlsen, J. Sørensen, J. Dairy Sci. 86 (2003) 1101–1107. [12] C.M. Andersen, R. Bro, J. Chemom. 17 (2003) 200–215. [13] R. Bro, Chemom. Intell. Lab. Syst. 46 (1999) 133–147. [14] R.G. Zeppa, W.M. Sheldonb, M.A. Moran, Mar. Chem. 89 (2004) 15–36. [15] S. Elcoroaristizabal, et al., Chemom. Intell. Lab. Syst. 142 (2015) 124–130. [16] C.A. Andersson, R. Bro Chemom, Intell. Lab. Syst. 52 (1) (2000) 1–4.