Statistical inference

Here we will briefly cover multiple concepts of inferential statistics in an introductory manner, and demonstrate how to use some MNE statistical functions. :depth: 3

Hypothesis testing

Null hypothesis ^^^^^^^^^^^^^^^ From Wikipedia <https://en.wikipedia.org/wiki/Null_hypothesis>__:

In inferential statistics, a general statement or default position that
there is no relationship between two measured phenomena, or no
association among groups.

We typically want to reject a null hypothesis with some probability (e.g., p < 0.05). This probability is also called the significance level $\alpha$. To think about what this means, let's follow the illustrative example from [1]_ and construct a toy dataset consisting of a 40 x 40 square with a "signal" present in the center with white noise added and a Gaussian smoothing kernel applied.

The data averaged over all subjects looks like this:

In this case, a null hypothesis we could test for each voxel is:

There is no difference between the mean value and zero
($H_0 \colon \mu = 0$).

The alternative hypothesis, then, is that the voxel has a non-zero mean ($H_1 \colon \mu \neq 0$). This is a two-tailed test because the mean could be less than or greater than zero, whereas a one-tailed test would test only one of these possibilities, i.e. $H_1 \colon \mu \geq 0$ or $H_1 \colon \mu \leq 0$.

Parametric tests ^^^^^^^^^^^^^^^^ Let's start with a paired t-test, which is a standard test for differences in paired samples. Mathematically, it is equivalent to a 1-sample t-test on the difference between the samples in each condition. The paired t-test is parametric because it assumes that the underlying sample distribution is Gaussian, and is only valid in this case. This happens to be satisfied by our toy dataset, but is not always satisfied for neuroimaging data.

In the context of our toy dataset, which has many voxels ($40 \cdot 40 = 1600$), applying the paired t-test is called a mass-univariate approach as it treats each voxel independently.

"Hat" variance adjustment ~~~~~ The "hat" technique regularizes the variance values used in the t-test calculation [1]_ to compensate for implausibly small variances.

Non-parametric tests ^^^^^^^^^^^^^^^^^^^^ Instead of assuming an underlying Gaussian distribution, we could instead use a non-parametric resampling method. In the case of a paired t-test between two conditions A and B, which is mathematically equivalent to a one-sample t-test between the difference in the conditions A-B, under the null hypothesis we have the principle of exchangeability. This means that, if the null is true, we can exchange conditions and not change the distribution of the test statistic.

When using a paired t-test, exchangeability thus means that we can flip the signs of the difference between A and B. Therefore, we can construct the null distribution values for each voxel by taking random subsets of samples (subjects), flipping the sign of their difference, and recording the absolute value of the resulting statistic (we record the absolute value because we conduct a two-tailed test). The absolute value of the statistic evaluated on the veridical data can then be compared to this distribution, and the p-value is simply the proportion of null distribution values that are smaller.

In the case where n_permutations is large enough (or "all") so that the complete set of unique resampling exchanges can be done (which is $2^{N_{samp}}-1$ for a one-tailed and $2^{N_{samp}-1}-1$ for a two-tailed test, not counting the veridical distribution), instead of randomly exchanging conditions the null is formed from using all possible exchanges. This is known as a permutation test (or exact test).

Multiple comparisons

So far, we have done no correction for multiple comparisons. This is potentially problematic for these data because there are $40 \cdot 40 = 1600$ tests being performed. If we use a threshold p < 0.05 for each individual test, we would expect many voxels to be declared significant even if there were no true effect. In other words, we would make many type I errors (adapted from here <errors_>_):

.. rst-class:: skinnytable

+----------+--------+------------------+------------------+ | | Null hypothesis | | +------------------+------------------+ | | True | False | +==========+========+==================+==================+ | | | Type I error | Correct | | | Yes | False positive | True positive |

• Reject +--------+------------------+------------------+ | | | Correct | Type II error | | | No | True Negative | False negative | +----------+--------+------------------+------------------+

To see why, consider a standard $\alpha = 0.05$. For a single test, our probability of making a type I error is 0.05. The probability of making at least one type I error in $N_{\mathrm{test}}$ independent tests is then given by $1 - (1 - \alpha)^{N_{\mathrm{test}}}$:

To combat this problem, several methods exist. Typically these provide control over either one of the following two measures:

1. Familywise error rate (FWER) <fwer_>_ The probability of making one or more type I errors:

   .. math::

   \mathrm{P}(N_{\mathrm{type\ I}} >= 1 \mid H_0)

2. False discovery rate (FDR) <fdr_>_ The expected proportion of rejected null hypotheses that are actually true:

   .. math::

   \mathrm{E}(\frac{N_{\mathrm{type\ I}}}{N_{\mathrm{reject}}}
   \mid N_{\mathrm{reject}} > 0) \cdot
   \mathrm{P}(N_{\mathrm{reject}} > 0 \mid H_0)

We cover some techniques that control FWER and FDR below.

Bonferroni correction ^^^^^^^^^^^^^^^^^^^^^ Perhaps the simplest way to deal with multiple comparisons, Bonferroni correction <https://en.wikipedia.org/wiki/Bonferroni_correction>__ conservatively multiplies the p-values by the number of comparisons to control the FWER.

False discovery rate (FDR) correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Typically FDR is performed with the Benjamini-Hochberg procedure, which is less restrictive than Bonferroni correction for large numbers of comparisons (fewer type II errors), but provides less strict control of type I errors.

Non-parametric resampling test with a maximum statistic ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Non-parametric resampling tests can also be used to correct for multiple comparisons. In its simplest form, we again do permutations using exchangeability under the null hypothesis, but this time we take the maximum statistic across all voxels in each permutation to form the null distribution. The p-value for each voxel from the veridical data is then given by the proportion of null distribution values that were smaller.

This method has two important features:

1. It controls FWER.
2. It is non-parametric. Even though our initial test statistic (here a 1-sample t-test) is parametric, the null distribution for the null hypothesis rejection (the mean value across subjects is indistinguishable from zero) is obtained by permutations. This means that it makes no assumptions of Gaussianity (which do hold for this example, but do not in general for some types of processed neuroimaging data).

Clustering ^^^^^^^^^^ Each of the aforementioned multiple comparisons corrections have the disadvantage of not fully incorporating the correlation structure of the data, namely that points close to one another (e.g., in space or time) tend to be correlated. However, by defining the connectivity/adjacency/neighbor structure in our data, we can use clustering to compensate.

To use this, we need to rethink our null hypothesis. Instead of thinking about a null hypothesis about means per voxel (with one independent test per voxel), we consider a null hypothesis about sizes of clusters in our data, which could be stated like:

The distribution of spatial cluster sizes observed in two experimental
conditions are drawn from the same probability distribution.

Here we only have a single condition and we contrast to zero, which can be thought of as:

The distribution of spatial cluster sizes is independent of the sign
of the data.

In this case, we again do permutations with a maximum statistic, but, under each permutation, we:

1. Compute the test statistic for each voxel individually.
2. Threshold the test statistic values.
3. Cluster voxels that exceed this threshold (with the same sign) based on adjacency.
4. Retain the size of the largest cluster (measured, e.g., by a simple voxel count, or by the sum of voxel t-values within the cluster) to build the null distribution.

After doing these permutations, the cluster sizes in our veridical data are compared to this null distribution. The p-value associated with each cluster is again given by the proportion of smaller null distribution values. This can then be subjected to a standard p-value threshold (e.g., p < 0.05) to reject the null hypothesis (i.e., find an effect of interest).

This reframing to consider cluster sizes rather than individual means maintains the advantages of the standard non-parametric permutation test -- namely controlling FWER and making no assumptions of parametric data distribution. Critically, though, it also accounts for the correlation structure in the data -- which in this toy case is spatial but in general can be multidimensional (e.g., spatio-temporal) -- because the null distribution will be derived from data in a way that preserves these correlations.

.. sidebar:: Effect size

For a nice description of how to compute the effect size obtained
in a cluster test, see this
`FieldTrip mailing list discussion <ft_cluster_effect_size_>`_.

However, there is a drawback. If a cluster significantly deviates from the null, no further inference on the cluster (e.g., peak location) can be made, as the entire cluster as a whole is used to reject the null. Moreover, because the test statistic concerns the full data, the null hypothesis (and our rejection of it) refers to the structure of the full data. For more information, see also the comprehensive FieldTrip tutorial <ft_cluster_>_.

Defining the connectivity/neighbor/adjacency matrix ~~~~~~~~~~~ First we need to define our connectivity/neighbor/adjacency matrix. This is a square array (or sparse matrix) of shape (n_src, n_src) that contains zeros and ones to define which spatial points are connected, i.e., which voxels are adjacent to each other. In our case this is quite simple, as our data are aligned on a rectangular grid.

Let's pretend that our data were smaller -- a 3 x 3 grid. Thinking about each voxel as being connected to the other voxels it touches, we would need a 9 x 9 connectivity matrix. The first row of this matrix contains the voxels in the flattened data that the first voxel touches. Since it touches the second element in the first row and the first element in the second row (and is also a neighbor to itself), this would be::

[1, 1, 0, 1, 0, 0, 0, 0, 0]

:mod:sklearn.feature_extraction provides a convenient function for this:

In general the connectivity between voxels can be more complex, such as those between sensors in 3D space, or time-varying activation at brain vertices on a cortical surface. MNE provides several convenience functions for computing connectivity/neighbor/adjacency matrices (see the Statistics API <api_reference_statistics>).

Standard clustering ~~~~~~~ Here, since our data are on a grid, we can use connectivity=None to trigger optimized grid-based code, and run the clustering algorithm.

"Hat" variance adjustment ~~~~~ This method can also be used in this context to correct for small variances [1]_:

Threshold-free cluster enhancement (TFCE) ~~~~~~~~~ TFCE eliminates the free parameter initial threshold value that determines which points are included in clustering by approximating a continuous integration across possible threshold values with a standard Riemann sum <https://en.wikipedia.org/wiki/Riemann_sum>_ [2]. This requires giving a starting threshold start and a step size step, which in MNE is supplied as a dict. The smaller the step and closer to 0 the start value, the better the approximation, but the longer it takes.

A significant advantage of TFCE is that, rather than modifying the statistical null hypothesis under test (from one about individual voxels to one about the distribution of clusters in the data), it modifies the data under test while still controlling for multiple comparisons. The statistical test is then done at the level of individual voxels rather than clusters. This allows for evaluation of each point independently for significance rather than only as cluster groups.

We can also combine TFCE and the "hat" correction:

Visualize and compare methods

Let's take a look at these statistics. The top row shows each test statistic, and the bottom shows p-values for various statistical tests, with the ones with proper control over FWER or FDR with bold titles.

The first three columns show the parametric and non-parametric statistics that are not corrected for multiple comparisons:

• Mass univariate t-tests result in jagged edges.
• "Hat" variance correction of the t-tests produces less peaky edges, correcting for sharpness in the statistic driven by low-variance voxels.
• Non-parametric resampling tests are very similar to t-tests. This is to be expected: the data are drawn from a Gaussian distribution, and thus satisfy parametric assumptions.

The next three columns show multiple comparison corrections of the mass univariate tests (parametric and non-parametric). These too conservatively correct for multiple comparisons because neighboring voxels in our data are correlated:

• Bonferroni correction eliminates any significant activity.
• FDR correction is less conservative than Bonferroni.
• A permutation test with a maximum statistic also eliminates any significant activity.

The final four columns show the non-parametric cluster-based permutation tests with a maximum statistic:

• Standard clustering identifies the correct region. However, the whole area must be declared significant, so no peak analysis can be done. Also, the peak is broad.
• Clustering with "hat" variance adjustment tightens the estimate of significant activity.
• Clustering with TFCE allows analyzing each significant point independently, but still has a broadened estimate.
• Clustering with TFCE and "hat" variance adjustment tightens the area declared significant (again FWER corrected).

Statistical functions in MNE

The complete listing of statistical functions provided by MNE are in the Statistics API list <api_reference_statistics>, but we will give a brief overview here.

MNE provides several convenience parametric testing functions that can be used in conjunction with the non-parametric clustering methods. However, the set of functions we provide is not meant to be exhaustive.

If the univariate statistical contrast of interest is not listed here (e.g., interaction term in an unbalanced ANOVA), consider checking out the :mod:statsmodels package. It offers many functions for computing statistical contrasts, e.g., :func:statsmodels.stats.anova.anova_lm. To use these functions in clustering:

1. Determine which test statistic (e.g., t-value, F-value) you would use in a univariate context to compute your contrast of interest. In other words, if there were only a single output such as reaction times, what test statistic might you compute on the data?
2. Wrap the call to that function within a function that takes an input of the same shape that is expected by your clustering function, and returns an array of the same shape without the "samples" dimension (e.g., :func:mne.stats.permutation_cluster_1samp_test takes an array of shape (n_samples, p, q) and returns an array of shape (p, q)).
3. Pass this wrapped function to the stat_fun argument to the clustering function.
4. Set an appropriate threshold value (float or dict) based on the values your statistical contrast function returns.

Parametric methods provided by MNE ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

• :func:mne.stats.ttest_1samp_no_p Paired t-test, optionally with hat adjustment. This is used by default for contrast enhancement in paired cluster tests.

• :func:mne.stats.f_oneway One-way ANOVA for independent samples. This can be used to compute various F-contrasts. It is used by default for contrast enhancement in non-paired cluster tests.

• :func:mne.stats.f_mway_rm M-way ANOVA for repeated measures and balanced designs. This returns F-statistics and p-values. The associated helper function :func:mne.stats.f_threshold_mway_rm can be used to determine the F-threshold at a given significance level.

• :func:mne.stats.linear_regression Compute ordinary least square regressions on multiple targets, e.g., sensors, time points across trials (samples). For each regressor it returns the beta value, t-statistic, and uncorrected p-value. While it can be used as a test, it is particularly useful to compute weighted averages or deal with continuous predictors.

Non-parametric methods ^^^^^^^^^^^^^^^^^^^^^^

• :func:mne.stats.permutation_cluster_test Unpaired contrasts with connectivity.

• :func:mne.stats.spatio_temporal_cluster_test Unpaired contrasts with spatio-temporal connectivity.

• :func:mne.stats.permutation_t_test Paired contrast with no connectivity.

• :func:mne.stats.permutation_cluster_1samp_test Paired contrasts with connectivity.

• :func:mne.stats.spatio_temporal_cluster_1samp_test Paired contrasts with spatio-temporal connectivity.

References

.. [1] Ridgway et al. 2012, "The problem of low variance voxels in statistical parametric mapping; a new hat avoids a 'haircut'", NeuroImage. 2012 Feb 1;59(3):2131-41.

.. [2] Smith and Nichols 2009, "Threshold-free cluster enhancement: addressing problems of smoothing, threshold dependence, and localisation in cluster inference", NeuroImage 44 (2009) 83-98.

.. include:: ../../links.inc