effect_size.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Effect Size\n",
    "===\n",
    "\n",
    "Examples and exercises for a tutorial on statistical inference.\n",
    "\n",
    "Copyright 2016 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "import numpy\n",
    "import scipy.stats\n",
    "\n",
    "import matplotlib.pyplot as pyplot\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# seed the random number generator so we all get the same results\n",
    "numpy.random.seed(17)\n",
    "\n",
    "# some nice colors from http://colorbrewer2.org/\n",
    "COLOR1 = '#7fc97f'\n",
    "COLOR2 = '#beaed4'\n",
    "COLOR3 = '#fdc086'\n",
    "COLOR4 = '#ffff99'\n",
    "COLOR5 = '#386cb0'\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part One\n",
    "\n",
    "To explore statistics that quantify effect size, we'll look at the difference in height between men and women.  I used data from the Behavioral Risk Factor Surveillance System (BRFSS) to estimate the mean and standard deviation of height in cm for adult women and men in the U.S.\n",
    "\n",
    "I'll use `scipy.stats.norm` to represent the distributions.  The result is an `rv` object (which stands for random variable)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu1, sig1 = 178, 7.7\n",
    "male_height = scipy.stats.norm(mu1, sig1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu2, sig2 = 163, 7.3\n",
    "female_height = scipy.stats.norm(mu2, sig2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function evaluates the normal (Gaussian) probability density function (PDF) within 4 standard deviations of the mean.  It takes and rv object and returns a pair of NumPy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eval_pdf(rv, num=4):\n",
    "    mean, std = rv.mean(), rv.std()\n",
    "    xs = numpy.linspace(mean - num*std, mean + num*std, 100)\n",
    "    ys = rv.pdf(xs)\n",
    "    return xs, ys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the two distributions look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs, ys = eval_pdf(male_height)\n",
    "pyplot.plot(xs, ys, label='male', linewidth=4, color=COLOR2)\n",
    "\n",
    "xs, ys = eval_pdf(female_height)\n",
    "pyplot.plot(xs, ys, label='female', linewidth=4, color=COLOR3)\n",
    "pyplot.xlabel('height (cm)')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assume for now that those are the true distributions for the population.\n",
    "\n",
    "I'll use `rvs` to generate random samples from the population distributions.  Note that these are totally random, totally representative samples, with no measurement error!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "male_sample = male_height.rvs(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "female_sample = female_height.rvs(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both samples are NumPy arrays.  Now we can compute sample statistics like the mean and standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean1, std1 = male_sample.mean(), male_sample.std()\n",
    "mean1, std1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sample mean is close to the population mean, but not exact, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean2, std2 = female_sample.mean(), female_sample.std()\n",
    "mean2, std2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the results are similar for the female sample.\n",
    "\n",
    "Now, there are many ways to describe the magnitude of the difference between these distributions.  An obvious one is the difference in the means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "difference_in_means = male_sample.mean() - female_sample.mean()\n",
    "difference_in_means # in cm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On average, men are 14--15 centimeters taller.  For some applications, that would be a good way to describe the difference, but there are a few problems:\n",
    "\n",
    "* Without knowing more about the distributions (like the standard deviations) it's hard to interpret whether a difference like 15 cm is a lot or not.\n",
    "\n",
    "* The magnitude of the difference depends on the units of measure, making it hard to compare across different studies.\n",
    "\n",
    "There are a number of ways to quantify the difference between distributions.  A simple option is to express the difference as a percentage of the mean.\n",
    "\n",
    "**Exercise 1**: what is the relative difference in means, expressed as a percentage?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**STOP HERE**: We'll regroup and discuss before you move on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part Two\n",
    "\n",
    "An alternative way to express the difference between distributions is to see how much they overlap.  To define overlap, we choose a threshold between the two means.  The simple threshold is the midpoint between the means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_thresh = (mean1 + mean2) / 2\n",
    "simple_thresh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A better, but slightly more complicated threshold is the place where the PDFs cross."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh = (std1 * mean2 + std2 * mean1) / (std1 + std2)\n",
    "thresh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, there's not much difference between the two thresholds.\n",
    "\n",
    "Now we can count how many men are below the threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "male_below_thresh = sum(male_sample < thresh)\n",
    "male_below_thresh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And how many women are above it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "female_above_thresh = sum(female_sample > thresh)\n",
    "female_above_thresh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"overlap\" is the area under the curves that ends up on the wrong side of the threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "male_overlap = male_below_thresh / len(male_sample)\n",
    "female_overlap = female_above_thresh / len(female_sample)\n",
    "male_overlap, female_overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practical terms, you might report the fraction of people who would be misclassified if you tried to use height to guess sex, which is the average of the male and female overlap rates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "misclassification_rate = (male_overlap + female_overlap) / 2\n",
    "misclassification_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to quantify the difference between distributions is what's called \"probability of superiority\", which is a problematic term, but in this context it's the probability that a randomly-chosen man is taller than a randomly-chosen woman.\n",
    "\n",
    "**Exercise 2**:  Suppose I choose a man and a woman at random.  What is the probability that the man is taller?\n",
    "\n",
    "HINT: You can `zip` the two samples together and count the number of pairs where the male is taller, or use NumPy array operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overlap (or misclassification rate) and \"probability of superiority\" have two good properties:\n",
    "\n",
    "* As probabilities, they don't depend on units of measure, so they are comparable between studies.\n",
    "\n",
    "* They are expressed in operational terms, so a reader has a sense of what practical effect the difference makes.\n",
    "\n",
    "### Cohen's effect size\n",
    "\n",
    "There is one other common way to express the difference between distributions.  Cohen's $d$ is the difference in means, standardized by dividing by the standard deviation.  Here's the math notation:\n",
    "\n",
    "$ d = \\frac{\\bar{x}_1 - \\bar{x}_2} s $\n",
    "\n",
    "where $s$ is the pooled standard deviation:\n",
    "\n",
    "$s = \\sqrt{\\frac{n_1 s^2_1 + n_2 s^2_2}{n_1+n_2}}$\n",
    "\n",
    "Here's a function that computes it:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CohenEffectSize(group1, group2):\n",
    "    \"\"\"Compute Cohen's d.\n",
    "\n",
    "    group1: Series or NumPy array\n",
    "    group2: Series or NumPy array\n",
    "\n",
    "    returns: float\n",
    "    \"\"\"\n",
    "    diff = group1.mean() - group2.mean()\n",
    "\n",
    "    n1, n2 = len(group1), len(group2)\n",
    "    var1 = group1.var()\n",
    "    var2 = group2.var()\n",
    "\n",
    "    pooled_var = (n1 * var1 + n2 * var2) / (n1 + n2)\n",
    "    d = diff / numpy.sqrt(pooled_var)\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the denominator is a little complicated; in fact, people have proposed several ways to do it.  This implementation uses the \"pooled standard deviation\", which is a weighted average of the standard deviations of the two groups.\n",
    "\n",
    "And here's the result for the difference in height between men and women."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CohenEffectSize(male_sample, female_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most people don't have a good sense of how big $d=1.9$ is, so let's make a visualization to get calibrated.\n",
    "\n",
    "Here's a function that encapsulates the code we already saw for computing overlap and probability of superiority."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def overlap_superiority(control, treatment, n=1000):\n",
    "    \"\"\"Estimates overlap and superiority based on a sample.\n",
    "    \n",
    "    control: scipy.stats rv object\n",
    "    treatment: scipy.stats rv object\n",
    "    n: sample size\n",
    "    \"\"\"\n",
    "    control_sample = control.rvs(n)\n",
    "    treatment_sample = treatment.rvs(n)\n",
    "    thresh = (control.mean() + treatment.mean()) / 2\n",
    "    \n",
    "    control_above = sum(control_sample > thresh)\n",
    "    treatment_below = sum(treatment_sample < thresh)\n",
    "    overlap = (control_above + treatment_below) / n\n",
    "    \n",
    "    superiority = (treatment_sample > control_sample).mean()\n",
    "    return overlap, superiority"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the function that takes Cohen's $d$, plots normal distributions with the given effect size, and prints their overlap and superiority."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_pdfs(cohen_d=2):\n",
    "    \"\"\"Plot PDFs for distributions that differ by some number of stds.\n",
    "    \n",
    "    cohen_d: number of standard deviations between the means\n",
    "    \"\"\"\n",
    "    control = scipy.stats.norm(0, 1)\n",
    "    treatment = scipy.stats.norm(cohen_d, 1)\n",
    "    xs, ys = eval_pdf(control)\n",
    "    pyplot.fill_between(xs, ys, label='control', color=COLOR3, alpha=0.7)\n",
    "\n",
    "    xs, ys = eval_pdf(treatment)\n",
    "    pyplot.fill_between(xs, ys, label='treatment', color=COLOR2, alpha=0.7)\n",
    "    \n",
    "    o, s = overlap_superiority(control, treatment)\n",
    "    pyplot.text(0, 0.05, 'overlap ' + str(o))\n",
    "    pyplot.text(0, 0.15, 'superiority ' + str(s))\n",
    "    pyplot.show()\n",
    "    #print('overlap', o)\n",
    "    #print('superiority', s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example that demonstrates the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_pdfs(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And an interactive widget you can use to visualize what different values of $d$ mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slider = widgets.FloatSlider(min=0, max=4, value=2)\n",
    "interact(plot_pdfs, cohen_d=slider)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cohen's $d$ has a few nice properties:\n",
    "\n",
    "* Because mean and standard deviation have the same units, their ratio is dimensionless, so we can compare $d$ across different studies.\n",
    "\n",
    "* In fields that commonly use $d$, people are calibrated to know what values should be considered big, surprising, or important.\n",
    "\n",
    "* Given $d$ (and the assumption that the distributions are normal), you can compute overlap, superiority, and related statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, the best way to report effect size depends on the audience and your goals.  There is often a tradeoff between summary statistics that have good technical properties and statistics that are meaningful to a general audience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}