Beyond linear neural envelope tracking: a mutual information approach

The most common approach to estimating probability values is based on the histogram-based algorithm. Here, the data are first partitioned into a set of discrete bins. The probability of observing a discrete bin value is then obtained by dividing the number of occurrences in a single bin by the total number of data points. The probability values can then be plugged into the marginal and joint entropy formulas.

Binning most commonly occurs according to the equipartition principle, where bins have equal occupancy (e.g. if using ten bins, each bin contains 10% of the data). An important consideration is how many bins one should choose for constructing the histogram. In theory, one should not have too few bins, otherwise, there will be little sensitivity to the distribution, and the relationship between two variables can be lost. On the other hand, if one has too many bins, estimation errors for the joint distribution can be made hence generality is lost. The latter case can be seen as a form of overfitting, as the distribution will become sensitive to noise components in the data. Therefore, the bin parameter can induce either an under- (too few bins) or overestimation (too many bins) to the true MI. Previous neural tracking studies have chosen the bin parameter rather arbitrarily with little justification or empirical data to support their choice (Kaufeld et al 2020, Zan et al 2020).

Instead of partitioning the data into discrete bins, several continuous alternatives exist to estimate the probability distribution. In the kernel density estimation (KDE) approach, each datapoint is convolved with a (most commonly Gaussian) kernel. The probability density function is then obtained by summing all kernel values and normalizing so the integral equals one. In the KDE approach, researchers need to pre-specify the width of the kernel (i.e. the bandwidth), which is the equivalent critical choice to the number of bins in the histogram-based approach (Venelli 2010). The bandwidth parameter determines the smoothness of the obtained probability distribution. A bandwidth that is too large leads to over-smoothing, hence specificity is lost, while a bandwidth that is too small leads to under-smoothing, hence generality is lost. We currently are not aware of any study applying the KDE-based method to neural envelope tracking. The KDE method does not solve the issue of bin-selection, and additionally suffers from a high computational cost.

Correctly choosing the number of bins or the bandwidth parameter is critical for constructing the probability distribution. Several rules of thumb exist, such as the Freedman–Diaconis rule (Freedlan and Diaconis 1981) for histogram-based distributions and the Silverman's rule (Silverman 1986) for KDE estimations. Cross-validation approaches have also been proposed in the past (Sain et al 1994). With these approaches, the obtained parameter value is data-specific. Yet, it is crucial to fix this parameter value to enable comparisons across participants and channels since different parameter values introduce different baseline MI values (Vinh et al 2009). An alternative approach that bypasses the need to set the bins or bandwidth parameter is to assume that the data are normally distributed. Under this assumption, the MI between two Gaussian variables is determined based on the covariance matrices:

$$\begin{equation} I(X;Y) = \frac{1}{2ln2}\Biggl[\frac{|\sum_{X}||\sum_{Y}|}{|\sum_{XY}|}\Biggr] \end{equation} \tag{ 3 }$$

where $\sum_{X}$ and $\sum_{Y}$ are the covariance matrices of variables X and Y, and $\sum_{XY}$ is the covariance matrix for the joint variable. Straight lines indicate that the determinant of the covariance matrices is calculated. The resulting algorithm does not require any parameter values. However, we must rely on the strong assumption that the data are normally distributed.

Recently, Ince et al (2017) introduced the concept of Gaussian copula MI for neuroimaging research. The method uses the Gaussian definition of MI but transforms the data so it meets the normality assumption. The procedure first implies that values for the individual variables are ranked on a scale from 0 to 1, obtaining the cumulative density functions (CDF's). By computing the inverse standard normal CDF of the variables separately, the data distributions are transformed to perfect standard Gaussians. With these steps, the statistical relationship between both variables (i.e. the 'copula') is preserved. After data transformations to perfect standard Gaussians, the parametric MI estimate in equation (3) can be applied. Crucial, the MI in this method provides a lower bound to the true MI, i.e. erroneous high values cannot occur. This is because the joint Gaussian distribution has the maximum entropy for a variable with a given mean and covariance matrix (Cover and Thomas 1991). For a more in-depth explanation, we refer to appendix. Over recent years, the Gaussian copula MI has been applied in numerous neural envelope tracking papers (Giordano et al 2017, Daube et al 2019, Coopmans et al 2022, Perez et al 2022).

The Gaussian copula has a particular advantage over the other MI derivations in multidimensional spaces. While the histogram- and KDE-based approaches severely suffer from the curse of dimensionality and computational complexity, the Gaussian copula method can efficiently handle multidimensional variables. After transforming each individual dimension to a standard Gaussian, the formula in equation (3) can be applied. For the context of neural envelope tracking, this may open perspectives to calculating the MI between the speech envelope and multiple channels combined:

$$\begin{align*} I(Envelope;\;Channel\;1,\;Channel\;2,\;\ldots\;Channel\;N). \end{align*}$$

This could be considered an information-theoretic analog to the linear backward model, as we combine multiple EEG channels to determine the statistical relationship with the envelope.

In the present study, we compare the Gaussian copula approach with the probability density estimation approach in the context of neural envelope tracking. We restrict the comparison to the most popular approach, i.e. the histogram-based method, to demonstrate the effect of arbitrarily setting a parameter (i.e. the number of bins) to estimate the probability distribution. Next, we validate the most suitable MI technique against the linear forward and backward model. To this end, we apply the measures to two datasets containing EEG data of participants listening to continuous speech. The first dataset includes a large group study of healthy young participants listening to a natural story. With this dataset, we compare the temporal and spatial patterns of the TMIF and TRF in-depth, discuss the properties of multivariate MI in relation to the linear backward model, and investigate whether non-linear components are present in the neural response to the speech envelope. The second dataset includes normal hearing younger and older adults listening to continuous speech to demonstrate the application of MI in a group study and highlight certain interesting differences between MI and its linear counterpart.

For both datasets, the speech envelope was extracted using a gammatone filter bank (Sø ndergaard et al 2012), with 28 channels spaced by 1 equivalent rectangular bandwidth and center frequencies from 50 Hz until 5000 Hz. The envelopes were extracted from each subband by taking each sample's absolute value and raising it to the power of 0.6 (Stevens 1955, Biesmans et al 2015). The resulting 28 subband envelopes were averaged to obtain one single envelope. The envelope was then downsampled to 512 Hz to decrease processing time. Next, the envelope was highpass-filtered at 0.5 Hz (transition band 0.45–0.5 Hz) and lowpass-filtered at 8 Hz (transition band 8–8.8 Hz) using a least squares filter of order 2000 and compensated for the group delay. After filtering, the envelope was normalized and further downsampled to 128 Hz.

The EEG signals were first downsampled to 512 Hz to decrease processing time. Eye-blink artifacts were removed using a multi-channel Wiener filter (Somers et al 2018) and the EEG was referenced to common average. The data were then highpass-filtered at 0.5 Hz and lowpass-filtered at 8 Hz using the same least squares filter used in the envelope extraction method and compensated for the group delay. Next, normalization and additional downsampling to 128 Hz were applied.

In the linear forward model, the EEG is predicted as a function of the speech envelope. This approach estimates the TRF, a kernel that describes the neural response at different latencies to the speech envelope. We used the Eelbrain toolbox (Brodbeck 2020) to estimate the TRF for each EEG channel separately using the boosting algorithm (David et al 2007). To avoid overfitting, a ten-fold cross-validation was used (ten equally long folds, eight for training, one for validation and one for testing). The TRF covered an integration window from −200 to 600 ms, with a 50 ms Hamming window basis. Training stopped when the mean squared error on the validation data stopped decreasing. Next, the TRFs were used to predict the EEG-response in the test fold by convolving it with the speech envelope. We then calculated the Spearman correlation between the predicted EEG and the actual EEG per channel, yielding a measure of neural tracking. The TRFs were averaged across all folds to interpret the temporal and spatial patterns.

In the linear backward model, the speech envelope is reconstructed from the neural data. A linear decoder attributes weights to each EEG channel and a time-shifted response of each channel. The reconstruction of the speech envelope $\hat{s}(t)$ is obtained through:

$$\begin{equation} \hat{s}(t) = \sum_{n=1}^{N}\sum_{\tau}g(n, \tau)R(t + \tau, n) \end{equation} \tag{ 4 }$$

where g is the linear decoder as a function of the number of electrodes N and the integration window τ, (−200 to 600 ms). Matrix R contains the shifted neural responses of all channels (with time index t). The decoder was calculated by solving $g = (RR^t)^{-1}(Rs^t)$ with s the speech envelope. To prevent overfitting, ridge regularization and ten-fold cross-validation (ten equally long folds, nine used to train the decoder and one used for testing) were applied. In the test fold, the reconstructed envelope was correlated with the actual envelope yielding a measure of neural tracking. Results were averaged across all folds. For a more in-depth explanation of linear models and ridge regression, we refer to Crosse et al (2021).

The MI between the envelope and the EEG was calculated using the histogram-based and Gaussian copula methods. The MI in the histogram-based approach is obtained by applying equation (2) in section 1.1. Histograms were constructed according to the equipartition principle, where bins have varying widths yet equal observation probabilities (hence, we obtain uniform distributions). We partitioned the data into several numbers of bins to demonstrate the effect of this parameter on the results. The Gaussian copula MI was calculated by applying the summarized data transformations in section 1.1 and the Gaussian definition of MI (see equation (3)). The analyses for the Gaussian copula MI were performed with the 'GCMI'-toolbox (Ince et al 2017). For both MI measures, the temporal evolution of the MI (i.e. the TMIF) is obtained by shifting the EEG relative to the speech envelope from 200 ms prestimulus to 600 ms poststimulus and recalculating the MI at each sample.

The MI between all individual channels and the envelope was calculated to derive spatial interpretations. This measure will be compared to the linear forward model. The Gaussian copula approach further allows calculating the MI between the speech envelope and all channels combined by computing the covariance matrices of multidimensional data. This multivariate MI will be compared to the analogous linear counterpart, i.e. the backward model, where the envelope is reconstructed as a weighted sum of multiple channels.

All methods were subject to permutation testing to quantify the meaningfulness of the derived values. We created stationary noise that matched the spectrum of the EEG responses. This was repeated 1000 times for each channel and each participant individually. For each repetition, we calculated the MI, the prediction accuracy (linear forward model) or the reconstruction accuracy (linear backward model). The significance level was then determined as the 95th percentile of the neural tracking result of all repetitions.

For the multivariate MI, permutations were created from noise matching the spectrum of the envelope rather than noise matching the spectrum of the EEG. This was done to prevent dependencies between EEG channels from disappearing.

We investigated whether nonlinear relationships are captured by MI in the neural response to the speech envelope, both for individual channels and for the multivariate MI (all channels combined). To this end, we first applied linear models to remove all linear components in the data. The predicted EEG was subtracted from the actual EEG for the forward model. Next, the MI was calculated between the speech envelope and the residual EEG. For the multivariate MI, the linear backward model reconstructed the envelope, and this was subtracted from the actual envelope. The MI between the residual envelope and the 64-dimensional EEG data (64 channels) was then calculated. The results were compared to the significance level.

We used ordinary least squares without regularization to ensure that all linear components that the linear models could find were removed entirely from the data. Weights and correlations were then exactly equal to 0 (within machine precision) in a second linear fit, performed as a sanity check.

We will refer to a central and a parieto-occipital channel selection in the results section to plot the results in channels that exhibit the same polarity and that are known to contribute to auditory speech perception (Lesenfants et al 2019). These respectively include channels [FC5, FC3, FC1, C5, C3, C1, FC6, FC4, FC2, C6, C4, C2, FCz, Cz] and [TP7, CP5, P5, P7, P9, PO7, PO3, O1, TP8, CP6, P6, P8, P10, PO8, PO4, O2, Iz, Oz, POz]. Since MI measures are strictly positive, we combine both clusters into a single large cluster. We will also refer to a control region (i.e. frontal channels) where no envelope tracking is expected. This control region includes channels [Fpz, Fp1, AF7, F7, Fp2, AF8, F8, AFz, AF3, AF4]. The channel selections are depicted in supplementary figure 1.

3.1.1.1. For different channels

The histogram-based approach requires a prespecified number of bins. Figure 2 illustrates the effect of this parameter value on the TMIF for Dataset 1. When taking two bins (figure 2(a)), the variability across channels is minimal, but sensitivity to the full relationship is lost. For example, the left-lateralized temporal peak at 380 ms, present with four bins (figure 2(b)), ten bins (figure 1(a)) and the linear TRF (figure 1(c)), is largely suppressed. With four bins, specificity is preserved compared to the pattern that results from a larger number of bins (10, figure 1(a)). But when the number of bins is too high (20), the spatial pattern changes drastically: the MI in frontal channels dramatically increases compared to the other channels (figure 2(c)). Higher MI values for the frontal channels are even present before 0 ms, where the brain cannot yet process speech. We further observe that the topography of the significance level of frontal channels, obtained through calculating the MI between the envelope and spectrally matched EEG-shaped noise, is higher than the significance level of our channel selection (see figure 2(d)), irrespective of the number of bins. The histogram-based method thus introduces a bias in the results: some channels systematically exhibit higher MI values than others by chance. The enhanced envelope tracking in frontal channels, observed when using too many bins, does not reflect a true effect.

Zoom In Zoom Out Reset image size

The issue of unequal baseline levels across channels has implications, for example, when interpreting the spatial distribution of MI. We demonstrate this in figure 2(e). The relative MI values for frontal channels and our channel selection were calculated by subtracting its average from the global average (all channels). We plot these relative measures as a function of the number of bins. For illustration, we interpret the case of four and ten bins, two parameter values previously set in the field of envelope tracking (Gross et al 2013, Zan et al 2020). If four bins were selected, one might conclude that there is relatively higher neural tracking in our channel selection. Conversely, MI is higher for the frontal channels when using ten bins. Thus, the number of bins that the researcher arbitrarily sets influences the spatial distribution of MI. Compared to significance-level, figure 2(e) shows that from eight bins on, the envelope is no longer tracked significantly for the mean of all channels (gray line crosses the mean MI value).

Unequal baseline MI across channels originates from differences in power in the lowest frequencies. We illustrate this bias with 30 s of EEG data of one participant in supplementary figure 2. Channels with relatively high power in the low frequencies (0.5–1 Hz) have a more slowly oscillating behavior and remain in the same bin for a more extended consecutive period of time (see Panels (a) and (b)). As a result, when no true relationship is present, the joint entropy with a random variable has a lower expected value: periods of high occurrences for certain bin combinations alternate with low to zero occurrences for other combinations in the joint distribution (Panel (c)). Lower joint entropy is directly related to higher MI (see equation (2)). Therefore, channels with high power in the low frequencies have a higher expected MI by chance. Across channels, we report a strong correlation between power in the low frequencies (0.5–1 Hz) and the significance level (Spearman's r = 0.75, p $\lt$ 0.001).

The MI based on the Gaussian copula, on the other hand, entails a transformation of the marginal distributions to standard Gaussians. The calculation of MI is then based on the covariance matrices. Hence, it does not require prespecified parameter values. For the Gaussian copula TMIF, we observe a pattern similar to the histogram-based approach, but with the additional benefit that the frontal channels, where no true relationship between the EEG and the envelope is expected, are effectively suppressed (see figure 1(b)).

3.1.1.2. For different stimuli

The histogram-based MI thus suffers from a bias induced by power in the low frequencies. Yet, this bias arises in the stimulus as well. We illustrate this with the two stimuli we consider in this paper: the story 'De Kleine Zeemeermin' ('DKZ') for Dataset 1, and 'Milan' for Dataset 2. The envelope of the latter story has higher power in the low frequencies, visualized by the power spectrum (figure 3(a)). As a result, the latter story has higher baseline MI values (obtained through calculating the MI with EEG-shaped noise). This bias again exacerbates with a higher number of bins. The Gaussian copula MI, on the other hand, is much more robust to this issue: baseline neural tracking measures are approximately equal for both stimuli (figure 3(c)).

Zoom In Zoom Out Reset image size

The linear TRF models the brain's response to the stimulus, hence polarity naturally occurs. The TMIF, on the other hand, is a non-negative metric of neural tracking at each timestep. Nonetheless, the temporal and spatial patterns of both measures correspond well (see figure 1). A negative peak in the TRF is also represented as a positive peak in the TMIF. Hence, TRF topographies with high positive (red) and high negative (blue) values are both high positive (red) in the TMIF topographies. The prediction accuracy in the linear model, which can be obtained by correlating the actual and the predicted EEG per channel, reflects neural tracking in the linear approach. The topography of these correlations resembles the MI (see figure 1(c)).

We quantified the degree of correspondence between the two measures for the temporal patterns and neural tracking. For the temporal patterns, we computed the correlation between each corresponding channel in the TMIF and the absolute value of the TRF at the single-subject level. We used an integration window of 0–500 ms, including the peaks and latencies of interest. The group results are visualized in figure 4(a). The sample represents the mean Spearman correlation across channels per participant. Across participants, correlations are moderate to strong, with an average of r = 0.57.

Zoom In Zoom Out Reset image size

The prediction accuracy quantifies neural tracking in the linear forward model, resulting in a single value per participant. To compare this measure to MI, we extracted the mean MI over the same integration window as the linear model, i.e. 0–500 ms. Both neural tracking measures were averaged over frontal and temporal channels. We fitted a regression over all participants (see figure 4(b)) and report a Spearman correlation of r = 0.94 (p $\lt$ 0.001). In conclusion, the MI and the linear forward model exhibit great correspondence.

To check whether MI captures nonlinear components in the single-channel response to the speech envelope, we first filtered out the linear components in the data by subtracting the predicted from the actual EEG (see Methods). The TMIF between the residual EEG and the speech envelope was then calculated and interpreted against significance.

The TMIF for our channel selection is depicted in figure 5(a). Although the temporal pattern exhibits a reasonable evolution over time (i.e. with peak MI around 120 ms post-onset and a left-lateralized parietal topography), the TMIF did not reach significance on the group level. Not a single participant exhibited significant tracking of nonlinear components for the mean in a 0–400 ms integration window (figure 5(b)). For the peak MI value, the TMIF reached significance in 27 out of 64 participants. The magnitude of MI for solely nonlinear components dropped with a factor of ≈ 1000 compared to the analyses including both linear and nonlinear components. This shows that, although there is significant envelope tracking for some participants for nonlinear components in the data, the added value of MI over the linear forward model is small and, in this case, not significant on the group level.

Zoom In Zoom Out Reset image size

Finally, we compare the linear forward model to MI in a group study. In a recent comparison, Gillis et al (2023) found higher neural envelope tracking for younger listeners compared to older adults. We replicated these results (Wilcoxon signed-rank test, W = 95, p = 0.039, figure 6(a)). Yet, enhanced neural tracking for younger listeners was not represented in the linear TRF amplitudes (a non-parametric cluster-based permutation test found no significant higher amplitudes for younger adults, figure 6(a)). This is in contrast to the TMIF, where group differences do emerge in the spatio-temporal domain (a left parieto-occipital, p = 0.007, and a right temporo-central, p = 0.019) and on average over channels and latencies (W = 107, p = 0.005, figure 6(b)).

Zoom In Zoom Out Reset image size

A popular feature of the linear TRF is that it allows for peak latency analyses. In this dataset, a shift in the latency of the first peak was observed as age increased (Gillis et al 2023). We plotted this trend for the parieto-occipital (positive peak) and central (negative peak) channels in figure 6(c) (Pearson's r = −0.27, p = 0.048; r = −0.39, p = 0.004, respectively). The correlation in peak latency is also present for the TMIF (channels combined; Pearson's r = −0.31, p = 0.020). In conclusion, we show that the TMIF can handle peak latency analyses with statistical relationships comparable to the TRF.

Section 3.2 reported that the single-channel MI analysis did not capture significant nonlinear relationships with the envelope. Here, we investigate whether such relationships do arise for the multivariate MI. Linear components captured by the backward model were removed by subtracting the reconstructed envelope from the actual envelope. Next, we calculated the TMIF between the residual envelope and the 64-dimensional EEG data.

The results are depicted in figure 8. The temporal pattern reaches significance at the group level and peaks at 120 ms post envelope (figure 8(a)). The mean MI (integration window 0–400 ms) is significant for 56 out of 64 participants, and the peak value reaches significance for 62 participants. Compared to the analysis that included linear components, MI dropped by a factor 10. This drop is largely reduced compared to the MI for single channels, where a factor of 1000 was observed (section 3.2). We conclude that multivariate MI has significant added value over the linear backward model. Furthermore, we conclude that multivariate MI analysis is a more powerful tool for detecting nonlinearities in the data than single-channel MI analyses.

Zoom In Zoom Out Reset image size

An interaction with linear relations can partially explain the spread across participants in the nonlinear MI. Figure 8(c) illustrates a negative correlation between linear and nonlinear relations in the data. Lower reconstruction accuracies for the linear backward model are strongly related to higher nonlinear MI values (Spearman's r = −0.56, p $\lt$ 0.001). This shows that nonlinear relations between the EEG and the envelope tend to be more strongly represented when linear relations are less prominent and demonstrates the added value of a model that captures all kinds of variance in the data: there is no loss of information with nonlinear methods.

The MI between the EEG and the envelope is a relative measure and requires permutation testing. Although the bits scale provides a meaningful effect size (Friston 2012), many factors can impact the magnitude of MI. These include the length of the stimulus and the EEG recording (shorter length results in higher MI), the sampling rate (lower sampling results in higher MI), and the power spectrum (higher power in lowest frequency bands result in higher MI by chance). Permutation testing must be performed by computing a null distribution that controls for these covariates. Random permutations of EEG- or envelope-shaped noise (shaped with the same power spectrum as the original signal) are both valid choices (Nichols and Holmes 2002, Crosse et al 2021). The 95th percentile of permutations can be taken as the alpha = 0.05 significance level. With the TMIF, individual time points can then be compared to the significance level. In line with our expectations, brain latencies following the envelope were significantly tracked while latencies preceding the envelope were not (see figures 1, 6 and 7). In conclusion, neural tracking in terms of bits should not be interpreted as an absolute value, but rather as a relative value compared to a significance level. Note that permutation tests are also required for linear methods (Combrisson and Jerbi 2015, Crosse et al 2021).

As discussed, the Gaussian copula MI can be applied to either individual EEG channels or a combination of channels (i.e. the multivariate MI). The former approach has the advantage that it allows for an interpretation of the spatial distribution of MI. In contrast, the latter method is more informative as it determines the multivariate relationship between the envelope and all EEG signals combined. Note, however, that there will be an upper limit to the number of dimensions (i.e. channels) that one can use in the multivariate MI, which will strongly depend on the amount of data. Future work needs to address the interplay between the amount of data and the maximum number of dimensions. In the present study, the envelope was significantly tracked for a combination of 64 channels with 15 min of data. We suggest that researchers relate the obtained outcome to the significance level. When neural tracking is not significant with a certain number of dimensions, researchers can consider applying data reduction techniques such as denoising source separation (Särelä and Valpola 2005) or a multiway canonical correlation analysis (de Cheveigné et al 2019). A deep canonical correlation analysis (Andrew et al 2013) may be the preferred option, as it preserves the nonlinear nature of the data.

The most straightforward approach to estimating the probability distribution is based on the histogram-based algorithm. Here, the data are first partitioned into a set of discrete bins. The probability of observing a discrete value is obtained by dividing the number of occurrences in a single bin by the total number of data points. The probability values can then be plugged into the formulas of marginal and joint entropy.

The binned approach is appealing due to its simplicity. However, there are several considerations to take into account. First, the data discretization can be performed via two main strategies: either bins are set with equal bin width ('equidistant'), or bins are set to have equal occupancy ('equipartition', e.g. if using ten bins, each bin contains 10% of the data). The equidistant approach estimates the true probability distribution, while the equipartition binning transforms the data to a uniform distribution.

The equipartition approach is the preferred option for two main reasons. First, it is robust to sampling errors. In the equidistant approach, bins with low occupancy can receive too low probability due to insufficient sampling, not reflecting the true probability curve. In contrast, the equipartition binning assigns equal observation probabilities to all bins. And second, the equidistant approach facilitates MI comparisons over variables since the entropy for all individual variables is equal for a certain number of bins (i.e. all variables have uniform distributions). Hence, the MI will not depend on the marginal entropies, but will be fully encapsulated by the joint entropy (see equation (1)). This is preferred as the joint entropy determines the relationship between both variables. By contrast, the distributions and marginal entropies differ across variables in the equidistant approach. As a result, the MI is a complex interplay between the marginal and the joint entropy.

A second consideration is the number of bins to choose for constructing the histogram. In theory, one should not have too few bins, otherwise, there will be little sensitivity to the distribution, and the relationship between two variables can be lost. On the other hand, if one has too many bins, estimation errors for the joint distribution can be made hence generality is lost. The latter case can be seen as a form of overfitting, as the distribution will become sensitive to noise components in the data. Therefore, the bin parameter can induce either an under- (too few bins) or overestimation (too many bins) to the true MI. Previous neural tracking studies have chosen the bin parameter rather arbitrarily with little justification or empirical data to support their choice.

Instead of partitioning the data into discrete bins, several continuous alternatives exist to estimate the probability distribution. In the KDE approach, each datapoint is convolved with a (most commonly Gaussian) kernel. The probability density function is then obtained by summing all kernel values and normalizing so the integral equals one. In the KDE approach, researchers need to pre-specify the width of the kernel (i.e. the bandwidth), which is the equivalent critical choice to the number of bins in the histogram-based approach (Venelli 2010). The bandwidth parameter determines the smoothness of the obtained probability distribution. A bandwidth set too large leads to over smoothing, hence specificity is lost, while a bandwidth too small leads to under smoothing, hence generality is lost. A recent neural envelope tracking paper calculated MI based on KDE estimations (Keshavarzi et al 2021), but the bandwidth parameter was not specified.

Correctly choosing the number of bins or the bandwidth parameter is critical for constructing the probability distribution. Several rules of thumb exist, such as the Freedman–Diaconis rule (Freedlan and Diaconis 1981) for histogram-based distributions and the Silverman's rule (Silverman 1986) for KDE estimations. Cross-validation approaches have also been proposed in the past (Sain et al 1994). With these approaches, the obtained parameter value is data-specific. Yet, it is crucial to fix this parameter value to enable comparisons across participants and channels since different parameter values introduce different baseline MI values (Vinh et al 2009). An alternative approach that bypasses the need to set the bins or bandwidth parameter is to assume that the data are normally distributed. Under this assumption, the MI between two Gaussian variables is determined based on the covariance matrices:

$$\begin{equation} I(X;Y) = \frac{1}{2ln2}\Biggl[\frac{|\sum_{X}||\sum_{Y}|}{|\sum_{XY}|}\Biggr] \end{equation} \tag{ 4 }$$

where $\sum_{X}$ and $\sum_{Y}$ are the covariance matrices of variables X and Y, and $\sum_{XY}$ is the covariance matrix for the joint variable. Straight lines indicate that the determinant of the covariance matrices is calculated. The resulting algorithm does not require any parameter values. However, we must rely on the strong assumption that the data is normally distributed.

Recently, Ince et al (2017) introduced the concept of Gaussian Copula MI for neuroimaging research. The method uses the Gaussian definition of MI but transforms the data, so it meets the normality assumption. The procedure first implies that the individual variables are ranked on a scale from 0 to 1, obtaining the CDF's. These ranked variables are then uniformly distributed. The joint distribution of these two variables, i.e. the joint CDF, is expressed in function of the marginal CDF's and a function that links these two variables (Sklar 1959). This function is called the copula, and it describes the statistical relationship between the two variables. The entropy of the joint CDF is equal to the sum of the entropy of the marginal CDF's plus the entropy of the copula (Ma and Sun 2011):

$$\begin{equation} H(X,Y) = H(X) + H(Y) + H(c) \end{equation} \tag{ 5 }$$

where H(c) is the entropy of the copula. When plugging this formula for the joint entropy into the general formula of MI (equation (1)), we obtain:

$$\begin{equation} I(X;Y) = -H(c). \end{equation} \tag{ 6 }$$

MI is equal to the negative entropy of the copula (Ma and Sun 2011). We thus no longer require the entropy of the marginals, as the MI is fully captured by the copula (i.e. the statistical relationship between both variables). As a result, the MI can be calculated by estimating the entropy of the copula. Yet, this again results in the problem of accurately estimating probability values for the copula. Instead, Ince et al (2017) suggest exploiting the fact that the MI does not depend on the individual variables (see equation (6)). Therefore, we can transform the marginal distributions in any way we see fit as long as the copula is preserved. By computing the inverse standard normal CDF of the variables separately, the data distributions are transformed to perfect standard Gaussians, and the copula is preserved. As a result, the parametric MI estimate in equation (4) can be applied.

In summary, the Gaussian Copula TMIF in the context of neural envelope tracking is obtained by applying the following steps:

• (1)  
  Rank the speech envelope and the neural time series of a single channel on a 0–1 scale, obtaining the marginal CDFs as uniform distributions.
• (2)  
  Transform the marginal CDFs to standard Gaussians by computing the inverse standard normal CDF.
• (3)  
  Create the joint variable by stacking the envelope and the neural time series as a matrix sized time by 2.
• (4)  
  Compute the covariance matrices and calculate MI according to the formula in equation (4).
• (5)  
  Repeat for all channels and time lags.

The transformations of the data distributions are depicted in figure 2.

Zoom In Zoom Out Reset image size

This algorithm is computationally efficient and makes no assumptions on the distributions of the variables. Crucial to mention, the MI in this method cannot be overestimated, i.e. erroneous high values cannot occur. This is because the joint Gaussian distribution has the maximum entropy for a variable with a given mean and covariance matrix (Cover and Thomas 1991). As MI is the negative copula entropy, the MI derived from a Gaussian copula is a lower bound to the true MI (Calsaverini and Vicente 2009). For a more in-depth explanation of Gaussian copula MI, we refer to Ince et al (2017). Over the past recent years, the Gaussian copula MI has been applied in numerous neural envelope tracking papers (Giordano et al 2017, Daube et al 2019, Coopmans et al 2022, Perez et al 2022).

The Gaussian copula has a particular advantage over the other MI derivations in multidimensional spaces. While the histogram- and KDE-based approaches severely suffer from the curse of dimensionality and computational complexity, the Gaussian copula method can efficiently handle multidimensional variables. After transforming each individual dimension to a standard Gaussian, the formula in equation (4) can be applied. For the context of neural envelope tracking, this may open perspectives to calculating the MI between the speech envelope and multiple channels combined: I(Envelope; Channel 1, Channel 2, $\ldots$ Channel N). This could be considered an information-theoretic analog to the linear backward model.