{ "cells": [ { "cell_type": "markdown", "id": "0bc94814", "metadata": {}, "source": [ "# How bad is the 4th wave going to be in Germany?" ] }, { "cell_type": "markdown", "id": "9dbeaf47", "metadata": {}, "source": [ "We explore two scenarios for the next month:\n", "\n", "- Option 1: The R-value stays at 1.25 (as it has been for a few days now). This is exponential growth.\n", "\n", "- Option 2: We see a linear increase in the *daily* new infections for the past week: we assume this trend continues - this is quadratic growth in the total number of infections. (This would reflect an R value that decreases over time.) \n", "\n", " This assumption is based on the good fit (see relevant plot below), and the experience that a similar fit was a fairly good prediction in the early stages of the outbreak in China. There is no other good reason for investigating this assumption.\n" ] }, { "cell_type": "markdown", "id": "2825976b", "metadata": {}, "source": [ "Outline of document:\n", "\n", "- we compute the daily new infections as an average over the last week (to get rid of reporting fluctuations)\n", "\n", "- we define some convenience functions\n", "\n", "- we show data for Option 1 and Option 2\n", "\n", "- we close with a discussion" ] }, { "cell_type": "code", "execution_count": 1, "id": "77fc39f5", "metadata": {}, "outputs": [ { "data": { "text/markdown": [ "[Execute this notebook with Binder](https://mybinder.org/v2/gh/oscovida/binder/master?filepath=ipynb/2021-11-12-germany-onset-wave4.ipynb)" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib.dates as mdates\n", "import pandas as pd\n", "import numpy as np\n", "import oscovida as ov\n", "%config InlineBackend.figure_formats = ['svg']\n", "# ov.clear_cache()\n", "ov.display_binder_link(\"2021-11-12-germany-onset-wave4.ipynb\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "de3a8995", "metadata": {}, "outputs": [], "source": [ "cases_jh, deaths_jh = ov.get_country_data(\"Germany\")\n", "\n", "# compute weekly mean (average over the last 7 days) - \n", "# this is proportional to the incidence\n", "c = cases_jh.rolling(window=7, center=False).mean()\n", "\n", "# infections per day (averaged over 7 days), for the last 42 days only\n", "daily_cases = c.diff()[-42:]" ] }, { "cell_type": "code", "execution_count": 3, "id": "bacbfbb7", "metadata": {}, "outputs": [], "source": [ "# We define one more conveniene function to plot our \n", "# results (please ignor the code)\n", "\n", "def daily_growth_to_incidence(g, country='Germany'):\n", " \"\"\"Given a time series g of new infections per day, \n", " return the corresponding incidence (for country)\"\"\"\n", " pop = ov.population(country)\n", " \n", " # 7 days * g new cases\n", " new_cases = 7 * g\n", " # divide by population and multiply with 100k\n", " return new_cases / pop * 1e5\n", "\n", "assert abs(daily_growth_to_incidence(33000) - 278) < 1\n", "def plot_predictions(series :pd.Series, title):\n", " \"\"\"Expect series to be series with daily \n", " new cases and date as index.\"\"\"\n", " assert isinstance(series, pd.Series)\n", " \n", " fig, ax = plt.subplots()\n", " incidence = daily_growth_to_incidence(series)\n", " \n", " ax.bar(series.index, series)\n", " fig.autofmt_xdate()\n", " \n", " # hacks to show two y-scales\n", " ax2 = ax.twinx()\n", " ax2.plot(series.index, incidence, alpha=0.5)\n", " ax2.set_ylim(bottom=0)\n", " ax.set_ylim(top=series.max()*1.05)\n", " ax2.set_ylim(top=incidence.max()*1.05)\n", " ax2.grid(False)\n", " \n", " # hacks for show second x-axis at the top\n", " ax3 = ax.twiny()\n", " ax3.grid(False)\n", " ax3.plot(mdates.date2num(series.index - series.index[0]), alpha=0)\n", "\n", " # labels\n", " ax.set_ylabel('daily new cases')\n", " ax2.set_ylabel('incidence')\n", " ax3.set_xlabel('days in the future')\n", " \n", " ax.set_title(title)\n", " \n", " print(f\"Total infections in this data set: {round(series.sum())}\")" ] }, { "cell_type": "markdown", "id": "bb6068e1", "metadata": {}, "source": [ "# Option 1: Assuming exponential growth with constant R" ] }, { "cell_type": "code", "execution_count": 4, "id": "e9b2992f", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-11-14T23:36:15.511107\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.4.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-08\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-22\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-08\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0.8\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0.9\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1.0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1.1\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1.2\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1.3\n", " \n", " \n", " \n", " R & growth factor\n", " (based on )\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " daily growth factor \n", " \n", " \n", " \n", " \n", " \n", " \n", " daily growth factor (rolling mean)\n", " \n", " \n", " \n", " \n", " \n", " \n", " estimated R (using )\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ov.plot_reproduction_number(ax=ax, series=cases_jh[-42:])\n", "fig.autofmt_xdate()" ] }, { "cell_type": "markdown", "id": "95f68607", "metadata": {}, "source": [ "Looking at the yellow line, we see an R value of approximately 1.25 since the beginning of November. To extract this number from the reported infections per day, we have assumed reproduction interval of 4 days: i.e. every 4 days each currently infected person will infect 1.25 people. (See https://oscovida.github.io/r-value.html for more details.)" ] }, { "cell_type": "markdown", "id": "2fceaf01", "metadata": {}, "source": [ "Current infection numbers per day (based on data averaged over the last week) are" ] }, { "cell_type": "code", "execution_count": 5, "id": "a9299552", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "36766" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "daily_cases[-1].astype('int')" ] }, { "cell_type": "markdown", "id": "094928f7", "metadata": {}, "source": [ "If the current exponential growth does not change, we expect the following number of daily infections for the next 4 weeks:\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "d42ac3ec", "metadata": {}, "outputs": [], "source": [ "r = 1.25\n", "int = 4 # serial interval (time re-infect next generation)\n", "n0 = float(daily_cases[-1])\n", "days = np.linspace(1, 28, 28)\n", "new_cases = n0*r**(days/4)\n", "index = days.astype('timedelta64[D]') + daily_cases.index[-1]\n", "option1 = pd.DataFrame(index=index, data={'new cases': new_cases, \n", " 'days' : days})\n", " " ] }, { "cell_type": "code", "execution_count": 7, "id": "cc93ce51", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total infections in this data set: 2553504\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-11-14T23:36:15.868132\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.4.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-13\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-17\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-21\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-25\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-29\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-05\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-09\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 25000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 50000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 75000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 100000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 125000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 150000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 175000\n", " \n", " \n", " \n", " daily new cases\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " Prediction Option 1:\n", " Exponential growth: assume R stays at 1.25\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 200\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 400\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 600\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 800\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1200\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 1400\n", " \n", " \n", " \n", " incidence\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 5\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 10\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 20\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 25\n", " \n", " \n", " \n", " days in the future\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_predictions(\n", " option1['new cases'], \n", " \"Prediction Option 1:\\n Exponential growth: assume R stays at 1.25\")" ] }, { "cell_type": "markdown", "id": "a7808326", "metadata": {}, "source": [ "So for ongoing exponential growth, we expect up to 175,000 new infections per day in Germany in 28 days - corresponding to an incidence of 1400. Over 28 days, we look at more than 2.5 million new infections (See discussion below.) " ] }, { "cell_type": "markdown", "id": "15f55724", "metadata": {}, "source": [ "# Option 2: assume linear growth in daily new infections" ] }, { "cell_type": "code", "execution_count": 8, "id": "e487318c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'Last data point is from 2021-11-13.'" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-11-14T23:36:16.122944\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.4.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-08\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-22\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-08\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 5000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 10000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 15000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 20000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 25000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 30000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 35000\n", " \n", " \n", " \n", " new daily cases\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " new cases per day (averaged over last week)\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.bar(daily_cases.index, daily_cases, color='C6')\n", "fig.autofmt_xdate()\n", "ax.set_title(\"new cases per day (averaged over last week)\")\n", "ax.set_ylabel(\"new daily cases\")\n", "f\"Last data point is from {daily_cases.index[-1].date()}.\"" ] }, { "cell_type": "markdown", "id": "633eb0f5", "metadata": {}, "source": [ "We see somewhat constant daily new infections up to the middle October 2021, then the daily new infections increase significantly. (The incidence is proportional to this.)" ] }, { "cell_type": "markdown", "id": "6677abac", "metadata": {}, "source": [ "### Estimate growth rate of new daily infections by using linear fit for individual weeks" ] }, { "cell_type": "code", "execution_count": 9, "id": "e0414904", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-11-14T23:36:16.345447\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.4.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-08\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-10-22\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-08\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 5000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 10000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 15000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 20000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 25000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 30000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 35000\n", " \n", " \n", " \n", " daily new cases\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " week \"-1\"\n", " \n", " \n", " \n", " \n", " \n", " \n", " week \"-2\"\n", " \n", " \n", " \n", " \n", " \n", " \n", " week \"-3\"\n", " \n", " \n", " \n", " \n", " \n", " \n", " week \"-4\"\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.bar(daily_cases.index, daily_cases, color='C6', alpha=1)\n", "fit_coef = {}\n", "fit_func = {}\n", "nweeks = 4\n", "for i in range(nweeks):\n", " # indices for days to look at one week at a time\n", " start = -7*(i+1)\n", " stop = -7*i\n", " if stop == 0: stop = None\n", " # convert dates to number of days for fitting routine\n", " X = mdates.date2num(daily_cases.index[start:stop]) \n", " Y = daily_cases[start:stop] \n", " # fit affine linear function to one week of data points\n", " ps = np.polyfit(X, Y, 1)\n", " m, a = ps # extract slope m and offset a \n", " fit_coef[i] = (m, a) # to display later\n", " fit = np.poly1d(ps)\n", " fit_func[i] = fit # re-use later in the notebook\n", "\n", " ax.plot(Y.index, fit(X), linewidth=5, alpha=0.5, label=f'week \"-{i+1}\"')\n", "\n", "fig.autofmt_xdate()\n", "ax.legend()\n", "ax.set_ylabel(\"daily new cases\");\n" ] }, { "cell_type": "markdown", "id": "4d18e709", "metadata": {}, "source": [ "We can see that the slope of increase tends to increase from week to week. \n", "\n", "In the following, we will use the linear trend shown with the red thick line (for week \"-1\": that is 1 week in the past) as the **growth model for Option 2**.\n", "\n", "The only justification for this choice is that the fit looks good (and that a similar approach worked well for the prediction of new infections in China late 2019/early 2020)." ] }, { "cell_type": "code", "execution_count": 10, "id": "985e7ee6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{0: (1766.4693877537343, -33427013.714259814),\n", " 1: (1095.1836734680578, -20714946.816301342),\n", " 2: (1069.4030612242636, -20224457.71938348),\n", " 3: (790.515306122628, -14944854.306125835)}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fit_coef" ] }, { "cell_type": "markdown", "id": "8dd5bc3a", "metadata": {}, "source": [ "# Assume number of *daily new infections* stays constant" ] }, { "cell_type": "markdown", "id": "780945f7", "metadata": {}, "source": [ "Based on the last week, we have new infections per day:" ] }, { "cell_type": "code", "execution_count": 11, "id": "2092d4a3", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
new casesincidence
2021-11-1438,749326
2021-11-2151,114430
2021-11-2863,479534
2021-12-0575,845638
2021-12-1288,210743
\n", "
" ], "text/plain": [ " new cases incidence\n", "2021-11-14 38,749 326\n", "2021-11-21 51,114 430\n", "2021-11-28 63,479 534\n", "2021-12-05 75,845 638\n", "2021-12-12 88,210 743" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# create data frame for convencience\n", "nc1 = [] # new cases\n", "model = fit_func[0]\n", "futuredays = np.arange(1, 30, dtype='int') + mdates.date2num(daily_cases.index[-1])\n", "new_cases = model(futuredays)\n", "incidence = daily_growth_to_incidence(new_cases)\n", "index = pd.to_datetime(mdates.num2date(futuredays)).date\n", "option2 = pd.DataFrame(index=index, data={'new cases' : new_cases,\n", " 'incidence' : incidence})\n", "\n", "pd.options.display.float_format = '{:,.0f}'.format\n", "option2.iloc[::7, :]" ] }, { "cell_type": "code", "execution_count": 12, "id": "4820ccbc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total infections in this data set: 1840903\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-11-14T23:36:16.759948\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.4.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-13\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-17\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-21\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-25\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-11-29\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-01\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-05\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-09\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2021-12-13\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 20000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 40000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 60000\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 80000\n", " \n", " \n", " \n", " daily new cases\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " Prediction:\n", " Option 2\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 100\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 200\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 300\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 400\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 500\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 600\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 700\n", " \n", " \n", " \n", " incidence\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 0\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 5\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 10\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 15\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 20\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 25\n", " \n", " \n", " \n", " days in the future\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_predictions(option2['new cases'], \n", " \"Prediction:\\n Option 2\")" ] }, { "cell_type": "markdown", "id": "f5cc4c99", "metadata": {}, "source": [ "With this assumption, we would see up to 88,000 new infections per day and an incidence of 700. That is a total of 1.8million new infections over 28 days." ] }, { "cell_type": "markdown", "id": "d6781fe7", "metadata": {}, "source": [ "# Discussion" ] }, { "cell_type": "markdown", "id": "db38d75d", "metadata": {}, "source": [ "### Summary data\n", "\n", "Using exponetial growth with R=1.25, we find an incidence up to 1400, up to 175,000 new infections per day, and 2.5million of infections over 28 days.\n", "\n", "Using Option2, we find an incidence up to 700, up to 88,000 new infections per day, and 1.8million of infections over 28 days.\n" ] }, { "cell_type": "markdown", "id": "ab6c1e53", "metadata": {}, "source": [ "### Importance" ] }, { "cell_type": "markdown", "id": "f94e1556", "metadata": {}, "source": [ "The number of infections will lead to hospital cases and deaths. For Germany, there were about 5 million infections so far with 100,000 deaths, leading to a rate of about 2% of those infected dying:" ] }, { "cell_type": "code", "execution_count": 13, "id": "49a7dce0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cases : 5037039\n", "deaths : 97677\n", "death rate : 1.94%\n" ] } ], "source": [ "print(f\"cases : {cases_jh[-1]}\\n\"\n", " f\"deaths : {deaths_jh[-1]}\\n\"\n", " f\"death rate : {deaths_jh[-1]/cases_jh[-1]*100:.3}%\")" ] }, { "cell_type": "markdown", "id": "3447b9f8", "metadata": {}, "source": [ "This death rate of 2% is likely to go down now that large parts of the population have been vaccinated, and thus have a much reduced risk of sever illness or death. On the other hand, there is still a very large population of not vaccinated people (1/4 to 1/3 of the population - so still representing a smaller country) which will experience the death rate as before when they get infected.\n", "\n", "\n", "**We can estimate an order of magnitude of the impact.**\n", "\n", "Let's assume a death rate of 0.4% only [1]. In that case, we \n", "\n", "- expect 200 new deaths for every 50,000 new infections \n", "\n", "- the worst case scenario of unlimited exponentional growth suggests 2.5million infections: this would correspond to 2.5m*0.4% = 10,000 new deaths in total." ] }, { "cell_type": "markdown", "id": "011d5ae1", "metadata": {}, "source": [ "About half of the patients going to an intensive care unit will eventually die. We can thus estimate the number of intensive care patients as twice that of the deaths (0.7% [1]): \n", "\n", "- With _every day_ of 50,000 new infections, we would expect 50,000*0.7% = 350 _new_ intensive care patients _every day_\n", "\n", "- with the exponential growth scenario (Option 1), we expect 0.7% * 2.5million = 17,500 new intensive care patients over the period of 28 days.\n", "\n", "[1] Source: [Lothar Wieler, Pressekonferenz zur Corona Lage, 12 November 2021](https://www.youtube.com/watch?v=7gK9YO3LUQI): for 50,000 new infections, expect at least 3000 hospital patients, 350 in need of intensive care, and 200 are expected to die." ] }, { "cell_type": "markdown", "id": "a0adfd2b", "metadata": {}, "source": [ "## Limits of growth" ] }, { "cell_type": "markdown", "id": "b0c80aed", "metadata": {}, "source": [ "The assumed growth assumptions are of course a huge simplification. Reasons that may slow down the growth include:\n", "\n", "- governements imposing new restrictions\n", "- people changing their behaviour voluntarily\n", "\n", "The big question is if and when these (or other factors) come into play to change our path. If there is no change to behaviour, no change to environment conditions, the virus, the accuracy and completeness of the testing, then reality may well follow a path between Option 1 and 2. \n", "\n", "[Austria](https://oscovida.github.io/html/Austria.html) seems to be on a similar trajectory, but is about 3 weeks ahead of Germany. (And has just [decided to have a lock down](https://www.theguardian.com/world/2021/nov/14/austria-orders-lockdown-those-not-fully-vaccinated-against-covid).)" ] } ], "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.1" } }, "nbformat": 4, "nbformat_minor": 5 }