{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# 2 samples permutation test on source data with spatio-temporal clustering\n\nTests if the source space data are significantly different between\n2 groups of subjects (simulated here 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# 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 scipy import stats as stats\n\nimport mne\nfrom mne import spatial_src_adjacency\nfrom mne.datasets import sample\nfrom mne.stats import spatio_temporal_cluster_test, summarize_clusters_stc\n\nprint(__doc__)"
      ]
    },
    {
      "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\"\nstc_fname = meg_path / \"sample_audvis-meg-lh.stc\"\nsubjects_dir = data_path / \"subjects\"\nsrc_fname = subjects_dir / \"fsaverage\" / \"bem\" / \"fsaverage-ico-5-src.fif\"\n\n# Load stc to in common cortical space (fsaverage)\nstc = mne.read_source_estimate(stc_fname)\nstc.resample(50, npad=\"auto\")\n\n# Read the source space we are morphing to\nsrc = mne.read_source_spaces(src_fname)\nfsave_vertices = [s[\"vertno\"] for s in src]\nmorph = mne.compute_source_morph(\n    stc,\n    \"sample\",\n    \"fsaverage\",\n    spacing=fsave_vertices,\n    smooth=20,\n    subjects_dir=subjects_dir,\n)\nstc = morph.apply(stc)\nn_vertices_fsave, n_times = stc.data.shape\ntstep = stc.tstep * 1000  # convert to milliseconds\n\nn_subjects1, n_subjects2 = 6, 7\nprint(f\"Simulating data for {n_subjects1} and {n_subjects2} subjects.\")\n\n#    Let's make sure our results replicate, so set the seed.\nnp.random.seed(0)\nX1 = np.random.randn(n_vertices_fsave, n_times, n_subjects1) * 10\nX2 = np.random.randn(n_vertices_fsave, n_times, n_subjects2) * 10\nX1[:, :, :] += stc.data[:, :, np.newaxis]\n# make the activity bigger for the second set of subjects\nX2[:, :, :] += 3 * stc.data[:, :, np.newaxis]\n\n#    We want to compare the overall activity levels for each subject\nX1 = np.abs(X1)  # only magnitude\nX2 = np.abs(X2)  # only magnitude"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute statistic\n\nTo use an algorithm optimized for spatio-temporal clustering, we\njust pass the spatial adjacency matrix (instead of spatio-temporal)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(\"Computing adjacency.\")\nadjacency = spatial_src_adjacency(src)\n\n#    Note that X needs to be a list of multi-dimensional array of shape\n#    samples (subjects_k) \u00d7 time \u00d7 space, so we permute dimensions\nX1 = np.transpose(X1, [2, 1, 0])\nX2 = np.transpose(X2, [2, 1, 0])\nX = [X1, X2]\n\n# Now let's actually do the clustering. This can take a long time...\n# Here we set the threshold quite high to reduce computation,\n# and use a very low number of permutations for the same reason.\nn_permutations = 50\np_threshold = 0.001\nf_threshold = stats.distributions.f.ppf(\n    1.0 - p_threshold / 2.0, n_subjects1 - 1, n_subjects2 - 1\n)\nprint(\"Clustering.\")\nF_obs, clusters, cluster_p_values, H0 = clu = spatio_temporal_cluster_test(\n    X,\n    adjacency=adjacency,\n    n_jobs=None,\n    n_permutations=n_permutations,\n    threshold=f_threshold,\n    buffer_size=None,\n)\n#    Now select the clusters that are sig. at p < 0.05 (note that this value\n#    is multiple-comparisons corrected).\ngood_cluster_inds = np.where(cluster_p_values < 0.05)[0]"
      ]
    },
    {
      "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 each cluster, where each\n#    cluster becomes a \"time point\" in the SourceEstimate\nfsave_vertices = [np.arange(10242), np.arange(10242)]\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\nbrain = stc_all_cluster_vis.plot(\n    \"fsaverage\",\n    hemi=\"both\",\n    views=\"lateral\",\n    subjects_dir=subjects_dir,\n    time_label=\"temporal extent (ms)\",\n    clim=dict(kind=\"value\", lims=[0, 1, 40]),\n)"
      ]
    }
  ],
  "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
}