intro_ml.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import random\n",
    "import pandas as pd\n",
    "import itertools\n",
    "import importlib\n",
    "import mlintro_min as mli"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### To anticipate section 2 on nilearn: please download the data below if not done already"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Note:** please run these commands from a different notebook / terminal to avoid having your Jupyter kernel busy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn import datasets\n",
    "# Replace \"datasets/age_netmats\" with whatever path you would like\n",
    "data = datasets.fetch_development_fmri(data_dir='datasets/age_netmats')\n",
    "atlases = datasets.fetch_atlas_basc_multiscale_2015(version='sym', data_dir='datasets/age_netmats')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Simple regression (OLS) model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Range of cannon angle psi to investigate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "psi_min=20\n",
    "psi_max=60"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Read pre-generated source of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ref_df = pd.read_csv('datasets/offline/refdata_100K.csv')\n",
    "# Get a very small subset of inter-spread data for plotting\n",
    "ref_df_light = mli.get_ref_light(ref_df, psi_min=psi_min, psi_max=psi_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Generate a single dataset of 50 observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "ds = mli.get_datasets(ref_df, n_datasets=1, sample_size=50, psi_min=psi_min, psi_max=psi_max)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Explore dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We are only interested in experimental measures of angle and range, `exp_angle` and `exp_range` respectively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Plot dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4, 4))\n",
    "sns.scatterplot(\"exp_angle\", \"exp_range\", data=ds, color='blue')\n",
    "plt.ylim([50, 80]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 1: model selection with train / test split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We will compare three models of different complexity (different number of parameters):\n",
    "   * Standard linear model (2 parameters: $\\beta_0$, $\\beta_\\psi$)\n",
    "   * Linear model with the feature squared (3 parameters: $\\beta_0$, $\\beta_\\psi$, $\\beta_{\\psi^2}$)\n",
    "   * Linear model with up to $5^{th}$ power of feature (6 parameters: $\\beta_0$, $\\beta_\\psi$, $\\beta_{\\psi^2}$, $\\beta_{\\psi^3}$, $\\beta_{\\psi^4}$, $\\beta_{\\psi^5}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's divide our dataset into training and test set !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X = ds[['exp_angle']]\n",
    "y = ds['exp_range']\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=42)\n",
    "print(\"Number of observations: {} for training and {} for testing\".format(len(X_train), len(X_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-14T20:06:42.259063Z",
     "start_time": "2020-06-14T20:06:42.249379Z"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X_train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Linear (OLS) model with single feature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "##### Let's train the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error, r2_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Create new estimator object (each model is an object in `sklearn`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "lm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('beta_psi: {}, beta_0: {}'.format(lm.coef_, lm.intercept_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Predict, and compute performance scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm.predict([[30], [45], [60]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The training score can provide an idea of best possible performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# We score the predictions made on the same data we trained on\n",
    "y_pred = lm.predict(X_train)\n",
    "\n",
    "R2_train = r2_score(y_train, y_pred)\n",
    "MSE_train = mean_squared_error(y_train, y_pred)\n",
    "print('lm training performance is R2: {:0.2f} and MSE: {:0.2f}'.format(R2_train, MSE_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "`R2` can be obtained directly from the model / estimator object as it is the default for `LinearRegression`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm.score(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's save all our models scores to compare models later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results = []\n",
    "train_test_results.append({'model': 'lm', 'stage': 'train', 'scorer': 'r2', 'val': R2_train})\n",
    "train_test_results.append({'model': 'lm', 'stage': 'train', 'scorer': 'MSE', 'val': -MSE_train})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Let's test on unseen data and get testing score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "See possible scores: https://scikit-learn.org/stable/modules/model_evaluation.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The testing score can provide an idea of performance for unseen data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# We score the predictions made on unseen data\n",
    "y_pred = lm.predict(X_test)\n",
    "\n",
    "R2_test = r2_score(y_test, y_pred)\n",
    "MSE_test = mean_squared_error(y_test, y_pred)\n",
    "print('lm testing performance is R2: {:0.2f} and MSE: {:0.2f}'.format(R2_test, MSE_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Again let's save our scores for comparison later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results.append({'model': 'lm', 'stage': 'test', 'scorer': 'r2', 'val': R2_test})\n",
    "train_test_results.append({'model': 'lm', 'stage': 'test', 'scorer': 'MSE', 'val': -MSE_test})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Polynomial degree 2 model (polynomial use feature powers as additional features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "##### Train and get training score of degree 2 polynomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The polynomial terms are created and then added to a standard `LinearRegression`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "poly_transformer = PolynomialFeatures(degree=2)\n",
    "lm_deg2 = LinearRegression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A polynomial model is fitted by first transforming the features and then training a `LinearRegression` on the transformed features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X_train_deg2 = poly_transformer.fit_transform(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "split",
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "split",
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X_train_deg2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "lm_deg2.fit(X_train_deg2, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm_deg2.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's compute training score  \n",
    "**Warning**: `LinearRegression` needs to be applied to the *transformed* (polynomial) features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "y_pred = lm_deg2.predict(X_train_deg2)\n",
    "\n",
    "R2_train = r2_score(y_train, y_pred)\n",
    "MSE_train = mean_squared_error(y_train, y_pred)\n",
    "print('poly deg2 training performance is R2: {:0.2f} and MSE: {:0.2f}'.format(R2_train, MSE_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results.append({'model': 'lm_deg2', 'stage': 'train', 'scorer': 'r2', 'val': R2_train})\n",
    "train_test_results.append({'model': 'lm_deg2', 'stage': 'train', 'scorer': 'MSE', 'val': -MSE_train})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Let's test on unseen data and get testing score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Again, should not forget to transform the test data to get polynomial features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X_test_deg2 = poly_transformer.fit_transform(X_test) \n",
    "\n",
    "y_pred = lm_deg2.predict(X_test_deg2)\n",
    "\n",
    "R2_test = r2_score(y_test, y_pred)\n",
    "MSE_test = mean_squared_error(y_test, y_pred)\n",
    "print('poly deg2 testing performance is R2: {:0.2f} and MSE: {:0.2f}'.format(R2_test, MSE_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results.append({'model': 'lm_deg2', 'stage': 'test', 'scorer': 'r2', 'val': R2_test})\n",
    "train_test_results.append({'model': 'lm_deg2', 'stage': 'test', 'scorer': 'MSE', 'val': -MSE_test})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Polynomial degree 5 model, introducing pipeline object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "It is common to use preprocessing steps such as features transformation. To avoid repeating these processing steps (e.g. when testing the model on new data) the `Pipeline` object is very useful. It builds a workflow which can be called with a single command. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Train and get training score of degree 5 polynomial with *pipeline*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The creation of the polynomial terms are embed in the `Pipeline` object which represents our new estimator (i.e. the one to use with `fit`, `predict`, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm_deg5 = Pipeline([('poly_transformer', PolynomialFeatures(degree=5)),\n",
    "                    ('lm', LinearRegression())])\n",
    "lm_deg5.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm_deg5['lm'].coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's compute training score\n",
    "y_pred = lm_deg5.predict(X_train)\n",
    "\n",
    "R2_train = r2_score(y_train, y_pred)\n",
    "MSE_train = mean_squared_error(y_train, y_pred)\n",
    "print('poly deg5 training performance is R2: {:0.2f} and MSE: {:0.2f}'.format(R2_train, MSE_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results.append({'model': 'lm_deg5', 'stage': 'train', 'scorer': 'r2', 'val': R2_train})\n",
    "train_test_results.append({'model': 'lm_deg5', 'stage': 'train', 'scorer': 'MSE', 'val': -MSE_train})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "##### Get score on test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "No more transformation to do manually as it is embed in the pipeline!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "y_pred = lm_deg5.predict(X_test)\n",
    "R2_test = r2_score(y_test, y_pred)\n",
    "MSE_test = mean_squared_error(y_test, y_pred)\n",
    "print('poly deg5 testing performance is R2: {:0.2f} and MSE: {:0.2f}'.format(R2_test, MSE_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results.append({'model': 'lm_deg5', 'stage': 'test', 'scorer': 'r2', 'val': R2_test})\n",
    "train_test_results.append({'model': 'lm_deg5', 'stage': 'test', 'scorer': 'MSE', 'val': -MSE_test})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "train_test_results_df = pd.DataFrame(train_test_results)\n",
    "r2_results = train_test_results_df.loc[train_test_results_df['scorer'] == 'r2']\n",
    "MSE_results = train_test_results_df.loc[train_test_results_df['scorer'] == 'MSE']\n",
    "with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "    g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=train_test_results_df, \n",
    "                    kind=\"bar\", sharey=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For many other regression models implemented in sklearn, cf https://scikit-learn.org/stable/supervised_learning.html#supervised-learning "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 2: Model selection with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "https://scikit-learn.org/stable/modules/cross_validation.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "KFold?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's see what K fold is actually doing on a 10-observation dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "center",
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X_example = pd.DataFrame({'obs': \n",
    "                          np.random.randint(0, 1000, 10)})\n",
    "X_example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "kf = KFold(n_splits=5)\n",
    "for (ix_train, ix_test) in kf.split(X_example):\n",
    "    print(\"train ix:\", ix_train, \"testix:\", ix_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "kfs = KFold(n_splits=5, shuffle=True, random_state=42)\n",
    "for (ix_train, ix_test) in kfs.split(X_example):\n",
    "    print(\"train ix:\", ix_train, \"test ix:\", ix_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "kf = KFold(n_splits=5, shuffle=True, random_state=42)\n",
    "print(\"Size X: {}, with X train: {} and X test: {}\".format(len(X), len(X_train), len(X_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ml_models = {'lm': LinearRegression(),\n",
    "             'lm_deg2': Pipeline([('poly_transformer', PolynomialFeatures(degree=2)),\n",
    "                                  ('lm', LinearRegression())]),\n",
    "             'lm_deg5': Pipeline([('poly_transformer', PolynomialFeatures(degree=5)),\n",
    "                                  ('lm', LinearRegression())])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "#from sklearn.model_selection import StratifiedKFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "kf_results = []\n",
    "\n",
    "kfs = KFold(n_splits=10, shuffle=True, random_state=42)\n",
    "for i_f, (ix_train, ix_test) in enumerate(kfs.split(X_train)):\n",
    "    # Loop over models\n",
    "    for mod_name in ml_models.keys():\n",
    "        # Fit the model on the training folds\n",
    "        ml_models[mod_name].fit(X_train.iloc[ix_train], y_train.iloc[ix_train])\n",
    "        # Test on both the training and testing folds to check for over-/under-fitting\n",
    "        y_pred_train = ml_models[mod_name].predict(X_train.iloc[ix_train])\n",
    "        y_pred_test = ml_models[mod_name].predict(X_train.iloc[ix_test])\n",
    "        # R2\n",
    "        kf_results.append({'model': mod_name, 'fold': i_f, 'stage': 'train', 'scorer': 'r2', \n",
    "                           'val': r2_score(y_train.iloc[ix_train], y_pred_train)})\n",
    "        kf_results.append({'model': mod_name, 'fold': i_f, 'stage': 'test', 'scorer': 'r2', \n",
    "                           'val': r2_score(y_train.iloc[ix_test], y_pred_test)})\n",
    "        # MSE\n",
    "        kf_results.append({'model': mod_name, 'fold': i_f, 'stage': 'train', 'scorer': 'MSE', \n",
    "                           'val': -mean_squared_error(y_train.iloc[ix_train], y_pred_train)})\n",
    "        kf_results.append({'model': mod_name, 'fold': i_f, 'stage': 'test', 'scorer': 'MSE', \n",
    "                           'val': -mean_squared_error(y_train.iloc[ix_test], y_pred_test)})\n",
    "kf_results_df = pd.DataFrame(kf_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "for mod_name in ml_models.keys():\n",
    "    kf_df = kf_results_df.loc[kf_results_df['model'] == mod_name]\n",
    "    with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "        g = sns.catplot(x=\"fold\", y=\"val\", hue=\"stage\", col=\"scorer\", data=kf_df, kind=\"bar\", \n",
    "                        sharey=False, height=3, aspect=1.5)\n",
    "        g.axes[0,0].set_ylim(-0.2,1)\n",
    "        #g.axes[0,1].set_ylim(0, 45)\n",
    "        g.fig.suptitle(mod_name, y=1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "    g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=kf_results_df, kind=\"box\",\n",
    "                    sharey=False, height=3, aspect=1.5)\n",
    "    g.axes[0, 0].set_ylim(-0.2, 1)\n",
    "    g.axes[0, 1].set_ylim(-50, 0)\n",
    "    g.fig.suptitle(\"Score accross fold across models\", y=1.05)\n",
    "    g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=kf_results_df, kind=\"box\", \n",
    "                    sharey=False, height=3, aspect=1.5)\n",
    "    g.axes[0, 0].set_ylim(0.45, 1)\n",
    "    g.axes[0, 1].set_ylim(-18, 0)\n",
    "    g.fig.suptitle(\"Close up on polynomial models\", y=1.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### `sklearn` helps to automate common operations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "cross_val_score?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import fbeta_score, make_scorer\n",
    "mse_scorer = make_scorer(mean_squared_error, greater_is_better=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ml_models = {'lm': LinearRegression(),\n",
    "             'lm_deg2': Pipeline([('poly_transformer', PolynomialFeatures(degree=2)),\n",
    "                                  ('lm', LinearRegression())]),\n",
    "             'lm_deg5': Pipeline([('poly_transformer', PolynomialFeatures(degree=5)),\n",
    "                                  ('lm', LinearRegression())])}\n",
    "# Get cv train AND test scores\n",
    "cv_test_scores = {}\n",
    "for mod_name in ml_models.keys():\n",
    "    cv_test_scores[mod_name] = cross_val_score(ml_models[mod_name], X_train, y_train, cv=kfs,\n",
    "                                               scoring=mse_scorer, n_jobs=-1)\n",
    "cv_test_scores_df = pd.DataFrame(cv_test_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "cv_test_scores_df.boxplot()\n",
    "plt.title('Negative MSE for the three models tested (larger is better)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Want even more with even less code?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_validate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "`cross_validate` returns train *and* test scores on *several* metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "cross_validate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ml_models = {'lm': LinearRegression(),\n",
    "             'lm_deg2': Pipeline([('poly_transformer', PolynomialFeatures(degree=2)),\n",
    "                                  ('lm', LinearRegression())]),\n",
    "             'lm_deg5': Pipeline([('poly_transformer', PolynomialFeatures(degree=5)),\n",
    "                                  ('lm', LinearRegression())])}\n",
    "# Get cv train AND test scores\n",
    "cv_scores = {}\n",
    "for mod_name in ml_models.keys():\n",
    "    cv_scores[mod_name] = cross_validate(ml_models[mod_name], X_train, y_train, cv=kfs, \n",
    "                                         scoring=['r2', 'neg_mean_squared_error'], n_jobs=-1,\n",
    "                                         return_train_score=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def crossval_to_df(cv_dict):\n",
    "    crossval_results = []\n",
    "    for model in cv_dict.keys():\n",
    "        for scorer in cv_dict[model].keys():\n",
    "            if scorer.startswith('train_'):\n",
    "                score = scorer.replace('train_', '')\n",
    "                for i_val, val in enumerate(cv_dict[model][scorer]):\n",
    "                    crossval_results.append({'model': model, 'fold': i_val, 'stage': 'train', \n",
    "                                             'scorer': score, 'val': val})\n",
    "            elif scorer.startswith('test_'):\n",
    "                score = scorer.replace('test_', '')\n",
    "                for i_val, val in enumerate(cv_dict[model][scorer]):\n",
    "                    crossval_results.append({'model': model, 'fold': i_val, 'stage': 'test', \n",
    "                                             'scorer': score, 'val': val})\n",
    "    return pd.DataFrame(crossval_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     1
    ],
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "crossval_df = crossval_to_df(cv_scores)\n",
    "for mod_name in crossval_df['model'].unique():\n",
    "    kf_df = crossval_df.loc[crossval_df['model'] == mod_name]\n",
    "    with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "        g = sns.catplot(x=\"fold\", y=\"val\", hue=\"stage\", col=\"scorer\", data=crossval_df,\n",
    "                        kind=\"bar\", sharey=False, height=3, aspect=1.5, ci=None)\n",
    "        g.axes[0,0].set_ylim(-0.2,1)\n",
    "        g.axes[0,1].set_ylim(-45, 0)\n",
    "        g.fig.suptitle(mod_name, y=1.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Penalized regression and feature selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We will use our new skills to compare models for predicting age from brain functional connectivity matrices. We will use data ready to be directly analyzed with machine learning. We will see later how to generate these data with `nilearn`. We will compare in this example:\n",
    "* Feature selection + classic linear regression\n",
    "* Penalized regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Feature selection** is useful for:\n",
    "* Helping prevent overfitting\n",
    "* Obtaining parsimonious models\n",
    "* Speeding up calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Note:** This section is using the same MRI data example as in the notebook Jake Vogel presented at the 2020 Brainhack School: https://github.com/neurodatascience/course-materials-2020/blob/master/lectures/14-may/03-intro-to-machine-learning/ML_Regression_Tutorial.ipynb  \n",
    "\n",
    "The nilearn part here is very similar while the machine learning part is quite different (different subset of data, different models). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Get the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "datasets.fetch_development_fmri?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Better starting the download now in a different Jupyter notebook / terminal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "data = datasets.fetch_development_fmri(data_dir='datasets/age_netmats')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 1: individual model evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### ML data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Create the feature matrix `X` and labels `y`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "features_file = 'datasets/offline/MAIN_BASC064_subsamp_features.npz'\n",
    "X_all = np.load(features_file)['a']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "no_download = True\n",
    "if no_download:\n",
    "    pheno_file = 'datasets/offline/pheno.csv'\n",
    "    pheno_df=pd.read_csv(pheno_file)\n",
    "else:\n",
    "    pheno_df = pd.DataFrame(data.phenotypic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's have a look at `X_all` (we should always look at our data when we can)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,16))\n",
    "plt.imshow(X_all, aspect='auto')\n",
    "plt.colorbar()\n",
    "plt.title('feature matrix')\n",
    "plt.xlabel('features')\n",
    "plt.ylabel('subjects')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,16))\n",
    "g = sns.heatmap(X_all, vmin=-1, vmax=1, cmap='viridis')\n",
    "g.set(xlabel='Features', ylabel='Subjects', title='Feature matrix');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pheno_df.iloc[36:63]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,4))\n",
    "plt.plot(np.std(X_all, axis=1), label = 'std')\n",
    "median_s_std = np.median(np.std(X_all, axis=1))\n",
    "std_s_std = np.std(np.std(X_all, axis=1))\n",
    "min_s_std = np.min(np.std(X_all, axis=1))\n",
    "max_s_std = np.max(np.std(X_all, axis=1))\n",
    "plt.axhline(median_s_std, label='median std', color='g')\n",
    "plt.vlines(np.where(pheno_df[\"Age\"]<5)[0], ymin=min_s_std, ymax=max_s_std, color='r', alpha=0.4,\n",
    "           label='child less than 5 years old')\n",
    "plt.title('Random QA plot')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pheno_df['AgeGroup'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Plotting hist without kde\n",
    "ax = sns.distplot(pheno_df[['Age']], kde=False)\n",
    "ax2 = ax.twinx()\n",
    "sns.distplot(pheno_df[['Age']], ax=ax2, kde=True, hist=False)\n",
    "ax2.set_yticks([]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Number of subjects less than 13: {}\".format(pheno_df['Age'].lt(13).sum()))\n",
    "print(\"Number of subjects more than 13: {}\".format(pheno_df['Age'].ge(13).sum()))\n",
    "print(\"Total number of subjects: {}\".format(len(pheno_df)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Get data in shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's look at models predicting age only for children until 12 (could do same exercise with also removing children less than 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lt13_boolmask = (pheno_df['Age'] < 15)\n",
    "y = pheno_df[['Age']].loc[lt13_boolmask]\n",
    "X = X_all[lt13_boolmask, :]\n",
    "age_class = pheno_df[['AgeGroup']].loc[lt13_boolmask]\n",
    "print(\"Number of subjects is {} and number of features is {}\".format(*X.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Some serious overfitting ahead if we don't do anything!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's split our dataset to have round numbers (easier to conceptualize CV we just studied).Options to `train_test_split` are:\n",
    "* 100/22 train/test subjects \n",
    "* shuffle, very important as data has some order!\n",
    "* stratification to make sure under-represented 8-12 group balanced in all folds\n",
    "* random see to reproduce the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "train_test_split?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tts_ix = np.arange(len(X))\n",
    "X_train, X_test, y_train, y_test, age_class_train, age_class_test, X_ix, y_ix = train_test_split(\n",
    "    X, y, age_class, tts_ix, test_size=22, shuffle = True, stratify=age_class, random_state = 42)\n",
    "print(\"Number of subjects: training={}, testing={}\".format(len(y_train), len(y_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's make sure the stratified splitting worked out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "sns.distplot(y_train, label='train')\n",
    "sns.distplot(y_test, label='test')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Classic lineal (OLS) model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Just for the sake of confirming what we know: what should we get for training and testing performance if we fit a standard linear regression model with 100 datapoints and 2000+ features? Let's find out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lm = LinearRegression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Evaluate model with (stratified!) cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Option 1**: manually (as the stratification can only be done on label data `y`, not on `age_class_train`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "skf_results = []\n",
    "\n",
    "skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)\n",
    "for i_f, (ix_train, ix_test) in enumerate(skf.split(X_train, age_class_train)):\n",
    "        lm.fit(X_train[ix_train], y_train.iloc[ix_train])\n",
    "        # Test on both the testing fold and training folds to check for over-/under-fitting\n",
    "        y_pred_train = lm.predict(X_train[ix_train])\n",
    "        y_pred_test = lm.predict(X_train[ix_test])\n",
    "        # R2\n",
    "        skf_results.append({'model': 'lm', 'fold': i_f, 'stage': 'train', 'scorer': 'r2', \n",
    "                            'val': r2_score(y_train.iloc[ix_train], y_pred_train)})\n",
    "        skf_results.append({'model': 'lm', 'fold': i_f, 'stage': 'test', 'scorer': 'r2', \n",
    "                           'val': r2_score(y_train.iloc[ix_test], y_pred_test)})\n",
    "        # MSE\n",
    "        skf_results.append({'model': 'lm', 'fold': i_f, 'stage': 'train', 'scorer': 'MSE', \n",
    "                           'val': -mean_squared_error(y_train.iloc[ix_train], y_pred_train)})\n",
    "        skf_results.append({'model': 'lm', 'fold': i_f, 'stage': 'test', 'scorer': 'MSE', \n",
    "                           'val': -mean_squared_error(y_train.iloc[ix_test], y_pred_test)})\n",
    "skf_results_df = pd.DataFrame(skf_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Option 2**: create custom cv based on `age_class_train`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "class AgeStratifiedKFold(StratifiedKFold):\n",
    "    def __init__(self, n_splits=5, *, shuffle=False, random_state=None, age_grp=age_class_train):\n",
    "        super().__init__(n_splits=n_splits, shuffle=shuffle, random_state=random_state)\n",
    "        self.age_grp = age_grp\n",
    "    \n",
    "    def split(self, X, y, groups=None):\n",
    "        local_age_grp = self.age_grp.loc[y.index.to_list()]\n",
    "        return super().split(X, local_age_grp, groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "age_skf = AgeStratifiedKFold(n_splits=10, shuffle=True, random_state=42, age_grp=age_class_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "cv_scores = {}\n",
    "cv_scores['lm'] = cross_validate(lm, X_train, y_train, cv=age_skf, \n",
    "                                 scoring=['r2', 'neg_mean_squared_error'], n_jobs=-1,\n",
    "                                 return_train_score=True)\n",
    "cv_scores_df = pd.DataFrame(cv_scores)\n",
    "crossval_df = crossval_to_df(cv_scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### PLot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "skf_df = skf_results_df.loc[skf_results_df['model'] == 'lm']\n",
    "with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "    # Barploat\n",
    "    g = sns.catplot(x=\"fold\", y=\"val\", hue=\"stage\", col=\"scorer\", data=skf_df, kind=\"bar\", \n",
    "                    sharey=False, height=3, aspect=1.5)\n",
    "    g.axes[0,0].set_ylim(-0.2,1)\n",
    "    g.fig.suptitle('lm', y=1.05)\n",
    "    # Boxplot\n",
    "    g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=skf_results_df, kind=\"box\",\n",
    "                    sharey=False, height=3, aspect=1.5)\n",
    "    g.axes[0, 0].set_ylim(-0.2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "(Below are identical plots from the custom cv solution )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "for mod_name in crossval_df['model'].unique():\n",
    "    kf_df = crossval_df.loc[crossval_df['model'] == mod_name]\n",
    "    with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "        # Bar plot\n",
    "        g = sns.catplot(x=\"fold\", y=\"val\", hue=\"stage\", col=\"scorer\", data=crossval_df,\n",
    "                        kind=\"bar\", sharey=False, height=3, aspect=1.5, ci=None,\n",
    "                        hue_order =[\"train\", \"test\"])\n",
    "        g.axes[0,0].set_ylim(-0.2,1)\n",
    "        g.fig.suptitle(mod_name, y=1.05)\n",
    "        # Box plot\n",
    "        g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=crossval_df, \n",
    "                        kind=\"box\", sharey=False, height=3, aspect=1.5, hue_order =[\"train\", \"test\"])\n",
    "        g.axes[0, 0].set_ylim(-0.2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Adding feature selection to OLS model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "https://scikit-learn.org/stable/modules/feature_selection.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectKBest, f_regression, mutual_info_regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's look into two different types of feature selection and create the associated pipelines "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# F-score (correlation with target turned into F score then pvalue)\n",
    "fs_freg = SelectKBest(f_regression, k=10)\n",
    "lm_fs_freg = Pipeline([('fs_freg', fs_freg), \n",
    "                       ('lm', LinearRegression())])\n",
    "# Mutual information score\n",
    "fs_mireg = SelectKBest(mutual_info_regression, k=10)\n",
    "lm_fs_mireg = Pipeline([('fs_mireg', fs_mireg), \n",
    "                        ('lm', LinearRegression())])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "SelectKBest?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Let's do a simple application of one of our pipeline with the default hyper-parameter K for K-best (k=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "age_skf = AgeStratifiedKFold(n_splits=10, shuffle=True, random_state=42, age_grp=age_class_train)\n",
    "cv_scores = {}\n",
    "for ml_model in [lm_fs_freg, lm_fs_mireg]:\n",
    "    cv_scores[ml_model] = cross_validate(ml_model, X_train, y_train, cv=age_skf,\n",
    "                                         scoring=['r2', 'neg_mean_squared_error'],\n",
    "                                         return_train_score=True, n_jobs=-1)\n",
    "cv_scores_df = pd.DataFrame(cv_scores)\n",
    "crossval_df = crossval_to_df(cv_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     1
    ],
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# For visualization, rename model names appropriately\n",
    "def rename_model(model_name):\n",
    "    if str(model_name).startswith('Pipeline'):\n",
    "        return '___'.join([s[0] for s in model_name.steps])\n",
    "    else:\n",
    "        return model_name\n",
    "crossval_df['model'] = crossval_df['model'].apply(rename_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "###### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "for mod_name in crossval_df['model'].unique():\n",
    "    kf_df = crossval_df.loc[crossval_df['model'] == mod_name]\n",
    "    with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "        # Bar plot\n",
    "        g = sns.catplot(x=\"fold\", y=\"val\", hue=\"stage\", col=\"scorer\", data=kf_df,\n",
    "                        kind=\"bar\", sharey=False, height=3, aspect=1.5, ci=None,\n",
    "                        hue_order =[\"train\", \"test\"])\n",
    "        g.axes[0,0].set_ylim(-0.2,1)\n",
    "        g.fig.suptitle(mod_name, y=1.05)\n",
    "        # Box plot\n",
    "        g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=kf_df, \n",
    "                        kind=\"box\", sharey=False, height=3, aspect=1.5,\n",
    "                        hue_order =[\"train\", \"test\"])\n",
    "        g.axes[0,0].set_ylim(-0.2,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Now looking at range of hyper-parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import validation_curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "validation_curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Warning:** the cell below can take some time to run. Could be useful to reduce the range of the `k` hyperparameter (e.g. [5, 15, 25, 50, 75, 100]) if you don't have a lot of CPU resources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "k_range = [1, 5, 10, 15, 20, 25, 50, 75, 100, 150, 200, 250, 500]\n",
    "# Feature selection with f-score\n",
    "f_train_r2, f_test_r2 = validation_curve(lm_fs_freg, X_train, y_train, param_name=\"fs_freg__k\", \n",
    "                                         param_range=k_range, cv=age_skf, scoring='r2', n_jobs=-1)\n",
    "f_train_mse, f_test_mse = validation_curve(lm_fs_freg, X_train, y_train, param_name=\"fs_freg__k\", \n",
    "                                           param_range=k_range, cv=age_skf, \n",
    "                                           scoring='neg_mean_squared_error', n_jobs=-1)\n",
    "# Feature selection with MI-score\n",
    "mi_train_re, mi_test_r2 = validation_curve(lm_fs_mireg, X_train, y_train, param_name=\"fs_mireg__k\", \n",
    "                                           param_range=k_range, cv=age_skf, scoring='r2', n_jobs=-1)\n",
    "mi_train_mse, mi_test_mse = validation_curve(lm_fs_mireg, X_train, y_train, param_name=\"fs_mireg__k\", \n",
    "                                             param_range=k_range, cv=age_skf, \n",
    "                                             scoring='neg_mean_squared_error', n_jobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def valcurve_to_df(train_scores, test_scores, param_range, model, scorer):\n",
    "    valcurve_results = []\n",
    "    n_folds = train_scores.shape[1]\n",
    "    for i_param, param in enumerate(param_range):\n",
    "        valcurve_results.append(pd.DataFrame(\n",
    "            {'score': train_scores[i_param], 'stage': ['train']*n_folds,\n",
    "             'fold': np.arange(n_folds), 'param': np.array([param]*n_folds),\n",
    "             'model': [model]*n_folds, 'scorer': [scorer]*n_folds}))\n",
    "        valcurve_results.append(pd.DataFrame(\n",
    "            {'score': test_scores[i_param], 'stage': ['test']*n_folds,\n",
    "             'fold': np.arange(n_folds), 'param': np.array([param]*n_folds),\n",
    "             'model': [model]*n_folds, 'scorer': [scorer]*n_folds}))\n",
    "    return pd.concat(valcurve_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "valcurve_df = pd.concat([valcurve_to_df(f_train_r2, f_test_r2, k_range, 'lm_fs_freg', 'r2'),\n",
    "                         valcurve_to_df(f_train_mse, f_test_mse, k_range, 'lm_fs_freg', 'mse'),\n",
    "                         valcurve_to_df(mi_train_re, mi_test_r2, k_range, 'lm_fs_mireg', 'r2'),\n",
    "                         valcurve_to_df(mi_train_mse, mi_test_mse, k_range, 'lm_fs_mireg', 'mse')])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Hyper-parameter evaluation results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "for mod_name in valcurve_df['model'].unique():\n",
    "    kf_df = valcurve_df.loc[valcurve_df['model'] == mod_name]\n",
    "    with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "        g = sns.catplot(x='param', y='score', hue='stage', col='scorer', data=kf_df, kind='point',\n",
    "                        sharey=False, hue_order =[\"train\", \"test\"], height=3, aspect=1.5)\n",
    "        #plt.xticks(range(10))\n",
    "        g.set_xticklabels(k_range, rotation=90)\n",
    "        g.axes[0,0].set_ylim(-0.2, 1.1)\n",
    "        g.axes[0,1].set_ylim(-45, 1)\n",
    "        g.fig.suptitle(mod_name, y=1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_predict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Look at relationship between predicted and real with `cross_val_predict`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%%capture output --no-stdout --no-display\n",
    "# F-score + OLS pipeline\n",
    "fs_freg = SelectKBest(f_regression, k=10)\n",
    "lm_fs_freg = Pipeline([('fs_freg', fs_freg), \n",
    "                       ('lm', LinearRegression())])\n",
    "y_pred = cross_val_predict(lm_fs_freg, X_train, y_train, cv=age_skf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%capture output --no-stdout --no-display\n",
    "\n",
    "# Scores\n",
    "r2 = r2_score(y_train, y_pred)\n",
    "mse = mean_squared_error(y_train, y_pred)\n",
    "print(\"Scores from CV are R2: {}, MSE: {}\".format(r2, mse))\n",
    "\n",
    "sns.regplot(y_pred, y_train)\n",
    "plt.xlabel('Predicted Age with f-score feature elimination + OLS')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now transforming the predicted variable to account for the possible power law  \n",
    "https://scikit-learn.org/stable/modules/generated/sklearn.compose.TransformedTargetRegressor.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.compose import TransformedTargetRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# F-score + OLS pipeline + log transformer\n",
    "fs_freg = SelectKBest(f_regression, k=10)\n",
    "lm_fs_freg = Pipeline([('fs_freg', fs_freg), \n",
    "                       ('lm', LinearRegression())])\n",
    "ttr = TransformedTargetRegressor(regressor=lm_fs_freg, func=np.log, inverse_func=np.exp)\n",
    "\n",
    "y_pred = cross_val_predict(ttr, X_train, y_train, cv=age_skf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Scores\n",
    "r2 = r2_score(np.log(y_train), np.log(y_pred))\n",
    "mse = mean_squared_error(np.log(y_train), np.log(y_pred))\n",
    "print(\"Scores from CV are R2: {}, MSE: {}\".format(r2, mse))\n",
    "\n",
    "sns.regplot(np.log(y_pred), np.log(y_train))\n",
    "plt.xlabel('Predicted Age with f-score feature elimination + OLS')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Using penalized regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lasso = Lasso(random_state=42, max_iter=10000)\n",
    "alpha_range = np.logspace(-4, -0.5, 30)\n",
    "alpha_range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's investigate the change in performance score (R2, MSE) as we vary the penalty parameters alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Warning:** the cell below can take some time to run. Could be useful to reduce the range of the `alpha` hyperparameter (e.g. `alpha_range = np.logspace(-2, -0.5, 15)`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "l1_train_re, l1_test_r2 = validation_curve(lasso, X_train, y_train, param_name=\"alpha\", cv=age_skf, \n",
    "                                           param_range=alpha_range, scoring='r2', n_jobs=-1)\n",
    "l1_train_mse, l1_test_mse = validation_curve(lasso, X_train, y_train, param_name=\"alpha\", \n",
    "                                             cv=age_skf, param_range=alpha_range,  \n",
    "                                             scoring='neg_mean_squared_error', n_jobs=-1)\n",
    "# Save results in dataframe for display\n",
    "valcurve_df = pd.concat([valcurve_to_df(l1_train_re, l1_test_r2, alpha_range, 'lasso', 'r2'),\n",
    "                         valcurve_to_df(l1_train_mse, l1_test_mse, alpha_range, 'lasso', 'mse')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "for mod_name in valcurve_df['model'].unique():\n",
    "    kf_df = valcurve_df.loc[valcurve_df['model'] == mod_name]\n",
    "    with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "        g = sns.catplot(x='param', y='score', hue='stage', row='scorer', data=kf_df, kind='point',\n",
    "                        sharey=False, hue_order =[\"train\", \"test\"], height=3, aspect=3.0)\n",
    "        #plt.xticks(range(10))\n",
    "        g.set_xticklabels(alpha_range, rotation=90)\n",
    "        #g.axes[0,0].set_ylim(-0.2, 1.1)\n",
    "        #g.axes[0,1].set_ylim(-45, 1)\n",
    "        g.fig.suptitle(mod_name, y=1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lasso_choice = Lasso(alpha=0.045, random_state=42, max_iter=100000)\n",
    "y_pred = cross_val_predict(lasso_choice, X_train, y_train, cv=age_skf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Scores\n",
    "r2 = r2_score((y_train), (y_pred))\n",
    "mse = mean_squared_error((y_train), (y_pred))\n",
    "print(\"Scores from CV are R2: {}, MSE: {}\".format(r2, mse))\n",
    "\n",
    "sns.regplot((y_pred), (y_train))\n",
    "plt.xlabel('Predicted Age with lasso')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Better approach: nested cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LassoCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LassoCV?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lasso_cv = LassoCV(alphas=alpha_range, random_state=42, max_iter=10000, cv=10, n_jobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%capture output --no-stdout --no-display\n",
    "nestedcv_results = []\n",
    "best_alphas = []\n",
    "\n",
    "skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)\n",
    "for i_f, (ix_train, ix_test) in enumerate(skf.split(X_train, age_class_train)):\n",
    "        lasso_cv.fit(X_train[ix_train], y_train.iloc[ix_train])\n",
    "        best_alphas.append(lasso_cv.alpha_)\n",
    "        # Test on both the testing fold and training folds to check for over-/under-fitting\n",
    "        y_pred_train = lasso_cv.predict(X_train[ix_train])\n",
    "        y_pred_test = lasso_cv.predict(X_train[ix_test])\n",
    "        # R2\n",
    "        nestedcv_results.append({'model': 'lasso', 'fold': i_f, 'stage': 'train', 'scorer': 'r2', \n",
    "                                 'val': r2_score(y_train.iloc[ix_train], y_pred_train)})\n",
    "        nestedcv_results.append({'model': 'lasso', 'fold': i_f, 'stage': 'test', 'scorer': 'r2', \n",
    "                                 'val': r2_score(y_train.iloc[ix_test], y_pred_test)})\n",
    "        # MSE\n",
    "        nestedcv_results.append({'model': 'lasso', 'fold': i_f, 'stage': 'train', 'scorer': 'MSE', \n",
    "                                 'val': -mean_squared_error(y_train.iloc[ix_train], y_pred_train)})\n",
    "        nestedcv_results.append({'model': 'lasso', 'fold': i_f, 'stage': 'test', 'scorer': 'MSE', \n",
    "                                 'val': -mean_squared_error(y_train.iloc[ix_test], y_pred_test)})\n",
    "nestedcv_results_df = pd.DataFrame(nestedcv_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "mod_name = 'lasso'\n",
    "skf_df = nestedcv_results_df.loc[nestedcv_results_df['model'] == mod_name]\n",
    "with sns.plotting_context(\"notebook\", font_scale=1.2):\n",
    "    # Barploat\n",
    "    g = sns.catplot(x=\"fold\", y=\"val\", hue=\"stage\", col=\"scorer\", data=skf_df, kind=\"bar\", \n",
    "                    sharey=False, height=3, aspect=1.5)\n",
    "    g.axes[0,0].set_ylim(-0.2,1)\n",
    "    g.fig.suptitle(mod_name, y=1.05)\n",
    "    # Boxplot\n",
    "    g = sns.catplot(x=\"model\", y=\"val\", hue=\"stage\", col=\"scorer\", data=skf_results_df, kind=\"box\",\n",
    "                    sharey=False, height=3, aspect=1.5)\n",
    "    g.axes[0, 0].set_ylim(-0.2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "best_alphas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 2: analyzing our final model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Fit on whole (training) dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lasso_choice = Lasso(alpha=0.045, random_state=42, max_iter=100000)\n",
    "lasso_choice.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Test and score on hold-out test dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "y_pred = lasso_choice.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Scores\n",
    "r2 = r2_score(y_test, y_pred)\n",
    "mse = mean_squared_error(y_test, y_pred)\n",
    "print(\"Scores from CV are R2: {}, MSE: {}\".format(r2, mse))\n",
    "\n",
    "sns.regplot(y_pred, y_test)\n",
    "plt.xlabel('Predicted Age with final lasso')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Interpret model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lasso_choice.coef_.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.stem(np.arange(len(lasso_choice.coef_)), lasso_choice.coef_)\n",
    "plt.xlabel('Feature')\n",
    "plt.ylabel('Beta coefficient')\n",
    "plt.title(\"Final lasso model with {} parameters\".format(np.count_nonzero(lasso_choice.coef_)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Want to see a brain plot of the connections prediction age? Let's dive into Nilearn!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Nilearn usage example for resting-state fMRI connectivity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For this section we will follow closely the notebook created by Jake Vogel and presented at the 2020 Brainhack School: https://github.com/neurodatascience/course-materials-2020/blob/master/lectures/14-may/03-intro-to-machine-learning/ML_Regression_Tutorial.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Given $N$ regions of the brain, a functional connectivity matrix is typically created by:\n",
    "* choosing a partition of the brain into $N$ regions (often non-overlapping)\n",
    "* computing the correlation in time between each possible pair of regions $(i, j)$\n",
    "* created the associated matrix where each element $r(i(t),j(t))$ of the matrix is the correlation $r$ between:\n",
    "    * the signal $i(t)$ in region $i$\n",
    "    * the signal $j(t)$ in region $j$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Standard pearson correlation is often used, so since $r(i(t), j(t)) = r(j(t), i(t))$ the matrix is symmetric, and since $r(k(t), k(t)) = 1$, the diagonal is made of $1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The data required to generate the connectivity matrices are fMRI data already preprocessed for motion, etc. The data is typically 4D: a 3D brain volume within which each voxel (3D pixel) includes a time course of that voxel activity. The activity is typically \"rest\": subjects just have to lie in the scanner during the scan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "With `nilearn` we can use the preprocessed data and transform them into data ready to be used with machine learning (the feature matrix `X` and label vector `y` we used before) with the following steps:\n",
    "* Choose an atlas of $N$ regions (already aligned with all subjects)\n",
    "* For each subject:\n",
    "    * Extract the mean (or median) signal within each atlas region\n",
    "    * Compute the correlation between each pair of regions\n",
    "    * Transform the matrix into a vector: the features for that subject\n",
    "* Concatenate all vectors as rows of the feature matrix `X`\n",
    "* Create the label vector `y`, one value per subject:\n",
    "    * a continuous value for regression (as we did before)\n",
    "    * *or* a categorical value for classification (cf e.g. of classification in section 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We will use the 64-region version of the multi-resolution MIST atlas (cf. https://mniopenresearch.org/articles/1-3)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Get and view the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn import plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "atlases = datasets.fetch_atlas_basc_multiscale_2015(version='sym', \n",
    "                                                    data_dir='datasets/age_netmats')\n",
    "atlas64_path = atlases.scale064\n",
    "# Get the 3D coordinates of the atlas regions for future plotting\n",
    "atlas64_coords = plotting.find_parcellation_cut_coords(atlas64_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plotting.plot_roi(atlas64_path, draw_cross=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plotting.view_img(atlas64_path, cmap=plotting.cm.bwr, symmetric_cmap=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### MRI data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's have a look at the first subject"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "### Already done earlier:\n",
    "# datasets.fetch_development_fmri?\n",
    "# data = datasets.fetch_development_fmri(data_dir='datasets/age_netmats')\n",
    "s0_fmri_files = data.func[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn import image "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "s0_avg_vol = image.mean_img(s0_fmri_files)\n",
    "display = plotting.plot_epi(s0_avg_vol, cmap=\"gray\", vmin=0, vmax=1000)\n",
    "#display.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plotting.view_img?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%%capture output --no-stdout --no-display\n",
    "plotting.view_img(s0_avg_vol, cmap=\"binary\", symmetric_cmap=False, threshold=0, vmin=0, vmax=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Create data in machine learning format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Extract the time series from each brain region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiLabelsMasker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A `NiftiLabelsMasker` object can be used not only to extract the signal from each region but also to \"clean it\" by regressing out \"confounds\".  \n",
    "In this case a confound file needs to be passed to `NiftiLabelsMasker `."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "s0_confounds_file = data.confounds[0]\n",
    "s0_confounds_df = pd.read_csv(s0_confounds_file, sep='\\t')\n",
    "s0_confounds_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "masker = NiftiLabelsMasker(labels_img=atlas64_path, standardize=True, \n",
    "                           memory='nilearn_cache',  verbose=1)\n",
    "s0_time_series = masker.fit_transform(s0_fmri_files, confounds=s0_confounds_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Can you guess the size of the resulting `time_series` 2D array knowing that one dimension is the number of activity time points (168)?  \n",
    "(Hint: try to remember the number of regions in the atlas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "s0_time_series.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can visualize the mean signal over some ROIS for subject 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "plt.plot(s0_time_series[:, 0], alpha=0.7, label='region 0')\n",
    "plt.plot(s0_time_series[:, 1], alpha=0.7, label='region 1')\n",
    "plt.plot(s0_time_series[:, 2], alpha=0.7, label='region 2')\n",
    "plt.title('Voxel Time Series')\n",
    "plt.xlabel('Time points (one 3D volume per time point)')\n",
    "plt.ylabel('Normalized signal')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Are the time series correlated? Let's quantify this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.corrcoef(s0_time_series[:, :3].T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Compute correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from nilearn.connectome import ConnectivityMeasure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "correlation_measure = ConnectivityMeasure(kind='correlation')\n",
    "s0_conmat = correlation_measure.fit_transform([s0_time_series])[0]\n",
    "s0_conmat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plotting.plot_matrix(s0_conmat, figure=(10, 8), labels=range(s0_time_series.shape[-1]),\n",
    "                     vmax=0.8, vmin=-0.8, reorder=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plotting.plot_connectome(s0_conmat, atlas64_coords, title='Functional connectivity of subject 0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The next step is to process all the subjects. We will not do it as we already analyzed these data. The code is left below for reference.  \n",
    "Note: the option `vectorize` is set to `True` in order to have the whole correlation matrix as a single vector of features (for each subject)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "masker = NiftiLabelsMasker(labels_img=atlas64_path, standardize=True, \n",
    "                           memory='nilearn_cache', verbose=0)\n",
    "correlation_measure = ConnectivityMeasure(kind='correlation', vectorize=True, discard_diagonal=True)\n",
    "vectorized_conmats = []\n",
    "#for i_s, s in enumerate(data.func):\n",
    "#    print(\"=== Processing subject {:03d}\".format(i_s))\n",
    "#    # Get the mean time series within each brain region\n",
    "#    s_time_series = masker.fit_transform(s, confounds=data.confounds[i_s])\n",
    "#    # Compute the correlation between each pair of region to create the connectivity matrix\n",
    "#    s_conmat = correlation_measure.fit_transform([s_time_series])[0]\n",
    "#    # Add to the list of vectorized connectivity matrices\n",
    "#    vectorized_conmats.append(s_conmat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize our results from the previous section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We'll first transform the vectorized `lasso` coefficients back to the matrix format (this is the inverse transform of what we just did)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Refit correlation to make sure the mat --> vec transform is calculated\n",
    "correlation_measure = ConnectivityMeasure(kind='correlation', vectorize=True, discard_diagonal=True)\n",
    "s0_conmat_vec = correlation_measure.fit_transform([s0_time_series])[0]\n",
    "# Now we can apply the inverse transform\n",
    "lasso_coeff_mat = correlation_measure.inverse_transform([lasso_choice.coef_])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plotting.plot_matrix(lasso_coeff_mat, figure=(10, 8), labels=range(lasso_coeff_mat.shape[0]),\n",
    "                     reorder=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plotting.plot_connectome(lasso_coeff_mat, atlas64_coords, colorbar=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python [conda env:ni38]",
   "language": "python",
   "name": "conda-env-ni38-py"
  },
  "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.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 760,
   "position": {
    "height": "40px",
    "left": "1650px",
    "right": "20px",
    "top": "120px",
    "width": "250px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "none",
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}