BayesianCovid.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trying out Bayesian inference with PyMC3 on covid data\n",
    "\n",
    " _Disclaimer: this is in no way intended to be relied on!_\n",
    " _this was done purely for me to learn something_\n",
    " \n",
    "It doesn't respect reactions of the countries, it doesn't respect the testing capabilities / numbers in the countries, it doesn't respect real biological models and past research in the field of virology and pandemics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
     ]
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import theano\n",
    "%matplotlib inline\n",
    "import warnings\n",
    "from scipy.stats import halfnorm\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data based on a dump from a wiki page offering country specific infections.\n",
    "\n",
    "Data is a snapshot form Kaggle taken from around mid April 2020 and wasn't updated since!\n",
    "\n",
    "To make the data more representative, days before 2000 infections were reached were removed, since there might have been just single hotspots that were under control, also only those timeseries were looked at, that had in it's current state more than 20.000 infections counted.\n",
    "Furthermore the data was restricted to series of at least 10 days.\n",
    "These restrictions allow to look at a smaller set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0, 'Austria-'), (1, 'Belarus-'), (2, 'Belgium-'), (3, 'Brazil-'), (4, 'Canada-Ontario'), (5, 'Canada-Quebec'), (6, 'Chile-'), (7, 'China-Hubei'), (8, 'Ecuador-'), (9, 'France-'), (10, 'Germany-'), (11, 'India-'), (12, 'Iran-'), (13, 'Ireland-'), (14, 'Israel-'), (15, 'Italy-'), (16, 'Japan-'), (17, 'South Korea-'), (18, 'Mexico-'), (19, 'Netherlands-'), (20, 'Pakistan-'), (21, 'Peru-'), (22, 'Poland-'), (23, 'Portugal-'), (24, 'Qatar-'), (25, 'Romania-'), (26, 'Russia-'), (27, 'Saudi Arabia-'), (28, 'Singapore-'), (29, 'Spain-'), (30, 'Sweden-'), (31, 'Switzerland-'), (32, 'Turkey-'), (33, 'US-California'), (34, 'US-Colorado'), (35, 'US-Connecticut'), (36, 'US-Florida'), (37, 'US-Georgia'), (38, 'US-Illinois'), (39, 'US-Indiana'), (40, 'US-Louisiana'), (41, 'US-Maryland'), (42, 'US-Massachusetts'), (43, 'US-Michigan'), (44, 'US-New Jersey'), (45, 'US-New York'), (46, 'US-Ohio'), (47, 'US-Pennsylvania'), (48, 'US-Texas'), (49, 'US-Virginia'), (50, 'US-Washington'), (51, 'United Arab Emirates-'), (52, 'United Kingdom-')]\n"
     ]
    }
   ],
   "source": [
    "infections = []\n",
    "countries = {}\n",
    "\n",
    "MIN_DATES = 10\n",
    "\n",
    "with open('untitled1.txt', 'r') as csv:\n",
    "    intermediate = []\n",
    "    counter = 0\n",
    "    \n",
    "    for line in csv:\n",
    "        \n",
    "        line = line.strip().split(',')\n",
    "        country = line[2]+'-'+line[1]\n",
    "        infection = int(float(line[4]))\n",
    "        deaths = int(float(line[5]))\n",
    "#        print(line)\n",
    "                \n",
    "        if infection < 2000:\n",
    "            continue\n",
    "        if not country in countries:\n",
    "            countries[country] = 0\n",
    "            counter = 0\n",
    "            if len(intermediate) > MIN_DATES and intermediate[-1][2] > 10000:\n",
    "                for i in intermediate:\n",
    "                    infections.append(i)\n",
    "                \n",
    "            intermediate = []\n",
    "        counter += 1\n",
    "        intermediate.append([country, counter, infection, deaths])\n",
    "    if len(intermediate) > MIN_DATES:\n",
    "        for i in intermediate:\n",
    "            infections.append(i)\n",
    "\n",
    "\n",
    "        \n",
    "full_df = None\n",
    "full_df = pd.DataFrame(infections, columns=['country', 'day', 'infections', 'deaths'])\n",
    "full_df = full_df.astype({'day': 'int32', 'infections': 'int32', 'deaths': 'int32'})\n",
    "\n",
    "#filters = full_df.country.apply(lambda x: x in [\n",
    "#        'China', 'Germany', 'Japan', 'South Korea', 'France', 'Netherlands'])\n",
    "#full_df=full_df[filters]\n",
    "countries = full_df.country.values\n",
    "uniq_countries = full_df.country.unique()\n",
    "n_countries = len(uniq_countries)\n",
    "\n",
    "full_df['country_idx'] = [list(uniq_countries).index(x) for x in countries]\n",
    "#print(full_df.country_idx)\n",
    "#print(full_df)\n",
    "print(list(enumerate(uniq_countries)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### here is the modeling part\n",
    "\n",
    "the base idea is to fit a sigmoid like function to model the number of total infections. This assumption alone is probably already enough reason to not trust any output of this model. So _please don't trust_ the model.\n",
    "\n",
    "Instead of using the regular sigmoid, I chose the _Gompertz Function_:\n",
    "\n",
    "\\begin{equation}\n",
    "\\large{\n",
    "f(x) = a \\cdot e^{b \\cdot e^{c \\cdot x} }\n",
    "}\n",
    "\\end{equation}\n",
    "\n",
    "The reason for using the Gompertz function is it's assymmetrie, allowing to adjust for the exponential increase ans slow down phases.\n",
    "with $b, c < 0$ the value of $a$ determines the upper limit and therefore in our investigation the upper limit of infections.\n",
    "$b$ and $c$ determine the speeed and acceleration.\n",
    "\n",
    "To have some benefit from all the past countries, I tried to model $b$ and $c$ hierarchical, having a \"mean value\" across all time series, and the individual time series deviates from this according to a small normal distribution. The idea is, to have estimates for how things will develop even when very little hints are in the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [c, c_base, b, b_base, a]\n",
      "Sampling 2 chains, 0 divergences:   4%|▍         | 3538/81000 [02:59<1:04:25, 20.04draws/s]"
     ]
    }
   ],
   "source": [
    "from theano import shared\n",
    "predictors = full_df.day.values.copy()\n",
    "predictors_shared = shared(predictors)\n",
    "country_id = full_df.country_idx.values.copy()\n",
    "country_idx = shared(country_id)\n",
    "\n",
    "\n",
    "from theano import shared\n",
    "predictors = full_df.day.values.copy()\n",
    "predictors_shared = shared(predictors)\n",
    "\n",
    "\n",
    "\n",
    "import scipy\n",
    "with pm.Model() as model:\n",
    "    \n",
    "    a = pm.Uniform('a', lower=1000, upper=2000000, shape=n_countries)\n",
    "    b_base = pm.Normal('b_base', mu=-4.5, sigma=0.5)\n",
    "    b = pm.Normal('b', mu=b_base, sigma=0.5, shape=n_countries)\n",
    "    c_base = pm.Normal('c_base', mu=-0.075, sigma=0.03)\n",
    "    c = pm.Normal('c', mu=c_base, sigma=0.03, shape=n_countries)\n",
    "    \n",
    "    y = (a[country_idx] * pm.math.exp(b[country_idx] * pm.math.exp(c[country_idx] * (predictors_shared))))\n",
    "    obs = pm.Normal('obs', mu=y, sigma=15000, observed=full_df.infections.values)\n",
    "\n",
    "    trace = pm.sample(40000, cores=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now plotting the results of the fittings\n",
    "\n",
    "The fittings did not work out very well, we will see why when we look at the traces.\n",
    "\n",
    "We can see some pretty wide confidence intervals, so like the output suggested it didn't work out too well.\n",
    "Interestingly this is especially then the case, when the counts haven't turned into the slow down phase where the infections are under control. This also makes sense, because the model has to guess which kind of behavior it will see when the infections get under control, without having any hints on it.\n",
    "But here is the hierarchical model at least helping a bit, interpolating from overal behavior of the infections to the individual case.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from pymc3 import forestplot\n",
    "\n",
    "plt.figure(figsize=(20,20))\n",
    "forestplot(trace, var_names=['a'])\n",
    "forestplot(trace, var_names=['b'])\n",
    "forestplot(trace, var_names=['c'])\n",
    "\n",
    "pm.traceplot(trace)\n",
    "print(list(enumerate(uniq_countries)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### now predicting the future...\n",
    "\n",
    "the traceplot above show what we already assumed, had some issues, especially the base values of c and b didn't fully converge to a single distribution, normally you would do a reparametrization and probably increase tuning steps to fix this.\n",
    "But still let us try to now use the found model parameters to simulate how it's going to continue.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ppc = pm.sample_posterior_predictive(trace, samples=500, model=model)\n",
    "x = np.tile(np.linspace(1, 100, 100).astype('int32'), n_countries)\n",
    "print(len(x))\n",
    "predictors_shared.set_value(x)\n",
    "\n",
    "y = np.repeat(np.linspace(0,n_countries-1,n_countries).astype('int32'), 100)\n",
    "print(len(y))\n",
    "country_idx.set_value(y)\n",
    "\n",
    "with model:\n",
    "    post_pred = pm.sample_posterior_predictive(trace, samples=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### looking at fittings and predictions\n",
    "What we can actually see is that the model fitted the given points quite ok, but the predictions have quite a lot uncertainty. Especially in those cases, where there is little hint as to how much the region was able to slow down. \n",
    "So again don't rely on this model for anything.\n",
    "This was done purely as an educational exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "means = post_pred['obs'].mean(axis=0, keepdims=False).copy()\n",
    "stds = post_pred['obs'].std(axis=0)\n",
    "for i in range(n_countries):\n",
    "    choice = y==i\n",
    "    old_choice = full_df.country_idx==i\n",
    "\n",
    "    plt.figure(figsize=(10,10))\n",
    "    plt.errorbar(np.linspace(1,100,100),\n",
    "                 means[choice],\n",
    "                 stds[choice],\n",
    "                 linestyle='None',\n",
    "                 marker='.')\n",
    "    plt.plot(np.linspace(1,len(full_df[old_choice]), len(full_df[old_choice])),\n",
    "             full_df.infections[old_choice],\n",
    "             marker='o')\n",
    "    plt.title(uniq_countries[i])\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}