hacker_stats_in_python_part_2.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "hacker_stats_in_python_part_2.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "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.7"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Amandianeze/go/blob/master/hacker_stats_in_python_part_2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "6Ijg5wUCTQYG"
      },
      "source": [
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/datacamp/python-live-training-template/blob/master/assets/datacamp.svg?raw=True\" alt = \"DataCamp icon\" width=\"50%\">\n",
        "</p>\n",
        "<br><br>\n",
        "\n",
        "\n",
        "This is part two of a pair of live training sessions to review the concepts of statistical inference laid out in the DataCamp courses [Statistical Thinking in Python I](https://learn.datacamp.com/courses/statistical-thinking-in-python-part-1) and [Statistical Thinking in Python II](https://learn.datacamp.com/courses/statistical-thinking-in-python-part-2). The concepts of those courses were reviewed and fortified in [Case Studies in Statistical Thinking](https://learn.datacamp.com/courses/case-studies-in-statistical-thinking). You can link to those courses by clicking on the badges below.\n",
        "\n",
        "<div style=\"margin: auto; width: 400px;\">\n",
        "\n",
        "<a href=\"https://learn.datacamp.com/courses/statistical-thinking-in-python-part-1\"><img src = \"https://assets.datacamp.com/production/course_1549/shields/original/shield_image_course_1549_20191010-1-13inj9n?1570728356\" width=\"100px\"></a>\n",
        "<a href=\"https://learn.datacamp.com/courses/statistical-thinking-in-python-part-2\"><img src = \"https://assets.datacamp.com/production/course_1550/shields/original/shield_image_course_1550_20180911-13-4iztc?1536680442\" width=\"100px\"></a>\n",
        "<a href=\"https://learn.datacamp.com/courses/case-studies-in-statistical-thinking\"><img src = \"https://assets.datacamp.com/production/course_4674/shields/original/shield_image_course_4674_20190318-13-1c5k8sf?1552948884\" width=\"100px\"></a> \n",
        "\n",
        "</div>\n",
        "\n",
        "In the Hacker Stats in Python Part 1 live training, we learned some basic ideas in nonparametric statistical inference using the plug-in principle. We applied those concepts to computing plug-in estimates and confidence intervals. Here, we will reinforce and extend those skills to:\n",
        "\n",
        "- Construct confidence intervals and perform NHSTs comparing two treatments in a data set.\n",
        "- Compute confidence intervals on two-dimensional data sets, investigating correlations and linear regressions.\n",
        "- Review concepts behind null hypothesis significance tests (NHSTs).\n",
        "\n",
        "While this live training can stand alone, I strongly encourage you to review the materials from Part 1."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jFQM2hT_xTgX",
        "colab_type": "text"
      },
      "source": [
        "## **Necessary packages**\n",
        "\n",
        "I always import my packages at the top of a notebook or `.py` file. It is good style."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "EMQfyC7GUNhT",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "sns.set()\n",
        "\n",
        "# Seed RNG\n",
        "np.random.seed(3252)\n",
        "\n",
        "# I want crisp graphics, so we'll output SVG\n",
        "%config InlineBackend.figure_format = 'svg'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "97mgKFaFxTgb"
      },
      "source": [
        "## **The Dataset**\n",
        "\n",
        "We will use the same dataset as in Part 1, but will explore different aspects of it. As a reminder, the data set we will be working with comes from a study by [Beattie, et al.](https://doi.org/10.1098/rsos.160321) in which they used the [Glasgow Facial Matching Test](https://en.wikipedia.org/wiki/Glasgow_Face_Matching_Test) (GFMT, [original paper](https://doi.org/10.3758/BRM.42.1.286)) to investigate how sleep deprivation affects a human subject's ability to match faces, as well as the confidence the subject has in those matches. Briefly, the test works by having subjects look at a pair of faces. Two such pairs are shown below.\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/datacamp/Hacker-Stats-in-Python-Live-Training/blob/master/assets/gfmt_faces.png?raw=True\" alt=\"GFMT faces\" width=\"400px\">\n",
        "</p>\n",
        "<br>\n",
        "\n",
        "\n",
        "For each of 40 pairs of faces, the subject gets as much time as he or she needs and then says whether or not they are the same person. The subject then rates his or her confidence in the choice.\n",
        "\n",
        "In this study, subjects also took surveys to determine properties about their sleep. The surveys provide three different metrics of sleep quality and wakefulness. \n",
        "\n",
        "- The Sleep Condition Indicator (SCI) is a measure of insomnia disorder over the past month. High scores indicate better sleep and scores of 16 and below indicate insomnia. \n",
        "- The Pittsburgh Sleep Quality Index (PSQI) quantifies how well a subject sleeps in terms of interruptions, latency, etc. A higher score indicates poorer sleep. \n",
        "- The Epworth Sleepiness Scale (ESS) assesses daytime drowsiness. Higher scores indicate greater drowsiness.\n",
        "\n",
        "We will explore how the various sleep metrics are related to each other and how sleep disorders affect subjects' ability to discern faces and their confidence in doing so."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "BMYfcKeDY85K"
      },
      "source": [
        "### Loading in the data set\n",
        "\n",
        "We will load in the data set as in the previous live sesson on hacker stats."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "l8t_EwRNZPLB",
        "colab": {}
      },
      "source": [
        "# Read in the dataset\n",
        "df = pd.read_csv(\n",
        "    \"https://github.com/datacamp/Hacker-Stats-in-Python-Live-Training/blob/master/data/gfmt_sleep.csv?raw=True\",\n",
        "    na_values=\"*\",\n",
        ")\n",
        "\n",
        "# Print header\n",
        "df.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s2NkyAMSxTgf",
        "colab_type": "text"
      },
      "source": [
        "For this session, we are particularly interested in the following columns:\n",
        "\n",
        "- `percent correct`: Percentage of correct responses among all trials\n",
        "- `confidence when correct`: Average confidence when the subject gave a correct response for for all trials\n",
        "- `confidence when incorrect`: Average confidence when the subject gave an incorrect response for for all trials\n",
        "- `sci`: The subject's Sleep Condition Indicator.\n",
        "- `psqi`: The subject's Pittsburgh Sleep Quality Index.\n",
        "- `ess`: The subject's Epworth Sleepiness Scale.\n",
        "\n",
        "Going forward, it will be useful to separate the subjects into two groups, insomniacs and normal sleepers. We will therefore add an `'insomnia'` column to the DataFrame with True/False entries. Recall that a person is deemed an insomniac if their SCI is 16 or below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fuf1ED8qxTgf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Add a column to the data frame for insomnia\n",
        "df['insomnia'] = df['sci'] <= 16"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1vWyV6gAxTgh",
        "colab_type": "text"
      },
      "source": [
        "With our data set in place, we can get moving with statistical inference."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eP4kzJ8FxTgi",
        "colab_type": "text"
      },
      "source": [
        "## **Review of principles from last time**\n",
        "\n",
        "In Part 1, we defined the following terms.\n",
        "\n",
        "- **Cumulative distribution function (CDF)**: A probability distribution can be defined by its CDF, given by CDF(*x*) = probability of observing a value ≤ to *x*.\n",
        "- **Generative distribution**: The unknown distribution from which the data are generated.\n",
        "- **Empirical distribution**: The distribution defined entirely by observed data.\n",
        "- **Empirical cumulative distribution function (ECDF)**: The cumulative distribution function (CDF) of the empirical distribution. It is defined as ECDF(*x*) = fraction of measured data points ≤ *x*.\n",
        "- **Plug-in principle**: Approximate the generative distribution by the empirical distribution. That is, we make the approximation CDF(*x*) ≈ ECDF(*x*).\n",
        "- **Confidence interval**: A *p*% confidence interval may be defined as follows. If an experiment is repeated over and over again, the estimate I compute will lie between the bounds of the *p*% confidence interval for *p*% of the experiments.\n",
        "- **Bootstrap sample**: If I obtained *n* measurements, a bootstrap sample is constructed by drawing *n* values out of the empirical distribution. In practice, this amounts to randomly drawing data points out of the measured data set with replacement.\n",
        "- **Bootstrap replicate**: A bootstrap replicate is the value of an estimate (e.g., a plug-in estimate) computed from a bootstrap sample.\n",
        "\n",
        "We showed that we can compute many bootstrap replicates and then compute a confidence interval for an estimate by computing percentiles of the replicates."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q2MbWXBrxTgi",
        "colab_type": "text"
      },
      "source": [
        "### **Computing plug-in estimates**\n",
        "\n",
        "The table below shows how to compute various plug-in estimates from a univariate set of measurements using NumPy.\n",
        "\n",
        "\n",
        "| Property           | Plug-in estimate      |\n",
        "|--------------------|-----------------------|\n",
        "| mean               | `np.mean(x)`          |\n",
        "| median             | `np.median(x)`        |\n",
        "| variance           | `np.var(x)`           |\n",
        "| standard deviation | `np.std(x)`           |\n",
        "| _p_ percentile     | `np.percentile(x, p)` |\n",
        "\n",
        "\n",
        "\n",
        "Two-dimensional distributions additionally have properties like covariances and correlations. For two-dimensional distributions, the plug-in estimates are:\n",
        "\n",
        "| Property           | Plug-in estimate              |\n",
        "|--------------------|-------------------------------|\n",
        "| covariance         | `np.cov(x, y, ddof=0)[0, 1]`  |\n",
        "| correlation        | `np.corrcoef(x, y)[0, 1]`     |"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FJsaDuYzxTgj",
        "colab_type": "text"
      },
      "source": [
        "### **Functions from the Part 1**\n",
        "\n",
        "We laid a lot of the groundwork we need in this live session in the previous one. We wrote functions that are useful for our hacker stats approaches to inference. First, we wrote a function to compute the *x*- and *y*-values for plotting ECDFs."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Yr76uLsRxTgk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def ecdf(data):\n",
        "    \"\"\"Compute ECDF for a one-dimensional array of measurements.\"\"\"\n",
        "    # Number of data points\n",
        "    n = len(data)\n",
        "\n",
        "    # x-data for the ECDF\n",
        "    x = np.sort(data)\n",
        "\n",
        "    # y-data for the ECDF\n",
        "    y = np.arange(1, n + 1) / n\n",
        "\n",
        "    return x, y"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-E6W9AaKxTgm",
        "colab_type": "text"
      },
      "source": [
        "We also wrote a pair of functions that enable us to draw many bootstrap replicates for a plug-in estimate."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4F61dIwixTgn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def bootstrap_replicate_1d(data, func):\n",
        "    \"\"\"Generate bootstrap replicate of 1D data.\"\"\"\n",
        "    bs_sample = np.random.choice(data, len(data))\n",
        "\n",
        "    return func(bs_sample)\n",
        "\n",
        "\n",
        "def draw_bs_reps(data, func, size=1):\n",
        "    \"\"\"Draw `size` bootstrap replicates.\"\"\"\n",
        "    # Initialize array of replicates\n",
        "    bs_replicates = np.empty(size)\n",
        "\n",
        "    # Generate replicates\n",
        "    for i in range(size):\n",
        "        bs_replicates[i] = bootstrap_replicate_1d(data, func)\n",
        "\n",
        "    return bs_replicates"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y2_1xgpdxTgp",
        "colab_type": "text"
      },
      "source": [
        "### **Automating and organizing confidence intervals**\n",
        "\n",
        "In our analysis in Part 1, we took an array of data, computed the x- and y-values for an ECDF, defined a function to compute a plug-in estimate (in our case, the mean), used the function to compute the estimate, generated bootstrap replicates of the estimate, and computed a confidence interval. We did this for two categories, normal sleepers and insomniacs.\n",
        "\n",
        "At various points in our analyses, we might like to access these results. It is therefore useful to define a class to compute and store the results of an ECDF, plug-in estimate, and confidence interval calculation for a one-dimensional data set. Many participants in this live training may not be familiar with Python's classes, so we will instead mimic a class by writing a function that does the necessary computations and returns the results in a dictionary."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SSFIPmKKxTgq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def plugin_summary(\n",
        "    data, func, ptiles=(2.5, 97.5), n_bs_reps=10000, label=None\n",
        "):\n",
        "    \"\"\"Compute and store ECDF, plug-in estimate, and confidence\n",
        "    intervals in a dictionary.\"\"\"\n",
        "    # Initialize output dictionary\n",
        "\n",
        "    # Store data and settings\n",
        "\n",
        "    # Compute ECDF x and y values\n",
        "\n",
        "    # Compute plug-in estimate\n",
        "\n",
        "    # Compute bootstrap confidence interval\n",
        "\n",
        "    return summary"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lXF4zDixxTgs",
        "colab_type": "text"
      },
      "source": [
        "Importantly, whenever we instantiate this class, we get an object that has the values for plotting an ECDF, a plug-in estimate, and a bootstrap confidence interval for that measurement. For example, to repeat the analyses from last time, we can pull out the percent correct for normal sleepers and for insomniacs and instantiate the class."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0CDUa_6QxTgs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Initialize list of plug-in summaries\n",
        "\n",
        "# Iterate through groups and instantiate conf intervals\n",
        "\n",
        "# Adjust label names to be descriptive\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zGoTkCZ2xTgv",
        "colab_type": "text"
      },
      "source": [
        "We now have a list of dictionaries that contain our analyses. Conveniently, we have the data set, what function was used to compute the plug-in estimate, the percentiles used for the confidence interval, the number of bootstrap replicates we used, which category we are considering (insomniacs vs. normal sleepers), the *x*- and *y*-values for the ECDF, the plug-in estimate for the mean, the bootstrap replicates we drew and the confidence interval. A treasure trove!\n",
        "\n",
        "Now that we have this convenient way of storing our analysis results, we should modify our plotting function from last time to plot confidence intervals that are stored in these dictionaries."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RXC8W4JtxTgv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def plot_conf_ints(summaries, palette=None):\n",
        "    \"\"\"Plot confidence intervals with estimates.\"\"\"\n",
        "    # Set a nice color palette\n",
        "    if palette is None:\n",
        "        palette = [\n",
        "            \"#1f77b4\",\n",
        "            \"#ff7f0e\",\n",
        "            \"#2ca02c\",\n",
        "            \"#d62728\",\n",
        "            \"#9467bd\",\n",
        "            \"#8c564b\",\n",
        "            \"#e377c2\",\n",
        "            \"#7f7f7f\",\n",
        "            \"#bcbd22\",\n",
        "            \"#17becf\",\n",
        "        ]\n",
        "    elif type(palette) == str:\n",
        "        palette = [palette]\n",
        "\n",
        "    labels = [ci[\"label\"] for ci in summaries][::-1]\n",
        "    estimates = [ci[\"estimate\"] for ci in summaries][::-1]\n",
        "    conf_intervals = [ci[\"conf_int\"] for ci in summaries][::-1]\n",
        "    palette = palette[: len(labels)][::-1]\n",
        "\n",
        "    # Set up axes for plot\n",
        "    fig, ax = plt.subplots(figsize=(5, len(labels) / 2))\n",
        "\n",
        "    # Plot estimates as dots and confidence intervals as lines\n",
        "    for i, (label, est, conf_int) in enumerate(\n",
        "        zip(labels, estimates, conf_intervals)\n",
        "    ):\n",
        "        color = palette[i % len(palette)]\n",
        "        ax.plot(\n",
        "            [est],\n",
        "            [label],\n",
        "            marker=\".\",\n",
        "            linestyle=\"none\",\n",
        "            markersize=10,\n",
        "            color=color,\n",
        "        )\n",
        "\n",
        "        ax.plot(conf_int, [label] * 2, linewidth=3, color=color)\n",
        "\n",
        "    # Make sure margins look ok\n",
        "    ax.margins(y=0.25 if len(labels) < 3 else 0.125)\n",
        "\n",
        "    return ax"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I-dzOTjPxTgx",
        "colab_type": "text"
      },
      "source": [
        "We can now put this class to use to make two of our plots from Part 1. We'll start with the ECDFs."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "l9dS4ogNxTgy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Dictionary connecting category to color\n",
        "colors = {\"normal\": \"#1f77b3\", \"insomniac\": \"#ff7f0e\"}\n",
        "\n",
        "# Plot the ECDFs"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ihcus2xTxTg1",
        "colab_type": "text"
      },
      "source": [
        "And now the confidence intervals."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fu4dCTufxTg1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Make plot\n",
        "\n",
        "# Set the axis label\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n6UECMqcxTg3",
        "colab_type": "text"
      },
      "source": [
        "## **Our approach**\n",
        "\n",
        "We are now all set to continue with statistical inference on this data set. In this live training, we will address the following questions (with question (1) answered in Part 1 and reviewed above).\n",
        "\n",
        "2. How different is *confidence* in facial matching for insomniacs versus normal sleepers?\n",
        "3. How are the different sleep metrics correlated?\n",
        "4. How do sleep metrics influence facial matching performance?\n",
        "\n",
        "Along the way, we will introduce the necessary new theoretical and technical concepts."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qdHPSQ4FxTg4",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "\n",
        "<center><h1> Q&A 1</h1> </center>\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bf1n2wyRxTg4",
        "colab_type": "text"
      },
      "source": [
        "## **2. Confidence of insomniacs versus normal sleepers**\n",
        "\n",
        "We have seen that there seems to be a difference in performance between normal sleepers and insomniacs. Now, let's compare how well the subjects *think* they are doing in their responses. We will investigate their confidence when correct and incorrect in the face matching tasks for both insomniacs and normal sleepers. We start by extracting the measurements from the DataFrame as NumPy arrays."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ny-4aAa-xTg5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Extract confidences for normal sleepers\n",
        "conf_corr_normal = (\n",
        "    df.loc[~df[\"insomnia\"], \"confidence when correct\"].dropna().values\n",
        ")\n",
        "conf_incorr_normal = (\n",
        "    df.loc[~df[\"insomnia\"], \"confidence when incorrect\"].dropna().values\n",
        ")\n",
        "\n",
        "# Extract confidences for insomniacs\n",
        "conf_corr_insom = (\n",
        "    df.loc[df[\"insomnia\"], \"confidence when correct\"].dropna().values\n",
        ")\n",
        "conf_incorr_insom = (\n",
        "    df.loc[df[\"insomnia\"], \"confidence when incorrect\"].dropna().values\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V4sZxFSKxTg7",
        "colab_type": "text"
      },
      "source": [
        "### **Performing the analysis (the automated way)**\n",
        "\n",
        "The analysis pipeline is much the same as we did for looking at the percent correct in the first lesson. We will again compute means and confidence intervals of the mean. We can conveniently use our `plugin_summary()` functions to do all of the calculations and give us a dictionary with the results."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mxIxMtPgxTg7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Put into dictionaries\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xIO5Zi2_xTg-",
        "colab_type": "text"
      },
      "source": [
        "### **EDA**\n",
        "\n",
        "While we did just compute confidence intervals, we should always do a bit of EDA to check things out. We will plot the ECDFs of all four classes (confidence when correct and incorrect for normal sleepers and insomniacs) and overlay them on the same plot."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "68gir3KzxTg-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Dictionary of kwargs for convenience\n",
        "\n",
        "# Plot ECDFs\n",
        "\n",
        "# Adjust labels\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RV11ed-wxThB",
        "colab_type": "text"
      },
      "source": [
        "We see a few striking trends directly from the ECDF.\n",
        "\n",
        "1. Regardless of sleep quality, subjects are less confident in their responses when they are wrong.\n",
        "2. When they are right, insomniacs and normal sleepers seem to have the same level of confidence in their responses.\n",
        "3. When they are wrong, insomniacs have more confidence in their responses than normal sleepers.\n",
        "\n",
        "This last point is really interesting. The insomniacs seem less likely to *know when they are wrong.* This kind of cognitive impairment could have real repercussions in how they go about their lives. We will give this last difference special attention."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bSzm2gnCxThB",
        "colab_type": "text"
      },
      "source": [
        "### **Plug-in estimates and confidence intervals**\n",
        "\n",
        "We have already computed the plug-in estimates for the mean and confidence intervals when we instantiated the classes. Let's compare them."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eYpmwGswxThC",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Make plot of confidence intervals\n",
        "\n",
        "# Set axis label\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eA23pR4dxThF",
        "colab_type": "text"
      },
      "source": [
        "The evidence from this plot is clear. Insomniacs are more confident that they have provided a correct answer when they have not. The difference is about 10% on the confidence scale, which is about 20% of the entire range of responses from subjects about their confidence. If we were to perform the experiment again, we would very likely continue to see this discrepancy."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gGAbg_iTxThG",
        "colab_type": "text"
      },
      "source": [
        "### Defining null hypothesis significance tests\n",
        "\n",
        "We have already drawn telling conclusions from the data from the analysis above. Considering the confidence when incorrect, the rightward shift of the ECDF for insomniacs relative to normal sleepers strongly suggests the overconfidence of insomniacs. The confidence intervals make clear that this effect is not due to the inherent variability among participants and the small sample size. We really could conclude our analysis there.\n",
        "\n",
        "Nonetheless, we will proceed to do a null hypothesis significance test. In some fields, this is referred to as an A/B test. I do so primarily because they are widely used and are important to know how to do and to interpret. I will give further commentary about why the results of NHST should not carry as much weight as conclusions derived from graphical EDA with ECDFs and calculation of confidence intervals after we perform the analysis.\n",
        "\n",
        "I begin by reviewing what a NHST is. A typical hypothesis test consists of these steps.\n",
        "\n",
        "1. Clearly state the hypothesis being considered, referred to as the **null hypothesis**.\n",
        "2. Define a **test statistic**, a scalar value that you can compute from data, almost always a statistical functional of the empirical distribution. Compute it directly from your measured data.\n",
        "3. *Simulate* data acquisition for the scenario where the null hypothesis is true. Do this many times, computing and storing the value of the test statistic each time.\n",
        "4. The fraction of simulations for which the test statistic is at least as extreme as the test statistic computed from the measured data is called the **p-value**, which is what you report.\n",
        "\n",
        "We need to be clear on our definition here. **The p-value is the probability of observing a test statistic being at least as extreme as what was measured if the null hypothesis is true.** It is exactly that, and nothing else. It is not the probability that the null hypothesis is true. In the frequentist interpretation of probability, we cannot assign a probability to the truth of a hypothesis.\n",
        "\n",
        "A complete definition of a hypothesis test then has the following components.\n",
        "\n",
        "- The null hypothesis.\n",
        "- The test statistic.\n",
        "- What it means to be at least as extreme. \n",
        "\n",
        "All of the named hypothesis tests you may have heard of, like the Student-t test, the Mann-Whitney U-test, Welch’s t-test, etc., describe a specific hypothesis with a specific test statistic with a specific definition of what it means to be at least as extreme (e.g., one-tailed or two-tailed). I can never remember what these are, nor do I encourage you to; you can always look them up. (I like [Kanji's *100 Statistical Tests*](https://dx.doi.org/10.4135/9781849208499) for this purpose.) Rather, you should just clearly write out what your test is in terms of the hypothesis, test statistic, and definition of extremeness."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tUmxw8Z_xThG",
        "colab_type": "text"
      },
      "source": [
        "### **Performing a NHST**\n",
        "\n",
        "To perform a hypothesis test using hacker stats, we will follow the prescription laid out above. We of course have to start by specifying our test.\n",
        "\n",
        "- Null hypothesis: The generative distributions for reported confidence is the same for insomniacs as it is for normal sleepers.\n",
        "- Test statistic: Mean confidence of normal sleepers minus that of insomniacs.\n",
        "- At least as extreme as: The test statistic is less than or equal to the observed difference of means.\n",
        "\n",
        "The bread and butter of conducting a hypothesis test is simulating data acquisition under the null hypothesis. For this specific hypothesis, that the distributions are the same, there is a very straight-forward way of simulating it. Let `x` and `y` be our two arrays, with `m` entries in `x` and `n` entries in `y`.\n",
        "\n",
        "1. Concatenate the two data sets into one. This is done with `np.concatenate((x, y))`.\n",
        "2. Randomly scramble the order of the combined data set. This is accomplished using `np.random.permutation()`.\n",
        "3. Designate the first `m` entries in this scrambled array to be \"`x`\" and the remaining to be \"`y`.\"\n",
        "\n",
        "Having simulated the data acquisition process, we then use these new \"`x`\" and \"`y`\" arrays to compute the test statistic.\n",
        "\n",
        "This simulation is exact; it is as if the labels of the data set (in this case \"normal sleeper\" and \"insomniac\") have no meaning; hence the distributions of the two data sets are entirely equal. Note that we are *not* say what the distribution is, only that the two data sets have the *same* distribution. \n",
        "\n",
        "A test done in this way is referred to as a **permutation test**. A **permutation sample** is akin to a bootstrap sample; it is a new pair of data sets generated after scrambling the concatenated data set. A **permutation replicate** is a value of the test statistic computed from a permutation sample, in this case the difference of means.\n",
        "\n",
        "Let's code up the functions to do the test, starting with generating the permutation sample."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yW9D8108xThG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def permutation_sample(data1, data2):\n",
        "    \"\"\"Generate a permutation sample from two data sets.\"\"\"\n",
        "    # Concatenate the data sets\n",
        "\n",
        "    # Permute the concatenated array\n",
        "\n",
        "    # Split the permuted array into two\n",
        "\n",
        "    return perm_sample_1, perm_sample_2"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WXzlxetuxThJ",
        "colab_type": "text"
      },
      "source": [
        "We need to provide a function that takes in two data sets (a permutation sample or the original data sets) and computes the test statistic. In our case, this is the difference of means."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ynHZCbCWxThK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def diff_of_means(data_1, data_2):\n",
        "    \"\"\"Difference in means of two arrays.\"\"\""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8v1WWJYqxThM",
        "colab_type": "text"
      },
      "source": [
        "Finally, we can write a function to generate multiple permutation replicates."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Pmr-q0MSxThM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def draw_perm_reps(data_1, data_2, func, size=1):\n",
        "    \"\"\"Generate multiple permutation replicates.\"\"\"\n",
        "    # Initialize array of replicates: perm_replicates\n",
        "\n",
        "    for i in range(size):\n",
        "        # Generate permutation sample\n",
        "\n",
        "        # Compute the test statistic\n",
        "\n",
        "    return perm_replicates"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5qYxo7ErxThO",
        "colab_type": "text"
      },
      "source": [
        "We'll start with the confidence when giving correct answers to the face matching test for insomniacs versus normal sleepers. We start by computing the test statistic from the actual measurements."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kQ_vzpUfxThO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Compute difference of means from original data\n",
        "\n",
        "# Take a look"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jZ5uDC5HxThQ",
        "colab_type": "text"
      },
      "source": [
        "Now, we'll draw permutation replicates."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lzcE0lb0xThQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Acquire permutation replicates\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7s9SzzpbxThT",
        "colab_type": "text"
      },
      "source": [
        "Now, we'll see what fraction were less than or equal to the test statistic calculated from the measured data sets."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Xgc7R_TpxThU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Compute p-value\n",
        "\n",
        "# Take a look\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BxDootnhxThW",
        "colab_type": "text"
      },
      "source": [
        "In this case, the p-value is about 0.15, which means that about 15% of the time, the null hypothesis that the confidence when correct for insomniacs and normal sleepers is the same will generate a difference of means at least as different in magnitude than what was observed. The data are then commensurate with the null hypothesis; they do not serve to rule out that hypothesis.\n",
        "\n",
        "Now, let's do the same analysis for the case where the subjects responded incorrectly to the face matching test."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rvJfDztaxThX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Compute difference of means from original data\n",
        "\n",
        "# Acquire permutation replicates\n",
        "\n",
        "# Compute p-value\n",
        "\n",
        "# Take a look\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gz15tZe5xThY",
        "colab_type": "text"
      },
      "source": [
        "This p-value is much smaller. Only 0.1 or 0.2% of the permutation samples had a difference of means less than or equal to what was observed. According to this result, it is unlikely that the observed data set could have been generated by identical distributions for normal sleepers and insomniacs. We already knew this from our analysis of the confidence intervals, and this is the icing on the cake, to borrow a phrase from [Allen Downey](http://www.allendowney.com/wp/)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "et05XUExxThZ",
        "colab_type": "text"
      },
      "source": [
        "### **Do not overemphasize NHSTs**\n",
        "\n",
        "It is important to know exactly what an NHST gives you, a p-value. It is the probability of observing a test statistic being at least as extreme as what was measured if the null hypothesis were true. To put interpretation of the p-value in perspective, we can compare to the interpretation of plug-in estimates and confidence intervals. Let's consider three questions addressed by these calculations.\n",
        "\n",
        "1. **Plug-in estimate**: How different are the means of the two samples?\n",
        "2. **Confidence interval**: If we were to do the experiment again, how similar would our plug-in estimates be?\n",
        "3. **NHST**: What is the probability of observing a difference in means of the two samples at least as large as the observed difference in means, if the two samples in fact have the same generative distribution?\n",
        "\n",
        "I contend that a person thinking about their dataset might naturally ask the first two questions, but the third question is convoluted and of little practical interest.\n",
        "\n",
        "To add further perspective, say we made trillions of measurements of two different samples and their mean differs by one part per million. This difference, though tiny, would still give a low p-value, and therefore often be deemed \"statistically significant.\" But, ultimately, it is the size of the difference, or the **effect size**, relative to the experiment-to-experiment variation, that we care about.\n",
        "\n",
        "Put succinctly: **Statistical significance and practical significance are not the same thing.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QatPBzdkxThZ",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "\n",
        "<center><h1> Q&A 2</h1> </center>\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y_isJGHBxThZ",
        "colab_type": "text"
      },
      "source": [
        "## **3. Correlation of sleep metrics**\n",
        "\n",
        "For the next study of the data set, we will explore the correlations among the three sleep/drowsiness metrics. Remember that two variables are **correlated** if one is high while the other is high. They are anticorrelated if one is high while the other is low. And they are said to be uncorrelated if the level of one is not related to the level of the other.\n",
        "\n",
        "We can do a quick exploration of the relationship between each pair of sleep indexes by making scatter plots of each pair of indexes. This is conveniently accomplished using seaborn's `pairplot()` function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1T58276uxTha",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "_ = sns.pairplot(df[['sci', 'psqi', 'ess']], diag_kind=None, height=2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8EmO_QE7xThb",
        "colab_type": "text"
      },
      "source": [
        "In glancing at these plots, the SCI and PSQI are clearly anticorrelated. This makes sense; they are both assess a similar feature of the subjects, the quality of nighttime sleep. Their scales are such that quality sleep is a high number in the SCI and a low number in the PSQI, so we would expect them to be anticorrelated.\n",
        "\n",
        "However, visually at least, the ESS does not seem to correlate with SCI or PQSI. Recall that ESS is a measure of daytime drowsiness. We might expect that to correlate with sleep quality, but it does not seem to, at least graphically (though the analysis we will soon do may show that it does). Previous studies (e.g., [this one](http://dx.doi.org/10.1016/j.smrv.2009.04.002)) have shown that insomniacs experience hyperarousal, meaning that they are not drowsy as a result of their poor sleep. This would suggest that we should not be surprised by a lack of correlation between SCI and ESS."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bFdE8qUExThc",
        "colab_type": "text"
      },
      "source": [
        "### **Plug-in estimate for correlation**\n",
        "\n",
        "With the graphical analysis out of the way, let's investigate the Pearson correlation between the pairs of metrics. We can compute the plug-in estimate using `np.corrcoef()`, first extracting the data as NumPy arrays. Remember that the Pearson correlation is 1 for perfect correlation, –1 for perfect anticorrelation, and zero for complete absence of correlation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SQd4LR0dxThc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Extract all three sleep/drowsiness metrics as NumPy arrays\n",
        "\n",
        "# Compute Pearson correlation between each pair\n",
        "\n",
        "# Print the result\n",
        "print('plug-in SCI-PSQI correlation:', rho_sci_psqi)\n",
        "print('plug-in SCI-ESS correlation: ', rho_sci_ess)\n",
        "print('plug-in PSQI-ESS correlation:', rho_psqi_ess)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vbMWucuQxThe",
        "colab_type": "text"
      },
      "source": [
        "Indeed, the plug-in estimates suggest that SCI and PSQI are strongly anticorrelated, but the SCI and PSQI are not really correlated with ESS."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OCsezjjGxThe",
        "colab_type": "text"
      },
      "source": [
        "### **Bootstrap confidence intervals for correlation**\n",
        "\n",
        "At face, it may seem more difficult to compute bootstrap confidence intervals for correlation, we are going to have to sample out of *two* arrays. However, we can again directly apply the plug-in principle and rely on our coding skills to do it.\n",
        "\n",
        "Consider for a moment the SCI and PSQI. Each individual subject was tested for these two sleep metrics. We cannot sample out of the array of SCI's and PSQI's independently, since they were measured *together* for each subject. Instead, we will choose a *subject* at random, and then record *both* his or her SCI and PSQI score as we build our bootstrap sample. Stated another way, we randomly choose *pairs* of scores.\n",
        "\n",
        "This procedure is called **pairs bootstrap**. To implement pairs bootstrap on a pair of arrays of equal length, we randomly sample an *index* of the arrays, and then select the pair of values that share that index. \n",
        "\n",
        "Let's write a function to draw bootstrap pairs. As with our function for drawing bootstrap replicates for one-dimensional data, we will supply a function to compute the estimate as an argument. It must have the call signature `func(x, y)`, taking two arrays as arguments and returning a scalar result."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yTrWSpnuxThf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def draw_bs_pairs(x, y, func, size=1):\n",
        "    \"\"\"Perform pairs bootstrap for single statistic.\"\"\"\n",
        "    # Set up array of indices to sample from: inds\n",
        "\n",
        "    # Initialize replicates: bs_replicates\n",
        "\n",
        "    # Generate replicates\n",
        "\n",
        "    return bs_replicates"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BDR0GhxfxThh",
        "colab_type": "text"
      },
      "source": [
        "Now that we have this function in place, we can write our function to compute the Pearson correlation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O2wBfHzbxThh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def pearson_r(x, y):\n",
        "    \"\"\"Compute Pearson correlation coefficient between two arrays.\"\"\"\n",
        "    # Compute correlation matrix\n",
        "\n",
        "    return corr_mat[0,1]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uGedb0UdxThj",
        "colab_type": "text"
      },
      "source": [
        "Ok, let's draw our bootstrap replicates! We will do it for each pair of sleep indexes, SCI–PSQI, SCI–ESS, and PSQI–ESS. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_mdFFnd7xThj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Pairs bootstrap replicates of correlation\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "npzKYDi3xThm",
        "colab_type": "text"
      },
      "source": [
        "Now that we have our replicates, we can compute the confidence intervals, again by taking percentiles as we have done in the one-dimensional case."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rQPraGF8xThm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Pairs bootstrap confidence intervals\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D336iUM0xTho",
        "colab_type": "text"
      },
      "source": [
        "To plot the results, we would like to use our `plot_conf_ints()` function. It expects a dictionary has `label`, `estimate`, and `conf_int` keys."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BVO58xjMxTho",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "sci_psqi = dict(\n",
        "    label=\"SCI–PSQI\", estimate=rho_sci_psqi, conf_int=sci_psqi_conf_int\n",
        ")\n",
        "sci_ess = dict(\n",
        "    label=\"SCI–ESS\", estimate=rho_sci_ess, conf_int=sci_ess_conf_int,\n",
        ")\n",
        "psqi_ess = dict(\n",
        "    label=\"PSQI–ESS\", estimate=rho_psqi_ess, conf_int=psqi_ess_conf_int,\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uWzn5yH2xThq",
        "colab_type": "text"
      },
      "source": [
        "Now that we have them in convenient dictionaries, we can make our plot of the confidence intervals of the correlations."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qhDa0vcNxThr",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Make plot\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bW_YHKW9xThu",
        "colab_type": "text"
      },
      "source": [
        "There is very clearly strong correlation between SCI and PSQI, but no real correlation between those metrics and the ESS metric for daytime drowsiness."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QD-_HFHCxThu",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "\n",
        "<center><h1> Q&A 3</h1> </center>\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mL-dvO8axThu",
        "colab_type": "text"
      },
      "source": [
        "## **4. Effect of sleep quality on performance**\n",
        "\n",
        "We have just investigate how the different sleep metrics are correlated. Now, we turn our attention to how the sleep metrics affect performance in the face matching test. We will start by considering the SCI. As usual, we'll do a little EDA and make a scatter plot of the percent correct versus SCI."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yjIAPwPAxThv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "_ = plt.plot(\n",
        "    df[\"sci\"], df[\"percent correct\"], marker=\".\", linestyle=\"none\", alpha=0.5,\n",
        ")\n",
        "_ = plt.xlabel(\"SCI\")\n",
        "_ = plt.ylabel(\"percent correct\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZNZnxIB-xThx",
        "colab_type": "text"
      },
      "source": [
        "If we squint, we can see there may be better scores for higher SCI. We will attempt to quantify this.\n",
        "\n",
        "We could compute the correlation between SCI and percent correct in the test, similarly to how we computed correlation between SCI and PSQI. But this is not really what we want to do. We want to know how SCI *explains* performance in the test. Contrast that with the analysis we just did in computing correlation of SCI and PSQI. We did not try to assert that one influences the other, but that the same thing, namely nighttime sleep quality, may influence them both, and hence make them correlated. In that case, a correlation was appropriate, but now, we are interested in how SCI *influences* performance. So here, a linear regression is appropriate.\n",
        "\n",
        "This subtle difference may be explained by investigating the mathematical expressions for correlation and for the slope computed from a linear regression by least squares.\n",
        "\n",
        "\\begin{align}\n",
        "\\text{Pearson correlation} = \\frac{\\sigma_{xy}}{\\sigma_x\\,\\sigma_y}, \\\\[1em]\n",
        "\\text{slope of lin. reg.} = \\frac{\\sigma_{xy}}{\\sigma_x^2}.\n",
        "\\end{align}\n",
        "\n",
        "Here, $\\sigma_{xy}$ is the covariance of $x$ and $y$, $\\sigma_x^2$ is the variance of $x$, and $\\sigma_y^2$ is the variance of $y$. The Pearson correlation is a comparison of covariance to the geometric mean of the variances of the individual variables. The slope computed from linear regression by least squares is a comparison of the covariance to the variance in the x-variable, the variable that may explain the variation in y-variable. The slope of the linear regression is thus asymmetric; we are looking at the effects of *x* on *y*, in this case the effects of SCI on percent correct."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MmDWLbkuxThx",
        "colab_type": "text"
      },
      "source": [
        "### **Performing a linear regression**\n",
        "\n",
        "When we perform a linear regression, we seek to find the line that is closest to the data in a least squares sense. This means that we want to find the line that passes through the data such that the sum of the square of the distance between the line and each data point is as small as possible.\n",
        "\n",
        "This calculation involves some beautiful mathematics, and the result is given above. The slope of the line is given by the ratio of the covariance to the variance in x. The intercept is given by the mean of the y-values minus the slope times the mean of the x-values. We can directly compute the slope and intercept using these formulas."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pmYuDyWpxThy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Extract SCI and percent correct from data frame\n",
        "sci = df['sci'].dropna().values\n",
        "pcorr = df['percent correct'].dropna().values\n",
        "\n",
        "# Compute slope\n",
        "\n",
        "# Compute intecept\n",
        "\n",
        "# Print results\n",
        "print(\"slope:     \", slope)\n",
        "print(\"intercept: \", intercept)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VzD24_6PxTh0",
        "colab_type": "text"
      },
      "source": [
        "Indeed, we see what our squinting suggested; there is a slight positive slope. For every unit increase in SCI, users get an improvement of 0.32% in their performance on the facial matching test.\n",
        "\n",
        "I pause to note here that the slope (which is of primary interest) and the intercept are numbers that are computed directly from the data. *They are therefore plug-in estimates.* All of the same principles we have been using still apply.\n",
        "\n",
        "I wanted to show you the respective formulas for the slope and intercept to make that point clear; we are \"just\" computing variances and means directly from the data and taking some ratios. In practice, we can more conveniently compute the slope and intercepts for a linear regression using the `np.polyfit()` function. It fits (*x*, *y*) data with a polynomial of arbitrary degree. A line is a polynomial of degree 1, so we can use it!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TaIEFYZQxTh1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Use polyfit to compute slope and intercept\n",
        "\n",
        "# Show that it's the same as our formulas\n",
        "print(\"slope:     \", slope)\n",
        "print(\"intercept: \", intercept)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "89wAP953xTh4",
        "colab_type": "text"
      },
      "source": [
        "### **Bootstrap replicates of linear regression parameters**\n",
        "\n",
        "Since the linear regression parameters are computed directly from data, we can again apply pairs bootstrap and compute the parameters. Since we might also be interested in the intercept, we can adapt our `draw_bs_pairs()` function to return two arrays of bootstrap replicates (one for slope and one for intercepts)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VkWM3vsAxTh4",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def draw_bs_pairs_linreg(x, y, size=1):\n",
        "    \"\"\"Perform pairs bootstrap for linear regression.\"\"\"\n",
        "    # Set up array of indices to sample from: inds\n",
        "\n",
        "    # Initialize replicates: bs_slope_reps, bs_intercept_reps\n",
        "\n",
        "    # Generate replicates\n",
        "\n",
        "    return bs_slope_reps, bs_intercept_reps"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ly1pjk6zxTh7",
        "colab_type": "text"
      },
      "source": [
        "Let's use this function to get some bootstrap replicates for the slope and intercepts!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vyFsmM0axTh8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "bs_slope_reps, bs_intercept_reps = draw_bs_pairs_linreg(sci, pcorr, size=10000)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nbcamKh5xTh_",
        "colab_type": "text"
      },
      "source": [
        "As before, we can compute the confidence intervals using `np.percentile()`. The slope is of primary interest, so we can compute its confidence interval."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QG-q0WB6xTiA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Compute confidence interval from percentiles\n",
        "\n",
        "# Take a look\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_2a-vfd6xTiC",
        "colab_type": "text"
      },
      "source": [
        "The confidence interval juuuust crosses zero, suggesting that most likely there is a real positive slope, meaning that SCI does inform performance on the test.\n",
        "\n",
        "I will visualize the linear regression on the scatter plot. To do so, I like to plot a bunch of the regression lines on the same plot as the scatter plot so I can see how the line might change if we were to do the experiment again."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IJ_HiCK5xTiC",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# x-values for lines\n",
        "\n",
        "# Plot 300 lines (thin with a bit of transparency)\n",
        "\n",
        "# Plot data\n",
        "\n",
        "# Label axes\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-aGGPrh9xTiF",
        "colab_type": "text"
      },
      "source": [
        "Indeed in this graphic, we see that on rare occasion the slope is close to zero or slightly negative, but most of the lines show positive slope."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2Zq1QOPExTiF",
        "colab_type": "text"
      },
      "source": [
        "## A NHST for a linear regression\n",
        "\n",
        "The above analysis, in which we computed a plug-in estimate for the slope, its confidence interval, and created a visualization, is sufficient, in my opinion, for this analysis. Nonetheless, we might want to perform a null hypothesis test. Remember, we have to specify three things; here is a reasonable test.\n",
        "\n",
        "1. Null hypothesis: The SCI is completely inconsequential when it comes to performance on the facial matching test.\n",
        "2. Test statistic: The slope of the linear regression line.\n",
        "3. At least as extreme as: The slope is greater than or equal to what was acquired from the measured data.\n",
        "\n",
        "Again, the trick is in simulating the null hypothesis. If the SCI has no bearing on the performance of the facial matching test, we can scramble the subjects' SCI scores, thereby decoupling them from the test scores. We can then compute the slope from the scrambled dataset. This is another version of a permutation test.\n",
        "\n",
        "Let's code up our function to draw permutation replicates of the slope."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "f0X0nn4BxTiF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def draw_perm_reps_slope(x, y, size=1):\n",
        "    \"\"\"Draw permuation replicates of the slope from lin. reg.\"\"\"\n",
        "    # Initialize array of replicates: perm_replicates\n",
        "\n",
        "    for i in range(size):\n",
        "        # Scramble x-values\n",
        "\n",
        "        # Compute the slope\n",
        "\n",
        "    return perm_replicates"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ImOctYqexTiI",
        "colab_type": "text"
      },
      "source": [
        "Now we can draw our replicates (we'll take 10,000 of them) and compute the p-value."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "d2vjJ23vxTiI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Draw permutation replicates\n",
        "perm_reps = draw_perm_reps_slope(sci, pcorr, size=10000)\n",
        "\n",
        "# p-value\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fm6duMKbxTiK",
        "colab_type": "text"
      },
      "source": [
        "Our p-value is about 0.02, suggesting that the measured data may not be commensurate with the null hypothesis being true. This again confirms what we already saw in our analysis of the plug-in estimate with confidence interval."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "poMaadwaxTiK",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "\n",
        "<center><h1> Q&A 4</h1> </center>\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MaJL2sT8xTiK",
        "colab_type": "text"
      },
      "source": [
        "## **Conclusions**\n",
        "\n",
        "We have seen that even when we have multidimensional data, the principles we apply to do nonparametric statistical inference are the same: We directly apply the plug-in principle and the frequentist interpretation of probability. We sometimes have to do bootstrapping with pairs of data, but the central ideas are the same.\n",
        "\n",
        "Null hypothesis significance testing is also possible with hacker stats. We simply directly apply the plug-in principle to simulate data acquisition under the null hypothesis. While we can perform NHSTs, they often add little insight, and we should be very careful to not over-interpret them."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r1JuRmalxTiK",
        "colab_type": "text"
      },
      "source": [
        "## **Take-home questions**\n",
        "\n",
        "There is much to still explore in this dataset, and I encourage you to do so. In particular, you might want to investigate the following:\n",
        "\n",
        "- Do the other sleep metrics, particularly the ESS which measured drowsiness, inform performance?\n",
        "- Does performance in the tests correlate with age? How about confidence?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vBYtQRz9xTiL",
        "colab_type": "text"
      },
      "source": [
        "Finally, I'll note that the [dc_stat_think package](https://github.com/justinbois/dc_stat_think) contains all of the functions that learners wrote in the Statistical Thinking courses (including Case Studies in Statistical Thinking). It is important that you know how to write all of your own functions, but this package can help you apply hacker stats in your analyses."
      ]
    }
  ]
}