2_regression_short.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### SHORT VERSION\n",
    "# Introduction to data science: regression\n",
    "\n",
    "[**This notebook is available on Google Colab.**](https://colab.research.google.com/drive/1d1tKJPLF8EXh0LU7OK-VHiFpiSygQ_M8)\n",
    "\n",
    "Classification is all about predicting discrete classes, but sometimes we want to predict quantities &mdash; continuous numerical properties with magnitude, like temperature or time.\n",
    "\n",
    "As before, we want to select the best (i.e. optimal) model -- and, as before, this means having an objective measure of 'best' and a way to prove that you have found it. \n",
    "\n",
    "First we'll import some data. Again, I'm using an extract from [the Rock Property Catalog](https://subsurfwiki.org/wiki/Rock_Property_Catalog)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_excel('https://github.com/scienxlab/datasets/raw/refs/heads/main/rpc/3-lithologies-imputed.xlsx')\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's select a column as `X` (input) and a `y` (target)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "X = df[['Vp']].values\n",
    "y = df['Vs'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "## Ridge regression\n",
    "\n",
    "**Linear regression** is a reliable regressor. The 'vanilla' algorithm does not include any regularization, which is sometimes what you want. But we will use regularization from the start, because you will often want to apply it &mdash; and it gives me the chance to talk about regularization and hyperparameter tuning in general."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.scatter(X, y)\n",
    "plt.scatter(X, y_pred, c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot a line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "X_model = np.linspace(1400, 6000).reshape(-1, 1)\n",
    "y_model = est.predict(X_model)\n",
    "\n",
    "plt.scatter(X, y)\n",
    "plt.plot(X_model, y_model, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we are fitting a model like:\n",
    "\n",
    "$$ \\Large \\hat{y} = w x + b $$\n",
    "\n",
    "We can ask for the gradient ($w$) and the intercept ($b$) of the line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est.coef_, est.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**❓ How can we measure the quality of this prediction?**\n",
    "\n",
    "<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**❓ What do we think? Are we satisfied?**\n",
    "\n",
    "<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Scoring\n",
    "\n",
    "Scores matter in all machine learning tasks. It is very common to see people reporting only **R<sup>2</sup>** for regression tasks, or only **accuracy** for classification tasks. But it is almost never enough to only look at (or report) the 'obvious' score."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **mean squared error** was used for the optimization, but its square root, **RMS error**, is often easier to interpret because it has the same units as the original quantity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error, root_mean_squared_error\n",
    "\n",
    "mean_squared_error(y, y_pred), root_mean_squared_error(y, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **mean absolute error** may be even more intuitive, especially if we feel like negative and positive errors cancel each other out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "mean_absolute_error(y, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make ourselves a small report."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def report(y, y_pred):\n",
    "    print(f\"R²: {r2_score(y, y_pred):.4f}\", end='    ')\n",
    "    print(f\"RMSE: {root_mean_squared_error(y, y_pred):.1f}\", end='    ')\n",
    "    print(f\"MAE: {mean_absolute_error(y, y_pred):.1f}\")\n",
    "    return\n",
    "\n",
    "report(y, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Assumptions of regression\n",
    "\n",
    "Along with these metrics, it's essential to inspect the residuals to ensure that they are:\n",
    "\n",
    "- Approximately normally distributed.\n",
    "- Not correlated with the inputs or output.\n",
    "- Homoskedastic (the variance is not correlated with inputs or output).\n",
    "\n",
    "These conditions are assumptions of Gauss–Markov theorem, which underlies linear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "sns.kdeplot(y - y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X, y - y_pred, 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "\n",
    "## Fair evaluation\n",
    "\n",
    "Now that we have an evaluation criterion, we need to address another issue: fairness. Our test was not fair.\n",
    "\n",
    "We should not train the model then check its accuracy only on that same training dataset. It's cheating, because in the future we'd like to predict on data that the model has never seen. So we should test the model on data it has never seen.\n",
    "\n",
    "Let's hold out some validation data, or 'blind' data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "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)\n",
    "\n",
    "X_train.shape, y_train.shape, X_test.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### EXERCISE\n",
    "\n",
    "Now let's train a model _on only the training data_ and validate it properly _on only the test data_.\n",
    "\n",
    "**❓ Do we think the score will be better or worse than before?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "est = Ridge()\n",
    "\n",
    "est.fit( ... )  # Replace the dots with the training data.\n",
    "\n",
    "y_pred = est.predict( ... )  # What will we predict on?\n",
    "\n",
    "report( ... )  # What do we give the report function?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Was: `R²: 0.8801    RMSE: 215.6    MAE: 172.9`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**❓ What do we need to think about when splitting? In other words: what is the most important thing about the test data?**\n",
    "\n",
    "&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />\n",
    "\n",
    "First and foremost, the test data must look like data we expect in the future.\n",
    "\n",
    "We also need to think about:\n",
    "\n",
    "- Independence: can you shuffle the data wthout losing information?\n",
    "- Identical distributions: are all the data from the same distribution?\n",
    "- Reproducibility: what can we do to make this reproducible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.scatter(X_train, y_train)\n",
    "plt.scatter(X_test, y_test, c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How are we doing?\n",
    "\n",
    "Let's check the model against the test data only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_model = est.predict(X_model)\n",
    "\n",
    "plt.scatter(X_train, y_train, c='r', alpha=0.1)\n",
    "plt.scatter(X_test, y_test)\n",
    "plt.plot(X_model, y_model, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**❓ What could we do to improve this prediction?**\n",
    "\n",
    "&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />\n",
    "\n",
    "Quite a few things!\n",
    "\n",
    "- We have another log: we could use the density as well.\n",
    "- We have lithology, we could use that as well.\n",
    "- We could think about whether we need to preprocess the data in any way.\n",
    "- We could add nonlinear transformations and combinations of the features.\n",
    "- We could tune the hyperparameter of the model, `alpha`.\n",
    "- We could try other models.\n",
    "\n",
    "🕑 If we have time, we can try some of these things."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## More features\n",
    "\n",
    "We also have the rock density, `Rho`, let's use that. Now the model will be like:\n",
    "\n",
    "$$ \\Large \\hat{y} = w_0 x_0 + w_1 x_1 + b $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df[['Vp', 'Rho']].values\n",
    "y = df['Vs'].values\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)\n",
    "\n",
    "est = Ridge()\n",
    "est.fit(X_train, y_train)\n",
    "y_pred = est.predict(X_test)\n",
    "\n",
    "report(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Was originally:\n",
    "\n",
    "`R²: 0.8801    RMSE: 215.6    MAE: 172.9`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One downside is that it is now a bit harder to draw than it was. but we won't fix it as it's just about to get worse..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "## Non-linearity\n",
    "\n",
    "Linear regression is linear **in the parameters**, but we can model non-linear relationships in the data by adding non-linear transformations of the data. In particular, we will add the squares of the features, and the **interactions** of the features (their products, basically).\n",
    "\n",
    "We can do this with **polynomial expansion** as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "np.set_printoptions(precision=3)\n",
    "\n",
    "poly = PolynomialFeatures(2, include_bias=False)\n",
    "X_ = poly.fit_transform(X)\n",
    "X_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have:\n",
    "\n",
    "$$ \\Large \\hat{y} = w_0 x_0 + w_1 x_1 + w_2 x_0^2 + w_3 x_0 x_1 + w_4 x_1^2 + b $$\n",
    "\n",
    "So the model will have six parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X_, y, random_state=42)\n",
    "\n",
    "est = Ridge()\n",
    "est.fit(X_train, y_train)\n",
    "y_pred = est.predict(X_test)\n",
    "\n",
    "report(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Was previously:\n",
    "\n",
    "`R²: 0.8616    RMSE: 210.5    MAE: 165.7`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the coefficients now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est.coef_, est.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**❓ Big coefficient means big effect, right? Can we interpret these as importances?**\n",
    "\n",
    "<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "\n",
    "## Scale matters\n",
    "\n",
    "Let's look at the first few rows of the data (before expansion):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "X[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The units of these columns are m/s and g/cm<sup>3</sup>, which have very different magnitudes for rocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (axl, axr) = plt.subplots(figsize=(12, 4), ncols=2, sharey=True)\n",
    "axl.plot(X_train[:, 0], y_train, 'o')\n",
    "axl.plot(X_test[:, 0], y_test, 'ro')\n",
    "axr.plot(X_train[:, 1], y_train, 'o')\n",
    "axr.plot(X_test[:, 1], y_test, 'ro')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're really only interested in the _distributions_ of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "sns.kdeplot(X[:, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Standardize the data\n",
    "\n",
    "It's essential to train some machine learning algorithms on scaled data, for example on the Z-scores of your data, i.e. zero mean, unit variance. This ensures that the different scales of the features is not causing a problem.\n",
    "\n",
    "For features distributed uniformly, or with strong min/max constraints, another strategy like normalization (e.g. in (0, 1) or (-1, 1)) might be better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**❓ Standardization needs the mean and variance of the data... which dataset shall we measure these stats on?**\n",
    "\n",
    "<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;<br />&nbsp;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "scaler = StandardScaler()\n",
    "\n",
    "scaler.fit(X_train)  # Important!\n",
    "\n",
    "X_train_sc = scaler.transform(X_train)\n",
    "X_test_sc = scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This doesn't change how the data are distributed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, (axl, axr) = plt.subplots(figsize=(12, 4), ncols=2, sharey=True)\n",
    "axl.plot(X_train_sc[:, 0], y_train, 'o')\n",
    "axl.plot(X_test_sc[:, 0], y_test, 'ro')\n",
    "axr.plot(X_train_sc[:, 1], y_train, 'o')\n",
    "axr.plot(X_test_sc[:, 1], y_test, 'ro')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "sns.kdeplot(X_train_sc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Now we can re-fit the model and look at the scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "est = Ridge()\n",
    "est.fit(X_train_sc, y_train)\n",
    "y_pred = est.predict(X_test_sc)\n",
    "\n",
    "report(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Was:\n",
    "\n",
    "`R²: 0.8655    RMSE: 207.5    MAE: 163.1`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quality of the prediction does not change much, but the coefficients are more interesting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est.coef_, est.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that some people think you should not interpret the coefficients of the regularized model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Should you scale the target?\n",
    "\n",
    "I have not been able to get a clear answer on this, but in my experience it makes no difference except with neural networks (the reason being that large error values may not backpropagate properly during training)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Be careful!\n",
    "\n",
    "Solving one problem gives us a new one. It is now essential to scale the data now before inference -- although the model will happily make (terrible) predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "est.predict([X_train[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This is known as an \"out of distribution\" or OOD error, and it's a classic pitfall in machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "\n",
    "## Put everything in a pipeline\n",
    "\n",
    "This is the most flexible way to compose data pipelines in `sklearn`. It is better than implementing everything individually in a stepwise manner.\n",
    "\n",
    "For now, it won't change anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "pipe = make_pipeline(StandardScaler(), Ridge())\n",
    "\n",
    "pipe.fit(X_train, y_train)\n",
    "\n",
    "y_pred = pipe.predict(X_test)\n",
    "\n",
    "report(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hyperparameter tuning\n",
    "\n",
    "It is sensible to use cross-validation when tuning hyperparameters. Using the test data for tuning will overfit to the train-test split."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_validate, StratifiedKFold\n",
    "from sklearn.metrics import make_scorer\n",
    "\n",
    "df_train, df_test = train_test_split(df, stratify=df['Lithology'], random_state=42)\n",
    "\n",
    "pipe = make_pipeline(StandardScaler(), Ridge())\n",
    "skf = StratifiedKFold().split(X_train, df_train['Lithology'])\n",
    "\n",
    "cross_validate(pipe,\n",
    "               X_train, y_train,\n",
    "               scoring='neg_root_mean_squared_error',\n",
    "               cv=skf,\n",
    "               return_train_score=True,\n",
    "              )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use a loop to try lots of different values of `alpha`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = np.logspace(-4, 2, 13)\n",
    "alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_val, mean_train = [], []\n",
    "for alpha in alphas:\n",
    "    pipe = make_pipeline(StandardScaler(), Ridge(alpha=alpha))\n",
    "    skf = StratifiedKFold().split(X_train, df_train['Lithology'])\n",
    "    cross_val = cross_validate(pipe, X_train, y_train, scoring='neg_root_mean_squared_error', cv=skf, return_train_score=True)\n",
    "    mean_val.append(-np.mean(cross_val['test_score']))\n",
    "    mean_train.append(-np.mean(cross_val['train_score']))\n",
    "\n",
    "alpha_opt = alphas[np.argmin(mean_val)]\n",
    "\n",
    "plt.plot(alphas, mean_val, label='Validation')\n",
    "plt.plot(alphas, mean_train, c='r', label='Training')\n",
    "plt.fill_between(alphas, mean_val, mean_train, alpha=0.2)\n",
    "plt.xlabel('alpha (complexity)')\n",
    "plt.ylabel('RMSE (bias)')\n",
    "plt.xscale('log')\n",
    "plt.title('Bias–variance tradeoff curve')\n",
    "plt.axvline(alpha_opt, c='k', lw=0.67)\n",
    "plt.text(0.9*alpha_opt, 200, f\"← alpha={alpha_opt:.3f}\")\n",
    "plt.text(1.2*alpha_opt, 203, f\"UNDERFITTING\", horizontalalignment='right', c='gray')\n",
    "plt.text(0.8*alpha_opt, 203, f\"OVERFITTING\", horizontalalignment='left', c='gray')\n",
    "plt.text(1e-2, 195.5, f\"variance\", verticalalignment='center', c='C0')\n",
    "plt.legend()\n",
    "plt.gca().invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This curve is well known in the field of machine learning. It's called the **bias–variance trade-off curve** and I think it makes the most sense when the axes are labeled as shown. _Bias_ is basically 'error', and the shaded area is called _variance_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Categorical features\n",
    "\n",
    "Turns out the different lithologies have quite different rock physics characteristics. It makes a lot of sense to account for this in the model. You might think you want to train a different model for each lithology, but this is not necessary. It is sufficient to introduce the lithology variable into `X`, using 'dummy encoding'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "LITHS = ['limestone', 'dolomite', 'shale']\n",
    "\n",
    "def lith_index(y):\n",
    "    return [LITHS.index(lith) for lith in y]\n",
    "\n",
    "plt.scatter(*X.T, c=lith_index(df['Lithology']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more on how to do this, check out [this notebook](https://github.com/agilescientific/geocomputing/blob/develop/prod/Linear_regression.ipynb), and look for the **Take lithology into account** section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "\n",
    "## Carry on exploring!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "- How can you add Lithology to the features? Would it improve the prediction quality?\n",
    "- Choose another algorithm to try a prediction with, and implement it in a pipeline. For example, try KNN.\n",
    "- Choose a hyperparameter of the new algorithm and tune it. (If you have done this kind of thing before, try tuning 2 or 3 hyperparameters with grid or random search.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Test\n",
    "\n",
    "When you have tuned the predictor and are satisfied that it is as good as it can be, you can test against the holdout set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "alpha_opt = alphas[np.argmin(mean_val)]\n",
    "\n",
    "pipe = make_pipeline(StandardScaler(), Ridge(alpha=alpha_opt))\n",
    "\n",
    "pipe.fit(X_train, y_train)\n",
    "                    \n",
    "y_pred = pipe.predict(X_test)\n",
    "\n",
    "report(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Was:\n",
    "\n",
    "`R²: 0.8822    RMSE: 214.1    MAE: 161.9`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember we said the residuals should meet certain conditions? Let's check them again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(y_test - y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X_test[:, 0], y_test - y_pred, 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X_test[:, 1], y_test - y_pred, 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "If you are satisfied (think hard about what this means... you really have to decide before you start the model fitting process) then you are ready to fit the final model. If not, you must start all over again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Using this model\n",
    "\n",
    "We do not want to use this model &mdash; if we like its performance then we should now retrain it on all the data. Presumably, this new model will be at least as good as the one trained on the training set, we just don't have a way to check it now 😬"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "scaler = StandardScaler().fit(X_)\n",
    "X_sc = scaler.transform(X_)\n",
    "est = Ridge(alpha=alpha_opt).fit(X_sc, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "There is no way for us to test this model, but we should monitor it in production."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Takeaways\n",
    "\n",
    "There are lots, but here are a few.\n",
    "\n",
    "- **Machine learning is programming.** For now, for science and engineering applications, there is no way around this. Automatic machine learning is not a thing, at least not yet, not for us.\n",
    "- **Learn about the algorithms you apply.** For linear regression, you need to know about its cost function, its assumptions, and regularization.\n",
    "- **Give a lot of thought to how to fairly test your model.** Start by thinking about how it will be applied\n",
    "- **Scale your features.** It never hurts and sometimes it's essential.\n",
    "- **Use polynomial expansion on regression tasks.** It will lmost always improve your model.\n",
    "- **Use pipelines in `sklearn`.** You will avoid a lot of headaches and gotchas with preprocessing your data.\n",
    "- **Use regularization.** It usually makes sense for predictive applications, but make sure you understand how it works and note that L2 ('ridge') regularization may not be the best strategy for your application.\n",
    "- **Use appropriate measures of performance.** For example, make sure your metrics are compatible with the loss function.\n",
    "- **Check the assumptions of linear regression hold for your solution.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "&copy; 2025 Matt Hall / Equinor CC BY"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py311",
   "language": "python",
   "name": "py311"
  },
  "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}