2.1.

Functional Near-Infrared Spectroscopy Linear Statistical Models

In fNIRS studies, measurements are made between a series of discrete light emitter and detector positions on the head. The intensity of the light transmitted from an emitter to a detector position is sensitive to changes in the optical absorption of the underlying tissue along the light’s path through the tissue. This results in a time course of optical transmission measurements. Concurrent and colocalized measurements at two or more optical wavelengths are then used for estimating the changes in the concentration of ${\mathrm{HbO}}_2$ and Hb using the modified Beer–Lambert law¹ and given by

In many fNIRS studies, these changes in hemoglobin are recorded over time during some variation of a repeated cognitive task(s). The fNIRS signals between each source-to-detector pair are analyzed using a general linear regression model to test for statistical differences between the baseline and the task conditions for each scan. This approach is similar to fMRI, although several differences in the structure of noise in fNIRS compared to other modalities should be noted (see Ref. 9, for review). First-level statistical models for examining evoked signal changes are given by a regression model described by the equation

In general, statistical testing of the regression coefficients ($\beta$) to either compare these to zero (baseline) or between two task conditions is used for inferring the location and level of brain activity using a Student’s $t$-test. While there are a few approaches to solving Eq. (2), which have been used in fNIRS analysis including generalized linear models (e.g., Refs. 9, 16, 18, and ordinary least-squares regression¹⁹), in this paper, we will use the NIRS-specific generalized linear model formulation proposed by Barker et al.¹⁵

2.2.

Prewhitened and Robust Linear Model

In this work, we used a NIRS-adopted version of the general linear model, as previously described in Barker et al.¹⁵ This approach uses an autoregressive prewhitening method and an iteratively reweighted least squares (AR-IRLS) to control type-I errors in the fNIRS statistical model. This approach has been previously shown to have excellent sensitivity-specificity properties in comparison to ordinary least squares and other general linear regression models for fNIRS data containing physiological noise (serial-correlated noise structures) and motion-related artifacts (heavy-tailed outliers).¹⁵ In brief, this regression model uses an $n$’th order autoregressive filter ($W_{AR}$) determined by an Akaike model-order (AIC) selection to whiten both sides of this expression, e.g.,

As described in Barker et al.,¹⁵ the regression model is first solved using robust regression and the residual noise is then fit to an AR model. This filter ($W_{AR}$) is applied to both sides of the original model and then resolved and repeated until convergence. This AR filter removes serially correlated errors in the data that result from physiological noise and/or motion artifacts. AR whitening, however, does not address the heavy-tailed noise from motion artifacts. To do this, the AR-whitened model is solved using robust weighted regression, which is a procedure to iteratively downweight outliers, such as motion artifacts

Using this model, the regression coefficients ($\beta$) and their error covariance is estimated, which is used in defining the statistical tests between task conditions and baseline. The regression model is solved sequentially for each data file for each subject. All source–detector pairs within a file are solved concurrently yielding a full covariance model of the noise, which is used in group-level and region-of-interest analyses. The estimate of $\beta$ and its covariance matrix is given by the expressions

2.3.

Linear Regression Design Model ($X$)

In the fNIRS field, there is a controversy on what type of design model should be used for linear regression analysis. Namely, it is not clear what effect the assumption of a fixed temporal shape of the response or the use of specific time window for computing the contrast from a block average or deconvolution model has on the statistical power for brain signals from special populations (e.g., infants) that may not have the “typical” adult-like response. To examine this question, we note that in Eq. (2), the design model ($X$) describes the timing of the experimental task and conditions in the experiment and sets up the hypothesis that can be tested from the regression coefficients. This design can take the form of either an unrestricted (full deconvolution) model or an impulse response model. In the impulse response model, the response is assumed to be linearly additive such that (for example) the brain response to a 30-s task is the same as 30 repeated 1-s duration tasks building up on top of each other. In this model, only the hemodynamic impulse response function is estimated, which generally is thought to peak around 6 to 8 s and recover to baseline around 12 to 15 s. This impulse response is assumed to be the same for tasks of any duration. Thus, an advantage of this model is that it can be used even if the task blocks have differing durations, as is often the case in self-paced paradigms. In comparison, for the full deconvolution model, the entire time course of the block is estimated. For example, in a 30-s duration task, the response over a window of about 45 s (task plus 15-s recovery to baseline) would be used. Because the full deconvolution model does not assume the response to be linearly additive over the task block, it can be used for examining the saturation or habituation effects of a prolonged task duration. However, in this 30-s duration example, the full deconvolution model would require three times more regression terms compared to the impulse response version (modeling 45 s of data instead of just the 15-s impulse window). Note that this discussion of additive linear assumptions applies to blocked design or slow event-related designed experiments. In rapid event-related designs, both the full deconvolution and impulse response models assume linearity. Examples of somewhat simplified full and impulse response models are given in Eqs. (6a) and (6b).

Full deconvolution model

Impulse response model

The full model design matrix is constructed where each column is a vector of 1s at the task onset of each block. Each column shifts the onset by one entry [Eq. (6b)]. If the vectors of coefficients ($\beta$) were plotted down the columns, this would look like the time course of either the full or the impulse response. Note that the model described in Eq. (6a) reduces to the block-averaging equation in the limit that the events are nonoverlapping (specifically, there are no off-diagonal terms in the matrix $X^{T.}X$ in the Gauss–Markov equation).²¹ Under this condition, the regression model reduces to a weighted backprojection ($\beta =1/n^{.}{X}^{T.}Y$), which is the block-averaging operation in terms of the estimation of the beta weights. In block-averaging, the variance of the coefficients is estimated separately for each point, whereas in regression its assumed to be normally distributed with parameter ${\text{sigma}}^2$ [Eq. (5b)]. However, this distinction is moot when computing the contrast over a time window, as the parametric t-test, which is often used in these block-averaging results, reintroduces the normally distributed assumption. Thus, as commonly implemented by many literatures, e.g., in the HOMER NIRS software, the Student’s $t$-test contrast from block-averaging and deconvolution are equivalent (see Ref. 19, for more details).

For either the full deconvolution model or the impulse response model, a basis set (denoted as matrix $\mathrm{\Gamma }$) can be used to reparameterize the model such that

This models the time course (column of the coefficient vector) as a linear combination of the basis. This basis set matrix can be either a complete (lossless) or a restricted (e.g., smoothing or canonical) basis. In this paper, we have explored several of the commonly used basis sets (e.g., canonical HRF; see Secs. 2.3.1–2.3.5, for details).

For either the original Eq. (2) or the basis-set reparametrized Eq. (7) models, then, the estimated model coefficient vector ${\beta }_i(k)$ is used to calculate the $t$-value for a two-tailed $t$-test to test the null hypothesis $c^T{\beta }_i(k)=0$.²² In this study, the $t$-statistics of channel $i$ at time step $k$ are obtained using

4.1.

Sensitivity of Regression to Systematic Bias in the Hemodynamic Response Function Shape

Using the range of derived hemodynamic shapes (shown in Fig. 1), we have examined the sensitivity and specificity of the various regression bases to systematic bias in the shape of the underlying hemodynamic response compared to the basis used in recovery. For consistency across simulations, the bias is defined relative to the canonical HRF model as the R-squared fit between the simulated and the canonical models, which we denote as the “percent-overlap”. In Figs. 2(a)–2(d), we show the AUC from ROC analyses for the (i) canonical (blue solid line), (ii) canonical+RE(deriv) (red dashed line), (iii) FIR-IRF (yellow dashed line), (iv) FIR (dark-green and dotted line), and (v) nonlinear (light-blue and dash-dot line) deconvolution models as a function of the match between the generative response shape and the idea canonical HRF model. In Figs. 2(e)–2(f), we show the sensitivity of the ROC models (at $p=0.05$). The results are presented at $\mathrm{SNR}=0.3$, 0.5, 1.0, and 2.0 levels as shown in panels (a) and (e), (b) and (f), (c) and (g), (d) and (h), respectively. In all panels, the $x$-axis is the overlap between these models (100% overlap equals an exact match to the canonical HRF). A task duration of 5 s is used in these regression models. The rationale for showing this cross section at 5-s duration of the total results will be more apparent from the further results sections.

Fig. 2

Comparison of AUC and sensitivity at $p=0.05$ from ROC analysis at various mismatch levels (relative bias with the canonical HRF model) between 20% and 100% overlap for 5-s task duration. An overlap of 100% means a perfect match to the canonical HRF. Overlap is mathematically defined here as the $R^2$ fit. In each panel, it shows the mismatch or (overlap) for both AUC (first row) and sensitivity (second row) using five selected regression models: canonical, canonical+RE(deriv), FIR-IRF, FIR, and nonlinear. First to fourth columns show the performance at SNR levels of 0.3, 0.5, 1.0, and 2.0, respectively.

As shown in Fig. 2, as expected, the canonical HRF model’s sensitivity falls off as mismatch between the generative and recovery responses is increased, especially for simulations with $<50\%$ overlap with the canonical model [see panel (e) and (f) (SNR: 0.3, 0.5, 1.0, and 2.0), for detail]. Below about a 50% overlap, the canonical model fails, but adding derivatives to this model (as random effects) recovers some of this loss in sensitivity and an increase of the AUC. The nonlinear model does not ever perform better than the canonical model with derivatives and, thus, is not recommended. At high SNR levels (1:0 and 2:0), the FIR (full deconvolution) and FIR-IRF (impulse response model) are able to handle large mismatches. However, this is not substantially better than the canonical model with derivatives at this 5-s task duration. As the SNR is lowered, the FIR and FIR-IRF models fail faster than the canonical models. This is due to the higher degrees of freedom in the FIR models. Thus, even though there is a mismatch between the simulated and the recovered HRF shapes, the performance of canonical models is better than the FIR models at lower SNR levels. The FIR-IRF response model is more stable than the FIR (full deconvolution/block-averaging) model at moderate SNR levels.

We note that a ROC analysis quantifies the sensitivity and specificity of the model to test null hypothesis. In the regression model of fNIRS, the null hypothesis is that the magnitude of the evoked signal change during the task is not statistically different from the baseline period. As demonstrated by these results, the rejection of this null hypothesis does not require a “perfect” alternative hypothesis. That is, it is not required to get the time window of the contrast or the shape of the response exactly right to still have the power to reject the null hypothesis.

4.2.

Sensitivity of Linear Models to Hemodynamic Response Function Mismatch per Task Duration

Although adding derivatives (as random effects) to the canonical HRF model recovers some of the loss of sensitivity due to systemic mismatch between the underlying and the recovered hemodynamic shapes at the short task duration (5 s), as shown in Fig. 2, this effect is negated as the duration of the task increases. Figure 3 shows the comparison of the area under the ROC curve (AUC) for all (448) possible configurations from eight different regression models, seven different SNR levels, and eight different task durations. It is noted the ROC analysis is repeated 900 times for every configuration. This figure has eight panels whereas each of these panels represents eight different task durations, that is, 1-, 2-, 5-, 10-, 15-, 20-, 25-, and 30-s durations, as shown in Figs. 3(a)–3(h), respectively. For every panel, it describes the AUC values ($y$-axis) for seven different SNR levels of 0.01, 0.1, 0.3, 0.5, 1.0, 2.0, and 5.0 ($x$-axis). As mentioned earlier, the range of AUC is 0.5 (worst) to 1 (best). Eight different regression models have been compared in every panel: (i) canonical (blue solid line), (ii) canonical+RE(deriv) (red dashed line), (iii) canonical + FE(deriv) (green dotted line), (iv) gamma function (black and dash-dot line), (v) boxcar function (magenta solid line), (vi) FIR-IRF (yellow dashed line), (vii) FIR (dark-green and dotted line), and (viii) nonlinear (light-blue and dash-dot line).

Fig. 3

Comparison of AUC from ROC analysis for all (eight) different regression models at various SNR levels ($x$-axis) from 900 repeated simulations. Eight different activation or task periods of 1, 2, 5, 10, 15, 20, 25, and 30 s are investigated in panels (a)–(h), respectively. In every panel, it shows the comparison of AUC values at seven SNR levels of 0.01, 0.1, 0.3, 0.5, 1.0, 2.0, and 5.0 for all (eight) different HRF models (i.e., canonical, canonical+RE(deriv), canonical+FE(deriv), gamma, boxcar, FIR-IRF, IRF, and nonlinear).

As the task duration increases, the exact timing of the transients of the response at the onset and recover are washed out by the steady-state behavior of the response during the duration. We found that above $\sim 10\text{-}\mathrm{s}$ task durations, the models all had about the same performance even when there was a mismatch between the true underlying HRF shape and the regression basis. For task durations of 2 to 5 s, the performance of the models differed mainly at low SNR levels ($<\mathrm{SNR}1.0$). Similar to the findings in Fig. 2, the FIR and FIR-IRF models lost their sensitivity quickest as SNR decreased, which was most pronounced at the shortest duration tasks (1 to 2 s). We noted that the boxcar model had a bit of odd behavior with a peak in its sensitivity at about SNR = 0.5–2 but then a loss in sensitivity at high SNRs. This was actually because this boxcar basis set had sharp edges to the response (e.g., a binary regressor). This resulted in discontinuities in the residual of the model after fitting and much of the response signal at the edges of the response are considered noise, particularly at higher SNR levels. As described previously, the boxcar model, which can only model an instantaneous change in the brain response between the baseline and the task period, was an unrealistic model. These discontinuities in the model residual at the transition points were more pronounced at higher SNRs.

4.3.

Receiver-Operatingor Characteristic Analysis for Short Task Duration

Figure 4 shows the comparison of the ROC analysis for 5-s task duration using six different SNR levels, that is, 0.1, 0.3, 0.5, 1.0, 2.0, and 5.0, as shown in Figs. 4(a)–4(f), respectively. Similar to the previous figure, eight different regression models have been compared in every panel. In detail [SNR: 0.1 in Fig. 4(a)], the ROC performance shows the AUC in descending order: gamma (AUC: 0.86), canonical+RE(deriv) (AUC: 0.82), FIR-IRF (AUC: 0.81), canonical (AUC: 0.80), nonlinear (AUC: 0.76), FIR (AUC: 0.70), boxcar (AUC: 0.63), and canonical+FE(deriv) (AUC: 0.60). As SNR levels increase ($>0.5$), the choice of basis used in the regression models has little effect on the ROC analysis [except boxcar function for reasons previously noted; see Figs. 4(c)–4(f)].

Fig. 4

ROC curves [TPR or sensitivity against FPR or (1-specificity)] for all (eight) different regression models at various SNR levels for 5-s task duration. Six different SNR levels of 0.1, 0.3, 0.5, 1.0, 2.0, and 5.0 are investigated in panels (a)–(f), respectively. In every panel, it shows the comparison of all (eight) different regression models: canonical, canonical+RE(deriv), canonical+FE(deriv), gamma, boxcar, FIR-IRF, IRF, and nonlinear.

5.1.

Response Magnitude versus Timing Hypotheses

We note that these findings only apply if the primary hypothesis of the study is testing the magnitude of the hemodynamic response, that is, if the hypothesis is comparing the evoked brain signal during a task to baseline or is comparing the amplitude of two tasks to each other. In this work, we did not look at how statistical tests can be used to examine if two tasks or groups have different evoked response timings from each other.

If two responses had different underlying timings, the use of a single basis set (including the FIR model using the same temporal window to define the contrast in both responses) would provide differing sensitivities to the two responses. Thus, rejection of the null hypothesis that the two responses had the same amplitude cannot distinguish between two responses with the same timing but different amplitudes or two responses with the same amplitude but differing timings. When viewed as a first-order Taylor-series expansion on the timing of the canonical basis set, the inclusion of derivatives to the canonical model can be used to test amplitude (testing the main coefficients) and first-order differences to the timing (using tests on the coefficients assigned to the derivatives). This however, would only test a null hypothesis that the two responses were not statistically different to the first order in terms of the onset or dispersion of the canonical shape. In this work, we had proposed a nonlinear regression model but found that these models required a minimum level of SNR to be useful as currently implemented. When testing magnitudes of the evoked response, we could not recommend the use of the nonlinear model in this current work. However, future work applying more robust nonlinear estimation methods, such as regularized methods, may improve the utility of these nonlinear models and allow for better testing of timing differences.

5.2.

Limitations of this Study

In this work, our results and conclusions are based on the testing of the null statistical hypothesis that the magnitude of the brain response during the task does not differ from the magnitude of the baseline period. In both our statistical tests and our simulations to generate ROC curves, we have assumed a linearly additive hemodynamic model. It is known that the hemodynamic response can be nonlinear, such as observed during long task durations demonstrating habituation effects or tasks with very short interstimulus intervals. In these cases where nonlinearities may exist, using a linear model to test and reject the null hypothesis is still statistically valid as a first-order test but will lose sensitivity (increase false-negatives) as the nonlinearity becomes more pronounced. Alternative methods such as the Volterra series³⁰ can be used to model such nonlinearities as better alternative hypothesis in these scenarios. An interesting future extension of this work would be to examine how these linear regressions begin to break down as the assumption of a linear additive hemodynamic response is violated. However, in this current work, our simulations are restricted only to the linear additive case.

A further conclusion of this work is that simpler (lower degrees of freedom) regression models are favorable particularly at lower SNR. We found that the increased sensitivity of these simple models often outweighs the loss in specificity due to the usage of a canonical model that is not a perfect match to the data. The results show that the canonical model performed well down to about a 50% overlap, which covers a fairly large variation in responses. However, this conclusion is not intended to say that the canonical model is necessarily the “best” for all data. For example, if there is a systematic bias in the shape of the brain response (e.g., infants may have different shape than adults), then using a slightly modified canonical model will recover some of the loss in specificity due to this mismatch. However, the conclusion of the paper is that, based on the AUC of the ROC analysis, it is better to have a slightly mismatched canonical shape than it is to use a model with high degrees of freedom.