incubator/IC50_synthetic-bayesian.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pymc3 as pm\n",
    "import pandas as pd\n",
    "import theano.tensor as tt\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Type\n",
    "\n",
    "IC<sub>50</sub> determination is a commonplace problem type for experimental biochemists. The goal is to determine the concentration of an inhibitor that causes a protein to have 1/2 as much activity as without the inhibitor.\n",
    "\n",
    "Classically, this is done by fitting a logistic function of some sort, where the IC<sub>50</sub> is a parameter in the model. We then want to determine IC<sub>50</sub>, given the data at hand.\n",
    "\n",
    "In the implementation here, we allow for multi-drug IC<sub>50</sub> determination. Additionally, I also show how we can have a mix of noisy data and have non-uniform coverage over the drug concentration space, and still come to a fairly precise determination of the IC<sub>50</sub> values as long as we have dense enough data around the true IC50 area. The simulated data below assumes that a first-pass experiment was done over a wide range of data points, followed by a second-pass set of experiments where a range of concentrations was filled in with more resolution.\n",
    "\n",
    "## Data structure\n",
    "\n",
    "To use it with this model, the data should be structured as such:\n",
    "\n",
    "- Each row is one measurement for one concentration of drug.\n",
    "- The columns should indicate, at the minimum:\n",
    "    - What treatment group the sample belonged to.\n",
    "    - The measured value.\n",
    "\n",
    "## Extensions to the model\n",
    "\n",
    "None\n",
    "\n",
    "## Reporting summarized findings\n",
    "\n",
    "Here are examples of how to summarize the findings.\n",
    "\n",
    "> The IC<sub>50</sub> of drug X was `mean` (95% HPD: [`lower`, `upper`]).\n",
    "\n",
    "## Other notes\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate some fake data.\n",
    "x_concs = np.concatenate(\n",
    "    [\n",
    "        np.arange(0, 100, 10),\n",
    "        np.arange(5, 21, 1),\n",
    "        np.arange(40, 51, 1),\n",
    "        np.arange(85, 96, 1),\n",
    "    ]\n",
    ").reshape(-1, 1)\n",
    "log_xconcs = np.log(x_concs)\n",
    "ic50_true = np.array([42, 13, 88])\n",
    "ic50_true = ic50_true.reshape(-1, ic50_true.shape[0])\n",
    "beta_true = 1\n",
    "slope_true = 1\n",
    "intercept_true = 150\n",
    "# y_true = slope_true * x_concs + intercept_true\n",
    "\n",
    "y_true = beta_true / (1 + np.exp(x_concs - ic50_true))\n",
    "\n",
    "y_noisy = y_true + np.random.normal(0, 0.15, size=y_true.shape)  # homoskedastic error\n",
    "y_noisy[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_concs, y_noisy[:, 0])\n",
    "plt.scatter(x_concs, y_noisy[:, 1])\n",
    "plt.scatter(x_concs, y_noisy[:, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how there's sparse coverage in the regions outside of the true IC<sub>50</sub>s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "concentrations = np.concatenate([x_concs] * ic50_true.shape[1])\n",
    "concentrations[::10]  # print every 10th"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "concentrations.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_noisy.flatten(order=\"F\").shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_concs.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.DataFrame()\n",
    "data[\"concentrations\"] = concentrations.reshape(concentrations.shape[0],)\n",
    "data[\"measurements\"] = y_noisy.flatten(order=\"F\")\n",
    "\n",
    "drugs = []\n",
    "for i in range(ic50_true.shape[1]):\n",
    "    drugs.extend([i] * x_concs.shape[0])\n",
    "data[\"drug\"] = drugs\n",
    "\n",
    "data = pm.floatX(data)\n",
    "\n",
    "\n",
    "from sklearn.preprocessing import LabelEncoder, MinMaxScaler\n",
    "\n",
    "le = LabelEncoder()\n",
    "data[\"idxs\"] = le.fit_transform(data[\"drug\"]).astype(\"int32\")\n",
    "\n",
    "# Normalize data['measurements'] to 0-1\n",
    "# mms = MinMaxScaler()\n",
    "# data['measurements'] = mms.fit_transform(data['measurements'].values.reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head().dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    datasize = len(set(data[\"idxs\"]))\n",
    "    beta = pm.Normal(\"beta\", mu=0, sd=100 ** 2, shape=datasize)\n",
    "    noise = pm.HalfCauchy(\"noise\", beta=100 ** 2, shape=datasize)\n",
    "    ic50 = pm.Normal(\"IC50\", sd=100 ** 2, shape=datasize)\n",
    "    measurements = pm.Deterministic(\n",
    "        \"measurements\",\n",
    "        beta[data[\"idxs\"].values]\n",
    "        / (1 + tt.exp(data[\"concentrations\"].values - ic50[data[\"idxs\"].values])),\n",
    "    )\n",
    "\n",
    "    y_like = pm.Normal(\n",
    "        \"y_like\",\n",
    "        mu=measurements,\n",
    "        sd=noise[data[\"idxs\"].values],\n",
    "        observed=data[\"measurements\"],\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    trace = pm.sample(draws=10000, step=pm.Metropolis(), start=pm.find_MAP())\n",
    "    # trace = pm.sample(draws=2000)  # ADVI init is fast, but NUTS sampling is slow later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.traceplot(trace[5000:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the posterior distributions around the IC50s are fairly nice and tight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.summary(trace[5000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bayesian",
   "language": "python",
   "name": "bayesian"
  },
  "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.6"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "102px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}