bonus-1-network-statistical-inference-student.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Load the data\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import numpy.random as npr\n",
    "from scipy.stats import norm, ks_2samp  # no scipy - comment out\n",
    "from custom import load_data as cf\n",
    "from custom import ecdf\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction\n",
    "\n",
    "In this notebook, we will walk through a hacker's approach to statistical thinking, as applied to network analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Statistics in a Nutshell\n",
    "\n",
    "All of statistics can be broken down into two activities:\n",
    "\n",
    "- Descriptively summarizing data. (a.k.a. **\"descriptive statistics\"**)\n",
    "- Figuring out whether something happened by random chance. (a.k.a. **\"inferential statistics\"**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Descriptive Statistics\n",
    "\n",
    "- Centrality measures: mean, median, mode\n",
    "- Variance measures: inter-quartile range (IQR), variance and standard deviation\n",
    "\n",
    "### Inferential Statistics\n",
    "\n",
    "- Models of Randomness (see below)\n",
    "- Hypothesis Testing\n",
    "- Fitting Statistical Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Load Data\n",
    "\n",
    "Let's load a [protein-protein interaction network dataset](http://konect.uni-koblenz.de/networks/moreno_propro).\n",
    "\n",
    "> This undirected network contains protein interactions contained in yeast. Research showed that proteins with a high degree were more important for the surivial of the yeast than others. A node represents a protein and an edge represents a metabolic interaction between two proteins. The network contains loops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Read in the data.\n",
    "# Note from above that we have to skip the first two rows, and that there's no header column,and that the edges are\n",
    "# delimited by spaces in between the nodes. Hence the syntax below:\n",
    "G = cf.load_propro_network()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Compute some basic descriptive statistics about the graph, namely:\n",
    "\n",
    "- the number of nodes,\n",
    "- the number of edges,\n",
    "- the graph density,\n",
    "- the distribution of degree centralities in the graph,\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Number of nodes:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Number of edges:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Graph density:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Degree centrality distribution:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "How are protein-protein networks formed? Are they formed by an [Erdos-Renyi](https://en.wikipedia.org/wiki/Erdős–Rényi_model) process, or something else?\n",
    "\n",
    ">In the G(n, p) model, a graph is constructed by connecting nodes randomly. Each edge is included in the graph with probability p independent from every other edge.\n",
    "\n",
    "If protein-protein networks are formed by an E-R process, then we would expect that properties of the protein-protein graph would look statistically similar to those of an actual E-R graph."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Make an ECDF of the degree centralities for the protein-protein interaction graph, and the E-R graph.\n",
    "- The construction of an E-R graph requires a value for `n` and `p`. \n",
    "- A reasonable number for `n` is the number of nodes in our protein-protein graph.\n",
    "- A reasonable value for `p` might be the [density](https://www.quora.com/What-is-graph-density) of the protein-protein graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ppG_deg_centralities = _______\n",
    "plt.plot(*ecdf(__________))\n",
    "\n",
    "erG = nx.erdos_renyi_graph(n=_______, p=________)\n",
    "erG_deg_centralities = __________\n",
    "plt.plot(*ecdf(erG_deg_centralities))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "From visualizing these two distributions, it is clear that they look very different. How do we quantify this difference, and statistically test whether the protein-protein graph could have arisen under an Erdos-Renyi model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One thing we might observe is that the variance, that is the \"spread\" around the mean, differs between the E-R model compared to our data. Therefore, we can compare variance of the data to the distribtion of variances under an E-R model.\n",
    "\n",
    "This is essentially following the logic of statistical inference by 'hacking' (not to be confused with the statistical bad practice of p-hacking)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Fill in the skeleton code below to simulate 100 E-R graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# 1. Generate 100 E-R graph degree centrality variance measurements and store them.\n",
    "# Takes ~50 seconds or so.\n",
    "n_sims = ______\n",
    "er_vars = np.zeros(________)  # variances for n simulaed E-R graphs.\n",
    "for i in range(n_sims):\n",
    "    erG = nx.erdos_renyi_graph(n=____________, p=____________)\n",
    "    erG_deg_centralities = __________\n",
    "    er_vars[i] = np.var(__________)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# 2. Compute the test statistic that is going to be used for the hypothesis test.\n",
    "ppG_var = np.var(______________)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Do a quick visual check\n",
    "n, bins, patches = plt.hist(er_vars)\n",
    "plt.vlines(ppG_var, ymin=0, ymax=max(n), color='red', lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Visually, it should be quite evident that the protein-protein graph did not come from an E-R distribution. Statistically, we can also use the hypothesis test procedure to quantitatively test this, using our simulated E-R data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Conduct the hypothesis test.\n",
    "ppG_var > np.percentile(er_vars, 99)  # we can only use the 99th percentile, because there are only 100 data points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Another way to do this is to use the 2-sample Kolmogorov-Smirnov test implemented in the `scipy.stats` module. From the docs:\n",
    "\n",
    "> This tests whether 2 samples are drawn from the same distribution. Note\n",
    "that, like in the case of the one-sample K-S test, the distribution is\n",
    "assumed to be continuous.\n",
    ">\n",
    "> This is the two-sided test, one-sided tests are not implemented.\n",
    "The test uses the two-sided asymptotic Kolmogorov-Smirnov distribution.\n",
    ">\n",
    "> If the K-S statistic is small or the p-value is high, then we cannot\n",
    "reject the hypothesis that the distributions of the two samples\n",
    "are the same.\n",
    "\n",
    "As an example to convince yourself that this test works, run the synthetic examples below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Scenario 1: Data come from the same distributions.\n",
    "# Notice the size of the p-value.\n",
    "dist1 = npr.random(size=(100))\n",
    "dist2 = npr.random(size=(100))\n",
    "\n",
    "ks_2samp(dist1, dist2)\n",
    "# Note how the p-value, which ranges between 0 and 1, is likely to be greater than a commonly-accepted\n",
    "# threshold of 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Scenario 2: Data come from different distributions. \n",
    "# Note the size of the KS statistic, and the p-value.\n",
    "\n",
    "dist1 = norm(3, 1).rvs(100)\n",
    "dist2 = norm(5, 1).rvs(100)\n",
    "\n",
    "ks_2samp(dist1, dist2)\n",
    "# Note how the p-value is likely to be less than 0.05, and even more stringent cut-offs of 0.01 or 0.001."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Now, conduct the K-S test for one synthetic graph and the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Now try it on the data distribution\n",
    "ks_2samp(___________________, ___________________)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Networks may be high-dimensional objects, but the logic for inference on network data essentially follows the same logic as for 'regular' data:\n",
    "\n",
    "- Identify a model of 'randomness' that may model how your data may have been generated.\n",
    "- Compute a \"test statistic\" for your data and the model.\n",
    "- Compute the probability of observing the data's test statistic under the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Further Reading\n",
    "\n",
    "Jake Vanderplas' \"Statistics for Hackers\" slides: https://speakerdeck.com/jakevdp/statistics-for-hackers\n",
    "\n",
    "Allen Downey's \"There is Only One Test\": http://allendowney.blogspot.com/2011/05/there-is-only-one-test.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Environment (conda_nams)",
   "language": "python",
   "name": "conda_nams"
  },
  "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"
  },
  "widgets": {
   "state": {},
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}