Linear regression models and K-means clustering for statistical analysis of fNIRS data

We propose a new algorithm, based on a linear regression model, to statistically estimate the hemodynamic activations in fNIRS data sets. The main concern guiding the algorithm development was the minimization of assumptions and approximations made on the dataset for the application of statistical tests. Further, we propose a K-means method to cluster fNIRS data (i.e. channels) as activated or not activated. The methods were validated both on simulated and in vivo fNIRS data. A time domain (TD) fNIRS technique was preferred because of its high performances in discriminating cortical activation and superficial physiological changes. However, the proposed method is also applicable to continuous wave or frequency domain fNIRS datasets. 2014 Optical Society of America OCIS codes: (000.5490) Probability theory, stochastic processes, and statistics (170.2655) Functional monitoring and imaging; (170.6920) Time-resolved imaging; (170.1470) Blood or tissue constituent monitoring. References and links 1. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” Neuroimage 63, 921–935 (2012). 2. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. Mata Pavia, U. Wolf and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” Neuroimage 85(1), 6-27 (2014). 3. T. J. Huppert, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “HomER: a review of time-series analysis methods for near-infrared spectroscopy of the brain,” Appl. Opt. 48(10), D280–D298 (2009). 4. P.H. Koh, D.E. Glaser, G. Flandin, B. Butterworth, A. Maki, D. Delpy and C.E. Elwell, “Functional Optical Signal Analysis (fOSA): A Software Tool for NIRS Data Processing Incorporating Statistical Parametric Mapping (SPM),” Journal of Biomedical Optics 12(6), 064010 (2007). 5. J. Ye, S. Tak, K. Jang, J. Jung and J. Jang, “NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy,” Neuroimage 44, 428–447 (2009). 6. H. Li, S. Tak and J.C. Ye, “Lipschitz–Killing curvature based expected Euler characteristics for p-value correction in fNIRS,” J. Neurosci. Methods 204(1), 61–67 (2012). 7. T. Fekete, D. Rubin, J. M. Carlson, and L. R. Mujica-Parodi, “The NIRS Analysis Package: noise reduction and statistical inference,” PLoS One 6, 24322 (2011). 8. http://www.biopac.com/fNIR-Software-Professional-Edition. 9. http://www.hitachi.co.jp/products/ot/analyze/kaiseki_en.html. 10. http://fnirs.org/software. 11. S. Tak and J.C. Ye, “Statistical analysis of fNIRS data: A comprehensive review,” Neuroimage 85(1), 72-91 (2014). 12. K. Friston, J. Ashburner, S. Kiebel, T. Nichols and W. Penny, “Statistical Parametric Mapping: The Analysis of Functional Brain Images: The Analysis of Functional Brain Images,” Academic Press (2011). 13. M.L. Schroeter, M.M. Bucheler, K. Muller, K. Uludag, H. Obrig, G. Lohmann, M. Tittgemeyer, A. Villringer and D.Y. von Cramon, “Towards a standard analysis for functional near-infrared imaging,” Neuroimage 21, 283–290 (2004). 14. J. Ye, S. Tak, K. Jang, J. Jung and J. Jang, “NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy,” Neuroimage 44, 428–447 (2009). 15. L. Gagnon, C. Gauthier, R. D. Hoge and F. Lesage, “Double-layer estimation of intraand extracerebral hemoglobin concentration with a time-resolved system,” J. Biomed. Opt. 13(5) 054019 (2008). 16. F. Scarpa, S. Brigadoi, S. Cutini, P. Scatturin, M. Zorzi, R. Dell'Acqua and G. Sparacino, “A reference-channel based methodology to improve estimation of event-related hemodynamic response from fNIRS measurements,” Neuroimage 72, 106-119 (2013). 17. G.H. Glover, “Deconvolution of impulse response in event-related BOLD fMRI,” NeuroImage 9(4), 416–429 (1999). 18. S. Prahl, “Optical absorption of hemoglobin,” (2013) http://omlc.ogi.edu/spectra/hemoglobin/index.html. 19. J.R. Mourant, T. Fuselier, J. Boyer, T.M. Johnson and I.J. Bigio, “ Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms,” Applied Optics 36(4), 949-957 (1997). 20. A. Torricelli, A. Pifferi, P. Taroni, E. Giambattistelli and R. Cubeddu, “In vivo optical characterization of human tissues from 610 to 1010 nm by time-resolved reflectance spectroscopy,” Physics in medicine and biology 46(8), 2227 (2001). 21. F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, “Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions and Software” (SPIE Press, 2010). 22. L. Zucchelli, D. Contini, R. Re, A. Torricelli and L. Spinelli, “Method for the discrimination of superficial and deep absorption variations by time domain fNIRS,” Biomedical Optics Express 4(12), 2893-2910 (2013). 23. D. Contini, L. Spinelli, M. Caffini, R. Cubeddu, and A. Torricelli, “A multichannel time-domain brain oximeter for clinical studies,” Proceedings of SPIE-OSA Biomedical Optics 7369 1D (2009). 24. E. Visani, L. Canafoglia, I. Gilioli, D.R. Sebastiano, V.E. Contarino, D. Duran, F. Panzica, R. Cubeddu, D. Contini, L. Zucchelli, L. Spinelli, M. Caffini, E. Molteni, A.M. Bianchi, S. Cerutti, S. Franceschett and A. Torricelli, “Hemodynamic and EEG Time-Courses During Unilateral Hand Movement in Patients with Cortical Myoclonus. An EEG-fMRI and EEG-TD-fNIRS Study,” Brain Topography (in press). 25. “R: A language and environment for statistical computing. R Foundation for Statistical Computing,” Vienna, Austria. http://www.R-project.org/. 26. A. Struyf, M. Hubert, P. Rousseeuw, “Clustering in an Object-Oriented Environment,” Journal of Statistical Software 1(4) 10630 (1997). 27. Terry Bo-Jau Kuo, Chang-Ming Chern, Wen-Yung Sheng, Wen-Jang Wong and Han-HwaHu, “Frequency Domain Analysis of Cerebral Blood Flow Velocity and Its Correlation With Arterial Blood Pressure,” J. of Cerebral Blood Flow and Metabolism 18, 311-318 (1998). 28. H. Obrig, M. Neufang, R. Wenzel, M. Kohl, J. Steinbrink, K. Einhaupl, and A. Villringer, “Spontaneous Low Frequency Oscillations of Cerebral Hemodynamics and Metabolism in Human Adults,” NeuroImage 12, 623– 639 (2000). 29. Y. Tong and B.D. Frederick, “Time lag dependent multimodal processing of concurrent fMRI and near-infrared spectroscopy (NIRS) data suggests a global circulatory origin for low-frequency oscillation signals in human brain,” NeuroImage 53, 553–564 (2010). 30. D. Canova, S. Roatta, D. Bosone and G. Micieli, “Inconsistent detection of changes in cerebral blood volume by near infrared spectroscopy in standard clinical tests,” J. Appl. Physiol. 110, 1646-1655 (2011). 31. L. Gagnon, M.A. Yücel, M. Dehaes, R.J. Cooper, K.L. Perdue, J. Selb, T.J. Huppert, R.D. Hoge and D.A. Boas, “Quantification of the cortical contribution to the NIRS signal over the motor cortex using concurrent NIRSfMRI measurements,” NeuroImage 59(4), 3933-3940 (2012). 32. T. Takahashi, Y. Takikawa, R. Kawagoe, S. Shibuya, T. Iwano and S. Kitazawa, “Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task,” NeuroImage 57, 991–1002 (2011). 33. B. Hallacoglu, A. Sassaroli, M. Wysocki, E. Guerrero-Berroa, M. Schnaider Beeri, V. Haroutunian, M. Shaul, I.H. Rosenberg, A.M. Troen and S. Fantini, “Absolute measurement of cerebral optical coefficients, hemoglobin concentration and oxygen saturation in old and young adults with near-infrared spectroscopy,” Journ. Biomed. Opt. 17(8), 081406 (2012).


Introduction
Functional near-infrared spectroscopy (fNIRS) is an optical technique able to noninvasively monitor the cerebral hemodynamic at cortical level [1,2].Exploiting the relatively low absorption of biological tissues, light in the red and near-infrared wavelength range can penetrate the human head down to some centimeters and reach the cerebral cortex.Therefore, fNIRS can provide a measure of oxy-and deoxy-hemoglobin (O2Hb and HHb, respectively), the main chromophores contributing to light absorption at this wavelength range.
Various techniques are proposed in literature to enhance and detect the fNIRS signal, irrespectively of the specific hardware implementation (e.g.continuous wave, frequency domain, or time domain regime) of the fNIRS technique.However, because of the relative novelty of the fNIRS technique, a common accepted framework for data analysis has not yet been agreed.The first public domain software package for fNIRS data analysis was Homer (acronym for Hemodynamic Evoked Response) [3].The software provided a graphical user interface and Matlab scripts for both the preprocessing and the standard statistics on fNIRS data.Homer has been upgraded and the new release Homer2 more easily supports group analyses and re-configuration of the processing stream, and it integrates users algorithms into the processing stream.Another free software is Functional Optical Signal Analysis (fOSA) [4], which offers Matlab based functions for a basic analysis of fNIRS data, incorporating several filters for signal denoising and providing also the Statistical Parametric Mapping (SPM) methodology for statistical analysis based on the general linear model (GLM) approach.More focused on the development of SPM routines is the non-commercial software NIRS-SPM [5].A novelty introduced by this program is represented by a voxel based alignment between interpolated maps instead of an inter-subjects realignment of optodes, in order to facilitate the group analysis.Further they proposed a new theory to deal with the multiple comparison problem for the p-value correction [6].Another software is NIRS analysis package (NAP) [7], which allows noise removal and GLM analysis, as well as anatomical registration of the measurements.fNIRSOFT is a stand-alone software to process, analyse and visualize fNIRS signals through a graphical user interface and/or scripting distributed by BIOPAC Systems, Inc. [8].Finally POTATo (Platform for Optical Topography Analysis Tools) is a software package for fNIRS signal processing and analysis, developed by Hitachi, Ltd. [9].A comprehensive list of software can be found in the website of the fNIRS Society [10].
The above mentioned tools share some common procedures for preprocessing the fNIRS data and for extracting the features of interest.Given that the functional studies are usually performed by repeating a particular task during several temporal slots, a preliminary way to inspect the fNIRS results is to detrend the signal by subtracting a mean value registered in a rest period before each task repetition and then to visualize the time series of O2Hb and HHb signals averaged over all the trials.The expected hemodynamic response is identified where the O2Hb concentration increases and simultaneously the HHb decreases.However the functional activity could be not easily detected because of the simultaneous presence of confounding effects like hemodynamic changes in the superficial layer (either systemic or task related), and movement artefacts.Therefore a careful data analysis and statistical inference must be considered to properly detect the signals related to a neuronal activation [11].
For statistical inference, the GLM approach has been adopted in most of the fNIRS software tools.GLM is a regression model assuming that a functional of a given signal can be modeled as a linear combination of known regressors, usually consisting of task-related boxcar functions.GLM was originally adopted for fMRI data analysis [12], and later used in fNIRS, taking advantage of the similarity between fMRI and fNIRS experiments in terms of designs and hypotheses [4,13].Anyway, because of the substantial differences between fMRI and fNIRS techniques, above all the lower spatial resolution of fNIRS due to the sparse optodes distribution over the head, there are a number of limitations, assumptions and specific issues that have to be considered when applying GLM on optical data [14].
In this paper we propose a new statistical method for the analysis of fNIRS data to statistically discriminate the hemodynamic activations registered by fNIRS.The proposed method is based on a linear regression model, and the main concern guiding the algorithm development was the minimization of assumptions and approximations made on the datasets for the application of statistical tests, such as the assumptions of independence, homoscedasticity and Gaussianity of residuals in order to guarantee the Gaussianity of the coefficients estimators.Furthermore, we propose a clustering algorithm aiming at better localizing the activated vs. not activated measured channels.The method is tested on simulated data mimicking real fNIRS measurements, and then applied on in-vivo measurements.In particular, we focused on time domain (TD) fNIRS data to have the best discrimination between cortical activation and superficial physiological changes.However, the proposed method is also suitable for continuous wave or frequency-domain fNIRS datasets.

Synthetic fNIRS data
A synthetic fNIRS dataset has been created to mimic real multichannel TD fNIRS measurements on a healthy adult during a motor task (handgrip experiment).A forward procedure and an inverse procedure were the main steps to obtain the synthetic dataset.
The forward procedure consisted in: 1) defining the geometry of the head; 2) assigning values for the hemodynamic parameters so as to calculate the hemodynamic response functions during a protocol; 3) converting hemodynamic parameters into absorption coefficients; 4) adding the information on the reduced scattering coefficients; 5) using a photon diffusion model to generate distributions of photon time-of-flight; 6) adding measurement noise to mimic a real TD fNIRS measurement.
The head has been modeled as a bilayered medium, where the upper layer is 1 cm thick, and the lower layer is ideally a semi-infinite medium.To a first approximation in fact this geometry can be used to simulate fNIRS measurements on the head of an adult, where an extra-cerebral layer (composed by scalp, skull and cerebrospinal fluid) overlays the intracerebral one (gray and white matter).
Hemoglobin concentrations were simulated by considering reference values of 12 μM for the O2Hb and 7 μM for the HHb in the superficial layer, and reference values of 30 μM for the O2Hb and 20 μM for the HHb in the lower layer [15].A dataset for 30 independent channels was generated, considering an optodes distribution over both central hemispheres (see Fig. 1).Two different experiments were simulated.The first experiment (EXP1) considers an ideal situation where a neuronal activation is generated in the cortical layer by a motor task, and no physiological oscillations occur in the superficial layer.This experiment consists of 10 repetitions (also called trials) of 10 s of baseline, 20 s of right handgrip movement and 10 s of recovery, for an overall duration of 400 s.The O2Hb and HHb concentrations in the upper layer were simulated to be constant at reference values during the whole experiment, while the concentrations in the lower layer were perturbed in some channels so as to mimic a hemodynamic response in correspondence to the task periods.This superimposed response profile was calculated as a convolution of a boxcar function, representing the task and rest alternation, with the Hemodynamic Response Function (HRF) evoked by a single stimulus.By following the method proposed by Scarpa et al. [16], the HRF was modeled as a linear combination of two different gamma-variant timedependent functions Γn: with: where α regulates the amplitude, τ1 and τ2 regulate the HRF shape, ρ1 and ρ2 tune its scale, while β determinates the undershoot.p coefficient was set to 5 as suggested by Glover et al. [17].Variability in amplitude of about 5% was considered among the different repetitions to account for possible differences in the execution of the task and/or in the functional response (e.g.habituation effects).The peak of the HRF for the O2Hb was chosen to be around 1555±75 nM; the HRF for the HHb was inverted and with a maximum set at -1/3 with respect to the O2Hb response.The free parameters have been chosen so as to create a HRF similar to one expected for the motor task of interest (α =1282, β =0.17, τ1=1, τ2=1, ρ1=-0.5, ρ2=3.5).To simulate an actual neuronal activation localized around the central positions of the hemisphere contralateral to the movement, some channels were considered activated with different intensities: channel number 16 (25% HRF), 17 (100% HRF), 18 (50% HRF), 21 (50% HRF), 28 (25% HRF), 29 (50% HRF); in the other channels no hemodynamic response was added.
In a second experiment (EXP2) the reference values for both layers and for the hemodynamic changes happening in the lower layer were identical to EXP1.However a physiological noise was added in the superficial layer by following the procedure reported by Scarpa et al. [16].An oscillation was built for O2Hb signal for one channel as a superposition of sinusoidal functions at different mean frequencies and amplitudes, as reported in Table 1: Amplitudes   and frequencies fi vary in the channel for every repetition in the range described in Table 1, while phases i  are equally distributed between 0 and 2π and are different for each trial.The generated signal was then replicated for all the channels by modifying the oscillation amplitude of a random value between ±10%.The HHb variations have been generated by threefold reducing the magnitude of the physiological noise simulated for the O2Hb.
The absorption coefficients at two wavelengths (690 and 820 nm) for both layers were computed from these hemoglobin concentration changes by exploiting the Lambert Beer law and the a priori knowledge of the specific absorption of O2Hb and HHb [18].
The scattering coefficients at the same wavelengths were derived from a simple approximation of the Mie theory [19]: by setting the amplitude scattering a and the power scattering b fixed respectively at 12 cm -1 and 0.5 for the upper layer, and at 12 cm -1 and 1 for the lower layer, for a reference λ0 at 660 nm [20].A forward model for photon diffusion in a bilayered geometry [21] was used generate synthetic time-resolved reflectance (TRR) curves for each channel by using as input parameters the optical properties and the source detector distance (fixed at 3 cm).A count rate of 5•10 5 ph/s was considered, the integration time was set at 1 s, and Poisson noise was added to the simulated curves to mimic real measurements.
The inverse procedure involved the following steps: 1) estimating the baseline optical properties and the absorption changes in the upper and lower layer; 2) calculating the hemodynamic parameters from the absorption coefficients.
The absolute values of μa and μs' have been recovered by TD fNIRS data by fitting the curves of the baseline period preceding each task period with a physical model for reflectance geometry in a homogeneous medium [21].Then absorption changes have been computed by means of the method proposed by Zucchelli et al. [22].The method allows the discrimination of superficial and deep absorption variations from TRR curves.It takes into account the effect of system set-up, as described by the Instrument Response Function (IRF) and the heterogeneous structure of the human head for the refined computation of the photon timedependent pathlength within each layer the tissue is composed of.It makes use of an approach based on time-gating of the photon distribution of time-of-flights.The fNIRS signal coming from the deep regions, more likely involved in the cerebral activity, can be cleaned from the superficial variations of the absorption properties, mainly due to systemic hemodynamics changes.
Fig. 2 and Fig. 3 show the time courses (folding average) for changes in O2Hb and HHb for the EXP2 in the upper layer and in the lower layer, respectively.

In vivo fNIRS data
The proposed method has been applied on in vivo data, acquired by a multi-channel dualwavelength TD fNIRS medical device developed by the Physics Department of the Politecnico of Milan [23] to preliminary validate its performance in real life settings.One right-handed healthy subject (male, 44 years old) underwent the experiment, consisting in a motor task (i.e.squeezing a soft ball in the right hand) at a rate of 2 Hz guided by a metronome.The same protocol simulated for synthetic data was maintained for in vivo experiment (10 repetitions of 10 s baseline, 20 s task and 10 s recovery, total duration 400 s).Instructions about the movement and rest were given by presenting a picture on a screen, which always had a fixation cross in the center.A total of fifteen detection bundles and eight light sources were positioned over the sensorimotor areas centered on C3 and C4 positions of the 10/20 standard system, following the configuration represented in Fig. 1.Pairs of light sources were sequentially illuminated in the left and right hemisphere every 0.25 s allowing the acquisition of 30 measurement points (channels) with an overall acquisition time of 1 s.Hemodynamic parameters were estimated by following the same steps previously described as "inverse procedure".The experiment was part of another study [24] that was reviewed and approved by the local ethics committee and it was conducted in compliance with the Declaration of Helsinki; the subject provided written informed consent to their participation in the study.
Folding average results for changes in O2Hb and HHb in the superficial and lower layer are shown in Fig. 4 and Fig. 5 respectively.

Statistical analysis
The proposed regression model consists in the following steps: 1) pre-processing; 2) estimate of the coefficients of the regressors for each trial; 3) inference test on the coefficients; 4) Kmeans algorithm for cluster analysis of activated channels.The method has been implemented in R Core Team (2014) [25].
A pre-processing algorithm is initially applied to the fNIRS datasets.The sample mean on the first 10 s of each trial is subtracted to the related trial, in order to detrend data.Then a smoothing spline algorithm is applied to the whole signal.If t y is hemoglobin concentration at time 1: 400 t  s, the algorithm calculates the curve   ŷt that minimize (in the class of twice differentiable functions) the following quantity: The first term represents the estimated Mean Squared Error (MSE) when using   summarizes the information about activation in each sub-interval.We consequently have a sample of size 10 of these quantities for each channel, on which we can make inference.This approach aims at limiting the most critical problems that are found in analysing fNIRS data through a linear (or a generalized linear) regression model: correlation, heteroscedasticity and not Gaussianity of residuals, that question the Gaussianity of the coefficients estimates.
If i indicates the sub-interval and k indicates the channel, we build for each channel 1: 30 k  the 10 following linear regression models: ,, , 1 , 2, 1 0, where , ik y is the vector of data, X is the design matrix, k  the vector of coefficients and , ik  a term of error.They are defined as follows: 1 Under the hypothesis that O2Hb increases during the task period, the regressors for O2Hb, rest and task, are obtained through a convolution between the HRF and a step-function equal to 0 in the first and last 10 s and 1 elsewhere for task (the opposite for rest) (see Fig. 6 left).On the contrary, given that HHb is expected to decrease during the task, the regressor for HHb is built as a convolution between the HRF and a step-function equal to 0 in the first and last 10 s, -1 elsewhere for task, the opposite for rest (see Fig. 6 right).For each channel k and sub-interval i, we then calculate the Ordinary Least Squares estimators for k  , as: Following this procedure we obtain fitted values , ˆ. ik y more similar to , .ik y than the ones found through a single linear regression model.An example is reported in Fig. 7, where fitted values and real data are plotted for regression models on trials (Fig. 7 left) and for a unique regression model (Fig. 7 right).In the first situation the Squared Error (SE) is lower than in the second one (34.1 and 42.4 μM 2 respectively), in which the obtained coefficients are the same for all trials, preventing any variability.For this reason the approach with regression models on trials well suits also experiments where the intensity of the HRF varies in time.
Moreover, through this approach, we have a population of

 
. In order to discriminate between activated/not-activated channels we focus on the contrast of the coefficients, ,, ˆî k i k task rest   , coherently with the current literature [12].We use the 10 found linear combinations ,,   and their Gaussianity to implement an inference test.We then build a map that shows the degree of activation of each channel.We conduct a hypothesis test for each channel k on the expected value of ,,   and we use the P-values to draw the map.
In particular, for fixed channel k, the test will be: A one-tailed test is chosen due to the shape of the regressors and the HRF.
Theoretically the decision of the test would be acceptance of the null hypothesis for the not-activated channels, and the rejection for the activated ones.In fact if a channel is activated we expect that the linear combination of the regressors is significant, and the coefficients related to it have expected value higher than 0. Due to the Gaussianity of , where X is a random variable from a t-student distribution with 1 n  degrees of freedom.We finally plot a map (named activation map, see Fig. 9 and following) in which the color of each channel is proportional to its P-value.Colors vary from white (activation) to black (no activation).In the following section we refer to this map as activation map.

K-means clustering algorithm
The localization of an activated area through a statistical analysis can be confirmed through K-means clustering algorithm.If k indicates the channel and i the trial, the following vector is considered for each channel:  10 .This clustering algorithm separates the 30 vectors in m groups, finding clusters that minimize the Euclidean distance within clusters and maximizes the one between clusters.We set m equal to 2, because we expect to observe two clusters: one with activated channels, one with not-activated ones.
The algorithm consists of 3 steps: 1. Initialization, in which the initial centers are randomly fixed.
2. The Euclidean distances from m centers are calculated for each vector, then vectors are assigned to the cluster with the nearest center.3. Updating of centers: the center of each cluster is calculated as the mean between the vectors belonging to the cluster.Steps 2) and 3) are repeated until convergence.
The choice of 2 m is confirmed also by the average silhouette width, a quality index allowing select the number of clusters [26].For each fixed m , the index varies from -1 to 1, increasing if the algorithm well classifies the vectors, decreasing if they are badly classified.An index is calculated for K-means with different number of clusters.Then K-means with the highest average silhouette width is chosen.If 2 m the index is equal to 0.85 for O2Hb, while for 2 7 m  we obtain indices lower than 0.40.Similar results are obtained for EXP2.This confirms our choice.

Synthetic fNIRS data
Activation maps for O2Hb and HHb are reported in Fig. 9 for EXP1 and in Fig. 10 for EXP2.
For dataset EXP1, obtained simulating ideal hemodynamic evolution without confounding oscillations, it is clear how the proposed method can discriminate the activated channels in the deeper region.P-values related to activated channels are in fact close to 0 (equal to 0 rounding to the third decimal place), creating a sharp division in the map between white and colored channels.At first sight, channel 30 might seem activated.Its P-value (P = 0.004) is however bigger (by more than a factor 10) than the highest activated-channel Pvalue (P = 0.0003 for channel 28).
If we consider as activated a channel with P-value lower than 0.001, we note that the map can discriminate exactly the activated channels.More specifically, the lowest P-value (less than 10 -7 ) is the one referred to channel 17, the most activated one (100% HRF), while channels with the lowest activation intensity (ch.16 and 28, 25% HRF) have the highest Pvalues among the activated ones (less than 0.0005).
In the upper layer we don't find any activation, as expected, and P-values are all higher than 0.01.Also in the activation map for the deeper layer of dataset EXP2 activated channels are instantly detectable.In the superficial layer P-values are very high (higher than 0.5) and very uniform, confirming that the simulated O2Hb is noise.
We can notice in the HHb activation maps the same trends found for O2Hb.Revealing activation in this situation is more difficult, because activation amplitude is lower than in O2Hb measures.Nevertheless, this method performs a good channels classification on these measures too: the only channel with a slightly high P-value is 28 (25% HRF).

In vivo data
The method proposed for activated channels detection was tested also on in vivo data.In vivo dataset was treated in the same way as synthetic ones: after a baseline subtraction and smoothing spline algorithm application, each trial was used for a linear regression model.Coefficients estimates were then tested in order to evaluate the dependence between data and regressors.Activation maps were calculated in the same way.This procedure produced good results on in vivo data, and it proved to be suitable in revealing activated channels.Activation maps for in vivo subject are reported in Fig. 11.In the superficial layer the subject shows high P-values in most of channels, while active areas can be clearly identified in the left hemisphere of deep layer.This proves the efficacy of the proposed procedure for activation detection, that can be successfully applied also on in vivo data.

K-means
K-means activation maps for synthetic data are reported in Fig. 12 for EXP1 and EXP2.In these images channels from different clusters are represented in different colors.Channels simulated as active (in deep layer) are circled by green.
The K-means algorithm is able to identify most of the activated channels for both O2Hb and HHb in the deep layer.The channels that present a higher activation intensity (channels 17, 18, 21, 29) are precisely clustered, for both EXP1 and EXP2, O2Hb and HHb measures.Conversely, channels 16 and 28, that present a low intensity of activation (25% HRF), are assigned to the cluster of not-activated ones.If there is no activation, as happens in the upper layer, the algorithm is not suitable.In fact it is forced to separate channels in two groups, even if vectors shouldn't be divided in two clusters.Thus results are unpredictable.Fig. 13 shows K-means maps for in vivo subject.The algorithm performs a classification coherent with the activation maps, always assigning the activated channels to the same cluster.In some situations, if there are one or more highly noisy channels, the algorithm separates them from the others.Thus highly noisy channels must be excluded from the K-means analysis, not to have clusters composed only by one or two anomalous channels.Excluded channels are reported in pink (their irregularity can be verified looking at folded averages in Fig. 4 and 5).

Discussion
In this paper a new method for analysis of fNIRS data has been introduced.It is based on linear regression models using, as regressors, convolutions between scale functions and HRF, as in current literature.
The main novelty is the splitting of data in trial or sub-intervals, each one representing a realization of the same event (the elementary activation sequence), and the application of a linear regression model on each of them.This way, an investigable sample of coefficients is obtained for each channel.In agreement with the literature we specifically focused on a linear combination of the coefficients.We verified the Gaussianity of the linear combinations through Shapiro-Wilks tests.Thus, rather than a unique linear combination for each channel, that is difficult to analyse because of heteroscedasticity, not-normality and correlation of residuals, we have a normal population of linear combinations for each channel, which can be easily investigated through an inference test.A one-tailed hypothesis test on the expected values of the linear combinations is built for each channel, expecting a rejection (acceptance) of the null hypothesis for activated (not-activated) channels.The related P-values are used to draw the activation maps.This procedure was succesfully applied on simulated data.
A further novelty is the use of a clustering algorithm (K-means) as a useful additional instrument in activation detection.K-means algorithm separates channels in two sharp groups: the information on the degree of activation of each channel is then lost.The output of the algorithm is a binary assignment of each channel, that is simply labelled as "active/notactive".Thus it can happen that low activated channels are assigned to the cluster of notactivated ones.This happens for example for channels 16 and 28 in Fig. 12, that present a low intensity of activation (25% HRF).The channels that present a higher activation intensity (channels 17, 18, 21, 29) are precisely clustered, for both EXP1 and EXP2, O2Hb and HHb measures.The clustering algorithm can be an important instrument of control: its right clustering is a further confirmation of the accuracy of the previous analysis and of the calculated activation maps.In fact it follows a different procedure compared to the linear regression model.It aims to detect the same channels applying another kind of analysis, that doesn't use statistical tests and is based different hypothesis.The clustering algorithm works well if there are channels with "similar" vectors (thus with a similar evolution in time.This is expected to happen with activated channels).If there is no activation the algorithm is forced to separate channels in two groups, thus results are unpredictable and with no sense.Some problems can also arise if one (or more) channels are particularly noisy and product vectors much different from the others: in this situation the K-means algorithm can be inaccurate, separating noisy channels from the others.For this reason highly noisy channels should be excluded before an analysis with K-means clustering algorithm.
A limitation of this study could be the accuracy of the simulated dataset.We should in fact observe that the creation of synthetic data well reproducing the superficial and deep O2Hb and HHb concentration changes happening during a succession of rest and task periods is not a trivial issue.As a matter of fact the actual behavior (particularly the exact magnitude and frequency components) of O2Hb and HHb physiological oscillations occurring concurrently within the superficial layers of the head and in the cerebral cortex, during rest or task periods, is never fully reported in literature with adequate precision, but often only for a range of frequencies, for partial regions of the head, or for physiological parameters different than hemoglobin species (blood flow, pressure or volume components) [27][28][29].Moreover, the quantitative definition of the amplitude variation within the brain of both species of hemoglobin following a neuronal activation, i.e. the hemodynamic response, is still an open issue in the scientific community.In fact reported fNIRS data present inaccuracies related to different factors: for instance most of the retrieved concentration values are obtained with CW fNIRS instruments, and thus present an intrinsic measurement error due to the poor depth resolution of the technique: the obtained cerebral signal is inevitably affected by extracerebral concentration variations [30][31][32].Further, the in vivo optical (absorption and scattering) properties of biological tissues are hardly measurable and data in literature present a high variability in the results [33].Finally, the different anatomical characteristics within and between subjects produce unavoidable analysis errors.Despite all these limitations the simulated datasets presented in this paper can be effectively used to test the performances of data analysis procedures.

Conclusion
The present study proposes a new procedure for statistical analysis of activated channels in fNIRS data.The introduced method minimizes the hypothesis made on data for the application of statistical tests.Thus it can be employed on a wide class of datasets, without losing validity even if high correlation and heteroscedasticity of residuals are proved or their Gaussianity is not verified.All these assumptions are relaxed through a model that provides a sample of activation-related quantities for each channel.The unique required hypothesis is Gaussianity of the activation-related quantities, and this hypothesis was always confirmed by statistical tests.This procedure was validated on a synthetic dataset and then on in vivo data from TD fNIRS, and it produced good results in both situations, detecting activated channels with precision.A clustering algorithm (K-means) is also proposed in the present study as an useful additional tool for activated channels detection.This clustering algorithm doesn't require any statistical hypothesis.It divides channels in two groups trough geometric considerations on activation-related quantities of each channel.Because of the different proceedings compared to the previous algorithm, K-means can be used as a reinforcing control instrument after the proposed method execution.To our knowledge the proposed method can complement the procedures contained in the most used software for fNIRS data analysis (e.g Homer2 and NIRS-SPM).However, a thorough comparison of the outcomes of different software tools is out of the scope of this paper.

Fig. 1 .
Fig. 1.Position of sources (red circles) and detectors (green circles) according to the 10/20 International System.The position of the 30 simulated measurement channels (purple diamonds) is also highlighted.
second term penalizes the curvature of   ŷt .The parameter λ controls the trade-off between the accuracy of   ŷt (for λ=0 it corresponds to the original data) and how it is smoothed.Thus data are estimated through the smoothing spline   ŷt that minimizes a weighted sum of MSE and the average curvature.Instead of using all the 400 s for a single linear regression model, we divide the time series of each channel in 10 sub-intervals (i.e.repetitions or trials) lasting 40 s (made of 10 s rest, 20 s task, and 10 s rest).We then individually apply a linear regression model to each repetition.Each trial is the elementary sequence revealing activation.It represents a realization of the same phenomenon.It can therefore be interesting to first analyse each trial independently, writing a regression model for each of them.Then we look for an appropriate quantity that

Fig. 7 .
Fig. 7. Comparison between O2Hb fitted values and real data (channel 18, EXP1) from regression models on sub-intervals (left) and from a unique regression model on all the 400 measurements (right).

Fig
Fig. 8. P-values of Gaussianity for , ˆik rest  is the expected value of the linear combination of the coefficients.

Fig. 9 .
Fig. 9. Statistical detection of channels' activation for O2Hb (left column) and HHb (right column) in the upper layer (top row) and lower layer (bottom row) for EXP1.The numbers inside the circles are the channel numbers while the numbers outside the circles are the P-values.Channels simulated as active are circled in green.

Fig. 10 .
Fig. 10.Statistical detection of channels' activation for O2Hb (left column) and HHb (right column) in the upper layer (top row) and lower layer (bottom row) for EXP2.The numbers inside the circles are the channel numbers while the numbers outside the circles are the P-values.Channels simulated as active are circled in green.

Fig. 11 .
Fig. 11.Statistical detection of channels' activation for O2Hb (left column) and HHb (right column) in the superficial layer (top row) and lower layer (bottom row) for in vivo data.The numbers inside the circles are the channel numbers while the numbers outside the circles are the P-values.

Fig. 12 .
Fig. 12. K-means maps for O2Hb (left column) and HHb (right column) in the superficial layer (top row) and lower layer (bottom row) for EXP1 (a) and EXP2 (b).Channels simulated as

Fig. 13 .
Fig. 13.K-means maps for O2Hb (left column) and HHb (right column) in the superficial layer (top row) and lower layer (bottom row) for in vivo subject.Channels simulated as active are circled in green.The numbers inside the circles are the channel numbers while colors identify the two clusters.Channels that have to be excluded are in pink.

Table 1 . Mean and standard deviation of frequency and amplitude of each physiological component a
a From Scarpa et al.Neuroimage 72, 106-119 (2013).