2.1.1.

SIM data

Semisimulated data were generated based on resting-state measurements from $N=30$ individuals that originated from an intermittent NFB study (unpublished). One participant was excluded because the rest data was missing. For the remaining 29 subjects, $\sim 3\min$ of fNIRS data were recorded while subjects were looking at a black fixation cross on a gray background displayed on a computer screen. After applying the channel-based exclusion criteria of the present study (cf. Sec. 2.3.1 paragraph channel pruning for details), $N=23$ older adults (11 females, 12 males; age [mean ± SD]: $[65.52\pm 4.13]$ years; range: 58 to 71 years) were included for the final analyses.

2.1.2.

ME and MI data

The real data set comprised ME and MI data. Details of the study design can be found in Klein et al.²⁸ In short, $N=32$ healthy older adults performed five motor tasks, i.e., a self-paced sequential finger tapping task that was physically performed with either the left (ME LEFT) or the right (ME RIGHT) hand or kinesthetically imagined with either the left or the right hand. The fifth task was the kinesthetic imagination of an arbitrarily chosen whole-body movement involving both arms and legs. Each task consisted of 12 trials and the trial structure was a blocked design with 15 s of task period preceded and followed by 18 to 22 s of rest period. The order of tasks was pseudorandomized, with the restriction that not all of the three MI tasks appeared consecutively. Furthermore, a 15 s break after six trials of the same task and between two tasks was included in which subjects were allowed to slightly move their head, shoulders, and/or limbs if they felt that this helps for relaxation. All three MI tasks were combined into one MI task set (36 trials). After applying study specific exclusion criteria (cf. Sec. 2.3.1 paragraph channel pruning), $N=24$ subjects (13 females, 11 males; age [mean ± SD]: $[63.63\pm 4.73]$ years; range: 56 to 71 years) remained for analysis.

Participants of both studies gave written informed consent and were paid 10 Euros/h as reimbursement. General inclusion criteria for both studies were right-handedness, $>4\mathrm{h}$ of sleep, and no history of neurological, psychiatric, or psychological diseases (cf. Klein et al.,²⁸ for more details and all exclusion criteria). Both studies were approved by either the regular (resting-state data) or the medical (real data) Ethics Committee of the University of Oldenburg.

2.2.1.

Additional unused data

Optode positions and fiducial points were digitized with an optical digitizer (Xensor, ANT Neuro, The Netherlands). Because most of the digitizer data of the SIM data set was corrupted due to missing or inaccurate optode locations, digitzer data was not further considered. Moreover, the data sets contain electromyography data that were used in a previous publication²⁸ but are of no relevance for the present study.

2.3.1.

fNIRS analysis

FNIRS data processing was conducted with a combination of the NIRS Brain AnalyzIR toolbox³³ ( https://github.com/huppertt/nirs-toolbox) and custom made scripts. For a better overview, in Fig. 1 the workflow from data types (A, E) over the preprocessing pipeline (B), the applied correction methods (F) to the generated quality measures (G) is visualized.

Quality check and channel pruning. Channel quality was assessed by means of the qt-nirs toolbox that uses the scalp coupling index (SCI^34,35) and the peak spectral power (PSP³⁵) for evaluating signal quality (SCI threshold = 0.6, $Q$ threshold = 0.65, and PSP threshold = 0.1; https://github.com/lpollonini/qt-nirs). Note that the GLM-based correction methods are using the AR-IRLS model³⁶ (cf. Sec. 2.3.2), which is relatively robust against channels with poor signal quality and thus these channels do not necessarily need to be removed. However, as the quality assessment is important for other correction methods it was however performed as a general step in the processing pipeline. The first inclusion criterion was that all short-distance channels remained for analyses, which was the case for 25 participants of the SIM data set and for 24 of the ME/MI data set. For the SIM data, the second inclusion criterion concerned the consistent number of identical channels for adding the simulated HRF across subjects. Here, the aim was to find the same two channels across subjects that remain within an ROI after channel pruning. This was the case for channels 6 and 8 of ROI M1 LEFT in $N=23$ subjects [cf. Figs. 1(d) and 1(e)]. Across the final sample of the SIM data, on average (±SD) 2.26 (± 2.56; ranging from 0 to 8 channels) regular channels were pruned. Pruned channels were replaced by NaNs in the respective data set.

For the ME/MI data, the second inclusion criterion was that in each of the three ROIs at least two channels remained for analysis after pruning. Based on this, all N = 24 subjects were considered for analysis and on average across subjects 2.71 of all regular channels (SD± 3.20; ranging from 0 to 9 channels) were pruned. Data of the pruned channels were replaced by NaNs in the respective data set. A visualization of the frequency of the remaining channels of both data sets after pruning can be found in Fig. 1(c) and exact numbers are documented in Tables S21–S26 in the Supplementary Material.

Preprocessing. Preprocessing was applied to all channels, including both short and regular-distance channels. Raw data were transformed into optical density changes and were then corrected for possible motion artifacts by means of the temporal derivative distribution repair (TDDR) method.³⁷ Afterwards, data were converted into hemoglobin concentration changes by means of the modified Beer–Lambert law (mBLL). Therefore, the individual age-related³⁸ differential pathlength factor (DPF) was used in order to calculate the partial pathlength factor [PPF = DPF × partial volume factor (PVF = 1/60)], which was used for the mBLL. To generate the SIM data, the canonical HRF was hereupon added to the rest data [see following paragraph Generation of semisimulated (SIM) data for details]. Then, data were first low-pass filtered with a zero-phase digital Butterworth filter (cut-off = 0.09 Hz; order = 2) and then high-pass filtered with a zero-phase digital Butterworth filter (cut-off = 0.01 Hz; order = 2). The cut-off frequencies of [0.01, 0.09] Hz were selected based on the recommendations by Pinti et al.;³⁹ however, departing from this, in the present data analysis, no finite impulse response (FIR) filter with a high filter order ($>1000$) was applied because the SIM data had no sufficient length for this filter order ($\sim 3\min$ of rest data) and we wanted to keep the preprocessing the same across data sets.

Generation of semisimulated (SIM) data. Resting-state data of the first data set were used to generate the SIM data. Canonical HRFs (peak at 6 s) with different amplitudes were generated and added to the rest data of neighboring channels 6 and 8, located in the ROI M1 LEFT [cf. Fig. 1(c)]. Semisimulated data were generated for two channels and with different amplitudes to simulate a spatial pattern with a signal strength gradient (cf. Sec. 2.4.3). The amplitude for channel 6 was $10\mu \mathrm{M}$ for $\mathrm{\Delta }[\mathrm{HbO}]$ and $10/3\mu \mathrm{M}$ for $\mathrm{\Delta }[\mathrm{HbR}]$ (100% amplitude) and $5\mu \mathrm{M}$ for $\mathrm{\Delta }[\mathrm{HbO}]$ and $5/3\mu \mathrm{M}$ for $\mathrm{\Delta }[\mathrm{HbR}]$ (50% amplitude) for channel 8. For each subject, a total of five trials were simulated as a blocked design with 15 s of task period preceded and followed by 18 to 22 s of rest period.

2.3.2.

Systemic activity correction (SAC) methods

No systemic activity correction (NO SAC). Data after the band-pass filtered preprocessing stage served as a reference measure and are referred to as no systemic activity corrected (NO SAC) data.

Common average reference (CAR). A very simple correction method without the use of SDCs is the common average reference (CAR) spatial filter, which has its origin in EEG analysis.²² Based on the idea that a global signal is present in all channels, a spatial filter can be generated by simply subtracting the average time signal across all $N$ channels from each single channel ${\mathbf{x}}_i$, $(i=1,\dots ,N)$ to reduce the global signal influence, resulting in the corrected channels ${\mathbf{x}}_{{\mathrm{CAR}}_i}$ [cf. Eq. (1)]²²

Global Component Removal (GCR). The second correction method without using SDCs is the global component removal (GCR) method.^23,24 The GCR is a more advanced spatial filter based on Gaussian spatial filtering and singular value decomposition (SVD). The Gaussian kernel smoother $G$ is a two-dimensional kernel smoother that is applied to remove detail and noise from a data set.²³ It requires the Montreal Neurological Institute (MNI) coordinates of each channel that are stored in a $N\times N$ distance matrix $D$

Short Separation Regression (SSR). A simple and easy to implement SAC method that involves SDCs is the short separation regression (SSR).^11,19 By simply applying Eq. (7), it was shown that SA can be reduced

Here, $\alpha$ is defined as the quotient of two scalar products ($⟨·,·⟩$)

In this method, ${\mathbf{y}}_{\mathrm{SDC}}$ can be either the closest SDC, the SDC showing the strongest correlation with the regular channel ${\mathbf{y}}_{\mathrm{RC}}$ or an average across all available SDCs. In the present study ${\mathbf{y}}_{\mathrm{SDC}}$ represents the signal of the closest SDC to each regular channel determined by the Euclidian distance.

General Linear Model with all SDCs as Regressors (GLM ALL). The most recommended SAC method when using SDCs is to incorporate the SDCs as regressors within the general linear model (GLM).^13,17 With this approach one can add the SDC regressors $X_{\mathrm{SDC}}$ in addition to the task-related regressors $X_{\text{task}}$ directly to the GLM analysis

In this model, each SDC time course in $X_{\mathrm{SDC}}$ results in a beta value in ${\mathbf{\beta }}_{\mathrm{SDC}}$ estimating its contribution to the time signals in $Y$ and the part that can not be explained by the model is included in the residual term $ϵ$. By subtracting the product of design matrix $X_{\mathrm{SDC}}$ and beta values ${\mathbf{\beta }}_{\mathrm{SDC}}$ from the data $Y$, $ϵ$ results in the cleaned data set $Y_{\mathrm{GLM}\mathrm{ALL}}$

As recommended by Santosa and colleagues,¹⁷ for the GLM ALL method all eight SDCs of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ were included in the present study as regressors for running the GLM. For solving the GLM, an autoregressive iterative least square (AR-IRLS³⁶) algorithm implemented in the Brain AnalyzIR toolbox³³ ( https://github.com/huppertt/nirs-toolbox) was applied.

General Linear Model with one SDC on Both Hemispheres as Regressors (GLM BH). Similar to the GLM ALL method, the GLM BH correction method uses SDCs as regressors within the GLM in order to correct for SA. However, here only two SDCs were used for correction to pretend to be in a scenario where only a small number of SDCs are available.^41,42 As for GLM ALL, by subtracting the product of design matrix $X_{{\mathrm{SDC}}_{\mathrm{B}\mathrm{H}}}$ and beta values ${\mathbf{\beta }}_{{\mathrm{SDC}}_{\mathrm{B}\mathrm{H}}}$ from the uncorrected data $Y$, $ϵ$ results in the cleaned data set $Y_{\mathrm{GLM}\mathrm{BH}}$

For the GLM BH method, SDCs 4 and 7 were chosen because they cover both hemispheres (BH) and were included in all of the three ROIs [cf. Fig. 1(d)]. Also for this GLM method, regressors of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ of the respective channels were included in the model.

2.4.1.

Channel selection

The majority of quality measures were calculated for selected channels. To enhance comparability across SDC methods, the same channel set was used for all correction methods. Channels were not selected from the whole channel pool but selection was restricted to the task-specific ROI. The task-specific ROI was defined a priori as the ROI that most likely shows the strongest activation for a given task. For ME LEFT the task-specific ROI was M1 RIGHT, for ME RIGHT and SIM it was M1 LEFT for MI it was SMA.

For each subject, the best channel was selected separately for each task (SIM, ME LEFT, ME RIGHT, and MI) and signal type ($\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$), resulting in four tasks × two signal types = eight channels per subject. The channel selection was based on the individual output of two AR-ILS³⁶ GLMs (one for SIM data and one for real data) run on the GLM ALL corrected data because it was expected that this correction method results in the cleanest data.^13,17 For tasks ME LEFT, ME RIGHT, and MI, besides the three task-related regressors, regressors of the instruction and the breaks were added to the GLM. Additionally, first and second derivatives of each regressor were included in the model. For SIM data, the GLM included one task-related regressor consisting of all simulated onsets. From the GLM outputs, the individual channels with the highest ($\mathrm{\Delta }[\mathrm{HbO}]$) and lowest ($\mathrm{\Delta }[\mathrm{HbR}]$) beta values were selected from within the task-specific ROIs. The resulting channel frequencies across tasks and signal types are visualized in Fig. 1(h).

2.4.2.

Quality measures for signal improvement

Signal improvement was assessed with three different measures expected to be reflective of single trial noise reduction after applying SAC. It was expected that for SAC methods that effectively remove SA, quality measures for signal improvement should be enhanced relative to uncorrected data. It was also expected that the measures demonstrate a superiority of SDC methods compared with methods without SDCs.

For the individually selected channels (see Sec. 2.4.1), measures were calculated for data epoched $[-5,25]\mathrm{s}$ around stimulus onset and baseline corrected with the average of the signal from five seconds prior to stimulus onset. For ME LEFT and ME RIGHT, quality measures were first calculated separately and then averaged for each subject to generate a single value termed ME.

Scaled root mean square error (sRMSE). For this quality measure, a linear regression was fitted (using the function lm() of the ‘stats’ package v4.0.2⁴³ and rescale() of the ‘scales’ package v1.1.1 [ https://github.com/r-lib/scales] in R⁴³ in combination with RStudio⁴⁴) between a single trial epoch of either $\mathrm{\Delta }[\mathrm{HbO}]$ or $\mathrm{\Delta }[\mathrm{HbR}]$ scaled to [0, 1] ${\mathbf{y}}_s$ and a single trial canonical HRF scaled to [0, 1] ${\mathbf{x}}_s$

Then, the residuals ${\mathbf{ϵ}}_s$ of the model fit were extracted and the sRMSE was calculated as follows:

For each subject and selected channel, the single trial sRMSE values were then averaged across trials to have a single value per subject. When comparing different models resulting from different SAC methods, a low sRMSE value indicates low error of the corresponding model and is therefore interpreted as stronger signal improvement as compared with a model with a higher sRMSE value. Note that the RMSE is an important and often applied measure to validate the performance of a linear regression model. However, since the RMSE is a scale-dependent measure, it is only meaningful when comparing models based on the same data.^45–47 If it is of interest to compare several models based on different data sets, as in the present study, a scaling or normalization step should be performed.⁴⁷

Correlation (COR). For this quality measure, the Spearman rank correlation coefficient $\rho$ between a single trial epoch of either $\mathrm{\Delta }[\mathrm{HbO}]$ or $\mathrm{\Delta }[\mathrm{HbR}]$ $y$ and the corresponding canonical HRF $x$ was calculated (using the function cor() of the ‘stats’ package v4.0.2⁴³ in R⁴³ in combination with RStudio⁴⁴). To use this value for statistical analysis, each single trial $\rho$ was then Fisher’s z-transformed (function FisherZ() of the ‘DescTools’ package v0.99.38 [ https://github.com/AndriSignorell/DescTools] in R⁴³ in combination with RStudio⁴⁴) and for each subject and selected channel averaged across trials. With respect to the COR quality measure, when comparing two methods, the method with the higher COR value indicates a stronger relationship between HRF and trial and therefore a stronger signal improvement.

Contrast-to-noise ratio (CNR). The third single trial quality measure for signal improvement is the contrast-to-noise ratio^25,48,49

2.4.3.

Quality measures for spatial specificity

Spatial specificity characterizes how well defined a signal is in space. To prepare the data, the time series signal of all regular channels were epoched $[-\mathrm{5,25}]\mathrm{s}$ around stimulus onset and baseline corrected with the average of the signal from 5 s prior stimulus onset. In addition, the beta values resulting from the GLMs of each correction method were extracted for both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$.

Correlation matrices with SDCs (SDC CORMAT). For each single trial, a correlation matrix¹⁶ [Spearman correlation using the MATLAB function corr()] was calculated based on the data of each of the 16 regular channels concatenated to the data of the eight SDCs from the NO SAC data. Matrices consisted thus of correlations between all regular channels as well as between SDCs and all regular channels. For pruned channels, correlations were replaced by a NaNs. The single trial matrices were then averaged [using function nanmean()] first within a subject and then across subjects separately for each method (NO SAC, CAR, GCR, SSR, GLM ALL, and GLM BH), data type ($\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$) and task (SIM, ME LEFT, ME RIGHT, and MI), resulting in a total of $6\times 2\times 4=48$ correlation matrices with dimensions $24\times 24$. A more spatially specific signal would, on a descriptive level, be indicated by (i) reduced correlations between regular channels and SDCs as well as (ii) exclusively high correlations in the respective task-specific ROIs.

Topographical beta maps (BETA MAPS). For this quality measure no single trial data were considered, instead, the whole time series was of interest. For both data sets (SIM data and real data) and each correction method a separate AR-ILS³⁶ GLM was applied, resulting in two data sets × six correction methods = 12 GLMs per subject. For the SIM data one task-related regressor consisting of all simulated onsets was added to the model. For the real (ME/MI) data, besides the three task-related regressors (ME LEFT, ME RIGHT, and MI) also regressors of the instruction and the breaks were added. Additionally, first and second derivatives of each regressor were included in the model.

Beta maps²⁸ were generated based on the GLM outputs by averaging beta values [using the MATLAB function nanmean()] for each regular channel across subjects separately for SDC method (NO SAC, CAR, GCR, SSR, GLM ALL, and GLM BH), data type ($\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$), and task (SIM, ME LEFT, ME RIGHT, and MI), resulting in a total of $6\times 2\times 4=48$ beta maps. Beta values of pruned regular channels were replaced by an NaNs. For the topographical beta maps, spatial specificity would be judged as poor if most of the channels show strong activation and as good if mainly channels of the expected task-specific ROI show strong activation. That is, for the SIM data, only channels 6 and 8 should be activated, for ME LEFT channels of ROI M1 RIGHT, for ME RIGHT channels of ROI M1 LEFT, and for MI channels in ROI SMA should be activated.

2.5.1.

Signal improvement

For the comparison of quality measures sRMSE, COR, and CNR for each task (SIM, ME, and MI) and signal type ($\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$) a one factorial Bayesian repeated-measures ANOVA (rmBANOVA) with the factor METHOD (NO SAC, CAR, GCR, SSR, GLM ALL, and GLM BH) compared with the null model was conducted and Bayes factors ${\mathrm{BF}}_{10}$ were reported. Prior odds $P(M)$ were set as the default of $P(M)=0.5$. Posthoc analyses were only conducted if the rmBANOVA showed at least moderate (${\mathrm{BF}}_{10}>3$) evidence favoring the alternative hypothesis ($H_1$). As implemented by default in JASP, the posthoc analyses are based on the default $t$-test with a Cauchy (0, $\frac{1}{\sqrt{2}}$) prior, i.e., prior odds were set as the default value of $P(M)=0.26$ and posterior odds $P(M|\text{data})$ were corrected for multiple testing by fixing to 0.5 the prior probability that the $H_0$ holds across all comparisons.⁵⁰

For all three measures the general expectation was that (i) all correction methods show lower (sRMSE) or higher (COR and CNR) values, respectively, as compared with NO SAC and that (ii) values resulting from correction methods with SDCs (SSR, GLM ALL, and GLM BH) are lower (sRMSE) or higher (COR and CNR), respectively, as compared with values resulting from correction methods without SDCs (CAR and GCR).

2.5.2.

Spatial specificity

Correlation Matrices with SDCs (SDC CORMAT). On a purely descriptive level, it was expected that data from corrected channels show less correlation with the SDCs and that matrices show higher spatial specificity by showing stronger correlations within ROIs as compared with the NO SAC data.

For statistical analysis, a Bayesian paired $t$-test was performed on the lower triangular part of the correlation matrices for any two pairs of correction methods (total of 15 pairs) within task (SIM, ME LEFT, ME RIGHT, and MI) and signal type ($\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$). In total, $15\times 4$ tasks × 2 signal types = 120 $t$-tests were conducted for which Bayes factors ${\mathrm{BF}}_{10}$ were reported. The default of a medium prior with $P(M)=\frac{\sqrt{2}}{2}$ was applied. Across tasks and signal types, we expected (i) that matrices resulting from any correction method will differ from the correlation matrix of NO SAC as well as (ii) differences between correlation matrices resulting from methods with SDC (SSR, GLM ALL, and GLM BH) and those resulting from methods without SDCs (CAR and GCR).

Topographical beta maps (BETA MAPS). For visual inspection and qualitative comparisons, the beta values of the GLMs were used to generate average topographical beta maps for each method (NO SAC, CAR, GCR, SSR, GLM ALL, and GLM BH), signal type ($\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$), and task (SIM, ME LEFT, ME RIGHT, and MI), resulting in a total of six methods × two signal types × four tasks = 48 topographical beta maps. For statistical comparisons, for each of the 16 regular channels of any of these beta maps an one-sample Bayesian $t$-test was conducted by using the default of a medium prior with $P(M)=\frac{\sqrt{2}}{2}$ and Bayes factors ${\mathrm{BF}}_{10}$ are reported.

On a purely descriptive level, the general expectation within each task and signal type was that (i) spatial specificity increases after applying any SAC methods as compared with NO SAC and that (ii) methods with SDC correction show a higher spatial specificity as compared with methods without SDCs. Statistically, this would be reflected by the strongest Bayes factors in channels 6 and 8 of the SIM DATA and for ME LEFT and ME RIGHT by the strongest Bayes factors resulting from channels on the hemisphere contralateral to the performing hand. For MI data, the expectation is less strong, but the general expectation would be to get the strongest Bayes factors for ROI SMA.

3.1.1.

Scaled root mean square error (sRMSE)

The mean sRMSE values ($\pm SEM$) are visualized in Fig. 2. Descriptively, the sRMSE values of $\mathrm{\Delta }[\mathrm{HbO}]$ show only for the SIM data a clear difference between data of all SAC methods and the data of NO SAC [cf. Fig. 2(a)]. For $\mathrm{\Delta }[\mathrm{HbR}]$, the pattern for the correction methods of the SIM data is less clear. However, mean sRMSE values of all correction methods seem to decrease as compared with NO SAC [cf. Fig. 2(b)]. For both data types of ME and MI data, descriptively, the mean sRMSE values either increase for the correction methods or seem to be comparable to those of the NO SAC data.

Fig. 2

Illustration of the sRMSE results including all three data types (columns; SIM, ME, and MI) of (a) $\mathrm{\Delta }[\mathrm{HbO}]$ and (b) $\mathrm{\Delta }[\mathrm{HbR}]$. Boxes represent mean sRMSE values and error bars represent standard error of the mean (SEM) across participants. Below each subfigure statistical results are displayed by color coded stars representing Bayes factors ${\mathrm{BF}}_{10}>3$, which indicate evidence favoring the alternative hypothesis $H_1$, i.e., that there are differences between two correction methods. For details of statistical results including ${\mathrm{BF}}_{10}<3$, see Tables S1–S6 in the Supplementary Material.

SIM data. Detailed statistical results are reported in Tables S1 and S2 in the Supplementary Material. With respect to the $\mathrm{\Delta }[\mathrm{HbO}]$ data, the rmBANOVA revealed extreme evidence favoring the $H_1$. Post hoc tests revealed very strong to extreme evidence that all of the correction methods differed from NO SAC in that their sRMSE values were lower than that of NO SAC (all < NO SAC). Additionally, the comparisons between GLM ALL and SSR, as well as GLM BH, indicated strong evidence and between GLM ALL and GCR moderate evidence for a difference in favor of GLM ALL (GLM ALL < SSR, GLM BH, GCR).

For $\mathrm{\Delta }[\mathrm{HbR}]$, data did not provide evidence for a difference in sRMSE between SAC methods. Instead, rmBANOVA indicated anecdotal evidence favoring the $H_0$.

ME Data. Detailed statistical results are reported in Tables S3 and S4 in the Supplementary Material. For the $\mathrm{\Delta }[\mathrm{HbO}]$ data, the rmBANOVA revealed extreme evidence favoring the $H_1$. Posthoc tests indicated that evidence was strong to very strong evidence for a raised sRMSE for GCR. This included the comparisons between GCR and NO SAC, GCR and GLM ALL, and GCR and GLM BH (GCR > NO SAC, GLM ALL, GLM BH). Similarly, there was strong evidence for a higher sRMSE for CAR compared with NO SAC, CAR compared with GLM BH, and extreme evidence for CAR compared with GLM ALL (CAR > NO SAC, GLM ALL, GLM BH). Note that for the comparisons with NO SAC the opposite direction was expected, that is, expected were lower sRMSE values for GCR and CAR as compared with NO SAC.

For $\mathrm{\Delta }[\mathrm{HbR}]$, data did not provide evidence for a difference in sRMSE between SAC methods. Instead, rmBANOVA indicated anecdotal evidence favoring the $H_0$.

MI Data. Detailed statistical results are reported in Tables S5 and S6 in the Supplementary Material. For the $\mathrm{\Delta }[\mathrm{HbO}]$ data, the rmBANOVA revealed anecdotal evidence favoring the $H_0$. Regarding the $\mathrm{\Delta }[\mathrm{HbR}]$ data, the rmBANOVA revealed moderate evidence favoring the $H_0$.

3.1.2.

Correlation (COR)

The means of the Fisher’s $z$-transformed Spearman correlations COR ($\pm \mathrm{SEM}$) are illustrated in Fig. 3. Descriptively, for both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ SIM data mean COR values were higher for all SAC methods as compared with NO SAC. For ME $\mathrm{\Delta }[\mathrm{HbO}]$ data, only the SAC methods with SDCs showed an increase in COR as compared with NO SAC, whereas both methods without SDCs showed a decrease in COR. For $\mathrm{\Delta }[\mathrm{HbR}]$ of ME data and both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ of the MI data, COR values of the SAC methods either decreased or were comparable to NO SAC.

Fig. 3

Illustration of the Fisher’s $z$-transformed Spearman correlation (COR) results including all three data types (columns: SIM, ME, and MI) of (a) $\mathrm{\Delta }[\mathrm{HbO}]$ and (b) $\mathrm{\Delta }[\mathrm{HbR}]$. Circles represent mean COR values and error bars represent standard error of the mean across participants. Below each subfigure statistical results are displayed by color coded stars representing Bayes factors ${\mathrm{BF}}_{10}>3$ which indicate evidence favoring the alternative hypothesis $H_1$, i.e., that there are differences between two correction methods. For details of statistical results including ${\mathrm{BF}}_{10}<3$, see Tables S7–S12 in the Supplementary Material.

SIM data. Detailed statistical results are reported in Tables S7 and S8 in the Supplementary Material. For $\mathrm{\Delta }[\mathrm{HbO}]$, the rmBANOVA revealed extreme evidence favoring the $H_1$. Posthoc tests indicated extreme evidence for the differences between NO SAC and all other correction methods, resulting from higher COR values of the correction methods as compared with NO SAC. For GLM ALL, there was also extreme evidence for higher COR values compared with GCR, SSR, and GLM BH and very strong evidence for higher COR values compared with CAR (GLM ALL > GCR, SSR, GLM BH, and CAR).

For $\mathrm{\Delta }[\mathrm{HbR}]$, the rmBANOVA revealed moderate evidence favoring the $H_1$. Posthoc tests indicated higher COR values for GLM ALL and GLM BH compared with NO SAC, with extreme evidence for the difference between NO SAC and GLM ALL and strong evidence for differences between GLM BH and NO SAC (GLM ALL and GLM BH > NO SAC).

ME data. Detailed statistical results are reported in Tables S9 and S10 in the Supplementary Material. For the $\mathrm{\Delta }[\mathrm{HbO}]$ data, the rmBANOVA revealed extreme evidence favoring the $H_1$. Posthoc tests indicated extreme evidence for differences between GLM ALL and and GCR and between GLM ALL and CAR with higher COR values for GLM ALL in both comparisons (GLM ALL > GCR, CAR). Moreover, there was evidence for particularly low COR values for GCR, with moderate evidence for differences between GCR and SSR and between GCR and GLM BH (SSR, GLM BH > GCR). There was also moderate evidence for the difference between CAR and GLM BH, with higher COR values for GLM BH (GLM BH > CAR).

With respect to the $\mathrm{\Delta }[\mathrm{HbR}]$ data, the rmBANOVA revealed anecdotal evidence favoring the $H_1$.

MI data. Detailed statistical results are reported in Tables S11 and S12 of the Supplementary Material. For both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$, the rmBANOVAs revealed anecdotal evidence favoring the $H_1$.

3.1.3.

Contrast-to-noise ratio (CNR)

The mean CNRs ($\pm \mathrm{SEM}$) are illustrated in Figs. 4(a) and 4(b). Descriptively, with respect to the SIM data of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ all SAC methods show increased CNR values as compared with NO SAC. For $\mathrm{\Delta }[\mathrm{HbO}]$ of the ME data, all SAC methods with SDCs show increased CNR values as compared with NO SAC and for $\mathrm{\Delta }[\mathrm{HbR}]$ of the ME data, only both GLM methods show increased CNR values. In contrast, methods without SDCs show slightly decreased CNR values as compared with NO SAC. With respect to both data types of the MI data, all methods seem to be comparable with NO SAC.

Fig. 4

Illustration of the contrast-to-noise ratio (CNR) results including all three data types (columns: SIM, ME, and MI) of (a) $\mathrm{\Delta }[\mathrm{HbO}]$ and (b) $\mathrm{\Delta }[\mathrm{HbR}]$. Diamonds represent mean CNR values and error bars represent standard error of the mean across participants. Below each subfigure statistical results are displayed by color coded stars representing Bayes factors ${\mathrm{BF}}_{10}>3$, which indicate evidence favoring the alternative hypothesis $H_1$, i.e., that there are differences between two correction methods. For details of statistical results including ${\mathrm{BF}}_{10}<3$, see Tables S13–S18 in the Supplementary Material.

SIM data. Detailed statistical results are reported in Tables S13 and S14 in the Supplementary Material. With respect to the $\mathrm{\Delta }[\mathrm{HbO}]$ data, the rmBANOVA revealed extreme evidence favoring the $H_1$. As compared with NO SAC, posthoc tests revealed extreme evidence for differences to GLM ALL, very strong evidence for differences to GLM BH, strong evidence for differences to CAR and SSR and moderate evidence for differences to GCR, all resulting from higher CNR values in favor of the correction methods as compared with NO SAC (NO SAC < CAR, GCR, SSR, GLM BH, and GLM ALL). Moreover, there was very strong evidence for differences between GLM ALL and GCR and between GLM ALL and GLM BH and moderate evidence for differences between GLM ALL and CAR as well as between GLM ALL and SSR, reflecting the overall highest CNR for GLM ALL (GLM ALL > GCR, SSR, GLM BH, and CAR).

Regarding the $\mathrm{\Delta }[\mathrm{HbR}]$ data, the rmBANOVA revealed anecdotal evidence favoring the $H_0$.

ME data. Detailed statistical results are reported in Tables S15 and S16 in the Supplementary Material. For $\mathrm{\Delta }[\mathrm{HbO}]$, the rmBANOVA revealed extreme evidence favoring the $H_1$. Posthoc tests revealed extreme evidence for differences between GLM ALL and NO SAC, CAR as well as GCR, strong evidence for differences between GLM ALL and GLM BH, and moderate evidence for differences between GLM ALL and SSR, all resulting from higher CNR values of GLM ALL as compared with all other methods (GLM ALL > CAR, GCR, SSR, GLM BH, and NO SAC). Evidence was strong for a difference between GCR and SSR and extreme evidence was found for differences between GCR and GLM BH, with reduced CNR values for GCR (SSR, GLM BH > GCR). Moreover, evidence was strong for a smaller CNR for CAR compared with GLM BH and moderate for a smaller CNR for CAR compared with SSR (GLM BH, SSR > CAR). There was extreme evidence for differences between GLM BH and GCR, and strong evidence was found for the difference between NO SAC and GLM BH as well as CAR and GLM BH with higher CNR for GLM BH (GLM BH > NO SAC).

For $\mathrm{\Delta }[\mathrm{HbR}]$, the rmBANOVA indicated extreme evidence favoring the $H_1$. Posthoc tests indicated that this was related to reduced CNR values for CAR compared with GLM ALL, for which strong evidence was found (GLM ALL > CAR). Similarly, there was strong evidence for differences between GCR and GLM ALL and moderate evidence between GCR and GLM BH reflecting smaller CNR values for GCR (GLM ALL, GLM BH > GCR).

MI data. Detailed statistical results are reported in Tables S17 and S18 in the Supplementary Material. Regarding both the $\mathrm{\Delta }[\mathrm{HbO}]$ and the $\mathrm{\Delta }[\mathrm{HbR}]$ data, the rmBANOVA revealed anecdotal evidence favoring the $H_0$.

3.2.1.

Correlation matrices with SDCs (SDC CORMAT)

Overall and across tasks, descriptively, the correlation matrices of $\mathrm{\Delta }[\mathrm{HbO}]$ for NO SAC showed a strong influence of the SA across channels, indicated by high correlation coefficients between regular channels as well as between regular channels and SDCs. The high correlations are already attenuated but not absent when considering $\mathrm{\Delta }[\mathrm{HbR}]$ [cf. Figs. 5(a), 6(a), and 7(a)]. For both signal types correlations with SDCs are highly reduced after applying correction methods irrespective of using SDCs for correction or not. However, negative correlations between regular channels as well as between regular channels and SDCs resulted from correction with methods CAR and GCR.

Statistical analyses were run with Bayesian paired $t$-tests, performed on the lower triangular part of the correlation matrices for any two pairs of correction methods. The results of all $t$-tests are summarized in Tables S19 and S20 in the Supplementary Material.

SIM data. As visualized in Fig. 5(a), for both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$, nearly all SDC CORMAT comparisons indicated extreme evidence for differences between correction methods. Exceptions were comparisons between CAR and GCR with respect to $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ as well as between methods GLM ALL and GLM BH for the $\mathrm{\Delta }[\mathrm{HbR}]$ data.

Fig. 5

Illustration of the spatial specificity results for SIM data. (a) The averaged correlation matrices of $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbO}]$. Colored stars below a correlation matrix represent Bayes factors of the statistical comparison between the respective matrix and the matrix belonging to the correction method that is indicated by the colored polygon above the Bayes factor. Note that the order of the regular channels in the matrix is arbitrarily chosen in that the order always starts in the lower left corner with channels outside the ROI (in ascending order, indicated by white circles) and are followed by channels inside the ROI (in ascending order, indicated by gray circles and corresponding to ROI M1 left). (b) Visualizes the topographical beta maps (BETA MAPS) of $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$. Each circle represents a channel with its corresponding mean beta value. The star within a circle represents the Bayes factor resulting from the Bayesian $t$-test of the beta values of the respective channel. The black oval marks the ROI, M1 LEFT.

Fig. 6

Illustration of the spatial specificity results of ME LEFT and ME RIGHT data. (a) The averaged correlation matrices of $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$. Colored stars below a correlation matrix represent Bayes factors of the statistical comparison between the respective matrix and the matrix belonging to the correction method that is indicated by the colored polygon above the Bayes factor. Note that the order of the regular channels in the matrix is arbitrarily chosen in that the order always starts in the lower left corner with channels outside the ROI (in ascending order, indicated by white circles) and are followed by channels inside the ROI (in ascending order, indicated by gray circles and corresponding to ROI M1 left). (b) shows the topographical beta maps (BETA MAPS) of $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$. Each circle represents a channel with its corresponding mean beta value. The star within a circle represents the Bayes factor resulting from the Bayesian $t$-test of the beta values of the respective channel. The black ovals mark the ROIs, M1 LEFT and RIGHT.

Fig. 7

Illustration of the spatial specificity results for SIM data. (a) The averaged correlation matrices of $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$. Colored stars below a correlation matrix represent Bayes factors of the statistical comparison between the respective matrix and the matrix belonging to the correction method that is indicated by the colored polygon above the Bayes factor. Note that the order of the regular channels in the matrix is arbitrarily chosen in that the order always starts in the lower left corner with channels outside the ROI (in ascending order, indicated by white circles) and are followed by channels inside the ROI (in ascending order, indicated by gray circles and corresponding to ROI M1 left). (b) The topographical beta maps (BETA MAPS) of $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$. Each circle represents a channel with its corresponding mean beta value. The star within a circle represents the Bayes factor resulting from the Bayesian $t$-test of the beta values of the respective channel. The black oval marks the ROI, SMA.

ME data. For both ME LEFT and ME RIGHT [cf. Fig. 6(a)], results indicated extreme evidence for differences between correlation matrices. Exceptions were the comparisons between CAR and GCR for both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ of ME LEFT and ME RIGHT and between SSR and GLM BH for both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ of ME LEFT and for $\mathrm{\Delta }[\mathrm{HbR}]$ of ME RIGHT.

MI data. As illustrated in Fig. 7(a), all comparisons between correlation matrices resulted in extreme evidence for differences. Sole exception was the comparison between CAR and GCR. This pattern was identical for $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$.

3.2.2.

Topographical beta maps (BETA MAPS)

The beta maps for each method are visualized for SIM in Fig. 5(b), ME in Fig. 6(b), and MI in Fig. 7(b). As expected, descriptively, the beta maps of the SIM data showed the strongest activation in channel 6 (100% amplitude) directly followed by channel 8 (50% amplitude) across methods and data types. For the real data, visual inspection shows all SAC methods differ clearly from the NO SAC data for which beta values indicated comparable activation across channels. For instance, a clearer lateralization in the ME beta maps can be seen after correction [cf. Fig. 6(b)]. In addition, clear differences between SAC methods with and without SDCs are visible across tasks. That is, for SAC methods without SDCs but not for SAC methods with SDCs beta maps show either unexpectedly negative ($\mathrm{\Delta }[\mathrm{HbO}]$) or positive ($\mathrm{\Delta }[\mathrm{HbR}]$) beta values for some channels.

SIM data. With respect to the $\mathrm{\Delta }[\mathrm{HbO}]$ BETA MAPS [cf. Fig. 5(b) top], for all methods the one sample Bayesian $t$-tests of channels 6 and 8 revealed extreme evidence for differing positively from zero. This was expected because these channels contain the simulated hemodynamic response. Two additional channels differed positively from zero for method GLM ALL. For methods, CAR and GCR one sample Bayesian $t$-tests indicated at least moderate evidence for a difference from zero at five and seven additional channels, respectively. Differences were mostly but not exclusively due to negative beta values (cf. Table 21 in the Supplementary Material for exact values).

Regarding the $\mathrm{\Delta }[\mathrm{HbR}]$ BETA MAPS [cf. Fig. 5(b) bottom], results indicated as expected extreme evidence for a negative difference from zero for data of simulated channels 6 and 8 of all methods. Only methods CAR, GCR, and SSR resulted in additional channels mostly with evidence for positive differences from zero (cf. Table 22 in the Supplementary Material for exact values).

ME data. Regarding the $\mathrm{\Delta }[\mathrm{HbO}]$ data of both ME LEFT and ME RIGHT [cf. Fig. 6(b) top two rows], BETA MAPS for method NO SAC showed moderate to extreme evidence for differing positively from zero across all channels (16/16). For CAR one sample Bayesian t-tests indicated that a small number of channels differed positively (ME LEFT: 1/6 in ROI M1 RIGHT; ME RIGHT: 2/6 in ROI M1 LEFT) or for negatively (both ME LEFT and ME RIGHT: 2/16 in ROI SMA) from zero. Results were similar for GCR (positive differences ME LEFT: 3/16 in all ROIs, 1/16 in M1 RIGHT; ME RIGHT: 1/16 in ROI SMA; negative differences ME LEFT: 2/16 in all ROIs; ME RIGHT: 1/16 in ROI SMA). Statistical analysis of BETA MAPS of method SSR showed for most channels moderate to extreme evidence for differing positively from zero (ME LEFT: 12/16 in all ROIs, 5/6 in M1 RIGHT; ME RIGHT: 12/16 in all ROIs, 5/6 in M1 LEFT). BETA MAPS for methods GLM ALL and GLM BH differed with moderate to extreme evidence positively from zero in a medium number of channels (GLM ALL: ME LEFT: 10/16 in all ROIs, 5/6 in M1 RIGHT; ME RIGHT: 7/16 in all ROIs, 5/6 in M1 LEFT; GLM BH: ME LEFT: 8/16 in all ROIs, 4/6 in M1 RIGHT; ME RIGHT: 7/16 in all ROIs, 5/6 in M1 LEFT) (cf. Table S23 in the Supplementary Material for exact values).

For $\mathrm{\Delta }[\mathrm{HbR}]$ of NO SAC, for ME LEFT data 14 of 16, and for ME RIGHT data of 12 of 16 channels provided evidence for moderate to extreme evidence for negatively differing from zero as indicated by the results of one sample Bayesian $t$-tests [cf. Fig. 6(b) bottom two rows]. Statistical results for methods CAR and GCR revealed in more channels evidence for differing positively from zero (CAR: ME LEFT: 3/16 in all ROIs; ME RIGHT: 5/16 in all ROIs; GCR: ME LEFT: 3/7 in all ROIs; ME RIGHT: 4/16 in all ROIs) than for differing negatively from zero (CAR: ME LEFT: 1/6 in ROI M1 RIGHT; ME RIGHT: 0/16; GCR: ME LEFT: 3/16 in all ROIs; ME RIGHT: 1/8 in ROI SMA), which would however be the expected direction. Results for SSR showed moderate to extreme evidence for most channels to differ negatively from zero (ME LEFT: 14/16 in all ROIs, 6/6 in M1 RIGHT; ME RIGHT: 10/16 in all ROIs, 4/6 in M1 LEFT). For GLM ALL and GLM BH, the number of channels with evidence for differing negatively from zero was similar or slightly reduced compared with SSR (GLM ALL: ME LEFT: 13/16 in all ROIs, 5/6 in M1 RIGHT; ME RIGHT: 8/16 in all ROIs, 4/6 in M1 LEFT; GLM BH: ME LEFT: 13/16 in all ROIs, 5/6 in M1 RIGHT; ME RIGHT: 11/16 in all ROIs, 4/6 in M1 LEFT). For more details see Table S24 in the Supplementary Material.

MI data. Regarding the $\mathrm{\Delta }[\mathrm{HbO}]$ data of NO SAC, data of 15 out of 16 channels across all ROIs provided moderate to extreme evidence for differing positively from zero [cf. Fig. 7(b) top]. Regarding CAR and GCR, there was evidence for one channel to differ positively (ROI M1 LEFT) and for one channel to differ negatively (SMA) from zero. For SSR, 12 out of 16 channels across all ROIs differed positively from zero with moderate to extreme evidence. For GLM ALL, data of 8 out of 16 and for GLM BH 10 out of 16 channels across all ROIs provided moderate to very strong evidence for differing positively from zero (cf. Table S25 in the Supplementary Material for exact values).

With respect to $\mathrm{\Delta }[\mathrm{HbR}]$, one sample Bayesian $t$-tests provided moderate to extreme evidence for beta values to differ negatively from zero in 11 out of 16 channels across all ROIs for NO SAC [cf. Fig. 7(b) bottom]. For method CAR, for 5 out of 16 channels across all ROIs evidence was at least moderate for a positive difference from zero, whereas for three further channels evidence was at least moderate for a negative difference (ROIs M1 LEFT and M1 RIGHT). For GCR, there was evidence for differing positively from zero for 5 out of 16 channels across all ROIs and for differing negatively from zero for 4 out of 16 channels across all ROIs. Method SSR was linked to evidence indicating a negative difference from zero for 10 out of 16 channels across all ROIs. For GLM ALL, this was only the case for 7 and for GLM BH for 8 out of 16 channels and across all ROIs (cf. Table S26 in the Supplementary Material for exact values).

4.3.1.

Common average reference (CAR)

The CAR method²² is a spatial filter for which in the present study the averaged signal across all regular channels was subtracted from each individual channel. Of the two methods without SDCs, CAR resulted in a greater signal improvement. Surprisingly, for the SIM data, signal improvement following CAR was larger as compared with SSR and GLM BH across all three quality measures. However, for the real data CAR seems to overcorrect the data. Overcorrection was also suggested by spatial specificity measures with negative correlations (SDC CORMAT) and negative beta values (BETA MAPS) across all tasks. On the other hand, channels indicated as active by BETA MAPS were comparable to the methods with SDCs, as for instance reflected in the lateralization pattern for the ME DATA [cf. Fig. 6(b)]. This indicates a good spatial specificity of CAR for strongly activating tasks.

The negativity in the signal following CAR is not surprising given that strongly activated channels affect the amplitude of the averaged signal and channels with no or lower amplitudes get pushed into the reversed direction. A workaround would be to omit the channels in the critical ROI from the averaged signal. However, this only works if the ROI covers the peak activation and when channels outside of the ROI are not expected to be active. This is rather unlikely for any but purely perceptual tasks, and even here it is doubtful. A more advanced version of this method could be to regress out the average across all channels, e.g., by using the SSR¹⁹ method,¹⁵ instead of simply subtracting it. Moreover, for the implemented version of the CAR method or any modification thereof, a proper channel quality check and channel pruning are important in order to avoid to subtract data of poor quality into a given channel.

Overall, although it can result in artificial negativity of inactive channels, CAR seems to be a good and easy-to-implement choice if no SDCs are available and if only data of one channel or a very circumscribed ROI covering the peak activation will be used for (online) processing.

4.3.2.

Global component removal (GCR)

The GCR method²³ is a spatial filter that combines a Gaussian kernel smoothing by taking into account the MNE coordinates of the channel positions and a singular value decomposition. Regarding the signal improvement measures, for SIM the results indicated a signal improvement across all quality measures. For the real data the opposite was evident, that is, results indicated a signal decline. Overcorrection was indicated by strong negative correlations between regular channels and between regular channels and SDCs, and by negative beta activation across tasks and data types. The reason for the artificial negative activation is probably the combination of the relatively low spatial coverage and the large-sized kernel.²³ Although the probe in the present study falls into the recommended minimum area²³ it seems that this method needs a much higher head coverage. In a recent study,²⁴ it was shown that higher-density data corrected by means of the GCR was comparable with that corrected by a regression-based approach using SDCs. However, this approach considered only the first two principal components²⁴ resulting from a principal component analysis (PCA) of the SDC data and therefore missed potentially important spatial information that could have made a difference to the results.

On the whole, the GCR method performed well for SIM data but not for the real data. Therefore, it should be used with care and should only be applied with sufficiently high spatial coverage to avoid overcorrection. To decide whether the number of channels is sufficient, a comparison of the beta output of a GLM (as done for BETA MAPS) before and after correction might be helpful.

4.3.3.

Short separation regression (SSR)

The SSR method¹⁹ is a simple regression method. In the present study, we applied SSR-based correction using the spatially closest SDC to a given regular channel. For both signal improvement and spatial specificity measures, out of the correction methods with SDCs the SSR performed poorest and in some cases performed worse than methods without SDCs. Moreover, for real data, the spatial specificity measures suggested a stronger impact of SA after correction compared with the two GLM methods (cf. Figs. 6Fig. 7Fig. 8–9). One explanation for this limited performance might be that the closest SDC is not necessarily the best choice for SSR. Alternatively, the SDC with the highest correlation, with the average across SDCs or with the first $n$ principal components resulting from a PCA based on all SDCs²⁴ could be used.

Fig. 8

Illustration of the averaged time series data of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ including all channels for all correction methods. For method NO SAC, additionally, the averaged time series data of the SDCs is visualized (top left panel). The black ovals mark the ROI, M1 LEFT and RIGHT.

Fig. 9

Illustration of the averaged time series data of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ including all channels for all correction methods. For method NO SAC, additionally, the averaged time series data of the SDCs is visualized (top left panel). The black oval marks the ROI, SMA.

In spite of these shortcomings, given the simplicity of the SSR method especially regarding online implementation, it can still be recommended. However, one important step before applying this method, especially when considering the closest SDC for correction, is the signal quality check. If SDCs or regular channels of poor signal quality are not properly pruned before correction, regular channels might either be corrupted, remain uncorrected or SA will get amplified.

4.3.4.

General linear model with all SDCs as regressors (GLM ALL)

In the present study, of the selected SAC methods, the GLM ALL method^13,17 was overall performing best. For this correction method, all eight SDCs of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ were included as regressors to the SAC GLM, because previous work provided evidence that the SA is not a homogeneous but rather a heterogeneous artifact^13,57 and therefore, a higher number of SDCs covering a large portion of the head improves correction.¹⁷

That the GLM ALL method seems to be the most accurate method for SA correction is in line with previous findings.^13,17 Although the AR-ILS³⁶ GLM is generally computational more intensive as compared with the ordinary least squares (OLS) GLM, its feasibility for online application has been shown⁵⁸ and it has been implemented in online preprocessing tools such as TurboSatori.⁵⁹ Moreover, AR-ILS can handle the presence of motion artifacts very well.³⁶ This is, however, not the case for the other SA correction methods and hence, could have been a gain for the GLM-based approaches. However, if motion artifacts were present in our data sets then usually at the beginning/end of the recording or during the breaks and only in very rare cases during the task blocks. Hence, even if some motion artifacts survived the TDDR correction, the advantage of the AR-ILS GLM over the other correction methods should be low.

In spite of the good performance of GLM ALL, it is likely that results can still be improved by adding more information to the GLM. For instance, gains can be expected by adding additional regressors based on auxiliary recorded physiological measures like blood pressure and oxygen saturation,⁵⁵ information extracted from the fNIRS measurement itself like Mayer waves,¹³ or a global signal regressor generated from the SDCs.¹³ Another promising way would be to increase the number of SDCs as it was shown that the number of SDCs has a strong effect on the correction performance of the GLM, i.e., the more information in terms of SDCs is available the better the correction performance.¹⁷ Another promising approach would be to cover blind spots by, e.g., adding short separation sources to the detectors the same way as for the available hardware solutions short separation detectors are added to the sources.¹³ Moreover, in the present study for method GLM ALL but also for method GLM BH, two GLMs were applied consecutively. That is, first a GLM was applied only for SA correction by including all SDCs as regressors. From this GLM, SAC corrected data was extracted for further evaluation, among others based on a second GLM in which task-induced activation was modeled. This two-step procedure was chosen because it provided the SAC corrected data required as input for all quality measures. In principle, if only the task-related activation is of interest, it would also be possible to run a single GLM with both the SDC regressors and the task-related regressors. Indeed, for BCI and NFB online scenarios in which the feedback is based on the GLM output this would seem like the natural choice. One major issue with this one-step approach is, however, the risk of multicollinearity due to the strong task-evoked activity in the SDCs, which can lead to unpredictable results.¹⁷ To prevent this problem a two-step GLM approach as applied in the present study should be considered, even if the main focus of the GLM analysis is the extraction of task-related activation.

In sum, the GLM ALL method strongly increased both signal quality and spatial specificity in the present data sets. Accordingly, of the evaluated methods, it is the best choice for SA correction if SDCs are available. Within a real-time GLM calculation, it can be easily applied for online SA correction either by implementing it independently or by using available toolboxes.⁵⁹

4.3.5.

General linear model with one SDC on both hemispheres as regressors (GLM BH)

For the GLM BH method only two out of the eight available SDCs (SDC 4 and 7 of both $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$) were included in the GLM for correction. The intention was to evaluate GLM correction performance for constellations in which only a small number of SDCs are available. GLM BH produced overall satisfying results that were, however, inferior to the results from GLM ALL. The finding of overall inferiority to a GLM relying on larger numbers of SDCs is in line with previous findings.¹⁷ On closer inspection, results for GLM BH were comparable to those of the SSR method. For spatial specificity, however, GLM BH results were, descriptively, more comparable to GLM ALL. As a workaround for unavailable SDCs, some researchers help themselves by sacrificing one or two regular channels and placing the vacated sources and detectors as close as possible to each other to create a self-made SDC.^41,42 This approach should be considered with care because depending on the manufacturer the resulting distance between source and detector is likely $\ge 10\mathrm{mm}$. SDCs placed at distances $\ge 10\mathrm{mm}$ may capture different signals than SDCs placed at the recommended 8 mm as used in the present study because of NIR light penetrating to cortical brain tissue.¹⁸ This will likely influence the recovery of the true brain signal and bears the risk of overcorrection.¹⁸ This is, however, similarly true for all correction methods using SDCs, regardless of the number of SDCs used. Moreover, if not handled correctly, e.g., by the recording software, self-made SDCs can possibly lead to a saturation of the self-made SDCs. Nevertheless, if properly used, self-made SDCs are should still be better than using no SDCs.

In sum, with a reduced set of 8-mm SDCs, GLM BH achieved satisfying results for GLM-based correction, but results were somewhat inferior to GLM ALL. Similar to GLM-ALL, this method can easily be implemented online as part of a real-time GLM calculation.

4.3.6.

Recommendations

Based on the present results and in line with previous findings,^13,17 the best performing SA correction method is the GLM ALL method directly followed by the GLM BH and the SSR method. For the GLM ALL method all eight SDCs were used as regressors, for the GLM BH method two predefined SDCs were used and for SSR only the spatially closest SDC was used for correction. These results indicate that correction performance increases with the number of SDCs used for the correction, which is in line with findings of Santosa and colleagues¹⁷ who investigated the effect of the number of SDCs for a GLM-based correction approach. In the study by Santosa et al.¹⁷ the maximum number of available SDCs was likewise eight, however, the authors predicted that the performance will likely increase with more SDCs. The methods without SDCs showed a lower performance compared with the methods with SDCs and tended to overcorrect the data. However, based on the performance results of CAR and GCR and given the strong influence of the SA artifact, the application of either of these or another correction method^17,27 should be preferred over the option to neglect the SA correction step. Based on the present results and especially if only a limited optode coverage is available (cf. next Sec. 4.4), we would recommend the CAR method over the GCR method because the overcorrection seems slightly lower and CAR results in a higher spatial specificity as compared with the GCR method.