{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Repeated measures ANOVA on source data with spatio-temporal clustering\n\nThis example illustrates how to make use of the clustering functions\nfor arbitrary, self-defined contrasts beyond standard t-tests. In this\ncase we will tests if the differences in evoked responses between\nstimulation modality (visual VS auditory) depend on the stimulus\nlocation (left vs right) for a group of subjects (simulated here\nusing one subject's data). For this purpose we will compute an\ninteraction effect using a repeated measures ANOVA. The multiple\ncomparisons problem is addressed with a cluster-level permutation test\nacross 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#          Denis Engemannn <denis.engemann@gmail.com>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom numpy.random import randn\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.minimum_norm import apply_inverse, read_inverse_operator\nfrom mne.stats import (\n    f_mway_rm,\n    f_threshold_mway_rm,\n    spatio_temporal_cluster_test,\n    summarize_clusters_stc,\n)\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\"\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\")\n# we'll load all four conditions that make up the 'two ways' of our ANOVA\n\nevent_id = dict(l_aud=1, r_aud=2, l_vis=3, r_vis=4)\nreject = dict(grad=1000e-13, mag=4000e-15, eog=150e-6)\nepochs = 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)\nepochs.equalize_event_counts(event_id)"
      ]
    },
    {
      "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)\n\n# we'll only use one hemisphere to speed up this example\n# instead of a second vertex array we'll pass an empty array\nsample_vertices = [inverse_operator[\"src\"][0][\"vertno\"], np.array([], int)]\n\n#    Let's average and compute inverse, then resample to speed things up\nconditions = []\nfor cond in [\"l_aud\", \"r_aud\", \"l_vis\", \"r_vis\"]:  # order is important\n    evoked = epochs[cond].average()\n    evoked.resample(30).crop(0.0, None)\n    condition = apply_inverse(evoked, inverse_operator, lambda2, method)\n    #    Let's only deal with t > 0, cropping to reduce multiple comparisons\n    condition.crop(0, None)\n    conditions.append(condition)\n\ntmin = conditions[0].tmin\ntstep = conditions[0].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 them\nto the same cortical space (e.g. fsaverage). For example purposes, we will\nsimulate this by just having each \"subject\" have the same response (just\nnoisy in source space) here.\n\nWe'll only consider the left hemisphere in this tutorial.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_vertices_sample, n_times = conditions[0].lh_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, 4) * 10\nfor ii, condition in enumerate(conditions):\n    X[:, :, :, ii] += condition.lh_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 ICO source space\nwith vertices 0:10242 for each hemisphere. Usually you'd have to morph\neach subject's data separately, but here since all estimates are on\n'sample' we can use one morph 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 (just left hemisphere)\nsrc = mne.read_source_spaces(src_fname)\nfsave_vertices = [src[0][\"vertno\"], []]\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    smooth=20,\n).morph_mat\nmorph_mat = morph_mat[:, :n_vertices_sample]  # just left hemi from src\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 * 4)\nprint(\"Morphing data.\")\nX = morph_mat.dot(X)  # morph_mat is a sparse matrix\nX = X.reshape(n_vertices_fsave, n_times, n_subjects, 4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we need to prepare the group matrix for the ANOVA statistic. To make the\nclustering function work correctly with the ANOVA function X needs to be a\nlist of multi-dimensional arrays (one per condition) of shape: samples\n(subjects) \u00d7 time \u00d7 space.\n\nFirst we permute dimensions, then split the array into a list of conditions\nand discard the empty dimension resulting from the split using numpy squeeze.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "X = np.transpose(X, [2, 1, 0, 3])  #\nX = [np.squeeze(x) for x in np.split(X, 4, axis=-1)]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Prepare function for arbitrary contrast\nAs our ANOVA function is a multi-purpose tool we need to apply a few\nmodifications to integrate it with the clustering function. This\nincludes reshaping data, setting default arguments and processing\nthe return values. For this reason we'll write a tiny dummy function.\n\nWe will tell the ANOVA how to interpret the data matrix in terms of\nfactors. This is done via the factor levels argument which is a list\nof the number factor levels for each factor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "factor_levels = [2, 2]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we will pick the interaction effect by passing 'A:B'.\n(this notation is borrowed from the R formula language).\nAs an aside, note that in this particular example, we cannot use the A*B\nnotation which return both the main and the interaction effect. The reason\nis that the clustering function expects ``stat_fun`` to return a 1-D array.\nTo get clusters for both, you must create a loop.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "effects = \"A:B\"\n# Tell the ANOVA not to compute p-values which we don't need for clustering\nreturn_pvals = False\n\n# a few more convenient bindings\nn_times = X[0].shape[1]\nn_conditions = 4"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A ``stat_fun`` must deal with a variable number of input arguments.\n\nInside the clustering function each condition will be passed as flattened\narray, necessitated by the clustering procedure. The ANOVA however expects an\ninput array of dimensions: subjects \u00d7 conditions \u00d7 observations (optional).\n\nThe following function catches the list input and swaps the first and the\nsecond dimension, and finally calls ANOVA.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>For further details on this ANOVA function consider the\n          corresponding\n          `time-frequency tutorial <tut-timefreq-twoway-anova>`.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def stat_fun(*args):\n    # get f-values only.\n    return f_mway_rm(\n        np.swapaxes(args, 1, 0),\n        factor_levels=factor_levels,\n        effects=effects,\n        return_pvals=return_pvals,\n    )[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute clustering 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": [
        "# as we only have one hemisphere we need only need half the adjacency\nprint(\"Computing adjacency.\")\nadjacency = mne.spatial_src_adjacency(src[:1])\n\n# Now let's actually do the clustering. Please relax, on a small\n# notebook and one single thread only this will take a couple of minutes ...\npthresh = 0.005\nf_thresh = f_threshold_mway_rm(n_subjects, factor_levels, effects, pthresh)\n\n# To speed things up a bit we will ...\nn_permutations = 50  # ... run way fewer permutations (reduces sensitivity)\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    threshold=f_thresh,\n    stat_fun=stat_fun,\n    n_permutations=n_permutations,\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\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\nsubjects_dir = data_path / \"subjects\"\n# The brighter the color, the stronger the interaction between\n# stimulus modality and stimulus location\n\nbrain = stc_all_cluster_vis.plot(\n    subjects_dir=subjects_dir,\n    views=\"lat\",\n    time_label=\"temporal extent (ms)\",\n    clim=dict(kind=\"value\", lims=[0, 1, 40]),\n)\nbrain.save_image(\"cluster-lh.png\")\nbrain.show_view(\"medial\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, let's investigate interaction effect by reconstructing the time\ncourses:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inds_t, inds_v = [\n    (clusters[cluster_ind]) for ii, cluster_ind in enumerate(good_cluster_inds)\n][0]  # first cluster\n\ntimes = np.arange(X[0].shape[1]) * tstep * 1e3\n\nplt.figure()\ncolors = [\"y\", \"b\", \"g\", \"purple\"]\nevent_ids = [\"l_aud\", \"r_aud\", \"l_vis\", \"r_vis\"]\n\nfor ii, (condition, color, eve_id) in enumerate(zip(X, colors, event_ids)):\n    # extract time course at cluster vertices\n    condition = condition[:, :, inds_v]\n    # normally we would normalize values across subjects but\n    # here we use data from the same subject so we're good to just\n    # create average time series across subjects and vertices.\n    mean_tc = condition.mean(axis=2).mean(axis=0)\n    std_tc = condition.std(axis=2).std(axis=0)\n    plt.plot(times, mean_tc.T, color=color, label=eve_id)\n    plt.fill_between(\n        times, mean_tc + std_tc, mean_tc - std_tc, color=\"gray\", alpha=0.5, label=\"\"\n    )\n\nymin, ymax = mean_tc.min() - 5, mean_tc.max() + 5\nplt.xlabel(\"Time (ms)\")\nplt.ylabel(\"Activation (F-values)\")\nplt.xlim(times[[0, -1]])\nplt.ylim(ymin, ymax)\nplt.fill_betweenx(\n    (ymin, ymax), times[inds_t[0]], times[inds_t[-1]], color=\"orange\", alpha=0.3\n)\nplt.legend()\nplt.title(\"Interaction between stimulus-modality and location.\")\nplt.show()"
      ]
    }
  ],
  "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
}