02_Estimating_parameters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Welcome back!  In the [first post](http://www.davidketcheson.info/2020/03/17/SIR_model.html) of this series, we learned about the SIR model, which consists of three differential equations describing the rate of change of susceptible (S), infected (I), and recovered (R) populations:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{dS}{dt} & = -\\beta I \\frac{S}{N} \\\\\n",
    "\\frac{dI}{dt} & = \\beta I \\frac{S}{N}-\\gamma I \\\\\n",
    "\\frac{dR}{dt} & = \\gamma I\n",
    "\\end{align*}\n",
    "\n",
    "As we discussed, the model contains two key parameters ($\\beta$ and $\\gamma$) that influence the spread of a disease.  In this second post on modeling the COVID-19 outbreak, we will take the existing data and use it to estimate the values of those parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters we want to estimate are:\n",
    "\n",
    " - $\\beta$: The average number of people that come in close contact with a given infected individual, per day\n",
    " - $\\gamma$: The reciprocal of the average duration of the disease (in days)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating $\\gamma$\n",
    "\n",
    "A rough estimate of $\\gamma$ is available directly from medical sources.  Most cases of COVID-19 are mild and recovery occurs after about two weeks, which would give $\\gamma = 1/14 \\approx 0.07$.  However, a smaller portion of cases are more severe and can last for several weeks, so $\\gamma$ will be somewhat smaller than this value.  Estimates I have seen put the value in the range $0.03 - 0.06$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating $\\beta$\n",
    "\n",
    "It's much more difficult to get a good estimate of $\\beta$.  To be clear, we are trying to estimate, for an infected individual, the average number of other individuals with whom they have close contact per day.  Here *close contact* means contact that would lead to infection of the other individual (if that individual is still susceptible).\n",
    "\n",
    "As we discussed earlier, this number is affected by many factors.  It will also be affected by mitigation strategies implemented to reduce human contact.  For now, we want to estimate the value of $\\beta$ *in the absence of mitigation*.  Later, we will try to take mitigation into account."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that our equation for the number of infected is\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = \\left(\\beta \\frac{S}{N}-\\gamma \\right) I(t)\n",
    "$$\n",
    "\n",
    "Very early in an outbreak, the ratio $S/N \\approx 1$ since hardly anyone has\n",
    "been infected.  Also, at extremely early times, we can ignore $\\gamma$ because the\n",
    "disease is so new that nobody has been sick for long enough to recover.  For COVID-19,\n",
    "this is true for about the first two weeks of the disease' spread in a new population.\n",
    "During that time we have simply\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = \\beta I(t)\n",
    "$$\n",
    "\n",
    "This is one of the simplest differential equations, and its solution is just a growing exponential:\n",
    "\n",
    "$$\n",
    "I(t) = e^{\\beta t} I(0).\n",
    "$$\n",
    "\n",
    "Here $I(0)$ is of course the number of initially infected individuals.  Thus we can try to estimate $\\beta$ by fitting\n",
    "an exponential curve to the initial two weeks of spread.  This is not the only way to estimate $\\beta$; using this approach is the first of several choices that we'll make, and those choices will influence the our eventual predictions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting the data\n",
    "Fortunately for us, comprehensive data on the spread of COVID-19 is available from [this Github repository](https://github.com/CSSEGISandData/COVID-19) provided by the Johns Hopkins University Center for Systems Science and Engineering.  Specifically, I'll be using the data in [this file](https://github.com/CSSEGISandData/COVID-19/blob/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.csv).  Note that the file gets updated daily; as I write it is March 17th."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'm using Python and Pandas to work with the data.  For this blog post, I have removed most of the computer code, but you can [download the Jupyter notebook](https://github.com/ketch/covid-blog-posts/blob/master/02_Estimating_parameters.ipynb) and play with the code and data yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from scipy import optimize\n",
    "from ipywidgets import interact, widgets\n",
    "plt.style.use('seaborn-poster')\n",
    "matplotlib.rcParams['figure.figsize'] = (10., 6.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv\"\n",
    "data = pd.read_csv(url)\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "today = '3/17/20'  # Update this to include more data\n",
    "days = pd.date_range(start='1/22/20',end=today)\n",
    "dd = np.arange(len(days))\n",
    "\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')].sum() for day in days]\n",
    "print(total_cases[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate $\\beta$, we just pick a particular country from this dataset, plot the number of cases over time, and fit an exponential function to it.  We can use a standard mathematical tool called *least squares fitting* to find a reasonable value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the data from Italy\n",
    "For instance, here is the data from Italy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row=data['Country/Region'].tolist().index('Italy')\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][row] for day in days]\n",
    "plt.plot(total_cases);\n",
    "plt.xlabel('Days'); plt.ylabel('Total cases');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this data starts back in January, before the virus reached Italy, the number of cases at the beginning is zero.  We can use the interval from day 30 to day 43 (inclusive) to try to fit $\\beta$, since this seems to be when the outbreak began to take off.  Here it must be emphasized that the choice of this particular interval is somewhat arbitrary; different choices will give somewhat different values for $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exponential_fit(cases,start,length):\n",
    "\n",
    "    def resid(beta):\n",
    "        prediction = cases[start]*np.exp(beta*(dd-start))\n",
    "        return prediction[start:start+length]-cases[start:start+length]\n",
    "\n",
    "    soln = optimize.least_squares(resid,0.2)\n",
    "    beta = soln.x[0]\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how well this value predicts the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_fit(cases,start,end=56):\n",
    "    length=end-start\n",
    "    plt.plot(cases)\n",
    "    beta = exponential_fit(cases,start,length)\n",
    "    prediction = cases[start]*np.exp(beta*(dd-start))\n",
    "    plt.plot(dd[start:start+length],prediction[start:start+length],'--k');\n",
    "    plt.legend(['Data',r'fit with $\\beta$={:.3f}'.format(beta)]);\n",
    "    plt.xlabel('Days'); plt.ylabel('Total cases');\n",
    "    return beta\n",
    "    \n",
    "beta = plot_fit(total_cases,start=35,end=49)\n",
    "plt.title('Italy');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit seems reasonably good, over the interval we used.  How well does it match if we plot the fit over the whole time interval?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start=35\n",
    "plt.plot(total_cases)\n",
    "dd = np.arange(len(days))\n",
    "prediction = total_cases[start]*np.exp(beta*(dd-start))\n",
    "plt.plot(dd[start:],prediction[start:],'--k');\n",
    "plt.legend(['Data','fit']);\n",
    "plt.xlabel('Days'); plt.ylabel('Total cases');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the prediction is not accurate at later times.  There are two main reasons for this:\n",
    "\n",
    "  - Our assumption of exponential growth was based on other assumptions that are only valid at the very start of the outbreak;\n",
    "  - Italian society has taken measures to combat the spread of the virus, effectively reducing $\\beta$ at later times.\n",
    "\n",
    "We can resolve the first issue by using the full SIR model (instead of just exponential growth) to make predictions.  The second issue is more complicated; we will try to deal with it in a later blog post."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting to data from other regions\n",
    "\n",
    "To have more confidence in our value of $\\beta$, we can perform a similar fit with data from other regions, and see if we get a similar value.  Next, let's try fitting the data from the USA.  Here's the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'US'\n",
    "tophalf = data.iloc[:200]  # In lower part of file, duplicate data is given for cities\n",
    "\n",
    "rows = tophalf['Country/Region'].isin([country])\n",
    "total_cases = [tophalf[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plt.plot(dd,total_cases);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we only have about 1 week of meaningful data.  Let's try to fit an exponential to it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_fit(total_cases,start=48);\n",
    "plt.xlim(30,56); plt.title('USA');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a fairly similar value for $\\beta$.  Furthermore, the fit using this value seems to be pretty good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'Spain'\n",
    "rows = data['Country/Region'].isin([country])\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plot_fit(total_cases,start=48)\n",
    "plt.xlim(30,56); plt.title('Spain');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'United Kingdom'\n",
    "rows = data['Country/Region'].isin([country])\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plot_fit(total_cases,start=48)\n",
    "plt.xlim(30,56); plt.title('UK');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### France"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'France'\n",
    "rows = data['Country/Region'].isin([country])\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plot_fit(total_cases,start=44)\n",
    "plt.xlim(30,56); plt.title(country);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hubei Province, China\n",
    "Let's look at the data from where it all started: Hubei province, China.  Here it makes sense to start the fit from day zero of the JHU data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'Hubei'\n",
    "rows = data['Province/State'].isin([country])\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plot_fit(total_cases,start=0,end=14)\n",
    "plt.xlim(0,30); plt.title(country);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of these countries seems to fit the model reasonably well and to give a more or less similar value of $\\beta$, in the range $0.22$ to $0.29$.  It would be wrong to feel completely confident about this value, or to try to extrapolate too much from such a short time interval of data, but the consistency of these results does seem to suggest that our estimate is meaningful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at some countries that don't fit this pattern."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iran and South Korea\n",
    "Here is the number of confirmed cases for Iran:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'Iran'\n",
    "rows = data['Country/Region'].isin([country])\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plt.plot(dd[-30:],total_cases[-30:]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here is Korea:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'Korea, South'\n",
    "rows = data['Country/Region'].isin([country])\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "plt.plot(total_cases[:]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At a glance we can see that this data doesn't follow the pattern of the previous countries.  In Iran, after the first week, the growth seems to be linear.  In Korea, the initial exponential growth eventually slows down drastically and is beginning to level off.  This tells us that something we left out of our model must be at play.\n",
    "\n",
    "In the case of Korea, it seems straightforward to understand what is going on.  Korea has deployed the most extensive COVID-19 testing system in the world, with over 270,000 people tested to date.  This is combined with an extensive effort to isolate infected people and those they have been in recent contact with.  Essentially, South Korea has reduced the value of $\\beta$.  Based on our earlier analysis, to prevent future exponential growth, they will need to keep $\\beta$ down to approximately the value of $\\gamma$ or less.  If we believe that $\\gamma \\approx 0.05$ and $\\beta \\approx 0.25$, this means reducing the amount of human contact by infected people by five times. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iran's case is at first more puzzling, since the testing and quarantine measures there have not been exceptional compared to countries like Italy and Spain.  Instead, there are [strong](https://www.theatlantic.com/ideas/archive/2020/03/irans-coronavirus-problem-lot-worse-it-seems/607663/) [suspicions](https://www.nytimes.com/2020/02/28/world/middleeast/coronavirus-iran-confusion.html) that [the official numbers from Iran are wildly inaccurate](https://medicalxpress.com/news/2020-03-covid-outbreak-iran-larger.html) and the real number of cases (and deaths) is [drastically higher than what is reported](https://www.washingtonpost.com/world/middle_east/coronavirus-pummels-iran-leadership-as-data-show-spread-is-far-worse-than-reported/2020/03/04/7b1196ae-5c9f-11ea-ac50-18701e14e06d_story.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problems with the our approach\n",
    "\n",
    "Before we go finish, it's important to understand the limitations of the data we're working with and the technique we have used.  Most importantly, the numbers we have certainly **do not represent the real number of infected individuals**.  That's because many infected individuals are never tested for the virus.  This is especially true for diseases like COVID-19 in which the majority of cases are mild and do not require professional medical care.  Estimates I have seen claim that only about 10-20% of all cases are detected.\n",
    "\n",
    "If we assume that the fraction of cases that are actually detected is constant over time, then this discrepancy does not hinder our ability to estimate $\\beta$, since dividing the initial and final number of infected by the same constant will lead to the same estimate of $\\beta$ that would be obtained if we counted all the cases.  However, it's clear that in many places this factor changes over time as a country starts doing more and more testing.  This would cause the number of reported cases to grow even faster than the real number.  This is most likely occurring, for instance, in the US where previously many individuals with symptoms were not tested due to a lack of test availability.\n",
    "\n",
    "Another issue is that in some cases governments may be intentionally hiding the true number of infections.  As we have seen, this is likely the case in Iran.\n",
    "\n",
    "Finally, mitigation strategies may already be in place and influencing the rate of spread in some countries, even in the early days of outbreak.  This would lead to us underestimating the natural value of $\\beta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "\n",
    "What we can take away from this analysis are the rough estimates for the SIR parameters:\n",
    "\n",
    "$$\\gamma \\approx 0.05$$\n",
    "$$\\beta \\approx 0.25.$$\n",
    "\n",
    "Notice that the behavior in this initial phase of the epidemic that we have\n",
    "focused on is very similar to the simple behavior we considered at the start of\n",
    "the first post.  There, the number of infected individuals doubled each day, but we\n",
    "knew that was unrealistic.  Here, the number of infected individuals doubles every few days.\n",
    "How many days does it take for the number to double?\n",
    "If it takes $m$ days for the number of cases to double, then we have\n",
    "\n",
    "$$\n",
    "e^{\\beta m} = 2\n",
    "$$\n",
    "\n",
    "so $m = \\log(2)/\\beta$ where $\\log$ means the natural logarithm.  For $\\beta=0.25$, this gives a\n",
    "doubling time of 2.8 days.  This growth will slow down somewhat after the first\n",
    "couple of weeks for reasons we have already discussed.\n",
    "It should be emphasized that the value of $\\beta$ here is what we expect **in the absence of mitigation strategies**.   In later posts, we'll look at what these values mean for the future spread of the epidemic, and what the potential effect of mitigation may be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [Jupyter notebook for this post](https://github.com/ketch/covid-blog-posts/blob/master/02_Estimating_parameters.ipynb) there is an interactive setup where you can make your own fits to the data from a variety of regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_beta(region,interval):\n",
    "    start = interval[0]\n",
    "    end = interval[1]\n",
    "    if region in ['Hubei']:\n",
    "        rows = data['Province/State'].isin([region])\n",
    "    else:\n",
    "        rows = data['Country/Region'].isin([region])\n",
    "    total_cases = [data[day.strftime('%-m/%-d/%y')][rows].sum() for day in days]\n",
    "    beta = plot_fit(total_cases,start=start,end=end)\n",
    "    plt.xlim(start-10,end+10); plt.title(region);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(fit_beta,region=widgets.Dropdown(options=['Italy','Spain','France','Germany','United Kingdom','Hubei'],value='Italy'),\n",
    "                  interval=widgets.IntRangeSlider(min=0,max=len(days),value=(35,49)));"
   ]
  }
 ],
 "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.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}