{ "cells": [ { "cell_type": "markdown", "id": "f92d65af-a7b7-4a5c-8166-1a547773dca1", "metadata": {}, "source": [ " Kepler Mears (2026)\n", "\n", " This notebook goes through \n", " 1. Notes on scoring data and p-value definitions: How we go from the definition of a p-value in class to something we can use on real data\n", " 2. The distribution of p-values under the null: This gives a concrete example of what was said in lecture and why we shouldn't trust a p-value on its own\n", " 3. The distribution of p-values under a true effect: This provides a simulation to tweak parameters to generate data and see how the distribution changes as a result\n", " 4. P-hacking: A cool meta example of p-values in the wild and how to detect potential cheating" ] }, { "cell_type": "code", "execution_count": null, "id": "55297221-c6bb-4676-bf7c-3afc77f05647", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import random" ] }, { "cell_type": "markdown", "id": "7f04d46c-1975-4910-be7c-55747151bf7e", "metadata": {}, "source": [ "# The set up\n", "\n", "So somehow you hacked into canvas and saw the grades from the a recent pset from this year and the previous time the course was taught (pset 4 to be specific). You are curious what the data looks like so you plot it and notice that the most recent pset grades (i.e. yours) are lower. You want to prove to Sean that something is fishy here and want to make some claim about the likelihood of observing such a difference between data, so naturally you reach for p-values." ] }, { "cell_type": "code", "execution_count": null, "id": "65990e17-69eb-4b55-95d1-1cfd9d077c0c", "metadata": {}, "outputs": [], "source": [ "pset4_2024 = np.array([9.5, 9.1, 9.1, 9.5, 8.6, 9. , 9.5, 8.7, 9.1, 9.1, 8.7, 8.2, 8.7,\n", " 9.3, 9.4, 9.1, 9.6, 9.4, 9.4, 9.4, 8.5, 9.1, 8.8, 8.8, 9.2, 8.1,\n", " 8.7, 8.5, 8.5, 8.3])\n", "\n", "pset4_2026 = np.array([9.3, 8.8, 8.2, 8.3, 9.5, 8.6, 9.3, 9. , 8.8, 9.3, 8.3, 9.2, 8.2,\n", " 8.3, 7.7, 8.3, 9. , 9.2, 8.6, 8.5, 9.2, 9.1, 7.8, 9.2, 8.5, 8.1,\n", " 9. , 9.1, 8.2, 8.7])\n", "\n", "\n", "plt.figure(figsize=(4,5))\n", "\n", "# Add jitter to make it easier to see the points\n", "x1 = np.random.normal(1-0.15,0.03,len(pset4_2024))\n", "x2 = np.random.normal(1+0.15,0.03,len(pset4_2026))\n", "\n", "plt.scatter(x1, pset4_2024, alpha=0.25, color='red')\n", "plt.scatter(x2, pset4_2026, alpha=0.25, color='blue')\n", "\n", "# Set x-ticks to label cohorts\n", "plt.xticks([1-0.15, 1+0.15], [\"2024\", \"2026\"])\n", "\n", "plt.xlabel(\"Cohort\")\n", "plt.ylabel(\"Score\")\n", "plt.title(\"Raw Student Scores by Cohort\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "56ceab55-5a9f-4b7f-a81a-99817f2ba228", "metadata": {}, "source": [ "# 1. Notes on scoring data and p-value definitions\n", "\n", "As a reminder **A p-value is the probability that we would have gotten a result at least this extreme, if the null hypothesis is true.** We can determine this by calculating the CDF, which conveniently is already included in scipy stats. However there is no one CDF to rule them all and one p-value calculation to rule them all, it depends on **what you are measuring** and **what the null is**. Not to mention whatever distribution you are assuming your data is under. Here we are treating the data a Guassian but the same definition applies. Check the notes from last year to see what a binomial looks like!" ] }, { "cell_type": "markdown", "id": "44334590-757f-42fb-b50d-a73a10e7bca8", "metadata": {}, "source": [ "## 1. Z-Score: Single Observation\n", "\n", "From lecture, we decsribed finding p-values for Gaussian's by calculating the Z-score which is a way to normalize Gaussian's:\n", "\n", "$$Z = \\frac{x - \\mu}{\\sigma}$$\n", "\n", "p-value calculation where z is a specific observed score:\n", "\n", "$$P(Z \\geq z) = 1 - \\text{CDF}\\left(z\\right) \\quad$$\n", "\n", "Two-tailed p-value:\n", "\n", "$$P(|Z| \\geq z) = 2\\left(1 - \\text{CDF}\\left(z\\right)\\right)$$\n", "\n", "CDF is the standard normal $N(0,1)$.\n", "\n", "| Variable | Meaning |\n", "|---|---|\n", "| $x$ | A single observed value |\n", "| $\\mu$ | Known population mean |\n", "| $\\sigma$ | Known population standard deviation |\n", "\n", "**Use when:** You have one observation and you know the true population parameters (which we never really do)." ] }, { "cell_type": "markdown", "id": "d324bc0d-dfa8-48b4-ad7c-24e8575e58ec", "metadata": {}, "source": [ "## 2. Z-Score: Sample Mean\n", "\n", "The previous Z score is idealized when we know $$\\mu$$ and $$\\sigma$$. But we often have to estimate these parameters as well. Estimmating means is easy so we can redefinte our score where the observed value is a mean, which will adjust the variance:\n", "\n", "$$Z = \\frac{\\bar{X} - \\mu_0}{\\sigma / \\sqrt{n}}$$\n", "\n", "p-value calculation is the same except z would be some observed mean:\n", "\n", "$$P(Z \\geq z) = 1 - \\text{CDF}\\left(z\\right) \\quad$$\n", "\n", "Two-tailed p-value:\n", "\n", "$$P(|Z| \\geq z) = 2\\left(1 - \\text{CDF}\\left(z\\right)\\right)$$\n", "\n", "CDF is the standard normal $N(0,1)$. Requires known $\\sigma$.\n", "\n", "| Variable | Meaning |\n", "|---|---|\n", "| $\\bar{X}$ | Sample mean (average of $n$ observations) |\n", "| $\\mu_0$ | Hypothesized population mean |\n", "| $\\sigma / \\sqrt{n}$ | SE of the mean — $\\sigma$ is still known, but shrinks by $\\sqrt{n}$ |\n", "\n", "**Use when:** You have a sample mean and you know the true $\\sigma$ (which we also rarely do)." ] }, { "cell_type": "markdown", "id": "83c06ca4-e4c1-4704-94bb-088f61f3658b", "metadata": {}, "source": [ "## 3. T-Score: Sample Mean (Unknown $\\sigma$)\n", "\n", "However we also usually have to estimate the varaince as well. To do this we change our metric to a T score which instead of the normal/Gaussian distribution, follows a T distribution which is the subject of the p-set. The T distribution is similar and approximates a normal with enough data. But it is not the same as the process of estimating $$\\sigma$$ alter the distribution. I have outlined the changes but to see a full derivation check the links on the notes from the previous year.\n", "\n", "$$T = \\frac{\\bar{X} - \\mu_0}{S / \\sqrt{n}}$$\n", "\n", "Because we don't know $\\sigma$ we use the sample standard deviation $S$ which you can look up in greater detail in the previous year's section notes:\n", "\n", "$$S = \\sqrt{\\frac{1}{n-1}\\sum_{i=1}^n (X_i - \\bar{X})^2}$$\n", "\n", "p-value calculation is the same concept but the CDF is for a different distribution. The T distribution is also sensitive to the degrees of freedom, which itself is also commonly misintepreted. The official definition is how many independent pieces of information you have left to estimate variability after accounting for constrains. What this means for us is that for every variable we estimate, we lose a degree of freedom as we have constrained our final analysis:\n", "\n", "$$P(T \\geq t) = 1 - \\text{CDF}_{t(n-1)}\\left(t\\right) \\quad \\text{where} \\quad t = \\frac{\\bar{X} - \\mu_0}{S / \\sqrt{n}}$$\n", "\n", "Two-tailed:\n", "\n", "$$P(|T| \\geq t) = 2\\left(1 - \\text{CDF}_{t(n-1)}\\left(t\\right)\\right)$$\n", "\n", "CDF is the $t$-distribution with $n - 1$ degrees of freedom.\n", "\n", "| Slot | Meaning |\n", "|---|---|\n", "| $\\bar{X}$ | Sample mean |\n", "| $\\mu_0$ | Hypothesized population mean |\n", "| $S$ | **Estimated** standard deviation — replaces known $\\sigma$ |\n", "| $S / \\sqrt{n}$ | SE of the mean, estimated from data |\n", "\n", "Follows a $t$-distribution with $n - 1$ degrees of freedom.\n", "\n", "**Use when:** You have a sample mean but $\\sigma$ is unknown and must be estimated from the data." ] }, { "cell_type": "markdown", "id": "04adad11-be6d-4aed-87d5-65bf38639ea1", "metadata": {}, "source": [ "## 4. T-Score: Difference of Two Means\n", "\n", "Now we have to adjust our score to deal with the difference of means as our metric of interest. I am re-creating a t-test here but obviuosly omitting some details as this is not intended to be a complete guide. I largely consulted wikipedia to make sure I am feeding you correct information so don't get mad at me if it is wrong (https://en.wikipedia.org/wiki/Student%27s_t-test).\n", "\n", "$$T = \\frac{(\\bar{X}_A - \\bar{X}_B) - 0}{S_{pooled}\\sqrt{\\frac{1}{n_A} + \\frac{1}{n_B}}}$$\n", "\n", "Because we are subtracting two values, our null hypothesis is that the true means of both groups is the same. So when we great the difference statistic ($(\\bar{X}_A - \\mu_0) - (\\bar{X}_B - \\mu_0)$ the $\\mu_0$ drops out. Because under $H_0$** there is no difference between the groups the \"mean\" term should be zero! So at least here we have one aspect where the calculation is simpler and we don't have to worry about our estimatation of $\\mu_0$.\n", "\n", "Again we don't know $\\sigma$ so we have to use an estimator. This derivation is also a bit more complex and requires knowledge of variance rules so again you will have to trust me. For the difference in mean between two sample we use the Standard Error or $SE$ which one could straightforwardly derive from variance rules: \n", "\n", "$$SE = \\sqrt{\\frac{S_A^2}{n_A} + \\frac{S_B^2}{n_B}}$$\n", "\n", "However when we know (or assume) that the variance between the two sample populations is the same we use yet a more specified version:\n", "\n", "$$SE = S_{pooled}\\sqrt{\\frac{1}{n_A} + \\frac{1}{n_B}}$$\n", "\n", "Where:\n", "\n", "$$S_{pooled}^2 = \\frac{(n_A - 1)S_A^2 + (n_B - 1)S_B^2}{n_A + n_B - 2}$$\n", "\n", "\n", "When we assume the variances are equal between populations we use the \n", "\n", "p-value calculation\n", "\n", "$$P(T \\geq t) = 1 - \\text{CDF}_{t(n_A + n_B - 2)}\\left(t\\right) \\quad \\text{where} \\quad t = \\frac{(\\bar{X}_A - \\bar{X}_B)}{S_{pooled}\\sqrt{\\frac{1}{n_A} + \\frac{1}{n_B}}}$$\n", "\n", "Two-tailed:\n", "\n", "$$P(|T| \\geq t) = 2\\left(1 - \\text{CDF}_{t(n_A + n_B - 2)}\\left(t\\right)\\right)$$\n", "\n", "CDF is the $t$-distribution with $n_A + n_B - 2$ degrees of freedom.\n", "\n", "| Slot | Meaning |\n", "|---|---|\n", "| $\\bar{X}_A - \\bar{X}_B$ | Observed difference in sample means |\n", "| $0$ | Expected difference under $H_0$ (no effect) |\n", "| $S_{pooled}$ | Pooled estimate of shared $\\sigma$ across both cohorts |\n", "| $S_{pooled}\\sqrt{\\frac{1}{n_A} + \\frac{1}{n_B}}$ | SE of the difference — uncertainty accumulates from both cohorts |\n", "\n", "Follows a $t$-distribution with $n_A + n_B - 2$ degrees of freedom.\n", "\n", "**Use when:** You want to test whether two independent samples have different means, and $\\sigma$ is unknown but assumed equal across groups." ] }, { "cell_type": "markdown", "id": "d0fd10bc-5b17-48ff-847f-a877ca8a0b06", "metadata": {}, "source": [ "## What Changes Across All Four\n", "\n", "| | Score | CDF |\n", "|---|---|---|\n", "| Z, single obs | $\\frac{x - \\mu}{\\sigma}$ | $N(0,1)$ |\n", "| Z, sample mean | $\\frac{\\bar{X} - \\mu_0}{\\sigma/\\sqrt{n}}$ | $N(0,1)$ |\n", "| T, sample mean | $\\frac{\\bar{X} - \\mu_0}{S/\\sqrt{n}}$ | $t(n-1)$ |\n", "| T, diff of means | $\\frac{\\bar{X}_A - \\bar{X}_B}{S_{pooled}\\sqrt{1/n_A + 1/n_B}}$ | $t(n_A + n_B - 2)$ |\n", "\n", "The numerator is always **observed minus expected under $H_0$**. \n", "The denominator is always the **standard error of whatever is in the numerator**.\n", "\n", "In every case the p-value has the same form: $2(1 - \\text{CDF}(|\\text{score}|))$ for a two-tailed test. \n", "What changes is how the score is computed and which distribution's CDF you evaluate it against.\n", "\n", "All of these details are conveniently wrapped in scipy stats, but I wanted to be slightly more specific than just saying \"trust me bro\" and now we have everything we need to actually calculate some p-values!\n", "\n", "You assume that the null would be the same score on a pset, so the difference of means being 0. You also assume that the variance in grades to be relatively similar. Just to be sure you decide to also calculate p-values with a permutation test (also known as order statistics from lecture). Order statistics is how we get the analogy that p-values are false positive rates. The code for this is quite straightforward and doesn't require any fancy math (which the probably only did because they didn't have computers back then). Under the null hypothesis there is no difference between the data sets, i.e. the labels are meaningless, so to get a p-value you merge the sets and shuffle them to mix the labels creating two new datasets. Then you compute the difference of means and do this over and over and at the end ask how many are greater than the one you observed. This is often called a permutation test and it is my favorite way to calculate p-value like statistics because I find it intuitive and broadly applicable, this approach doesn't depend on the underlying distribution." ] }, { "cell_type": "code", "execution_count": null, "id": "191dcc8f-513e-4ae3-82dd-646bba803843", "metadata": {}, "outputs": [], "source": [ "def pvalue_analytical(a, b):\n", " \"\"\"\n", " Calculates a p-value using a two tailed T-test.\n", "\n", " Steps:\n", " 1. Compute the observed means for each group\n", " 2. Compute the Bessel-corrected sample variance for each group\n", " 3. Compute the pooled-variance and standard error of the difference in means\n", " 4. Generate the T-score\n", " 5. Use scipy stats to determine the p-value\n", " \"\"\"\n", "\n", " n1 = len(a)\n", " n2 = len(b)\n", "\n", " # Step 1: sample means\n", " mean1 = np.mean(a)\n", " mean2 = np.mean(b)\n", "\n", " # Step 2: Individual Bessel-corrected sample variances\n", " s1_sq = np.sum((a - mean1)**2) / (n1 - 1)\n", " s2_sq = np.sum((b - mean2)**2) / (n2 - 1)\n", " \n", " # Step 3: Pooled-variance (assuming equal variance)\n", " sp_sq = ((n1 - 1)*s1_sq + (n2 - 1)*s2_sq) / (n1 + n2 - 2)\n", " # standard error of difference in means\n", " se = np.sqrt(sp_sq * (1/n1 + 1/n2))\n", "\n", " # Step 4: t score - making it more explicit to show it is a difference of means compared to the null\n", " t = ( (mean1 - mean2) - 0)/ se\n", " # degrees of freedom\n", " df = n1 + n2 - 2\n", "\n", " # Step 5: two-tailed p-value\n", " p = 2 * stats.t.sf(abs(t), df)\n", "\n", " return p" ] }, { "cell_type": "code", "execution_count": null, "id": "7d95170d-569d-4877-a180-64e5e8c03b0f", "metadata": {}, "outputs": [], "source": [ "def pvalue_order_statistics(a, b, n_synthetic=5000, random_seed=None):\n", " \"\"\"\n", " Calculates a p-value using order statistics (permutation test).\n", "\n", " Core idea: if there is truly no difference between cohort A and cohort B\n", " (the null hypothesis), then the labels \"A\" and \"B\" are arbitrary -- we\n", " can shuffle them and the distribution of the test statistic shouldn't\n", " change. We use this to empirically build the null distribution, rather\n", " than assuming it follows a t-distribution analytically.\n", "\n", " Steps:\n", " 1. Compute the observed test statistic on the real data\n", " 2. Pool all observations together (labels are meaningless under H0)\n", " 3. Repeatedly shuffle the labels and recompute the statistic\n", " 4. The p-value is the fraction of shuffles that produced a statistic at least as extreme as the one we observed\n", " \"\"\"\n", " \n", " rng = np.random.default_rng(random_seed)\n", "\n", " n1 = len(a)\n", "\n", " # Step 1: Observed test statistic\n", " observed_stat = np.mean(a) - np.mean(b)\n", "\n", " # Step 2: Pool all observation\n", " # Under H0, the group labels carry no information so all observations are drawn from the same distribution.\n", " pooled = np.concatenate([a, b])\n", "\n", " # Step 3: Shuffle labels and recompute statistic a whole bunch of times (why we need computer for this)\n", " null_stats = np.zeros(n_synthetic)\n", " for i in range(n_synthetic):\n", " shuffled = rng.permutation(pooled)\n", " synthetic_a = shuffled[:n1]\n", " synthetic_b = shuffled[n1:]\n", " null_stats[i] = np.mean(synthetic_a) - np.mean(synthetic_b)\n", "\n", " # Step 4: P-value calculation. For a two-tailed count how often the null statistic is at least as extreme as the observed statistic in either direction.\n", " # You can easilt change this to a one tailed by removing the absolute value.\n", " n_as_extreme = np.sum(np.abs(null_stats) >= np.abs(observed_stat))\n", " p = n_as_extreme / n_synthetic\n", "\n", " return p" ] }, { "cell_type": "code", "execution_count": null, "id": "b20fcb8d-0ad4-49ba-8daa-e52640951d92", "metadata": {}, "outputs": [], "source": [ "print(pvalue_analytical(pset4_2024, pset4_2026))\n", "print(pvalue_order_statistics(pset4_2024, pset4_2026))" ] }, { "cell_type": "markdown", "id": "a749215d-cee7-41e5-9eb4-8c57d9f33576", "metadata": {}, "source": [ "HAH so the p-value is less than 0.05! So we rejct the null hypothesis and thus Sean (and the TFs) are definitely up to something! \n", "\n", "You present your data to Sean to which he says \"well this is just ONE pset out of the 12, there is a 5% of getting a pvalue <0.05 under the null so one out of 12 scoring significant is reasonable to expect. We teach the same material each year and draw from a similar population for TFs and students, each pset is essentially a repeated experiment and subject to random variation. You will have to do better to convince me of bias.\"\n", "\n", "So you patiently wait until the end of the year and again hack into canvas to rip the scores from 2024 and all the scores from 2026. You expect so see that a bunch of the psets the 2024 cohort scores higher on with a p-value <0.05" ] }, { "cell_type": "markdown", "id": "1b85b6c8-831c-4298-8949-5063860e451a", "metadata": {}, "source": [ "# 2. The distribution of p-values under the null" ] }, { "cell_type": "code", "execution_count": null, "id": "8d84240f-d4a5-4a8d-9a32-3c28decd69e5", "metadata": {}, "outputs": [], "source": [ "# Read in previous years data\n", "scores_2024 = np.loadtxt(\"2024_pset_scores.txt\", delimiter=',')\n", "scores_2026 = np.loadtxt(\"2026_pset_scores.txt\", delimiter=',')" ] }, { "cell_type": "code", "execution_count": null, "id": "87108cc8-caac-4b42-860f-8321c6e0fe01", "metadata": {}, "outputs": [], "source": [ "# ---------- PLOT RAW SCORES ----------\n", "plt.figure(figsize=(12,5))\n", "\n", "for i,(c1,c2) in enumerate(zip(scores_2024, scores_2026)):\n", " x1 = np.random.normal(i+1-0.15,0.03,len(c1))\n", " x2 = np.random.normal(i+1+0.15,0.03,len(c2))\n", "\n", " plt.scatter(x1,c1,alpha=0.25,label=\"2024 scores\" if i==0 else \"\", color='red')\n", " plt.scatter(x2,c2,alpha=0.25,label=\"2026 scores\" if i==0 else \"\", color='blue')\n", "\n", " # Compute means for each cohort\n", " mean1 = np.mean(c1)\n", " mean2 = np.mean(c2)\n", " \n", " # Plot horizontal mean lines for each cohort individually\n", " plt.hlines(mean1, i+1-0.25, i+1-0.05, colors='black', linestyles='dashed', linewidth=2)\n", " plt.hlines(mean2, i+1+0.05, i+1+0.25, colors='black', linestyles='dashed', linewidth=2)\n", "\n", "plt.xlabel(\"Problem Set\")\n", "plt.ylabel(\"Score\")\n", "plt.title(\"Raw Student Scores by Cohort\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "10b6dde6-df8a-4270-b648-83e521340699", "metadata": {}, "outputs": [], "source": [ "# Calculate p-values for difference in means for each p-set\n", "p_values_analytical = []\n", "p_values_order = []\n", "p_values_z = []\n", "\n", "for c1, c2 in zip(scores_2024, scores_2026):\n", " p_values_analytical.append(pvalue_analytical(c1, c2))\n", " p_values_order.append(pvalue_order_statistics(c1, c2))" ] }, { "cell_type": "code", "execution_count": null, "id": "685fb176-eb85-4c9c-8348-54a7020a60d5", "metadata": {}, "outputs": [], "source": [ "# Plot histograms\n", "plt.figure(figsize=(8,5))\n", "\n", "# Plot three separate histograms\n", "fig, axs = plt.subplots(1, 2, figsize=(15, 4))\n", "\n", "axs[0].hist(p_values_analytical, bins=20, edgecolor='black')\n", "axs[0].set_title('Analytical t-test p-values')\n", "axs[0].set_xlabel('p-value')\n", "axs[0].set_ylabel('Frequency')\n", "axs[0].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[0].legend()\n", "\n", "axs[1].hist(p_values_order, bins=20, edgecolor='black')\n", "axs[1].set_title('Order-statistic p-values')\n", "axs[1].set_xlabel('p-value')\n", "axs[1].set_ylabel('Frequency')\n", "axs[1].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[1].legend()\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "7fdd3c8f-089e-4cf3-9b7a-50bfe5b021e3", "metadata": {}, "source": [ "Rats, it appears Sean was right and there is likely no systematic bias in pset grading. And week 4 was probably just a fluke random week your class did better. Now you owe HIM a bag of donuts." ] }, { "cell_type": "markdown", "id": "677eeab7-893e-4950-bf5e-74c80bb88084", "metadata": {}, "source": [ "# 3. The distribution of p-values under an effect\n", "\n", "## Simulation time!\n", "\n", "Hopefully you are now convinced that under the null p-values are uniformly distributed. But what about when there is a real effect?\n", "\n", "Below is the for generating synthetic p-set data. Lets try a few tests. Also, because we know the variance we input we can also calculate the p value with the official z score statistic.\n", "\n", "First lets add more data to generate the distribution under the null. Next try modulating cohort_diff that will start to acctually systematically shift the p-values. See what happens as the true difference increases!\n", "\n", "I also encourage you to try changing the ammount of data. That also influence p-value calculations in interesting ways!" ] }, { "cell_type": "code", "execution_count": null, "id": "94220a7b-2c87-4311-a77b-0c23eac9b495", "metadata": {}, "outputs": [], "source": [ "def generate_pset_data(\n", " n_students=30,\n", " n_psets=12,\n", " std_dev=0.5,\n", " cohort_diff=0,\n", " lower_bound=7,\n", " upper_bound=9.5\n", "):\n", " \"\"\"\n", " Generate synthetic problem set scores for two cohorts.\n", "\n", " Parameters\n", " ----------\n", " n_students : int Number of students per cohort.\n", " n_psets : int Number of problem sets.\n", " std_dev : float Standard deviation of scores for each problem set.\n", " cohort_diff : float Amount by which cohort 2 is lower than cohort 1.\n", " lower_bound : float Minimum possible mean for problem sets.\n", " upper_bound : float Maximum possible mean for problem sets.\n", "\n", " Returns\n", " -------\n", " scores_cohort1 : list of np.ndarray List of arrays, one per problem set, for cohort 1.\n", " scores_cohort2 : list of np.ndarray List of arrays, one per problem set, for cohort 2.\n", " mean_psets : list of floas The generated mean grades for each problem set.\n", " \"\"\"\n", "\n", "\n", " # Generate random mean for each problem set\n", " mean_psets = [round(random.uniform(lower_bound, upper_bound), 1) for _ in range(n_psets)]\n", "\n", " scores_cohort1 = []\n", " scores_cohort2 = []\n", "\n", " for mean_grade in mean_psets:\n", " # Generate scores for cohort 1\n", " cohort1 = np.round(\n", " np.random.normal(loc=mean_grade, scale=std_dev, size=n_students),\n", " decimals=1\n", " )\n", "\n", " # Generate scores for cohort 2 (shifted by cohort_diff)\n", " cohort2 = np.round(\n", " np.random.normal(loc=mean_grade - cohort_diff, scale=std_dev, size=n_students),\n", " decimals=1\n", " )\n", "\n", " scores_cohort1.append(cohort1)\n", " scores_cohort2.append(cohort2)\n", "\n", " return scores_cohort1, scores_cohort2, mean_psets" ] }, { "cell_type": "code", "execution_count": null, "id": "560c39a6-476b-46ee-a196-01be83705730", "metadata": {}, "outputs": [], "source": [ "def pvalue_known_variance(data1, data2, stdev):\n", " \"\"\"\n", " Compute a p-value for the difference of two means when the\n", " population variance is known (two-sided z-test).\n", " I don't feel like annotating this at the moment super good.\n", " \"\"\"\n", "\n", " # sample sizes\n", " n1 = len(data1)\n", " n2 = len(data2)\n", "\n", " # sample means\n", " mean1 = np.mean(data1)\n", " mean2 = np.mean(data2)\n", "\n", " # difference in means\n", " diff = mean1 - mean2\n", "\n", " # standard error of the difference\n", " true_variance = stdev**2\n", " se = np.sqrt(true_variance/n1 + true_variance/n2)\n", "\n", " # z statistic\n", " z = diff / se\n", "\n", " # two-sided p-value\n", " p_value = 2 * stats.norm.sf(abs(z))\n", "\n", " return p_value" ] }, { "cell_type": "code", "execution_count": null, "id": "80e825dc-9dd4-4074-b36b-226c447b2dcf", "metadata": {}, "outputs": [], "source": [ "# No actual difference\n", "scores_cohort1, scores_cohort2, _ = generate_pset_data(n_students = 30, n_psets=500)\n", "\n", "# Calculate p-values for difference in means for each p-set\n", "p_values_analytical = []\n", "p_values_order = []\n", "p_values_z = []\n", "\n", "for c1, c2 in zip(scores_cohort1, scores_cohort2):\n", " p_values_analytical.append(pvalue_analytical(c1, c2))\n", " p_values_order.append(pvalue_order_statistics(c1, c2))\n", " p_values_z.append(pvalue_known_variance(c1, c2, 0.5))" ] }, { "cell_type": "code", "execution_count": null, "id": "b705a177-1098-4bd0-adae-05a5ece435e2", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(8,5))\n", "\n", "# Plot three separate histograms for the different types of p-value derivations\n", "fig, axs = plt.subplots(1, 3, figsize=(15, 4))\n", "\n", "axs[0].hist(p_values_analytical, bins=20, edgecolor='black')\n", "axs[0].set_title('Analytical t-test p-values')\n", "axs[0].set_xlabel('p-value')\n", "axs[0].set_ylabel('Frequency')\n", "axs[0].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[0].legend()\n", "\n", "axs[1].hist(p_values_order, bins=20, edgecolor='black')\n", "axs[1].set_title('Order-statistic p-values')\n", "axs[1].set_xlabel('p-value')\n", "axs[1].set_ylabel('Frequency')\n", "axs[1].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[1].legend()\n", "\n", "axs[2].hist(p_values_z, bins=20, edgecolor='black')\n", "axs[2].set_title('Z-test (known variance) p-values')\n", "axs[2].set_xlabel('p-value')\n", "axs[2].set_ylabel('Frequency')\n", "axs[2].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[2].legend()\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "90aa8f14-6ff1-4be3-9331-342efc5c211c", "metadata": {}, "source": [ "Looks like all three of our methods for calculating are pretty similar, which is good! And indeed with more data it does appear that the distribution is quite uniform. You will note this takes a long time, that is because the permutation test takes a long time unfortunately (even though it is my favorite).\n", "\n", "Now try messing with the generation of the data, specifically the cohort_diff which will change the means between the two groups. Also try the n_students. This does take a long time so I suggest you do something to save the outputs so you can compare them. In the future I will code this better but for now this is the best I got." ] }, { "cell_type": "code", "execution_count": null, "id": "eea3a0cd-f8b1-4150-a066-506cf54b6f47", "metadata": {}, "outputs": [], "source": [ "# Meaningful difference\n", "scores_cohort1_dif, scores_cohort2_dif, _ = generate_pset_data(n_students = 30, n_psets=500, cohort_diff=0.1)\n", "\n", "# Calculate p-values for difference in means for each p-set\n", "p_values_analytical_dif = []\n", "p_values_order_dif = []\n", "p_values_z_dif = []\n", "\n", "for c1, c2 in zip(scores_cohort1_dif, scores_cohort2_dif):\n", " p_values_analytical_dif.append(pvalue_analytical(c1, c2))\n", " p_values_order_dif.append(pvalue_order_statistics(c1, c2))\n", " p_values_z_dif.append(pvalue_known_variance(c1, c2, 0.5))" ] }, { "cell_type": "code", "execution_count": null, "id": "d3048f6b-f41b-4682-baf9-4e82b0bdbb02", "metadata": {}, "outputs": [], "source": [ "# Plot histograms\n", "plt.figure(figsize=(8,5))\n", "\n", "# Plot three separate histograms\n", "fig, axs = plt.subplots(1, 3, figsize=(15, 4))\n", "\n", "axs[0].hist(p_values_analytical_dif, bins=20, edgecolor='black')\n", "axs[0].set_title('Analytical t-test p-values')\n", "axs[0].set_xlabel('p-value')\n", "axs[0].set_ylabel('Frequency')\n", "axs[0].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[0].set_xlim(0, 1)\n", "axs[0].legend()\n", "\n", "axs[1].hist(p_values_order_dif, bins=20, edgecolor='black')\n", "axs[1].set_title('Order-statistic p-values')\n", "axs[1].set_xlabel('p-value')\n", "axs[1].set_ylabel('Frequency')\n", "axs[1].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[1].set_xlim(0, 1)\n", "axs[1].legend()\n", "\n", "axs[2].hist(p_values_z_dif, bins=20, edgecolor='black')\n", "axs[2].set_title('Z-test (known variance) p-values')\n", "axs[2].set_xlabel('p-value')\n", "axs[2].set_ylabel('Frequency')\n", "axs[2].axvline(0.05, color='red', linestyle='--', label='0.05 threshold')\n", "axs[2].set_xlim(0, 1)\n", "axs[2].legend()\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "53f1e980-4591-4116-a48a-6d3dab70322a", "metadata": {}, "source": [ "As you increase the actualy difference you start to shift the p-value curve towards 0 (annoying called right shifted). However it is still possible to get high p-values under a meaningful difference, just like you can get low p-values with a non-meaningful difference." ] }, { "cell_type": "markdown", "id": "6f40e530-2f29-4f82-925a-c9f96e41771e", "metadata": {}, "source": [ "# 4. P-hacking \n", "# Lets talk about data manipulation and how to catch it\n", "\n", "There has been a very successfull psy-opp to convince everyone that p < 0.05 is necessary to publish. You may have heard Goodhardt's law \"Any observed statistical regularity will tend to collapse once pressure is palced upon it for control purposes\". As a result, people will artifically try to get specific p-values (whether concious or not). Can you think of ways to do this? Maybe try them out above!\n", "\n", "But we can use the same statistics to catch these potential cheaters, or at least identify it is happening.\n", "\n", "Under a real difference from the null p-values quickly get very small, even with sparse data. We can use this knowledge to detect so called \"p-hacking\", researchers manipulating their data to get a \"significant\" result. I am drawing data here from a very good paper that explains this phenomena The Extent and Consequences of P-Hacking in Science (https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.1002106). In this paper the authors collected a bunch of reported p-values and looking to see if it deviates from expectation in a very clever way. Now since people only usually publish meaningful p-values, they are really looking at a conditional probablity i.e. P(p | p < 0.05). But the ideas behind the p-curve still apply." ] }, { "cell_type": "code", "execution_count": null, "id": "61468cb5-2433-46b3-9e00-d7a0db328c95", "metadata": {}, "outputs": [], "source": [ "path = 'p_values_meta.csv'\n", "df = pd.read_csv(path)" ] }, { "cell_type": "code", "execution_count": null, "id": "f27f39e7-f300-4282-bf0d-e768757b91fa", "metadata": {}, "outputs": [], "source": [ "plt.hist( df[df['P-VALUES FOR P-CURVES'] > 0 ]['P-VALUES FOR P-CURVES'] )" ] }, { "cell_type": "markdown", "id": "9498b3c0-2497-440a-8ed9-0b4e7118cfe6", "metadata": {}, "source": [ "You can see there is a very clear bump right below 0.05. The p-curve definintely appears to be exponentially decreasing, so a significant increase as p increases should sound alarms.\n", "\n", "In the paper they do some more fancy statistics than what we are doing. However as a more simple version we can make the assumption that because hte p-curve decays so exponentially, the number of p-values in the bins 0.04 0.045) & (p_values <= 0.05)\n", "bin_lower = (p_values > 0.04) & (p_values <= 0.045)\n", "\n", "count_upper = bin_upper.sum()\n", "count_lower = bin_lower.sum()\n", "\n", "print(f\"Count 0.045–0.05: {count_upper}\")\n", "print(f\"Count 0.04–0.045: {count_lower}\")\n", "\n", "# Total in both bins\n", "n_total = count_upper + count_lower\n", "\n", "# Binomial test: probability of being in upper bin if chance = 0.5\n", "res = binomtest(count_upper, n=n_total, p=0.5, alternative='greater')\n", "\n", "# Proportion in upper bin\n", "prop_upper = count_upper / n_total if n_total > 0 else np.nan\n", "\n", "print(f\"Proportion in 0.045–0.05: {prop_upper:.3f}\")\n", "print(f\"Binomial test p-value: {res.pvalue:.4f}\")" ] }, { "cell_type": "markdown", "id": "ae6bffc4-9448-41ea-926e-b770d3998df3", "metadata": {}, "source": [ "So it is quite surprising that we get this many p-values close to 0.05, suggestive that people are manipulating their results to get what they want! But be careful not to misinterpret the p-value!" ] }, { "cell_type": "code", "execution_count": null, "id": "a2bda3be-1e03-4bc2-b7a8-a16712ae5f11", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "seq_analysis", "language": "python", "name": "seq_analysis" }, "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.12.8" } }, "nbformat": 4, "nbformat_minor": 5 }