generate_survey_materials.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49efb72c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn import preprocessing, svm, metrics, tree, decomposition, svm\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier, GradientBoostingClassifier, AdaBoostClassifier\n",
    "from sklearn.linear_model import LogisticRegression, LinearRegression, Lasso, Ridge, Perceptron, SGDClassifier, OrthogonalMatchingPursuit\n",
    "from sklearn.naive_bayes import GaussianNB, MultinomialNB, BernoulliNB\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.metrics import plot_confusion_matrix, accuracy_score, roc_auc_score, precision_score, recall_score, f1_score, roc_curve, auc\n",
    "from sklearn.model_selection import train_test_split, ParameterGrid, KFold, StratifiedKFold\n",
    "from sklearn.dummy import DummyClassifier\n",
    "from sklearn import tree\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import yaml\n",
    "import pickle\n",
    "from statistics import mean\n",
    "import random\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "import shap\n",
    "\n",
    "pd.set_option('display.max_colwidth', None)\n",
    "pd.set_option('display.max_columns', None)\n",
    "pd.set_option('display.max_rows', None)\n",
    "\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc0b992d",
   "metadata": {},
   "source": [
    "## Inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e17106e",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = pd.read_csv(\"matrices/por_X.csv\",index_col=0)\n",
    "y = pd.read_csv(\"matrices/por_y.csv\",index_col=0,squeeze = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e18cce0",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1eb26833",
   "metadata": {},
   "source": [
    "# Tree-based Models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2f9ab5",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b46bb2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_params = {'max_depth': 5, 'max_features': 'sqrt', 'min_samples_split': 2, 'n_estimators': 100, 'n_jobs': -1}\t\n",
    "clf = RandomForestClassifier(**clf_params,random_state=42).fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2146000a",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_params = {'max_depth': 3}\n",
    "clf = DecisionTreeClassifier(**clf_params,random_state=42).fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7070f0a0",
   "metadata": {},
   "source": [
    "## Feature importances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a685a4a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "explainer = shap.TreeExplainer(clf)\n",
    "shap_values = explainer.shap_values(X_test)\n",
    "shap.summary_plot(shap_values[1], X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef6fe5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows = 1,ncols = 1,figsize = (4,4), dpi=300)\n",
    "\n",
    "tree.plot_tree(clf,\n",
    "               feature_names = X_test.columns, \n",
    "               class_names=['Pass','Fail'],\n",
    "               filled = True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d5942b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = pd.DataFrame()\n",
    "temp['names'] = X_test.columns\n",
    "a = np.array(shap_values[1])\n",
    "a = np.abs(a)\n",
    "temp['abs_vals'] = np.mean(a, axis=0)\n",
    "temp = temp.sort_values(by=['abs_vals'],ascending=False)\n",
    "temp.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eed88aeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "top_features = temp['names'].to_list()[0:5]\n",
    "top_features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "863ad3b9",
   "metadata": {},
   "source": [
    "## Get profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f9205f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_random_index(np_arr):\n",
    "    passed = 0\n",
    "    while passed == 0:\n",
    "        i = random.randint(0,len(np_arr)-1)\n",
    "        if np_arr[i] == 1:\n",
    "            return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a4484df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get profiles\n",
    "\n",
    "# Get threshold\n",
    "y_pred_proba = clf.predict_proba(X_test)[:,1]\n",
    "k = 25\n",
    "threshold = np.percentile(y_pred_proba,(100-k))\n",
    "print(threshold)\n",
    "threshold_high = np.percentile(y_pred_proba,90)\n",
    "print(threshold_high)\n",
    "threshold_low = np.percentile(y_pred_proba,1)\n",
    "print(threshold_low)\n",
    "\n",
    "\n",
    "# Get index of 2 profiles close to threshold\n",
    "next_score = min(y_pred_proba[y_pred_proba > threshold])\n",
    "print(next_score)\n",
    "arr = np.where((y_pred_proba == next_score), 1, 0)\n",
    "i_above_threshold = get_random_index(arr)\n",
    "\n",
    "prev_score = max(y_pred_proba[y_pred_proba < threshold])\n",
    "print(prev_score)\n",
    "arr = np.where((y_pred_proba == prev_score), 1, 0)\n",
    "i_below_threshold = get_random_index(arr)\n",
    "\n",
    "# Get index of low profile\n",
    "arr = np.where((y_pred_proba <= threshold_low), 1, 0)\n",
    "i_low = get_random_index(arr)\n",
    "\n",
    "# Get index of high profile\n",
    "arr = np.where((y_pred_proba >= threshold_high), 1, 0)\n",
    "i_high = get_random_index(arr)\n",
    "\n",
    "print(i_above_threshold)\n",
    "print(i_below_threshold)\n",
    "print(i_low)\n",
    "print(i_high)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f68aff",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Slightly above threshold | {i_above_threshold}\")\n",
    "X_test.filter(items=top_features).iloc[i_above_threshold]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21ab2213",
   "metadata": {},
   "outputs": [],
   "source": [
    "feat_names_long = ['Number of absences',\n",
    " 'Number of failures',\n",
    " 'First year grade',\n",
    " 'traveltime',\n",
    " 'Study time',\n",
    " 'Quality of family relationships',\n",
    " 'freetime',\n",
    " 'Frequency of going out with friends',\n",
    " 'Workday alcohol consumption',\n",
    " 'Weekly alcoholic drinks',\n",
    " 'Health status',\n",
    " 'activities',\n",
    " 'nursery',\n",
    " 'Has internet',\n",
    " 'Has romantic relationship',\n",
    " 'Attends school Gabriel Pereira',\n",
    " 'address_urban',\n",
    " 'Family size greater than 3',\n",
    " 'Pstatus_togethor',\n",
    " 'sex_female',\n",
    " 'health_mjob',\n",
    " 'other_mjob',\n",
    " 'services_mjob',\n",
    " 'teacher_mjob',\n",
    " 'health_fjob',\n",
    " \"Father's job not listed on survey\",\n",
    " \"Father's job in services\",\n",
    " \"Father's job is teacher\",\n",
    " 'home_reason',\n",
    " 'other_reason',\n",
    " 'School chosen based on reputation',\n",
    " 'mother_guardian',\n",
    " 'Guardian not mother or father']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5efac8b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.waterfall_plot(shap.Explanation(values=shap_values[1][i_above_threshold], \n",
    "                                     base_values=explainer.expected_value[1], \n",
    "                                     data=X_test.iloc[i_above_threshold],\n",
    "                                     feature_names=feat_names_long),max_display=(len(top_features)+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68002c0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Slightly below threshold | ID {i_below_threshold}\")\n",
    "X_test.filter(items=top_features).iloc[i_below_threshold]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8331422",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.waterfall_plot(shap.Explanation(values=shap_values[1][i_below_threshold], \n",
    "                                     base_values=explainer.expected_value[1], \n",
    "                                     data=X_test.iloc[i_below_threshold],\n",
    "                                     feature_names=feat_names_long),max_display=(len(top_features)+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef1a6a2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"High likelihood of failure | ID {i_high}\")\n",
    "X_test.filter(items=top_features).iloc[i_high]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06e20578",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.waterfall_plot(shap.Explanation(values=shap_values[1][i_high], \n",
    "                                     base_values=explainer.expected_value[1], \n",
    "                                     data=X_test.iloc[i_high],\n",
    "                                     feature_names=feat_names_long),max_display=(len(top_features)+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39d52122",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"High likelihood of passing | ID {i_low}\")\n",
    "X_test.filter(items=top_features).iloc[i_low]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1f7b3a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.waterfall_plot(shap.Explanation(values=shap_values[1][i_low], \n",
    "                                     base_values=explainer.expected_value[1], \n",
    "                                     data=X_test.iloc[i_low],\n",
    "                                     feature_names=feat_names_long),max_display=(len(top_features)+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf476218",
   "metadata": {},
   "source": [
    "# Linear model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2be855f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_params = {'alpha': 0.01}\n",
    "clf = Lasso(**clf_params,random_state=42).fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be61c5b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f05f1c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = pd.DataFrame()\n",
    "temp['coef'] = clf.coef_\n",
    "temp['names'] = X_train.columns\n",
    "temp['coef_abs'] = abs(clf.coef_)\n",
    "temp = temp.sort_values(by=['coef_abs'],ascending=False)\n",
    "temp = temp[temp.coef != 0]\n",
    "top_features = temp['names'].to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aa3c7c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp.head(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b38e055",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.barplot(x='coef',y='names',data=temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18f653db",
   "metadata": {},
   "source": [
    "## Get profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db93c1f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get profiles\n",
    "\n",
    "# Get threshold\n",
    "y_pred_proba = clf.predict(X_test)\n",
    "k = 25\n",
    "threshold = np.percentile(y_pred_proba,(100-k))\n",
    "print(threshold)\n",
    "threshold_high = np.percentile(y_pred_proba,90)\n",
    "print(threshold_high)\n",
    "threshold_low = np.percentile(y_pred_proba,1)\n",
    "print(threshold_low)\n",
    "\n",
    "\n",
    "# Get index of 2 profiles close to threshold\n",
    "next_score = min(y_pred_proba[y_pred_proba > threshold])\n",
    "print(next_score)\n",
    "arr = np.where((y_pred_proba == next_score), 1, 0)\n",
    "i_above_threshold = get_random_index(arr)\n",
    "\n",
    "prev_score = max(y_pred_proba[y_pred_proba < threshold])\n",
    "print(prev_score)\n",
    "arr = np.where((y_pred_proba == prev_score), 1, 0)\n",
    "i_below_threshold = get_random_index(arr)\n",
    "\n",
    "# Get index of low profile\n",
    "arr = np.where((y_pred_proba <= threshold_low), 1, 0)\n",
    "i_low = get_random_index(arr)\n",
    "\n",
    "# Get index of high profile\n",
    "arr = np.where((y_pred_proba >= threshold_high), 1, 0)\n",
    "i_high = get_random_index(arr)\n",
    "\n",
    "print(i_above_threshold)\n",
    "print(i_below_threshold)\n",
    "print(i_low)\n",
    "print(i_high)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "533e0e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Slightly above threshold | {i_above_threshold}\")\n",
    "X_test.filter(items=top_features).iloc[i_above_threshold]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ecc2455",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Slightly below threshold | ID {i_below_threshold}\")\n",
    "X_test.filter(items=top_features).iloc[i_below_threshold]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f19feef",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"High likelihood of failure | ID {i_high}\")\n",
    "X_test.filter(items=top_features).iloc[i_high]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b3df8d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"High likelihood of failure | ID {i_low}\")\n",
    "X_test.filter(items=top_features).iloc[i_low]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f73ad5d6",
   "metadata": {},
   "source": [
    "## Deprecated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25b1e2cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# profiles:\n",
    "# i_above_threshold = 28\n",
    "# i_below_threshold = 142\n",
    "# i_low = 14\n",
    "# i_high = 190"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7132aeec",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_proba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "405a4464",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test.iloc[i_below_threshold]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5f015c",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = get_random_index(arr)\n",
    "print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30e0d457",
   "metadata": {},
   "outputs": [],
   "source": [
    "explainer = shap.TreeExplainer(rf)\n",
    "shap_values = explainer.shap_values(X_test)\n",
    "shap.summary_plot(shap_values[1], X_test, plot_type='bar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3e9eb7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = pd.DataFrame()\n",
    "temp['names'] = X_test.columns\n",
    "a = np.array(shap_values[1])\n",
    "a = np.abs(a)\n",
    "temp['abs_vals'] = np.mean(a, axis=0)\n",
    "temp = temp.sort_values(by=['abs_vals'],ascending=False)\n",
    "temp.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "563bd888",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "y_pred_proba = rf.predict_proba(X_test)[:,1]\n",
    "k = 25\n",
    "y_pred = np.where(y_pred_proba > np.percentile(y_pred_proba,(100-k)), 1, 0)\n",
    "i = np.where(y_pred == 1)\n",
    "\n",
    "temp = pd.DataFrame()\n",
    "temp['names'] = X_train.columns\n",
    "a = np.array(shap_values[1][i])\n",
    "temp['vals'] = np.mean(a, axis=0)\n",
    "temp['abs_vals'] = temp['vals'].abs()\n",
    "temp = temp.sort_values(by=['abs_vals'],ascending=False)\n",
    "temp.head()\n",
    "\n",
    "sns.barplot(x='vals', y='names',data=temp[0:10])\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79e553aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.waterfall_plot(shap.Explanation(values=shap_values[1][i_low], \n",
    "                                     base_values=explainer.expected_value[1], \n",
    "                                     data=X_test.iloc[i_below_threshold],\n",
    "                                     feature_names=X_train.columns.tolist()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f66b9fc4",
   "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.9.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}