DSP_IHC_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import scipy.stats as ss\n",
    "# import cufflinks as cf\n",
    "import pickle\n",
    "#import pystan\n",
    "#import arviz as az\n",
    "# cf.go_offline()\n",
    "import os\n",
    "from itertools import combinations\n",
    "#%reload_ext blackcellmagic\n",
    "from joblib import Parallel, delayed\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV, RidgeClassifierCV\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.base import TransformerMixin, BaseEstimator\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import (\n",
    "    GroupKFold,\n",
    "    KFold,\n",
    "    cross_val_score,\n",
    "    cross_validate,\n",
    "    LeaveOneOut,\n",
    "    StratifiedKFold,\n",
    "    RepeatedStratifiedKFold,\n",
    "    permutation_test_score\n",
    ")\n",
    "from sklearn.metrics import (\n",
    "    accuracy_score,\n",
    "    make_scorer,\n",
    "    log_loss,\n",
    "    brier_score_loss,\n",
    "    roc_auc_score,\n",
    "    average_precision_score,\n",
    "    roc_curve,\n",
    "    precision_recall_curve\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set random state\n",
    "random_state = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "nanostring_dat = \"cd45_R_data.txt\"\n",
    "\n",
    "raw_data = pd.read_csv(nanostring_dat, sep=\"\\t\", index_col=0)\n",
    "data=raw_data\n",
    "\n",
    "markers = [\"DSP_CD45R\",\"qpath_CD45R\"]\n",
    "\n",
    "data[markers] = StandardScaler().fit_transform(data[markers])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>DSP_CD45R</th>\n",
       "      <th>qpath_CD45R</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timepoint</th>\n",
       "      <th>R</th>\n",
       "      <th>R</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>patient</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-0.346514</td>\n",
       "      <td>-0.838139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>-0.756862</td>\n",
       "      <td>-0.796544</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>1.284794</td>\n",
       "      <td>0.997820</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>0.421587</td>\n",
       "      <td>1.266739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>-0.050441</td>\n",
       "      <td>-0.659771</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          DSP_CD45R qpath_CD45R\n",
       "timepoint         R           R\n",
       "patient                        \n",
       "6         -0.346514   -0.838139\n",
       "22        -0.756862   -0.796544\n",
       "30         1.284794    0.997820\n",
       "36         0.421587    1.266739\n",
       "44        -0.050441   -0.659771"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Reshape to wide format\n",
    "data = data.reset_index().pivot_table(\n",
    "    index=\"patient\", values=markers, columns=\"timepoint\"\n",
    ")\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>DSP_CD45R</th>\n",
       "      <th>qpath_CD45R</th>\n",
       "      <th>pCR</th>\n",
       "      <th>ER</th>\n",
       "      <th>pam50</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timepoint</th>\n",
       "      <th>R</th>\n",
       "      <th>R</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>patient</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>6</th>\n",
       "      <td>-0.346514</td>\n",
       "      <td>-0.838139</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>-0.756862</td>\n",
       "      <td>-0.796544</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>1.284794</td>\n",
       "      <td>0.997820</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>0.421587</td>\n",
       "      <td>1.266739</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>-0.050441</td>\n",
       "      <td>-0.659771</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          DSP_CD45R qpath_CD45R pCR ER pam50\n",
       "timepoint         R           R             \n",
       "patient                                     \n",
       "6         -0.346514   -0.838139   0  1     1\n",
       "22        -0.756862   -0.796544   0  1     0\n",
       "30         1.284794    0.997820   1  0     0\n",
       "36         0.421587    1.266739   0  0     1\n",
       "44        -0.050441   -0.659771   0  1     1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pCRs = raw_data.groupby(\"patient\")[\"pCR\"].first()\n",
    "pCRs = pCRs.sort_index()\n",
    "data[\"pCR\"]=pCRs\n",
    "\n",
    "ERs = raw_data.groupby(\"patient\")[\"ER\"].first()\n",
    "ERs = ERs.sort_index()\n",
    "data[\"ER\"]=ERs\n",
    "\n",
    "pam50s = raw_data.groupby(\"patient\")[\"pam50\"].first()\n",
    "pam50s = pam50s.sort_index()\n",
    "data[\"pam50\"]=pam50s\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt('norm_data.txt', data.values, fmt='%f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# log loss scorer for leave one out - must set the possible label lavels\n",
    "log_loss_with_labels = make_scorer(lambda y_true, y_pred: log_loss(y_true, y_pred, labels=(0, 1)), needs_proba=True, greater_is_better=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FeatureSelector(TransformerMixin, BaseEstimator):\n",
    "    \"\"\"\n",
    "    Subclass to select features for introduction into pipeline.\n",
    "    Setup to answer questions above \n",
    "    Select features based on timepoints, change in timepoints, and different markers\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, marker_names, timepoints):\n",
    "        \"\"\"\n",
    "        :param list names: one or more marker names\n",
    "        :param list timepoints: one or more of [\"B\", \"R\", \"RminusB\"]\n",
    "        \"\"\"\n",
    "\n",
    "        super().__init__()\n",
    "\n",
    "        self.marker_names = marker_names\n",
    "        self.timepoints = timepoints\n",
    "\n",
    "    def transform(self, X):\n",
    "        \"\"\"\n",
    "        X is assumed to have multilevel columns, marker names at level 0 and marker timepoints at level 1\n",
    "        \"\"\"\n",
    "        return X.loc[:, pd.IndexSlice[self.marker_names, self.timepoints]]\n",
    "\n",
    "    def fit(self, X, y=None):\n",
    "        return self\n",
    "\n",
    "    def fit_transform(self, X, y=None):\n",
    "        return X.loc[:, pd.IndexSlice[self.marker_names, self.timepoints]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The folds to use in the inner CV loop used by LogisticRegressionCV\n",
    "kf = StratifiedKFold(n_splits=5, random_state=100, shuffle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_cross_validation(dataframe, features_to_test, scoring_names):\n",
    "    \"\"\"\n",
    "    Run cross validation and return scores\n",
    "    :param pd.DataFrame dataframe: a wide dataframe, with single patient per row, multiindex columns w/ level 0 marker, level 1 timepoint. also must contain \"pCR\" as a column\n",
    "    :param dict features_to_test: dictionary with label as keys, and tuple as values. First index of tuple is markers, Second index of tuple is timepoint\n",
    "    :param list scoring_names: list of sklearn metrics to use for scoring\n",
    "    :return pd.DataFrame scores: dataframe of scores\n",
    "    \"\"\"\n",
    "    \n",
    "    pipelines_to_compare = dict(\n",
    "        [\n",
    "            (\n",
    "                k,\n",
    "                make_pipeline(\n",
    "                    FeatureSelector(*v),\n",
    "                    StandardScaler(),\n",
    "                    LogisticRegressionCV(\n",
    "                        Cs=100,\n",
    "                        penalty=\"l2\",\n",
    "                        solver=\"liblinear\",\n",
    "                        class_weight=\"balanced\",\n",
    "                        scoring=\"accuracy\",\n",
    "                        random_state=random_state,\n",
    "                        cv=kf,\n",
    "                    )\n",
    "                ),\n",
    "            )\n",
    "            for k, v in features_to_test.items()\n",
    "        ]\n",
    "    )\n",
    "    ## We'll do 100x repeated max-fold stratified cross-validation, then unpaired t-testing for significance\n",
    "    scores = {}\n",
    "\n",
    "    for feature, pipeline in pipelines_to_compare.items():\n",
    "        scores[feature] = pd.DataFrame.from_records(\n",
    "            cross_validate(\n",
    "                pipeline,\n",
    "                dataframe,\n",
    "                dataframe.pCR,\n",
    "                cv=RepeatedStratifiedKFold(random_state=random_state,\n",
    "                    n_splits=dataframe.pCR.value_counts().min()//2, # max-folds for stratified is the minimum of number of samples in class\n",
    "                    n_repeats=100 # TODO: since small dataset, probably should check what the max possible combinations are\n",
    "                ),\n",
    "                scoring=scoring_names,\n",
    "                n_jobs=36,\n",
    "            )\n",
    "        )\n",
    "    scores = pd.concat(scores.values(), keys=scores.keys(), names=[\"marker_group\", \"fold\"])\n",
    "    \n",
    "    return scores,pipelines_to_compare\n",
    "\n",
    "\n",
    "def get_model_weights(dataframe, model_pipeline, features_to_use):\n",
    "    \"\"\"\n",
    "    Train a model and then return model weights\n",
    "    :param pd.DataFrame dataframe:  a wide dataframe, with single patient per row, multiindex columns w/ level 0 marker, level 1 timepoint. also must contain \"pCR\" as a column\n",
    "    \"\"\"\n",
    "    model_pipeline.fit(dataframe, dataframe.pCR)\n",
    "    \n",
    "    weights = pd.DataFrame(model_pipeline[-1].coef_.T, index=dataframe.loc[:, pd.IndexSlice[features_to_use[0], features_to_use[1]]].columns, columns=[\"weight\"])\n",
    "    weights = weights.append(pd.DataFrame([[model_pipeline[-1].intercept_[0]]], index=pd.MultiIndex.from_tuples([(\"intercept\", \"NA\")]), columns=[\"weight\"])).rename_axis([\"marker\", \"timepoint\"])\n",
    "    return weights\n",
    "    \n",
    "\n",
    "def calculate_significance(scores):\n",
    "    \"\"\"\n",
    "    Calculates significance with unpaired t-test and holm-bonferroni correction\n",
    "    :param pd.DataFrame scores: dataframe output by sklearn cross_validate function\n",
    "    \"\"\"\n",
    "    \n",
    "    # significance testing\n",
    "    cc = [list(a) for a in combinations(scores.index.levels[0],2)]\n",
    "    metric_p_values = {}\n",
    "    for metric in scores.columns:\n",
    "        p_values = {}\n",
    "        for combination in cc:\n",
    "            p_values[tuple(combination)] = ss.ttest_ind(scores.loc[combination[0], metric], scores.loc[combination[1], metric])[1]\n",
    "\n",
    "        metric_p_values[metric] = p_values\n",
    "    metric_p_values = pd.DataFrame.from_records(metric_p_values)\n",
    "    significances = pd.concat([\n",
    "        metric_p_values,\n",
    "        metric_p_values.apply(lambda x: sm.stats.multipletests(x, alpha=0.05, method=\"holm\")[1], axis=0)\n",
    "    ], axis=1, keys=[\"p-value\", \"adj_p-value\"])\n",
    "    return significances\n",
    "\n",
    "def do_permutation_test_vs_shuffled_null(dataframe, features_to_use):\n",
    "    \"\"\"\n",
    "    Train a model and compare using permutation test to null\n",
    "    :param pd.DataFrame dataframe:  a wide dataframe, with single patient per row, multiindex columns w/ level 0 marker, level 1 timepoint. also must contain \"pCR\" as a column\n",
    "    :param tuple features_to_use: Tuple of features to use. First index of tuple is markers, Second index of tuple is timepoint\n",
    "    :return The true score without permuting targets, The scores obtained for each permutations, The p-value, which approximates the probability that the score would be obtained by chance. \n",
    "    \"\"\"\n",
    "    \n",
    "    model_pipeline = make_pipeline(\n",
    "                    FeatureSelector(*features_to_use),\n",
    "                    StandardScaler(),\n",
    "                    LogisticRegressionCV(\n",
    "                        Cs=100,\n",
    "                        penalty=\"l2\",\n",
    "                        solver=\"liblinear\",\n",
    "                        class_weight=\"balanced\",\n",
    "                        scoring=\"accuracy\",\n",
    "                        random_state=random_state,\n",
    "                        cv=kf,\n",
    "                    )\n",
    "                )\n",
    "    score, permutation_scores, pvalue = permutation_test_score(model_pipeline, dataframe, dataframe.pCR, scoring=\"roc_auc\", cv=RepeatedStratifiedKFold(random_state=random_state,n_repeats=10,n_splits=dataframe.pCR.value_counts().min()//2), n_permutations=1000, n_jobs=36, )\n",
    "    \n",
    "    return score, permutation_scores, pvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_cross_validation_plots(dataframe, features_to_test):\n",
    "    \"\"\"\n",
    "    Run cross validation and return scores for plotting\n",
    "    :param pd.DataFrame dataframe: a wide dataframe, with single patient per row, multiindex columns w/ level 0 marker, level 1 timepoint. also must contain \"pCR\" as a column\n",
    "    :param dict features_to_test: dictionary with label as keys, and tuple as values. First index of tuple is markers, Second index of tuple is timepoint\n",
    "    :return tuple (pd.DataFrame roc_scores, pd.DataFrame pr_scores)\n",
    "    \"\"\"\n",
    "    pipelines_to_compare = dict(\n",
    "        [\n",
    "            (\n",
    "                k,\n",
    "                make_pipeline(\n",
    "                    FeatureSelector(*v),\n",
    "                    StandardScaler(),\n",
    "                    LogisticRegressionCV(\n",
    "                        Cs=100,\n",
    "                        penalty=\"l2\",\n",
    "                        solver=\"liblinear\",\n",
    "                        class_weight=\"balanced\",\n",
    "                        scoring=\"accuracy\",\n",
    "                        random_state=random_state,\n",
    "                        cv=kf,\n",
    "                    ),\n",
    "                ),\n",
    "            )\n",
    "            for k, v in features_to_test.items()\n",
    "        ]\n",
    "    )\n",
    "    ## We'll do 100x repeated max-fold stratified cross-validation, then unpaired t-testing for significance\n",
    "    roc_scores = {}  # used to store fpr and tpr\n",
    "    pr_scores = {}\n",
    "\n",
    "    def train_and_score_model(pipeline, train_X, test_X, train_y, test_y):\n",
    "        pipeline.fit(train_X, train_y)\n",
    "        preds = pipeline.decision_function(test_X)\n",
    "        fpr, tpr, roc_threshs = roc_curve(test_y, preds)\n",
    "        precision, recall, pr_threshs = precision_recall_curve(test_y, preds)\n",
    "        return (\n",
    "            pd.DataFrame([fpr, tpr], index=[\"fpr\", \"tpr\"], columns=roc_threshs),\n",
    "            pd.DataFrame(\n",
    "                [recall, precision],\n",
    "                index=[\"recall\", \"precision\"],\n",
    "                columns=list(pr_threshs) + [1000],\n",
    "            ),\n",
    "        )\n",
    "    \n",
    "    for feature, pipeline in pipelines_to_compare.items():\n",
    "        rskf = RepeatedStratifiedKFold(\n",
    "            random_state=random_state,\n",
    "            n_splits=dataframe.pCR.value_counts().min()\n",
    "            // 2,  # max-folds for stratified is the minimum of number of samples in class, but we want more than one positive/negative class in each test fold\n",
    "            n_repeats=100,\n",
    "        )\n",
    "        pr_data = []\n",
    "        data = Parallel(n_jobs=32)(\n",
    "            delayed(train_and_score_model)(\n",
    "                pipeline,\n",
    "                dataframe.iloc[train],\n",
    "                dataframe.iloc[test],\n",
    "                dataframe.pCR.iloc[train],\n",
    "                dataframe.pCR.iloc[test],\n",
    "            )\n",
    "            for train, test in rskf.split(dataframe, dataframe.pCR)\n",
    "        )\n",
    "        roc_data, pr_data = list(zip(*data))\n",
    "        roc_scores[feature] = pd.concat(\n",
    "            roc_data, keys=np.arange(len(roc_data)), names=[\"fold\", \"metric\"]\n",
    "        )\n",
    "        pr_scores[feature] = pd.concat(\n",
    "            pr_data, keys=np.arange(len(pr_data)), names=[\"fold\", \"metric\"]\n",
    "        )\n",
    "    roc_scores = pd.concat(roc_scores, keys=features_to_test, names=[\"feature\"])\n",
    "    pr_scores = pd.concat(pr_scores, keys=features_to_test, names=[\"feature\"])\n",
    "    return roc_scores, pr_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "scoring_names = [\n",
    "    \"neg_log_loss\",\n",
    "    \"accuracy\",\n",
    "    \"brier_score_loss\",\n",
    "    \"roc_auc\",\n",
    "    \"average_precision\",\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. how does DSP protein data (CD45) compare to IHC data (CD45) in terms of predictive ability? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_to_test = {\n",
    "    \"qpath_CD45R\": ([\"qpath_CD45R\"],[\"R\"]),\n",
    "    \"DSP_CD45R\": ([\"DSP_CD45R\"],[\"R\"])\n",
    "}\n",
    "\n",
    "scores,pipelines = run_cross_validation(data, features_to_test, scoring_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"3\" halign=\"left\">test_accuracy</th>\n",
       "      <th colspan=\"3\" halign=\"left\">test_average_precision</th>\n",
       "      <th colspan=\"3\" halign=\"left\">test_brier_score_loss</th>\n",
       "      <th colspan=\"3\" halign=\"left\">test_neg_log_loss</th>\n",
       "      <th colspan=\"3\" halign=\"left\">test_roc_auc</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>50%</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>50%</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>50%</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>50%</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>50%</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>marker_group</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>DSP_CD45R</th>\n",
       "      <td>0.7232</td>\n",
       "      <td>0.168492</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>0.823100</td>\n",
       "      <td>0.202447</td>\n",
       "      <td>0.833333</td>\n",
       "      <td>-0.215919</td>\n",
       "      <td>0.076017</td>\n",
       "      <td>-0.209955</td>\n",
       "      <td>-0.703278</td>\n",
       "      <td>1.024404</td>\n",
       "      <td>-0.611456</td>\n",
       "      <td>0.810667</td>\n",
       "      <td>0.226427</td>\n",
       "      <td>0.875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>qpath_CD45R</th>\n",
       "      <td>0.8856</td>\n",
       "      <td>0.123850</td>\n",
       "      <td>0.833333</td>\n",
       "      <td>0.865967</td>\n",
       "      <td>0.171419</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.242682</td>\n",
       "      <td>0.028705</td>\n",
       "      <td>-0.249897</td>\n",
       "      <td>-0.677861</td>\n",
       "      <td>0.064997</td>\n",
       "      <td>-0.692943</td>\n",
       "      <td>0.871250</td>\n",
       "      <td>0.171917</td>\n",
       "      <td>1.000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             test_accuracy                     test_average_precision  \\\n",
       "                      mean       std       50%                   mean   \n",
       "marker_group                                                            \n",
       "DSP_CD45R           0.7232  0.168492  0.800000               0.823100   \n",
       "qpath_CD45R         0.8856  0.123850  0.833333               0.865967   \n",
       "\n",
       "                                 test_brier_score_loss                      \\\n",
       "                   std       50%                  mean       std       50%   \n",
       "marker_group                                                                 \n",
       "DSP_CD45R     0.202447  0.833333             -0.215919  0.076017 -0.209955   \n",
       "qpath_CD45R   0.171419  1.000000             -0.242682  0.028705 -0.249897   \n",
       "\n",
       "             test_neg_log_loss                     test_roc_auc            \\\n",
       "                          mean       std       50%         mean       std   \n",
       "marker_group                                                                \n",
       "DSP_CD45R            -0.703278  1.024404 -0.611456     0.810667  0.226427   \n",
       "qpath_CD45R          -0.677861  0.064997 -0.692943     0.871250  0.171917   \n",
       "\n",
       "                     \n",
       "                50%  \n",
       "marker_group         \n",
       "DSP_CD45R     0.875  \n",
       "qpath_CD45R   1.000  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# summary statistics\n",
    "summary_scores = scores.groupby(level=0).describe()\n",
    "summary_scores.loc[:, pd.IndexSlice[map(lambda x: \"test_\" + x, scoring_names), [\"mean\", \"std\", \"50%\"]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "roc_scores, pr_scores = run_cross_validation_plots(data, features_to_test)\n",
    "roc_scores = roc_scores.fillna(axis=1, method=\"bfill\").fillna(axis=1, method=\"ffill\")\n",
    "roc_scores = roc_scores.groupby([\"feature\", \"metric\"]).apply(lambda x: pd.DataFrame([x.mean(axis=0), x.quantile(0.40, axis=0), x.quantile(0.66, axis=0)], index=[\"mean\", \"lower\", \"upper\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'ROC')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEjCAYAAAAsbUY2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3iUVfrw8e9JIxASaqihBwyEQKiC0rsiiAiisoq6gK6goqsI9kVsgCKKjfUVlGUBwRX5YUMEBOmo9BZKCIFAAoT0OnPeP04yJCSZTCCTSbk/15Ur85R55k4Ic89zyn2U1hohhBCiIG6uDkAIIUTpJolCCCGEXZIohBBC2CWJQgghhF2SKIQQQtgliUIIIYRdkiiEEELYJYlCCEApFa6USlFKJSqlziulFimlquY4fotSar1SKkEpFaeU+j+lVJtrruGnlHpfKRWRdZ3jWdu1S/4nEqL4SKIQ4qphWuuqQCjQAZgOoJTqDqwFvgMaAM2AvcAWpVTzrHO8gF+BYGAI4AfcAlwCupbsjyFE8VIyM1sIc0cBjNdar8vangUEa62HKqU2A/u11o9f85wfgRit9YNKqfHAG0ALrXViCYcvhFPJHYUQ11BKBQC3AceVUlUwdwYr8jn1a2Bg1uMBwE+SJER5JIlCiKtWKaUSgDNANPAqUBPz/yQqn/OjgOz+h1oFnCNEmSeJQoirRmitfYE+QBAmCcQCVqB+PufXBy5mPb5UwDlClHmSKIS4htb6N2ARMEdrnQRsA0bnc+o9mA5sgHXAYKWUT4kEKUQJkkQhRP7eBwYqpUKBacA4pdSTSilfpVQNpdRMoDvwr6zzF2OarL5RSgUppdyUUrWUUi8opW53zY8gRPGQRCFEPrTWMcBXwMta69+BwcBITD/Eaczw2R5a67Cs89MwHdpHgF+AeGAnpvlqR4n/AEIUIxkeK4QQwi65oxBCCGGXJAohhBB2SaIQQghhlyQKIYQQdkmiEEIIYZckCiGEEHZJohBCCGGXJAohhBB2SaIQQghhlyQKIYQQdkmiEEIIYZckCiGEEHZJohBCCGGXJAohhBB2OS1RKKW+UEpFK6UOFHBcKaU+UEodV0rtU0p1dFYsQgghrp8z7ygWAUPsHL8NaJn1NRH4xImxCCGEuE5OSxRa603AZTun3Al8pY3tQHWllCxOL4QQpYyHC1+7IWaN4WyRWfuirj1RKTURc9eBj49Pp6CgoBIJUAghSiOtwZLP6qTpmVYyrVaS0y2kpFuonBJPjbgYKpHBH2lc1Fr7X8/ruTJRqHz25bsuq9Z6AbAAoHPnznr37t3OjEsIIW6I1ao5fD6e9ExroedeiE9jV/hlTl9KorJX3rfknacu4a4U7u5X3zLPXE6xe82aqYlM2PktE8NXUammBfw9UfsyThf9JzFcmSgigUY5tgOAcy6KRQhRgUXHp5KWaeXM5WRiEtM4czmZ8/GpeLqb1vltJy5RycONvZFxTo2jnp83Vbzcc+2r7OnO5aR0BrSua9vXsbHGp5IHrev55jrXatW0PbiDNp+/TeXI/eChoG9TmPMT1Gl13XG5MlGsBiYrpZYBNwNxWus8zU5CCHG90jItWKxXGypSM6y88+MRLiWl2fb9FXGFS0np+T7ft5IHKMi0aFIyLHRvXosaPp60rOOb7/k5JaVlcmtg7fzbTq5Rs4oXN9XzxdvTvfCT7Vm/Hj59G87vhwbu0KsS/HMF+Le8ocs6LVEopZYCfYDaSqlI4FXAE0Br/SnwA3A7cBxIBh52VixCiIohKS2TyNgU0jItPLVsD6cuJhV4bnADPwDq+nlTycONf/QNxNvDDX/fSjSqWYVaPl5Ur+JVUqHfOEsGpP0K7odMgpjwAvR/qVgu7bREobW+r5DjGpjkrNcXQpR/47/czfojF2yfxJPTLXnOmdw3EF/v3G91993cGD9vzxKJ0akiI+HdWdC3Bpz8H1yJgNu9oVYgdHyg2F7GlU1PQghRJFprjpxPYMfJS/x3ZwTHLiQCMPbmxlnHoZKnG8ENquHuphjYui5ubg60/ZQ1+1bCu1PhtyhIT4e/PKCPtzn29AGo3rhYX04ShRDCpTItVqLiUgE4dTGJKykZBZ67Jewiy3efybXv13/2poV/VafG6HKZ6bD2JVBuEHYU/t/3EGMx27eEwD29YORs8KjklJeXRCGEKHGnLiZx9Hw8B8/F8+H640V+/pt3hdAvqA7+vpVwL493DABWK+xbBnuXQuQfkJoIezJgjwUsFmjaBmZ9Aj16OD0USRRCiGIXn5pB1JVU2/axCwks3n6as7Ep+FRytzUZZfPz9uCVYcForWlQvTJ1/bwLvLaftwd17Bwv9SK2Q8xRUAUkOEs6rJsBadcMxXVvCeeBugpGjYInngAfH6eHC5IohBDFKNNi5fPfT/H2j0cKPKdny9o0rlmFWwNr071FLZrXroqnu0IV9MZZ2l06Yd78d34GidGQEAXuXuCWT2d5RsGjsPLVtD/cMQtqB5rtlkugdWvoWLI1VCVRCCGKLD3Tyq7wy6w9eJ5zcakcjorHYtW2vgaAzk1q8EiPZrbtRjWqEBJQzRXhOofVAqe3wpd35N5fvQm0GgLuBYyqSk+EloOhfruCr+3hDYdOwcyZEHgeemQlirFjiyf2IpJEIYRw2KFz8RyPSeTJpX/l2l/Pz5v41AxGdmhIZS93nht8U9mag5AtPRniIiFqj+kozhb+O5zaBJVrmCajyF25n3fTULjt7eIZbZSQAHPnwapVZnvFihLph7BHEoUQolCpGRY2h11kwle566x9/Wh32gVUu/EZxa5gyTBfl8Lgl1fg5MbCn9OgA3jXgOZ9wZoJDTtCi/7Q5JaC7yCKYtMmeOstiIkBT08YPx7Gjbvx694gSRRCiDziUzM4n9WMtC8yjmdX7LUduzO0Ac8MbEXjmlVKT7+CJQPCfoGkmKudxFrDvuVwMQx8apOrlkZKLCTkU1qu0c3QsDPUa2u+51TV39xROEN8PLz9Nqxda7ZDQuDll6F5c+e8XhFJohBCABCXnMHPB88z9Zt9BZ7z41M9aV3fr+SCSo03o4ASL0Ds6atJIC4Stn8CsadME5E10/51qjUEv4Y5djSDtECo1w58/KFuWwjsX/BIJGdzd4c9e8DbGyZNgjFjwK30rFQtiUKICs5i1bR66cdcxfP63uRPu4DqtKrri0Vr2jWsRtPaThiKabXAhQNgyflGr00n8Y5PIf6s/ed7eEONphA4ALQV2twJ1QKuHnf3gqp1ij/u4nDhAlSrZpKDj4+5o6hdGxo0cHVkeUiiEKICOXYhgQWbTuLn7UnE5WTWHb6Q6/iMO4NpH1Cd9o2qOy+IzHQzhDT5Ivy7X+Hn93kBqtSESn7gn6NUdtV64FcGF8W0WuHbb2HePBg5EqZMMfvb2RkF5WKSKIQow9IyLfwVcQWr1uw9E4c1a9Wznw6c5+C5OKpV9sQ9RxPGxcSr5bW9Pc3+kR0bUtfPm8l9A/GpVExvCVrDptmw4Q2oXBPcclw3KTr3ub71YdgHufcpBQ06gk+t4omntIiIMENe//zTbEdFmcRRipqZ8iOJQohSJC4lg7UHz9ve8A+diyc1w5rv+8jJmCR2nLK3LD10alKTOn5X6/9oDR0bV2d050Z2nnWDTm+Dr+4ES3ZS0hA0NMcJ2iSOBh2hki8E3VHq3yhvmMUCS5bAp5+aIn41a8Lzz0O/fq7rFykCSRRCOFl6ppVPNp4gw5L/spjpFisbjkQTFp2Y73GAOr55i71lJ5Perfz5R58WALSu54e3l3nT9XJ3K9lRSTFH4aObybWi8bNhpbePoKQkJcFjj8Hhw2Z76FB45hnTP1FGSKIowMo/Itl4NLr8FhwTJea7PVeHYeb395TdiVzLx4uGNSozsHVd7u50tUO2po9X6ZynYMmEC/tNCYtt8+Fcjkl4Y5ZA6zsKfm5F4uMD9epBbCy88ALccourIyoySRQF+H+/n+JwVDxNa1VxdSiijGtcswo+lTxY+Vj34usDKAlam2Go2mqGn343GTJTTLOR1nB2d97nDJsHHceVieYUp9q/H6pUgRbmTo8XXwQvL7OvDCpDf7Uly8vDjVta1OK/E7q5OhQhnMdqzdGXkOXKGdj2Ifz5Vf7PaZE1UqlpT/CsAl3Gm9FHtVs5bT2EMiMlBT7+GJYtM8X7Fi40cySqO3EUWQmQRFEABXi4l/MONlExpSXCX4shYhsc+s7+ufXawc2PmceelU3Hs0cZrOFUEnbuNCOazp0znfM332w6sd1LYbNhEUmiEKKiyEg1lUtnt8i9v0EHM1Etm7aaCqht7iye+kXlXXw8vP8+rF5ttlu1gldegaAg18ZVjCRRCFGeaQ3rZ8Kl43BoVY4DCv55BHzruSy0csFigYceMvMjvLxgwgR44AHwKF9vreXrpxFCmH6HrfPMIjqXT8Kxn8x+vwBo1MX0LQTfZWY7ixvj7g733gs//WTuIpo2dXVETiGJQoiyzJIBf35piueBqZt04Jurxz0qm+8PrIIWfUs+vvJGa/jxR3MnMWyY2TdqlPkqx5MGJVEIURZdiTAd0Wtfyv940B1mqKpP7ZKNqzyLijJrRWzdaoa5du9uiviV4wSRTRKFEKXdn4shNhxObzHVUGOOQuL5q8dbDoIRn4BXVbPt5gHu8l+72FitsHIlzJ8Pycng52dmVtcqZ3Wo7JC/JiFKgtYQcwQyknPvt2TAmR2mrHYl39zHoo+YTujMlNz7A7pC5epm/kKrwcWz/KbI3+nT8PrrZq0IgP79YerUCpUkQBKFEDfu5EZIvmQeR+zIXSkVYPtHjl+rZo4VzbTVfN00FPq9CHWDbzhUUUQzZsDevaaI37RppohfBSSJQgh7UmLh/IGr2xePQmqcefzHl3DldP7P88pxd+DmAf6toemt0KxX3kTi5mHmMsgopNJB66slSJ5/3syynjLFNDlVUJIohLBaTRPPldNmxnLM0atv5hcO2H8umHWW+75g1lUAqNlC+gjKovR0+PxziIyEN980+7Inz1Vw8tcsyqeUWJMAYg5DUozpAzi1Kf8JZic35t0XlFX5tEZTqFrXzDvIVqMJ+NQxazVLOYvyYe9e08x0Omtd7oceMklCAJIoRFl2ajOseMiUmXD3Mm/cALGnCn5OfBTUuaa0QkBXsGbALU+Cjz807AReZbPKpyii5GQzmmnFCtPk1LSpuYOQJJGLJApReqTGQ+oV81hbYdnYHE0/15at1rk3g0debS4K6AIpl6HlYLCkQ902UK0x1A50ZvSirNm2Dd54A86fNzOsx42D8eNNKQ6RiyQKUfwy08wbvT3n98OJ9eZOQFtMPaKCBHSB5vnMKk5PhNbDoEnZWwhGlALbtpkkERQkdxGFkEQhHKO1adKJjzJt+lciwNPbHDu+3swBUMp0BFszru81qtY1i97UaHJ1X7t7pWNYFJ/YWKhRwzz+xz8gIADuvrtclAJ3Jqf+D1RKDQHmAe7A51rrt6853hj4Eqiedc40rfUPzoxJXCMzHdIScu+LOQLH18HRH0zn7+lteRe3AfCubhaqsWZC8kWzoE2NppB4AVoNATc7//m01QwZzV4ERzqGhTNdvAizZsGhQ/D116YER+XKcM89ro6sTHBaolBKuQMfAQOBSGCXUmq11vpQjtNeAr7WWn+ilGoD/AA0dVZMZU7UPtj2UeHNONcjbC14Vyt4HkC2+Cio3960+Te6GQL7QyU/aNTVPF+I0kxrWLMG3nsPEhJMgjhyBDp2dHVkZYoz7yi6Ase11icBlFLLgDuBnIlCA9mzWKoB56hoMlLNp/ewn83i9B7e2DpuI3ea7zWaFf/releDzFRoO8qsXFav3dVjlnQzC7jJLbK0pSi7zp0zndU7dpjtW26BF16AerIGR1E5M1E0BM7k2I4Ebr7mnNeAtUqpJwAfYIAT47k+Z3bClnnOufaJ9Xlr/8DVjtvmfc0b9uA3nPP6QpRX338Pb79t1rD284Nnn4Xbbrs641oUiTMTRX7/IteMaeQ+YJHW+l2lVHdgsVKqrda521qUUhOBiQCNG5dwAbQ1T5shmnXbFv+1azQz/QOB/aH7ZKjVQv6QhSgOvr4mSQwaZJJETSmPciOcmSgigUY5tgPI27T0d2AIgNZ6m1LKG6gNROc8SWu9AFgA0Llz52uTjXNVqWUWf/nHlhJ9WSFEEWRmwl9/QZcuZrtXL/jqK2jTxrVxlRPOXHFjF9BSKdVMKeUF3AusvuacCKA/gFKqNeANxDgxpqI79Zv5pC+EKJ2OHIEHH4RJk8yopmySJIqN0+4otNaZSqnJwM+Yoa9faK0PKqVmALu11quBfwL/Vko9jWmWekhrXbJ3DIXxrX91QRghROmRlgYLFsDixaauV4MG5s5CFDunzqPImhPxwzX7Xsnx+BBwqzNjuGHKHWpJ6QchSpW//jILCkVEmH69++83E+gqV3Z1ZOWSTHkVQpQtq1bBzKySL82bw8svQ0iIa2Mq5yRRFCah4k3tEKJUu/VWqF4dRo+Ghx+WIn4lQBKFHZ463cyKTipd/etCVChxcbB8uans6uYG/v6werWZZS1KhCQKO9x1VsdYo66uDUSIikhrWLfO1GiKjQUfHxg71hyTJFGiJFE4QspYCFGyYmLgnXdg40az3bGjmRshXEIShRCi9NDaNCvNnQuJiebOYcoUGDHCNDsJl5BEYYet6UkIUTJ+/dUMewXo0cMU8atTx7UxCUkU9twZv8Q8KGVzAIUot/r1M01MgwebOk1S+6xUkHs5O3R2XcObH3NtIEKUVydPmtIb0Vnl3dzczNoRgwdLkihFJFEUxquqrLwmRHHLyIDPPzczqnfsgE8+cXVEwg5pehJClKxDh2DGDDh+3GyPHAlPPunamIRdkiiEECUjNRU++wyWLDFF/AIC4KWXoHNnV0cmCiGJQghRMiIiTJIAeOABePRR8PZ2bUzCIZIohBDOk5YGlbImrLZqZVabCw42X6LMkM5sIYRz/P473HUX/Pbb1X333CNJogySOwohRPGKjYV334WffjLb//d/0Lu3a2MSN0QShRCieGgNv/xiivhduWKanCZNgnvvdXVk4gZJohBC3LjYWFN6Y9Mms92lixnR1LCha+MSxUIShRDixlWqBGFhULUqPP00DB8uM6vLEUkUQojrc+YM1KplKrxWqWLKgvv7my9RrsioJyFE0VitsHgxjBkDH398dX+bNpIkyimH7yiUUj5Aqtba4sR4hBCl2YkT8K9/mTIcAAkJJnHIWhHlWoGJQinlBtwLjAW6AGlAJaVUDPADsEBrHVYiUQohXCsjAxYuhC++gMxMs0bEiy/Crbe6OjJRAuzdUWwA1gHTgQNaayuAUqom0Bd4Wyn1rdb6P84P0zXctAUsGa4OQwjXSkyERx4xJcEBRo2CJ54wa1iLCsFeohigtc7zLqm1vgx8A3yjlPJ0WmSlQJu0PWBJd3UYQrhW1arQooW5k3jpJbN+tahQCkwU2UlCKTUHWKi1PljQOeWRl06lUUY4VK7h6lCEKHm7dkG1aqY+E8D06WYIbHbdJlGhONKZfQRYoJTyABYCS7XWcc4Ny8WWjObri2vN4zpSl0ZUIAkJMG8erFplksRXX4GHB/j5uToy4UKFDlXQWn+utb4VeBBoCuxTSv1XKdXX2cG5xLk9EGaSxGq/++D+5S4OSIgSsmmTKdq3ahV4esKAAa6OSJQSDg2PVUq5A0FZXxeBvcAzSqlHtdblq5DLAlO87MOqU9hVfSjDK1V1cUBCONnlyzBnDqzNuotu1w5efhmaNXNtXKLUKDRRKKXeA4YB64E3tdY7sw69o5Q66szgStyprDo1NVvwq/sg5GZblHsWCzz8MJw9axYRmjzZ3FXIvAiRgyN/DQeA9lrrR3MkiWxdnRCTa1gt8OUw87j380RcTnZtPEKUBHd3GDcOunaFr782lV4lSYhrOPIXMVZrnetdUyn1K0C56tQ++K353rwvut09XE5K50JcqmtjEqK4Wa3wzTfwv/9d3XfXXfDRR9CggeviEqWavZnZ3kAVoLZSqgaQXQrSDyh/f1GnzCpceuC/eHjRLgCCG0jjkyhHIiJg5kz480/TzNS7tynqJ1VeRSHs9VE8CkzBJIU/c+yPBz5yZlCuoDNTUUCreeFkZP1axt3S1KUxCVEsLBZYsgQ+/RTS06FmTZg2zSQJIRxgb8LdPGCeUuoJrfWHJRiTSySmZBCna5OBB31v8ue9e0Kp4ePl6rCEuDHHjsGMGXDkiNm+4w545hmZFyGKxF7TUz+t9XrgrFJq5LXHtdb/y+dp115jCDAPcAc+11q/nc859wCvARrYq7W+3/Hwi8fl2MvUDPuWZKoz795Q7gyVVblEOaC1WSPiyBGoV88U8eve3dVRiTLIXtNTb8yQ2GH5HNOA3USRNffiI2AgEAnsUkqt1lofynFOS0zRwVu11rFKqTpFjP+Gaa15fs7H/NsTTnm2ZGCbuiUdghDFK7vst1Lwwgum43rSJLO4kBDXwV7T06tZD8df5xoUXYHjWuuTAEqpZcCdwKEc50wAPtJax2a9ZvR1vM4N2Rx2EaspjEvI396iipcs+ifKqORk+OQTOH8eZs0yiaJFC3juOVdHJso4R4bHnlJKLVBK9VeqSMMjGgJncmxHZu3LqRXQSim1RSm1PaupKg+l1ESl1G6l1O6YmJgihFC4+RuO00SZ/OQjSUKUVTt2mDkQS5fCb7+ZBYaEKCaOJIqbMOtSTMIkjflKqR4OPC+/pKKv2fYAWgJ9gPuAz5VS1fM8SesFWuvOWuvO/sW81KKftwfTPP6btSF9E6KMiY83ndWTJsG5c6aQ3+LFEBjo6shEOVLoR2itdQrwNfB11nyKecBvmA5qeyKBRjm2A4Bz+ZyzPatc+amskiAtgV2OhX/jlLbgpSxQuSZULfEuEiGu38aN8NZbcOkSeHnBxInwt7+Zaq9CFCOH5uorpXorpT7GzKfwBu5x4Gm7gJZKqWZKKS/MsqqrrzlnFWa1PJRStTFNUScdjL1Y+Fz4wzxoKks6ijJm3z6TJEJDTZPTQw9JkhBO4UhRwFPAHsxdxXNa6yRHLqy1zlRKTQZ+xtx9fKG1PqiUmgHs1lqvzjo2SCl1CLBkXf/Sdf4sRbb+yAVU3BnwArpOLKmXFeL6aA0xMWa9ajB3EE2bmrkRUp9JOJEjHz/aa63jr+fiWusfgB+u2fdKjscaeCbrq8R98+dZBrvvMRuVa7oiBCEcExUFb75pOqm//tosT+rtDcOHuzoyUQHYm3A3VWs9C3hDKXVtJzRa6yedGlkJ2B52gY/ct5mN2i1dG4wQ+bFaYeVK+PBDSEkxM6pPnjRrRghRQuzdURzO+r67JAIpaRarpnfaRtPs1KADeMhawKKUOX0aXn8d9mTd9fbvD88/b2o1CVGC7E24+7+sh8la6xU5jymlRjs1qhKQYbHip7Kqp9+z2LXBCHGt//3PrDqXnm6K9z3/PPTr5+qoRAXlSA/YdAf3lTkTPdaYB7LcqSht6tUzSWL4cFixQpKEcCl7fRS3AbcDDZVSH+Q45AdkOjuwknBF+1LbIxWvyjVcHYqo6NLTYdcuuDVrmPYtt8Dy5aYEhxAuZu+O4hymfyIV+CPH12pgsPNDc67UDAsaOFu9i6tDERXd3r1w330wZQrs3391vyQJUUrY66PYC+xVSi3RWpeLO4icIi4nFzq1XAinSk6G+fNN05LWZk6Eu/xVitLHXtPT11rre4C/rhkeqzBTIMr0+Dz3jGSC3U4TU0lq4ggX2LYN3njDVHp1dzezqv/+d1OKQ4hSxt7w2Keyvt9REoGUNOv6mQBY3L1dHImocFauhLez1vBq3RpeftkU8xOilCqwj0JrHZX18CJwRmt9GqgEtCdvcb8yJyrcLA2ZPDDPontCOFffvuDvD08+CYsWSZIQpZ4jw2M3Ad5KqYbAr8DDwCJnBlUSaql4ot38ad64satDEeXdxYvwwQdgyVr/q1Yt+O47ePBB6ZMQZYIjiUJprZOBkcCHWuu7gDbODcvJMtPp5BZGVbd0V0ciyjOtYfVqGD0avvoKliy5ekz6IkQZ4khRQKWU6g6MBf5ehOeVXodNtfOTNXrQ1sWhiHLq3DnTWb1jh9m+5RYYXOZHlYsKypE3/KcwM7G/zSoT3hzY4NywnMv603TcgM11/iaJQhQvq9VUd50/H1JToVo1ePZZGDLErGEtRBnkyAp3mzD9FNnbJ4GyWzk2JRa3JLNGtme9IBcHI8qddetMjSaAQYNMkpAifqKMc2TholbAs0DTnOdrrctm8ZmkiwAszhyAv69UjBXFbMAA2LDB3EH07u3qaIQoFo40Pa0APgU+x6xCVy7ssgbRV+dZZkOIojl82NxBzJwJ9eublebeesvVUQlRrBxJFJla60+cHklJuXwKADesdGgkxQDFdUpLg88+g//8x/RL/Pvf8MorhT9PiDLIkUTxf0qpx4FvgbTsnVrry06LypkyUwG4QA3c3aRzUVyHP/80dxAREeYOYuxYeOwxV0clhNM4kijGZX1/Lsc+DTQv/nBKTqz2dXUIoqxJSjJLkq5cababNzd3EW1l7Jwo3xwZ9dSsJAIRotQ7dw6+/RY8POCRR+Dhh8HT09VRCeF0jox6qgI8AzTWWk9USrUEbtJar3F6dEK4WlIS+PiYxy1bwgsvQHAwBErVYVFxOFLCYyGQDtyStR0JzHRaREKUBlrD2rUwYgSsX391/513SpIQFY4jiaKF1noWkAGgtU7BrElRNsUccXUEorSLiYF//tPcPcTGmkl0QlRgjnRmpyulKmM6sFFKtSDH6Kcy54hpMYvSMltWXENrU9X1/fchMdE0OT31lLmrEKICcyRRvAr8BDRSSi0BbgUecmZQTpOWAFF70SjiqOrqaERpcukSvPQS7Npltnv2hOnToU4d18YlRCngyKinX5RSfwLdME1OT2mtLzo9MmfYaBYpCm90F4S5OBZRuvj4mFFN1avDc8+ZOk1SxE8IwP6a2U2AK1rrOK31JaVUMjACaKWUmq+1LnuLOURsA+BIywkQFuviYITLnTwJdeuaJOHtDbNnm5XnasiMfSFysteZ/dnMiXEAACAASURBVDXgA6CUCsXUfIrALIX6sfNDcwLlDs37kuQjq9pVaBkZpuTG/febCXTZWrWSJCFEPuw1PVXWWmevjf034Aut9btKKTdgj/NDE8IJDh2CGTPg+HGzrbWp1eTmyABAISome4kiZwNtP8ziRWitrUrabkVZk5pqivgtWWISQ0AAvPwydOrk6siEKPXsJYr1SqmvgSigBrAeQClVHzMBT4iyISEBHngAIiPNncMDD8Cjj5p+CSFEoewliinAGKA+0ENrnZG1vx7worMDE6LY+Pqawn3e3qaIX5s2ro5IiDLF7vBYrfWyfPb9lf1YKaW0LkOr/0TuhOZ9XR2FKAmbN0Pt2tC6tdmePh28vKSInxDXwV4P3gal1BNKqVxDhJRSXkqpfkqpL7lagjxfSqkhSqmjSqnjSqlpds4bpZTSSqnORQv/OqTFO/0lhAvFxsKLL8LTT8O//mVGOIEZAitJQojrYu+OYgjwCLBUKdUMuAJ4A+7AWmCu1rrA0U9KKXfgI2AgppDgLqXUaq31oWvO8wWeBHbcyA9SKGvWKq6BA5z6MsJFtIaffzZzIeLiTDPT8OHg7u7qyIQo8wpMFFrrVMx8iY+VUp5AbSBFa33FwWt3BY5rrU8CKKWWAXcCh64573VgFvBsEWMvmqi95nvWCneiHImONutUb95strt2NXcVDRu6Ni4hyglHaj2R1ZEdVcRrNwTO5NiOBG7OeYJSqgPQSGu9Rinl3EQRnZWfmvSABKe+kihJmZlmEaHz56FqVdPkNHy4lN8Qohg5lCiuU37/U20d31kT9+biQIFBpdREYCJA48bXOas6fIv5XquFJIryxMMDJkyATZtg2jRTgkMIUaycOR01EmiUYzsAOJdj2xdoC2xUSoVjig6uzq9DW2u9QGvdWWvd2f963wgqZ5VmqNXi+p4vSgeLBRYvhq+/vrpv+HCYM0eShBBOYveOQik1AggE9mutfy7itXcBLbM6ws8C9wL3Zx/UWsdh+j2yX2sj8KzWencRX8dxXr5Ou7QoAWFh8PrrpgxHpUowYADUrCnNTEI4mb3qsR8DwcBW4HWlVFet9euOXlhrnamUmgz8jBkp9YXW+qBSagawW2u9+gZjFxVFejosXAhffGHuKOrWNZ3VNWXxKSFKgr07il5Ae621RSlVBdiMGaHkMK31D8AP1+x7pYBz+xTl2kV28ShYMwo/T5Qu+/ebu4iTJ8326NEwebKZFyGEKBH2EkW61toCoLVOVmW9EmBmmm1obGKqJIwyQWuYN88kicaNTRG/Dh1cHZUQFY69RBGklNqX9VgBLbK2FaC11u2cHl1xcveChqZS6Ce/nQCgkqeUli6VMjPNaCalTBPT99+bkU2VKrk6MiEqJHuJonWJRVFizE1R9cpeuCtFHV+pHlqqJCTA+++b9avnzjWJolkz09QkhHAZezOzT+e3P6s0x71AvsfLAqUgJKCaq8MQOf32m5ldffGiqcl06hQ0b+7qqIQQ2JlHoZTyU0pNV0rNV0oNUsYTwEngnpILUZRrly+byq7//KdJEu3awX//K0lCiFLEXtPTYiAW2AaMB54DvIA77RUDFMJhP/0Es2ZBfDxUrmyamEaPlmVJhShl7CWK5lrrEACl1OfARaCx1loKYIjicfKkSRI332w6rRs0cHVEQoh82EsUtjGkWXMpTpXtJKHJUWpKuILVClFRV6u6jh8PgYEwcKDMrhaiFLN3jx+qlIrP+koA2mU/VkqVvdV/TqwHi8yfcJmICLNO9SOPmLsIMCvODRokSUKIUs5eotirtfbL+vLVWnvkeOxXYhEWFw9v8JaRTiXOYoGvvoJ774W//jKT6M6cKfx5QohSw17TU/lqp3HzhHpla45gmXfsGMyYAUeOmO1hw8x6EX5l73OGEBWZvURRRyn1TEEHtdbvOSEeUV4sW2YmzVksUL++6azu1s3VUQkhroO9ROEOVCX/BYiEsK95c9N5PWYMTJoEVaq4OiIhxHWylyiitNYzSiwSUbYlJ8P27dCvn9nu2hX+9z9o1Mj+84QQpZ69zmy5kxCO2b7d3Dk8/zzsyTEXU5KEEOWCvTuK/iUWhSgTMjIyiIyMJDXVlGvHajVDXVNSTCe1p6eZVX34sGsDFcIJvL29CQgIwNPT09WhlDh7RQEvl2QgovSLjIzE19eXpk2bohIS4Px5qFYNqlc361XXqiVzIkS5pLXm0qVLREZG0qxZM1eHU+KkqI5wWGpqKrVq1UJdvgyRkWbdiCpVTMd17dqSJES5pZSiVq1aV++mKxh7TU9C5KGUMncRly+bO4gaNSRBiAqhrC/yeSMq3B3FXxGxHDlfhktWucK5c6bKq86ag+nhYWo01awpSUKICqBCJYo/I2K56+OtAHh5uLs4mjLAaoXly82Ipq+/hqSkq8dclCCqVq0KQHh4OG3bts117LXXXmPOnDm27Tlz5hAUFETbtm1p3749X331Vb7XLOi8Pn36cNNNN9GuXTuCgoKYPHkyV65cyfVci8VChw4duOOOO2z7HnroIZo1a0ZoaCihoaHsyRoJtmjRIvz9/QkNDSUoKIi5c+cW+HOuWrWKGTNm2H4upRTHjx+3HZ87dy5KKXbv3l3o7+xGzZ07F29vb+Li4mz7Fi1axORrVh7s06ePLZ7ExEQeffRRWrRoQXBwML169WLHjh0Ovd7ly5cZOHAgLVu2ZODAgcTGxuZ73tSpUwkODqZ169Y8+eSTaK1JSEiw/d5DQ0OpXbs2U6ZMAWDTpk107NgRDw8PVq5cabtOTEwMQ4YMKdLvpKIp94li3row3vnpCAlpmfwZYf7gPhnbkQ/v6+DiyEq58HCzTvXs2WZU04ABZWrS3Keffsovv/zCzp07OXDgAJs2bULrvFVpCjtvyZIl7Nu3j3379lGpUiXuvPPOXM+fN28erVvnXTV49uzZ7Nmzhz179hAaGmrbP2bMGPbs2cOWLVt44403OFNA3atZs2bx+OOP27ZDQkJYtmyZbXvlypW0adPG8V/IDVi6dCldunTh22+/dfg548ePp2bNmoSFhXHw4EEWLVrExYsXHXru22+/Tf/+/QkLC6N///68/fbbec7ZunUrW7ZsYd++fRw4cIBdu3bx22+/4evra/u979mzhyZNmjBy5EgAGjduzKJFi7j//vtzXcvf35/69euzZcsWh3++iqbcJgqrVXPr2+uZu+4Yn2w8Ydv/6d86cltIfRdGVsplZsIXX8B998HevaYfYvZsePvtMrWg0JtvvsnHH3+MX1ZdqWrVqjFu3LjrPs/Ly4tZs2YRERHB3r17ATMK7Pvvv2f8+PFFjq9WrVoEBgYSFRWV59ixY8eoVKkStWvXtu0bMWIE3333HQAnT56kWrVq+Pv7246vXbuW7t2707FjR0aPHk1iYiIAM2bMoEuXLrRt25aJEyfakmCfPn14/vnn6dq1K61atWLz5s35xnnixAkSExOZOXMmS5cudehnO3HiBDt27GDmzJm4Zf3NNG/enKFDhzr0/O+++872bzBu3DhWrVqV5xylFKmpqaSnp5OWlkZGRgZ169bNdU5YWBjR0dH07NkTgKZNm9KuXTtbTDmNGDGCJUuWOBRfRVTmO7O11kTGppCSYWFfZBzLdkaQnG7hUNTVSujrnumN7+cejO/YHNpKkrBr/Xr4+GPz+M474amn8i3i96//O8ihc8Vbbb5NAz9eHRbs8PknTpzI9Wn9/PnzPPvssyQkJJCQkECLFi3sPt/R87K5u7vTvn17jhw5Qvv27ZkyZQqzZs0iISFvn9eLL77IjBkzbJ+IK1WqlOt4REQEqamptGuXt1Dlli1b6NixY659fn5+NGrUiAMHDvDdd98xZswYFi5cCMDFixeZOXMm69atw8fHh3feeYf33nuPV155hcmTJ/PKK68A8MADD7BmzRqGDRsGQGZmJjt37uSHH37gX//6F+vWrcsTy9KlS7nvvvvo2bMnR48eJTo6mjp16tj9PR08eJDQ0FDc3fNv3u3Zs2e+v7M5c+YwYMAALly4QP365v9p/fr1iY6OznNu9+7d6du3L/Xr10drzeTJk/Pc2S1dupQxY8Y41AnduXNnXnrppULPq6jKfKLoPHMdl5LS8+zvc5M/iamZLH+0O+5u0uFql9ZX+xwGDICtW+H2200ZjlKsRYsWtvZ/MG35YD48OPLm4Oh51z4HYM2aNdSpU4dOnTqxcePGXOe89dZb1KtXj/T0dCZOnMg777xje7Nevnw5GzZs4OjRo/z73//G29s7z2tERUXlulvIdu+997Js2TJ+/vlnfv31V1ui2L59O4cOHeLWW28FID09ne7duwOwYcMGZs2aRXJyMpcvXyY4ONiWKLKbZDp16kR4eHi+P++yZcv49ttvcXNzY+TIkaxYsYJJkyYV+Htz5PdZ0N1LURw/fpzDhw8TGRkJwMCBA9m0aRO9evXKFfvixYsdul6dOnU4d+7cDcdVXpX5RJGSYQHgvXvaoxSENqpBg+reVJLOasfs2WOalmbNMivPublB1huuPUX55F/S/Pz88PHx4eTJkzRv3vyGz8tmsVjYv38/rVu3ZunSpaxevZoffviB1NRU4uPj+dvf/sZ//vMf26fhSpUq8fDDD+fqYB8zZgzz589n27ZtDB06lNtuu4169erlep3KlSvn6jjONmzYMJ577jk6d+5sayoDk7wGDhyYp2koNTWVxx9/nN27d9OoUSNee+21XPMAsu9y3N3dyczMzPN6+/btIywsjIEDBwImATVv3pxJkyZRq1atPJ3Mly9fpnbt2lSvXp29e/ditVrzbeYp7I6ibt26REVFUb9+faKiovK9g/n222/p1q2bbXDDbbfdxvbt222JYu/evWRmZtKpU6c8z81PamoqlStXdujciqjsNDoXwMNN8dAtTRnZMYC7OgTQrLZP3iShNaTLkNhckpNNchg/Ho4ehUWLXB1RsZo+fTqTJk0iPms1vfj4eBYsWHDd52VkZDB9+nQaNWpEu3bteOutt4iMjCQ8PJxly5bRr18//vOf/wDY+h201qxatSrP6CwwTScPPPAA8+bNy3OsdevWuUY4ZatcuTLvvPMOL774Yq793bp1Y8uWLbbnJCcnc+zYMVtSqF27NomJiblG+jhi6dKlvPbaa4SHhxMeHs65c+c4e/Ysp0+fpkuXLmzZsoXz588DsHv3btLS0mjUqBEtWrSgc+fOvPrqq7Y7sLCwMFsfy+bNm3N1OGd/DRgwAIDhw4fz5ZdfAvDll1/mGUAApmP6t99+IzMzk4yMDH777bdcTU/ZTWaOOnbsWL7/TsIo84nCIQnmj5m0vJ/SKqStW+Gee8yQV3d3kyyee87VURWrf/zjH/Tt29fWkdu7d2+q5DNqq7Dzxo4dS7t27Wjbti1JSUm2Nzt7xo4dS0hICCEhIVy8eLHAtu/nn3+ehQsX5vl03atXL/766698R2nde++9efov/P39WbRoEffddx/t2rWjW7duHDlyhOrVqzNhwgRCQkIYMWIEXbp0KTT2nJYtW8Zdd92Va99dd93FsmXLqFu3LvPmzeP2228nNDSUKVOmsHTpUtsdxOeff8758+cJDAwkJCSECRMm0KBBA4ded9q0afzyyy+0bNmSX375hWnTpgEmGWUPHBg1ahQtWrQgJCSE9u3b0759e1uTGsDXX3+dJ1Hs2rWLgIAAVqxYwaOPPkpw8NW74g0bNjjc2V4Rqfz+GEuzzp0765xjx9u99jMjOwbw2nA7TSFxZ2FuGxj2AXTKO6KlwoiPh3ffhe+/N9utW8Mrr0DLlg49/fDhw/kOBRXF76mnnmLYsGG2T9nCuXr16sV3331HjRo17J5Xlv8PKKX+0Fp3vp7nVow7CmFcvAg//wxeXvDkk6a5ycEkIUrWCy+8QHJysqvDqBBiYmJ45plnCk0SFVmZ78wWhYiLM8NblTLF+155Bdq2hcaNXR2ZsKNu3boMHz7c1WFUCP7+/owYMcLVYZRqFeOOYsv75ruqGD8uYDrwV6+GESPgl1+u7r/9dkkSQogiqRjvnPFZ46Nvus21cZSUc+fMOtUzZkBCAkhpAiHEDagYTU9u7uAfBD61Cz+3LMsu4vfRR5CaasqBP/ssSMEzIcQNqBiJImI7eFd3dRTOFR0N06bBvn1me9AgkyRq1nRtXEKIMs+pTU9KqSFKqaNKqeNKqWn5HH9GKXVIKbVPKfWrUqpJUa6fYbESn5p3Rmkelfwg+VJRLl32VKsGV66YJUnfew/efLNcJgl3d3dCQ0MJDg6mffv2vPfee1itVsBMNMuew9C2bVt69OhhK46X/by2bdsyevTofEcU2SuNbe91s0VERFC1atVcM7GbNm1KSEgIoaGhdO58dWRizlLk7du359dffy3wZ54yZQqbNm0CTDG/xo0b55pjMWLECNsMZWe78847beVBsj300EN5JvPljOfYsWPcfvvtBAYG0rp1a+655x4uXLjg0Ov98ccfhISEEBgYaCslfq24uDiGDRtG+/btCQ4OtpU2ARgyZAjVq1fPVQYeYP78+QQGBqKUylXVds2aNbz66qsOxVahaK2d8gW4AyeA5oAXsBdoc805fYEqWY//ASwv7LqdOnXS2b7ceko3eX6Nnrpir7brwy5aL3/Q/jll0aFDWsfHX90+cSL3drG/3CGnXdtRPj4+tscXLlzQ/fv316+88orWWus333xTP/3007bjR44c0ampqXmed//99+t33303z7XHjBmjp02bpi0Wi9Za6xMnTug1a9YU+rrZRo4cqUeNGqVnz55t29ekSRMdExOT57XGjRunV6xYobXWev369TowMDDfn/fSpUv65ptvtm337t1bh4SE6M2bN2uttY6NjdVdu3bNFZ+zxMbG6oCAAB0UFKRPnjxp25/zZ8mWHU9KSooODAzUq1evth1bv3693r9/v0Ov2aVLF71161ZttVr1kCFD9A8//JDnnDfeeENPnTpVa611dHS0rlGjhk5LS9Naa71u3Tq9evVqPXTo0FzP+fPPP/WpU6fy/PtYrVYdGhqqk5KS8o2nNPwfuF7Abn2d7+fOvKPoChzXWp/UWqcDy4Bcc/G11hu01tkf7bYDAUV5gTX7TKmER3sXXqenXElLgw8+gHHjzPdszZuDr6/r4iphderUYcGCBcyfPx+tNVFRUTRs2NB2/KabbspTtRVMraFrS2QUpTT2ta8LZqGh5s2b55rt66ju3btz9uzZfI+tXLkyz6I62cUBAf73v//Zivtlmz17Nl26dKFdu3a5Ph2PGDGCTp06ERwcnKtMSdWqVXnxxRdp37493bp1K/DT/jfffMOwYcNyvX5h/vvf/9K9e/dcs6b79u3rULmMqKgo4uPj6d69O0opHnzwwQJLjickJKC1JjExkZo1a+LhYVrV+/fvj28+/yc6dOhA06ZN871Wnz59WLNmjUM/X0XhzD6KhkDOVVkigZvtnP934Mf8DiilJgITwdR4ATgRk8jOU5fx9nSjuX/J3HaXCn/+Ca+/DmfOmAJ+Varkrv5aUn6cBuf3F+8164XAbXkXqbGnefPmWK1WoqOjeeSRRxg0aBArV66kf//+jBs3jpbXTCjMzMzkxx9/zPPmW1hpbHuvW7VqVd555x1++eWXXM1OYN54Bg0ahFKKRx99lIkTJ+a51k8//VTgOP4tW7YwatSoXPv69+/PhAkTsFgsLFu2jAULFvD6668DZl2KsLAwdu7cidaa4cOH26qqfvHFF9SsWZOUlBS6dOnC3XffTa1atUhKSqJbt2688cYbTJ06lX//+9/5lh1ZunQpr776KnXr1mXUqFFMnz690N/TgQMHCizMd/ToUcaMGZPvsY0bN3L27FkCAq5+dgwICMg3oU6ePJnhw4fToEEDEhISWL58eb7FCB3VuXNnNm/ezD333HPd1yhvnJko8nvnyrdeiFLqb0BnoHd+x7XWC4AFYEp4ABw4a+o2Pd4nsBhCLQOSkszdwzffmO2ck+cquOxP9aGhoZw8eZK1a9eybt06unTpwrZt22jdujUpKSm2tSt69uzJ3//+92J73VdffZWnn346336CLVu20KBBA6Kjoxk4cCBBQUG2CqfPPfccU6dOJTo6mu3bt+f7GvmVHHd3d6dHjx4sX76clJSUXJ+M165dy9q1a+nQwazgmJiYSFhYGL169eKDDz6wrVJ35swZwsLCqFWrFl5eXrY2/E6dOvFLznk3WS5cuMDx48fp0aMHSik8PDw4cOAAbdu2zbe0uCPlxm+66aZcZeKvlf37Ley6P//8M6Ghoaxfv54TJ04wcOBAevbsmavCblFIyfG8nJkoIoFGObYDgDy/faXUAOBFoLfWOs3Ri0fFmcqYg4PrFXJmORAfb1acu3ABPDzgkUfg4YfB09N1MRXxk7+znDx5End3d1sp6qpVqzJy5EhGjhyJm5sbP/zwA61bt6Zy5cp235SCg4Ptlsa297o7duxg5cqVTJ06lStXruDm5oa3tzeTJ0+2FcKrU6cOd911Fzt37rQlitmzZzNy5Eg++OADxo0bxx9//JHndSpXrpyrNHi2e++9l7vuusu2Bkc2rTXTp0/n0UcfzbV/48aNrFu3jm3btlGlShX69Olju66np6ftDbigkuPLly8nNjaWZs2aAabK7rJly5g5c2aekuPZ5cazf6+//fZbvr/Dwu4oAgICbOtNgFlRML/CggsXLmTatGkopQgMDKRZs2YcOXKErte5noqUHM/LmX0Uu4CWSqlmSikv4F5gdc4TlFIdgM+A4VrrvMtY2eHpbkKvVy3vwi95XDxalEuXPn5+0KULtGkD//kPTJzo2iRRSsTExPDYY48xefJklFJs2bLF9oaVnp7OoUOHaNLEsYF0hZXGtve6mzdvtpXinjJlCi+88AKTJ08mKSnJVhk2KSmJtWvX5mmbd3Nz46mnnsJqtfLzzz/nea2CSo737NmT6dOn56mQOnjwYL744gvbaK+zZ88SHR1NXFwcNWrUoEqVKhw5cqTAO5iCLF26lJ9++sn2c/7xxx+2foo+ffqwfPly0tPNAmKLFi2ib9++ANx///1s3bqV77MLUWKa2vbv32+7o8jvq3r16tSvXx9fX1+2b9+O1pqvvvqqwJLj2aPGLly4wNGjRx1aX6QgUnI8L6clCq11JjAZ+Bk4DHyttT6olJqhlMouYjMbqAqsUErtUUqtLuByNxKI+V6WhsdqDWvXwv4cfQDPP2+K+AVWkKa2AmQ3IQUHBzNgwAAGDRpk67A9ceIEvXv3JiQkhA4dOtC5c2fuvvtuh69trzS2vdctyIULF+jRowft27ena9euDB06NE/fCJjmlJdeeolZs2blOTZ06NA8K+hlP+fZZ5/Nta42wKBBg7j//vvp3r07ISEhjBo1ioSEBIYMGUJmZibt2rXj5Zdfplu3bg7/XsLDw4mIiMj1nGbNmuHn58eOHTu444476NmzJ506dSI0NJQtW7bwzjvvAOaOaM2aNXz44Ye0bNmSNm3asGjRokKXU832ySefMH78eAIDA2nRogW33WaqK3z66ad8+umnALz88sts3bqVkJAQ+vfvzzvvvGP7vfTs2ZPRo0fz66+/EhAQYEvGH3zwge2OpV27drnWPZeS4/m43uFSrvrKHh77+eaTusnza/SV5HT7Y8KSLmn9qp/W69+0f15pceGC1k8/rXWnTlqPGqV11jC/0qAsDw0sy2699VYdGxvr6jAqhPPnz+t+/foVeLws/x/gBobHlv+Z2RfDzPdKpXxklNawahW8/77puPbxMf0SHuX/n0jY9+677xIREUH16uW8ukApEBERwbvvvuvqMEqdivMuVKeNqyMoWGQkzJwJ2Qsy9ewJ06eDg7fnony7+WZ7o8pFcSrqKoAVRflPFDGHzffSupJfZqbpnI6OhurVYepUGDiw5OdFCCFEAcp/oshWo6mrI8ifh4cpCb59O/zznyZZCCFEKVJxEoVXFVdHYGRkwMKFpg9i7Fizb+hQ8yWEEKVQxUkUpcHBg2YxoRMnzLrVt98Osk6vEKKUK/8r3EUfdnUEZhGh9983s6lPnIBGjeDDDyVJXIec5cKHDRvGlStXiuW64eHhTptk9dVXX9G2bVuCg4Np06aNrR5Unz592J09gEGIUqx8J4r0ZNhhJuXg5aLhsbt3w5gxZkY1wIMPwtKlUEChNGFfdimOAwcOULNmTT766CNXh2TXjz/+yPvvv8/atWs5ePAgf/75J9WqVXN1WEIUSflOFKfMYi80vgW8r69A2A3RGhYsgLNnzYzqRYvgySfB24GyI6JQOctzJyYm0r9/fzp27EhISIit9EZ4eDitW7dmwoQJBAcHM2jQIFJSUgCzKE779u3p3r17roSTmprKww8/bJvhvWHDBsCUphgxYgTDhg2jWbNmzJ8/n/fee48OHTrQrVs3Ll++nCfGt956izlz5thmeHt7ezNhwgTb8RUrVtC1a1datWrF5s2bbTH37NmTjh070rFjR7Zu3QqY+kd9+vRh1KhRBAUFMXbsWFu5kV27dnHLLbfYZoEnJCRgsVh47rnnbCXHP/vss2L9/YsK5Hpn6rnqq0gzs3f9PzMrO3xLUSYw3rics6lPn9b688+1Ti9kBnkZkGdWaqdOBX99883V8775xv65RZC9IE5mZqYeNWqU/vHHH7XWWmdkZOi4uDittdYxMTG6RYsW2mq16lOnTml3d3f9119/aa21Hj16tF68eLHWWuuQkBC9ceNGrbXWzz77rA4ODtZaaz1nzhz90EMPaa21Pnz4sG7UqJFOSUnRCxcu1C1atNDx8fE6Ojpa+/n56U8++URrrfWUKVP03Llz88Rbo0YNfeXKlXx/lt69e+tnnnlGa631999/r/v376+11jopKUmnpKRorbU+duyYzv6b37Bhg/bz89NnzpzRFotFd+vWTW/evFmnpaXpZs2aBfAVBAAAD5tJREFU6Z07d2qttY6Li9MZGRn6s88+06+//rrWWuvU1FTdqVOnXAsOiaKTmdnlUtZchBrNSublYmNhzhyzJOn8+WYuROPGUAwlrYWRXXMpPDycTp06MXDgQMB84HnhhRfYtGkTbm5unD171rYAT/aSo2DKaIeHhxMXF8eVK1fo3dtUtn/ggQf48UezHMrvv//OE088AUBQUBBNmjTh2LFjgFl0x9fXF19fX6pVq2ZbkCckJIR92euVF0H2okPZcQFkZGQwefJk9uzZg7u7u+21Abp27WpboyH791CtWjXq169vmyyWXV577dq17Nu3z7ZMaVxcHGFhYbYKsEI4qpwnihKiNfz8M8yeDXFxpmkpPBzK+39IRztiR440X8Ugu48iLi6OO+64g48++ognn3ySJUuWEBMTwx9//IGnpydNmza1ldHOucqdu7s7KSkpaK0LXDNB25mcmfNabm5utm03N7d8y3MHBwfzxx9/0K9fP7vXy1nee+7cudStW9dW9tw7R1PltT9LZmZmgT+L1poPP/yQwYMHF/jzCOGI8t1HcWan81/jwgV4+ml46SWTJLp2heXLy3+ScLFq1arxwQcfMGfOHDIyMoiLi6NOnTp4enqyYcMGTp8+bff51atXp1q1avz+++8ALFmyxHasV69etu1jx44RERHBTTfddF1xTp8+nalTp3L+/HkA0tLS+CDn8rX5iIuLo379+ri5ubF48WIsFovd84OCgjh37hy7du0CICEhgczMTAYPHswnn3xCRkaG7WdJSkq6rp9DVGzl+47iUlZBQB9/++ddr+++g3ffheRkqFoVnnkGhg2T8hslpEOHDrRv355ly5YxduxYhg0bRufOnQkNDSUoKKjQ5y9cuJBHHnmEKlWq5PrU/fjjj/PYY48REhKCh4cHixYtynftbUfcfvvtXLhwgQEDBtg++T/yyCN2n/P4449z9913s2LFCvr27YuPj4/d8728vFi+fDlPPPEEKSkpVK5cmXXr1jF+/HjCw8Pp2LEjWmv8/f3zXXNaiMIoe7fZpVHnzp317t27+X+/n+L1NYfY++ogqlUuYBGf16oBCl4rnrH2eSxYYL5694Zp08DfSQmplDh8+DCtW7d2dRhCuExZ/j+glPpDa935ep5bvu8oUFCvGCdRWSxw5gxkr1H88MMQFGSqvcpdhBCinCrffRRoaN6neC4VFmYSw8SJpi8CzHKkvXpJkhBClGvl947iyhnzPf0GO+/S0+GLL0whP4sF6taFqCiooLNr7Y0WEqI8K2vN9MWpzCaKtEz7I0FIzfrU36Dj9b/I/v3w+utw8qTZHj0aJk82lV8rIG9vby5dukStWrUkWYgKRWv9/9u78yApizOO498fCAKRSxHFIyIoKhIORcV4oaKFiEAlGqFKI8ZgSmNiaTxi0GiZVEk0R0WlJIiIEg+EhIQQLcWDgAoiyBFBvFARMYIgR+ICKk/+6B5mHGZnX7aYc59P1da+87799vZ07c6z3f2+z8u6deu+dqlyQ1KxgWLeeyFdQtPGtcyevXxP+N60nh/qDz8cEveZhZvmbrkFevWqX11VIvUw+rVr15a6Kc4VXbNmzXbc7NjQVGygePuT/wLQvGnj3AW2heP1XqM4+mho1AguvhhGjIB6Xh5ZTZo0aeJ39TrXAFVsoGjfak8a5VuKX78C2h8NLfZOVuHmzTB7dnhGBITsrtOmhTUJ55xrwCoyUCxdvZGFKzdwYqd9ai9UswG+rElW4cyZMGoUfPop7L8/HBPXNTxIOOdcZQaKp/4d0iGc1bWWD/IV/4LNq+GoQfkrWr8e7rwTnn02vO7eHfZOOAJxzrkGoiIDhREuU/vBybXMl38Uk9X1uriWCgyeeipket20CZo3D1czXXAB+eeznHOu4am4QGHA6BfezV9o9u/D9059cx9/5JHwaFKAE06AkSMhPljGOefc11VcoFizaQuNgVbNamn6l9vSVzzt0TR3mYEDYerUcKf1uef6ndXOOZdHxc2zrNm8FYDnr+u788Ht2+HB/mG74ynp/R98ALffDjHdMm3awOTJIWB4kHDOubwqbkQBcPoR+9Jurxz3NTxxMXy0IGwPuS+k3Jg4MWR43bYt3Dg3fHg47msRzjmXSEUGih4Ht9l55zM3w/LpYfvqJbCmBq66BJYvD/vOO2+3PWXNOecakooMFDnNfzB8HzwWHvsnTJgQpqI6dAiL1X36lLR5zjlXqSoyUOzfKisx14fzwgL2cSNgbVsYf1dYexg6FK68Elq0KE1DnXOuClRkoOjdsW36xbp3YVy/EBgOPBZ6nAULFoRUHD16lK6RzjlXJSpyRfegtnGE8PYMuLE7TKmBjsOg57AQMG66yYOEc87tJgUNFJL6S3pT0juSfp7j+J6SJsXjr0jqmKTeZk0aw/JZMOI8eHILfLYdPj5kdzffOeccBQwUkhoDo4FzgK7AMElds4pdBnxmZocBfwB+k6jyfzwB554Bb30B+3aCO8bDDTfsxtY755xLKeQaxfHAO2a2AkDS48BgYFlGmcHAbXF7CnCvJFmeZw4etPkT+NGw8KLLIfDQs3CIjyacc65QCjn1dCDwYcbrVXFfzjJm9iWwEciTOxzabtsMTQQn7QkzlnuQcM65AivkiCJXbozskUKSMki6HLg8vtyqlV+9zsoamNLgnzrXDvi01I0oE94Xad4Xad4XaUfU98RCBopVwMEZrw8CVtdSZpWkPYDWwPrsisxsLDAWQNJ8M+tdkBZXGO+LNO+LNO+LNO+LNEnz63tuIaeeXgUOl3SopKbAUGBaVplpwCVx+3zg+XzrE84554qvYCMKM/tS0lXA00BjYLyZLZV0OzDfzKYBDwATJb1DGEkMLVR7nHPO1U9B78w2syeBJ7P2/TJjewtwwS5WO3Y3NK1aeF+keV+keV+keV+k1bsv5DM9zjnn8qnIFB7OOeeKp2wDRaHSf1SiBH1xraRlkpZIek5S1d5cUldfZJQ7X5JJqtorXpL0haTvxd+NpZIeLXYbiyXB38g3Jb0gaWH8OxlQinYWmqTxktZIer2W45J0d+ynJZKOSVSxmZXdF2Hx+12gE9AUWAx0zSpzJTAmbg8FJpW63SXsi9OBFnH7iobcF7FcS2AWMBfoXep2l/D34nBgIdA2vm5f6naXsC/GAlfE7a7A+6Vud4H64lTgGOD1Wo4PAJ4i3MPWB3glSb3lOqLYkf7DzLYBqfQfmQYDD8XtKcCZUlU+ALvOvjCzF8zs8/hyLuGelWqU5PcC4FfAncCWYjauyJL0xQhgtJl9BmBma4rcxmJJ0hcGtIrbrdn5nq6qYGazyHEvWobBwMMWzAXaSOpQV73lGigKkv6jQiXpi0yXEf5jqEZ19oWkXsDBZja9mA0rgSS/F12ALpJekjRXUv+ita64kvTFbcBFklYRrsT8SXGaVnZ29fMEKN8HF+229B9VIPH7lHQR0Bs4raAtKp28fSGpESEL8fBiNaiEkvxe7EGYfupLGGXOltTNzDYUuG3FlqQvhgETzOx3kk4k3L/Vzcy2F755ZaVen5vlOqLYlfQf5Ev/UQWS9AWS+gEjgUFmtrVIbSu2uvqiJdANmCnpfcIc7LQqXdBO+jfydzP7wszeA94kBI5qk6QvLgOeADCzOUAzQh6ohibR50m2cg0Unv4jrc6+iNMtfyIEiWqdh4Y6+sLMNppZOzPraGYdCes1g8ys3jluyliSv5G/ES50QFI7wlTUiqK2sjiS9MVK4EwASUcRAsXaorayPEwDvh+vfuoDbDSzj+s6qSynnszTf+yQsC/uAvYCJsf1/JVmNqhkjS6QhH3RICTsi6eBsyUtA74CrjezdaVrdWEk7IufAfdLuoYw1TK8Gv+xlPQYYaqxXVyPuRVoAmBmYwjrMwOAd4DPgUsT1VuFfeWcc243KtepJ+ecc2XCA4Vzzrm8PFA455zLywOFc865vDxQOOecy8sDhStbkr6StCjjq6OkvpI2xiygb0i6NZbN3L9c0m/z1NtL0ri4faSkOZK2Srouzzn7SZouaXHMxvpkbWXr+V57S7o7bu8p6dn4ni+UNE5S1zznDkplTJU0JF/ZjHMel1SNN9+5AijL+yici2rMrGfmDoV08rPNbKCkbwCLJKXyOqX2NwcWSppqZi/lqPcXwK/j9nrgp8CQOtpyOzDDzP4Y29G9Xu+oFvGmwNSNgb2AJhnvfVId504jfYPZEGA6sKyOH3kfcAMhcaBzefmIwlUsM/sfsADonLW/BlhEjmRnkloC3c1scSy7xsxeBb6o48d1IKQ/SP2MJbG+vpJmSZoaRxpjYs4pJJ0dRyuvSZosaa+4/zhJL8fRyTxJLWM90yW1B/4M9Iwjis6SZqbSkCg8d+G1eO5zcd9wSfdK+jYwCLgr49zXMt774ZIWxJezgX4x/Y1zeXmgcOWseca009Tsg5L2IeRzWpq1vy0hp9GsHHX2BnI+1KUOo4EHFB5+M1LSARnHjifc+fstQtD6TkyZcTPQz8yOIYwWro0pJiYBV5tZD6AfUJOqKKZg+SFhdNTTzN7NeF/7AvcD343nfu1582b2MmFkcX3GuRslpUYmlwITYtnthLtze9SjL1wD4/9NuHK209RTdIqkhcB2YFRM19A37l8CHBH3/yfHuR2oR44fM3taUiegP3AOYWqrWzw8z8xWwI4UCicTnoXRFXgpplVpCsyJbfs4jmIws03xvCTN6APMign+MLMkSTDHAZdKuha4kBDUUtYABxBGZc7VygOFq0SzzWxgbfsldQFejGsUi7LK1BASwuUl6cek5+8HmNnq+MH8KPBoXBc5FVjHzmmajZDOeYaZDcuqt3uO8kmpHuf+hZDv53lgQVaup2ZkjGacq41PPbmqY2ZvAXcAN+Y4/AZwWII6Rsfpm55mtlrSGZJawI51js6EjKQAxytkLm1E+K/9RULm2pMkHRbPaRED2HLgAEnHperahXWCOcBpkg6N5+6do8xmQrr11PvYQkiWdx/wYFbZLmRN2zmXiwcKV63GAKemPlRTzGw50Dp+2CNp/5hl81rgZkmrJLXauTqOBebHqa05wLjU9FF8PYqw9vEeMNXM1hIeoPRYPGcucGR8VOeFwD2SFgMzSDDCiW1fC1wO/DWem+tqqMeB6+NlwqlF/kcII5FnUoUk7UeY2qszxbRznj3WNTgKqaY3m9m43VBXX+C6WqbCykK8P6S1md2Sse8aYJOZPVC6lrlK4WsUriG6j6wrhqpVvFqsM3BG1qENwMTit8hVIh9ROOecy8vXKJxzzuXlgcI551xeHiicc87l5YHCOedcXh4onHPO5eWBwjnnXF7/ByNdWMbA+7QeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot roc\n",
    "f, ax = plt.subplots()\n",
    "ax.step(roc_scores.loc[pd.IndexSlice[\"qpath_CD45R\", \"fpr\", \"mean\"], :],\n",
    "        roc_scores.loc[pd.IndexSlice[\"qpath_CD45R\", \"tpr\", \"mean\"], :],\n",
    "        label=\"IHC CD45R (Mean AUC=0.871)\",)\n",
    "ax.step(roc_scores.loc[pd.IndexSlice[\"DSP_CD45R\", \"fpr\", \"mean\"], :],\n",
    "        roc_scores.loc[pd.IndexSlice[\"DSP_CD45R\", \"tpr\", \"mean\"], :],\n",
    "        label=\"DSP CD45R (Mean AUC=0.811)\",)\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=2, color=\"r\", label=\"Random Chance\", alpha=0.8)\n",
    "ax.legend(bbox_to_anchor=(.27, .33))\n",
    "ax.set_xlim([0, 1])\n",
    "ax.set_xlabel(\"FPR (1-Specificity)\")\n",
    "ax.set_ylim([0, 1])\n",
    "ax.set_ylabel(\"TRP (Sensitivity)\")\n",
    "f.suptitle(\"ROC\")\n",
    "#plt.savefig(\"CD45R_DSP_qpath_L2.pdf\", transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th colspan=\"7\" halign=\"left\">p-value</th>\n",
       "      <th colspan=\"7\" halign=\"left\">adj_p-value</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>fit_time</th>\n",
       "      <th>score_time</th>\n",
       "      <th>test_accuracy</th>\n",
       "      <th>test_average_precision</th>\n",
       "      <th>test_brier_score_loss</th>\n",
       "      <th>test_neg_log_loss</th>\n",
       "      <th>test_roc_auc</th>\n",
       "      <th>fit_time</th>\n",
       "      <th>score_time</th>\n",
       "      <th>test_accuracy</th>\n",
       "      <th>test_average_precision</th>\n",
       "      <th>test_brier_score_loss</th>\n",
       "      <th>test_neg_log_loss</th>\n",
       "      <th>test_roc_auc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>qpath_CD45R</th>\n",
       "      <th>DSP_CD45R</th>\n",
       "      <td>0.0005</td>\n",
       "      <td>0.2983</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0003</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.5799</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0005</td>\n",
       "      <td>0.2983</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0003</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.5799</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       p-value                           \\\n",
       "                      fit_time score_time test_accuracy   \n",
       "qpath_CD45R DSP_CD45R   0.0005     0.2983           0.0   \n",
       "\n",
       "                                                                    \\\n",
       "                      test_average_precision test_brier_score_loss   \n",
       "qpath_CD45R DSP_CD45R                 0.0003                   0.0   \n",
       "\n",
       "                                                     adj_p-value             \\\n",
       "                      test_neg_log_loss test_roc_auc    fit_time score_time   \n",
       "qpath_CD45R DSP_CD45R            0.5799          0.0      0.0005     0.2983   \n",
       "\n",
       "                                                            \\\n",
       "                      test_accuracy test_average_precision   \n",
       "qpath_CD45R DSP_CD45R           0.0                 0.0003   \n",
       "\n",
       "                                                                            \n",
       "                      test_brier_score_loss test_neg_log_loss test_roc_auc  \n",
       "qpath_CD45R DSP_CD45R                   0.0            0.5799          0.0  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calculate_significance(scores).round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}