2.1.

Solving the Linear Model

The linear regression model [Eq. (1)] is typically solved by use of the Gauss–Markov expression, which yields the solution to the least-squares fit of the data ($Y$) under the assumption of spherical, uncorrelated, and homoscedastic noise. These assumptions will be discussed in Sec. 3. Under these assumptions, the best linear unbiased estimator of the coefficients is given by

The uncertainty (covariance) of these coefficients is the given by the expression

Based on the value of the coefficients and their covariance, statistical $T$-tests can be used to test if the coefficients are nonzero or if multiple coefficients differ from each other (e.g., the difference is nonzero). The $T$-test on these coefficients is then given by the expression

3.1.

The Noise in Functional Near-Infrared Spectroscopy Is Correlated

Typical fNIRS data have strong physiological noise due to signals such as cardiac, respiratory, and blood pressure changes. These physiological signals are typically much slower than the sample rate of an fNIRS system. This implies that the noise in fNIRS data has a colored structure and serial correlations. Colored noise refers to the fact that specific temporal (or spatial) frequencies are present in the noise spectrum. In contrast, white noise is characterized by a flat (uniform) noise spectrum. Colored noise results in serial correlations in the error term $ϵ$, which means that each sample point is not independent from the others and depends on its self history. Thus, the effective degrees of freedom of the linear model are often much less than the number of sample points. This can also result in a potential bias in the estimate of the coefficients if the noise is not mean zero. While fMRI also has colored noise due to the same systemic physiology, since the sample rate of fNIRS (up to tens to hundreds of Hz) is typically higher than fMRI (generally around 0.5 Hz), and since fNIRS is more sensitive to superficial physiology in the scalp, this noise tends to affect fNIRS measurements more strongly. Figure 1(a) shows an example of an fNIRS-recorded oxy-hemoglobin signal demonstrating physiological noise from both cardiac and blood pressure oscillations. In Fig. 1(c), the autocorrelation of this data is demonstrated. As can be seen, the autocorrelation of the raw signal remains significantly above chance well beyond the 20 s of lag time shown in the figure. This shows that two data points sampled 10 s apart in this file still share on average 60% of the information (down to 20% after 20 s) and cannot be considered independent measurements. Thus, in this sample data in Fig. 1, although there are over 7000 time points (700 s at 10 Hz) in this data trace, there are far fewer than this number of independent degrees of freedom, thus corrections must be made to account for this in the statistical model and estimate of the false-discovery rates for reports of brain activity, which we will further detail in Sec. 6 of this work.

Fig. 1

Example of noise in fNIRS data. This figure shows an example of an fNIRS time course of oxy-hemoglobin. The zoomed view (a, inset) shows physiological noise due to cardiac and blood pressure changes. Typical motion artifacts (a) are clearly visible on this data as well. (b) A histogram noise distribution of the data and the best-fit ideal distribution (dotted red line). This is shown in logarithm scale for the $y$-axis to highlight the deviations at the tails of this distribution (highlighted in red), which constitute heavy-tailed noise associated with the motion artifacts. (c) The temporal autocorrelation of this raw data along with the autocorrelation after filtering and AR prewhitening (described in Sec. 4), which are used to reduce these serially correlated errors.

3.2.

Functional Near-Infrared Spectroscopy Data Is Not Independent Across Measurement Channels

In fNIRS, correlations are also found between different source–detector pairs due to the typically low spatial frequencies of systemic and superficial physiological signals. Motion artifacts can also give rise to strong spatial noise structures. In addition, measurements are correlated between optical wavelengths and estimates of oxy- and deoxy-hemoglobin. Oxy- and deoxy-hemoglobin are partially correlated variables, which are together related to underlying changes in neural activity and the resulting blood flow and oxygen metabolism changes in the brain. Thus, unlike fMRI, fNIRS analysis generally involves a multivariate problem with partially correlated variables. More formally, oxy- and deoxy-hemoglobin are nonlinear correlated as related to the underlying dynamics of blood vessel dilation and oxygen transport. However, this topic is considered beyond the scope of this discussion. Spatial covariance and spatially structured noise have an impact on the assumptions of group level, region-of-interest, and image reconstruction procedures used in fNIRS.

3.3.

Functional Near-Infrared Spectroscopy Data Exhibits Heteroscedasticity

Noise is said to have homoscedasticity when all the samples arise from the same noise distribution. In other words, the variance of the noise is uniform. In general, this is not true for fNIRS data, which is therefore said to be heteroscedastic. Over the time dimension, motion artifacts often confound fNIRS data. A motion artifact, which may often be considerably larger than physiological noise, constitutes a heavy-tail to the noise distribution. These artifacts arise from the movement of the fNIRS sensors on the head and transient changes in the contact of these sensors with the scalp. A heavy-tailed noise distribution is one in which there are a few strong noise components compared to the normal noise background. These strong noise terms are outliers in the noise model. Such noise is heteroscedastic because these outlier samples are from a different distribution than the rest of the noise. The example data trace shown in Fig. 1(a) shows several motion artifacts. The variance of this oxy-hemoglobin signal evaluated over the time window free of motion-artifacts from 200 to 260 s [inset; (${\sigma }^2=3.5\mu {\mathrm{M}}^2$)] is a factor of 25 less than the period from 300 to 360 s (${\sigma }^2=95\mu {\mathrm{M}}^2$), which has several large motion artifacts. Thus, this data demonstrates noise from at least two distinct distributions. In Fig. 1(b), the noise distribution of this data is shown on a logarithm scale. The best fit of the normal distribution for this data is shown in red. Heavy-tailed noise elevates the tails of this distribution and is highlighted in red/shaded on this figure, which deviates from the assumed normal statistical model. We note that not all motion artifacts are necessarily outliers from the normal distribution of noise, however, and since these specific types of artifacts do not necessarily violate the homoscedastic assumption of the noise model, these are not of specific concern to the solution to the linear model. Although such noise would lower effect sizes and reduce the power to detect brain activity, homoscedastic noise would not necessarily result in an inaccurate statistical model or uncontrolled false-discovery.

In addition, the noise over the spatial dimension across fNIRS channels (source–detector pairs) is also typically nonuniform. Specific source–detector pairs or measurements at a particular wavelength may vary considerably in the level of noise compared with the other measurements. This noise is largely determined by how well the sensors are placed on the head, contact with the skin, the presence of hair under the fNIRS probe, and so on. In practice, it is not uncommon to have an order of magnitude difference or more in the signal-to-noise ratio between channels with good and poor contract on the head. Thus, the noise across fNIRS channels should not be assumed to be normal and exhibits heteroscedasticity. Due to the small number of fNIRS measurement channels compared with voxels in an fMRI scan, the central limit theorem generally does not apply in the fNIRS spatial domain. In the analysis methods used in the program NIRS-SPM,¹⁴ noise covariance estimates are performed by pooling the data across channels with the iterative restricted maximal likelihood (ReML) method. The ReML method is used to estimate both the noise covariance and the autoregressive (AR) coefficients within the covariance structure used in the prewhitening procedure (discussed in Sec. 4 of this work). Because this method obtains only one set of hyperparameters for oxy- and deoxy-hemoglobin that define the covariance structure applied to the whole channel set (e.g., across channels with varied qualities of optode-head contact), this estimate can be skewed by having one or a few channels that have very different noise from the other channels, as demonstrated in Sec. 6 of this work. As a result, the use of ReML in the NIRS-SPM software requires an additional preprocessing step, described in the user manual, to remove “bad” channels prior to computing the ReML covariance models. This step is not generally required in fMRI analysis, since the spatial noise of fMRI is more uniform and can otherwise be spatially smoothed (spatial precoloring). However, in fNIRS data, the noise per optical source–detector pair can vary significantly, thus the same noise model is often not sufficient for analysis (we will return to this concept in Sec. 6 of this paper). In addition, spatial heteroscedastic noise across the fNIRS measurements is of concern in group-level analyses, where corresponding channels are analyzed across subjects, and in analysis pooling data from multiple channels on the probe such as region-of-interest or image reconstruction methods.

3.4.

Functional Near-Infrared Spectroscopy Data May Be Nonergodic and Nonstationary

fNIRS noise due to systemic physiology may also appear to change over the length of a scan or study as a reflection of very slow underlying physiological changes (e.g., due to posture, blood pressure, and so on). Of more concern in analysis, however, systemic noise may also change corresponding to the stimulus task. During physical tasks such as motor movement (e.g., finger tapping) or tasks such as walking or exercise, changes in respiration, heart-rate, and blood pressure may be altered in a task-dependent manner. Thus, the noise during these time blocks may have different properties (mean, variance, or spectral content). This type of noise is said to be nonstationary or nonergodic. Specifically, because systemic physiology can change in a task-dependent manner, the noise structure of the “baseline” can be different from the “task,” which violates assumptions of the $T$-task and the estimation of brain activity under these conditions. The issue of nonstationary noise in fNIRS has been addressed in a number of publications in the context of dynamic estimation models such as a Kalman filter,^15–18 but since there is no simple solution to dealing with these noise errors, we will simply note that this is a potential source of error that could cause fNIRS signals to violate the common assumptions in most statistical models used to date in analysis, and very little is known about how much of an effect this has on the sensitivity and specificity of fNIRS methods.

In Fig. 1(a), we show an example of an fNIRS data trace of oxy-hemoglobin from an experimental study (unpublished). This is a fairly representative trace of the quality of data that we often see in our own studies. This particular time course was from a child imaging study and shows a few more motion artifacts than one might see in an adult fNIRS study since often for adults, the head cap can be placed tighter and there is a bit better subject compliance. While a zoomed view (inset figure) shows clear physiological signals related to cardiac and respiratory fluctuations, the full time course also shows larger but less frequent noise fluctuations from the motion of the fNIRS head cap relative to the head (indicated in the figure). These motion artifacts are almost an order of magnitude larger than the variance of other parts of the data (e.g., physiological noise shown in the zoomed view) and much greater than expected from physiological fluctuations. These infrequent, larger fluctuations deviate from uniform normal statistical distribution and create a heavy-tailed noise distribution, as shown in Fig. 1(b), which shows the best fit of a normal distribution to this data and highlights the deviations due to these outliers at the tail of this distribution (red/shaded). A metric such as the Tukey’s bisquare weight (also shown in panel B as the dotted blue line) can be used to estimate the probability of these values at the tails of the normal distribution as outliers, which is iteratively used in robust regression methods, discussed in Sec. 5 of this paper. As further shown in Fig. 1(c), the structures in the noise also generate serially correlated errors in which the current sample point depends on previous samples and the assumption of independent noise is violated. In this example trace, the raw data maintains autocorrelation for more than 20 s of lag (200 sample points at 10 Hz).

4.1.

Noise Prewhitening

Prewhitening of a regression model is a process in which the data and model are reweighted to remove colored noise from the residual $ϵ$. The prewhitening expression is given by

As a further note, prewhitening is not the same as prefiltering. Prefiltering is a preprocessing step in which a high-pass filter is applied only to the data. For example, a $1/60\mathrm{Hz}$ high-pass filter is often used in fNIRS to remove slow drifts from measurements. As demonstrated in Fig. 1, some of the observed temporal autocorrelation can be reduced by high-pass filtering, but this may introduce anticorrelations in the data, as shown in Fig. 1(c) by the negative correlation at a 5-s lag and beyond. In comparison to prewhitening, prefiltering is applied only to the measurements ($Y$ and $ϵ$ terms), but not the design matrix $X$. Since this filtering is not applied to $X$, the estimate of $\beta$ is no longer unbiased. Thus, prefiltering should only be used if one has confidence that the frequencies being removed were not relevant to model of the data (e.g., $X\beta$). Since this is not the case in prewhitening, which is applied to both sides of the expression, prewhitening is generally recommended over prefiltering to avoid a biased estimate of $\beta$.

The whitening matrix $W$ can be specified either based on prior expectations about the frequency content of the noise in the data or empirically by iteratively fitting the model and examining the residual. The prior case requires the user to specify a specific frequency filter. For example, a 0 to $1/60\mathrm{Hz}$ discrete cosine transform (DCT) might be used to whiten slow drifts in the data akin to a high-pass filter. In this case, $W$ is constructed from the residual forming matrix of the filtering regressors [e.g., $W=I-D·{(D^T·D)}^{-1}·D^T$, where $I$ is an identity matrix and $D$ is the matrix constructed from the DCT regressors]. The advantage of using filters based on prior expectations on the noise structure is that it is fast and computationally simple. In fMRI, e.g., iterative estimation of a whitening filter on a per-voxel basis would be time consuming. Prior defined whitening filters such as the DCT filter also have the advantage of having a more straightforward interpretation in terms of control over exactly which frequencies are affected in the model. Prewhitening using a DCT filter [or similar, such as the minimum descriptive length (MDL) wavelet method²⁵] is used in programs such as SPM¹² and NIRS-SPM.¹⁴ Predefined whitening filters, however, may not be entirely efficient. For example, if there were no low frequencies in the noise, one is unnecessarily discarding information from the data, but in principle does not change the estimate of $\beta$, although its uncertainty will likely increase because less information is being used in the estimate.

An alternative is an iterative approach that estimates $W$ directly from the residual noise in the data, which can offer a more efficiently whitened noise spectrum, but is more computationally intensive. In the iterative methods, the linear model is first fit using no whitening (e.g., $W=\text{identity}$). The residual noise is then estimated from the model fit via Eq. (6), and the whitening matrix is found based on the current residual. The whitened model is refit and continued until convergence. A common approach to estimating the whitening filter is to fit an autoregressive model to the residual.^19,20,22,24 An autoregressive model will remove periodic noise structures (e.g., the dependence of a current time point to its self history). Additional components such as moving-average (ARMA) or integrative (ARIMA) models can also be used to additionally model drift and nonstationarities in the noise. The whitening model (AR or other) is fit to the residual of the linear regression, then used to construct the whitening matrix $W$. As a slight variation of this approach, in the SPM¹² and NIRS-SPM¹⁴ programs, the AR model is currently implemented as an off-diagonal structure in the covariance components used in the iterative ReML method and estimated using an expectation-maximization algorithm. Here, the same AR model is applied across all spatial channels and is often used in addition to the DCT or MDL prewhitening methods.

Figure 1(c) shows the effect of prewhitening on the temporal autocorrelation of a sample NIRS data file. As demonstrated in this figure, AR models can further reduce this autocorrelation, but require a fairly high model order to completely whiten the noise. In this case of 10 Hz data, a model order of $n=33$ was determined by Bayesian information criteria for optimal whitening. This model order is notably larger than typically used for similar corrections in fMRI research. In comparison, the default AR model order used in the ReML prewhitening step of NIRS-SPM¹⁴ is $n=1$, which, as we will show in Sec. 5 of this paper, corrects some but often not all of these serially correlated errors. The first order AR model ($n=1$), shown in Fig. 1(c), reduces the significant autocorrelation to below a 5-s lag but introduces some ringing in the autocorrelation and negative correlations in the structure of the noise.

4.2.

Noise Precoloring

Noise precoloring is an alternative to prewhitening that has been used in fMRI analysis including SPM¹² and NIRS-SPM.¹⁴ In noise precoloring, a known noise color spectrum is applied to the data to mask over the natural but unknown noise.²⁶ Although there is now more colored structure in the noise, it has a known profile and can be accounted for by changing the statistical model. Similar to prewhitening, the temporal precoloring matrix is applied on both sides of Eq. (1) to yield an expression similar to Eq. (5). While prewhitening is akin to a high-pass filter, precoloring is similar to a low-pass filter and imposes additional autocorrelation structure on the noise that overpowers the intrinsic (unknown) structure.

Precoloring assumes homoscedastic noise and is adversely affected by the presence of heavy-tailed and outlier noise points. By imposing strong autocorrelation structures on the noise, outlier noise points can be smeared, affecting large portions of the data. This is similar to the effect a low-pass filter would have on an outlier spike in the data (e.g., due to motion). Thus, precoloring in the presence of heavy-tailed noise distributions can increase the overall noise of the data and be counterproductive. This blurring also reduces the effectiveness of outlier regression (robust regression) methods, as described in Sec. 5. Thus, in general, prewhitening is recommended over precoloring for fNIRS data.

7.1.

The Deconvolution Finite-Impulse Response Model

In a deconvolution model, the design matrix encodes the timing of the stimulus events but makes no prior assumptions about the shape of the response. The design matrix is constructed from shifted columns of a binary vector (e.g., 0’s and 1’s) constructed from the onset of the stimulus. Each column of $X$ contains a shifted version of the stimulus onsets. For example, the first (unshifted) column models the response at $\text{time}=0$ from when the stimulus onset occurred. The second column (shifted by $+1$) models the response one time point after the onset of the stimulus. Columns shifted backward (e.g., $\text{shift}=-1$) model the response before the stimulus onset (e.g., prebaseline periods), and so on. Thus, in this design matrix, there is one regressor for each time point over the window in which one is trying to estimate the response. As an example, if the data was sampled twice a second (2 Hz), the design model for a 30-s-long task with a 14-s-long poststimulus (recovery) period would have 88 ($30\mathrm{s}\times 2\mathrm{Hz}+14\mathrm{s}\times 2\mathrm{Hz}$) columns. When the model was estimated, there would be 88 corresponding coefficients $\beta$ that model the amplitude of the response over this entire time period from stimulus onset until 14 s after the end of the task. If these coefficients were plotted as a time course, one would hope to see a response rising with the stimulus, staying elevated for around 30 s, then falling back to baseline around 6 to 10 s after the task ended. In this model, the entire time course for the response is estimated, and the only assumption that is made about the response is that it has recovered before the end of the estimation window. One limitation of the full deconvolution model, however, is that all events of the same type must have the same duration since the same estimation time window is used for all events of the same type.

In order to define the statistical contrast, a contrast vector must then be specified to compute the average contrast over some time window via Eq. (4). Continuing the same example, if the estimated response were 44 s long and estimated at 2 Hz, then the $T$-statistic for the window from 10 to 30 s could be found by setting the contrast vector $c$ to 1’s for the 20th to 60th entries and 0’s for everything else, as defined in Eq. (3).

A variation on the full deconvolution model is the impulse response model. In the impulse response model, the design matrix is constructed with knowledge of the duration of events, and the model of the response is the convolution of the task duration with an impulse response. An impulse response is the theoretical response to a short task only one time point in duration. To construct this design matrix, each column of the matrix is a binary vector of 1’s during the duration of the stimulus and 0’s for baseline periods. In this case, the columns are shifted to model the impulse response. In the same example discussed above, each single column of the design matrix would have 30-s blocks ($30\mathrm{s}\times 2\mathrm{Hz}=60\text{points}$) for each stimulus event. This would be shifted 28 times to model the 14-s window of the impulse response ($14\mathrm{s}\times 2\mathrm{Hz}=28$). This model has fewer coefficients than the full deconvolution model, but makes the assumption that the response is sustained at the same level for the example 30-s duration of the task. The advantage of this model, however, is that it can be used to describe experiments where stimulus blocks of the same type have different durations (e.g., in the case of a self-paced study) since the duration of each block is encoded in the design matrix and the impulse response function does not vary with the stimulus duration.

7.2.

The Canonical Model

The canonical model is widely used in fMRI studies. Instead of modeling the shape of the response by allowing each time point in the response to be described by a separate coefficient $\beta$, the canonical model uses a prior hypothesis on the shape of the response. This idealized response shape is then used to construct the design matrix such that only one coefficient (or a few) is needed to model the amplitude of the response. In this model, the shape of the curve is fixed. This shape represents a prior on the expected response. A number of different variations of this model have been proposed.³⁰