05_Regression_select_importance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df1  = pd.read_csv('QSAR/ER_alpha_Fingerprinter.csv'               , header = 0)\n",
    "df2  = pd.read_csv('QSAR/ER_alpha_ExtendedFingerprinter.csv'       , header = 0)\n",
    "df3  = pd.read_csv('QSAR/ER_alpha_EStateFingerprinter.csv'         , header = 0)\n",
    "df4  = pd.read_csv('QSAR/ER_alpha_GraphOnlyFingerprinter.csv'      , header = 0)\n",
    "df5  = pd.read_csv('QSAR/ER_alpha_MACCSFingerprinter.csv'          , header = 0)\n",
    "df6  = pd.read_csv('QSAR/ER_alpha_PubchemFingerprinter.csv'        , header = 0)\n",
    "df7  = pd.read_csv('QSAR/ER_alpha_SubstructureFingerprinter.csv'   , header = 0)\n",
    "df8  = pd.read_csv('QSAR/ER_alpha_SubstructureFingerprintCount.csv', header = 0)\n",
    "df9  = pd.read_csv('QSAR/ER_alpha_KlekotaRothFingerprinter.csv'    , header = 0)\n",
    "df10 = pd.read_csv('QSAR/ER_alpha_KlekotaRothFingerprintCount.csv' , header = 0)\n",
    "df11 = pd.read_csv('QSAR/ER_alpha_AtomPairs2DFingerprinter.csv'    , header = 0)\n",
    "df12 = pd.read_csv('QSAR/ER_alpha_AtomPairs2DFingerprintCount.csv' , header = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chemblId</th>\n",
       "      <th>pIC50</th>\n",
       "      <th>FP1</th>\n",
       "      <th>FP2</th>\n",
       "      <th>FP3</th>\n",
       "      <th>FP4</th>\n",
       "      <th>FP5</th>\n",
       "      <th>FP6</th>\n",
       "      <th>FP7</th>\n",
       "      <th>FP8</th>\n",
       "      <th>...</th>\n",
       "      <th>FP1015</th>\n",
       "      <th>FP1016</th>\n",
       "      <th>FP1017</th>\n",
       "      <th>FP1018</th>\n",
       "      <th>FP1019</th>\n",
       "      <th>FP1020</th>\n",
       "      <th>FP1021</th>\n",
       "      <th>FP1022</th>\n",
       "      <th>FP1023</th>\n",
       "      <th>FP1024</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHEMBL370037</td>\n",
       "      <td>8.221849</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CHEMBL189073</td>\n",
       "      <td>5.762708</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 1026 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       chemblId     pIC50  FP1  FP2  FP3  FP4  FP5  FP6  FP7  FP8   ...    \\\n",
       "0  CHEMBL370037  8.221849  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0   ...     \n",
       "1  CHEMBL189073  5.762708  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ...     \n",
       "\n",
       "   FP1015  FP1016  FP1017  FP1018  FP1019  FP1020  FP1021  FP1022  FP1023  \\\n",
       "0     1.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0     1.0   \n",
       "1     0.0     0.0     0.0     0.0     0.0     0.0     1.0     0.0     0.0   \n",
       "\n",
       "   FP1024  \n",
       "0     0.0  \n",
       "1     1.0  \n",
       "\n",
       "[2 rows x 1026 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#select feature\n",
    "\n",
    "df1  = df1 [['chemblId', 'pIC50'\n",
    "             ,'FP135', 'FP97', 'FP933', 'FP639', 'FP75', 'FP28', 'FP239', 'FP571', 'FP404', 'FP373', 'FP137', 'FP243', 'FP362', 'FP676', 'FP343', 'FP133', 'FP237', 'FP78', 'FP206', 'FP385']]\n",
    "df2  = df2 [['chemblId', 'pIC50'\n",
    "             ,'ExtFP286', 'ExtFP914', 'ExtFP315', 'ExtFP767', 'ExtFP338', 'ExtFP343', 'ExtFP106', 'ExtFP111', 'ExtFP726', 'ExtFP711', 'ExtFP868', 'ExtFP811', 'ExtFP783', 'ExtFP423', 'ExtFP212', 'ExtFP718', 'ExtFP840', 'ExtFP467', 'ExtFP31', 'ExtFP560']]\n",
    "df3  = df3 [['chemblId', 'pIC50'\n",
    "             ,'EStateFP34', 'EStateFP30', 'EStateFP29', 'EStateFP35', 'EStateFP36', 'EStateFP37', 'EStateFP24', 'EStateFP50', 'EStateFP13', 'EStateFP16', 'EStateFP18', 'EStateFP51', 'EStateFP9', 'EStateFP32', 'EStateFP38', 'EStateFP7', 'EStateFP25', 'EStateFP19', 'EStateFP11', 'EStateFP54']]\n",
    "df4  = df4 [['chemblId', 'pIC50'\n",
    "             ,'GraphFP60', 'GraphFP382', 'GraphFP83', 'GraphFP21', 'GraphFP342', 'GraphFP228', 'GraphFP101', 'GraphFP140', 'GraphFP16', 'GraphFP224', 'GraphFP394', 'GraphFP163', 'GraphFP329', 'GraphFP72', 'GraphFP3', 'GraphFP499', 'GraphFP59', 'GraphFP270', 'GraphFP462', 'GraphFP89']]\n",
    "df5  = df5 [['chemblId', 'pIC50'\n",
    "             ,'MACCSFP86', 'MACCSFP139', 'MACCSFP92', 'MACCSFP97', 'MACCSFP77', 'MACCSFP36', 'MACCSFP131', 'MACCSFP75', 'MACCSFP154', 'MACCSFP118', 'MACCSFP113', 'MACCSFP109', 'MACCSFP80', 'MACCSFP133', 'MACCSFP95', 'MACCSFP136', 'MACCSFP62', 'MACCSFP98', 'MACCSFP127', 'MACCSFP72']]\n",
    "df6  = df6 [['chemblId', 'pIC50'\n",
    "             ,'PubchemFP193', 'PubchemFP308', 'PubchemFP375', 'PubchemFP439', 'PubchemFP696', 'PubchemFP391', 'PubchemFP420', 'PubchemFP534', 'PubchemFP199', 'PubchemFP493', 'PubchemFP345', 'PubchemFP186', 'PubchemFP566', 'PubchemFP580', 'PubchemFP258', 'PubchemFP405', 'PubchemFP346', 'PubchemFP681', 'PubchemFP33', 'PubchemFP757']]\n",
    "df7  = df7 [['chemblId', 'pIC50'\n",
    "             ,'SubFP23', 'SubFP169', 'SubFP88', 'SubFP18', 'SubFP181', 'SubFP137', 'SubFP20', 'SubFP2', 'SubFP184', 'SubFP5', 'SubFP183', 'SubFP1', 'SubFP135', 'SubFP182', 'SubFP172', 'SubFP49', 'SubFP3', 'SubFP275', 'SubFP100', 'SubFP33']]\n",
    "df8  = df8 [['chemblId', 'pIC50'\n",
    "             ,'SubFPC169', 'SubFPC23', 'SubFPC274', 'SubFPC302', 'SubFPC88', 'SubFPC181', 'SubFPC2', 'SubFPC18', 'SubFPC287', 'SubFPC137', 'SubFPC1', 'SubFPC20', 'SubFPC5', 'SubFPC172', 'SubFPC135', 'SubFPC182', 'SubFPC100', 'SubFPC3', 'SubFPC183', 'SubFPC33']]\n",
    "df9  = df9 [['chemblId', 'pIC50'\n",
    "             ,'KRFP1148', 'KRFP1642', 'KRFP1193', 'KRFP3560', 'KRFP433', 'KRFP4736', 'KRFP341', 'KRFP1645', 'KRFP2981', 'KRFP677', 'KRFP45', 'KRFP2242', 'KRFP2856', 'KRFP2', 'KRFP4820', 'KRFP2979', 'KRFP3788', 'KRFP2855', 'KRFP344', 'KRFP3681']]\n",
    "df10 = df10[['chemblId', 'pIC50'\n",
    "             ,'KRFPC1148', 'KRFPC1642', 'KRFPC297', 'KRFPC2711', 'KRFPC1193', 'KRFPC45', 'KRFPC3560', 'KRFPC444', 'KRFPC1645', 'KRFPC433', 'KRFPC1', 'KRFPC2981', 'KRFPC2262', 'KRFPC677', 'KRFPC4331', 'KRFPC2547', 'KRFPC669', 'KRFPC2856', 'KRFPC2979', 'KRFPC2438']]\n",
    "df11 = df11[['chemblId', 'pIC50'\n",
    "             ,'AD2D92', 'AD2D726', 'AD2D248', 'AD2D170', 'AD2D91', 'AD2D638', 'AD2D180', 'AD2D415', 'AD2D102', 'AD2D482', 'AD2D717', 'AD2D414', 'AD2D492', 'AD2D258', 'AD2D326', 'AD2D703', 'AD2D648', 'AD2D93', 'AD2D336', 'AD2D247']]\n",
    "df12 = df12[['chemblId', 'pIC50'\n",
    "             ,'APC2D1_C_C', 'APC2D2_N_O', 'APC2D10_O_O', 'APC2D1_C_N', 'APC2D1_C_O', 'APC2D4_N_O', 'APC2D2_N_N', 'APC2D3_N_O', 'APC2D9_N_O', 'APC2D7_O_O', 'APC2D6_O_S', 'APC2D3_O_O', 'APC2D5_N_O', 'APC2D7_N_O', 'APC2D5_O_O', 'APC2D10_N_S', 'APC2D6_O_O', 'APC2D1_C_S', 'APC2D4_O_O', 'APC2D2_N_S']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chemblId</th>\n",
       "      <th>pIC50</th>\n",
       "      <th>FP135</th>\n",
       "      <th>FP97</th>\n",
       "      <th>FP933</th>\n",
       "      <th>FP639</th>\n",
       "      <th>FP75</th>\n",
       "      <th>FP28</th>\n",
       "      <th>FP239</th>\n",
       "      <th>FP571</th>\n",
       "      <th>...</th>\n",
       "      <th>FP137</th>\n",
       "      <th>FP243</th>\n",
       "      <th>FP362</th>\n",
       "      <th>FP676</th>\n",
       "      <th>FP343</th>\n",
       "      <th>FP133</th>\n",
       "      <th>FP237</th>\n",
       "      <th>FP78</th>\n",
       "      <th>FP206</th>\n",
       "      <th>FP385</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHEMBL370037</td>\n",
       "      <td>8.221849</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CHEMBL189073</td>\n",
       "      <td>5.762708</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 22 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       chemblId     pIC50  FP135  FP97  FP933  FP639  FP75  FP28  FP239  \\\n",
       "0  CHEMBL370037  8.221849    1.0   1.0    1.0    0.0   0.0   0.0    1.0   \n",
       "1  CHEMBL189073  5.762708    1.0   1.0    1.0    0.0   0.0   0.0    1.0   \n",
       "\n",
       "   FP571  ...    FP137  FP243  FP362  FP676  FP343  FP133  FP237  FP78  FP206  \\\n",
       "0    0.0  ...      0.0    0.0    0.0    0.0    0.0    1.0    0.0   0.0    1.0   \n",
       "1    0.0  ...      0.0    0.0    0.0    0.0    0.0    1.0    0.0   0.0    1.0   \n",
       "\n",
       "   FP385  \n",
       "0    1.0  \n",
       "1    1.0  \n",
       "\n",
       "[2 rows x 22 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df1 .to_csv('QSAR_select/Fingerprinter.csv'               , sep=',' ,index=False)\n",
    "df2 .to_csv('QSAR_select/ExtendedFingerprinter.csv'       , sep=',' ,index=False)\n",
    "df3 .to_csv('QSAR_select/EStateFingerprinter.csv'         , sep=',' ,index=False)\n",
    "df4 .to_csv('QSAR_select/GraphOnlyFingerprinter.csv'      , sep=',' ,index=False)\n",
    "df5 .to_csv('QSAR_select/MACCSFingerprinter.csv'          , sep=',' ,index=False)\n",
    "df6 .to_csv('QSAR_select/PubchemFingerprinter.csv'        , sep=',' ,index=False)\n",
    "df7 .to_csv('QSAR_select/SubstructureFingerprinter.csv'   , sep=',' ,index=False)\n",
    "df8 .to_csv('QSAR_select/SubstructureFingerprintCount.csv', sep=',' ,index=False)\n",
    "df9 .to_csv('QSAR_select/KlekotaRothFingerprinter.csv'    , sep=',' ,index=False)\n",
    "df10.to_csv('QSAR_select/KlekotaRothFingerprintCount.csv' , sep=',' ,index=False)\n",
    "df11.to_csv('QSAR_select/AtomPairs2DFingerprinter.csv'    , sep=',' ,index=False)\n",
    "df12.to_csv('QSAR_select/AtomPairs2DFingerprintCount.csv' , sep=',' ,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(df1 ),len(df2 ),len(df3 ),len(df4 ),len(df5 )\\\n",
    ",len(df6 ),len(df7 ),len(df8 ),len(df9 ),len(df10)\\\n",
    ",len(df11),len(df12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zeromtmu/anaconda2/envs/ER_alpha/lib/python2.7/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.cross_validation import StratifiedShuffleSplit\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn import cross_validation\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from collections import defaultdict\n",
    "import os.path\n",
    "import shutil\n",
    "import pickle\n",
    "\n",
    "def build_model(X, Y, seed, hx, f):\n",
    "    \n",
    "    \n",
    "    #Data split using 70/30 ratio\n",
    "    X_internal, X_external, Y_internal, Y_external = train_test_split(X,Y,\n",
    "                                                        test_size=0.3,\n",
    "                                                        random_state=seed)\n",
    "\n",
    "    # Training set\n",
    "    rf = RandomForestRegressor(n_estimators=300, \n",
    "                               max_features='sqrt', \n",
    "                               min_samples_leaf = 1,\n",
    "                               random_state = 13,\n",
    "                               n_jobs=-1)\n",
    "    \n",
    "    rf.fit(X_internal,Y_internal)\n",
    "    prediction = rf.predict(X_internal)\n",
    "    \n",
    "    # save the model to disk\n",
    "    filename = f + '_model' + str(seed) + '.sav'\n",
    "    pickle.dump(rf, open(filename, 'wb'))\n",
    "          \n",
    "    # Cross-validation\n",
    "    cv = cross_validation.cross_val_predict(rf, X_internal, Y_internal, cv=10, n_jobs=-1)\n",
    "    \n",
    "    # External set  \n",
    "    prediction_external = rf.predict(X_external)\n",
    " \n",
    "    #print result from each seed    \n",
    "    R2_train.append(r2_score(Y_internal, prediction))\n",
    "    RMSE_train.append((mean_squared_error(Y_internal, prediction)))\n",
    "    Q2_CV.append(r2_score(Y_internal, cv))\n",
    "    RMSE_CV.append((mean_squared_error(Y_internal, cv)))\n",
    "    Q2_External.append(r2_score(Y_external, prediction_external))\n",
    "    RMSE_External.append((mean_squared_error(Y_external, prediction_external)))\n",
    "    \n",
    "    #Feature Importance\n",
    "    Feature = hx[:]\n",
    "    feature_importance = rf.feature_importances_\n",
    "    importances = 100.0 * (feature_importance / feature_importance.max()) #index\n",
    "    \n",
    "    for i, fx in enumerate(Feature):\n",
    "        importances_dict[fx].append(importances[i])\n",
    "    \n",
    "    return R2_train, RMSE_train, Q2_CV, RMSE_CV, Q2_External, RMSE_External, Feature, \\\n",
    "           X_internal, X_external, Y_internal, Y_external, rf, prediction, importances_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "\n",
    "def Y_scrambling(X_internal, X_external, Y_internal, Y_external):\n",
    "    # Do the Y-scrambling. Loop over the actual learning for 100 times.\n",
    "    for randomseedcounter in range(1,101):\n",
    "        y_train_scrambled = deepcopy(Y_internal)\n",
    "        X_train_scrambled = deepcopy(X_internal)\n",
    "        np.random.shuffle(y_train_scrambled)\n",
    "        np.random.shuffle(X_train_scrambled)\n",
    "\n",
    "        # training was done on \"scrambled\" data - prediction on test set\n",
    "        RF_scrambled         = RandomForestRegressor()\n",
    "        RF_scrambled         = RF_scrambled.fit(X_internal,y_train_scrambled)\n",
    "        y_predict_scrambled  = RF_scrambled.predict(X_external)\n",
    "    \n",
    "        acclist_predictionOnTest_scrambledtrain.append((RF_scrambled.score(X_external,Y_external))**2)\n",
    "    \n",
    "        # training was done on \"scrambled\" data - prediction on train set\n",
    "        y_predict_scrambled_predictTrain  = RF_scrambled.predict(X_internal)\n",
    "    \n",
    "        acclist_predictionOnTrain_scrambledtrain.append((RF_scrambled.score(X_internal,Y_internal))**2)\n",
    "        \n",
    "        r2 = pd.DataFrame(acclist_predictionOnTrain_scrambledtrain, columns=['R2'])\n",
    "        q2 = pd.DataFrame(acclist_predictionOnTest_scrambledtrain, columns=['Q2'])        \n",
    "        \n",
    "    return acclist_predictionOnTest_scrambledtrain,acclist_predictionOnTrain_scrambledtrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mean(R2_train, RMSE_train, Q2_CV, RMSE_CV, Q2_External, RMSE_External, importances_dict):\n",
    "    R2_train_mean = np.mean(R2_train)\n",
    "    RMSE_train_mean = np.mean(RMSE_train)\n",
    "    Q2_CV_mean = np.mean(Q2_CV)\n",
    "    RMSE_CV_mean = np.mean(RMSE_CV)\n",
    "    Q2_External_mean = np.mean(Q2_External)\n",
    "    RMSE_External_mean = np.mean(RMSE_External)\n",
    "    importances_mean0 = {}\n",
    "    for fx in importances_dict:\n",
    "        importances_mean0[fx] = np.mean(importances_dict[fx])\n",
    "    importances_mean = sorted([(k,v) for k,v in importances_mean0.iteritems()],\n",
    "                                    key=lambda x: x[1], reverse=True)\n",
    "    \n",
    "    return R2_train_mean, RMSE_train_mean, Q2_CV_mean, RMSE_CV_mean, Q2_External_mean, \\\n",
    "            RMSE_External_mean, importances_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def std(R2_train, RMSE_train, Q2_CV, RMSE_CV, Q2_External, RMSE_External, importances_dict):\n",
    "    R2_train_std = np.std(R2_train)\n",
    "    RMSE_train_std = np.std(RMSE_train)\n",
    "    Q2_CV_std = np.std(Q2_CV)\n",
    "    RMSE_CV_std = np.std(RMSE_CV)\n",
    "    Q2_External_std = np.std(Q2_External)\n",
    "    RMSE_External_std = np.std(RMSE_External)\n",
    "    importances_std0 = {}\n",
    "    for fx in importances_dict:\n",
    "        importances_std0[fx] = np.std(importances_dict[fx])\n",
    "    importances_std = sorted([(k,v) for k,v in importances_std0.iteritems()],\n",
    "                                    key=lambda x: x[1], reverse=True)\n",
    "  \n",
    "    return R2_train_std, RMSE_train_std, Q2_CV_std, RMSE_CV_std, Q2_External_std,\\\n",
    "            RMSE_External_std, importances_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def print_output(R2_train_mean, RMSE_train_mean, Q2_CV_mean, RMSE_CV_mean, Q2_External_mean, \n",
    "           RMSE_External_mean, R2_train_std, RMSE_train_std, Q2_CV_std, RMSE_CV_std, Q2_External_std, \n",
    "           RMSE_External_std, X_internal, X_external, X, f ):\n",
    "    #outfile = open('all_output.csv', 'a')\n",
    "    print >> outfile, '%s,%d,%d,%0.4f,%0.4f,%0.4f,%0.4f,%d,%0.4f,%0.4f,%0.4f,%0.4f,%d,%0.4f,%0.4f,%0.4f,%0.4f' \\\n",
    "                                     % (f, len(X_internal),\n",
    "                                        len(X[0]),\n",
    "                                        R2_train_mean, R2_train_std,\n",
    "                                        RMSE_train_mean, RMSE_train_std,\n",
    "                                        len(X_internal),\n",
    "                                        Q2_CV_mean, Q2_CV_std,\n",
    "                                        RMSE_CV_mean, RMSE_CV_std,\n",
    "                                        len(X_external),\n",
    "                                        Q2_External_mean, Q2_External_std,\n",
    "                                        RMSE_External_mean, RMSE_External_std,\n",
    "                                       )  \n",
    "        \n",
    "        \n",
    "    \n",
    "    print '\\nTraining set\\n------------'\n",
    "    print 'N: ' + (str(len(X_internal)))\n",
    "    print 'R2: %0.4f'%(R2_train_mean)\n",
    "    print 'std_R2: %0.4f'%(R2_train_std)\n",
    "    print 'RMSE: %0.4f'%(RMSE_train_mean)\n",
    "    print 'std_RMSE: %0.4f'%(RMSE_train_std)\n",
    "\n",
    "    print '\\nCross-validation set\\n------------'\n",
    "    print 'N: ' + (str(len(X_internal)))\n",
    "    print 'Q2: %0.4f'%(Q2_CV_mean)\n",
    "    print 'std_Q2: %0.4f'%(Q2_CV_std)\n",
    "    print 'RMSE: %0.4f'%(RMSE_CV_mean)\n",
    "    print 'std_RMSE: %0.4f'%(RMSE_CV_std)\n",
    "\n",
    "    print '\\nExternal set\\n------------'\n",
    "    print 'N: ' + (str(len(X_external)))\n",
    "    print 'Q2_EXt: %0.4f'%(Q2_External_mean)\n",
    "    print 'std_Q2_EXt: %0.4f'%(Q2_External_std)\n",
    "    print 'RMSE: %0.4f'%(RMSE_External_mean)\n",
    "    print 'std_RMSE: %0.4f'%(RMSE_External_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zeromtmu/anaconda2/envs/ER_alpha/lib/python2.7/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "import pylab as py \n",
    "import statsmodels.api as sm\n",
    "from statsmodels.stats.outliers_influence import summary_table\n",
    "\n",
    "def plot_model (f, X_internal, X_external, Y_internal, Y_external,\n",
    "                R2_train_mean, Q2_External_mean,\n",
    "                importances_mean, importances_std, Feature, prediction, \n",
    "                acclist_predictionOnTest_scrambledtrain, acclist_predictionOnTrain_scrambledtrain):\n",
    "    \n",
    "    # Prepare plot\n",
    "    m  = rf.fit(X_internal,Y_internal)\n",
    "    cm = plt.cm.RdBu\n",
    "    cv = cross_validation.cross_val_predict(rf, X_internal, Y_internal, cv=10, n_jobs=-1)\n",
    "\n",
    "    fig_size = plt.rcParams[\"figure.figsize\"]\n",
    "    fig_size[0]=5\n",
    "    fig_size[1]=5\n",
    "\n",
    "    # Train Set\n",
    "    x_train = np.array(Y_internal)\n",
    "    y_train = m.predict(X_internal).flatten()\n",
    "    py.scatter(x_train, y_train, s=50, marker='.', alpha=0.3,\n",
    "            c='g', cmap=cm ,edgecolors='g')\n",
    "                 #label=r\"$R^{2}_{Tr}$ = %.4f\" % R2_train_mean)\n",
    "    \n",
    "    \n",
    "    # CV Set\n",
    "    np.array(cv)\n",
    "    x_test = np.array(Y_internal)\n",
    "    y_test = cv\n",
    "    py.scatter(x_test, y_test, s=50, marker='.', alpha=0.3,\n",
    "            c='b', cmap=cm ,edgecolors='b') \n",
    "                 #label=r\"$Q^{2}_{Ext}$ = %.4f\" % Q2_External_mean)\n",
    "        #2SD line\n",
    "    X = sm.add_constant(x_test)\n",
    "    res = sm.OLS(y_test, X).fit()\n",
    "\n",
    "    st, data, ss2 = summary_table(res, alpha=0.05)\n",
    "    fittedvalues = data[:,2]\n",
    "        #predict_mean_se  = data[:,3]\n",
    "    predict_mean_ci_low, predict_mean_ci_upp = data[:,4:6].T\n",
    "    predict_ci_low, predict_ci_upp = data[:,6:8].T\n",
    "        #2SD line\n",
    "    plt.plot(X, predict_ci_low, '--b', linewidth=0.5, alpha=0.5)\n",
    "    plt.plot(X, predict_ci_upp, '--b', linewidth=0.5, alpha=0.5)\n",
    "\n",
    "        \n",
    "    # External Set\n",
    "    x_test = np.array(Y_external)\n",
    "    y_test = m.predict(X_external).flatten()\n",
    "    py.scatter(x_test, y_test, s=50, marker='.', alpha=0.3,\n",
    "            c='r', cmap=cm ,edgecolors='r') \n",
    "    \n",
    "    plt.legend(loc=2,prop={'size':6})\n",
    "    plt.xlabel(\"Experimental $pIC_{50}$ values\", fontsize=10)\n",
    "    plt.ylabel(\"Predicted $pIC_{50}$ values\", fontsize=10)\n",
    "        \n",
    "    min_axis = np.min(np.concatenate([Y_internal, prediction], axis=0))\n",
    "    max_axis = np.max(np.concatenate([Y_internal, prediction], axis=0))\n",
    "    plt.xlim([(min_axis*0.9),(max_axis*1.05)])\n",
    "    plt.ylim([(min_axis*0.9),(max_axis*1.05)])\n",
    "    plt.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "    \n",
    "    # Save plot to file\n",
    "    plt.savefig(f+'_compair_Fp.pdf', dpi=300)\n",
    "    plt.show()\n",
    "    \n",
    "    # Y-scrambling plot\n",
    "    py.scatter(Q2_CV_mean, R2_train_mean, s=100, marker='.', alpha=0.3,\n",
    "            c='b', cmap=cm ,edgecolors='b') \n",
    "            \n",
    "    py.scatter(acclist_predictionOnTest_scrambledtrain,acclist_predictionOnTrain_scrambledtrain, \\\n",
    "               s=100, marker='.', alpha=0.3, c='r', cmap=cm ,edgecolors='r')\n",
    "\n",
    "    plt.legend(loc=2,prop={'size':6})\n",
    "    plt.xlabel(\"$Q^{2}$\", fontsize=10)\n",
    "    plt.ylabel(\"$R^{2}$\", fontsize=10)\n",
    "    plt.tick_params(axis='both', which='major', labelsize=14)\n",
    "    \n",
    "    plt.ylim([0, 1])\n",
    "    plt.xlim([0, 1])\n",
    "    plt.axhline(0.5, color='gray',linestyle='--', dashes=(5, 10), linewidth=0.5)\n",
    "    plt.axvline(0.5, color='gray',linestyle='--', dashes=(5, 10), linewidth=0.5)\n",
    "    \n",
    "    plt.savefig(f+'_Y-scrambling.pdf', dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "    #Feature Importance \n",
    "    fig_size = plt.rcParams[\"figure.figsize\"]\n",
    "    fig_size[0]=5\n",
    "    fig_size[1]=10\n",
    "\n",
    "    barlist = plt.barh(range(20), [x[1] for x in importances_mean[:20]], \n",
    "          color=\"g\", align=\"center\") # , xerr=[x[1] for x in importances_std[:20]], error_kw=dict(ecolor='gray', lw=2, capsize=5, capthick=2)\n",
    "\n",
    "    \n",
    "    plt.yticks(range(20), [x[0] for x in importances_mean[:20]])\n",
    "    plt.ylim([-1, 20])\n",
    "    \n",
    "    plt.xlabel(r\"$\\bf{Gini}$\" + \" \"+ r\"$\\bf{index}$\", fontsize=12)\n",
    "    ax = plt.gca()\n",
    "    ax.invert_yaxis()\n",
    "    plt.tight_layout(pad=2.0, w_pad=0.7, h_pad=2.0)\n",
    "    plt.savefig(f+'_Feature_importances.pdf', dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constructing the QSAR model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running....QSAR_select/AtomPairs2DFingerprintCount.csv\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zeromtmu/anaconda2/envs/ER_alpha/lib/python2.7/site-packages/ipykernel_launcher.py:25: DeprecationWarning: \n",
      ".ix is deprecated. Please use\n",
      ".loc for label based indexing or\n",
      ".iloc for positional indexing\n",
      "\n",
      "See the documentation here:\n",
      "http://pandas.pydata.org/pandas-docs/stable/indexing.html#deprecate_ix\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.9158\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.1599\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.7001\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.5696\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7219\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5326\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/SubstructureFingerprintCount.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.9340\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.1254\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.7011\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.5677\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7374\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5029\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/SubstructureFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8501\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.2846\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6585\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6486\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.6853\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.6027\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/PubchemFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8127\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3558\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6673\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6318\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7217\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5329\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/GraphOnlyFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.7961\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3872\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6455\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6733\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.6298\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.7091\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/KlekotaRothFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8077\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3651\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6422\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6796\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.6842\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.6049\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/ExtendedFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8252\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3319\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6575\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6505\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.6701\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.6318\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/AtomPairs2DFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8398\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3042\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6539\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6573\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.6735\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.6252\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/MACCSFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8391\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3056\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.7022\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.5655\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7063\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5624\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/KlekotaRothFingerprintCount.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.9017\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.1866\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6780\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6115\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7125\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5507\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/EStateFingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.8500\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.2849\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6460\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.6724\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7138\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5481\n",
      "std_RMSE: 0.0000\n",
      "running....QSAR_select/Fingerprinter.csv\n",
      "\n",
      "Training set\n",
      "------------\n",
      "N: 861\n",
      "R2: 0.7924\n",
      "std_R2: 0.0000\n",
      "RMSE: 0.3943\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "Cross-validation set\n",
      "------------\n",
      "N: 861\n",
      "Q2: 0.6281\n",
      "std_Q2: 0.0000\n",
      "RMSE: 0.7063\n",
      "std_RMSE: 0.0000\n",
      "\n",
      "External set\n",
      "------------\n",
      "N: 370\n",
      "Q2_EXt: 0.7001\n",
      "std_Q2_EXt: 0.0000\n",
      "RMSE: 0.5743\n",
      "std_RMSE: 0.0000\n"
     ]
    }
   ],
   "source": [
    "import glob #to read multiple files \n",
    "from rdkit import Chem, DataStructs\n",
    "from rdkit.Chem import AllChem\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "\n",
    "! rm Result/CHEMBL206_IC50_Revision2.csv\n",
    "\n",
    "outfile = open('Result/CHEMBL206_IC50_Revision2.csv', 'a')\n",
    "\n",
    "print >> outfile, 'Filename,N_train,Descriptors,R2_train,R2_train_std,' + \\\n",
    "                    'RMSE_train,RMSE_train_std,N_CV,Q2_CV,Q2_CV_std,RMSE_CV,' + \\\n",
    "                    'RMSE_CV_std,N_External,Q2_External,Q2_External_std,RMSE_External,RMSE_External_std'\n",
    "\n",
    "path = r'QSAR_select/'\n",
    "\n",
    "for f in glob.glob(path + '*.csv'):\n",
    "    print 'running....' + f\n",
    "    \n",
    "    df = pd.read_csv(f)\n",
    "    df = df.apply(lambda x: pd.to_numeric(x,errors='ignore'))\n",
    "    df = df.fillna(method='ffill')\n",
    "    \n",
    "    Y = df[\"pIC50\"].tolist()\n",
    "    \n",
    "    X = df.ix[:,2:]\n",
    "    hx = X.columns.tolist()\n",
    "    data = X.as_matrix().astype(np.float)\n",
    "    X = np.array(data)    \n",
    "    \n",
    "    # Prepare empty lists to plot QSAR model\n",
    "    R2_train = []\n",
    "    RMSE_train = []\n",
    "    Q2_CV = []\n",
    "    RMSE_CV = []\n",
    "    Q2_External = []\n",
    "    RMSE_External = []\n",
    "    importances_dict = defaultdict(list)\n",
    "    \n",
    "    # Prepare empty lists to plot the performance of accuracy.\n",
    "    acclist_realRF                          = []\n",
    "    acclist_realRF_predictTrain             = []\n",
    "    acclist_predictionOnTest_scrambledtrain  = []\n",
    "    acclist_predictionOnTrain_scrambledtrain = []\n",
    "    \n",
    "    seed = 0\n",
    "    i = seed \n",
    "    R2_train, RMSE_train, Q2_CV, RMSE_CV, Q2_External, RMSE_External, Feature, \\\n",
    "           X_internal, X_external, Y_internal, Y_external, rf, prediction, importances_dict = build_model(X, Y, i, hx, f)\n",
    "            \n",
    "    acclist_predictionOnTest_scrambledtrain, acclist_predictionOnTrain_scrambledtrain = Y_scrambling( \\\n",
    "                                                                        X_internal, X_external, Y_internal, Y_external)\n",
    "    R2_train_mean, RMSE_train_mean, Q2_CV_mean, RMSE_CV_mean, Q2_External_mean, \\\n",
    "        RMSE_External_mean, importances_mean = mean(R2_train, RMSE_train, Q2_CV, RMSE_CV, \\\n",
    "                                 Q2_External, RMSE_External, importances_dict)\n",
    "    R2_train_std, RMSE_train_std, Q2_CV_std, RMSE_CV_std, Q2_External_std,\\\n",
    "        RMSE_External_std, importances_std = std(R2_train, RMSE_train, Q2_CV, RMSE_CV, \\\n",
    "                                 Q2_External, RMSE_External, importances_dict)\n",
    "    print_output(R2_train_mean, RMSE_train_mean, Q2_CV_mean, RMSE_CV_mean, Q2_External_mean, \n",
    "           RMSE_External_mean, R2_train_std, RMSE_train_std, Q2_CV_std, RMSE_CV_std, Q2_External_std, \n",
    "           RMSE_External_std, X_internal, X_external, X, f)\n",
    "    \n",
    "outfile.close()\n",
    "    #END"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}