{
"nbformat_minor": 0,
"nbformat": 4,
"cells": [
{
"execution_count": null,
"cell_type": "code",
"source": [
"%matplotlib inline"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"\n\n# Permutation t-test on source data with spatio-temporal clustering\n\n\nTests if the evoked response is significantly different between\nconditions across 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\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n# Eric Larson <larson.eric.d@gmail.com>\n# License: BSD (3-clause)\n\n\nimport os.path as op\n\nimport numpy as np\nfrom numpy.random import randn\nfrom scipy import stats as stats\n\nimport mne\nfrom mne import (io, spatial_tris_connectivity, compute_morph_matrix,\n grade_to_tris)\nfrom mne.epochs import equalize_epoch_counts\nfrom mne.stats import (spatio_temporal_cluster_1samp_test,\n summarize_clusters_stc)\nfrom mne.minimum_norm import apply_inverse, read_inverse_operator\nfrom mne.datasets import sample\n\nprint(__doc__)"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"Set parameters\n--------------\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"data_path = sample.data_path()\nraw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'\nevent_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'\nsubjects_dir = data_path + '/subjects'\n\ntmin = -0.2\ntmax = 0.3 # Use a lower tmax to reduce multiple comparisons\n\n# Setup for reading the raw data\nraw = io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"Read epochs for all channels, removing a bad one\n------------------------------------------------\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"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(raw, events, event_id, tmin, tmax, picks=picks,\n baseline=(None, 0), reject=reject, preload=True)\n\nevent_id = 3 # L visual\nepochs2 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n baseline=(None, 0), reject=reject, preload=True)\n\n# Equalize trial counts to eliminate bias (which would otherwise be\n# introduced by the abs() performed below)\nequalize_epoch_counts([epochs1, epochs2])"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"Transform to source space\n-------------------------\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"fname_inv = data_path + '/MEG/sample/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 or sLORETA)\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"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"Transform to common cortical space\n----------------------------------\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 7 subjects with a two-sided statistical test, the minimum\n significance under a permutation test is only p = 1/(2 ** 6) = 0.015,\n which is large.</p></div>\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"n_vertices_sample, n_times = condition1.data.shape\nn_subjects = 7\nprint('Simulating data for %d subjects.' % n_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]"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"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": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"fsave_vertices = [np.arange(10242), np.arange(10242)]\nmorph_mat = compute_morph_matrix('sample', 'fsaverage', sample_vertices,\n fsave_vertices, 20, subjects_dir)\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)"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"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": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"X = np.abs(X) # only magnitude\nX = X[:, :, :, 0] - X[:, :, :, 1] # make paired contrast"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"Compute statistic\n-----------------\n\nTo use an algorithm optimized for spatio-temporal clustering, we\njust pass the spatial connectivity matrix (instead of spatio-temporal)\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"source": [
"print('Computing connectivity.')\nconnectivity = spatial_tris_connectivity(grade_to_tris(5))\n\n# Note that X needs to be a multi-dimensional array of shape\n# samples (subjects) x time x space, so we permute dimensions\nX = np.transpose(X, [2, 1, 0])\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.\np_threshold = 0.001\nt_threshold = -stats.distributions.t.ppf(p_threshold / 2., n_subjects - 1)\nprint('Clustering.')\nT_obs, clusters, cluster_p_values, H0 = clu = \\\n spatio_temporal_cluster_1samp_test(X, connectivity=connectivity, n_jobs=1,\n threshold=t_threshold)\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]"
],
"outputs": [],
"metadata": {
"collapsed": false
}
},
{
"source": [
"Visualize the clusters\n----------------------\n\n"
],
"cell_type": "markdown",
"metadata": {}
},
{
"execution_count": null,
"cell_type": "code",
"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(clu, tstep=tstep,\n vertices=fsave_vertices,\n subject='fsaverage')\n\n# Let's actually plot the first \"time point\" in the SourceEstimate, which\n# shows all the clusters, weighted by duration\nsubjects_dir = op.join(data_path, 'subjects')\n# blue blobs are for condition A < condition B, red for A > B\nbrain = stc_all_cluster_vis.plot(hemi='both', views='lateral',\n subjects_dir=subjects_dir,\n time_label='Duration significant (ms)')\nbrain.save_image('clusters.png')"
],
"outputs": [],
"metadata": {
"collapsed": false
}
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"name": "python2",
"language": "python"
},
"language_info": {
"mimetype": "text/x-python",
"nbconvert_exporter": "python",
"name": "python",
"file_extension": ".py",
"version": "2.7.13",
"pygments_lexer": "ipython2",
"codemirror_mode": {
"version": 2,
"name": "ipython"
}
}
}
}