{ "nbformat_minor": 0, "nbformat": 4, "cells": [ { "execution_count": null, "cell_type": "code", "source": [ "%matplotlib inline" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "\n\n# Spatiotemporal permutation F-test on full sensor data\n\n\nTests for differential evoked responses in at least\none condition using a permutation clustering test.\nThe FieldTrip neighbor templates will be used to determine\nthe adjacency between sensors. This serves as a spatial prior\nto the clustering. Significant spatiotemporal clusters will then\nbe visualized using custom matplotlib code.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# Authors: Denis Engemann \n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.axes_grid1 import make_axes_locatable\nfrom mne.viz import plot_topomap\n\nimport mne\nfrom mne.stats import spatio_temporal_cluster_test\nfrom mne.datasets import sample\nfrom mne.channels import read_ch_connectivity\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'\nevent_id = {'Aud_L': 1, 'Aud_R': 2, 'Vis_L': 3, 'Vis_R': 4}\ntmin = -0.2\ntmax = 0.5\n\n# Setup for reading the raw data\nraw = mne.io.read_raw_fif(raw_fname, preload=True)\nraw.filter(1, 30, l_trans_bandwidth='auto', h_trans_bandwidth='auto',\n filter_length='auto', phase='zero')\nevents = mne.read_events(event_fname)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Read epochs for the channel of interest\n---------------------------------------\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "picks = mne.pick_types(raw.info, meg='mag', eog=True)\n\nreject = dict(mag=4e-12, eog=150e-6)\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n baseline=None, reject=reject, preload=True)\n\nepochs.drop_channels(['EOG 061'])\nepochs.equalize_event_counts(event_id)\n\ncondition_names = 'Aud_L', 'Aud_R', 'Vis_L', 'Vis_R'\nX = [epochs[k].get_data() for k in condition_names] # as 3D matrix\nX = [np.transpose(x, (0, 2, 1)) for x in X] # transpose for clustering" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Load FieldTrip neighbor definition to setup sensor connectivity\n---------------------------------------------------------------\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "connectivity, ch_names = read_ch_connectivity('neuromag306mag')\n\nprint(type(connectivity)) # it's a sparse matrix!\n\nplt.imshow(connectivity.toarray(), cmap='gray', origin='lower',\n interpolation='nearest')\nplt.xlabel('{} Magnetometers'.format(len(ch_names)))\nplt.ylabel('{} Magnetometers'.format(len(ch_names)))\nplt.title('Between-sensor adjacency')" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Compute permutation statistic\n-----------------------------\n\nHow does it work? We use clustering to `bind` together features which are\nsimilar. Our features are the magnetic fields measured over our sensor\narray at different times. This reduces the multiple comparison problem.\nTo compute the actual test-statistic, we first sum all F-values in all\nclusters. We end up with one statistic for each cluster.\nThen we generate a distribution from the data by shuffling our conditions\nbetween our samples and recomputing our clusters and the test statistics.\nWe test for the significance of a given cluster by computing the probability\nof observing a cluster of that size. For more background read:\nMaris/Oostenveld (2007), \"Nonparametric statistical testing of EEG- and\nMEG-data\" Journal of Neuroscience Methods, Vol. 164, No. 1., pp. 177-190.\ndoi:10.1016/j.jneumeth.2007.03.024\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# set cluster threshold\nthreshold = 50.0 # very high, but the test is quite sensitive on this data\n# set family-wise p-value\np_accept = 0.001\n\ncluster_stats = spatio_temporal_cluster_test(X, n_permutations=1000,\n threshold=threshold, tail=1,\n n_jobs=1,\n connectivity=connectivity)\n\nT_obs, clusters, p_values, _ = cluster_stats\ngood_cluster_inds = np.where(p_values < p_accept)[0]" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Note. The same functions work with source estimate. The only differences\nare the origin of the data, the size, and the connectivity definition.\nIt can be used for single trials or for groups of subjects.\n\nVisualize clusters\n------------------\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# configure variables for visualization\ntimes = epochs.times * 1e3\ncolors = 'r', 'r', 'steelblue', 'steelblue'\nlinestyles = '-', '--', '-', '--'\n\n# grand average as numpy arrray\ngrand_ave = np.array(X).mean(axis=1)\n\n# get sensor positions via layout\npos = mne.find_layout(epochs.info).pos\n\n# loop over significant clusters\nfor i_clu, clu_idx in enumerate(good_cluster_inds):\n # unpack cluster information, get unique indices\n time_inds, space_inds = np.squeeze(clusters[clu_idx])\n ch_inds = np.unique(space_inds)\n time_inds = np.unique(time_inds)\n\n # get topography for F stat\n f_map = T_obs[time_inds, ...].mean(axis=0)\n\n # get signals at significant sensors\n signals = grand_ave[..., ch_inds].mean(axis=-1)\n sig_times = times[time_inds]\n\n # create spatial mask\n mask = np.zeros((f_map.shape[0], 1), dtype=bool)\n mask[ch_inds, :] = True\n\n # initialize figure\n fig, ax_topo = plt.subplots(1, 1, figsize=(10, 3))\n title = 'Cluster #{0}'.format(i_clu + 1)\n fig.suptitle(title, fontsize=14)\n\n # plot average test statistic and mark significant sensors\n image, _ = plot_topomap(f_map, pos, mask=mask, axes=ax_topo,\n cmap='Reds', vmin=np.min, vmax=np.max)\n\n # advanced matplotlib for showing image with figure and colorbar\n # in one plot\n divider = make_axes_locatable(ax_topo)\n\n # add axes for colorbar\n ax_colorbar = divider.append_axes('right', size='5%', pad=0.05)\n plt.colorbar(image, cax=ax_colorbar)\n ax_topo.set_xlabel('Averaged F-map ({:0.1f} - {:0.1f} ms)'.format(\n *sig_times[[0, -1]]\n ))\n\n # add new axis for time courses and plot time courses\n ax_signals = divider.append_axes('right', size='300%', pad=1.2)\n for signal, name, col, ls in zip(signals, condition_names, colors,\n linestyles):\n ax_signals.plot(times, signal, color=col, linestyle=ls, label=name)\n\n # add information\n ax_signals.axvline(0, color='k', linestyle=':', label='stimulus onset')\n ax_signals.set_xlim([times[0], times[-1]])\n ax_signals.set_xlabel('time [ms]')\n ax_signals.set_ylabel('evoked magnetic fields [fT]')\n\n # plot significant time range\n ymin, ymax = ax_signals.get_ylim()\n ax_signals.fill_betweenx((ymin, ymax), sig_times[0], sig_times[-1],\n color='orange', alpha=0.3)\n ax_signals.legend(loc='lower right')\n ax_signals.set_ylim(ymin, ymax)\n\n # clean up viz\n mne.viz.tight_layout(fig=fig)\n fig.subplots_adjust(bottom=.05)\n plt.show()" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Exercises\n----------\n\n- What is the smallest p-value you can obtain, given the finite number of\n permutations?\n- use an F distribution to compute the threshold by traditional significance\n levels. Hint: take a look at ``scipy.stats.distributions.f``\n\n" ], "cell_type": "markdown", "metadata": {} } ], "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" } } } }