examples/basic_usage.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MEEGnobis: basic usage\n",
    "\n",
    "MEEGnobis can be used to perform multivariate analyses on M/EEG datasets over time. It interfaces with [MNE-Python](https://www.martinos.org/mne/stable/index.html) and assumes that all the data have been already preprocessed according to your taste.\n",
    "\n",
    "The main function is `compute_temporal_rdm`, which we import in the following cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from meegnobis.rsa import compute_temporal_rdm, make_pseudotrials\n",
    "from meegnobis.testing import generate_epoch\n",
    "from meegnobis.log import log, logging\n",
    "# increase log level for this notebook\n",
    "log.setLevel(logging.ERROR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please take time to read the docstring of the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computes pairwise metric across time\n",
      "\n",
      "    Arguments\n",
      "    ---------\n",
      "    epoch : instance of mne.Epoch\n",
      "    targets : array (n_trials,)\n",
      "        target (condition) for each trials\n",
      "    metric: str | BaseMetric | sklearn estimator\n",
      "        type of metric to use, one of 'cityblock', 'correlation', 'cosine',\n",
      "        'dice', 'euclidean', 'hamming', 'sqeuclidean' for distances.\n",
      "        Alternatively, any object with attributes fit/score, similar to\n",
      "        sklearn estimators.\n",
      "    cv : instance of sklearn cross-validator\n",
      "        Cross-validator used to split the data into training/test datasets\n",
      "        (default StratifiedShuffleSplit(n_splits=10, test_size=0.5)\n",
      "    cv_normalize_noise : str | None (default None)\n",
      "        Multivariately normalize the trials before applying the metric\n",
      "        (distance or classification). Normalization is performed with\n",
      "        cross-validation, that is the covariance matrix is estimated in the\n",
      "        training set, and then applied to the test set. Normalization is\n",
      "        always performed on single trials, thus before averaging within each\n",
      "        class if `mean_groups` is set to True.\n",
      "        Valid values for `cv_normalize_noise` are 'epoch' | 'baseline' | None:\n",
      "            - 'epoch' computes the covariance matrix on the entire epoch;\n",
      "            - 'baseline' uses only the baseline condition; it requires `epoch`\n",
      "            to have a valid baseline (i.e., to have performed baseline\n",
      "            correction)\n",
      "    mean_groups : bool (default False)\n",
      "        Whether the trials belonging to each target should be averaged prior\n",
      "        to running the metric. This is useful if the metric is a distance\n",
      "        metric. Should be set to False for classification.\n",
      "    time_diag_only : bool (default False)\n",
      "        Whether to run only for the diagonal in time,\n",
      "        e.g. for train_time == test_time\n",
      "    n_jobs : int (default 1)\n",
      "    batch_size : int (default 200)\n",
      "        size of the batches for cross-validation. To be used if the\n",
      "        number of splits are very large in order to reduce the memory load\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    rdm: ndarray (n_pairwise_targets, n_pairwise_times) or\n",
      "                 (n_pairwise_targets, n_times * n_times)\n",
      "        the cross-validated RDM over time; if the metric is symmetric,\n",
      "        then for efficiency purposes only the upper triangular matrix over\n",
      "        time is returned, and the full matrix can be reconstructed with\n",
      "        `np.triu_indices_from`. Otherwise, the entire flattened matrix over\n",
      "        time is returned, and the full matrix can be reconstructed with\n",
      "        `np.reshape((n_times, n_times))`. The order is row first, then columns.\n",
      "    targets_pairs : list-like (n_pairwise_targets,)\n",
      "        the labels for each row of rdm\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "print(compute_temporal_rdm.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to run temporal RDM we need \n",
    "\n",
    "- An `mne.Epoch` object containing the data\n",
    "- A `targets` array containing the target/class information for each trial\n",
    "- A `metric` to be used on the data. We can pass any object that has scikit-learn's `fit/predict` constructs, such as any scikit-learn's classifier. Alternatively, one can pass a string to use a specific distance metric.\n",
    "- A cross-validation object to perform, well, cross-validation. By default we use `StratifiedShuffleSplit` from scikit-learn, with 10 splits and `test_size = 0.5`. This will split the data randomly for 10 times to generate training/testing sets of equal size.\n",
    "\n",
    "Let's see an example with made up data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n"
     ]
    }
   ],
   "source": [
    "epoch = generate_epoch(n_epochs_cond=10,\n",
    "                       n_conditions=5, \n",
    "                       n_channels=10,\n",
    "                       sfreq=200.0, \n",
    "                       n_times=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<EpochsArray  |  n_events : 50 (all good), tmin : 0.0 (s), tmax : 0.495 (s), baseline : None, ~417 kB, data loaded,\n",
      " '0': 10, '1': 10, '2': 10, '3': 10, '4': 10>\n"
     ]
    }
   ],
   "source": [
    "print(epoch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we generated 50 trials with 10 channels/sensors each, with 5 conditions total (thus 10 trials for each condition), and time between 0 and 0.495 s. We're going to use this dataset for different example analyses. First, we need to baseline correct our data (in this case, we do it by considering as baseline the first 100ms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Applying baseline correction (mode: mean)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<EpochsArray  |  n_events : 50 (all good), tmin : 0.0 (s), tmax : 0.495 (s), baseline : (0, 0.1), ~417 kB, data loaded,\n",
       " '0': 10, '1': 10, '2': 10, '3': 10, '4': 10>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epoch.apply_baseline((0, 0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, to make things more interesting, we can just rename the conditions with semanticly meaningful labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<EpochsArray  |  n_events : 50 (all good), tmin : 0.0 (s), tmax : 0.495 (s), baseline : (0, 0.1), ~417 kB, data loaded,\n",
      " 'cat': 10, 'dog': 10, 'face': 10, 'house': 10, 'shoe': 10>\n"
     ]
    }
   ],
   "source": [
    "labels = ['face', 'shoe', 'house', 'cat', 'dog']\n",
    "epoch.event_id = {lbl: idx for idx, lbl in enumerate(labels)}\n",
    "print(epoch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: cross-validated mahalanobis distance (crossnobis)\n",
    "\n",
    "We will run a form of RSA that is noise-normalized using the sensor-wise covariance estimates across the entire epoch using only the trials in the training set. We are going to use the Euclidean distance as a metric, since the Euclidean distance on such normalized data is in fact the Mahalanobis distance. We are going to use only one split in the training/test set for speed, but you might want to use more splits with real data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedShuffleSplit\n",
    "cv = StratifiedShuffleSplit(n_splits=1, test_size=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['face', 'face', 'face', 'face', 'face', 'face', 'face', 'face', 'face', 'face', 'shoe', 'shoe', 'shoe', 'shoe', 'shoe', 'shoe', 'shoe', 'shoe', 'shoe', 'shoe', 'house', 'house', 'house', 'house', 'house', 'house', 'house', 'house', 'house', 'house', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'dog', 'dog', 'dog', 'dog', 'dog', 'dog', 'dog', 'dog', 'dog', 'dog']\n"
     ]
    }
   ],
   "source": [
    "# we are going to take as targets the event labels\n",
    "targets = [labels[e] for e in epoch.events[:, -1]]\n",
    "print(targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we're running a distance metric, we need to set the argument `mean_groups` to `True`: this will average the trials for each class or group in both the training/test set and then compare the average patterns between sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "Estimating covariance using SHRUNK\n",
      "Done.\n",
      "Using cross-validation to select the best estimator.\n",
      "Number of samples used : 2500\n",
      "[done]\n",
      "log-likelihood on unseen data (descending order):\n",
      "   shrunk: -359.709\n",
      "selecting best estimator: shrunk\n",
      "estimated rank (mag): 10\n",
      "Setting small MEG eigenvalues to zero.\n",
      "Not doing PCA for MEG.\n"
     ]
    }
   ],
   "source": [
    "rdm, target_pairs = compute_temporal_rdm(\n",
    "    epoch, targets, cv=cv, metric='euclidean', cv_normalize_noise='epoch', mean_groups=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((15, 5050), 15)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rdm.shape, len(target_pairs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a result we got an array of shape `(15, 5050)`. Each row correspond to a pair of categories (labeled in `target_pairs`), whereas each column correspond to a pair of timepoints on which the RDM is computed on; because we used a symmetric metric, only 5050 (instead of 100 * 100) pairs of timepoints are returned, corresponding to the upper triangular matrix (diagonal included) over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['cat+cat', 'cat+dog', 'cat+face', 'cat+house', 'cat+shoe', 'dog+dog', 'dog+face', 'dog+house', 'dog+shoe', 'face+face', 'face+house', 'face+shoe', 'house+house', 'house+shoe', 'shoe+shoe']\n"
     ]
    }
   ],
   "source": [
    "print(target_pairs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: decoding\n",
    "\n",
    "We can use the same code to run decoding using any of scikit-learn's classifiers. We're going to use an SVM with linear kernel, without noise-normalization, and increasing the number of splits to 4 this time. Because we're using a classifier, we do not want to average the trials within each class, so we set `mean_groups` to `False`. The function will perform all pairwise binary classifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "clf = SVC(kernel='linear')\n",
    "cv = StratifiedShuffleSplit(n_splits=4, test_size=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "50 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "50 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "50 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n"
     ]
    }
   ],
   "source": [
    "rdm, target_pairs = compute_temporal_rdm(epoch, targets, cv=cv, metric=clf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((15, 10000), 15)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rdm.shape, len(target_pairs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the result is a `(15, 10000)` array. Because the classifier is not symmetric over time (e.g., training at time `t0` and testing at time `t1` is different than training at time `t1` and testing at time `t0`), the entire matrix is returned. The order is rows first, then columns, so that one can create a tensor `n_pairwise_targets X n_time X n_time` by simply using reshape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15, 100, 100)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rdm = rdm.reshape((-1, 100, 100))\n",
    "rdm.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: decoding with \"pseudotrials\"\n",
    "\n",
    "SNR can be increased by using \"pseudotrials\", generated by randomly averaging subsets of trials. An approach that is used in the literature (see e.g., Cichy et al., 2014) is to repeat this averaging 100 times, randomly selecting trials to average every time, and testing the classifiers using a leave-one-sample-out scheme. We can easily do this in `meegnobis` using `make_pseudotrials`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Create pseudotrials by averaging within each group defined in `targets`\n",
      "    a number `navg` of trials. The trials are randomly divided into groups of\n",
      "    size `navg` for each target. If the number of trials in a group is not\n",
      "    divisible by `navg`, one pseudotrials will be created by averaging the\n",
      "    remainder trials.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    epoch : mne.Epoch\n",
      "    targets : (n_epochs,)\n",
      "    navg : int\n",
      "        number of trials to average\n",
      "    rng : random number generator\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    avg_epoch : mne.Epoch\n",
      "        epoch with ceil(n_epochs/navg) trials\n",
      "    avg_targets : (ceil(n_epochs/navg),)\n",
      "        unique targets corresponding to the avg_trials\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "print(make_pseudotrials.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example dataset we have 10 trials for each condition, so we can average three trials at a time to generate four samples for each condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n"
     ]
    }
   ],
   "source": [
    "avg_epoch, avg_targets = make_pseudotrials(epoch, targets, navg=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<EpochsArray  |  n_events : 20 (all good), tmin : 0.495 (s), tmax : 0.99 (s), baseline : (0, 0.1), ~182 kB, data loaded,\n",
      " 'cat': 4, 'dog': 4, 'face': 4, 'house': 4, 'shoe': 4>\n"
     ]
    }
   ],
   "source": [
    "print(avg_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['cat', 'cat', 'cat', 'cat', 'dog', 'dog', 'dog', 'dog', 'face', 'face', 'face', 'face', 'house', 'house', 'house', 'house', 'shoe', 'shoe', 'shoe', 'shoe']\n"
     ]
    }
   ],
   "source": [
    "print(avg_targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can  repeat this 5 times, and perform classification every time using a leave-one-out-sample cross-validation scheme, and then averaging across the 5 folds. We're still going to use `StratifiedShuffleSplit`, but this time we ask a test set with size exactly equal to the number of unique targets. In this way we'll have a single sample for every category in the test set.\n",
    "\n",
    "We can write the code easily as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n",
      "20 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "0 bad epochs dropped\n"
     ]
    }
   ],
   "source": [
    "n_unique_targets = len(np.unique(targets))\n",
    "clf = SVC(kernel='linear')\n",
    "\n",
    "out = []\n",
    "for _ in range(5):\n",
    "    # we reset the cv inside the loop so that it gets random splits every time\n",
    "    cv = StratifiedShuffleSplit(n_splits=1, test_size=n_unique_targets)    \n",
    "    # get randomly generate pseudotrials\n",
    "    avg_epoch, avg_targets = make_pseudotrials(epoch, targets, navg=3)\n",
    "    # compute RDM over time\n",
    "    rdm, target_pairs = compute_temporal_rdm(avg_epoch, avg_targets, cv=cv, metric=clf)\n",
    "    out.append(rdm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, (15, 10000))"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(out), out[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can then finally average to obtain the final RDM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15, 10000)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "final_rdm = np.mean(out, axis=0)\n",
    "final_rdm.shape"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}