Statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Statistics for Biologists\n",
    "This notebook will discuss **descriptive** and **inferential** statistics, and introduce ways to implement them in Python.\n",
    "\n",
    "You should run the cells in this notebook step-by-step, completing necessary tasks and answering questions on your quiz along the way.\n",
    "\n",
    "### By the end of this notebook, you will be able to:\n",
    "* Identify when to use descriptive or inferential statistics\n",
    "* Apply the appropriate statistical tests to compare two groups\n",
    "* Use the stats package from SciPy to run simple tests in Python\n",
    "* Test differences between spiny and aspiny cells\n",
    "\n",
    "### Table of Contents\n",
    "1. [Population vs sample distributions](#one)\n",
    "2. [The Central Limit Theorem](#two)\n",
    "3. [Skewed Distributions](#three)\n",
    "4. [Hypothesis Testing](#four)\n",
    "5. [Testing the Difference Between Cell Types](#five)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "First things first, let's configure our coding environment with the packages we'll need. We'll use `import` to bring in additional code packages and nickname them. For example, we're importing a code package called `numpy` and nicknaming it `np`. In the future, when we refer to numpy, we'll write `np`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import our necessary toolboxes and tell matplotlib to plot inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "import pandas as pd\n",
    "from scipy import stats # new package for today, for statistical tools!\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina' # configure high-resolution plotting\n",
    "\n",
    "print('Packages imported.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<a id=\"one\"></a>\n",
    "## Part One: Population vs sample distributions\n",
    "**Descriptive statistics** summarize the main features of a data set.\n",
    "\n",
    "It's important to distinguish between the following:\n",
    "* **Observation**: result from one trial of an experiment\n",
    "* **Sample**: results from multiple independent trials\n",
    "* **Population**: the *ground truth*; all possible observations that could be seen\n",
    "\n",
    "Distributions differ in their **location** (mean, $\\mu$) and **spread** (standard deviation, $\\sigma$). Using `np.random.normal`, which draws a given number of data points from a defined population, we'll define a **population distribution** and plot it. [`np.random.normal`](https://numpy.org/doc/stable/reference/random/generated/numpy.random.normal.html) has three arguments: a mu, a sigma, and a size.\n",
    "\n",
    ">**Task**: Run the cell below. Change the mu and sigma to see how this changes the plotted distribution. When you're done observing how this changes your distribution, change the mu back to 3 and the sigma back to 2 and re-run the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Decide on a mean and a standard deviation\n",
    "mu = 3\n",
    "sigma = 2\n",
    "\n",
    "# Use np.random.normal to create a distribution of 10,000 points with our given mu & sigma\n",
    "pop = np.random.normal(mu, sigma, 10000)\n",
    "\n",
    "# Plot a histogram with 30 bins\n",
    "# Giving it the argument density=True will plot normalized counts\n",
    "# This will create a probability density (rather than raw counts)\n",
    "plt.hist(pop, 30, density=True)\n",
    "plt.axvline(mu,color='r')\n",
    "plt.title('Population distribution of 10,000 points')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are various ways we can describe the distribution of the dataset, beyond the standard deviation:\n",
    "* Range (minimum and maximum)\n",
    "* Variance ($\\sigma^2$)\n",
    "* Standard Error of the Mean (S.E.M., $\\sigma/\\sqrt{n}$)\n",
    "* Confidence Intervals\n",
    "\n",
    "We can easily get many of these descriptive statistics by using the `scipy stats` package method `describe()`. [Documentation here.](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.describe.html). Note that <code>stats.describe()</code> doesn't give us the median (annoyingly) but `np.median()` can!\n",
    "* The `stats.skewtest()` method ([documentation here]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "stats.describe(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We created a normal distribution from a mean of 3 but with limited data points, so these values are *just* slightly off from what we would expect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our variable `pop` is the \"ground truth\" population. However, we'll rarely have *10,000* datapoints in our sample. So, let's generate a more realistic sample, and see how the mean compares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Create a sample distribution with less data points\n",
    "sample_mean, sample_sigma = 3, 2\n",
    "sample = np.random.normal(sample_mean, sample_sigma, 20)\n",
    "\n",
    "# Look at the descriptive statistics of our sample\n",
    "print(stats.describe(sample))\n",
    "\n",
    "# Plot our histogram, with alpha to 0.5 which will make the chart slightly transparent\n",
    "plt.hist(pop, 30, alpha=0.5, density=True)\n",
    "plt.hist(sample, 30, alpha=0.5, color='r',density=True)\n",
    "plt.axvline(np.mean(pop),color='blue') # Take the mean and plot a vertical blue line \n",
    "plt.axvline(np.mean(sample),color='red') # Take the mean and plot a vertical red line \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">Answer Q1 on the quiz.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<a id=\"two\"></a>\n",
    "## Part Two: The Central Limit Theorem\n",
    "\n",
    "With fewer samples, the mean of the sample distribution tends to be further from the mean of the population distribution. This is known as the **central limit theorem**, which states that the distribution of sample means will become increasingly close to a normal distribution as the sample size increases, regardless of the shape of the population distribution.\n",
    "\n",
    ">**Task**: The cell below will illustrate the central limit theorem by plotting the distribution of sample means with different sample sizes. All you need to do is run the cell! This cell will take ~10 seconds to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,5,figsize=(20,5),sharey=True) # Set up the plot\n",
    "\n",
    "mu = 0 # These distributions will have a mean of zero\n",
    "sample_sizes = [5,10,25,100,250] # Set a list of sample sizes for our plots \n",
    "sample_means = []\n",
    "\n",
    "# For each subplot, create a plot.\n",
    "for a in range(len(ax)):   \n",
    "    \n",
    "    # Calculate the mean of sample of sample_size designated above, 1000 times\n",
    "    for x in range(1000):\n",
    "        sample_size = sample_sizes[a]\n",
    "        sample_dist = np.random.normal(mu, 10, sample_size) # Create a normal distribution with mu, sigma\n",
    "        sample_means.append(np.mean(sample_dist)) # Append the mean of this distribution\n",
    "        \n",
    "    ax[a].hist(sample_means,color='teal',alpha = .5) # Plot the distribution of means\n",
    "    ax[a].set_title('size='+str(sample_size)+',mean='+str(np.round(np.mean(sample_means),3))+',stdev='+str(np.round(np.std(sample_means),2)))\n",
    "    ax[a].set_xlim([-20,20])\n",
    "    sample_means = [] # Reset the sample means\n",
    "\n",
    "plt.suptitle('Distributions of 10,000 sample means',fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<a id=\"three\"></a>\n",
    "## Part Three: Skewed Distributions\n",
    "\n",
    "However, not every population in nature is **normally distributed**. In fact, most populations are slightly skewed. Let's demonstrate a population distribution and sample distribution that is drawn from a [gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution).\n",
    "\n",
    "We're going to use [`np.random.gamma`](https://numpy.org/doc/stable/reference/random/generated/numpy.random.gamma.html) to create this distribution. This function takes several parameters: the first two define the shape and scale, the third defines the size of the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a skewed distribution of 10,000 points with our given mu & sigma\n",
    "pop_size = 10000\n",
    "sample_size = 30\n",
    "\n",
    "skewed_pop = np.random.gamma(7.5,1,pop_size)\n",
    "skewed_sample = np.random.gamma(7.5,1,sample_size)\n",
    "\n",
    "pop_stats = stats.describe(skewed_sample)\n",
    "sample_stats = stats.describe(skewed_pop)\n",
    "\n",
    "plt.hist(skewed_pop, 30, alpha = .3, density=True)\n",
    "plt.hist(skewed_sample, 30, alpha = .3, density=True)\n",
    "plt.axvline(pop_stats.mean,color='blue')\n",
    "plt.axvline(sample_stats.mean,color='orange')\n",
    "plt.legend(['Population','Sample'])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You might notice that with this skewed population, the mean is a pretty poor descriptor of both distributions. **When the skew is bad (*statistically bad*), we should report the median.**\n",
    "\n",
    "><b>Task</b>: Rework our demonstration of the central limit theorem for a skewed, rather than a normal, population. You can do so by changing the type of histogram that the central limit theorem demonstration is drawing from. In other words, change the line that says `np.random.normal` to create a gamma distribution with the parameters above. **Be sure to keep the original parameters (7.5, 1) for the shape of the gamma distribution, but use the `sample_size` variable to define the size.**\n",
    "\n",
    "> Does the theorem still hold? Answer Q2 on the quiz."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the skewness of the Allen Institute time constant data\n",
    "\n",
    "As an example of how you'd implement these statistics, let's look at some data that characterizes the intrinsic electrophysiology in different cell types, and includes information about whether the cells have spiny or aspiny dendrites. Here, we'll test for significant differences in those data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('Data/mouse_cell_metrics.csv')\n",
    "spiny_data = data[data['dendrite_type']=='spiny']\n",
    "aspiny_data = data[data['dendrite_type']=='aspiny']\n",
    "\n",
    "spiny_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to know if our data is skewed or not. To test for skewness, we can use [`stats.skewtest`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skewtest.html#scipy.stats.skewtest) to test. This method implements the <a href=\"https://www.jstor.org/stable/2684359?seq=1\">D'Agostino-Pearson skewness test</a>, one of many different tests (e.g., the Kolmogorov-Smirov test) that can be used to check the normality of a distribution. **If the skew test gives us a p-value of less than 0.05, the population is skewed.**\n",
    "\n",
    ">**Task**: Use the cell below to check the skew of both the aspiny and spiny **tau** data. You'll need to change the sample that is created in the first line in order to change spiny to aspiny. Once you're done, answer Q3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sample = spiny_data['tau'] # Choose which data to use\n",
    "\n",
    "stat,pvalue = stats.skewtest(sample) # Run the skew test\n",
    "\n",
    "# Print the p value of the skew test up to 30 decimal points\n",
    "print('The skewtest p-value is ' + '%.30f' % pvalue) \n",
    "\n",
    "plt.hist(sample) # Create a histogram\n",
    "plt.ylabel('# of Cells')\n",
    "plt.xlabel('Rheobase')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<a id=\"four\"></a>\n",
    "## Part Four: Hypothesis Testing\n",
    "\n",
    "### *Inferential statistics* generalize from observed data to the world at large\n",
    "\n",
    "Most often, the goal of our hypothesis testing is to test whether or not two distributions are different, or if a distribution has a different mean than the underlying population distribution.\n",
    "\n",
    "The SciPy stats package has [many hypothesis testing tools](https://docs.scipy.org/doc/scipy/reference/stats.html) (see Statistical Tests). In most cases, we will be testing whether two or more distributions are different from eachother.\n",
    "\n",
    "If we know our distributions are normal (they're generated from a normal distribution!) we could use **parametric statistics** to test our hypothesis. To test for differences between normal populations, we can use the independent t-test in our stats package: `stats.ttest_ind()`. If we had paired samples, we would use a dependent t-test [as seen here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_rel.html#scipy.stats.ttest_rel).\n",
    "\n",
    "If one of our populations is skewed, however, we **cannot use a t-test**. A t-test assumes that the populations are normally distributed. For skewed populations, we can use either the [Mann-Whitney U](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html#scipy.stats.mannwhitneyu) (for independent samples, `stats.mannwhitneyu()`) or the [Wilcoxon Signed Rank Test](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wilcoxon.html#scipy.stats.wilcoxon) (for dependent/paired samples,`stats.wilcoxon()`). Below, we'll generate two random populations, one normal and one skewed, and see whether they are significantly different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "skewed_pop = np.random.gamma(8,1,100) # Create a skewed distribution with 100 points\n",
    "comparison_pop = np.random.normal(7.5,2,100) # Create a normal distribution with 100 points\n",
    "\n",
    "print(stats.ttest_ind(skewed_pop,comparison_pop)) # to run an independent t-test\n",
    "print(stats.mannwhitneyu(skewed_pop,comparison_pop)) # to run a mannwhitneyu\n",
    "print(stats.wilcoxon(skewed_pop,comparison_pop)) # to run a wilcoxon signed rank test\n",
    "\n",
    "# Plot the distributions we created above\n",
    "plt.hist(skewed_pop,alpha=.5)\n",
    "plt.hist(comparison_pop,alpha=.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"five\"></a>\n",
    "## Part 5. Do spiny and aspiny cells have significantly different time constants?\n",
    "\n",
    "> *Task*: Using the code provided above, and given what we know about the skewness of our data, choose and run a statistical test that compares `aspiny_data['tau']` to `spiny_data['tau']`. When you're done, answer Q4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run your code here!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<a id=\"refs\"></a>\n",
    "## References & resources\n",
    "I *strongly* recommend reading the [Points of significance](https://www.nature.com/collections/qghhqm/pointsofsignificance) series from *Nature* which covers many of these topics. This lecture specifically focuses on [The Importance of Being Uncertain](https://www.nature.com/articles/nmeth.2613).\n",
    "\n",
    "For a demonstration on how to animate the central limit theorem, see [this Github notebook](https://github.com/rohanjoseph93/Central-Limit-Theorem/blob/master/Central%20Limit%20Theorem.ipynb).\n",
    "\n",
    "Consider working through the examples in [Inferential thinking](https://www.inferentialthinking.com/chapters/11/Testing_Hypotheses.html).\n",
    "\n",
    "This notebook borrows code from [Hypothesis tests in Python](https://datasciencechalktalk.com/2019/09/02/hypothesis-tests-with-python/) by Valentina Alto and ideas from a variety of other sources, including [Towards Data Science](https://towardsdatascience.com/hypothesis-testing-in-machine-learning-using-python-a0dc89e169ce)."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}