logistic-regression-titanic-mit-shap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistische Regression auf dem Titanic Datensatz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hole Daten\n",
    "\n",
    "Der Datensatz stammt von [Kaggle](https://www.kaggle.com/sureshbhusare/titanic-dataset-from-kaggle).\n",
    "\n",
    "Er beinhaltet Daten zu den über 1.300 Passagieren der tragischen Jungfernfahrt.\n",
    "\n",
    "Die Daten der über 900 Mann starken Besatzung sind nicht enthalten.\n",
    "\n",
    "Die Spalten haben folgende Bedeutung:\n",
    "\n",
    "Variable\t|Definition\t|Key\n",
    ":---|:---|---\n",
    "Survival\t|Survival - Label\t|0 = No, 1 = Yes\n",
    "Pclass\t|Ticket class\t|1 = 1st, 2 = 2nd, 3 = 3rd\n",
    "Sex\t|Sex\t| male, female\n",
    "Age\t|Age in years|\t\n",
    "SibSp\t|# of siblings / spouses aboard the Titanic\t|\n",
    "Parch\t| # of parents / children aboard the Titanic\t|\n",
    "Ticket\t|Ticket number\t|\n",
    "Fare\t|Passenger fare\t|\n",
    "Cabin\t|Cabin number\t|\n",
    "Embarked\t|Port of Embarkation\t|C = Cherbourg, Q = Queenstown, S = Southampton\n",
    "Name|Name des Passagiers |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = pd.read_csv('./data/titanic/train.csv')\n",
    "test = pd.read_csv('./data/titanic/test.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ein erster Blick auf die Daten\n",
    "\n",
    "`shape` sagt uns die Dimensionen des Datensatzes: Wieviele Datensätze (Zeilen) und wieviele Features (Spalten) haben wir?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(train.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mit `head` oder `tail` können wir uns die ersten bzw. letzten paar Datensätze anschauen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`info` gibt uns ein paar Informationen zu den Datentypen der Spalten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.info();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`describe` liefert uns einige Kennzahlen zur statistischen Verteilung der Daten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.describe(exclude='O')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir konnten aus `train.info()` schon ablesen, dass zu \"Age\", \"Embarked\" und \"Cabin\" weniger non-null Werte als 891 angezeigt werden. Zu diesen Null-Werten müssen wir noch überlegen, wie wir damit umgehen.\n",
    "\n",
    "Hier können wir uns ein Bild davon machen, wo die Null-Werte auftreten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(data = train.isnull(), yticklabels=False, cbar=False, cmap='viridis');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploratory Data Analysis\n",
    "\n",
    "... wie sehen unsere Daten eigentlich aus?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotten wir mal eine Verteilung:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "f, axes = plt.subplots(2, 4, figsize=(20, 7), sharex=False)\n",
    "\n",
    "showColumns = ['Survived', 'Pclass', 'Sex', 'Age', 'SibSp', 'Parch', 'Fare', 'Embarked']\n",
    "\n",
    "r = 0\n",
    "c = 0\n",
    "for col in showColumns:\n",
    "    sns.countplot(train[col], ax=axes[r, c])\n",
    "    c = c + 1\n",
    "    if (c > 3):\n",
    "        r = 1\n",
    "        c = 0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Fare\" sieht seltsam aus, schauen wir uns mal näher an:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(1, 2, figsize=(20, 7), sharex=False)\n",
    "sns.distplot(train['Fare'], ax=axes[0])\n",
    "sns.countplot(train['Fare'], ax=axes[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation\n",
    "\n",
    "In der Statistik misst der Korrelationskoeffizient $\\rho$ die Stärke und Richtung einer linearen Beziehung zwischen zwei Variablen.\n",
    "\n",
    "Ab einem Absolutwert größer 0.5 sollte man die Korrelation betrachten."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "sns.heatmap(train.corr(), annot=True, center=0, square=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Daten Aufbereitung\n",
    "\n",
    "### Fehlende Werte: Cabin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Spalte \"Cabin\" ist nur zu einem geringen Teil gefüllt und es scheint nicht plausibel, dass aus der Kabinenbezeichnung \n",
    "ein Einfluss auf die Überlebenswahrscheinlichkeit ausgeht, daher droppen wir diese Spalte:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.drop(['Cabin'], axis=1, inplace=True)\n",
    "train.info();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fehlende Werte: Age\n",
    "\n",
    "Die Spalte \"Age\" hat einige fehlende Einträge. Wir haben die Alternativen entweder die Spalte auch zu droppen, \n",
    "dann verlieren wir allerdings einige Information, oder die fehlenden Werte aufzufüllen.\n",
    "\n",
    "Ein Ansatz dazu, ist zu schauen, wo die höchste Korrelation zu \"Age\" besteht und dies zum Auffüllen zu verwenden:\n",
    "\"Pclass\" ist am höchsten korreliert mit \"Age\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "age_group = train.groupby(\"Pclass\")[\"Age\"]\n",
    "print(age_group.median())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Der Altersmedian unterscheided sich zwischen den Klassen signifikant, also setzen wir diesen für die fehlenden Werte ein:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.loc[train.Age.isnull(), 'Age'] = train.groupby(\"Pclass\").Age.transform('median')\n",
    "print(train[\"Age\"].isnull().sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsere Altersverteilung sieht jetzt so aus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (16, 8))\n",
    "\n",
    "sns.distplot(train[\"Age\"])\n",
    "plt.title(\"Age Histogram\")\n",
    "plt.xlabel(\"Age\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fehlende Werte: Embarked\n",
    "\n",
    "Für die zwei fehlenden Ausgangshäfen nehmen wir einfach den häufigsten Wert an - \"S\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train[\"Embarked\"] = train[\"Embarked\"].fillna('S')\n",
    "train.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... keine fehlenden Werte mehr."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spalten droppen, die wir nicht weiter betrachten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.drop(['PassengerId', 'Name', 'Ticket', 'Fare'], axis = 1, inplace = True)\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One-hot Encoding \"Sex\" und \"Embarked\"\n",
    "\n",
    "Pandas hat eine sehr nützliche Funktion, die kategorische Variablen One-Hot encoded und den Dataframe gleich entsprechend umwandelt - diese heißt `get_dummies`.\n",
    "\n",
    "Der Name rührt daher, dass mit dem On-Hot Encoding neue Spalten (= neue Variablen) entstehen, die \"Dummy-Variablen\" genannt\n",
    "werden - da sie in gewissem Sinne ja keine \"echten\" Variablen sind.\n",
    "\n",
    "Bei Scikit-Learn leistet dasselbe der `OneHotEncoder`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_onehot = pd.get_dummies(train, drop_first=True)\n",
    "train_onehot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training vorbereiten\n",
    "\n",
    "`train` enthält ja sowohl unsere Features als auch das Label `Survived`. Dies teilen wir jetzt auf in X und y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = train_onehot['Survived']\n",
    "X = train_onehot.drop('Survived', axis=1)\n",
    "X, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noch mal eine Analyse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import feature_visualisation as fv\n",
    "fv.fc(X, 'clusters')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fv.fc(X, 'dendrogram')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fv.fc(X, 'punchcard')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fv.fc(X, 'heatmap')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schließlich splitten wir unsere Daten noch in ein Training- und ein Test-Set im Verhältnis 80:20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=101)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.shape, X_test.shape, y_train.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistische Regression\n",
    "\n",
    "Trainiere das Modell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "model = LogisticRegression()\n",
    "model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... und bewerte das Ergebnis ...\n",
    "\n",
    "**Confusion Matrix:** Die Confusion-Matrix $C$ ist definiert durch: $C_{i,j}$ ist die Anzahl der Beobachtungen der wahren Gruppe $i$, die als zur Gruppe $j$ vorhergesagt werden."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import classification_report, confusion_matrix\n",
    "\n",
    "predictions = model.predict(X_test)\n",
    "tn, fp, fn, tp = confusion_matrix(y_test, predictions).ravel()\n",
    "print(\"\\t'false'\\t'true'\")\n",
    "print(\"false\\t  {}\\t  {}\".format(tn, fp))\n",
    "print(\"true\\t  {}\\t  {}\".format(fn, tp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kennzahlen zur Bewertung\n",
    "\n",
    "**Precision** ist das Verhältnis $\\frac{tp}{tp + fp}$, wobei $tp$ die Anzahl der echten positiven und $fp$ die Anzahl der falschen positiven Werte ist. Intuitiv ist Precision die Fähigkeit des Klassifikators, eine negative Probe nicht als positiv zu kennzeichnen.\n",
    "\n",
    "**Recall** ist das Verhältnis $\\frac{tp}{tp + fn}$, wobei $fn$ die Anzahl der falschen Negativen ist. Intuitiv ist Recall die Fähigkeit des Klassifikators, alle positiven Proben zu finden.\n",
    "\n",
    "Der **F-Beta-Score** kann als ein gewichteter harmonischer Mittelwert von Precision und Recall interpretiert werden, wobei ein F-Beta-Score seinen besten Wert bei 1 und den schlechtesten Wert bei 0 erreicht.\n",
    "\n",
    "$f1=2\\times\\frac{Precision \\times Recall}{Precision + Recall}$\n",
    "\n",
    "Der F-Beta-Score gewichtet um den Faktor Beta Recall mehr als die Precision. Beta == 1,0 bedeutet, dass Recall und Precision gleich wichtig sind.\n",
    "\n",
    "**Support** ist die jeweilige Anzahl der Vorkommnisse der wahren Labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(classification_report(y_test, predictions, target_names=['Nicht überlebt', 'Überlebt']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation des Modells\n",
    "\n",
    "Schauen wir uns Intercept und Koeffizienten an"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Intercept/Bias: {}'.format(model.intercept_))\n",
    "\n",
    "coef_dict = sorted(list(zip(X.columns.tolist(), model.coef_.ravel())), key=lambda tup: tup[1])\n",
    "for tup in coef_dict:\n",
    "    print(tup)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Das Geschlecht ist mit Abstand der am stärksten eingehende Faktor, gefolgt von der Klasse.\n",
    "\n",
    "Da wir die Features nicht skaliert haben, stimmt die Reihenfolge bzgl SibSp, Parch und insbesondere Age so allerdings nicht."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X['Age'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "d.h. Age geht eigentlich mit einem fast 30-fachen Gewicht gegenüber den 0/1-Features ein und steht damit eher an dritter Stelle.\n",
    "\n",
    "Skalieren der Variable Age wäre sicherlich sinnvoll gewesen:\n",
    "\n",
    "${age}_{skaliert}=\\frac{age - mittleres Alter}{Ältestes Alter - Jüngstes Alter}=\\frac{age-29.066}{80}$\n",
    "\n",
    "SibSp und Parch könnte man vielleicht als Familie zusammenfassen (Familiengröße = 1 + SibSp + Parch) und dann trainieren.\n",
    "\n",
    "Das geht dann in Richtung \"Feature Engineering\". Ein interessanter Artikel dazu im Kontext des Titanic Datasets \n",
    "ist [hier](https://medium.com/i-like-big-data-and-i-cannot-lie/how-i-scored-in-the-top-9-of-kaggles-titanic-machine-learning-challenge-243b5f45c8e9) zu finden.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Erklärung des Modells mit Shap\n",
    "\n",
    "\n",
    "... Game theoretischer Ansatz zur Erklärung eines Modells\n",
    "\n",
    "![](Shap.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shap\n",
    "\n",
    "# load JS visualization code to notebook\n",
    "shap.initjs()\n",
    "\n",
    "# explain the model's predictions using SHAP\n",
    "explainer = shap.LinearExplainer(model, X, feature_perturbation=\"interventional\")\n",
    "shap_values = explainer.shap_values(X)\n",
    "#X_array = X.toarray() # we need to pass a dense version for the plotting functions\n",
    "shap.summary_plot(shap_values, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.force_plot(explainer.expected_value, shap_values, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.dependence_plot(\"rank(1)\", shap_values, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.decision_plot(explainer.expected_value, shap_values, X, link='logit')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Anhang: LogisticRegression with Cross Validation\n",
    "\n",
    "Abwandlung der \"normalen\" LogisticRegression mit automatischer Bestimmung der Hyperparameter.\n",
    "(Hier Regularisierung)\n",
    "\n",
    "Leicht bessere Ergebnisse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV\n",
    "\n",
    "model_cv = LogisticRegressionCV()\n",
    "model_cv.fit(X_train, y_train)\n",
    "\n",
    "predictions_cv = model_cv.predict(X_test)\n",
    "tn_cv, fp_cv, fn_cv, tp_cv = confusion_matrix(y_test, predictions_cv).ravel()\n",
    "print(\"\\t'false'\\t'true'\")\n",
    "print(\"false\\t  {}\\t  {}\".format(tn_cv, fp_cv))\n",
    "print(\"true\\t  {}\\t  {}\".format(fn_cv, tp_cv))\n",
    "\n",
    "print(classification_report(y_test, predictions_cv, target_names=['Nicht überlebt', 'Überlebt']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Übung\n",
    "\n",
    "Versuchen Sie das Modell zu verbessern, indem Sie das Alter skalieren und die Familiengröße verwenden."
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}