2b_ClassificationSubjects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification to identify subjects\n",
    "\n",
    "This notebook reproduces plots that appear in Fig. 4 of the paper:\n",
    "\n",
    "- Gilson M, Zamora-López G, Pallarés V, Adhikari MH, Senden M, Tauste Campo A, Mantini D, Corbetta M, Deco G, Insabato A (submitted) \"Model-based whole-brain effective connectivity to study distributed cognition in health and disease\", bioRxiv; https://doi.org/10.1101/531830.\n",
    "\n",
    "The goal of to compare the connectivity measures (in particular, effective versus functional connectivity) in identifying subjects.\n",
    "\n",
    "It uses the scikit-learn library, see https://scikit-learn.org/ for details and tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Toggle to True to create directory and store results there\n",
    "save_outputs = False\n",
    "if save_outputs:\n",
    "    import os\n",
    "    res_dir = 'classif/'\n",
    "    if not os.path.exists(res_dir):\n",
    "        os.mkdir(res_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load the data and calculate the features: 'EC' stands for effective connectivity (estimated for the MOU model, see the *MOU_EC_Estimation* notebook), 'FC' for the BOLD correlations, 'FC + mask' for the BOLD correlations only concerning the ROI connections that exist according to the structural data (diffusion tensor imaging), 'PC' for the partial correlations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import dependencies\n",
    "import numpy as np\n",
    "import scipy.stats as stt\n",
    "import sklearn.linear_model as skllm\n",
    "import sklearn.neighbors as sklnn\n",
    "import sklearn.discriminant_analysis as skda\n",
    "import sklearn.preprocessing as skppc\n",
    "import sklearn.pipeline as skppl\n",
    "import sklearn.metrics as skm\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set data parameters\n",
    "param_dir = 'model_param_movie/'\n",
    "n_sub = 22 # number of subjects\n",
    "n_run = 5  # number of sessions per subject (2 rest + 3 movie)\n",
    "N = 66     # number of ROIs\n",
    "\n",
    "# Load data previously generated in notebook 'MOU_EC_Estimation.ipynb'\n",
    "# Effective connectivity estimated using the MOU dynamic model\n",
    "EC = np.load(param_dir + 'J_mod.npy')\n",
    "# Mask of existing connections in SC (structural connectivity)\n",
    "mask_EC = np.load(param_dir + 'mask_EC.npy')\n",
    "# BOLD covariances (without time lag)\n",
    "FC0 = np.load(param_dir + 'FC_emp.npy')[:,:,0,:,:]\n",
    "# triangular mask to retain half of the matrix elemenets in symmetric matrices\n",
    "mask_tri = np.tri(N,N,-1, dtype=bool)\n",
    "\n",
    "# Calculate features\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        # apply z-scoring to EC\n",
    "        EC[i_sub,i_run,mask_EC] = stt.zscore(EC[i_sub,i_run,mask_EC])\n",
    "\n",
    "corr = np.copy(FC0)\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        corr[i_sub,i_run,:,:] /= np.sqrt( np.outer(corr[i_sub,i_run,:,:].diagonal(), \n",
    "                                                   corr[i_sub,i_run,:,:].diagonal()) )\n",
    "\n",
    "PC = np.copy(FC0)\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        PC[i_sub,i_run,:,:] = -np.linalg.pinv(PC[i_sub,i_run,:,:])\n",
    "        PC[i_sub,i_run,:,:] /= np.sqrt( np.outer(PC[i_sub,i_run,:,:].diagonal(),\n",
    "                                                PC[i_sub,i_run,:,:].diagonal()) )\n",
    "\n",
    "# Plotting labels for connectivity measures\n",
    "n_conn = 4\n",
    "label_conn = ['EC', 'FC', 'FC + mask', 'PC']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code repeats the classification by splitting the fMRI sessions into a train set and test set (`n_rep` times). Here 2 sessions are used for training and 3 for testing, the split being randomly performed for each subject.\n",
    "\n",
    "We compare 2 classifiers:\n",
    "\n",
    "- MLR is the multinomial logistic regression;\n",
    "- 1NN is the 1-nearest-neighbor.\n",
    "\n",
    "See the scikit-learn [website](https://scikit-learn.org/stable/) for details.\n",
    "\n",
    "Note that the chance level here is equal to `1 / n_sub`, since all subjects have the same number of sessions in the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Classifier and learning parameters\n",
    "c_MLR = skppl.make_pipeline( skppc.StandardScaler(), \n",
    "                            skllm.LogisticRegression(C=10, penalty='l2', \n",
    "                                                     multi_class='multinomial', \n",
    "                                                     solver='lbfgs', \n",
    "                                                     max_iter=500) )\n",
    "c_1NN = sklnn.KNeighborsClassifier(n_neighbors=1, \n",
    "                                   algorithm='brute', \n",
    "                                   metric='correlation')\n",
    "\n",
    "# Labels of sessions for classification\n",
    "sub_labels = np.repeat(np.arange(n_sub).reshape([-1,1]), n_run, axis=1)\n",
    "\n",
    "# Chance level\n",
    "chance_level = 1. / n_sub\n",
    "\n",
    "# Store results\n",
    "n_rep = 10\n",
    "# last index: MLR/1NN\n",
    "perf = np.zeros([n_rep,n_conn,2])\n",
    "# confusion matrices\n",
    "CM = np.zeros([n_rep,n_conn,n_sub,n_sub], dtype=float)\n",
    " \n",
    "\n",
    "# repeat classification\n",
    "for i_rep in range(n_rep):\n",
    "    # Split indices in train and test sets \n",
    "    # (2 sessions are randomly chosen for each subject)\n",
    "    train_ind = np.zeros([n_sub,n_run], dtype=bool)\n",
    "    for i_sub in range(n_sub):\n",
    "        while train_ind[i_sub,:].sum()<2:\n",
    "            train_ind[:,np.random.randint(n_run)] = True\n",
    "    test_ind = np.logical_not(train_ind)\n",
    "\n",
    "    # loop over connectivity measures\n",
    "    for i_conn in range(n_conn):\n",
    "\n",
    "        if i_conn==0:\n",
    "            # vectorized EC matrices (only retaining existing connections)\n",
    "            vect_features = EC[:,:,mask_EC]\n",
    "        elif i_conn==1:\n",
    "            # vectorized FC matrices (only retaining low triangle)\n",
    "            vect_features = corr[:,:,mask_tri]\n",
    "        elif i_conn==2:\n",
    "            # vectorized FC matrices (only retaining SC existing connections)\n",
    "            vect_features = corr[:,:,mask_EC] \n",
    "        else:\n",
    "            # vectorized PC  matrices (only retaining low triangle)\n",
    "            vect_features = PC[:,:,mask_tri]\n",
    "\n",
    "        # train and test classifiers with subject labels\n",
    "        c_MLR.fit(vect_features[train_ind,:], sub_labels[train_ind])\n",
    "        perf[i_rep,i_conn,0] = c_MLR.score(vect_features[test_ind,:], \n",
    "                                           sub_labels[test_ind])\n",
    "\n",
    "        c_1NN.fit(vect_features[train_ind,:], sub_labels[train_ind])\n",
    "        perf[i_rep,i_conn,1] = c_1NN.score(vect_features[test_ind,:], \n",
    "                                           sub_labels[test_ind])\n",
    "\n",
    "        # confusion matrix for MLR\n",
    "        CM[i_rep,i_conn,:,:] = \\\n",
    "            skm.confusion_matrix( y_true=sub_labels[test_ind], \n",
    "                                  y_pred=c_MLR.predict(vect_features[test_ind,:]) )\n",
    "\n",
    "# Save results\n",
    "if save_outputs:\n",
    "    np.save(res_dir + 'perf_subs.npy', perf)\n",
    "\n",
    "# Print and plot summary of performance\n",
    "for i_conn in range(n_conn):\n",
    "    print(label_conn[i_conn])\n",
    "    print('Mean/std performance MLR:', perf[:,i_conn,0].mean(), ',', perf[:,i_conn,0].std())\n",
    "    print('Mean/std performance 1NN:', perf[:,i_conn,1].mean(), ',', perf[:,i_conn,1].std())\n",
    "\n",
    "plt.figure()\n",
    "plt.violinplot( perf[:,:,0], positions=np.arange(n_conn)-0.2, widths=[0.25]*n_conn )\n",
    "plt.violinplot( perf[:,:,1], positions=np.arange(n_conn)+0.2, widths=[0.25]*n_conn )\n",
    "plt.plot([-1,n_conn],[chance_level]*2,'--k')\n",
    "plt.axis(xmin=-0.6, xmax=n_conn-0.4, ymin=0, ymax=1.01)\n",
    "plt.xticks(range(n_conn), label_conn, rotation=20)\n",
    "plt.ylabel('accuracy', fontsize=14)\n",
    "if save_outputs:\n",
    "    plt.savefig(res_dir + 'perf_sub.png', format='png')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Confusion matrices are useful to understand the errors performed by the train classifier on the test set. Each row represents the classification of a test sample, with the diagonal element corresponding to the correct classification and off-diagonal elements to errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the confusion matrices\n",
    "for i_conn in range(n_conn):\n",
    "    tag = label_conn[i_conn] + '_'\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.imshow(CM[:,i_conn,:,:].mean(0), origin='lower', cmap='Greens', vmin=0, vmax=3)\n",
    "    plt.xticks([0,9,19], [1,10,20])\n",
    "    plt.yticks([0,9,19], [1,10,20])\n",
    "    plt.xlabel('Predicted subject', fontsize=14)\n",
    "    plt.ylabel('True subject', fontsize=14)\n",
    "    plt.colorbar(ticks=range(4))\n",
    "    plt.title(tag[:-1], fontsize=14)\n",
    "    if save_outputs:\n",
    "        plt.savefig(res_dir + 'conf_mat_' + tag + 'sub.png', format='png')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}