7-Generative-Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "78849b21",
   "metadata": {},
   "source": [
    "*This jupyter notebook is part of Arizona State University's course CAS 523 (Methods for Complex Systems Science: Statistics and Dimensionality Reduction) and was written by Bryan Daniels.  It was last updated September 28, 2022.*\n",
    "\n",
    "*This assignment uses data gathered by Ying Wang and Robert E. Page, Jr. at Arizona State University.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc6619fc",
   "metadata": {},
   "source": [
    "# Performing model selection using the Bayesian information criterion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cb19bb3",
   "metadata": {},
   "source": [
    "In this assignment, we will look again at the honey bee gene expression data that we previously analyzed using PCA.  Now going beyond a descriptive analysis, we will attempt to find a generative model that can output data that matches the statistics we measure in the data.  \n",
    "\n",
    "In this case, the inferred model will take the form of a probability distribution of gene expression values, predicting which combinations of gene expression are more or less likely.  A typical challenge when inferring such models is selecting the best form of model—in this case, selecting the form of the model sets the allowed shapes of the distribution.  Will a simple model suffice, or do we need to include more detail?  Here, we will use the \"Bayesian information criterion\" as a measure to decide which model is best."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af821c96",
   "metadata": {},
   "source": [
    "## Get set up and load the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5765354",
   "metadata": {},
   "source": [
    "**NOTE:** As in the previous PCA assignment, the dataset that we use in this assignment is not yet publicly available, but I have permission to share it with students.  For this reason, I do not include the data file in the public github repository.  Instead, it is available for download on Canvas.  Please follow the link on the Canvas assignment page to download the **two** files `Wang_Page_nanostring_data_2016_Day_1.csv` and `Wang_Page_nanostring_data_2016_Day_15.csv` and place them in the folder `data/WangPage2016/` in your own copy of the github repository.  Please do not share these data outside of this class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fc28530",
   "metadata": {},
   "source": [
    "Let's first load some standard packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6baeeaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams.update({'font.size': 18}) # increases font size on plots\n",
    "from pathlib import Path # to handle file paths across all operating systems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03bb6d0c",
   "metadata": {},
   "source": [
    "And some more specialized functions that we'll make use of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b50fe3e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "from sklearn.decomposition import PCA\n",
    "from helpers.prettyPlotting import scatter1D # custom 1D scatter plot\n",
    "from helpers.landau import LandauDistributionPDF # form of bimodal distribution that we fit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dff7804b",
   "metadata": {},
   "source": [
    "And now we'll load the gene expression data.  In the previous PCA assignment, we only looked at gene expression of bees at age 15 days.  Here, we'll also separately load gene expression profiles for bees at age 1 day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "113ce9cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataPathDay1 = Path('data/WangPage2016/Wang_Page_nanostring_data_2016_Day_1.csv')\n",
    "expressionDataDay1 = pd.read_csv(dataPathDay1).drop(columns=['Age'])\n",
    "\n",
    "dataPathDay15 = Path('data/WangPage2016/Wang_Page_nanostring_data_2016_Day_15.csv')\n",
    "expressionDataDay15 = pd.read_csv(dataPathDay15).drop(columns=['Age'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34079874",
   "metadata": {},
   "source": [
    "To remind ourselves what these data look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a403f6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "expressionDataDay1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad8bef2c",
   "metadata": {},
   "source": [
    "For each age, we have measurements of the expression levels of 91 genes for each of 16 bees."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b07b5f",
   "metadata": {},
   "source": [
    "## Reduce dimensionality"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "871607fc",
   "metadata": {},
   "source": [
    "Recall that, in the PCA assignment, we found what looked like two different groups of bees when we looked along the first principal component of the gene expression data.  Let's first reduce the dimensionality of our data by focusing only on this first component.  We'll do the calculation separately for bees of age 1 day and 15 days:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85003e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcaProjectionsDay1 = PCA(n_components=1).fit_transform(expressionDataDay1)\n",
    "pcaProjectionsDay15 = PCA(n_components=1).fit_transform(expressionDataDay15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8858ffe5",
   "metadata": {},
   "source": [
    "Recall also *why* it makes sense to look only along the first principal component: If we think the bees' gene expression is separating out into two distinct clusters, then the genes involved in that separation will have larger variance.  Focusing along the dimension with maximal variance—the first principal component—will highlight any such differences."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71155f03",
   "metadata": {},
   "source": [
    "The following makes a one-dimensional scatter plot of the Day 15 data as in our previous assignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c11eef6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter1D(pcaProjectionsDay15)\n",
    "plt.xlabel('Distance along first principal component')\n",
    "plt.title('Day 15');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5edc08ce",
   "metadata": {},
   "source": [
    "❓ **Make a 1D scatter plot for bees of age 1 day.  How do the two ages of bees look different?  By eye, would you guess that one or the other, or both, ages display distinct clusters of gene expression?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71d692b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✳️ **Answer:**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c5a7487",
   "metadata": {},
   "source": [
    "In the rest of the exercise, we will use statistical model selection to more precisely check our intuitive guess."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "100bb897",
   "metadata": {},
   "source": [
    "## Find best-fit Gaussian (unimodal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff10a443",
   "metadata": {},
   "source": [
    "The simplest guess for the distribution of gene expression along the principal component would be a Gaussian—the classic bell-shaped curve corresponding to a single continuous cluster.  This is called a unimodal distribution because it only has one peak.\n",
    "\n",
    "Recall from our exercise on quantifying uncertainty how to fit a Gaussian distribution to data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62a3e3cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "paramsNormalDay15 = stats.norm.fit(pcaProjectionsDay15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db33e2c",
   "metadata": {},
   "source": [
    "We can use the fit parameters to plot the probability density alongside a histogram of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91490180",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the best-fit Gaussian\n",
    "plt.hist(pcaProjectionsDay15,bins=10,density=True)\n",
    "xs = np.linspace(-15,15,100)\n",
    "plt.plot(xs,[stats.norm.pdf(x,*paramsNormalDay15) for x in xs],lw=5)\n",
    "plt.xlabel('Distance along first component')\n",
    "plt.ylabel('Probability density')\n",
    "plt.title('Day 15, Gaussian fit');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88cbdfe0",
   "metadata": {},
   "source": [
    "❓ **Looking at this plot, how well would you say the Gaussian does at fitting the data?**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31e6ca9f",
   "metadata": {},
   "source": [
    "✳️ **Answer:** "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9a377ec",
   "metadata": {},
   "source": [
    "Recall that the function `nnlf` returns the negative log-likelihood of the data given the model, so `-nnlf` is the standard log-likelihood:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15651b08",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaussianLogLDay15 = -stats.norm.nnlf(paramsNormalDay15,pcaProjectionsDay15)[0]\n",
    "print(\"The log-likelihood for the Gaussian model on Day 15 data is {}.\".format(gaussianLogLDay15))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57f7671f",
   "metadata": {},
   "source": [
    "We will use this next in comparing to a more complicated distribution with two bumps."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a46d6a6",
   "metadata": {},
   "source": [
    "## Find best-fit \"Landau\" distribution (bimodal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bfa2a2b",
   "metadata": {},
   "source": [
    "In a particular mathematical sense that I won't describe in detail here, the simplest distribution with *two* peaks (bimodal) is what I will call the \"Landau\" distribution.  This distribution corresponds to a very simplistic model of phase transitions in statistical physics. Here's the form of the distribution in case you are interested:\n",
    "\n",
    "$$\n",
    "p(x) = \\frac{1}{Z} \\exp \\left( -\\frac{c}{2} (x-\\mu)^2 - \\frac{d}{4} (x-\\mu)^4 \\right),\n",
    "$$\n",
    "\n",
    "where $Z$ is a normalization constant.  \n",
    "\n",
    "The Landau distribution has three parameters ($\\mu$, $c$, and $d$) compared to the Gaussian model's two parameters ($\\mu$ and $\\sigma$).  Note that the Gaussian is a special case of the Landau distribution: setting $d=0$ gives the usual form of a Gaussian.  So this model can have a single peak or two peaks depending on the values of the parameters.\n",
    "\n",
    "The Landau distribution is not standard enough to be implemented in the packages we use here, so I've done the fitting of parameters for you.  We load these parameters here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ecae352",
   "metadata": {},
   "outputs": [],
   "source": [
    "landauFitParameters = pd.read_csv(Path('data/WangPage2016/landau_fit_parameters.csv'),index_col='age')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba29481d",
   "metadata": {},
   "source": [
    "We can plot the fit distribution using the function `LandauDistributionPDF` that I wrote:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a030e57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the best-fit Landau distribution\n",
    "age = 15 # age of bees in days\n",
    "params = landauFitParameters.loc[age]\n",
    "plt.hist(pcaProjectionsDay15,bins=10,density=True)\n",
    "xs = np.linspace(-15,15,100)\n",
    "plt.plot(xs,LandauDistributionPDF(xs,params['mu'],params['c'],params['d']),lw=5)\n",
    "plt.xlabel('Distance along first component')\n",
    "plt.ylabel('Probability density')\n",
    "plt.title('Day 15, Landau fit');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6255c79b",
   "metadata": {},
   "source": [
    "This looks like a much better fit!  Still, it's not exactly right—but then again, we only have 16 datapoints, so we don't expect the fit to be exactly right.  How do we decide if it's better enough compared to the unimodal distribution?  This is a job for statistical model selection!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a6c0fa4",
   "metadata": {},
   "source": [
    "## Compute and interpret the Bayesian information criterion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b393ec4",
   "metadata": {},
   "source": [
    "First, let's look at the log-likelihood of the data given the Landau model.  I also have this pre-saved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c61a158",
   "metadata": {},
   "outputs": [],
   "source": [
    "landauLogLDay15 = landauFitParameters.loc[15]['log-likelihood']\n",
    "print(\"The log-likelihood for the Landau model on Day 15 data is {}.\".format(landauLogLDay15))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65866b32",
   "metadata": {},
   "source": [
    "❓ **Which model, Gaussian or Landau, fits the data better in terms of the log-likelihood?**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51eefaeb",
   "metadata": {},
   "source": [
    "✳️ **Answer:** "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1cfe990",
   "metadata": {},
   "source": [
    "❓ **Given that the Gaussian model is a special case of the Landau distribution, why does your answer to the last question have to be the case, no matter what the data were?** *Hint: Imagine you find the parameters of a simple model that best fit your data, and then you make a more complicated model by adding an extra parameter.  You add the new parameter in such a way that one of its possible values produces the same fit as the simpler model.  What do you know about the new fit?*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d588bcff",
   "metadata": {},
   "source": [
    "✳️ **Answer:** "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcd08afb",
   "metadata": {},
   "source": [
    "This logic means that our better fit with the more complicated model is not so impressive.  In other language we have encountered in this class, it's possible that the new, more complicated model is simply *overfitting*: fitting noise in the data.\n",
    "\n",
    "To determine if the extra parameter is \"worth the trouble\", we will compute the Bayesian information criterion (BIC).  Recall that the BIC includes a penalty on more complicated models, essentially setting a bar for *how much* better a model with extra parameters must fit the data to be statistically favored.\n",
    "\n",
    "As we saw in lecture, the BIC is defined as\n",
    "$$\n",
    "\\textrm{BIC} = L - \\frac{1}{2} N_{params} \\log(N_{datapoints}),\n",
    "$$\n",
    "where $L$ is the log-likelihood, $N_{params}$ is the number of free parameters in the model, and $N_{datapoints}$ is the number of datapoints to which the model was fit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f05825f",
   "metadata": {},
   "source": [
    "The following code computes BIC for the Gaussian and Landau models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9485f18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute bic for gaussian model\n",
    "numParamsGaussian = 2\n",
    "numDatapoints = 16\n",
    "gaussianBICDay15 = gaussianLogLDay15 - numParamsGaussian/2 * np.log(numDatapoints)\n",
    "\n",
    "# compute bic for landau model\n",
    "numParamsLandau = 3\n",
    "numDatapoints = 16\n",
    "landauBICDay15 = landauLogLDay15 - numParamsLandau/2 * np.log(numDatapoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea1cff09",
   "metadata": {},
   "source": [
    "Using this definition of BIC, having a BIC that is larger is indicative of evidence favoring that hypothesis;\n",
    "a BIC difference of more than 3 is considered strong evidence in favor of that hypothesis;\n",
    "and a BIC difference of more than 5 is considered very strong evidence.\n",
    "\n",
    "(See for instance Table 6 of [this reference](https://www.jstor.org/stable/271063).  Note that some authors define BIC to be twice the value we consider here, and/or use the opposite sign.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "886cb1c5",
   "metadata": {},
   "source": [
    "❓ **Which model is selected by BIC for bees of age 15 days?  Is one model strongly favored over the other?  How do you interpret this in terms of the number of distinct groups of bees?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26310c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✳️ **Answer:**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c70ae54",
   "metadata": {},
   "source": [
    "❓ **Repeat the analysis, including plots of the fit distributions and model selection using BIC, on data from bees of age 1 day.  Interpret your results.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db1cd5be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✳️ **Answer:**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}