2.1.

Samples: Selection of Volunteers

For this study, volunteers with skin phenotypes I and II according to Fitzpatrick classification³³ were selected. Volunteers with presence or history of irritation, sensitivity to cosmetics, dermatological diseases, and those who used cream cosmetics over the past 48 h preceding this examination were excluded from the study. Prior approval from the Research Ethics Committee was obtained for this study (report number 132.812).

This study involves confocal Raman measurements of skin of 30 female participants segregated into three groups on the basis of age and presence of diabetes mellitus type II as young healthy (YH), elderly healthy (EH), and elderly diabetic (ED). YH contains the volunteers in the age group of 20 to 30 years, whereas EH and ED groups include the volunteers in the age group of 56 to 81 years. This article aims to determine the vibrational markers indicating intrinsic skin aging by comparing the sets of samples of young and elderly volunteers. The process of skin aging can be studied from the changes in the biochemical composition of the dermis layer of the skin, as it is mainly represented by the degradation of collagen present in this dermis layer.

The ED group was included in this study as diabetes mellitus causes an acceleration in the aging process due to the glycation process. Therefore, the changes due to intrinsic aging would be enhanced in this group. In this research group, there is an advanced glycation end products (AGEs) overload in skin tissue, degrading the framework formed by the amino acids of type I collagen such as proline (P) and hydroxyproline (HP).^34,35 These diabetic people have three times more Amadori products (precursors of AGEs) and twice more AGEs than normoglycemic people. Furthermore, independent authors use this group of volunteers to assess extreme conditions in aging of adipose tissue³⁶ and other diagnostics.³⁷ Therefore, in the volunteers with diabetes mellitus, the variability is expected to be high owing to the rapid rate of skin aging caused by the glycation process.

2.2.

Study Conditions

The instrumental technique described here is noninvasive and presents no risk of causing any damage to the volunteer. Before starting the skin measurements, the volunteers were acclimatized to $23\pm 2°\mathrm{C}$ temperature and $51\pm 5\%$ relative humidity for 60 min. After acclimatization, the volunteers were seated with the forearm positioned on an aluminum plate, which contains an optical window. Before collection of Raman spectra, the forearm sites used for the measurements were cleaned with cotton soaked in 1.0 mL of ethyl alcohol 97%.

2.3.

Instrumentation

The acquisition of confocal Raman spectra of human skin was performed using a confocal Raman system (River Diagnostics, Model 3510 skin composition analyzer) coupled to a laser with exciting radiation of wavelength 785 nm. The laser light is focused on the skin with a microscope lens ($40\times$) located under a quartz window. The Raman signal is collected by a charge coupled device (CCD) and recorded on a computer connected to this system. The calibration of wavelength was performed by using lines of a neon–argon lamp. The laser power utilized for skin analysis was 27 mW.

2.4.

Study Design

The measurements were performed on the forearm of the volunteers. The specific and selected regions of forearm were exposed to the laser light from the objective. The scattered laser light was collected by the same objective and sent to the spectrophotometer after being filtered. The data were collected in a spectral range 400 to $1800{\mathrm{cm}}^{-1}$ with spatial resolution of $2{\mathrm{cm}}^{-1}$ and two accumulations by spectra. The objective used is coupled to a piezoelectric system that controls the depth of the focal plane. Raman spectra were acquired from the surface of skin up to a depth of $122\mu \mathrm{m}$. Confocal Raman spectra were recorded in different tracks with varying step sizes. Steps of $2\mu \mathrm{m}$ for the first $20\text{-}\mu \mathrm{m}$ depth, then with a step size of $4\mu \mathrm{m}$ in the region from 20- to $60\text{-}\mu \mathrm{m}$ depth, and finally with a step size of $2\mu \mathrm{m}$ in the region from 60- to $122\text{-}\mu \mathrm{m}$ depth. In total, 52 spectra were recorded from each volunteer. The exposure time for the first two tracks, i.e., from surface to $20\mu \mathrm{m}$ and 20 to $60\mu \mathrm{m}$, was 10 s. For the higher depths greater than $60\mu \mathrm{m}$, the exposure time was increased to 60 s. This is due to the fact that the number of photons that reaches the detector will be reduced in these higher depths. In the preliminary tests, with this exposure time, spectra of the dermis region were obtained with good quality and high signal-to-noise ratio.

2.5.

Preprocessing of Data

To minimize the influence of noise, the spectra were smoothed by Labspec software (Horiba JobinYvon, France) using Savitzky–Golay filter (grade 2, size 3). Using the same software, fluorescence present in the data was eliminated by subtracting a baseline using the line segments between the wavenumbers: 400, 470, 552, 588, 633, 794, 798, 994, 1144, 1495, and $1717{\mathrm{cm}}^{-1}$. The specific noise coming from cosmic rays (CRs) was eliminated through a specific feature of this software. CRs occasionally affect CCD detectors, introducing large spikes with very narrow bandwidth in the spectrum and can occur randomly in the spectra at various wavenumber positions. These spikes should be removed in order to nullify its influence on the normalization and statistical analysis of data.³⁸ Spikes can be visually detected in the Raman spectrum because of its characteristic of having positive peaks with high intensity values with bandwidth much narrower than Raman bandwidth. Following this subtraction, the spectra were normalized by vector normalization, in which each Raman intensity was divided by square root of sum of the squared intensities calculated using the full spectrum.^39,40 After these procedures, the preprocessed spectra were subjected to statistical analysis.

3.1.

Cluster Analysis: Selecting the Correct Skin Layer

Before starting statistical analysis of the differences in confocal Raman spectra for the groups of samples, it should be ensured that the comparisons between the groups are made using the same layer of the skin.

Confocal Raman system acquires the spectra at different depths of the skin. Therefore, it is necessary to know which set of spectra corresponds to what layer of the skin. In order to identify and separate the different layers of the skin, an exploratory analysis by the HCA was performed on a single sample containing spectra varying in depth. This is necessary to select which skin depth will be used for the subsequent discussion of the study.

Cluster analysis is one of the important techniques for pattern recognition. It finds “natural” grouping (unsupervised) based on the similarity of the intensities of the dataset variables. It is applied when the spectra are to be interpreted as members of a category like skin depth. An excellent visualization of these clusters by similarity is dendrograms, also called the tree of clusters.

Figure 1 shows the dendrogram of the HCA for human skin according to the different depths obtained by confocal Raman equipment. By HCA, it is possible to classify the different layers of the skin into separate clusters. From these clusters, the set of spectra that corresponds to the skin layer to be monitored can be selected.

Fig. 1

Dendrogram obtained by HCA from the confocal Raman spectra of human skin at different depths.

In addition to the classification of spectra corresponding to various skin layers, cluster analysis is also useful to detect outliers. Forming a cluster with only one or two spectra or a distant cluster to others is a strong indication of anomalous data. In this case, the spectrum of this sample should be inspected to detect the possible error. Based on the inspection of the spectra, decision has to be made either to reprocess that spectrum or to repeat the experiment if possible or to reject that particular spectrum. There are some approaches for outlier detection such as Hotelling T, Mahalanobis distance, and Chi-square distribution.^39,41 However, it is highlighted that it is not acceptable to remove any sample simply by considering it as an outlier. Outliers can be legitimate observations and sometimes very interesting and important details can be revealed by them. Therefore, it is very important to investigate the nature of the outliers before making a decision on their exclusion.

From this study, the dermis region was chosen for further analysis as skin aging is mainly represented by the glycation and degradation of proteins like type I collagen present in the dermis region of the skin. Out of 52 spectra obtained per volunteer, spectra representing the dermis layer of each volunteer that were classified by cluster analysis were used for further analysis. The number of spectra of dermis region was found to be 30 per volunteer by cluster analysis. Average spectra calculated from these 30 spectra of dermis region for each individual were used in univariate and multivariate analyses.

3.2.

Correlation Matrix: Representativity and Reproducibility of the Study

An issue in academia is the representation of groups and volunteer numbers used for in vivo studies for the reproducibility of results. These issues must be addressed in studies using confocal Raman spectroscopy as these measurements can be questioned since each measurement is collected in a focal plane within specific surface area by the lens used. However, this issue of confocal Raman spectroscopy can be dealt with in many ways such as by assessing the changes in analytical information with depth and by determining the intragroup and intergroup correlations.

The Pearson correlation matrix provides an easy interpretation of correlation between two measurements. This matrix considers each spectrum as a variable and calculates the correlation coefficient on a scale ranging from $-1$ to $+1$. The data can be with high negative correlation ($-1$), with no correlation (zero), or with strong positive correlation ($+1$). This correlation coefficient is calculated from the mean and standard deviation of the observations given by the intensity at each wavenumber in the spectra.

Figure 2 shows an example in which the Pearson correlation matrix can be used. This array has arranged three spectra at different depths in the dermis of 70, 93, and $117\mu \mathrm{m}$ and compared them with the average spectrum calculated by taking the average of all the spectra that make up the layer of the dermis. It is very common to use the average of a set of spectra to plot a graph representing that set of spectra. The Pearson correlation array can be used to compare the correlation between two individual spectra as well as with the average of that set of spectra.

Fig. 2

Pearson correlation matrix. The range of the scale from high-negative correlation ($-1$ with red), uncorrelated (0 with white), to strong positive correlation ($+1$ with blue).

The analysis of Fig. 2 shows that the lowest correlation occurs for spectral comparison between the depths of 70 and $117\mu \mathrm{m}$ ($r=0.95$). This is expected because spectra are taken from the initial and terminal portions of the dermis, where there is high biochemical variability in the composition of the skin. However, since both these spectra belong to the basal layer, significant correlation has been observed even though there are certain differences due to variable skin composition. The same explanation applies for the comparison of the average of all the spectra from dermis with the individual spectra shown in Fig. 2. All the spectra show a strong positive correlation indicating that there is not much deviation from the average. Therefore, the average can be used to express the overall behavior of this layer.

3.3.

Univariate and Multivariate Statistical Approaches: Revealing Vibrational Markers

Once the adequate range of depth to analyze was assessed by HCA and correlation of individual spectra at different depths compared with the average spectrum used to represent the dermis by the Pearson correlation matrix, statistical data analysis starts the search for the spectral regions that can be used as an individual group with characteristics as relevant as possible with respect to the vibrational ways that discriminate the groups.

Figure 3(a) shows the average spectra of dermis of each volunteer (average between the spectra collected in the range from 60- to $122\text{-}\mu \mathrm{m}$ depth) of the total of 30 volunteers in this study, and Fig. 3(b) shows the average of the average spectra of the dermis region of the volunteers belonging to the three groups, namely YH, EH, and ED. Figure 3(b) also gives the information on the residual spectra obtained upon comparing YH, EH, and ED groups. From the visual inspection of the average confocal Raman spectra and the residual spectra in Figs. 3(a) and 3(b), it is quite evident that it is not possible to clearly distinguish a spectral region capable of discriminating the three groups of volunteers analyzed.

Fig. 3

Raman spectra of human skin: (a) average dermis spectra of 30 volunteers and (b) average dermis confocal Raman spectra of young healthy (YH), elderly healthy (EH), and elderly with diabetes mellitus (ED) with the residual spectra obtained for the comparisons: ED minus YH, ED minus EH, and EH minus YH.

Different statistical approaches can be used to seek the spectral regions that are significant enough to describe the differences among the groups analyzed. These spectral regions indicate the vibrational markers responsible for explaining skin aging. In this study, we describe the steps used by both univariate and multivariate statistical approaches to explore these vibrational markers.

These approaches can be used either alone or concomitantly as a result of their complementary character. In univariate statistical application on the set of spectra, each variable, i.e., each wavenumber, is analyzed separately. The sequence for this analysis is made variable by variable. In a multivariate approach, dimension reduction methods will be used and the spectral data are interpreted as a data matrix. In the case of nonsupervised analysis, the experiment corresponds to an array and by applying matrix algebra methods, the spectral regions that have greater variation in the data are determined, giving little importance to spectral overlap and noise regions. In the case of supervised analysis, only the regions correlated with a property of interest have more importance. Irrespective of the approach, all the wavenumbers that contribute are considered in a multivariate model with different weights.

Figure 4 represents the flowchart of steps in statistical data analysis. Univariate analysis is performed by applying a set of hypotheses tests. The hypotheses tests consist of a statistical procedure based on probability theory in which a parameter is tested on a set of values. Two hypotheses are considered: the null hypothesis $H_0$ (to be tested) and the alternative hypothesis $H_a$ that will be accepted if the group does not fall within the tested null hypothesis. The decision whether the test set belongs to $H_0$ or $H_a$ can be made through the analysis of $p$-values. It is important to note that for the correct application of hypothesis testing, the verification of the assumptions should be performed. If the assumptions are not checked properly, it may result in an incorrect inference about the dataset analyzed. For the univariate analysis, there is a great distinction between the sequence of processes to be employed if the data follow a normal distribution (parametric statistics) or not (nonparametric statistics).

Fig. 4

Flowchart of statistical data analysis.

The univariate approach begins with checking the normality of the data. If the data show the normal Gaussian distribution pattern, then it follows the path of parametric analysis; otherwise, it will follow nonparametric analysis. In the scheme shown in Fig. 4, Lilliefors test was used for the normality check.⁴² In this case, the default null hypothesis, which states that the sample follows normal distribution, is tested against the alternative hypothesis stating that it does not follow normal distribution. The nonparametric test used was the Kruskal–Wallis test.⁴³ This test compares the medians of the samples and gives the $p$-value for the hypothesis testing. Null hypothesis indicates that all the samples are drawn from the same population, whereas the alternative states that they are from different populations.

The test for the homogeneity of variances is the F-test.⁴⁴ An F-test of the null hypothesis checks that two independent samples coming from normal distributions are with the same variance against the alternative hypothesis, which states that they came from normal distributions with different variances. $t$-Test⁴⁴ was used to check if the average of the intensities of the tested groups is the same or different. $t$-Test of null hypothesis checks that the data in independent random samples from the normal distributions are with equal averages against the alternative hypothesis stating that the averages are not equal. The complete sequence shown in Fig. 4 can be implemented in free R software, which is presented in the Appendix and discussed in detail in various statistics books.^43,44

The $p$-value analysis for a given level of statistical significance is used to infer whether a variable is significant or not. This result is usually disposed in a table of $p$-values for each wavenumber. However, it is possible to make a $p$-value graph to facilitate the visualization of the results (as shown in Fig. 5).

Fig. 5

Results of univariate analysis for young healthy (YH), elderly healthy (EH), and elderly with diabetes mellitus (ED). $p\text{-}\text{values}<0.05$ are represented by black stars, and $p\text{-}\text{values}>0.05$ are represented by open circles. Horizontal line represents $p\text{-}\text{value}=0.05$.

Figure 5 shows the results by univariate statistical approach for the three comparisons (ED versus EH, ED versus YH, and EH versus YH). Black stars represent the variables with statistically significant differences at 95% confidence level. Therefore, $p\text{-}\text{values}<0.05$ are represented by black stars and $p\text{-}\text{values}>0.05$ are represented by open circles.

The region below the horizontal line ($p\text{-}\text{values}<0.05$) shows the candidates that can act as vibrational markers for each comparison. In order to be more specific, emphasis will be given to the spectral region between 800 and $1000{\mathrm{cm}}^{-1}$. This spectral region does not show any significant variables for comparing EH versus YH, but for comparing ED versus EH and ED versus YH, many variables appear as significant. In other words, the spectral region between 800 and $1000{\mathrm{cm}}^{-1}$ appears to be a good spectral marker to differentiate the group of ED to other groups, but it does not show significant difference in the comparison between the YH and EH groups. The vibrational modes in this spectral region 800 to $1000{\mathrm{cm}}^{-1}$ mainly correspond to the region of amino acids P and HP.^34,35

The considerations for a multivariate statistical approach start from the right side of the flowchart shown in Fig. 4. In case of spectral overlap [see Fig. 3(a)] or in case of very difficult analysis by visual inspection or in case of the presence of many vibrational modes, multivariate analysis provides good results by use of dimension reduction techniques.¹⁹ With the advantage of visualization of the behavior of the samples using a much smaller number of dimensions than spectral variables, this multivariate analysis is also useful in quantifying compounds. In multivariate analysis, we will discuss two different approaches: unsupervised analysis by PCA and supervised analysis by PLS-DA.

PCA is a method of dimension reduction, in which the new variables called principal components (PCs) are projected in hyperplanes in such a way that the dispersion between the original variables is maximum. This approach is extremely useful for the data whose matrices are in the order of tens to hundreds of columns and can therefore be visualized in diagrams in two or three dimensions.

The multivariate approach for the full Raman spectra of 30 samples after preprocessing steps is described in Sec. 2. The PCA model was fitted with three PCs with coefficient of explained variation of data matrix $R^{2}x=0.71$. The diagram of score plot for this PCA model is shown in Fig. 6(a). The PCA diagram shows evidence of the formation of groups of samples, which reinforces the importance of the proposed study as the PCA analysis is based only on the spectra of the samples without any additional information, i.e., the samples are grouped in accordance with their confocal Raman spectra profile. However, the PCA model is influenced by individual characteristics such as age, skin pigmentation, and any other variable that was not controlled, thus presenting a very good way to inspect the experiment and detect outliers. It is not mandatory to use unsupervised methods such as PCA before a supervised analysis. However, this step is recommended to detect the outliers where there is no provision of any detailed information.

Fig. 6

Multivariate statistical analysis: (a) nonsupervised PCA analysis and (b) supervised PLS-DA analysis for young healthy (YH), elderly healthy (EH), and elderly with diabetes mellitus (ED) groups.

Figure 6(a) shows the position of the samples in the space of PCs. PCA is an unsupervised method that is extremely useful in detecting outliers and natural groupings based solely on the information of spectra. In this figure, it is possible to identify the tendency to form clusters; however, this does not occur due to intrinsic variability of each volunteer. When there is no clear separation among the groups in the PCA, score and loading plot analysis cannot be performed as it does not explain any group behavior.

Since the goal of the study is to find characteristics of groups and not of individuals, supervised methods maximize the classification of samples owing to the characteristics of groups, reducing the effect of variables that contain the characteristics of individuals. A supervised model of PLS-DA was performed after the use of OSC on the sets of samples. For the PLS-DA approach, it is necessary to have some additional information along with the spectra of the sample like the origin of the samples or the group to which it belongs so that the algorithm tries to classify the samples into groups and consequently highlights the spectral regions responsible for their separation into groups.

Figure 6(b) shows the PLS-adjusted model with $R^{2}x=0.748$ and coefficient of variation explained of the discriminant class $R^{2}y=0.826$. In this model, the samples are grouped more efficiently than when compared with the PCA model, presenting complete separation of the groups. When this separation occurs, it is possible to ascertain which variables from the data matrix were uppermost, i.e., with higher weight for explaining the classification of the groups. These most important variables are then interpreted as the candidates for vibrational markers because they make the difference between these analyzed groups.

As the PLS-DA model was more robust and showed good results to discriminate the three groups of samples, the subsequent strategy was to perform comparisons of two groups with the PLS-DA approach and then to check which spectral regions of confocal Raman spectra were the most important variables to explain this separation. With this approach, it is possible to verify the contribution of each original variable and wavenumber to explain the separation of the data into groups (Fig. 7). The graph as shown in Fig. 7 was performed using a variable extraction method of variable importance in the projection (VIP),^45–47 where the higher the intensity of the peak, the greater the significance of that wavenumber in the separation of the groups. The list of VIP variables and the coefficient prediction capability $Q^2$, described in the Appendix, are obtained from rounds of cross-validation and have confidence intervals from Jack-knifing estimate. As the objective of this study is not to validate a new method, all other additional figures of merit that are generally used for various quantitative studies and comparisons with established techniques like limits of detection and quantification, root mean square, bias, sensibility and sensitivity, and curve ROC^48–50 were not discussed in this article.

Fig. 7

Graph of VIP for PLS-DA model of young healthy (YH), elderly healthy (EH), and elderly with diabetes mellitus (ED).

The results obtained by univariate and multivariate analyses can be complementary. This is clearly observed from the comparison of the results of univariate and multivariate analyses as shown in Figs. 5 and 7, respectively. In the present study, as per the results of PLS-DA model, spectral region between 800 and $1000{\mathrm{cm}}^{-1}$ on the graph appears to have important variables and this is the region that stands out most to explain the differences between ED and other groups in univariate analysis, thus corroborating and complementing the results obtained by univariate analysis.

One of the goals of this study is to disseminate the knowledge on proper application of various statistical tools to the academic community. For this purpose, a complete routine of analysis including all the statistical tests used in this study is presented in the Appendix using the free R software. R software is a language and an environment for statistical computing, providing a wide variety of statistical and graphical techniques. It runs on a wide variety of platforms like UNIX, FreeBSD, Linux, Windows, and MacOS.