regression_estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d558ea7",
   "metadata": {},
   "source": [
    "# Starting least-squares regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43f22604",
   "metadata": {},
   "source": [
    "## Introduction and definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90fee4a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#: Import numerical and plotting libraries\n",
    "import numpy as np\n",
    "# Print to four digits of precision\n",
    "np.set_printoptions(precision=4, suppress=True)\n",
    "import numpy.linalg as npl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9877385",
   "metadata": {},
   "source": [
    "This exercise is to practice thinking about how the regression estimation\n",
    "works, and the relationship of correlation and regression.\n",
    "\n",
    "Here are another couple of samples of the “psychopathy” and “clamminess”\n",
    "scores, of the same type that you will see in the [on_regression\n",
    "page](https://textbook.nipraxis.org/on_regression.html):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "057791b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#: The data, that we are trying to model.\n",
    "psychopathy = np.array([10.3512,  4.7685,  8.8875, 15.8477,\n",
    "                        7.9221,  5.275 , 16.7839, 11.2392,\n",
    "                        13.2752, 11.0242, 21.9108, 21.0303])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34fb7535",
   "metadata": {},
   "outputs": [],
   "source": [
    "#: The regressor that we will use to model the data.\n",
    "clammy = np.array([0.0552,  1.8306, -0.3534, -0.7319,\n",
    "                   5.1021,  1.1133,  0.8126,  3.861,\n",
    "                   6.7069,  4.3189,  4.3741,  7.2527])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b295991",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of values in psychopathy and clammy\n",
    "n = len(psychopathy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db9856eb",
   "metadata": {},
   "source": [
    "The `psychopathy` scores are scores on a psychopathy scale from a standardized\n",
    "questionnaire, from 12 students.\n",
    "\n",
    "The `clammy` scores are scores of palm skin conductance from the same 12 students.\n",
    "\n",
    "Higher `clammy` scores correspond to a sweatier handshake.\n",
    "\n",
    "We predict that students with higher `clammy` scores will also score higher on\n",
    "the psychopathy scale.  In particular, we think there will be a linear\n",
    "(straight-line) relationship between these two.\n",
    "\n",
    "We want to find the slope and intercept for the \"best\" straight line to predict\n",
    "`psychopathy` scores from `clammy` scores.\n",
    "\n",
    "First do a scatter plot with `clammy` on the x-axis and `psychopathy` on the y-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b37683c8",
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "#- Plot psychopathy as a function of clammy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d8fe943",
   "metadata": {},
   "source": [
    "Before looking for the best line relating `clammy` to `psychopathy`, we will\n",
    "look at the correlation coefficient between these values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61c8a6aa",
   "metadata": {},
   "source": [
    "## Correlation coefficient\n",
    "\n",
    "As you remember we calculate the correlation by first calculating standard\n",
    "scores (z-scores), then taking the mean of the multiplied standard scores.\n",
    "\n",
    "Make a function to calculate and return the standard scores from an input\n",
    "array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60977892",
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_scores(v):\n",
    "    \"\"\" Return standard scores from input array `v`\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    v : array, shape (n,)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    z : array, shape (n,)\n",
    "        Standard scores corresponding to input array `v`.\n",
    "    \"\"\"\n",
    "   ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbe2d920",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()\n",
    "# Make 20 random numbers.\n",
    "my_vector = rng.normal(10, 2, size=20)\n",
    "my_vector_standard = standard_scores(my_vector)\n",
    "# Mean should be very close to 0\n",
    "assert np.isclose(np.mean(my_vector_standard), 0)\n",
    "# Standard deviation should be very close to 1\n",
    "assert np.isclose(np.std(my_vector_standard), 1)\n",
    "assert np.isclose(\n",
    "    my_vector_standard[0],\n",
    "    (my_vector[0] - np.mean(my_vector)) / np.std(my_vector))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71500c4f",
   "metadata": {},
   "source": [
    "Make a function that uses your `standard_scores` function to calculate the\n",
    "correlation coefficient between two 1D arrays (vectors):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80e87f97",
   "metadata": {},
   "outputs": [],
   "source": [
    "def corr(x, y):\n",
    "    \"\"\" Pearson correlation coefficient between vectors `x`, `y`.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array, shape (n,)\n",
    "    y : array, shape (n,)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    r : float\n",
    "        Pearson product-moment correlation coefficient between `x` and `y`.\n",
    "    \"\"\"\n",
    "   ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0019d950",
   "metadata": {},
   "source": [
    "Calculate the correlation coefficient between `clammy` and `psychopathy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bca8296",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = ...\n",
    "# Show the result\n",
    "r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19b85848",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert -1 <= r <= 1\n",
    "assert np.isclose(r, np.corrcoef(clammy, psychopathy)[0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c18bbf7f",
   "metadata": {},
   "source": [
    "Make a standard scored version of `clammy`.  Call it `z_x`.\n",
    "\n",
    "Make a standard scored version of `psychopathy`.  Call it `z_y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7225f21f",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_x = ...\n",
    "z_y = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1c0ecdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(z_x[0], (clammy[0] - np.mean(clammy)) / np.std(clammy))\n",
    "assert np.isclose(z_y[0], (psychopathy[0] - np.mean(psychopathy)) /\n",
    "    np.std(psychopathy))\n",
    "for vec in z_x, z_y:\n",
    "    assert np.isclose(np.mean(vec), 0)\n",
    "    assert np.isclose(np.std(vec), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f38b596",
   "metadata": {},
   "source": [
    "Do a scatterplot with `z_x` on the x-axis and `z_y` on the y-axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c70217c",
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "#- Plot z_y as a function of z_x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb60750f",
   "metadata": {},
   "source": [
    "Here is a function to calculate the Root Mean Squared Error (RMSE) when using a given\n",
    "slope and intercept to estimate values in `y`, from values in `x`.  You will\n",
    "recognize it from [on regression](https://textbook.nipraxis.org/on_regression).\n",
    "The only variation is that the function takes the square root of the MSE to\n",
    "give the Root MSE (RMSE)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "315dbd5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_rmse(x, y, slope, inter):\n",
    "    predicted = inter + slope * x\n",
    "    errors = y - predicted\n",
    "    return np.sqrt(np.mean(errors ** 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d2e16b8",
   "metadata": {},
   "source": [
    "Calculate the RMSE for `z_x` vs `z_y`, using `r` above as the slope\n",
    "and 0 as the intercept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8290f93",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_rmse = ...\n",
    "# Show the result\n",
    "r_rmse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b35655b1",
   "metadata": {},
   "source": [
    "Here is a selection of slopes and intercepts to try, to predict `z_y` with `z_x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03fec50e",
   "metadata": {},
   "outputs": [],
   "source": [
    "slopes = np.arange(-1, 1, 0.01)\n",
    "intercepts = np.arange(-1, 1, 0.01)\n",
    "n_vals = len(slopes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b2b94fd",
   "metadata": {},
   "source": [
    "For every slope, intercept pair from `slopes`, `intercepts` above, calculate\n",
    "the RMSE value and store it in the array `rmses`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37dba15e",
   "metadata": {},
   "outputs": [],
   "source": [
    "rmses = np.zeros((n_vals, n_vals))\n",
    "for i in range(n_vals):  # For each slope value.\n",
    "    slope = slopes[i]\n",
    "    for j in range(n_vals):\n",
    "       ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1275342",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.all(rmses >= 0)\n",
    "assert np.isclose(\n",
    "    rmses[50, 61],\n",
    "    calc_rmse(z_x, z_y, slopes[50], intercepts[61]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fd2cd3e",
   "metadata": {},
   "source": [
    "Show the `rmses` array as a grayscale image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33ad042a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Show rmse array as an image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b3af13d",
   "metadata": {},
   "source": [
    "Show that all values in `rmses` are greater than the RMSE you calculated for\n",
    "slope=`r` and intercept of 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b34965b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba8d1ee",
   "metadata": {},
   "source": [
    "The fact that this is so suggests that slope=`r` and intercept of 0 are in fact\n",
    "the slope, intercept pair that minimize the RMSE, when predicting `z_y` from\n",
    "`z_x`.\n",
    "\n",
    "Run the cell below to generate a plot which shows the RMSE for a selection of\n",
    "the candidate slope/intercept pairs (in green) — as well as for the\n",
    "slope/intercept pair that gives the minimum RMSE (shown as a red X) - when\n",
    "predicting `z_y` from `z_x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b22efdd8",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Run this cell.  This cell has some advanced plotting\n",
    "# that we do not teach on the course.\n",
    "slopes_2 = np.arange(-1, 1, 0.1)\n",
    "intercepts_2 = np.arange(-1, 1, 0.1)\n",
    "less_vals = len(slopes_2)\n",
    "\n",
    "fig = plt.figure(figsize = (10, 7))\n",
    "ax = plt.axes(projection =\"3d\")\n",
    "ax.set_xlabel('Intercept')\n",
    "ax.set_ylabel('Slope')\n",
    "ax.set_zlabel('RMSE')\n",
    "ax.view_init(elev=15, azim=15)\n",
    "ax.scatter3D(0, r, calc_rmse(z_x, z_y, r, 0),\n",
    "             color = \"red\", marker = 'X', s = 200)\n",
    "for i in range(less_vals):  # For each slope value.\n",
    "    slope = slopes_2[i]\n",
    "    for j in range(less_vals):\n",
    "        intercept = intercepts_2[j]\n",
    "        ax.scatter3D(intercept, slope,\n",
    "                     calc_rmse(z_x, z_y, slope, intercept),\n",
    "                     color = \"green\", alpha = 0.2);\n",
    "ax.view_init(15, 45)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "714f6e38",
   "metadata": {},
   "source": [
    "## Regression and correlation\n",
    "\n",
    "We have found the best straight line for predicting the psychopathy values as\n",
    "standard scores, from the clammy values as standard scores.  That line has\n",
    "slope `r` and intercept 0.\n",
    "\n",
    "We now want to find the best slope and intercept to predict the original\n",
    "`psychopathy` values from the original `clammy` values.\n",
    "\n",
    "Here are some slopes and intercepts to try for the line predicting\n",
    "`psychopathy` from `clammy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac9f613c",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_slopes = 200\n",
    "orig_slopes = np.linspace(0, 1, n_slopes)\n",
    "n_inters = 300\n",
    "orig_intercepts = np.linspace(5, 15, n_inters) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87ac68da",
   "metadata": {},
   "source": [
    "For every slope, intercept pair from `orig_slopes`, `orig_intercepts` above, calculate\n",
    "the RMSE value and store it in the array `rmses_orig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f019db9",
   "metadata": {},
   "outputs": [],
   "source": [
    "rmses_orig = np.zeros((n_slopes, n_inters))\n",
    "for i in range(n_slopes):  # For each slope value.\n",
    "   ...\n",
    "   ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb7031e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.all(rmses_orig >= 0)\n",
    "assert np.isclose(\n",
    "    rmses_orig[50, 61],\n",
    "    calc_rmse(clammy, psychopathy, orig_slopes[50], orig_intercepts[61]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d93d8930",
   "metadata": {},
   "source": [
    "Show the `rmses_orig` array as a grayscale image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b99fda1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Show rmse array as an image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f166f95f",
   "metadata": {},
   "source": [
    "From the [on regression](https://textbook.nipraxis.org/on_regression) page, we can work out the best (least squares, least RMSE) slope from `r`, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3d80c0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_orig_slope = r * np.std(psychopathy) / np.std(clammy)\n",
    "best_orig_slope"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "783db9aa",
   "metadata": {},
   "source": [
    "We can reconstruct the best intercept with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef89444d",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_orig_inter = np.mean(psychopathy) - best_orig_slope * np.mean(clammy)\n",
    "best_orig_inter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffd39c13",
   "metadata": {},
   "source": [
    "This is also what we find using the Scipy linear regression function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2a85fc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as sps\n",
    "res = sps.linregress(clammy, psychopathy)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d37ac9e",
   "metadata": {},
   "source": [
    "Calculate the RMSE value for `clammy`, `psychopathy`, using the best slope and\n",
    "intercept from the calculation above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dce6ae1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_rmse = ...\n",
    "# Show the result\n",
    "best_rmse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03f1be06",
   "metadata": {},
   "source": [
    "Show that all values in `rmses_orig` are greater than the RMSE you calculated for\n",
    "slope=`best_slope_orig` and intercept of `best_slope_orig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5d9fbf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1692a9d",
   "metadata": {},
   "source": [
    "This implies that slope, intercept pair we calculated using the `r` value, or\n",
    "with `scipy.stats.linregress`, does give the line minimizing the squared error\n",
    "(and RMSE) when predicting `psychopathy` from `clammy`.\n",
    "\n",
    "Run the cell below to generate a plot which shows the RMSE for a selection of\n",
    "the candidate slope/intercept pairs (in green). The plot marks the\n",
    "slope/intercept pair that gives the minimum RMSE with a red X."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc09ef49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell.  It has some advanced plotting that you\n",
    "# will not otherwise need for the course.\n",
    "n_slopes = 20\n",
    "orig_slopes = np.linspace(0, 1, n_slopes)\n",
    "n_inters = 30\n",
    "orig_intercepts = np.linspace(5, 15, n_inters)\n",
    "\n",
    "fig = plt.figure(figsize = (10, 7))\n",
    "ax = plt.axes(projection =\"3d\")\n",
    "ax.set_xlabel('Intercept')\n",
    "ax.set_ylabel('Slope')\n",
    "ax.set_zlabel('RMSE')\n",
    "ax.view_init(elev=15, azim=75)\n",
    "ax.scatter3D(best_orig_inter, best_orig_slope,\n",
    "             calc_rmse(clammy, psychopathy, best_orig_slope, best_orig_inter),\n",
    "             color = \"red\", marker = 'X', s = 200)\n",
    "\n",
    "for i in range(n_slopes):\n",
    "    slope = orig_slopes[i]\n",
    "    for j in range(n_inters):\n",
    "        intercept = orig_intercepts[j]\n",
    "        ax.scatter3D(intercept, slope,\n",
    "                     calc_rmse(clammy, psychopathy, slope, intercept),\n",
    "                     color = \"green\", alpha = 0.2);\n",
    "ax.view_init(15, 25)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "notebook_metadata_filter": "all,-language_info",
   "split_at_heading": true,
   "text_representation": {
    "extension": ".Rmd",
    "format_name": "rmarkdown",
    "format_version": "1.2",
    "jupytext_version": "1.10.3"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}