2a - unit 07.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "edd76bd1",
   "metadata": {},
   "source": [
    "# Exercise: Improving a logistic regression model\n",
    "\n",
    "In the previous exercise, we fit a simple logistic regression model to predict the chance of an avalanche. This time, we will improve its performance by using multiple features intelligently.\n",
    "\n",
    "## Data visualisation\n",
    "\n",
    "Let's load our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eec3f24",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas\n",
    "import graphing # custom graphing code. See our GitHub repo for details\n",
    "\n",
    "#Import the data from the .csv file\n",
    "dataset = pandas.read_csv('Data/avalanche.csv', delimiter=\"\\t\", index_col=0)\n",
    "\n",
    "# Split our data into training and test\n",
    "import sklearn.model_selection\n",
    "train, test = sklearn.model_selection.train_test_split(dataset, test_size=0.25, random_state=10)\n",
    "\n",
    "print(\"Train size:\", train.shape[0])\n",
    "print(\"Test size:\", test.shape[0])\n",
    "\n",
    "#Let's have a look at the data\n",
    "print(train.head())"
   ]
  },
  {
   "source": [
    "We have numerous features available:\n",
    "\n",
    "* `surface_hoar` is how disturbed the surface of the snow is\n",
    "* `fresh_thickness` is how thick the top layer of snow is, or 0 if there is no fresh snow on top\n",
    "* `wind` is the top wind speed that day, in km/h\n",
    "* `weak_layers` is the number of layers of snow that are not well-bound to other layers\n",
    "* `no_visitors` is the number of hikers who were on the trail that day\n",
    "* `tracked_out` is a `1` or `0`. A `1` means that the snow has been trampled heavily by hikers \n",
    "\n",
    "## Simple logistic regression\n",
    "\n",
    "Let's make a simple logistic regression model and assess its performance with accuracy.\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "from sklearn.metrics import accuracy_score\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# Perform logistic regression.\n",
    "model = smf.logit(\"avalanche ~ weak_layers\", train).fit()\n",
    "\n",
    "# Calculate accuracy\n",
    "def calculate_accuracy(model):\n",
    "    '''\n",
    "    Calculates accuracy\n",
    "    '''\n",
    "    # Make estimations and convert to categories\n",
    "    avalanche_predicted = model.predict(test) > 0.5\n",
    "\n",
    "    # Calculate what proportion were predicted correctly\n",
    "    # We can use sklearn to calculate accuracy for us\n",
    "    print(\"Accuracy:\", accuracy_score(test.avalanche, avalanche_predicted))\n",
    "\n",
    "calculate_accuracy(model)"
   ]
  },
  {
   "source": [
    "Let's see how we can improve our model\n",
    "\n",
    "## Utilizing multiple features\n",
    "\n",
    "Most of our features seem like they could be useful, least in theory. Let's try a model with all features we have available."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "# Perform logistic regression.\n",
    "model_all_features = smf.logit(\"avalanche ~ weak_layers + surface_hoar + fresh_thickness + wind + no_visitors + tracked_out\", train).fit()\n",
    "calculate_accuracy(model_all_features)"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": []
  },
  {
   "source": [
    "That's a big improvement on the simpler model we've been working with. \n",
    "\n",
    "To understand why, we can look at the summary information"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_all_features.summary()"
   ]
  },
  {
   "source": [
    "Take a look at the P column, recalling that values less than 0.05 mean we can be confident that this parameter is helping the model make better predictions.\n",
    "\n",
    "Both `surface_hoar` and `wind` have very small values here, meaning they are useful predictors and probably explain why our model is working better. If we look at the `coef` (which states _parameters_) column we see that these have positive values. This means that higher winds, and greater amounts of surface hoar result in higher avalanche risk.\n",
    "\n",
    "## Simplifying our model\n",
    "\n",
    "Looking at the summary again, we can see that `tracked_out` (how trampled the snow is), and `fresh_thickness` have large p-values. This means they are not useful predictors. Let's see what happens if we remove them from our model:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform logistic regression.\n",
    "model_simplified = smf.logit(\"avalanche ~ weak_layers + surface_hoar + wind + no_visitors\", train).fit()\n",
    "calculate_accuracy(model_simplified)"
   ]
  },
  {
   "source": [
    "Our new model works very similarly to the old one! In some circumstances simplifying a model like this can even improve it, as it becomes less likely to overfit.\n",
    "\n",
    "## Careful feature selection\n",
    "\n",
    "Usually, we don't just pick features blindly. Let's think about what we've just done - we removed how much fresh snow was in a model, trying to predict avalanches. Something seems off. Surely avalanches are much _more_ likely after it has snowed? Similarly, the number of people on the track seems unrelated to how many avalanches there were, but we know that people often can trigger avalanches.\n",
    "\n",
    "Let's review our earlier model again:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "model_all_features.summary()"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": []
  },
  {
   "source": [
    "Look at the `fresh_thickness` row. We're told that it has a negative coefficient. This means that as thickness increases, avalanches decrease.\n",
    "\n",
    "Similarly, `no_visitors` has a negative coefficient, meaning that fewer hikers means more avalanches.\n",
    "\n",
    "How can this bes? Well, while visitors can cause avalanches if there is a lot of fresh snow, presumably they cannot do so easily if there is no fresh snow. This means that our features are not fully independent.\n",
    "\n",
    "We can tell the model to try to take into account that these features interact, using a multiply sign. Let's try that now."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a model with an interaction. Notice the end of the string where\n",
    "# we have a multiply sign between no_visitors and fresh_thickness\n",
    "formula = \"avalanche ~ weak_layers + surface_hoar + wind + no_visitors * fresh_thickness\"\n",
    "model_with_interaction = smf.logit(formula, train).fit()\n",
    "calculate_accuracy(model_with_interaction)"
   ]
  },
  {
   "source": [
    "The model has improved to 84% accuracy! Let's look at the summary information:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_with_interaction.summary()"
   ]
  },
  {
   "source": [
    "We can see that the interaction term is helpful - the p-value is less than 0.05. The model is also performing better than our previous attempts.\n",
    "\n",
    "## Making predictions with multiple features\n",
    "\n",
    "Very quickly, lets explore what this interaction means by looking at model predictions. \n",
    "\n",
    "We will first graph two independent features in 3D. Let's start with `weak_layers` and `wind`:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphing.model_to_surface_plot(model_with_interaction, [\"weak_layers\", \"wind\"], test)"
   ]
  },
  {
   "source": [
    "The graph is interactive - rotate it and explore how there is a clear s-shaped relationship between the features and probability.\n",
    "\n",
    "Let's now look at the features that we have said can interact:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphing.model_to_surface_plot(model_with_interaction, [\"no_visitors\", \"fresh_thickness\"], test)"
   ]
  },
  {
   "source": [
    "It looks quite different to the other! From any side, we can see an s-shape, but these combine in strange ways.\n",
    "\n",
    "We can see that the risk goes up on days with lots of visitors _and_ lots of snow. There is no real risk of avalanche when there is a lot of snow but no visitors, or when there are a lot of visitors but no snow.\n",
    "\n",
    "The fact that it shows high risk when there is no fresh snow and no visitors could be due to rain, which keeps visitors and snow clouds away but results in avalanches of the older snow. To confirm this, we'd need to explore the data in more depth, but we'll stop here for now.\n",
    "\n",
    "## Summary\n",
    "\n",
    "Well done! Let's recap. We've:\n",
    "\n",
    "* improved our simple model by adding more features. \n",
    "* practiced interpreting our model coefficients (parameters) from the model summary\n",
    "* eliminated unnecessary features\n",
    "* explored how sometimes it's important to think about what your data really mean\n",
    "* created a model that combined features to give superior result"
   ],
   "cell_type": "markdown",
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.7.10 64-bit ('ml_crash_course': conda)"
  },
  "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.10"
  },
  "interpreter": {
   "hash": "a970c7e51787c6ef7f3c805442fbb5a30ae2745635e1d2514e62d936ba85c55c"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}