multivariate_BWAS_replicability_analysis_CT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Analyzing replicability of structure-based multivariate BWAS on the Human Connectome Project dataset\n",
    "\n",
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-08-03T20:04:15.431840Z",
     "start_time": "2021-08-03T20:04:14.753565Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.model_selection import KFold, train_test_split, cross_val_predict, GridSearchCV\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.decomposition import PCA\n",
    "from joblib import Parallel, delayed\n",
    "from mlxtend.evaluate import permutation_test\n",
    "import fnmatch\n",
    "sns.set(rc={\"figure.figsize\":(4, 2)})\n",
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Load HCP data\n",
    "\n",
    "We load functional network matrices (netmats) from the HCP1200-release, as published on connectomeDB: https://db.humanconnectome.org/\n",
    "Due to licensoing issues, data is not supplied with the repository, but can be downloaded from the ConnectomeDB.\n",
    "See [hcp_data/readme.md](hcp_data/readme.md) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "data": {
      "text/plain": "        Release Acquisition Gender    Age  3T_Full_MR_Compl  T1_Count  \\\nSubject                                                                 \n100004     S900         Q06      M  22-25             False         0   \n100206     S900         Q11      M  26-30              True         1   \n100307       Q1         Q01      F  26-30              True         1   \n100408       Q3         Q03      M  31-35              True         1   \n100610     S900         Q08      M  26-30              True         2   \n...         ...         ...    ...    ...               ...       ...   \n992774       Q2         Q02      M  31-35              True         2   \n993675     S900         Q09      F  26-30              True         2   \n994273     S500         Q06      M  26-30              True         1   \n995174    S1200         Q13      M  22-25             False         1   \n996782     S900         Q08      F  26-30              True         2   \n\n         T2_Count  3T_RS-fMRI_Count  3T_RS-fMRI_PctCompl  3T_Full_Task_fMRI  \\\nSubject                                                                       \n100004          0                 0                  0.0              False   \n100206          1                 4                100.0               True   \n100307          1                 4                100.0               True   \n100408          1                 4                100.0               True   \n100610          1                 4                100.0               True   \n...           ...               ...                  ...                ...   \n992774          2                 4                100.0               True   \n993675          2                 4                100.0               True   \n994273          1                 4                100.0               True   \n995174          1                 2                  0.0               True   \n996782          2                 4                100.0               True   \n\n         ...  Odor_Unadj  Odor_AgeAdj  PainIntens_RawScore  PainInterf_Tscore  \\\nSubject  ...                                                                    \n100004   ...      101.12        86.45                  2.0               45.9   \n100206   ...      108.79        97.19                  1.0               49.7   \n100307   ...      101.12        86.45                  0.0               38.6   \n100408   ...      108.79        98.04                  2.0               52.6   \n100610   ...      122.25       110.45                  0.0               38.6   \n...      ...         ...          ...                  ...                ...   \n992774   ...      122.25       111.41                  4.0               50.1   \n993675   ...      122.25       110.45                  0.0               38.6   \n994273   ...      122.25       111.41                  7.0               63.8   \n995174   ...       88.61        64.58                  3.0               50.1   \n996782   ...      108.79        97.19                  0.0               38.6   \n\n         Taste_Unadj  Taste_AgeAdj  Mars_Log_Score  Mars_Errs  Mars_Final  \\\nSubject                                                                     \n100004        107.17        105.31            1.80        0.0        1.80   \n100206         72.63         72.03            1.84        0.0        1.84   \n100307         71.69         71.76            1.76        0.0        1.76   \n100408        114.01        113.59            1.76        2.0        1.68   \n100610         84.84         85.31            1.92        1.0        1.88   \n...              ...           ...             ...        ...         ...   \n992774        107.17        103.55            1.76        0.0        1.76   \n993675         84.07         84.25            1.80        1.0        1.76   \n994273        110.65        109.73            1.80        1.0        1.76   \n995174        117.16        117.40            1.80        0.0        1.80   \n996782         75.43         73.72            1.84        0.0        1.84   \n\n          age  \nSubject        \n100004   23.5  \n100206   28.0  \n100307   28.0  \n100408   33.0  \n100610   28.0  \n...       ...  \n992774   33.0  \n993675   28.0  \n994273   28.0  \n995174   23.5  \n996782   28.0  \n\n[1206 rows x 582 columns]",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Release</th>\n      <th>Acquisition</th>\n      <th>Gender</th>\n      <th>Age</th>\n      <th>3T_Full_MR_Compl</th>\n      <th>T1_Count</th>\n      <th>T2_Count</th>\n      <th>3T_RS-fMRI_Count</th>\n      <th>3T_RS-fMRI_PctCompl</th>\n      <th>3T_Full_Task_fMRI</th>\n      <th>...</th>\n      <th>Odor_Unadj</th>\n      <th>Odor_AgeAdj</th>\n      <th>PainIntens_RawScore</th>\n      <th>PainInterf_Tscore</th>\n      <th>Taste_Unadj</th>\n      <th>Taste_AgeAdj</th>\n      <th>Mars_Log_Score</th>\n      <th>Mars_Errs</th>\n      <th>Mars_Final</th>\n      <th>age</th>\n    </tr>\n    <tr>\n      <th>Subject</th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>100004</th>\n      <td>S900</td>\n      <td>Q06</td>\n      <td>M</td>\n      <td>22-25</td>\n      <td>False</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0.0</td>\n      <td>False</td>\n      <td>...</td>\n      <td>101.12</td>\n      <td>86.45</td>\n      <td>2.0</td>\n      <td>45.9</td>\n      <td>107.17</td>\n      <td>105.31</td>\n      <td>1.80</td>\n      <td>0.0</td>\n      <td>1.80</td>\n      <td>23.5</td>\n    </tr>\n    <tr>\n      <th>100206</th>\n      <td>S900</td>\n      <td>Q11</td>\n      <td>M</td>\n      <td>26-30</td>\n      <td>True</td>\n      <td>1</td>\n      <td>1</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>108.79</td>\n      <td>97.19</td>\n      <td>1.0</td>\n      <td>49.7</td>\n      <td>72.63</td>\n      <td>72.03</td>\n      <td>1.84</td>\n      <td>0.0</td>\n      <td>1.84</td>\n      <td>28.0</td>\n    </tr>\n    <tr>\n      <th>100307</th>\n      <td>Q1</td>\n      <td>Q01</td>\n      <td>F</td>\n      <td>26-30</td>\n      <td>True</td>\n      <td>1</td>\n      <td>1</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>101.12</td>\n      <td>86.45</td>\n      <td>0.0</td>\n      <td>38.6</td>\n      <td>71.69</td>\n      <td>71.76</td>\n      <td>1.76</td>\n      <td>0.0</td>\n      <td>1.76</td>\n      <td>28.0</td>\n    </tr>\n    <tr>\n      <th>100408</th>\n      <td>Q3</td>\n      <td>Q03</td>\n      <td>M</td>\n      <td>31-35</td>\n      <td>True</td>\n      <td>1</td>\n      <td>1</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>108.79</td>\n      <td>98.04</td>\n      <td>2.0</td>\n      <td>52.6</td>\n      <td>114.01</td>\n      <td>113.59</td>\n      <td>1.76</td>\n      <td>2.0</td>\n      <td>1.68</td>\n      <td>33.0</td>\n    </tr>\n    <tr>\n      <th>100610</th>\n      <td>S900</td>\n      <td>Q08</td>\n      <td>M</td>\n      <td>26-30</td>\n      <td>True</td>\n      <td>2</td>\n      <td>1</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>122.25</td>\n      <td>110.45</td>\n      <td>0.0</td>\n      <td>38.6</td>\n      <td>84.84</td>\n      <td>85.31</td>\n      <td>1.92</td>\n      <td>1.0</td>\n      <td>1.88</td>\n      <td>28.0</td>\n    </tr>\n    <tr>\n      <th>...</th>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n    </tr>\n    <tr>\n      <th>992774</th>\n      <td>Q2</td>\n      <td>Q02</td>\n      <td>M</td>\n      <td>31-35</td>\n      <td>True</td>\n      <td>2</td>\n      <td>2</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>122.25</td>\n      <td>111.41</td>\n      <td>4.0</td>\n      <td>50.1</td>\n      <td>107.17</td>\n      <td>103.55</td>\n      <td>1.76</td>\n      <td>0.0</td>\n      <td>1.76</td>\n      <td>33.0</td>\n    </tr>\n    <tr>\n      <th>993675</th>\n      <td>S900</td>\n      <td>Q09</td>\n      <td>F</td>\n      <td>26-30</td>\n      <td>True</td>\n      <td>2</td>\n      <td>2</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>122.25</td>\n      <td>110.45</td>\n      <td>0.0</td>\n      <td>38.6</td>\n      <td>84.07</td>\n      <td>84.25</td>\n      <td>1.80</td>\n      <td>1.0</td>\n      <td>1.76</td>\n      <td>28.0</td>\n    </tr>\n    <tr>\n      <th>994273</th>\n      <td>S500</td>\n      <td>Q06</td>\n      <td>M</td>\n      <td>26-30</td>\n      <td>True</td>\n      <td>1</td>\n      <td>1</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>122.25</td>\n      <td>111.41</td>\n      <td>7.0</td>\n      <td>63.8</td>\n      <td>110.65</td>\n      <td>109.73</td>\n      <td>1.80</td>\n      <td>1.0</td>\n      <td>1.76</td>\n      <td>28.0</td>\n    </tr>\n    <tr>\n      <th>995174</th>\n      <td>S1200</td>\n      <td>Q13</td>\n      <td>M</td>\n      <td>22-25</td>\n      <td>False</td>\n      <td>1</td>\n      <td>1</td>\n      <td>2</td>\n      <td>0.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>88.61</td>\n      <td>64.58</td>\n      <td>3.0</td>\n      <td>50.1</td>\n      <td>117.16</td>\n      <td>117.40</td>\n      <td>1.80</td>\n      <td>0.0</td>\n      <td>1.80</td>\n      <td>23.5</td>\n    </tr>\n    <tr>\n      <th>996782</th>\n      <td>S900</td>\n      <td>Q08</td>\n      <td>F</td>\n      <td>26-30</td>\n      <td>True</td>\n      <td>2</td>\n      <td>2</td>\n      <td>4</td>\n      <td>100.0</td>\n      <td>True</td>\n      <td>...</td>\n      <td>108.79</td>\n      <td>97.19</td>\n      <td>0.0</td>\n      <td>38.6</td>\n      <td>75.43</td>\n      <td>73.72</td>\n      <td>1.84</td>\n      <td>0.0</td>\n      <td>1.84</td>\n      <td>28.0</td>\n    </tr>\n  </tbody>\n</table>\n<p>1206 rows × 582 columns</p>\n</div>"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# HCP data can be obtained from the connectomeDB\n",
    "# data is not part of this repository\n",
    "behavior = pd.read_csv('hcp_data/hcp1200_behavioral_data.csv')\n",
    "behavior = behavior.set_index('Subject', drop=True)\n",
    "\n",
    "# convert age to numeric\n",
    "age = []\n",
    "for s in behavior['Age']:\n",
    "    if s == '36+':\n",
    "        age.append(36)\n",
    "    else:\n",
    "        split = s.split(sep='-')\n",
    "        age.append(np.mean((float(split[0]), float(split[1]))))\n",
    "\n",
    "behavior['age'] = age\n",
    "behavior"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Function to prepare target variable\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "def create_data(target='CogTotalComp_AgeAdj', feature_data=\"FS_*_Thck\"):\n",
    "    # it's a good practice to use pandas for merging, messing up subject order can be painful\n",
    "\n",
    "    features = fnmatch.filter(behavior.columns.values, feature_data)\n",
    "\n",
    "    df = behavior.dropna(subset = [target] + features)\n",
    "    y = df[target].values\n",
    "    X = df[features].values\n",
    "    return X, y"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [],
   "source": [
    "features = fnmatch.filter(behavior.columns.values, \"FS_*_Thck\")\n",
    "len(features)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Function implementing a single bootstrap iteration\n",
    "\n",
    "We define a workhorse function which:\n",
    "- randomly samples the discovery and the replication datasets,\n",
    "- creates cross-validated estimates of predictive performance within the discovery sample\n",
    "- finalizes the model by fitting it to the whole discovery sample (overfits the discovery but not the replication sample)\n",
    "- use it to predict the replication sample"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "outputs": [],
   "source": [
    "def bootstrap_workhorse(X, y, sample_size, model, random_state, shuffle_y=False):\n",
    "\n",
    "    #create discovery and replication samples by random sampling from the whole dataset (without replacement)\n",
    "\n",
    "    if shuffle_y:\n",
    "        rng = np.random.default_rng(random_state)\n",
    "        y = rng.permutation(y)\n",
    "\n",
    "    X_discovery, X_replication, y_discovery, y_replication = train_test_split(X, y, train_size=sample_size, test_size=sample_size, shuffle=True, random_state=random_state)\n",
    "\n",
    "    cv = KFold(10)\n",
    "    # obtain cross-validated predictions in the discovery sample\n",
    "\n",
    "    predicted_discovery_cv = np.zeros_like(y_discovery)\n",
    "    cor_per_fold = np.zeros(cv.n_splits)\n",
    "    i = 0\n",
    "    for train, test in cv.split(X=X_discovery, y=y_discovery):\n",
    "        model.fit(X=X_discovery[train], y=y_discovery[train])\n",
    "        predicted_discovery_cv[test] = model.predict(X=X_discovery[test])\n",
    "        cor_per_fold[i] = np.corrcoef(y_discovery[test], predicted_discovery_cv[test])[0,1]\n",
    "        i += 1\n",
    "    # correlation between the cross-validated predictions and observations in the discovery sample\n",
    "    # this is the correct, unbiased estimate!\n",
    "    # calculated as mean test performance across all folds\n",
    "    r_disc_cv = np.mean(cor_per_fold)\n",
    "    # finalize model by training it on the full discovery sample (without cross-validation)\n",
    "    final_model = model.fit(X=X_discovery, y=y_discovery)\n",
    "    # obtain predictions with the final model on the discovery sample, note that this model actually overfits this sample.\n",
    "    # we do this only to demonstrate biased estimates\n",
    "    predicted_discovery_overfit = final_model.predict(X=X_discovery)\n",
    "    # here we obtain the biased effect size (r) estimates for demonstrational purposes\n",
    "    r_disc_overfit = np.corrcoef(predicted_discovery_overfit, y_discovery)[0, 1]\n",
    "\n",
    "    # We use the final model to predict the replication sample\n",
    "    # This is correct (no overfitting here), the final model did not see this data during training\n",
    "    predicted_replication = final_model.predict(X=X_replication)\n",
    "    # we obtain the out-of-sample prediction performance estimates\n",
    "    r_rep = np.corrcoef(predicted_replication, y_replication)[0, 1]\n",
    "\n",
    "    # below we calculate permutation-based p-values for all three effect size estimates (in-sample unbiased, in-sample biased, out-of-sample)\n",
    "    # (one sided tests, testing for positive correlation)\n",
    "    p_disc_cv = permutation_test(predicted_discovery_cv, y_discovery, method='approximate', num_rounds=1000, func=lambda x, y: np.corrcoef(x, y)[1][0],seed=random_state)\n",
    "    p_disc_overfit = permutation_test(predicted_discovery_overfit, y_discovery, method='approximate', num_rounds=1000, func=lambda x, y: np.corrcoef(x, y)[1][0],seed=random_state)\n",
    "    p_rep = permutation_test(predicted_replication, y_replication, method='approximate', num_rounds=1000, func=lambda x, y: np.corrcoef(x, y)[1][0],seed=random_state)\n",
    "    # return results\n",
    "    return r_disc_cv, r_disc_overfit, r_rep, p_disc_cv, p_disc_overfit, p_rep"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "All set, now we start the analysis."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Replicability with sample sizes n=50, 100, 200, 300 and max\n",
    "Here we train a few different models on 100 bootstrap samples.\n",
    "\n",
    "We aggregate the results of our workhorse function in `n_bootstrap`=100 bootstrap cases (run in parallel).\n",
    "\n",
    "The whole process is repeated for all sample sizes, fetaure_sets and target variables.\n",
    "\n",
    "## Here we test age and 5 cognitive variables, including 'cognitive ability' (the 'main' target variable in the target paper)\n",
    "- age: age group of the participants\n",
    "- CogTotalComp_AgeAdj: total cognitive ability\n",
    "- PMAT24_A_CR, : Fluid Intelligence (Penn Progressive Matrices)\n",
    "- CardSort_AgeAdj: Executive Function/Cognitive Flexibility (Dimensional Change Card Sort)\n",
    "- Flanker_AgeAdj: Executive Function/Inhibition (Flanker Task)\n",
    "- PicSeq_AgeAdj: Episodic Memory (Picture Sequence Memory)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Reproducing the PCA+SVR-based model from the target paper\n",
    "### Like in the target paper:\n",
    "- Both PCA and SVR are done inside the cross-validation\n",
    "- PCA reatains the firts k principal components that together explain 50% of the variance\n",
    "- scikit-learn makes sure that PCA is only fit for the training samples\n",
    "- both for the test sets (in the cross-validation) and the replication sample PCA is not re-fit, bt features are simply transformed with the already fit PCA"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 50\n",
      "0.08589268853895356 0.11814430558345144\n",
      "Replicability at alpha = 0.05 : 30.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 100\n",
      "0.09976372318981333 0.1321490781638016\n",
      "Replicability at alpha = 0.05 : 41.17647058823529 %\n",
      "Replicability at alpha = 0.01 : 5.88235294117647 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 200\n",
      "0.1319123546390201 0.12993985943621775\n",
      "Replicability at alpha = 0.05 : 59.45945945945946 %\n",
      "Replicability at alpha = 0.01 : 24.324324324324326 %\n",
      "Replicability at alpha = 0.005 : 16.216216216216218 %\n",
      "Replicability at alpha = 0.001 : 2.7027027027027026 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 300\n",
      "0.12171574361938152 0.12269194343762159\n",
      "Replicability at alpha = 0.05 : 77.5 %\n",
      "Replicability at alpha = 0.01 : 40.0 %\n",
      "Replicability at alpha = 0.005 : 20.0 %\n",
      "Replicability at alpha = 0.001 : 7.5 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age max\n",
      "0.14052260175749642 0.14440243478574197\n",
      "Replicability at alpha = 0.05 : 98.73417721518987 %\n",
      "Replicability at alpha = 0.01 : 60.75949367088608 %\n",
      "Replicability at alpha = 0.005 : 45.56962025316456 %\n",
      "Replicability at alpha = 0.001 : 11.39240506329114 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 50\n",
      "0.027807350316040817 0.06308898890639722\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 100\n",
      "0.022757951091348018 0.061284458973813294\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 200\n",
      "0.04462184321671844 0.04577939992501347\n",
      "Replicability at alpha = 0.05 : 50.0 %\n",
      "Replicability at alpha = 0.01 : 25.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 300\n",
      "0.05439263443559869 0.06513144477011525\n",
      "Replicability at alpha = 0.05 : 66.66666666666666 %\n",
      "Replicability at alpha = 0.01 : 22.22222222222222 %\n",
      "Replicability at alpha = 0.005 : 22.22222222222222 %\n",
      "Replicability at alpha = 0.001 : 11.11111111111111 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj max\n",
      "0.08067948568531891 0.0815027616931159\n",
      "Replicability at alpha = 0.05 : 56.09756097560976 %\n",
      "Replicability at alpha = 0.01 : 12.195121951219512 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 50\n",
      "0.008305362063431911 0.02353136316127152\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 100\n",
      "0.009558285589515252 0.012546573293026489\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 200\n",
      "0.035337364949865675 0.02646112111538681\n",
      "Replicability at alpha = 0.05 : 12.5 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 300\n",
      "0.030473825143593255 0.01688841104240886\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR max\n",
      "0.00685802000666485 0.008780290087947814\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 50\n",
      "0.0037964289881357556 0.0036090118905972445\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 100\n",
      "0.008171689667888534 0.0026079995876381865\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 200\n",
      "-0.004929600225644181 -0.009921153702658813\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 300\n",
      "-0.015493604931604406 -0.029831004294138302\n",
      "Replicability at alpha = 0.05 : nan %\n",
      "Replicability at alpha = 0.01 : nan %\n",
      "Replicability at alpha = 0.005 : nan %\n",
      "Replicability at alpha = 0.001 : nan %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj max\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:59: RuntimeWarning: invalid value encountered in long_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.010603773264953613 -0.019130199826790952\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 50\n",
      "0.03207685702787072 0.030830317612604786\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 100\n",
      "0.04021060189802274 0.03077808783073594\n",
      "Replicability at alpha = 0.05 : 12.5 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 200\n",
      "0.04137443676857757 0.04091365957382529\n",
      "Replicability at alpha = 0.05 : 22.22222222222222 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 300\n",
      "0.06497921296029413 0.06333840296663894\n",
      "Replicability at alpha = 0.05 : 42.10526315789473 %\n",
      "Replicability at alpha = 0.01 : 10.526315789473683 %\n",
      "Replicability at alpha = 0.005 : 5.263157894736842 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj max\n",
      "0.07888863748370917 0.08736464060154424\n",
      "Replicability at alpha = 0.05 : 58.333333333333336 %\n",
      "Replicability at alpha = 0.01 : 2.7777777777777777 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 50\n",
      "-0.002662800743093452 -0.015009832164857642\n",
      "Replicability at alpha = 0.05 : nan %\n",
      "Replicability at alpha = 0.01 : nan %\n",
      "Replicability at alpha = 0.005 : nan %\n",
      "Replicability at alpha = 0.001 : nan %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 100\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:59: RuntimeWarning: invalid value encountered in long_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.011373542001246479 0.004689831017528047\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 200\n",
      "0.02151173146450509 0.022081798461224947\n",
      "Replicability at alpha = 0.05 : 25.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 300\n",
      "0.01946957724852659 0.030930201944994077\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj max\n",
      "0.03010521006039739 0.03567733413187921\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "CPU times: user 7.81 s, sys: 720 ms, total: 8.53 s\n",
      "Wall time: 3min 32s\n"
     ]
    },
    {
     "data": {
      "text/plain": "     connectivity    model         target    n r_discovery_cv  \\\n0           FS_CT  PCA_SVR            age   50       0.089136   \n1           FS_CT  PCA_SVR            age   50            NaN   \n2           FS_CT  PCA_SVR            age   50            NaN   \n3           FS_CT  PCA_SVR            age   50      -0.178271   \n4           FS_CT  PCA_SVR            age   50      -0.032969   \n...           ...      ...            ...  ...            ...   \n2995        FS_CT  PCA_SVR  PicSeq_AgeAdj  556       0.012511   \n2996        FS_CT  PCA_SVR  PicSeq_AgeAdj  556       0.083369   \n2997        FS_CT  PCA_SVR  PicSeq_AgeAdj  556       0.055206   \n2998        FS_CT  PCA_SVR  PicSeq_AgeAdj  556      -0.029972   \n2999        FS_CT  PCA_SVR  PicSeq_AgeAdj  556      -0.126373   \n\n     r_discovery_overfit r_replication p_discovery_cv p_discovery_overfit  \\\n0               0.468271      0.289675       0.612388            0.000999   \n1               0.551551      0.117138       0.655345            0.000999   \n2               0.546213      0.155825       0.959041            0.000999   \n3               0.614155     -0.093113       0.982018            0.000999   \n4               0.403477      0.113985       0.691309            0.003996   \n...                  ...           ...            ...                 ...   \n2995            0.214361      0.039045       0.706294            0.000999   \n2996            0.227902      0.041203       0.086913            0.000999   \n2997            0.225934      0.043133       0.120879            0.000999   \n2998            0.153215      0.032032       0.809191            0.000999   \n2999            0.163645      0.003725            1.0            0.000999   \n\n     p_replication  \n0         0.023976  \n1         0.215784  \n2          0.14985  \n3         0.727273  \n4         0.230769  \n...            ...  \n2995      0.184815  \n2996      0.160839  \n2997      0.148851  \n2998      0.227772  \n2999      0.484515  \n\n[3000 rows x 10 columns]",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>connectivity</th>\n      <th>model</th>\n      <th>target</th>\n      <th>n</th>\n      <th>r_discovery_cv</th>\n      <th>r_discovery_overfit</th>\n      <th>r_replication</th>\n      <th>p_discovery_cv</th>\n      <th>p_discovery_overfit</th>\n      <th>p_replication</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>0.089136</td>\n      <td>0.468271</td>\n      <td>0.289675</td>\n      <td>0.612388</td>\n      <td>0.000999</td>\n      <td>0.023976</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>NaN</td>\n      <td>0.551551</td>\n      <td>0.117138</td>\n      <td>0.655345</td>\n      <td>0.000999</td>\n      <td>0.215784</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>NaN</td>\n      <td>0.546213</td>\n      <td>0.155825</td>\n      <td>0.959041</td>\n      <td>0.000999</td>\n      <td>0.14985</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>-0.178271</td>\n      <td>0.614155</td>\n      <td>-0.093113</td>\n      <td>0.982018</td>\n      <td>0.000999</td>\n      <td>0.727273</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>-0.032969</td>\n      <td>0.403477</td>\n      <td>0.113985</td>\n      <td>0.691309</td>\n      <td>0.003996</td>\n      <td>0.230769</td>\n    </tr>\n    <tr>\n      <th>...</th>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n    </tr>\n    <tr>\n      <th>2995</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.012511</td>\n      <td>0.214361</td>\n      <td>0.039045</td>\n      <td>0.706294</td>\n      <td>0.000999</td>\n      <td>0.184815</td>\n    </tr>\n    <tr>\n      <th>2996</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.083369</td>\n      <td>0.227902</td>\n      <td>0.041203</td>\n      <td>0.086913</td>\n      <td>0.000999</td>\n      <td>0.160839</td>\n    </tr>\n    <tr>\n      <th>2997</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.055206</td>\n      <td>0.225934</td>\n      <td>0.043133</td>\n      <td>0.120879</td>\n      <td>0.000999</td>\n      <td>0.148851</td>\n    </tr>\n    <tr>\n      <th>2998</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>-0.029972</td>\n      <td>0.153215</td>\n      <td>0.032032</td>\n      <td>0.809191</td>\n      <td>0.000999</td>\n      <td>0.227772</td>\n    </tr>\n    <tr>\n      <th>2999</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>-0.126373</td>\n      <td>0.163645</td>\n      <td>0.003725</td>\n      <td>1.0</td>\n      <td>0.000999</td>\n      <td>0.484515</td>\n    </tr>\n  </tbody>\n</table>\n<p>3000 rows × 10 columns</p>\n</div>"
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "random_state = 42\n",
    "n_bootstrap = 100\n",
    "\n",
    "features = {\n",
    "    'FS_CT': \"FS_*_Thck\"\n",
    "}\n",
    "\n",
    "models = {\n",
    "    'PCA_SVR': Pipeline([('pca', PCA(n_components=0.5)),\n",
    "                         ('svr', SVR())])\n",
    "\n",
    "}\n",
    "\n",
    "# We aggregate all results here:\n",
    "df = pd.DataFrame(columns=['connectivity','model','target','n','r_discovery_cv','r_discovery_overfit','r_replication','p_discovery_cv','p_discovery_overfit','p_replication'])\n",
    "\n",
    "for feature_set in features:\n",
    "    for model in models:\n",
    "        for target_var in ['age', 'CogTotalComp_AgeAdj', 'PMAT24_A_CR', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'PicSeq_AgeAdj']:\n",
    "            for sample_size in [50, 100, 200, 300, 'max']:\n",
    "\n",
    "                print('*****************************************************************')\n",
    "                print(feature_set, model, target_var, sample_size)\n",
    "\n",
    "                X, y = create_data(target=target_var, feature_data=features[feature_set])\n",
    "\n",
    "                if sample_size=='max':\n",
    "                    sample_size = int(len(y)/2)\n",
    "\n",
    "                # create random seeds for each bootstrap iteration for reproducibility\n",
    "                rng = np.random.default_rng(random_state)\n",
    "                random_sates = rng.integers(np.iinfo(np.int32).max, size=n_bootstrap)\n",
    "\n",
    "                # run bootstrap iterations in parallel\n",
    "                r_discovery_cv, r_discovery_overfit, r_replication, p_discovery_cv, p_discovery_overfit, p_replication = zip(\n",
    "                    *Parallel(n_jobs=-1)(\n",
    "                    delayed(bootstrap_workhorse)(X, y, sample_size, models[model], seed) for seed in random_sates))\n",
    "\n",
    "                tmp_data_frame = pd.DataFrame({\n",
    "                    'connectivity' : feature_set,\n",
    "                    'model' : model,\n",
    "                    'target' : target_var,\n",
    "                    'n' : sample_size,\n",
    "                    'r_discovery_cv': r_discovery_cv,\n",
    "                    'r_discovery_overfit': r_discovery_overfit,\n",
    "                    'r_replication': r_replication,\n",
    "                    'p_discovery_cv': p_discovery_cv,\n",
    "                    'p_discovery_overfit': p_discovery_overfit,\n",
    "                    'p_replication': p_replication\n",
    "                })\n",
    "                #sns.scatterplot(x='r_replication', y='r_discovery_cv', data=tmp_data_frame)\n",
    "                #plt.ylabel('in-sample (r)')\n",
    "                #plt.xlabel('out-of-sample (r_pred)')\n",
    "                #plt.show()\n",
    "                print(tmp_data_frame.r_discovery_cv.mean(), tmp_data_frame.r_replication.mean())\n",
    "\n",
    "                for alpha in [0.05, 0.01, 0.005, 0.001]:\n",
    "                    print('Replicability at alpha =', alpha, ':',\n",
    "                          (tmp_data_frame.loc[tmp_data_frame['p_discovery_cv']<alpha,'p_replication']<alpha).sum() / (tmp_data_frame['p_discovery_cv']<0.05).sum() * 100, '%')\n",
    "\n",
    "                df = pd.concat((df, tmp_data_frame))\n",
    "                df.reset_index(drop=True, inplace=True)\n",
    "                df.to_csv('res/results_thck_PCA_SVR.csv')\n",
    "\n",
    "df"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Now we fit a simple Ridge regression\n",
    "(no feature selection, no hyperparameter optimization)\n",
    "This can be expected to perform better on low samples than SVR."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*****************************************************************\n",
      "FS_CT ridge age 50\n",
      "0.15322562774488466 0.16690836840736034\n",
      "Replicability at alpha = 0.05 : 55.00000000000001 %\n",
      "Replicability at alpha = 0.01 : 10.0 %\n",
      "Replicability at alpha = 0.005 : 5.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge age 100\n",
      "0.15399594921982113 0.1668875020529473\n",
      "Replicability at alpha = 0.05 : 58.97435897435898 %\n",
      "Replicability at alpha = 0.01 : 12.82051282051282 %\n",
      "Replicability at alpha = 0.005 : 2.564102564102564 %\n",
      "Replicability at alpha = 0.001 : 2.564102564102564 %\n",
      "*****************************************************************\n",
      "FS_CT ridge age 200\n",
      "0.2227858933204806 0.22609433664584924\n",
      "Replicability at alpha = 0.05 : 94.31818181818183 %\n",
      "Replicability at alpha = 0.01 : 65.9090909090909 %\n",
      "Replicability at alpha = 0.005 : 52.27272727272727 %\n",
      "Replicability at alpha = 0.001 : 26.136363636363637 %\n",
      "*****************************************************************\n",
      "FS_CT ridge age 300\n",
      "0.24120584624743915 0.2515387452025219\n",
      "Replicability at alpha = 0.05 : 98.93617021276596 %\n",
      "Replicability at alpha = 0.01 : 90.42553191489363 %\n",
      "Replicability at alpha = 0.005 : 87.2340425531915 %\n",
      "Replicability at alpha = 0.001 : 68.08510638297872 %\n",
      "*****************************************************************\n",
      "FS_CT ridge age max\n",
      "0.26225660935039835 0.2804729981031089\n",
      "Replicability at alpha = 0.05 : 100.0 %\n",
      "Replicability at alpha = 0.01 : 100.0 %\n",
      "Replicability at alpha = 0.005 : 100.0 %\n",
      "Replicability at alpha = 0.001 : 99.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CogTotalComp_AgeAdj 50\n",
      "0.050850174162353226 0.08536461668579186\n",
      "Replicability at alpha = 0.05 : 4.761904761904762 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CogTotalComp_AgeAdj 100\n",
      "0.10321259690671888 0.11834467765925837\n",
      "Replicability at alpha = 0.05 : 42.857142857142854 %\n",
      "Replicability at alpha = 0.01 : 14.285714285714285 %\n",
      "Replicability at alpha = 0.005 : 9.523809523809524 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CogTotalComp_AgeAdj 200\n",
      "0.1489345823380713 0.1531064890802809\n",
      "Replicability at alpha = 0.05 : 79.66101694915254 %\n",
      "Replicability at alpha = 0.01 : 23.728813559322035 %\n",
      "Replicability at alpha = 0.005 : 11.864406779661017 %\n",
      "Replicability at alpha = 0.001 : 1.694915254237288 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CogTotalComp_AgeAdj 300\n",
      "0.17062951787445102 0.17614166737215428\n",
      "Replicability at alpha = 0.05 : 94.04761904761905 %\n",
      "Replicability at alpha = 0.01 : 61.904761904761905 %\n",
      "Replicability at alpha = 0.005 : 48.80952380952381 %\n",
      "Replicability at alpha = 0.001 : 22.61904761904762 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CogTotalComp_AgeAdj max\n",
      "0.21520251020327877 0.2153618683467791\n",
      "Replicability at alpha = 0.05 : 100.0 %\n",
      "Replicability at alpha = 0.01 : 99.0 %\n",
      "Replicability at alpha = 0.005 : 99.0 %\n",
      "Replicability at alpha = 0.001 : 96.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PMAT24_A_CR 50\n",
      "0.06943463206291772 0.08859237659473558\n",
      "Replicability at alpha = 0.05 : 10.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PMAT24_A_CR 100\n",
      "0.08304969831389532 0.08786928956024742\n",
      "Replicability at alpha = 0.05 : 4.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PMAT24_A_CR 200\n",
      "0.12186914921832041 0.12159547700742503\n",
      "Replicability at alpha = 0.05 : 55.10204081632652 %\n",
      "Replicability at alpha = 0.01 : 12.244897959183673 %\n",
      "Replicability at alpha = 0.005 : 6.122448979591836 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PMAT24_A_CR 300\n",
      "0.14642306330817234 0.14732241003165472\n",
      "Replicability at alpha = 0.05 : 82.89473684210526 %\n",
      "Replicability at alpha = 0.01 : 31.57894736842105 %\n",
      "Replicability at alpha = 0.005 : 19.736842105263158 %\n",
      "Replicability at alpha = 0.001 : 5.263157894736842 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PMAT24_A_CR max\n",
      "0.18903810660571366 0.18340724009145545\n",
      "Replicability at alpha = 0.05 : 100.0 %\n",
      "Replicability at alpha = 0.01 : 98.0 %\n",
      "Replicability at alpha = 0.005 : 92.0 %\n",
      "Replicability at alpha = 0.001 : 72.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge Flanker_AgeAdj 50\n",
      "-0.016644932932071994 0.020506673063754265\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge Flanker_AgeAdj 100\n",
      "0.010556253902243375 0.013988843177383383\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge Flanker_AgeAdj 200\n",
      "0.02507157553877913 0.024872762330759916\n",
      "Replicability at alpha = 0.05 : 7.142857142857142 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge Flanker_AgeAdj 300\n",
      "0.04174560398029875 0.03948750867349983\n",
      "Replicability at alpha = 0.05 : 15.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge Flanker_AgeAdj max\n",
      "0.05566959147271744 0.0537568231850621\n",
      "Replicability at alpha = 0.05 : 3.4482758620689653 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CardSort_AgeAdj 50\n",
      "0.00825538705753114 0.020822821436505334\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CardSort_AgeAdj 100\n",
      "0.044479176202799627 0.04046146050797244\n",
      "Replicability at alpha = 0.05 : 8.333333333333332 %\n",
      "Replicability at alpha = 0.01 : 8.333333333333332 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CardSort_AgeAdj 200\n",
      "0.05030621786868544 0.043731296512729625\n",
      "Replicability at alpha = 0.05 : 4.545454545454546 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CardSort_AgeAdj 300\n",
      "0.06409981929476123 0.05981810878552007\n",
      "Replicability at alpha = 0.05 : 33.33333333333333 %\n",
      "Replicability at alpha = 0.01 : 4.166666666666666 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge CardSort_AgeAdj max\n",
      "0.09361181906608154 0.08859554355382677\n",
      "Replicability at alpha = 0.05 : 58.06451612903226 %\n",
      "Replicability at alpha = 0.01 : 11.29032258064516 %\n",
      "Replicability at alpha = 0.005 : 6.451612903225806 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PicSeq_AgeAdj 50\n",
      "0.003707916905543548 0.034294688035293515\n",
      "Replicability at alpha = 0.05 : 10.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PicSeq_AgeAdj 100\n",
      "0.045721320747586766 0.05414014050479514\n",
      "Replicability at alpha = 0.05 : 13.636363636363635 %\n",
      "Replicability at alpha = 0.01 : 4.545454545454546 %\n",
      "Replicability at alpha = 0.005 : 4.545454545454546 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PicSeq_AgeAdj 200\n",
      "0.06794095886658241 0.059806195071102956\n",
      "Replicability at alpha = 0.05 : 13.636363636363635 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PicSeq_AgeAdj 300\n",
      "0.08568195230918263 0.07890981970825868\n",
      "Replicability at alpha = 0.05 : 37.5 %\n",
      "Replicability at alpha = 0.01 : 5.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT ridge PicSeq_AgeAdj max\n",
      "0.1037082766166224 0.11518787327130221\n",
      "Replicability at alpha = 0.05 : 93.15068493150685 %\n",
      "Replicability at alpha = 0.01 : 34.24657534246575 %\n",
      "Replicability at alpha = 0.005 : 19.17808219178082 %\n",
      "Replicability at alpha = 0.001 : 2.73972602739726 %\n",
      "CPU times: user 7.79 s, sys: 922 ms, total: 8.71 s\n",
      "Wall time: 3min 19s\n"
     ]
    },
    {
     "data": {
      "text/plain": "     connectivity  model         target    n r_discovery_cv  \\\n0           FS_CT  ridge            age   50       0.218019   \n1           FS_CT  ridge            age   50            NaN   \n2           FS_CT  ridge            age   50            NaN   \n3           FS_CT  ridge            age   50        0.03394   \n4           FS_CT  ridge            age   50       0.153336   \n...           ...    ...            ...  ...            ...   \n2995        FS_CT  ridge  PicSeq_AgeAdj  556       0.206691   \n2996        FS_CT  ridge  PicSeq_AgeAdj  556       0.168531   \n2997        FS_CT  ridge  PicSeq_AgeAdj  556       0.038784   \n2998        FS_CT  ridge  PicSeq_AgeAdj  556       0.113999   \n2999        FS_CT  ridge  PicSeq_AgeAdj  556       0.049988   \n\n     r_discovery_overfit r_replication p_discovery_cv p_discovery_overfit  \\\n0               0.861656      0.176967       0.193806            0.000999   \n1                0.73807      0.016571       0.172827            0.000999   \n2               0.784162      0.073507       0.995005            0.000999   \n3               0.831648     -0.006591       0.682318            0.000999   \n4               0.823246      0.126225       0.406593            0.000999   \n...                  ...           ...            ...                 ...   \n2995            0.440568      0.073023       0.000999            0.000999   \n2996            0.408207      0.088871       0.000999            0.000999   \n2997            0.367894      0.098458       0.212787            0.000999   \n2998            0.397635      0.089869       0.004995            0.000999   \n2999            0.373858      0.061371       0.183816            0.000999   \n\n     p_replication  \n0         0.125874  \n1         0.492507  \n2         0.318681  \n3         0.511489  \n4           0.1998  \n...            ...  \n2995      0.041958  \n2996      0.021978  \n2997      0.004995  \n2998      0.021978  \n2999      0.074925  \n\n[3000 rows x 10 columns]",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>connectivity</th>\n      <th>model</th>\n      <th>target</th>\n      <th>n</th>\n      <th>r_discovery_cv</th>\n      <th>r_discovery_overfit</th>\n      <th>r_replication</th>\n      <th>p_discovery_cv</th>\n      <th>p_discovery_overfit</th>\n      <th>p_replication</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>age</td>\n      <td>50</td>\n      <td>0.218019</td>\n      <td>0.861656</td>\n      <td>0.176967</td>\n      <td>0.193806</td>\n      <td>0.000999</td>\n      <td>0.125874</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>age</td>\n      <td>50</td>\n      <td>NaN</td>\n      <td>0.73807</td>\n      <td>0.016571</td>\n      <td>0.172827</td>\n      <td>0.000999</td>\n      <td>0.492507</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>age</td>\n      <td>50</td>\n      <td>NaN</td>\n      <td>0.784162</td>\n      <td>0.073507</td>\n      <td>0.995005</td>\n      <td>0.000999</td>\n      <td>0.318681</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>age</td>\n      <td>50</td>\n      <td>0.03394</td>\n      <td>0.831648</td>\n      <td>-0.006591</td>\n      <td>0.682318</td>\n      <td>0.000999</td>\n      <td>0.511489</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>age</td>\n      <td>50</td>\n      <td>0.153336</td>\n      <td>0.823246</td>\n      <td>0.126225</td>\n      <td>0.406593</td>\n      <td>0.000999</td>\n      <td>0.1998</td>\n    </tr>\n    <tr>\n      <th>...</th>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n    </tr>\n    <tr>\n      <th>2995</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.206691</td>\n      <td>0.440568</td>\n      <td>0.073023</td>\n      <td>0.000999</td>\n      <td>0.000999</td>\n      <td>0.041958</td>\n    </tr>\n    <tr>\n      <th>2996</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.168531</td>\n      <td>0.408207</td>\n      <td>0.088871</td>\n      <td>0.000999</td>\n      <td>0.000999</td>\n      <td>0.021978</td>\n    </tr>\n    <tr>\n      <th>2997</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.038784</td>\n      <td>0.367894</td>\n      <td>0.098458</td>\n      <td>0.212787</td>\n      <td>0.000999</td>\n      <td>0.004995</td>\n    </tr>\n    <tr>\n      <th>2998</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.113999</td>\n      <td>0.397635</td>\n      <td>0.089869</td>\n      <td>0.004995</td>\n      <td>0.000999</td>\n      <td>0.021978</td>\n    </tr>\n    <tr>\n      <th>2999</th>\n      <td>FS_CT</td>\n      <td>ridge</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.049988</td>\n      <td>0.373858</td>\n      <td>0.061371</td>\n      <td>0.183816</td>\n      <td>0.000999</td>\n      <td>0.074925</td>\n    </tr>\n  </tbody>\n</table>\n<p>3000 rows × 10 columns</p>\n</div>"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "random_state = 42\n",
    "n_bootstrap = 100\n",
    "\n",
    "features = {\n",
    "    'FS_CT': \"FS_*_Thck\"\n",
    "}\n",
    "\n",
    "models = {\n",
    "    'ridge': Ridge()\n",
    "}\n",
    "\n",
    "# We aggregate all results here:\n",
    "df = pd.DataFrame(columns=['connectivity','model','target','n','r_discovery_cv','r_discovery_overfit','r_replication','p_discovery_cv','p_discovery_overfit','p_replication'])\n",
    "\n",
    "for feature_set in features:\n",
    "    for model in models:\n",
    "        for target_var in ['age', 'CogTotalComp_AgeAdj', 'PMAT24_A_CR', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'PicSeq_AgeAdj']:\n",
    "            for sample_size in [50, 100, 200, 300, 'max']:\n",
    "\n",
    "                print('*****************************************************************')\n",
    "                print(feature_set, model, target_var, sample_size)\n",
    "\n",
    "                X, y = create_data(target=target_var, feature_data=features[feature_set])\n",
    "\n",
    "                if sample_size=='max':\n",
    "                    sample_size = int(len(y)/2)\n",
    "\n",
    "                # create random seeds for each bootstrap iteration for reproducibility\n",
    "                rng = np.random.default_rng(random_state)\n",
    "                random_sates = rng.integers(np.iinfo(np.int32).max, size=n_bootstrap)\n",
    "\n",
    "                # run bootstrap iterations in parallel\n",
    "                r_discovery_cv, r_discovery_overfit, r_replication, p_discovery_cv, p_discovery_overfit, p_replication = zip(\n",
    "                    *Parallel(n_jobs=-1)(\n",
    "                    delayed(bootstrap_workhorse)(X, y, sample_size, models[model], seed) for seed in random_sates))\n",
    "\n",
    "                tmp_data_frame = pd.DataFrame({\n",
    "                    'connectivity' : feature_set,\n",
    "                    'model' : model,\n",
    "                    'target' : target_var,\n",
    "                    'n' : sample_size,\n",
    "                    'r_discovery_cv': r_discovery_cv,\n",
    "                    'r_discovery_overfit': r_discovery_overfit,\n",
    "                    'r_replication': r_replication,\n",
    "                    'p_discovery_cv': p_discovery_cv,\n",
    "                    'p_discovery_overfit': p_discovery_overfit,\n",
    "                    'p_replication': p_replication\n",
    "                })\n",
    "                print(tmp_data_frame.r_discovery_cv.mean(), tmp_data_frame.r_replication.mean())\n",
    "\n",
    "                for alpha in [0.05, 0.01, 0.005, 0.001]:\n",
    "                    print('Replicability at alpha =', alpha, ':',\n",
    "                          (tmp_data_frame.loc[tmp_data_frame['p_discovery_cv']<alpha,'p_replication']<alpha).sum() / (tmp_data_frame['p_discovery_cv']<0.05).sum() * 100, '%')\n",
    "\n",
    "                df = pd.concat((df, tmp_data_frame))\n",
    "                df.reset_index(drop=True, inplace=True)\n",
    "                df.to_csv('res/results_thck_Ridge.csv')\n",
    "\n",
    "df\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Null scenario with random target\n",
    "To evaluate false positives with biased estimates"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 50\n",
      "-0.006672762525657837 -0.007750463437365226\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 100\n",
      "-0.0054298235725445355 0.007488203892819826\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 200\n",
      "0.0006834925110597672 0.003009677488480728\n",
      "Replicability at alpha = 0.05 : 20.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age 300\n",
      "-0.0021544089614625975 -0.006674668177288535\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR age max\n",
      "-3.0507795103961733e-05 0.0034802531192872064\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 50\n",
      "-0.014336736049140757 0.018508549683804842\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 100\n",
      "-0.003689805439261738 -0.006670388528765944\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 200\n",
      "-0.011800723070629838 -0.014588490903023055\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj 300\n",
      "-0.021860033767778137 -0.01539863759125782\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CogTotalComp_AgeAdj max\n",
      "0.0024523435006187032 0.003429029401127996\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 50\n",
      "0.016860724627237315 0.0006253730701868871\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 100\n",
      "-0.00017436485006250778 -0.0010680882219589883\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 200\n",
      "0.009168087766612858 -0.0023315245712418327\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR 300\n",
      "-0.0006806813557972834 0.0035995761509331602\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PMAT24_A_CR max\n",
      "0.01329414262459903 -0.0016541675887437185\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 50\n",
      "0.01705340700162105 0.009926042029562062\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 100\n",
      "-0.01798737178052155 0.009547819778861162\n",
      "Replicability at alpha = 0.05 : nan %\n",
      "Replicability at alpha = 0.01 : nan %\n",
      "Replicability at alpha = 0.005 : nan %\n",
      "Replicability at alpha = 0.001 : nan %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 200\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:60: RuntimeWarning: invalid value encountered in long_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.0019681493786523753 0.01594666888263921\n",
      "Replicability at alpha = 0.05 : 66.66666666666666 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj 300\n",
      "0.0021855114322033033 0.0020920437201013766\n",
      "Replicability at alpha = 0.05 : 40.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR Flanker_AgeAdj max\n",
      "-0.002452461838318572 0.0021759495908249722\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 50\n",
      "0.010612277564530275 0.02340110949596856\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 100\n",
      "0.009198801492210631 -0.010635629438048064\n",
      "Replicability at alpha = 0.05 : 33.33333333333333 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 200\n",
      "-0.0007690192235844262 0.0005875610226140504\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj 300\n",
      "-0.008935067626667717 -0.0005087266895158477\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR CardSort_AgeAdj max\n",
      "-0.004974721352014822 -0.005315130990598048\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 50\n",
      "-0.03796633861357167 0.008561640161755903\n",
      "Replicability at alpha = 0.05 : nan %\n",
      "Replicability at alpha = 0.01 : nan %\n",
      "Replicability at alpha = 0.005 : nan %\n",
      "Replicability at alpha = 0.001 : nan %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 100\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:60: RuntimeWarning: invalid value encountered in long_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.010951429727356155 -0.007800783898997811\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 200\n",
      "0.003510008843708663 -0.0054599636733647415\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj 300\n",
      "0.008344196030507038 0.0011190906096572842\n",
      "Replicability at alpha = 0.05 : 33.33333333333333 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "*****************************************************************\n",
      "FS_CT PCA_SVR PicSeq_AgeAdj max\n",
      "-0.004087055146693218 -0.0036189329291476465\n",
      "Replicability at alpha = 0.05 : 0.0 %\n",
      "Replicability at alpha = 0.01 : 0.0 %\n",
      "Replicability at alpha = 0.005 : 0.0 %\n",
      "Replicability at alpha = 0.001 : 0.0 %\n",
      "CPU times: user 7.22 s, sys: 744 ms, total: 7.96 s\n",
      "Wall time: 3min 8s\n"
     ]
    },
    {
     "data": {
      "text/plain": "     connectivity    model         target    n r_discovery_cv  \\\n0           FS_CT  PCA_SVR            age   50            NaN   \n1           FS_CT  PCA_SVR            age   50       0.059573   \n2           FS_CT  PCA_SVR            age   50      -0.057247   \n3           FS_CT  PCA_SVR            age   50      -0.425656   \n4           FS_CT  PCA_SVR            age   50       0.406026   \n...           ...      ...            ...  ...            ...   \n2995        FS_CT  PCA_SVR  PicSeq_AgeAdj  556      -0.059353   \n2996        FS_CT  PCA_SVR  PicSeq_AgeAdj  556       -0.01275   \n2997        FS_CT  PCA_SVR  PicSeq_AgeAdj  556       0.023795   \n2998        FS_CT  PCA_SVR  PicSeq_AgeAdj  556      -0.056596   \n2999        FS_CT  PCA_SVR  PicSeq_AgeAdj  556       -0.07424   \n\n     r_discovery_overfit r_replication p_discovery_cv p_discovery_overfit  \\\n0               0.546397     -0.173968       0.925075            0.000999   \n1               0.512302      0.148245       0.523477            0.000999   \n2               0.287055     -0.117199       0.976024            0.022977   \n3               0.554909      0.021776       0.999001            0.000999   \n4               0.623122     -0.073094       0.064935            0.000999   \n...                  ...           ...            ...                 ...   \n2995            0.174407     -0.014806       0.992008            0.000999   \n2996            0.226049     -0.064682       0.921079            0.000999   \n2997            0.193551     -0.064718       0.436563            0.000999   \n2998            0.195581      0.066811       0.962038            0.000999   \n2999            0.144046      0.039266       0.987013            0.000999   \n\n     p_replication  \n0         0.906094  \n1         0.146853  \n2         0.791209  \n3          0.45954  \n4         0.695305  \n...            ...  \n2995      0.631369  \n2996      0.929071  \n2997      0.943057  \n2998      0.054945  \n2999      0.178821  \n\n[3000 rows x 10 columns]",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>connectivity</th>\n      <th>model</th>\n      <th>target</th>\n      <th>n</th>\n      <th>r_discovery_cv</th>\n      <th>r_discovery_overfit</th>\n      <th>r_replication</th>\n      <th>p_discovery_cv</th>\n      <th>p_discovery_overfit</th>\n      <th>p_replication</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>NaN</td>\n      <td>0.546397</td>\n      <td>-0.173968</td>\n      <td>0.925075</td>\n      <td>0.000999</td>\n      <td>0.906094</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>0.059573</td>\n      <td>0.512302</td>\n      <td>0.148245</td>\n      <td>0.523477</td>\n      <td>0.000999</td>\n      <td>0.146853</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>-0.057247</td>\n      <td>0.287055</td>\n      <td>-0.117199</td>\n      <td>0.976024</td>\n      <td>0.022977</td>\n      <td>0.791209</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>-0.425656</td>\n      <td>0.554909</td>\n      <td>0.021776</td>\n      <td>0.999001</td>\n      <td>0.000999</td>\n      <td>0.45954</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>age</td>\n      <td>50</td>\n      <td>0.406026</td>\n      <td>0.623122</td>\n      <td>-0.073094</td>\n      <td>0.064935</td>\n      <td>0.000999</td>\n      <td>0.695305</td>\n    </tr>\n    <tr>\n      <th>...</th>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n    </tr>\n    <tr>\n      <th>2995</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>-0.059353</td>\n      <td>0.174407</td>\n      <td>-0.014806</td>\n      <td>0.992008</td>\n      <td>0.000999</td>\n      <td>0.631369</td>\n    </tr>\n    <tr>\n      <th>2996</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>-0.01275</td>\n      <td>0.226049</td>\n      <td>-0.064682</td>\n      <td>0.921079</td>\n      <td>0.000999</td>\n      <td>0.929071</td>\n    </tr>\n    <tr>\n      <th>2997</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>0.023795</td>\n      <td>0.193551</td>\n      <td>-0.064718</td>\n      <td>0.436563</td>\n      <td>0.000999</td>\n      <td>0.943057</td>\n    </tr>\n    <tr>\n      <th>2998</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>-0.056596</td>\n      <td>0.195581</td>\n      <td>0.066811</td>\n      <td>0.962038</td>\n      <td>0.000999</td>\n      <td>0.054945</td>\n    </tr>\n    <tr>\n      <th>2999</th>\n      <td>FS_CT</td>\n      <td>PCA_SVR</td>\n      <td>PicSeq_AgeAdj</td>\n      <td>556</td>\n      <td>-0.07424</td>\n      <td>0.144046</td>\n      <td>0.039266</td>\n      <td>0.987013</td>\n      <td>0.000999</td>\n      <td>0.178821</td>\n    </tr>\n  </tbody>\n</table>\n<p>3000 rows × 10 columns</p>\n</div>"
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "random_state = 42\n",
    "n_bootstrap = 100\n",
    "\n",
    "features = {\n",
    "    'FS_CT': \"FS_*_Thck\"\n",
    "}\n",
    "\n",
    "models = {\n",
    "    'PCA_SVR': Pipeline([('pca', PCA(n_components=0.5)),\n",
    "                         ('svr', SVR())])\n",
    "\n",
    "}\n",
    "\n",
    "# We aggregate all results here:\n",
    "df = pd.DataFrame(columns=['connectivity','model','target','n','r_discovery_cv','r_discovery_overfit','r_replication','p_discovery_cv','p_discovery_overfit','p_replication'])\n",
    "\n",
    "for feature_set in features:\n",
    "    for model in models:\n",
    "        for target_var in ['age', 'CogTotalComp_AgeAdj', 'PMAT24_A_CR', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'PicSeq_AgeAdj']:\n",
    "            for sample_size in [50, 100, 200, 300, 'max']:\n",
    "\n",
    "                print('*****************************************************************')\n",
    "                print(feature_set, model, target_var, sample_size)\n",
    "\n",
    "                X, y = create_data(target=target_var, feature_data=features[feature_set]) # gives a random y when target is None\n",
    "\n",
    "                if sample_size=='max':\n",
    "                    sample_size = int(len(y)/2)\n",
    "\n",
    "                # create random seeds for each bootstrap iteration for reproducibility\n",
    "                rng = np.random.default_rng(random_state)\n",
    "                random_sates = rng.integers(np.iinfo(np.int32).max, size=n_bootstrap)\n",
    "\n",
    "                # run bootstrap iterations in parallel, with shuffle_y=True\n",
    "                r_discovery_cv, r_discovery_overfit, r_replication, p_discovery_cv, p_discovery_overfit, p_replication = zip(\n",
    "                    *Parallel(n_jobs=-1)(\n",
    "                    delayed(bootstrap_workhorse)(X, y, sample_size, models[model], seed, shuffle_y=True) for seed in random_sates))\n",
    "\n",
    "                tmp_data_frame = pd.DataFrame({\n",
    "                    'connectivity' : feature_set,\n",
    "                    'model' : model,\n",
    "                    'target' : target_var,\n",
    "                    'n' : sample_size,\n",
    "                    'r_discovery_cv': r_discovery_cv,\n",
    "                    'r_discovery_overfit': r_discovery_overfit,\n",
    "                    'r_replication': r_replication,\n",
    "                    'p_discovery_cv': p_discovery_cv,\n",
    "                    'p_discovery_overfit': p_discovery_overfit,\n",
    "                    'p_replication': p_replication\n",
    "                })\n",
    "\n",
    "                #sns.scatterplot(x='r_replication', y='r_discovery_cv', data=tmp_data_frame)\n",
    "                #plt.ylabel('in-sample (r)')\n",
    "                #plt.xlabel('out-of-sample (r_pred)')\n",
    "                #plt.show()\n",
    "                print(tmp_data_frame.r_discovery_cv.mean(), tmp_data_frame.r_replication.mean())\n",
    "\n",
    "                for alpha in [0.05, 0.01, 0.005, 0.001]:\n",
    "                    print('Replicability at alpha =', alpha, ':',\n",
    "                          (tmp_data_frame.loc[tmp_data_frame['p_discovery_cv']<alpha,'p_replication']<alpha).sum() / (tmp_data_frame['p_discovery_cv']<0.05).sum() * 100, '%')\n",
    "\n",
    "                df = pd.concat((df, tmp_data_frame))\n",
    "                df.reset_index(drop=True, inplace=True)\n",
    "                df.to_csv('res/results_thck_null_PCA_SVR.csv')\n",
    "\n",
    "df"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### *See the notebook called 'plot_results_CT.ipynb' for the results.*"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}