{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Permutation t-test on source data with spatio-temporal clustering\n\nThis example tests if the evoked response is significantly different between\ntwo conditions across subjects. Here just for demonstration purposes\nwe simulate data from multiple subjects using one subject's data.\nThe multiple comparisons problem is addressed with a cluster-level\npermutation test across space and time.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n#          Eric Larson <larson.eric.d@gmail.com>\n#          Stefan Appelhoff <stefan.appelhoff@mailbox.org>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom numpy.random import randn\nfrom scipy import stats as stats\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.epochs import equalize_epoch_counts\nfrom mne.minimum_norm import apply_inverse, read_inverse_operator\nfrom mne.stats import spatio_temporal_cluster_1samp_test, summarize_clusters_stc"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Set parameters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nraw_fname = meg_path / \"sample_audvis_filt-0-40_raw.fif\"\nevent_fname = meg_path / \"sample_audvis_filt-0-40_raw-eve.fif\"\nsubjects_dir = data_path / \"subjects\"\nsrc_fname = subjects_dir / \"fsaverage\" / \"bem\" / \"fsaverage-ico-5-src.fif\"\n\ntmin = -0.2\ntmax = 0.3  # Use a lower tmax to reduce multiple comparisons\n\n# Setup for reading the raw data\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read epochs for all channels, removing a bad one\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw.info[\"bads\"] += [\"MEG 2443\"]\npicks = mne.pick_types(raw.info, meg=True, eog=True, exclude=\"bads\")\nevent_id = 1  # L auditory\nreject = dict(grad=1000e-13, mag=4000e-15, eog=150e-6)\nepochs1 = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    picks=picks,\n    baseline=(None, 0),\n    reject=reject,\n    preload=True,\n)\n\nevent_id = 3  # L visual\nepochs2 = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    picks=picks,\n    baseline=(None, 0),\n    reject=reject,\n    preload=True,\n)\n\n# Equalize trial counts to eliminate bias (which would otherwise be\n# introduced by the abs() performed below)\nequalize_epoch_counts([epochs1, epochs2])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Transform to source space\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fname_inv = meg_path / \"sample_audvis-meg-oct-6-meg-inv.fif\"\nsnr = 3.0\nlambda2 = 1.0 / snr**2\nmethod = \"dSPM\"  # use dSPM method (could also be MNE, sLORETA, or eLORETA)\ninverse_operator = read_inverse_operator(fname_inv)\nsample_vertices = [s[\"vertno\"] for s in inverse_operator[\"src\"]]\n\n# Let's average and compute inverse, resampling to speed things up\nevoked1 = epochs1.average()\nevoked1.resample(50, npad=\"auto\")\ncondition1 = apply_inverse(evoked1, inverse_operator, lambda2, method)\nevoked2 = epochs2.average()\nevoked2.resample(50, npad=\"auto\")\ncondition2 = apply_inverse(evoked2, inverse_operator, lambda2, method)\n\n# Let's only deal with t > 0, cropping to reduce multiple comparisons\ncondition1.crop(0, None)\ncondition2.crop(0, None)\ntmin = condition1.tmin\ntstep = condition1.tstep * 1000  # convert to milliseconds"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Transform to common cortical space\n\nNormally you would read in estimates across several subjects and morph\nthem to the same cortical space (e.g., fsaverage). For example purposes,\nwe will simulate this by just having each \"subject\" have the same\nresponse (just noisy in source space) here.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>Note that for 6 subjects with a two-sided statistical test, the minimum\n    significance under a permutation test is only\n    ``p = 1/(2 ** 6) = 0.015``, which is large.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_vertices_sample, n_times = condition1.data.shape\nn_subjects = 6\nprint(f\"Simulating data for {n_subjects} subjects.\")\n\n# Let's make sure our results replicate, so set the seed.\nnp.random.seed(0)\nX = randn(n_vertices_sample, n_times, n_subjects, 2) * 10\nX[:, :, :, 0] += condition1.data[:, :, np.newaxis]\nX[:, :, :, 1] += condition2.data[:, :, np.newaxis]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It's a good idea to spatially smooth the data, and for visualization\npurposes, let's morph these to fsaverage, which is a grade 5 source space\nwith vertices 0:10242 for each hemisphere. Usually you'd have to morph\neach subject's data separately (and you might want to use morph_data\ninstead), but here since all estimates are on 'sample' we can use one\nmorph matrix for all the heavy lifting.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Read the source space we are morphing to\nsrc = mne.read_source_spaces(src_fname)\nfsave_vertices = [s[\"vertno\"] for s in src]\nmorph_mat = mne.compute_source_morph(\n    src=inverse_operator[\"src\"],\n    subject_to=\"fsaverage\",\n    spacing=fsave_vertices,\n    subjects_dir=subjects_dir,\n).morph_mat\n\nn_vertices_fsave = morph_mat.shape[0]\n\n# We have to change the shape for the dot() to work properly\nX = X.reshape(n_vertices_sample, n_times * n_subjects * 2)\nprint(\"Morphing data.\")\nX = morph_mat.dot(X)  # morph_mat is a sparse matrix\nX = X.reshape(n_vertices_fsave, n_times, n_subjects, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we want to compare the overall activity levels in each condition,\nthe diff is taken along the last axis (condition). The negative sign makes\nit so condition1 > condition2 shows up as \"red blobs\" (instead of blue).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "X = np.abs(X)  # only magnitude\nX = X[:, :, :, 0] - X[:, :, :, 1]  # make paired contrast"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Find adjacency matrix\n\nFor cluster-based permutation testing, we must define adjacency relations\nthat govern which points can become members of the same cluster. While\nthese relations are rather obvious for dimensions such as time or frequency\nthey require a bit more work for spatial dimension such as channels or\nsource vertices.\n\nHere, to use an algorithm optimized for spatio-temporal clustering, we\njust pass the spatial adjacency matrix (instead of spatio-temporal).\nBut note that clustering still takes place along the\ntemporal dimension and can be\ncontrolled via the ``max_step`` parameter in\n:func:`mne.stats.spatio_temporal_cluster_1samp_test`.\n\nIf we wanted to specify an adjacency matrix for both space and time\nexplicitly we would have to use :func:`mne.stats.combine_adjacency`,\nhowever for the present case, this is not needed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(\"Computing adjacency.\")\nadjacency = mne.spatial_src_adjacency(src)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute statistic\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Note that X needs to be a multi-dimensional array of shape\n# observations (subjects) \u00d7 time \u00d7 space, so we permute dimensions\nX = np.transpose(X, [2, 1, 0])\n\n# Here we set a cluster forming threshold based on a p-value for\n# the cluster based permutation test.\n# We use a two-tailed threshold, the \"1 - p_threshold\" is needed\n# because for two-tailed tests we must specify a positive threshold.\np_threshold = 0.001\ndf = n_subjects - 1  # degrees of freedom for the test\nt_threshold = stats.distributions.t.ppf(1 - p_threshold / 2, df=df)\n\n# Now let's actually do the clustering. This can take a long time...\nprint(\"Clustering.\")\nT_obs, clusters, cluster_p_values, H0 = clu = spatio_temporal_cluster_1samp_test(\n    X,\n    adjacency=adjacency,\n    n_jobs=None,\n    threshold=t_threshold,\n    buffer_size=None,\n    verbose=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Selecting \"significant\" clusters\nAfter performing the cluster-based permutationt test, you may wish to\nselect the observed clusters that can be considered statistically\nsignificant under the permutation distribution. This can easily be\ndone using the code snippet below.\n\nHowever, it is crucial to be aware that a statistically significant\nobserved cluster does not directly translate into statistical\nsignificance of the channels, time points, frequency bins, etc. that\nform the cluster!\n\nFor more information, see the [FieldTrip tutorial](ft_cluster_).\n\n.. include:: ../../links.inc\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Select the clusters that are statistically significant at p < 0.05\ngood_clusters_idx = np.where(cluster_p_values < 0.05)[0]\ngood_clusters = [clusters[idx] for idx in good_clusters_idx]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualize the clusters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(\"Visualizing clusters.\")\n\n# Now let's build a convenient representation of our results, where consecutive\n# cluster spatial maps are stacked in the time dimension of a SourceEstimate\n# object. This way by moving through the time dimension we will be able to see\n# subsequent cluster maps.\nstc_all_cluster_vis = summarize_clusters_stc(\n    clu, tstep=tstep, vertices=fsave_vertices, subject=\"fsaverage\"\n)\n\n# Let's actually plot the first \"time point\" in the SourceEstimate, which\n# shows all the clusters, weighted by duration.\n\n# blue blobs are for condition A < condition B, red for A > B\nbrain = stc_all_cluster_vis.plot(\n    hemi=\"both\",\n    views=\"lateral\",\n    subjects_dir=subjects_dir,\n    time_label=\"temporal extent (ms)\",\n    size=(800, 800),\n    smoothing_steps=5,\n    clim=dict(kind=\"value\", pos_lims=[0, 1, 40]),\n)\n\n# We could save this via the following:\n# brain.save_image('clusters.png')"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}