METER_Sample_Size-SHORT-Quarto.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "183f8782-7b20-4c50-9b24-330b84c51dee",
   "metadata": {},
   "source": [
    "---\n",
    "title: Sample Size - How much data is enough for your experiment?\n",
    "project:\n",
    "  type: website\n",
    "format:\n",
    "  html:\n",
    "    code-fold: true\n",
    "    code-tools: true\n",
    "jupyter: python 3\n",
    "number-sections: false\n",
    "filters:\n",
    "    - pyodide\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c27513a-502c-42e2-8a0a-a079e474ba3e",
   "metadata": {},
   "source": [
    "# 1 - Just Google it?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39749e00-c96f-423f-bc0c-06010e1f07bd",
   "metadata": {},
   "source": [
    "Based on the groundbreaking research previously conducted in your lab, you and your collaborators have formulated a compelling scientific hypothesis: substance $x$ could be a genetic biomarker for longevity, potentially influencing the age at which individuals pass away. This intriguing hypothesis opens up a new frontier in our understanding of genetics and lifespan, promising significant advancements in the field.\n",
    "\n",
    "Before we can embark on an experimental journey to test the predictive power of this novel biomarker, we must first tackle a critical step: determining the appropriate sample size for a follow-up research study. The sample size is not just a number; it is a cornerstone of experimental design that ensures our data will be robust enough to support or refute our hypothesis.\n",
    "\n",
    "To accurately compute this sample size, we need to consider our prior beliefs and existing knowledge about substance $x$ and its relationship to longevity. Let's delve into the specifics. Imagine we have the following limited yet crucial pieces of information:\n",
    "\n",
    "1. **Distribution of Substance $x$:** The expression levels of substance $x$ in people follow a normal distribution.\n",
    "\n",
    "\n",
    "2. **Impact on Longevity:** Individuals at the high end of the expression spectrum tend to live approximately 5 years longer than those at the low end.\n",
    "\n",
    "Given these insights, our task is to calculate a sample size that can yield statistically significant results. This endeavor will not only help us test our hypothesis with precision but also pave the way for future research that could revolutionize our understanding of genetic influences on lifespan. Let's proceed with this vital calculation, knowing that the outcomes will bring us one step closer to potentially groundbreaking discoveries in genetic biomarkers and longevity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f39fded2-9a36-466f-8db5-696cea8b430f",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Given this information, how many individuals should we include in our study to have a reasonable chance of demonstrating this hypothesis is correct? (I.e., What is the **sample size**?)\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad9dedf5-ea2d-44bf-99be-35ca9e0e6a69",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## <b>Wait, I have no idea how to answer this?</b>\n",
    "\n",
    "- Don't worry!\n",
    "- The goal of this unit is to teach you to tackle this problem.\n",
    "- Let’s first come up with **any approach** to compute a sample size, even if we’re not confident in the results.\n",
    "\n",
    "A few possible places to start:\n",
    "\n",
    "- *Take an educated guess*: Perhaps you have taken part in or [read about](https://journals.sagepub.com/doi/pdf/10.4103/0253-7176.116232) similar research before. What order of magnitude seems right for this sort of experiment?\n",
    "\n",
    "- *Find a source*: Sample size estimation is a common topic in introductory statistics textbooks. These often include formulas that students can use to compute sample size for specific categories of questions.\n",
    "- *Google it*: There are many web-based resources (including online calculators) that are designed to enable sample size calculations. Search engines provide a starting point for finding such resource Doing so, you might end up at a website [like this](https://researchmethodsresources.nih.gov/grt-calculator) or [like this](https://www.abs.gov.au/websitedbs/d3310114.nsf/home/sample+size+calculator).\n",
    "\n",
    "Or, **if you'd like to skip this step**, we'll suggest a sample size of 100.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a779ef0d-1da7-4353-b908-296a782731d7",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Now, try to estimate the sample size using one of these approaches. What obstacles did you encounter along the way?\n",
    "\n",
    "- Sample size calculations aren't always easy or obvious, even for veteran researchers!\n",
    "\n",
    "- Here's a [good video of the challenge](https://www.youtube.com/watch?v=Hz1fyhVOjr4).\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e934f9a5-daba-4e5d-96e5-ea50c1a7a20a",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Given the description of the scientific hypothesis and experiment, think about what data you would collect and what analyses you would perform to test the hypothesis.\n",
    "\n",
    "- What types of values do you expect for each variable? What are their distributions, do you think?\n",
    "- How do you expect the variables to be related?\n",
    "- Try drawing a sketch of what you imagine a successful result might look like?\n",
    "\n",
    "- (Text) For each participant, we will collect expression levels of substance $x$ and age at death.\n",
    "- (Text) I expect age at death to increase with $x$.\n",
    "- (Multiple Choice) Show different plots of $x$ versus age at death, and ask learner to select the plot most consistent with the hypothesis.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adad8561-5a97-4175-921b-1eef24c6f18e",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## We provided very little information and asked you to compute the sample size. What other information do you think would be helpful to estimate the sample size?\n",
    "\n",
    "- \n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a9712a9-ae79-4450-9f59-68966a7e2fad",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 2- Underpowered experiments are doomed to failure."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c0bb810-3a0e-4927-afe6-2439f483de9e",
   "metadata": {},
   "source": [
    "Now that you've determined (or guessed) the sample size `N` for your experiment, let's perform the experiment.\n",
    "\n",
    "You collect `N` samples of data, so that you receive from each individual:\n",
    "\n",
    "* `x` - a measure of the proposed biomarker for longevity,\n",
    "\n",
    "* `lifespan` - the individual's age at death.\n",
    "\n",
    "Let's start by defining our choice for `N`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59b8cbf9-dd06-45fb-a1bc-2b07a41c1010",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "N = 100\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecb42cee-e624-4669-bd0a-23900ca27a08",
   "metadata": {},
   "source": [
    "With your chosen sample size `N`, you do the experiment and collect both `x` and `lifespan` for each subject."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a25fc13a-acdd-44f4-b0c4-5fb9ebfa37ac",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Load modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import pearsonr\n",
    "import statsmodels.api as sm\n",
    "# Load custom functions\n",
    "import requests\n",
    "url = \"https://raw.githubusercontent.com/Mark-Kramer/METER-Units/main/sample_size_functions.py\"\n",
    "response = requests.get(url)\n",
    "exec(response.text)\n",
    "\n",
    "x,lifespan = load_data(N)             # ... and load the data with this many N.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5a0628f-693e-43b4-ae58-32acd1b18624",
   "metadata": {},
   "source": [
    "Let's start by plotting the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb19424e-1710-428b-8b62-f5e1fe71b5c1",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.scatter(x,lifespan)\n",
    "plt.xlabel('Genetic biomarker x')\n",
    "plt.ylabel('Lifespan (years)')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9855a6b9-dc64-415c-855e-96209ee2bbb5",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## What do you observe? Does the hypothesized relationship between the biomarker $x$ and lifespan appear present in the data?\n",
    "\n",
    "- The data look like a random cloud of points.\n",
    "- It's very difficult to see the hypothesized relationship between $x$ and lifespan.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab1e61e0-832d-4dd3-b169-76721af901f8",
   "metadata": {},
   "source": [
    "Let's assess the relationship between the biomarker `x` and `lifespan` beyond visual inspection.\n",
    "\n",
    "There are many ways to do so.\n",
    "\n",
    "Here, we'll fit a line to the data and compute the slope."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b2f8ef0-f86d-46e6-9b50-feed9e3636bc",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Estimate a line from the data.\n",
    "from statsmodels.formula.api import ols\n",
    "dat                = {\"x\": x, \"lifespan\": lifespan}\n",
    "regression_results = ols(\"lifespan ~ 1 + x\", data=dat).fit()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c22000f-d414-4617-ba28-da3fa045353a",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## If this code is new to you, don't worry. Can you find the equation for the line in the code above?\n",
    "\n",
    "- The general equation for a line is `y = b + mx` where `b` is the intercept and `m` is the slope.\n",
    "- Here, we're interested in the specific line `lifespan = b + mx`.\n",
    "- In the code above, we represent this equation with the notation `lifespan ~ 1 + x`. In this notation, we tell Python to estimate the outcome variable `lifespan` as a function of a constant (with label `1` in the code) and predictor `x`. Python then estimates the solution to `linespan = b + mx` by finding the best values for `b` (the intercept) and `m` (the slope).\n",
    "- In the code above, we estimate the slope `m`, which characterizes the relationship between `lifespan` and `x`.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f09478d-a504-4492-8189-73f8e1f56218",
   "metadata": {},
   "source": [
    "Now, with the line estimated, we can print the estimated slope,  and its p-value."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcc0a446-aeec-4a2e-b0c7-dff874a2a4bf",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "slope_estimate = regression_results.params.iloc[1]\n",
    "print('Slope estimate                   = {:.3f}'.format(slope_estimate))\n",
    "print('Standard error of slope estimate = {:.3f}'.format(regression_results.bse['x']))\n",
    "print('p-value                          = {:.3f}'.format(regression_results.pvalues.iloc[1]))\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5566da0f-8cce-450c-8511-d6aa313a77e4",
   "metadata": {},
   "source": [
    "Let's interpret these numbers:\n",
    "\n",
    "### Slope estimate\n",
    "\n",
    "**Meaning:** The slope estimate represents the change in the `lifespan` for a one-unit change in the genetic biomarker `x`.\n",
    "\n",
    "**Interpretation:** For every one-unit increase the genetic biomarker `x`, the lifespan is estimated to increase by the `slope_estimate` in years, on average.\n",
    "\n",
    "### Standard error of slope estimate\n",
    "\n",
    "**Meaning:** The standard error measures the average amount that the slope estimate varies from the true slope of the population regression line. It indicates the precision of the slope estimate.\n",
    "\n",
    "**Interpretation:** A standard error of $\\sigma$, for example, suggests that the slope estimate could vary by about $2 sigma$ units from the true slope. In our case, the standard error is relatively large compared to the slope estimate. This implies there is a considerable amount of uncertainty in the estimate.\n",
    "\n",
    "### p-value\n",
    "\n",
    "**Meaning:** The p-value is used to test the null hypothesis that the slope of the regression line is zero (no relationship between `x` and `lifespan`).\n",
    "\n",
    "**Interpretation:** The p-value describes the probability of seeing an effect at least this large if substance `x` had no relation to lifespan. A p-value of 0.3, for example, is much larger than commonly used thresholds to significance levels (e.g., 0.05). In this case, this means that there is not enough evidence to reject the null hypothesis. In other words, the data do not provide sufficient evidence to conclude that there is a statistically significant relationship between `x` and `lifespan`.\n",
    "\n",
    "There's much more to say about p-values. If you're curious, check out [this link](https://github.com/Mark-Kramer/METER-Units/blob/main/METER_P_Values.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2909e622-5e33-46b7-960a-a179742fc1dd",
   "metadata": {},
   "source": [
    "Let's also **visualize** the estimated line by plotting it with the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7227b3f6-f355-4846-b3e4-e9eee95ef4a7",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "pred   = regression_results.get_prediction().summary_frame()\n",
    "mn     = pred['mean']\n",
    "ci_low = pred['mean_ci_lower'] \n",
    "ci_upp = pred['mean_ci_upper']\n",
    "\n",
    "# And plot it.\n",
    "indices_sorted = np.argsort(x,0)\n",
    "plt.figure()\n",
    "plt.scatter(x,lifespan)\n",
    "plt.plot(x[indices_sorted],mn[indices_sorted], 'r')\n",
    "plt.plot(x[indices_sorted],ci_low[indices_sorted], ':r')\n",
    "plt.plot(x[indices_sorted],ci_upp[indices_sorted], ':r')\n",
    "plt.xlabel('Genetic biomarker x')\n",
    "plt.ylabel('Lifespan (years)')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "454310ca-3aee-4aff-9b5d-50937f094df3",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Do you find a significant relationship between the genetic biomarker `x` and `lifespan`?\n",
    "\n",
    "- No. The results suggest that while there is a positive slope, indicating a potential relationship between the genetic biomarker `x` and `lifespan`, the high standard error and non-significant p-value imply that this relationship is not statistically significant. Further investigation with more data or additional variables may be needed to draw more definitive conclusions.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d929d6-99ab-4e69-8fdb-6a5b011d23e1",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Wait, this doesn't make sense!\n",
    "\n",
    "- We've applied a standard approach to compute sample size `N` and performed the experiment using this sample size.\n",
    "\n",
    "- We see a trend supporting the hypothesized relationship, but it's not significant.\n",
    "\n",
    "- Why did we fail to detect a significant relationship?\n",
    "\n",
    "- What's going on?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24935a66-c91b-464a-bcbf-c4d2a4da1dec",
   "metadata": {},
   "source": [
    "# Choose your own adventure ...\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "<b>Alert to CENTER:</b> What follows is a \"choose your own adventure\" format, in the style of the old books. The learner picks an option below to continue (either 3A, 3B, or 3C). Not all choices lead to good outcomes. This supports \"replayability\" of the Mini (e.g., learner can return and play another path).\n",
    "</div>\n",
    "\n",
    "To make sense of these confusing results, we'll use the data collected to estimate the power and improve our choice of sample size.\n",
    "\n",
    "## Which path will you pursue to compute the sample size?\n",
    "\n",
    "## **3A-** I'll resample: [Turn To Page 3A](#3A).\n",
    "\n",
    "## **3B-** I'll build models [Turn To Page 3B](#3B).\n",
    "\n",
    "## **3C-** I'll use my current sample size choice [Turn to Page 3C](#3C).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6cf4ddc-77bb-42bc-8c41-17d1fabebbec",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 3A- With resampling you can compute the sample size! {#3A}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6221559-08a6-4f9f-8ada-19ff3f3fb0c7",
   "metadata": {},
   "source": [
    "The data provided in Mini 2 represent one instantiation of the experiment, conducted with a sample size `N`. While our analysis of these data did not yield evidence to support our hypothesis, they remain extremely useful for our continued investigation into sample size. Specifically, we can leverage these data to estimate the necessary sample size for a subsequent experiment. By implementing a resampling procedure, we will systematically examine how variations in sample size `N` influence our capacity to detect a significant result, thus optimizing our experimental design for future investigations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79fc0a44-15b2-485d-8428-d3db5f8933b4",
   "metadata": {},
   "source": [
    "### Resampling procedure (Introduction)\n",
    "We're going to attempt something that seems far-fetched and magical: we'll generate new data from our existing data. To do so, we'll implement a nonparametric bootstrap to generate new pseudodata from the observed data. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "597d8dfb-cecc-4040-b0f7-ee0f3a107dc8",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "## A note on the nonparametric bootstrap.\n",
    "\n",
    "Briefly, there is strong theoretical justification for the nonparametric bootstrap. The fundamental idea is that resampling the data with replacement is equivalent to sampling new pseudodata from the empirical cumulative distribution function (eCDF) of the observed data. For a large sample of independent, identically distributed random variables, the distribution of the pseudodata generated from the eCDF will be close to the true distribution of the data. Note the important caveat that the variables are independent, identically distributed; this assumption fails in many cases, such as for time series. Here, we assume that the genetic biomarker and lifespan from each subject are drawn independently from the same distribution (i.e., the values from a subject are independent, identically distributed variables).\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afd6082f-fd51-48a7-92e0-d9fb13bd4b44",
   "metadata": {},
   "source": [
    "### Resampling procedure (4 steps)\n",
    "\n",
    "Our resampling procedure consists of 4 steps:\n",
    "\n",
    "1) Choose a new sample size (call it `N_resampled`).\n",
    "2) Draw a new (random) set of `N_resampled` labels we can use to index our data (biomarker $x$ and lifespan).\n",
    "3) Use these indices to create new pseudodata: a resampled data set.\n",
    "4) Compute the relationship (and its statistical significance) between the biomarker $x$ and lifespan in our resampled data.\n",
    "\n",
    "We'll now describe each step. For a related example, [see this video](https://youtu.be/mqDEJyW_z4c?si=heigY8z5PqAjnwKZ)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8df1c1e-f7de-4daa-b141-ff7627d48c4e",
   "metadata": {},
   "source": [
    "### Resampling procedure: Step 1\n",
    "Our first step is to choose a new sample size. Let's call it `N_resampled`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5a625f9-f886-477e-996a-e8e28068c6c4",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## What will you choose for `N_resampled`?\n",
    "\n",
    "- Our original choice of sample size resulted in a positive slope estimate, indicating a potential relationship between the genetic biomarker `x` and `lifespan`. But, the high standard error and non-significant p-value imply that this relationship is not statistically significant. Further investigation with more data may draw more definitive conclusions. However, to start, let's fix `N_resampled = N`, the original sample size.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da96f309-28e1-4dff-bb56-d92a1f8ca160",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "N_resampled = N;\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8910dbae-6d9c-4abc-a7e9-6bac89e6fa2e",
   "metadata": {},
   "source": [
    "We'll start by fixing `N_resampled = N`, the original sample size. In what follows, we'll adjust this value and examine the impact.\n",
    "\n",
    "**Thus concludes Step 1 of our resampling procedure.**\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b54ff590-56e8-4713-baca-104b5714f07d",
   "metadata": {},
   "source": [
    "### Resampling procedure: Step 2\n",
    "Our second step is to draw a random set of `N_resampled` labels to index our data (biomarker `x` and `lifespan`).\n",
    "\n",
    "To visualize this procedure, imagine we assign each patient in the original data set a number, from $0$ up to `N`. We then write each number on a marble and place all `N` marbles in an opaque bag. Each marble is assigned a unique integer value from 0 to `N`-1. Now, reach your hand into the bag, grab a marble, record its number, and replace the marble in the bag. We assume that each marble is equally likely to be selected at each draw (i.e., there are no special features that allow some marbles to be drawn more often). Repeat this procedure `N_resampled` times to create a list of `N_resampled` integers. Notice that after recording the drawn marble’s number, we replace it in the bag. So, we could potentially draw the same marble `N_resampled` times, although that’s extremely unlikely.\n",
    "\n",
    "Performing this sampling with replacement procedure by hand would, of course, be extremely time consuming (e.g., who will paint integers on each marble?). Fortunately, this is the type of boring task where a computer excels.\n",
    "\n",
    "Let's have a look at this \"marble draw\":"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86205cd2-663a-4aaf-b80b-f897567058ce",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "ind = np.random.choice(np.size(x), N_resampled)\n",
    "print(ind)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d822ceb-678f-495c-8a43-77f0d4905d03",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Look at the values in `ind`. What do they mean?\n",
    "\n",
    "- There are `N_resampled` values in the vector `ind`. That's because we draw `N_resampled` marbles.\n",
    "- These are the indices to our original data set. You can think of these as numbers indicating participants in the study (e.g., Participant 1, Participant 10, Participating 102, ...)\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "403535ef-d037-4e42-942e-b17102f3d481",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Run the code to generate `ind` again. What do you find? (I.e., is it the same or different than the first time?)\n",
    "\n",
    "- Because we draw random sets of indices, the values in `ind` will differ each time we run the code.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bbc1a59-f991-410a-8ab5-f76bf4f6569d",
   "metadata": {},
   "source": [
    "We can now generate a set of random indices to our data. We'll use these to create pseudodata in the next step.\n",
    "\n",
    "**Thus concludes Step 2 of our resampling procedure.**\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfd0c1ab-17b8-453e-8932-9fd565fb6d2d",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Resampling procedure: Step 3\n",
    "\n",
    "Our third step is to use these indices to generate the resampled data. To do so, we'll draw data from the study participants using the indices in `ind`. Again, this is a task for a computer:  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "865dae98-64ba-40cc-8685-6a1fbc0c852b",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "x_resampled        = x[ind]\n",
    "lifespan_resampled = lifespan[ind]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ca4aeab-1afc-4ebb-a9a1-c94917239330",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## How many values are in the resampled data `x_sampled` and `lifespan_resampled`?\n",
    "\n",
    "- There are `N_resampled` values in resampled data. That's because we're usig the vector `ind` to resample the data, and we drew `N_resampled` marbles.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5cedb26-d94d-494f-8e7b-e7afce4fb459",
   "metadata": {},
   "source": [
    "Let's see what those values look like, compared to our original data set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4787f5d2-6dc0-4fe2-8fca-f3b70d39755b",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.scatter(x, lifespan, color='orange', alpha=0.25, label='Original Data')\n",
    "plt.scatter(x_resampled, lifespan_resampled, color='blue', alpha=0.25, label='Resampled Data')\n",
    "plt.xlabel('Genetic biomarker x')\n",
    "plt.ylabel('Lifespan (years)');\n",
    "plt.legend()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d49311e6-9aaf-4179-9d5d-2121355d11cc",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Compare the plots of the original data (`x` and `lifespan`) with the pseudodata (`x_resampled` and `lifespan_resampled`). What do you observe? Do the pseudodata \"look like\" the original data?\n",
    "\n",
    "- The pseudodata overlaps the original data. That makes sense because we draw the pseudodata from the original data.\n",
    "- By chance, we draw some of the original data multiple times, and other original data not at all.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85a5b687-9177-4611-a221-efe41ade507b",
   "metadata": {},
   "source": [
    "We can now generate pseudodata from our original data. We'll assess the relationship between these pseudodata in the next step.\n",
    "\n",
    "**Thus concludes Step 3 of our resampling procedure.**\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c881df81-4f0c-4597-a709-a370cfae4324",
   "metadata": {},
   "source": [
    "### Resampling procedure: Step 4\n",
    "\n",
    "Our fourth step is to compute the relationship (and its statistical significance) between the resampled biomarker $x$ and resampled lifespan.\n",
    "\n",
    "To do so, we'll follow the same approach as above. We'll fit the same line to new resampled data, and again compute the slope and significance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd173d33-94a4-47cc-afa0-e376565393d9",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Estimate a line from the resampled data.\n",
    "from statsmodels.formula.api import ols\n",
    "dat_resampled                = {\"x\": x_resampled, \"lifespan\": lifespan_resampled}\n",
    "regression_results_resampled = ols(\"lifespan ~ 1 + x\", data=dat_resampled).fit()\n",
    "\n",
    "print('Slope estimate (resampled data)                   = {:.3f}'.format(regression_results_resampled.params.iloc[1]))\n",
    "print('Standard error of slope estimate (resampled data) = {:.3f}'.format(regression_results_resampled.bse['x']))\n",
    "print('p-value (resampled data)                          = {:.3f}'.format(regression_results_resampled.pvalues.iloc[1]))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b980213-c5d7-47cb-b571-925e46f0de29",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Are the slope estimate in the original data and resampled data similar or different? What about the standard errors in the estimates?\n",
    "\n",
    "- The slope estimates are similar (near 1).\n",
    "- The standard error estimates are similar (near 1).\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8814f4c9-84e1-41a4-b287-c96d064fdda5",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## So far, we've fixed `N_resampled = N`, the original sample size. Change `N_resampled` and repeat Modeling Steps 2,3,4 to generate results from multiple \"experiments\". Do you ever find a significant result? How often do the p-values you find reach your desired level of statistical significance? How does this depend on the value `N_resampled`?\n",
    "\n",
    "- Yes, now we can sometimes find p<0.05 in the modeled data when `N_resampled` is large (e.g., 1000).\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "327ae7ac-0252-4037-94ad-87779f6ae08e",
   "metadata": {},
   "source": [
    "We've now marched through the entire modeling procedure.\n",
    "\n",
    "As a final step of this procedure, we'll use this modeling approach to estimate the statsitcal power of our original experiment and a good sample size for increased power.\n",
    "\n",
    "**Thus concludes Step 4 of our resampling procedure.**\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f80f32e-4626-430b-8262-c5d7d4664df5",
   "metadata": {},
   "source": [
    "### Putting it all together\n",
    "\n",
    "We'll now use resampling to determine a good sample size for our experiment.\n",
    "\n",
    "To do so, we'll first introduce the concept of **statistical power**.\n",
    "\n",
    "In the context of statistical analysis, power and sample size are closely interrelated concepts.\n",
    "\n",
    "**Statistical Power** is the probability that a test will correctly reject a false null hypothesis (i.e., detect an effect if there is one). Higher power reduces the risk of a Type II error, where a real effect is missed (failing to reject a false null hypothesis).\n",
    "\n",
    "Our initial challenge was to compute the **sample size**: the number of observations or data points included in a study. Our initial choice `N` was too small; with this choice, we did not detect a significant relationship between the biomarker `x` and lifespan, i.e., we did not have enough **statistical power**.\n",
    "\n",
    "Using resampling, we generated pseudodata with an increased sample size `N_resampled`. Doing so in the exercise above, you might have found (sometimes) a significant relationship between the siualted biomarker `x` and simulated lifespan; you might have (sometimes) found p<0.05, the arbitrary magical threshold often used to declare a significant effect. If you'd like to understand this magic, check out [LINKS TO OTHER METERS].\n",
    "\n",
    "\n",
    "We can use this same resampling procedure to compute the statistical power of our test given the sample size. We'll do so in a few steps:\n",
    "\n",
    "### Resampling procedure to estimate the statistical power (6 steps)\n",
    "\n",
    "1) Choose a new sample size (call it `N_resampled`).\n",
    "2) Draw a new (random) set of `N_resampled` labels we can use to index our data (biomarker $x$ and lifespan).\n",
    "3) Use these indices to create new pseudodata: a resampled data set.\n",
    "4) Compute  statistical significance (p-value) between the biomarker $x$ and lifespan in our resampled data.\n",
    "5) Repeat Steps 1-5 `K` times, saving the p-value each time.\n",
    "6) The **statistical power** is the proportion of p-values below a chosen threshold `alpha`.\n",
    "\n",
    "That's a lot of steps! Let's break them down:\n",
    "\n",
    "### Resampling procedure: Steps 1-4\n",
    "\n",
    "You've already done steps 1-4 when performing resampling to create pseudodata! Nothing new to see here.\n",
    "\n",
    "### Resampling procedure: Step 5\n",
    "\n",
    "We've added Step 5, in which we create `K` new instances of the pseudodata. For each instance, we calculate and save the p-value corresponding to the statistical significance of the relationship between the biomarker $x$ and lifespan in our resampled data.\n",
    "\n",
    "At the end of Step 5, we'll have created a vector of `K` p-values. Let's do so now:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b79cb0ea-905b-443c-b0fb-04f2492d2ecd",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "N_resampled = N\n",
    "K = 1000\n",
    "p_values = np.zeros(K)\n",
    "for k in np.arange(K):                                 # For each k,\n",
    "    ind = np.random.choice(np.size(x), N_resampled)    # ... get N_resampled indices,\n",
    "    x_resampled = x[ind]                               # ... to create the pseudodata.\n",
    "    lifespan_resampled = lifespan[ind]                 # Estimate the line,\n",
    "    dat                = {\"x\": x_resampled, \"lifespan\": lifespan_resampled}\n",
    "    regression_results = ols(\"lifespan ~ 1 + x\", data=dat).fit()\n",
    "    p_values[k] = regression_results.pvalues.iloc[1]         # ... and save the p-value of the slope.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9611dc3c-700b-4d43-8260-2927531f283b",
   "metadata": {},
   "source": [
    "Let's investigate this list of p-values by plotting a historgram:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0f749e6-cd0c-4c02-a36d-42c38475c35b",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.hist(p_values, bins=np.arange(0,1,0.05));\n",
    "plt.xlabel('p-values')\n",
    "plt.ylabel('Counts')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f64041a-c25d-4e80-876d-0f0869d1bed1",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## What p-values do you observe?\n",
    "\n",
    "- P-values extend from near 0 to near 1.\n",
    "- P-values are slightly more concentrated near 0, but extend to cover the entire range.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae1fcbda-ab04-4394-bd9f-496e55221ee1",
   "metadata": {},
   "source": [
    "### Resampling procedure: Step 6\n",
    "Our last step to compute the statistical power is the proportion of p-values below a chosen threshold `alpha`.\n",
    "\n",
    "The threshold `alpha` represents the threshold for rejecting the null hypothesis when it is actually true. It's conventional to set\n",
    "\n",
    "`alpha = 0.05`\n",
    "\n",
    "which means that there is a 5% chance of committing a Type I error, which is the error of incorrectly rejecting a true null hypothesis. This  value is not inherently magical or optimal in all circumstances. But, it has become a convention primarily because it offers a middle ground that has been deemed acceptable by the scientific community for controlling Type I errors.\n",
    "\n",
    "To implement Step 6, let's compute the `statistical_power` as the proportion of times that `p_values` is less than the threshold `alpha`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1954571b-fc65-4492-a2d8-83c5f97be2d4",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "alpha = 0.05;\n",
    "statistical_power = np.sum(p_values < 0.05)/K\n",
    "print(statistical_power)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9e94f92-138d-41ff-8931-cbd1b4652ecc",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Interpret the value in `statistical_power`. What does it mean?\n",
    "\n",
    "- This value represents the proportion of times we drew pseudodata and detected a significant relationship between the biomarker `x` and lifespan. \n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7faf35c-0896-4072-a39f-c03ae5b5786d",
   "metadata": {},
   "source": [
    "The value in `statistical_power` is the **statistical power** of our test. It represents the proportion of times we reject the null hypothesis and declare a significant relationship between the biomarker `x` and lifespan.\n",
    "\n",
    "To make this graphically explicit, let's replot the histogram of `p-values` with a line at our threshold `alpha`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22dfc378-c4ef-4b68-990e-f1b43f26b5e1",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.hist(p_values, bins=np.arange(0,1,0.025));\n",
    "plt.xlabel('p-values')\n",
    "plt.ylabel('Counts')\n",
    "plt.grid()\n",
    "plt.axvline(x=0.05, color='red', label='alpha')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ede1b51f-b3a9-404c-8e55-30acc227284b",
   "metadata": {},
   "source": [
    "In this plot, the **statistical power** is the proportion of values to the left (i.e., smaller than) the red line.\n",
    "\n",
    "And that's it!\n",
    "\n",
    "The **statistical power** is not a mystical quantity. It's the probability that a test will correctly reject a false null hypothesis. And, using the data we collected, we can compute this **statistical power** for different choices of sample size (`N_resample`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3a4a3b8-3fac-4efb-bbf1-e4bfdadf32a9",
   "metadata": {},
   "source": [
    "Let's collect all the code, and perform one more experiment:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee552979-7ed6-480c-bddb-03831ee7da7f",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "N_resampled = N\n",
    "alpha       = 0.05\n",
    "K           = 1000\n",
    "p_values = np.zeros(K)\n",
    "for k in np.arange(K):                                 # For each k,\n",
    "    ind = np.random.choice(np.size(x), N_resampled)    # ... get N_resampled indices,\n",
    "    x_resampled = x[ind]                               # ... to create the pseudodata.\n",
    "    lifespan_resampled = lifespan[ind]                 # Estimate the line,\n",
    "    dat                = {\"x\": x_resampled, \"lifespan\": lifespan_resampled}\n",
    "    regression_results = ols(\"lifespan ~ 1 + x\", data=dat).fit()\n",
    "    p_values[k] = regression_results.pvalues.iloc[1]         # ... and save the p-value of the slope.\n",
    "statistical_power = np.sum(p_values < 0.05)/K\n",
    "print('Statistical power = {:.3f}'.format(statistical_power), 'for N_resampled={:.0f}'.format(N_resampled), 'and alpha={:.2f}'.format(alpha))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6074146b-1223-4ec7-bb70-701746c3d1af",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## We've used our original sample size by setting `N_resampled = N` in the code above. What is the statistical power? Does this make sense with our original conclusion in Mini 2?\n",
    "\n",
    "- Using our original sample size (`N=100`), the statistical power is small, less than 0.15.\n",
    "- The probability that the test will correctly reject a false null hypothesis (here, no relationship between biomarker $x$ and longevitiy) is 0.15. That's not very high ... \n",
    "- Therefore, with this sample size, we do not expect enough power to detect a significant effect.\n",
    "- This is consistent with the results in Mini 2, in which we failed to detect a significant effect.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d70f6bf4-9d86-4ef7-a225-f48136f2576a",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Increase `N_resampled` in the code above. What happens to the statistical power?\n",
    "\n",
    "- Statistical power increases with `N_resampled`.\n",
    "- Larger samples provide more information about the population, leading to more precise estimates of the population parameters. This precision reduces the standard error and widens the gap between the null hypothesis and the alternative hypothesis if there is a true effect, making it easier to detect significant differences. Therefore, increasing the sample size typically increases the power of a statistical test.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7817e57c-0094-4022-ae8a-dc98cd3a99c8",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## At what value of `N_resampled` does the statistical power equal 0.80?\n",
    "\n",
    "- At approximately `N_resampled` = 1000, the statistical power equals 0.80.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "202bedd7-e87f-474a-b8b9-957be3dee2f7",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Why do we choose statistical power 0.8?\n",
    "\n",
    "Choosing a statistical power of 0.8, or 80%, is a common convention in many fields of research, particularly in the social and biomedical sciences. \n",
    "\n",
    "Statistical power is the probability of correctly rejecting a false null hypothesis, thus avoiding a Type II error. A power of 0.8 means there is a 20% chance of a Type II error (failing to detect a true effect). Setting the power at 0.8 provides a reasonable balance between the risks of Type I errors (false positives) and Type II errors (false negatives). Researchers often choose a 5% (`alpha=0.05`) significance level for Type I errors, aiming to maintain a pragmatic yet cautious approach to declaring findings.\n",
    "\n",
    "Increasing power beyond 0.8 generally requires larger sample sizes, which can escalate the costs and logistical complexity of a study. The choice of 0.8 is considered a good trade-off between increasing precision and controlling operational constraints.\n",
    "\n",
    "The 0.8 level has become somewhat of a standard through historical precedent and its endorsement in statistical texts and guidelines. Researchers often follow these conventions to align with accepted practices, making their studies comparable to others in the field.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3982ebe5-a6c3-49dd-a26b-852d84e678ce",
   "metadata": {},
   "source": [
    "## Turn to Page 4 [Summary](#Summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1806aa2-a4f9-404a-981f-28bc99b578a9",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 3B- Build models to compute the sample size! {#3B}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1cdb4d-c870-4e40-ac0b-65256a10fcae",
   "metadata": {},
   "source": [
    "The data provided in Mini 2 represent one instantiation of the experiment, conducted with a sample size `N`. While our analysis of these data did not yield evidence to support our hypothesis, they remain extremely useful for our continued investigation into sample size. Specifically, we can leverage these data to estimate the necessary sample size for a subsequent experiment. By estimating models of these data, we will systematically examine how variations in sample size `N` influence our capacity to detect a significant result, thus optimizing our experimental design for future investigations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11b34096-29ed-4f1e-bc60-396a6e5c6a04",
   "metadata": {},
   "source": [
    "### Modeling procedure (Introduction)\n",
    "We're going to attempt something that seems far-fetched and magical: we'll generate new data from our existing data. To do so, we'll estimate models from the existing data, and use those models to simulate new synthetic data.\n",
    "\n",
    "### Modeling procedure (4 steps)\n",
    "\n",
    "Our modeling procedure consists of 4 steps:\n",
    "\n",
    "1) Estimate a model for biomarker `x`.\n",
    "2) Estimate a model of the relationship between biomarker `x` and lifespan.\n",
    "3) Choose a new sample size (call it `N_modeled`) and simulate data from the models.\n",
    "4) Compute the relationship (and its statistical significance) between the simulated biomarker $x$ and lifespan.\n",
    "\n",
    "We'll now describe each step."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98429b4e-4e8e-44ae-b497-af1b242fb7cc",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Modeling procedure: Step 1\n",
    "Our first step is to estimate a model for biomarker `x`.\n",
    "\n",
    "To do so, let's begin by visualizing the biomarker `x` in a historgram:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e09390de-1b6d-48e5-a1ac-86cc17842636",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.hist(x)\n",
    "plt.xlabel('Biomarker x')\n",
    "plt.ylabel('Number of observations')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394f60f8-ec1f-461a-b488-057a36dccf72",
   "metadata": {
    "tags": []
   },
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Describe the histogram. What does it look like?\n",
    "\n",
    "- The histogram has most values near 0, with fewer values near +/- 2.\n",
    "- It's approximately bell-shaped.\n",
    "- It's looks \"normal\" or \"Gaussian\".\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23d1fb7b-b0c4-4226-9d2e-2637290971d5",
   "metadata": {},
   "source": [
    "We conclude that the values for biomarker `x` look approximately [normally distributed](https://en.wikipedia.org/wiki/Normal_distribution).\n",
    "\n",
    "That's very useful, because it suggests we can model the data as normally distributed.\n",
    "\n",
    "Normal distributions are nice because we only need to estimate two parameters: **the mean and standard deviation**.\n",
    "\n",
    "Let's compute those values in Python:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "536e0d08-cfd9-4017-8971-73b2523076e0",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "# Estimating the parameters of the normal distribution for the biomarker x\n",
    "mean_x = np.mean(x)\n",
    "std_x  = np.std(x)\n",
    "print('Mean of x               = {:.2f}'.format(mean_x))\n",
    "print('Standard deviation of x = {:.2f}'.format(std_x))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ebe9990-6fbe-4099-8759-d0cbf212ce10",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Interpert the values of `mean_x` and `std_x`. How do they compare to the histogram of biomarker `x` plotted above?\n",
    "\n",
    "- The mean of the distribution for biomarker `x` is very close to zero. This suggests that on average, the values in the biomarker center around zero, consistent with our visualization of the data in the histogram.\n",
    "- The standard deviation measures the dispersion or spread of the data points around the mean. A standard deviation of 0.71 indicates that the typical deviation from the mean biomarker values is about 0.71 units. For a normal distribution, we expect about 68% of the data to fall within one standard deviation of the mean (i.e., between -0.68 and 0.74), and about 95% of the data to fall within two standard deviations (i.e., between -1.39 and 1.45). Those ranges appear are consistent with our visualization of the data in the histogram.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd81b26-3ee0-417f-8144-137944be59fc",
   "metadata": {},
   "source": [
    "Now, with these two parameters estimated, we've completely specificed our model of biomarker `x`.\n",
    "\n",
    "In words, we'll model `x` as normally distributed with mean `mean_x` = 0.03 and standard deviation `std_x` = 0.71.\n",
    "\n",
    "In Python, it's **easy to simulate values from this model.**\n",
    "\n",
    "Let's do so, and simulate 100 values from the model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cde7f6e9-7491-4b87-9721-953323cbe149",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "N_modeled = 100;\n",
    "x_modeled = np.random.normal(loc=mean_x, scale=std_x, size=[N_modeled,1])\n",
    "print(x_modeled)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6160f323-5b0a-401d-8a4f-470756f37e5c",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Compared `x_modeled` to the original values `x`. Do the values appear consistent?\n",
    "\n",
    "- Yes. In `x_modeled` there are many values near 0, and few values near +/- 2.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e17261db-8bfe-431f-93ea-98dcc86690a7",
   "metadata": {},
   "source": [
    "Let's also plot histograms for the original data `x` and the modeled data `x_modeled`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa94fdb9-6953-470e-af47-d4b5fb4ab831",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.hist(x,         bins=30, alpha=0.5, label='Original Data (x)', color='blue')\n",
    "plt.hist(x_modeled, bins=30, alpha=0.5, label='Modeled data (x_modeled)', color='red')\n",
    "plt.title('Histogram of Original Data and Samples')\n",
    "plt.xlabel('Value')\n",
    "plt.ylabel('Frequency')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01bf6ef3-a384-43f7-94fa-dd743ef4a972",
   "metadata": {
    "tags": []
   },
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Compared the histograms of `x_modeled` and the original values `x`. Do the histograms appear consistent?\n",
    "\n",
    "- Yes. In `x_modeled` there are many values near 0, and few values near +/- 2.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96a6d368-59dd-49f0-aca5-1f0563ad88bb",
   "metadata": {},
   "source": [
    "**Thus concludes Step 1 of our modeling procedure.**\n",
    "\n",
    "We've modeled the biomarker `x` as normally distributed, with mean and standard deviation estimated from our observed data.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25641e27-418c-4833-aa47-c1926979da33",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Modeling procedure: Step 2\n",
    "Our second step is to estimate a model of the lifespan.\n",
    "\n",
    "Thankfully, we've already completed this step!\n",
    "\n",
    "In Mini 2, we modeled the relationship between biomarker $x$ and lifespan as a line.\n",
    "\n",
    "Here's that code again:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38d550e2-860c-4ca9-99b3-46e9b2757e72",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "# Estimate a line from the data.\n",
    "from statsmodels.formula.api import ols\n",
    "dat   = {\"x\": x, \"lifespan\": lifespan}\n",
    "model = ols(\"lifespan ~ 1 + x\", data=dat).fit()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5707ad81-ae95-4481-b51f-7e9d05b9e26b",
   "metadata": {},
   "source": [
    "As we discussed in Mini 2, the model consists of two parameters: the `slope` and `intercept`.\n",
    "\n",
    "Let's define and print those:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b42c3a-e53d-491c-b2d8-a45ef8a3afe2",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "slope     = model.params.iloc[1]\n",
    "intercept = model.params.iloc[0]\n",
    "\n",
    "print('Slope estimate     = {:.2f}'.format(slope))\n",
    "print('Intercept estimate = {:.2f}'.format(intercept))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe2e1c4d-d4e7-4b3a-92c0-917efc496ed3",
   "metadata": {},
   "source": [
    "Do those values make sense?\n",
    "\n",
    "Let's check by plotting the data with our estimated line:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb42e2f-e603-4e58-8d82-59881953cfd8",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "# Get the model prediciton.\n",
    "pred   = model.get_prediction().summary_frame();  mn = pred['mean']\n",
    "# And plot it.\n",
    "indices_sorted = np.argsort(x,0)\n",
    "plt.figure()\n",
    "plt.scatter(x,lifespan)\n",
    "plt.plot(x[indices_sorted],mn[indices_sorted], 'r')\n",
    "plt.xlabel('Genetic biomarker x')\n",
    "plt.ylabel('Lifespan (years)')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc1dbec5-f49c-4538-846a-7381de50c140",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Examine the plot above. Do the estimated values for `slope` and `intercept` make sense?\n",
    "\n",
    "- The intercept is the value of `lifespan` when `x=0`. Looking at the plot, this occurs near 73 years, consistent with the value of intercept.\n",
    "- The line has a small upward tilt. This indicates the lifespan increases a tiny bit for each unit increase in biomarker $x$, consistent with the slope of 0.91.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a93d393-5ee7-4d6e-8a7b-8807e8efacc6",
   "metadata": {},
   "source": [
    "To simulate the model of lifespan, we'll need to extract one more parameter from the estimated model: the **dispersion**.\n",
    "\n",
    "\n",
    "The **dispersion parameter** describes the amount of uncertainty in our ability to predinct each data point. In this case, it is the residual standard deviation of the lifespan after we have tried to predict it using the expression level of substance $x$.\n",
    "\n",
    "Let's get the dispersion from the estimated model:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71542449-3785-43cf-93b4-5ae73697b809",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "dispersion = np.sqrt(model.scale)\n",
    "print('Dispersion parameter = {:.2f}'.format(dispersion))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5de058d5-6338-4b06-91a9-6383568cb005",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "## *Programming aside*\n",
    "We calculate the dispersion parameter using the `np.sqrt(model.scale)` formula from the fitted Ordinary Least Squares (OLS) model. In the context of an OLS model, the `model.scale` attribute reflects the variance of the residuals (errors), and taking the square root of this variance gives you the standard deviation.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a943c540-dc0f-403c-b613-c6d792a40f06",
   "metadata": {},
   "source": [
    "In general the dispersion (or standard deviation) tells us how much the residuals (differences between observed and predicted values) are spread out around the mean of the residuals. A higher value indicates a greater spread, suggesting more variability in the errors.\n",
    "\n",
    "In this case, a dispersion value of 7.22 means that, on average, the actual `lifespan` values deviate from the values predicted by our model by about 7.22 units. This indicates an average error magnitude of approximately 7 years from the predicted lifespan based on your model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1911a2fb-8dc4-4588-b9b9-63b8c94a3f73",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## So, the dispersion parameter indicates how well my model fits the data. Don't I want the dispersion to be 0?\n",
    "\n",
    "- In general, we do *not* expect the dispersion parameter to be 0. Our model represents a simplification of the data: we're using a simple line to capture the relationship between biomarker $x$ and lifespan. We do *not* expect this line to capture every nuance of the relationship; these are complicated biological entities with complex relationships. So, we're happy with a non-zero dispersion.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7722a7c4-20fd-45d3-8f5c-00f067969c0e",
   "metadata": {},
   "source": [
    "**Thus concludes Step 2 of our modeling procedure.**\n",
    "\n",
    "With the 3 estimated model parameters, we can now simulate values of lifespan from the model. We'll do so in the next step.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd04ade6-20a0-48e3-93fc-5ddf27f1c26b",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Modeling procedure: Step 3\n",
    "Our next step is to choose a new sample size (call it `N_modeled`) and simulate data from our models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9ecb247-9c6e-45fc-b411-ba5b6e4a00cd",
   "metadata": {
    "tags": []
   },
   "source": [
    "With the our models estiamted from the original data, we can now simulate realizations of models. \n",
    "\n",
    "To do so, we'll evaluate these model:\n",
    "\n",
    "`x_modeled        = np.random.normal(loc=mean_x, scale=std_x, size=N_modeled)`\n",
    "\n",
    "`lifespan_modeled = intercept + slope * x_modeled + np.random.normal(loc=0.0, scale=dispersion, size=N_modeled)`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b839d36-f75c-405c-9f40-a283e5b1ff52",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Describe - in words - each term in the equations above.\n",
    "- What variables do you recognize?\n",
    "- What variables are now?\n",
    "- What is the equation doing?\n",
    "    \n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "124426e7-375d-4bf3-96db-f466397676fa",
   "metadata": {},
   "source": [
    "The first equation should look familiar ... it's our model of biomarker `x`, described in Step 2.\n",
    "\n",
    "The second equation might also look familar ... it's the linear model we originally estimated data, with an added random noise component to simulate data.\n",
    "\n",
    "Here's a breakdown of each term in the second equation:\n",
    "\n",
    "- **lifespan_modeled**: This is the dependent variable in the equation, representing the predicted values of lifespan based on the model. This variable is being assigned the values calculated by the formula on the right-hand side.\n",
    "\n",
    "- **intercept**: This is the intercept of the linear model. It represents the value of the dependent variable (`lifespan_modeled`) when the independent variable (`x`) is zero.\n",
    "\n",
    "- **slope**: This term is the coefficient of the independent variable `x` in the linear model. It measures the amount by which `lifespan_modeled` is expected to increase for a one-unit increase in `x`. It represents the steepness or incline of the regression line.\n",
    "\n",
    "- **x**: This is the independent variable or predictor in the model. Here we use biomarker 'x' to predict `lifespan_modeled`. We assume the relationship between `x` and `lifespan_modeled` is linear in this model.\n",
    "\n",
    "- **np.random.normal(loc=0.0, scale=dispersion, size=[N_modeled,1])**: This\n",
    "- function generates random noise added to the linear model, simulating variability in the data that is not explained by the independent variable `x` alone:\n",
    "\n",
    "    - **loc=0.0**:  This specifies the mean of the normal distribution from which the random noise is drawn. A mean of 0 indicates that the noise is centered around zero, adding no systematic bias to the predictions, just variability.\n",
    "    - **scale=dispersion**: This is the standard deviation of the normal distribution. It controls the variability of the noise added to the model. The term dispersion here is whate we calculated above, representing the standard deviation of the residuals in the linear model.\n",
    "    - **size=[N_modeled,1]**: This specifies the shape of the array of random values generated. `N_modeled` is the number of observations or samples for which you are modeling lifespan_modeled. The [N_modeled,1] format makes the output an N_modeled-by-1 array, where each element is a random noise value added to the corresponding model prediction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51c0386c-c22d-419a-ae64-32cf83a89983",
   "metadata": {},
   "source": [
    "With the estimated models in hand, it's now simple to simulate data from the models.\n",
    "\n",
    "To do so, we must choose a value for `N_modeled`.\n",
    "\n",
    "Let's start by setting `N_modeled = N`, our original sample size, and then simulate the models. To do so, "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e7e842d-2181-4558-8d16-99230446138d",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "N_modeled = N\n",
    "x_modeled        = np.random.normal(loc=mean_x, scale=std_x, size=[N_modeled,1])\n",
    "lifespan_modeled = intercept + slope * x_modeled + np.random.normal(loc=0.0, scale=dispersion, size=[N_modeled,1])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4b24db6-119c-4722-994c-e6bebc5a90fa",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## How many values are in the simulated data `x_sampled` and `lifespan_resampled`?\n",
    "\n",
    "- There are `N_resampled` values in resampled data. That's because we're usig the vector `ind` to resample the data, and we drew `N_resampled` marbles.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "defa462c-7cf2-4c3b-8d03-a34f603e3651",
   "metadata": {},
   "source": [
    "Let's see what those modeled values look like, compared to our original data set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e4f3f70-aa4d-4830-8622-8a6777df0413",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.scatter(x, lifespan, color='orange', alpha=0.25, label='Original Data')\n",
    "plt.scatter(x_modeled, lifespan_modeled, color='blue', alpha=0.25, label='Modeled Data')\n",
    "plt.xlabel('Genetic biomarker x')\n",
    "plt.ylabel('Lifespan (years)')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7fdb95a-fdb4-432a-8d78-6ad6537b4f32",
   "metadata": {
    "tags": []
   },
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Compare the plots of the original data (`x` and `lifespan`) with the modeled data (`x_modeled` and `lifespan_modeled`). What do you observe? Do the modeled data \"look like\" the original data?\n",
    "\n",
    "- The modeled data overlap the original data. The two sets of data have appoximately the same range of biomarker values (from -2 to 2) and lifespan values (from 55 to 90 years).\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deae8ccc-349c-4272-befb-a41e3779ab91",
   "metadata": {},
   "source": [
    "**Thus concludes Step 3 of our modeling procedure.**\n",
    "\n",
    "We can now simulate data from our models. We'll assess the relationship between these simulated values in the next step.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9774444-7641-453a-8e02-b19a3f5286cc",
   "metadata": {},
   "source": [
    "### Modeling procedure: Step 4\n",
    "\n",
    "Our fourth step is to compute the relationship (and its statistical significance) between the modeled biomarker $x$ and modeled lifespan.\n",
    "\n",
    "To do so, we'll follow the same approach as above. We'll fit the same line to new modeled data, and again compute the slope and significance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5a1818b-79d4-4695-a886-01f9e07f15ef",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "# Estimate a line from the modeled data.\n",
    "from statsmodels.formula.api import ols\n",
    "dat_modeled                = {\"x\": x_modeled, \"lifespan\": lifespan_modeled}\n",
    "regression_results_modeled = ols(\"lifespan ~ 1 + x\", data=dat_modeled).fit()\n",
    "\n",
    "print('Slope estimate (modeled data)                   = {:.3f}'.format(regression_results_modeled.params.iloc[1]))\n",
    "print('Standard error of slope estimate (modeled data) = {:.3f}'.format(regression_results_modeled.bse['x']))\n",
    "print('p-value (modeled data)                          = {:.3f}'.format(regression_results_modeled.pvalues.iloc[1]))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f8e88c2-2aad-427f-8e71-4b1e11af93e0",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Are the slope estimate in the original data and modeled data similar or different? What about the standard errors in the estimates?\n",
    "\n",
    "- The slope estimates are similar (near 1).\n",
    "- The standard error estimates are similar (also near 1).\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3ad03e0-dcbb-40a5-8cc5-f0205b7eaac3",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## So far, we've fixed `N_modeled = N`, the original sample size. Change `N_modeled` and repeat Modeling Steps 2,3,4 to generate results from multiple \"experiments\". Do you ever find a significant result? How often do the p-values you find reach your desired level of statistical significance? How does this depend on the value `N_modeled`?\n",
    "\n",
    "- Yes, now we can sometimes find p<0.05 in the modeled data when `N_modeled` is large (e.g., 1000).\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47fdb971-45bc-4542-a8de-6fbb0a011014",
   "metadata": {},
   "source": [
    "**Thus concludes Step 4 of our modeling procedure.**\n",
    "\n",
    "We've now marched through the entire modeling procedure.\n",
    "\n",
    "As a final step of this procedure, we'll use this modeling approach to estimate the statsitcal power of our original experiment and a good sample size for increased power.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8195ca25-f117-4d46-84e0-50391a4cdc7b",
   "metadata": {},
   "source": [
    "### Now, let's use this modeling approach to determine a good sample size for our experiment.\n",
    "\n",
    "To do so, we'll first introduce the concept of **statistical power**.\n",
    "\n",
    "In the context of statistical analysis, power and sample size are closely interrelated concepts.\n",
    "\n",
    "**Statistical Power** is the probability that a test will correctly reject a false null hypothesis (i.e., detect an effect if there is one). Higher power reduces the risk of a Type II error, where a real effect is missed (failing to reject a false null hypothesis).\n",
    "\n",
    "Our initial challenge was to compute the **sample size**: the number of observations or data points included in a study. Our initial choice `N` was too small; with this choice, we did not detect a significant relationship between the biomarker `x` and lifespan, i.e., we did not have enough **statistical power**.\n",
    "\n",
    "Using our modeling approach, we can generate simulated data with an increased sample size `N_modeled`. Doing so in the exercise above, you might have found (sometimes) a significant relationship between the siualted biomarker `x` and simulated lifespan; you might have (sometimes) found p<0.05, the arbitrary magical threshold often used to declare a significant effect. If you'd like to understand this magic, check out [LINKS TO OTHER METERS]. \n",
    "\n",
    "We can use this same modeling procedure to compute the statistical power of our test given the sample size. We'll do so in a few steps:\n",
    "\n",
    "### Modeling procedure to estimate the statistical power (6 steps)\n",
    "\n",
    "1) Estimate a model for biomarker `x`.\n",
    "2) Estimate a model of the relationship between biomarker `x` and lifespan.\n",
    "3) Choose a new sample size (call it `N_modeled`) and simulate data from the models.\n",
    "4) Compute the relationship (and its statistical significance) between the simulated biomarker $x$ and lifespan.\n",
    "5) Repeat Steps 1-4 `K` times, saving the p-value each time.\n",
    "6) The **statistical power** is the proportion of p-values below a chosen threshold `alpha`.\n",
    "\n",
    "That's a lot of steps! Let's break them down:\n",
    "\n",
    "### Modeling procedure: Steps 1-4\n",
    "\n",
    "You've already done Steps 1-4 when creating the models to generate simulated data. Nothing new to see here.\n",
    "\n",
    "### Modeling procedure: Step 5\n",
    "\n",
    "We've added Step 5, in which we create `K` new instances of the modeled data. For each instance, we calculate and save the p-value corresponding to the statistical significance of the relationship between the biomarker $x$ and lifespan in our modeled data.\n",
    "\n",
    "At the end of Step 5, we'll have created a vector of `K` p-values. Let's do so now, with `N_modeled=N`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e6041fb-4b75-4b4f-90f9-b7af9bd5ab72",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "N_modeled = N\n",
    "K = 1000\n",
    "p_values = np.zeros(K)\n",
    "for k in np.arange(K):                                 # For each k,\n",
    "    x_modeled          = np.random.normal(loc=mean_x, scale=std_x, size=[N_modeled,1])  # Simulate the model of biomarker x\n",
    "    lifespan_modeled   = intercept + slope * x_modeled + np.random.normal(loc=0.0, scale=dispersion, size=[N_modeled,1]) # Simulate the model of lifespan\n",
    "    dat                = {\"x\": x_modeled, \"lifespan\": lifespan_modeled}\n",
    "    regression_results = ols(\"lifespan ~ 1 + x\", data=dat).fit()\n",
    "    p_values[k] = regression_results.pvalues.iloc[1]         # ... and save the p-value of the slope.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03f6e463-3314-4251-b9a5-9d768d0ab7e6",
   "metadata": {},
   "source": [
    "Let's investigate this list of p-values by plotting a historgram:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa3347cb-7507-4e81-a850-8170cc1206f4",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.hist(p_values, bins=np.arange(0,1.05,0.05));\n",
    "plt.xlabel('p-values')\n",
    "plt.ylabel('Counts')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f94f8eb7-58a8-4f10-9277-f0c8f0f7301f",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## What p-values do you observe?\n",
    "\n",
    "- P-values extend from near 0 to near 1.\n",
    "- P-values are slightly more concentrated near 0, but extend to cover the entire range.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0e88cce-91a4-4987-afdf-5861811ebf05",
   "metadata": {},
   "source": [
    "### Modeling procedure: Step 6\n",
    "Our last step to compute the statistical power is the proportion of p-values below a chosen threshold `alpha`.\n",
    "\n",
    "The threshold `alpha` represents the threshold for rejecting the null hypothesis when it is actually true. It's conventional to set\n",
    "\n",
    "`alpha = 0.05`\n",
    "\n",
    "which means that there is a 5% chance of committing a Type I error, which is the error of incorrectly rejecting a true null hypothesis. This  value is not inherently magical or optimal in all circumstances. But, it has become a convention primarily because it offers a middle ground that has been deemed acceptable by the scientific community for controlling Type I errors.\n",
    "\n",
    "To implement Step 6, let's compute the `statistical_power` as the proportion of times that `p_values` is less than the threshold `alpha`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb8f9721-bd90-4c26-b5df-ff6456517ddd",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "alpha = 0.05;\n",
    "statistical_power = np.sum(p_values < 0.05)/K\n",
    "print(statistical_power)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd795c1e-6626-43c1-95bb-b591bda307a3",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## **Question :** Interpret the value in `statistical_power`. What does it mean?\n",
    "\n",
    "- This value represents the proportion of times we created simualted data and detected a significant relationship between the biomarker `x` and lifespan. \n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cb5ea0e-545e-4e1f-8f51-479b9577d0ad",
   "metadata": {},
   "source": [
    "The value in `statistical_power` is the **statistical power** of our test. It represents the proportion of times we reject the null hypothesis and declare a significant relationship between the biomarker `x` and lifespan.\n",
    "\n",
    "To make this graphically explicit, let's replot the histogram of `p-values` with a line at our threshold `alpha`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54be73d9-aee1-4d27-ae7c-f756ea94c36f",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.hist(p_values, bins=np.arange(0,1,0.025));\n",
    "plt.xlabel('p-values')\n",
    "plt.ylabel('Counts')\n",
    "plt.grid()\n",
    "plt.axvline(x=0.05, color='red', label='alpha')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7293c7a2-94d9-4d94-a4ee-b7579f77251b",
   "metadata": {},
   "source": [
    "In this plot, the **statistical power** is the proportion of values to the left (i.e., smaller than) the red line.\n",
    "\n",
    "And that's it!\n",
    "\n",
    "The **statistical power** is not a mystical quantity. It's the probability that a test will correctly reject a false null hypothesis.\n",
    "\n",
    "And, using the data we collected, we can compute how the **statistical power** depends on sample size by changing the value of (`N_modeled`).\n",
    "\n",
    "To that end, let's collect all the code, and perform one more experiment:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b107e2bf-928b-44d7-87d8-07f1c88690ed",
   "metadata": {
    "tags": []
   },
   "source": [
    "```{pyodide-python}\n",
    "N_modeled = 1000\n",
    "alpha     = 0.05\n",
    "K         = 1000\n",
    "p_values  = np.zeros(K)\n",
    "for k in np.arange(K):                                 # For each k,\n",
    "    x_modeled          = np.random.normal(loc=mean_x, scale=std_x, size=[N_modeled,1])  # Simulate the model of biomarker x\n",
    "    lifespan_modeled   = intercept + slope * x_modeled + np.random.normal(loc=0.0, scale=dispersion, size=[N_modeled,1]) # Simulate the model of lifespan\n",
    "    dat                = {\"x\": x_modeled, \"lifespan\": lifespan_modeled}\n",
    "    regression_results = ols(\"lifespan ~ 1 + x\", data=dat).fit()\n",
    "    p_values[k] = regression_results.pvalues.iloc[1]         # ... and save the p-value of the slope.\n",
    "statistical_power = np.sum(p_values < 0.05)/K\n",
    "print('Statistical power = {:.3f}'.format(statistical_power), 'for N_modeled={:.0f}'.format(N_modeled), 'and alpha={:.2f}'.format(alpha))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6808b51-9720-4167-9ebb-c4afa53efdae",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## We've used our original sample size by setting `N_modeled = N` in the code above. What is the statistical power? Does this make sense with our original conclusion in Mini 2?\n",
    "\n",
    "- Using our original sample size (N=100), the statistical power is small, near 0.15.\n",
    "- Therefore, with this sample size, we do not expect enough power to detect a significant effect.\n",
    "- This is consistent with the results in Mini 2, in which we failed to detect a significant effect.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b82e91fd-dac5-44b3-a784-c70553487ff2",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Now, use this code to determine the value at which `N_modeled` produces statistical power equal to 0.80?\n",
    "\n",
    "- At approximately `N_modeled` = 1000, the statistical power equals 0.80.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1577180f-7254-406d-898a-5cd6c632ffaf",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Why do we choose statistical power 0.8?\n",
    "\n",
    "Choosing a statistical power of 0.8, or 80%, is a common convention in many fields of research, particularly in the social and biomedical sciences. \n",
    "\n",
    "Statistical power is the probability of correctly rejecting a false null hypothesis, thus avoiding a Type II error. A power of 0.8 means there is a 20% chance of a Type II error (failing to detect a true effect). Setting the power at 0.8 provides a reasonable balance between the risks of Type I errors (false positives) and Type II errors (false negatives). Researchers often choose a 5% (`alpha=0.05`) significance level for Type I errors, aiming to maintain a pragmatic yet cautious approach to declaring findings.\n",
    "\n",
    "Increasing power beyond 0.8 generally requires larger sample sizes, which can escalate the costs and logistical complexity of a study. The choice of 0.8 is considered a good trade-off between increasing precision and controlling operational constraints.\n",
    "\n",
    "The 0.8 level has become somewhat of a standard through historical precedent and its endorsement in statistical texts and guidelines. Researchers often follow these conventions to align with accepted practices, making their studies comparable to others in the field.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b8b1a5f-50f2-4a99-9020-7b3b40cfa4f0",
   "metadata": {},
   "source": [
    "## Turn to Page 4 [Summary](#Summary).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d11fd60-6fdf-4d67-adc3-c32dca1da02e",
   "metadata": {},
   "source": [
    "# 3C- Do nothing, I'm happy with the current sample size choice {#3C}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0ed16ff-9727-475a-a9dc-353f63a1aed0",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Load modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import pearsonr\n",
    "import statsmodels.api as sm\n",
    "# Load custom functions\n",
    "import requests\n",
    "url = \"https://raw.githubusercontent.com/Mark-Kramer/METER-Units/main/sample_size_functions.py\"\n",
    "response = requests.get(url)\n",
    "exec(response.text)\n",
    "\n",
    "N = 100\n",
    "do_nothing_function(N)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44b00fc2-63fb-440c-b7cf-53cd36a52dd6",
   "metadata": {},
   "source": [
    "# 4- Summary {#Summary}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "354c3894-1abb-453d-ac24-6724d0894611",
   "metadata": {},
   "source": [
    "We've done something remarkable.\n",
    "\n",
    "We began with `N` observations in our original data set. Analyzing these data, we failed to detect a significant relationship between the biomarker $x$ and lifespan.\n",
    "\n",
    "Rather than abandoning these data, we instead repurposed these data to perform a sample size calculation.\n",
    "\n",
    "Depending on the adventure you choose, you (hopefully) used the original data to:\n",
    "\n",
    "- 3A: Resample the data to create pseudodata with different sample sizes `N_resample`.\n",
    "\n",
    "- 3B: Model the data to create simualted data with different sample sizes `N_modeled`.\n",
    "\n",
    "We repeated this procedure to compute the **statistical power**: the proportion of times we reject the null hypothesis.\n",
    "\n",
    "This procedure provides a more direct, intuitive approach to computing the **statistical power**.\n",
    "\n",
    "If you have some data, you can compute the sample size!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b00a662-ee04-4fda-9c52-64dd7aa1947b",
   "metadata": {},
   "source": [
    "::: {.callout-note}\n",
    "\n",
    "## Okay, I get it. When I have some data, I can use it to estimate a sample size for a future experiment But what if I do *not* have any data. Then how would I compute the sample size?\n",
    "\n",
    "- Good question. A meanigful sample size calculation requires:\n",
    "  1. A hypothesis, and\n",
    "  2. Knowledge of your data.\n",
    "\n",
    "Before planning your experiment, you'll certainly know *something* about your data. The challenge then is to use that *something* to estimate sample size.\n",
    "\n",
    "If you're interested in learning more about this scenario, check out the full Unit. \n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b36be1a-d4ed-443d-a883-7c8bdd24a479",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "*Conclusions*:\n",
    "\n",
    "- Beginning with some data, you can use resample or simulate the data to estimate the statistical power of your test. \n",
    "- You can apply this process to compute the statistical power for different sample sizes.\n",
    "- You can apply this process to determine the sample size that acheives a desired statistical power (typically 0.8).\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e35c7c2e-0d4a-4236-94d1-21ca374e6df2",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "<b>Warning:</b>\n",
    "</p>\n",
    "In this example, we were lucky that the initial draw of a small sample size produced the expected effect. An unlucky sample may have produced (by chance) an opposite effect. In that case, resampling will not produce meaningful power/sample size results. Preliminary data is often important for future experimental design, but it’s important to consider how variability in a small, preliminary dataset can influence power and sample size estimates.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}