Basic statistics in data analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CCP-BioSim Training Course:\n",
    "### Basic statistical methods in MD data analysis\n",
    "\n",
    "When you are analysing MD trajectories frequently you are dealing with lots of data, and the data is noisy. Showing that measurements you make and conclusions you come to are statistically significant (or not!) is important.\n",
    "\n",
    "Here you will look at the first 100 microseconds of the 1.03 millisecond ANTON BPTI dataset from Shaw et al. (http://science.sciencemag.org/content/330/6002/341.full). The trajectory file contains just the coordinates of the C-alpha atoms, and is sampled every ten nanoseconds.\n",
    "\n",
    "We begin by importing the key Python libraries: numpy, mdtraj, and matplotlib:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import mdtraj as mdt\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline \n",
    "plt.rcParams.update({'font.size': 15}) #This sets a better default label size for plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we load the trajectory data, and create a numpy array with the timepoint data which will be used later when we plot things:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajfile = 'data/statistical_analysis/bpti_ca_100us_dt10ns.xtc'\n",
    "topfile = 'data/statistical_analysis/bpti_ca.pdb'\n",
    "t_100 = mdt.load(trajfile, top=topfile)\n",
    "time_100 = np.arange(len(t_100)) * 0.01 # The time, in microseconds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's investigate a metric from this simulation: the average distance between the first and last C-alpha atoms.\n",
    "\n",
    "To begin with let's imagine we ran the simulation for just ten microseconds so only have the first 1000 data points. In the next cell we cut the trajectory down to size, and then use MDTraj's `compute_distances()` method to find the end-to-end distance in each frame: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_10 = t_100[:1000] # A new trajectory that just contains the first 1000 snapshots\n",
    "time_10 = time_100[:1000] # An array of time data to match t_10\n",
    "ends = np.array([[0, t_10.topology.n_atoms - 1]])\n",
    "dist = mdt.compute_distances(t_10, ends)[:,0]\n",
    "plt.plot(time_10, dist)\n",
    "plt.xlabel('time (micro_sec)')\n",
    "plt.ylabel('end-to-end distance (nm)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what's the mean end-to-end distance? Here's a very simple approach - calculate the mean, and because we know a little bit about statistics, we also calculate the standard error of the mean, and from this the classic 95% confidence interval. For this we make use of numpy's `mean()` and `std()` (standard deviation) methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following four lines create lists we will use to store values we calculate as we work through this notebook.\n",
    "d_means = [None] * 4\n",
    "d_c95lows = [None] * 4\n",
    "d_c95highs = [None] * 4\n",
    "methods = [None] * 4\n",
    "\n",
    "d_means[0] = dist.mean()\n",
    "stderr = dist.std()/np.sqrt(len(dist))\n",
    "d_c95lows[0] = d_means[0] - 1.96 * stderr # The lower bound of the 95% confidence interval\n",
    "d_c95highs[0] = d_means[0] + 1.96 * stderr # The upper bound of the 95% confidence interval\n",
    "\n",
    "text = \"{}:\\nMean end to end distance: {:3.3f} nm, 95% confidence interval {:3.3f} nm to {:3.3f} nm\\n\"\n",
    "methods[0] = \"Naive approach\"\n",
    "for i in range(1):\n",
    "    print(text.format(methods[i], d_means[i], d_c95lows[i], d_c95highs[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above uses every value for the end-to-end distance. But is this valid? Only if the individual values are uncorrelated (see https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation#Effect_of_autocorrelation_(serial_correlation)). Let's check for this.\n",
    "\n",
    "The cell below calculates the autocorrelation function for the distances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d2 = dist - dist.mean()\n",
    "c = np.correlate(d2, d2, 'full')\n",
    "c = c[int(c.size/2):]\n",
    "c /= (c.var() * np.arange(d2.size, 0, -1))\n",
    "plt.xlabel('sample interval')\n",
    "plt.ylabel('correlation')\n",
    "plt.plot(c[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The autocorellation function only reaches zero after a time lag of about 4 samples (equivalent to 40 nanoseconds), so we should subsample the distance data at this interval to get truely uncorrelated values from which to calculate statistics.\n",
    "\n",
    "This is what we do in the next cell, then recalculate the mean end-to-end distance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_uncorr = dist[::4]\n",
    "time_uncorr = time_10[::4]\n",
    "d_means[1] = d_uncorr.mean()\n",
    "stderr = d_uncorr.std()/np.sqrt(len(d_uncorr))\n",
    "d_c95lows[1] = d_means[1] - stderr * 1.96\n",
    "d_c95highs[1] = d_means[1] + stderr * 1.96\n",
    "methods[1] = \"After checking for correlation\"\n",
    "for i in range(2):\n",
    "    print(text.format(methods[i], d_means[i], d_c95lows[i], d_c95highs[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the new estimate of the mean end-to-end distance is considerably lower - below the 95% confidence limit in the original analysis.\n",
    "\n",
    "Up to now we have implicitly been working under the assumtion that the data come from a normal distribution - but is this likely? Not really, after all, the aim of MD is to generate a Boltzmann weighted ansemble, not a normal one. If our data is not normally distributed then exactly what a \"mean value\" means, and what the confidence limits on it are, will need to be reconsidered.\n",
    "\n",
    "In the cell below we plot a histogram of the distance values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H, edges = np.histogram(d_uncorr, bins='auto', density=True)\n",
    "midpoints = [edges[i:i+2].mean() for i in range(len(edges) - 1)]\n",
    "plt.xlabel('end-to-end distance')\n",
    "plt.ylabel('probability density')\n",
    "plt.plot(midpoints, H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yup -= not a normal distribution at all. In this situation it is better to use the method of bootstrapping https://en.wikipedia.org/wiki/Bootstrapping_(statistics) to estimate the confidence interval for the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_means = []\n",
    "for i in range(1000):\n",
    "    b_means.append(np.random.choice(d_uncorr, size=d_uncorr.size).mean())\n",
    "b_means.sort()\n",
    "d_means[2] = d_means[1] # our estimate of the mean is unchanged from last time, only the confidence interval\n",
    "d_c95lows[2] = b_means[50] # because only 50/1000 (5%) of the estimated means are less than this\n",
    "d_c95highs[2] = b_means[950] # because only 50/1000 (5%) of the estimated means are more than this\n",
    "methods[2] = \"With bootstrapping\"\n",
    "for i in range(3):\n",
    "    print(text.format(methods[i], d_means[i], d_c95lows[i], d_c95highs[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that by taking proper account of the non-normal distribution in the end-to-end distance samples, the 95% confidence range tightens up quite a bit.\n",
    "\n",
    "OK, we have a mean and some confidence intervals. But is this a converged value? Is the simulation long enough? \n",
    "\n",
    "One simple approach is to say that it may be, if the mean end to end distance calculated from the first half of the simulation is close to the value calculated from the second half. Below we plot the difference between the means for the first and second halves of the simulation as a function of increasing total simulation length. The expectation might be that as this increases, the difference between the means will tend to zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_length_difference = []\n",
    "stepsize = int(len(d_uncorr)/100)\n",
    "for i in range(stepsize, len(d_uncorr) + stepsize, stepsize):\n",
    "    j = int(i/2)\n",
    "    mean_length_difference.append(d_uncorr[:j].mean() - d_uncorr[j:i].mean())\n",
    "plt.xlabel('simulation length (microsec.)')\n",
    "plt.ylabel('difference in block mesns')\n",
    "plt.plot(time_uncorr[1::stepsize], mean_length_difference)\n",
    "plt.plot(time_uncorr[1::stepsize], np.zeros(len(mean_length_difference))) # add a line at y=0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm - worryingly, the difference is tending to a positive number, suggesting that the sampling is not converged.\n",
    "\n",
    "Let's go back to the full dataset (100 microseconds), and repeat the process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist = mdt.compute_distances(t_100, ends)[:,0]\n",
    "d_uncorr = dist[::4]\n",
    "time_uncorr = time_100[::4]\n",
    "\n",
    "mean_length_difference = []\n",
    "stepsize = int(len(d_uncorr)/100)\n",
    "for i in range(stepsize, len(d_uncorr) + stepsize, stepsize):\n",
    "    j = int(i/2)\n",
    "    mean_length_difference.append(d_uncorr[:j].mean() - d_uncorr[j:i].mean())\n",
    "plt.xlabel('simulation length (microsec.)')\n",
    "plt.ylabel('difference in block mesns')\n",
    "plt.plot(time_uncorr[1::stepsize], mean_length_difference)\n",
    "plt.plot(time_uncorr[1::stepsize], np.zeros(len(mean_length_difference)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's better. Now recalculate the mean end-to-end distance using our best method, and compare with the previous estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_means[3] = d_uncorr.mean()\n",
    "b_means = []\n",
    "for i in range(1000):\n",
    "    b_means.append(np.random.choice(d_uncorr, size=d_uncorr.size).mean())\n",
    "b_means.sort()\n",
    "d_c95lows[3] = b_means[50]\n",
    "d_c95highs[3] = b_means[950]\n",
    "methods[3] = \"With ten times more data\"\n",
    "for i in range(4):\n",
    "    print(text.format(methods[i], d_means[i], d_c95lows[i], d_c95highs[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By chance, the \"best\" estimate is very close to the \"naive\" one - but that's just a stroke of luck!\n",
    "\n",
    "### Summary\n",
    "MD generates lots of data, so there is every reason to take full advantage of it and put your data analysis on a sound statistical footing."
   ]
  }
 ],
 "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}