week_4/4.3_logistic_regression/logistic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 4.3: Logistic Regression\n",
    "\n",
    "This lecture, we are going to use logistic regression to predict the quality of wines.\n",
    "\n",
    "**Learning goals:**\n",
    "\n",
    "- solve a classification problem\n",
    "- fit a logistic regression model\n",
    "- predict using a logistic regression model\n",
    "- visualize the decision boundary of a logistic regression model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "Humans have been making wine for more than [8000 years](https://www.nationalgeographic.com/travel/destinations/asia/georgia/sponsor-content-secret-birthplace-of-wine/),  and that's a lot of bottles. 🍷 There's too many options! Let's use data science to help choose the _tastiest_ red wines. Or at least, to classify high rated wines vs low rated wines based on some chemical features. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Data Munging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the dataset into a `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('wine_quality.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each row of our dataset is a wine.\n",
    "\n",
    "The `tasty` column indicates if a wine is \"tasty\" or not. i.e if it has a high rating. Even if it is represented by an integer, it is a _categorical_ column, as values are either 0 or 1. This will be the _label_ we are aiming to _predict_ in this _classification_ task.\n",
    "\n",
    "The other columns are numerical features. Their means are close to 0, which suggests that this dataset was normalized, like for lecture 3.5.\n",
    "\n",
    "To keep things simple and easy to visualize, let's start by using only two features to predict the \"tastiness\" of these wines: `sulphates` and `alcohol`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df[['sulphates', 'alcohol']]\n",
    "y = df['tasty']\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our features are already a _matrix_, we do not need to convert our `X` and `y` to numpy arrays and reshape them. Which means we're ready for training!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Training\n",
    "\n",
    "We'll use sklearn's [`LogisticRegression`](https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression) model which, as always, will make things simple with the `.fit()` and `.predict()` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "clf = LogisticRegression(random_state=0).fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🧠 Can you recall all the optimization steps that sklearn takes care of in the `.fit()` method? Refer to lectures 4.1 & 4.2 if this isn't clear.\n",
    "\n",
    "Just like for linear regression, we can investigate the _parameters_ of our model, $\\theta$.\n",
    "They won't be _interpreted_ in the same way however, since our hypothesis uses the logistic function. Nevertheless, sklearn still calls these parameters `intercept_` and `coef_`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "theta = np.append(clf.intercept_, clf.coef_)\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that the optimization procedure found that the logistic hypothesis which best predits the `tasty` labels is:\n",
    "\n",
    "$$\\begin{align}\n",
    "h_{\\theta}(x_{1}, x_{2}) &= \\frac{1}{1 + e^{-(\\theta_{0} + \\theta_{1} x_{1} + \\theta_{2} x_{2})}} \\\\\n",
    "h_{\\theta}(x_{1}, x_{2}) &= \\frac{1}{1 + e^{-(0.244 + 0.451 x_{1} + 1.099 x_{2})}}\n",
    "\\end{align}$$\n",
    "\n",
    "where $x_{1}$ is our first feature, the `sulphates` column, and $x_{2}$ is our second feature, the `alcohol` column.\n",
    "\n",
    "This is a function with two variables that we can plot in 3D:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from mpl_toolkits import mplot3d\n",
    "import math\n",
    "\n",
    "sns.set()\n",
    "\n",
    "def hypothesis(x1, x2):\n",
    "    linear_term = theta[0] + theta[1] * x1 + theta[2] * x2\n",
    "    return 1 / (1 + math.exp(-linear_term))\n",
    "\n",
    "v_hypothesis = np.vectorize(hypothesis)\n",
    "    \n",
    "def create_mesh(theta):\n",
    "    x_values = np.linspace(-8, 8, 40)\n",
    "    y_values = np.linspace(-8, 8, 40)\n",
    "    x_mesh, y_mesh = np.meshgrid(x_values, y_values)\n",
    "    z_mesh = v_hypothesis(x_mesh, y_mesh)\n",
    "    return x_mesh, y_mesh, z_mesh\n",
    "    \n",
    "x_mesh, y_mesh, z_mesh = create_mesh(theta)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "ax.plot_wireframe(x_mesh, y_mesh, z_mesh, rstride=5, cstride=5, alpha=0.4, color='g')\n",
    "\n",
    "ax.set_title('Logistic Regression Model of Wine Tastiness')\n",
    "ax.set_xlabel('Sulphates (norm)')\n",
    "ax.set_ylabel('Alcohol (norm)')\n",
    "ax.set_zlabel('Probability of Being Tasty');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🧠 Is this the shape you expected? Can you relate this to the 1D plots of the logistic function we saw in the slides?\n",
    "\n",
    "We could add our data points to this 3D graph. However, things would get crowded, and _encoding_ a binary label on the z-axis of a 3D plot isn't very clear. We end up with some points near the \"roof\" of the volume, and others near the \"floor\", without a good way of comparing them.\n",
    "\n",
    "Instead, we can use _color_ to encode the labels of the data points. We'll visualize the wines with the `tasty=1` label in red, and the non `tasty` wines in blue. In short, we are using a _2D graph_, where:\n",
    "- on the x-axis, the `sulphates` feature\n",
    "- on the y-axis, the `alcohol` feature\n",
    "- as color, the `tasty` binary label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "colors = sns.color_palette('husl').as_hex()\n",
    "\n",
    "fig = plt.figure(dpi=120)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# plot points\n",
    "positives = X[y == 1]\n",
    "negatives = X[y == 0]\n",
    "ax.scatter(positives.iloc[:, 0], positives.iloc[:, 1], alpha=0.4, s=8,c=colors[0], label='tasty')\n",
    "ax.scatter(negatives.iloc[:, 0], negatives.iloc[:, 1], alpha=0.4, s=8, c=colors[4], label='not tasty')\n",
    "\n",
    "# formatting\n",
    "ax.set_ylim((-2, 4))\n",
    "ax.set_xlim((-2, 4))\n",
    "ax.set_xlabel('Sulphates (norm)')\n",
    "ax.set_ylabel('Alcohol (norm)')\n",
    "ax.legend()\n",
    "ax.set_title('Wine Tastiness vs Sulphates & Alcohol Contents');\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "😨 This data is messy! This is a good example of what _real world data_ looks like. No clean clusters, no clean cut classes. How easy it is to split the dataset into its categorical labels is often called the _separability_ of the data. A dataset that is _separable_ has clear optimal decision boundaries that would lead to few prediction errors. In this case, it's not even clear to a human where to draw the line between `tasty` and not `tasty`!\n",
    "\n",
    "Let's look at our logistic regression model's attempt at separating the labels... maybe it has tasted more bottles than us 😏. Here, the prediction probabilities we visualized in the 3D graph isn't really what we wish to see. We would like to compare the _binary predictions_ against the binary labels. \n",
    "\n",
    "For this, we have to show the _decision boundary_ of our model. But all we have are the model parameters, $\\theta$ ! How do we go from the shape of the _hypothesis_ to the shape of the _decision boundary_ ? 🤨\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Decision Boundary\n",
    "\n",
    "Remember we defined how to go from prediction probability to binary prediction in the slides. We choose a threshold such that:  \n",
    "\n",
    "if $h_{\\theta}(\\textbf{x}) >= 0.5$ :  \n",
    "    _predict 1_  \n",
    "else:   \n",
    "    _predict 0_\n",
    "\n",
    "The decision boundary is the set of points where predictions \"switch\" from 0 to 1. This means we can represent the it as all the points where:\n",
    "\n",
    "$$\n",
    "h_{\\theta}(x_{1}, x_{2}) = 0.5 \\label{eq:threshold} \\tag{1}\n",
    "$$\n",
    "\n",
    "However, this doesn't help us much with plotting it in matplotlib... Ideally, we expect a line equation of the form:\n",
    "\n",
    "$$\n",
    "y = b + ax\n",
    "$$\n",
    "\n",
    "And since we're using a 2D plot where our feature `sulphates` ($x_{1}$) is on the x-axis, and our feature `alcohol` ($x_{2}$) is on the y-axis, we can rewrite the above as:\n",
    "\n",
    "$$\n",
    "x_{2} = b + ax_{1}\n",
    "$$\n",
    "\n",
    "So this is our goal. Now, let's use a bit of mathematical trickery to get there. 🧙‍♀️  \n",
    "We can substitute the logistic hypothesis in equation $\\eqref{eq:threshold}$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{1}{1 + e^{-(\\theta_{0} + \\theta_{1} x_{1} + \\theta_{2} x_{2})}} &= 0.5 \\label{eq:logistic} \\tag{2} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "You might recall from the [plot](https://upload.wikimedia.org/wikipedia/commons/8/88/Logistic-curve.svg) of the logistic function that it is equal to $0.5$ when its input is $0$.  \n",
    "So we can rewrite equation $\\eqref{eq:logistic}$ as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\theta_{0} + \\theta_{1} x_{1} + \\theta_{2} x_{2} &= 0 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Finally, we can re-arrange the terms to get the line equation:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\theta_{0} + \\theta_{1} x_{1} + \\theta_{2} x_{2} &= 0 \\\\\n",
    "x_{2} &= \\frac{-(\\theta_{0} + \\theta_{1} x_{1})}{\\theta_{2}} \\\\\n",
    "x_{2} &= \\frac{-\\theta_{0}}{\\theta_{2}} + \\frac{ \\theta_{1}}{\\theta_{2}}x_{1} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "This is the equation of a _line_ with intercept $\\frac{-\\theta_{0}}{\\theta_{2}}$ and slope $\\frac{ \\theta_{1}}{\\theta_{2}}$. We can plot this in our 2D graph! 🎊 \n",
    "\n",
    "We'll add our logistic regression model's decision boundary as a green line:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "colors = sns.color_palette('husl').as_hex()\n",
    "\n",
    "fig = plt.figure(dpi=120)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# plot points\n",
    "positives = X[y == 1]\n",
    "negatives = X[y == 0]\n",
    "ax.scatter(positives.iloc[:, 0], positives.iloc[:, 1], alpha=0.4, s=8,c=colors[0], label='tasty')\n",
    "ax.scatter(negatives.iloc[:, 0], negatives.iloc[:, 1], alpha=0.4, s=8, c=colors[4], label='not tasty')\n",
    "\n",
    "\n",
    "def decision_boundary(x1, theta):\n",
    "    return -(x1 * theta[1] + theta[0])/theta[2]\n",
    "\n",
    "# plot decision boundary\n",
    "x1 = np.array(ax.get_xlim())\n",
    "x2 = decision_boundary(x1, theta)\n",
    "ax.plot(x1, x2, '--', c=colors[2]);\n",
    "\n",
    "# formatting\n",
    "ax.set_ylim((-2, 4))\n",
    "ax.set_xlim((-2, 4))\n",
    "ax.set_xlabel('Sulphates (norm)')\n",
    "ax.set_ylabel('Alcohol (norm)')\n",
    "ax.legend()\n",
    "ax.set_title('Wine Tastiness vs Sulphates & Alcohol Contents');\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The green line is all the points where $h_{\\theta}(x_{1}, x_{2}) = 0.5$. The points _above_ the line are the model's _positive_ predictions, i.e it thinks those wines are `tasty`, whilst the points _below_ the line are the model's _negative_ predictions, i.e it thinks those wines are not so `tasty`. Take you time to understand this graph, as it is typical of classification tasks. \n",
    "\n",
    "ℹ️ Notice that the predictions don't always match the true _label_ of the data. For example, take the blue point at $(-1.75, 2.4)$. The blue color shows that the label is 0, i.e a _not_ `tasty` wine. However, it is located _above_ the decision boundary, so the model think it is `tasty`. It's a classification error! Based on the _separability_ of the dataset, we can understand that the model makes a few mistakes... the goal of supervised learning algorithms isn't to perfectly predict a variable, but to find the model parameters that lead to the best _possible_ predictions.\n",
    "\n",
    "\n",
    "🧠 Can you match this 2D graph with the 3D representation of the prediction probabilities of section 3? Where would the decision boundary be on the 3D graph?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Prediction\n",
    "\n",
    "We are happy with our model, and want to try it out as soon as possible. During a visit to the supermarket, we find a strange looking wine bottle. It piques our curiosity, so we read the label: `sulphates=1.942` and `alcohol=-0.343` (in normalized units). Let's figure out if our model thinks this is a tasty bottle!\n",
    "\n",
    "🧠 Can you graphically read the model prediction for this strange wine using the graph above?\n",
    "\n",
    "We know by now that sklearn's `.predict()` method expects a matrix. Since we have only one example with two features, $x_{1}$ and $x_{2}$, it will be automatically instantiated as a 2D vector. Instead, let's reshape it to a ($1 \\times 2$) `ndarray` matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sulphates = 1.942\n",
    "alcohol = -0.343\n",
    "x_strange_wine = np.asarray([sulphates, alcohol]).reshape(1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then predict the wine's tastiness using the logistic regression model that we trained earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_strange_wine = clf.predict(x_strange_wine)\n",
    "y_strange_wine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prediction is positive, $1$, meaning our model thinks that this wine is `tasty`. Hopefully this means it doesn't taste as strange as it looks!\n",
    "\n",
    "We can compare this bottle to the rest of our dataset, by plotting it on the 2D graph from section 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "colors = sns.color_palette('husl').as_hex()\n",
    "\n",
    "fig = plt.figure(dpi=120)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# plot points\n",
    "positives = X[y == 1]\n",
    "negatives = X[y == 0]\n",
    "ax.scatter(positives.iloc[:, 0], positives.iloc[:, 1], alpha=0.4, s=8,c=colors[0], label='tasty')\n",
    "ax.scatter(negatives.iloc[:, 0], negatives.iloc[:, 1], alpha=0.4, s=8, c=colors[4], label='not tasty')\n",
    "\n",
    "\n",
    "def decision_boundary(x1, theta):\n",
    "    return -(x1 * theta[1] + theta[0])/theta[2]\n",
    "\n",
    "# plot decision boundary\n",
    "x1 = np.array(ax.get_xlim())\n",
    "x2 = decision_boundary(x1, theta)\n",
    "ax.plot(x1, x2, '--', c=colors[2]);\n",
    "\n",
    "# plot strange wine\n",
    "ax.plot(x_strange_wine[:, 0], x_strange_wine[:, 1], marker='x', markersize=20, markeredgewidth=4, color='r')\n",
    "\n",
    "# formatting\n",
    "ax.set_ylim((-2, 4))\n",
    "ax.set_xlim((-2, 4))\n",
    "ax.set_xlabel('Sulphates (norm)')\n",
    "ax.set_ylabel('Alcohol (norm)')\n",
    "ax.legend()\n",
    "ax.set_title('Wine Tastiness vs Sulphates & Alcohol Contents');\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can indeed see that this strange looking wine is placed _above_ the decision boundary, which is why the `.predict()` method returned a value of $1$. 🙌\n",
    "\n",
    "You might recall that we had many more features available in our dataset `df`. Going beyond two features means that we won't be able to plot the data or the decision boundary in a meaningful manner. However, the predictions are likely to be more accurate!\n",
    "\n",
    "💪💪 Train a new logistic regression model named `full_clf` using _all_ the features available in our dataset `df`. Then predict the tastiness of the `mystery_wine` 🕵️‍♀️, and store it in a variable called `y_mystery_wine`.  \n",
    "Pro-tip: don't use the `tasty` column in your feature matrix X! That should be reserved as label vector, y.  \n",
    "Pro-tip 2: please use the `random_state=0` argument when constructing your sklearn class, to make sure your results are reproducible. Else, the unit test won't pass!  \n",
    "Pro-tip 3: If you also don't like typing out a big list of column names, you might want to consider [`.drop()`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.drop.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mystery_wine = np.array([-0.12, 0.82, 1.43, -0.29, 0.03, -0.90, -0.55, 2.19, -0.62, 0.43, 1.85]).reshape(1, -1)\n",
    "\n",
    "# INSERT CODE\n",
    "\n",
    "def test_mystery_wine_prediction():\n",
    "    assert y_mystery_wine == 1, 'The mystery wine is actually tasty! But your prediction says otherwise 😢'\n",
    "    prob_mystery = full_clf.predict_proba(mystery_wine)[:, 1]\n",
    "    assert math.isclose(prob_mystery, 0.7615, rel_tol=1e-3), f'Your prediction probability for the mystery wine is wrong'\n",
    "    print('Success! 🎉')\n",
    "    \n",
    "test_mystery_wine_prediction()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Summary\n",
    "\n",
    "Today, we learned about our first classification method, **logistic regression**. We restated that classification tasks aim to predict **categorical** labels, and visualized this for one feature in a two dimensional plot. Linear models were a bad fit binary labels, so we discovered the **logistic function** that could predict **probabilities** between 0 & 1. We described how **thresholds** can turn these probabilities into a binary prediction. The logistic function can be **parameterized** by feeding a **linear model** as its argument, like a set of russian dolls. This means that we can **optimize** our logistic model using the **cross-entropy** cost function and **gradient descent**. We explained what **decision boundaries** are, and how important they are for **classifiers**. We also extended our logistic regression to the **multi-class** setting, where it could predict **multiple labels** with a **softmax** hypothesis. Finally, we put this knowledge to practice by predicting if wines were **tasty**, and visualizing the **decision boundary** of our trained **logistic regression** model.\n",
    "\n",
    "# Resources\n",
    "### Core Resources\n",
    "\n",
    "- [Coursera Machine Learning - Logistic Regression](https://www.coursera.org/lecture/machine-learning/classification-wlPeP)  \n",
    "Andrew Ng's always excellent course goes into detail for this section on Logistic Regression.\n",
    "- [Kaggle dataset - red wine quality](https://www.kaggle.com/uciml/red-wine-quality-cortez-et-al-2009)  \n",
    "The dataset used in this notebook.\n",
    "\n",
    "\n",
    "### Additional Resources\n",
    "\n",
    "- [ML cheatsheet cost functions](https://ml-cheatsheet.readthedocs.io/en/latest/loss_functions.html)  \n",
    "A handy and mathematical summary of common cost functions, including MSE and cross-entropy.\n",
    "- [Cross-entropy loss a visual explanation](https://towardsdatascience.com/understanding-binary-cross-entropy-log-loss-a-visual-explanation-a3ac6025181a)  \n",
    "A clear visual yet formal explanation of the cross-entropy cost function.\n",
    "- [Binary vs multi-class logistic regression](https://chrisyeh96.github.io/2018/06/11/logistic-regression.html#multinomial-logistic-regression-via-cross-entropy)  \n",
    "Comparison of logistic vs softmax functions for logistic regression models.\n",
    "- [Understanding the softmax function in minutes](https://medium.com/data-science-bootcamp/understand-the-softmax-function-in-minutes-f3a59641e86d)  \n",
    "Comprehensive blogpost defining the softmax function and its uses in Machine Learning models.\n"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}