Skip to content

Instantly share code, notes, and snippets.

@AlexRiina

AlexRiina/hospitalizations_and_viral_rna.csv

Last active May 27, 2020 05:07
Show Gist options
  • Save AlexRiina/7837ba1b74f8b9d8ea4c8a51a0cc6554 to your computer and use it in GitHub Desktop.
Save AlexRiina/7837ba1b74f8b9d8ea4c8a51a0cc6554 to your computer and use it in GitHub Desktop.
Download ZIP
sewer SARS-CoV-2 RNA vs hospital admissions
Raw
hospitalizations_and_viral_rna.csv
date viral_rna hospitalizations
19/03/2020 5 8
20/03/2020 5 10
21/03/2020 25 3
22/03/2020 0 8
23/03/2020 50 10
24/03/2020 30 12
25/03/2020 80 8
26/03/2020 80 10
27/03/2020 85 17
28/03/2020 75 14
29/03/2020 90 14
30/03/2020 95 26
31/03/2020 105 17
01/04/2020 45 12
02/04/2020 110 27
03/04/2020 70 21
04/04/2020 140 15
05/04/2020 115 17
06/04/2020 175 15
07/04/2020 105 21
08/04/2020 130 22
09/04/2020 125 21
10/04/2020 290 22
11/04/2020 120 25
12/04/2020 350 10
13/04/2020 185 29
14/04/2020 135 27
15/04/2020 205 20
16/04/2020 195 26
17/04/2020 80 25
18/04/2020 65 20
19/04/2020 100 5
20/04/2020 90 14
21/04/2020 120 12
22/04/2020 45 18
23/04/2020 55 23
24/04/2020 45 12
25/04/2020 90 17
26/04/2020 80 13
27/04/2020 85 15
28/04/2020 25 18
29/04/2020 55 11
30/04/2020 85 19
01/05/2020 90 19
Display the source blob
Display the rendered blob
Raw
sewer covid RNA vs hospitalizations.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas\n",
"import io\n",
"import seaborn\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.ticker import StrMethodFormatter\n",
"import matplotlib.dates as mdates\n",
"from datetime import date\n",
"\n",
"from statsmodels.nonparametric.smoothers_lowess import lowess\n",
"import statsmodels.formula.api as sm"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"df = pandas.read_csv(\"hospitalizations_and_viral_rna.csv\", \n",
" index_col=\"date\", \n",
" dayfirst=True, \n",
" parse_dates=[\"date\"])\n",
"df['viral_rna'] = df.viral_rna * 1000\n",
"\n",
"df['smoothed_viral_rna'] = lowess(df.viral_rna, (df.index - df.index.min()).days, return_sorted=False)\n",
"df['smoothed_hospitalizations'] = lowess(df.hospitalizations, (df.index - df.index.min()).days, return_sorted=False)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/alex/.pyenv/versions/3.7.4/envs/misc/lib/python3.7/site-packages/pandas/plotting/_matplotlib/converter.py:103: FutureWarning: Using an implicitly registered datetime converter for a matplotlib plotting method. The converter was registered by pandas on import. Future versions of pandas will require you to explicitly register matplotlib converters.\n",
"\n",
"To register the converters:\n",
"\t>>> from pandas.plotting import register_matplotlib_converters\n",
"\t>>> register_matplotlib_converters()\n",
" warnings.warn(msg, FutureWarning)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAEKCAYAAACfamUvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOyde3hU1bn/P28SICGQEC7hFpKAIBTBCyIXtYJSLN6KtdZKf1XbSu1NrW2fY7Xn+Z329LS/WntR7MWWorbaHtBWWz09WksV8Mb9GqGgCEm4hXANJIFAMu/vj7UnGYaZyUwyk5lJ3s/z7Gdmr73W3msmO/Pd77ve9S5RVQzDMAzDiExGsjtgGIZhGOmACaZhGIZhRIEJpmEYhmFEgQmmYRiGYUSBCaZhGIZhRIEJpmEYhmFEQcIFU0QyRWS9iPzN2x8uIitFZLuIPCsi3b3yHt7+du94aZjz3SEi73vbHQHlF4tImdf+MRERr7yviCz26i8WkYJEf2bDMAwjNNFqQirSERbm14B/Bez/CHhEVUcCR4A7vfI7gSNe+SNevTMQkb7Ad4DJwCTgOwEC+DjwBWCUt83yyh8AXlPVUcBr3r5hGIaRHKLVhJQjoYIpIkXAdcACb1+Aq4A/e1V+D9zovZ/t7eMdn+G3EgP4KLBYVQ+r6hFgMTBLRAYDeaq6Ql0mhqfDnDfweoZhGEYHEqMmpBxZCT7/o8D9QG9vvx9wVFUbvf3dwFDv/VBgF4CqNopIjVf/YMD5musEtR/qvQ8uBxioqvu891XAwFAdFZG7gLu83Yt79uwZ5Uc0DMMwAOrr6xVYF1A0X1XnB+zHogkpR8IEU0SuB6pVda2ITE/UdWJBVVVEQuYC9P6o8wFyc3O1rq6uQ/tmGIaR7ojICVWdGOZYymlCrCTSwrwM+JiIXAtkA3nAPKCPiGR5TxRFwB6v/h5gGLBbRLKAfOBQ0Dn3ANMD9ouApV55UVC5/7z7RWSwqu7zXLfV8fl4hmEYRgzEqgkpR8LGMFX1QVUtUtVS4FbgdVX9P8AS4Gav2h3Ai977l7x9vOOvexbhUBF5zSt/FbhaRAq8YJ+rgVc9l+sxEZni+cRvD3PewOsZhmEYHUQbNCHlSMY8zG8B3xCR7Tj/9RNe+RNAP6/8G7REsw4GGgFU9TDwX8Bqb/ueVwbwFdxA8nbgA+AVr/whYKaIvA98xNs3DMMwUoNwmpBySKov7yUidwOVqvpSR13TxjANwzBiR0TqVTU32f1IFCkvmMnABNMwDCN2OrtgWmo8wzAMw4gCE0zDMAzDiAITTMMwDMOIAhNMwzAMw4gCE0zDMAzDiAITTMMwDMOIAhNMwzAMw4gCE0zDMAzDiAITTMMwDMOIAhNMwzAMw4gCE0zDMAzDiAITTMMwDMOIAhNMwzAMw4gCE0zDMAzDiAITTMMwDMOIAhNMwzAMw4gCE0zDMAzDiIKECaaIZIvIKhHZKCKbReQ/vfLfichOEdngbRd65SIij4nIdhHZJCITwpx3lohs8+o9EFA+XERWeuXPikh3r7yHt7/dO16aqM9sGIZhhCdWXUg1EmlhNgBXqeoFwIXALBGZ4h37N1W90Ns2eGXXAKO87S7g8eATikgm8Euv7lhgjoiM9Q7/CHhEVUcCR4A7vfI7gSNe+SNePcPo2vh8UFUFFRXu1edLdo+MrkGsupBSJEww1VHr7XbzNo3QZDbwtNduBdBHRAYH1ZkEbFfVHap6ClgEzBYRAa4C/uzV+z1wY8B5f++9/zMww6tvGF0Tnw/KymDKFCgtda9lZSaaRsJpgy6kFAkdwxSRTBHZAFQDi1V1pXfoB57b9RER6eGVDQV2BTTf7ZUFEq5OP+CoqjaGaNvcxjte49U3jK5JdTXMnu2sS3Cvs2e7csNIMDHqQkqRUMFU1SZVvRAoAiaJyDjgQWAMcAnQF/hWIvsQLSJyl4isEZE1jY2NrTcwjHSloaFFLP1UVLhyw2gfWf7fUW+7K7hCOulCMB0SJauqR4ElwCxV3eeZ5Q3AUzg3K8AeYFhAsyKvLJBwdQ7hXLhZIdo2t/GO53v1g/s4X1UnqurErKys4MOG0Xno0QNKSs4sKylx5YbRPhr9v6PeNj9cxSh1IaVIZJTsABHp473PAWYCW/3jkt444o3Au16Tl4DbvWjZKUCNqu7z6m716qwGRnkRsd2BW4GXVFVxX/zNXr07gBcDznuH9/5m4HWvvmF0TQoL4cUXW0SzpMTtFxYmt19Gp6cNupBSJNKUGgz83otszQCeU9W/icjrIjIAEGAD8CWv/svAtcB2oB74HICI9PfqoqqNInI38CqQCTypqpu99t8CFonI94H1wBNe+RPAMyKyHTiME1nD6LpkZMD48bBihXPD9ujhxDLDpmUbCSdWXUgpJNWNLRG5Hhihqo911DVzc3O1rq6uoy5nGIbRKRCRelXNTXY/EkXKD9ap6t+S3QfDMAzDMB+MYRiGYUSBCaZhGIZhRIEJpmEYhmFEgQmmYRiGYUSBCaZhGIZhRIEJpmEYhmFEgQmmYRiGYUSBCaZhGIZhRIEJpmEYhmFEgQmmYRiGYUSBCaZhGIZhRIEJpmEYhmFEgQmmYRiGYUSBCaZhGIZhRIEJpmEYhmFEgQmmYRiGYUSBCaZhGIZhRIEJpmEYhmFEQcIEU0SyRWSViGwUkc0i8p9e+XARWSki20XkWRHp7pX38Pa3e8dLw5z3DhF539vuCCi/WETKvPaPiYh45X1FZLFXf7GIFCTqMxuGYRjhiVUXUo1EWpgNwFWqegFwITBLRKYAPwIeUdWRwBHgTq/+ncARr/wRr94ZiEhf4DvAZGAS8J0AAXwc+AIwyttmeeUPAK+p6ijgNW/fMAzD6Hhi1YWUImGCqY5ab7ebtylwFfBnr/z3wI3e+9nePt7xGX4rMYCPAotV9bCqHgEW477wwUCeqq5QVQWeDnPewOsZhmEYHUgbdCGlSOgYpohkisgGoBonbh8AR1W10auyGxjqvR8K7ALwjtcA/YJO2VwnqP1Q731wOcBAVd3nva8CBobp610iskZE1jQ2NoaqYhiGYUQmy/876m13BVeIURdSiqxEnlxVm4ALRaQP8BdgTCKvF0V/VEQ0zLH5wHyA3NzckHUMwzCMiDSq6sRIFVJNF2KhQ6JkVfUosASYCvQREb9QFwF7vPd7gGEA3vF84FDQqZrrBLXf470PLgfY77ls8V6r4/CRDCO5+HxQVQUVFe7V54vr6VWV2tpajh49Sm1tLW6kwzDiR5S6kFIkMkp2gPcEgYjkADOBf+G+oJu9ancAL3rvX/L28Y6/7lmEQ0XkNa/8VeBqESnwgn2uBl71XK7HRGSKN+55e5jzBl7PMNITnw/KymDKFCgtda9lZXETTVWlurqaBQsWMG/ePBYsWEB1dbWJptFu2qALKYUk6p9ARM7HDd5m4oT5OVX9noiMABYBfYH1wGdUtUFEsoFngIuAw8CtqrpDRCYCP1DVj3rn/Tzwbe8yP1DVp7zyicDvgBzgFeAeT3D7Ac8BxUAFcIuqHo7U99zcXK2rq4vXV2EY8aWqyolkRUVLWUkJrFgBgwa1+/S1tbUsWLCAmpqa5rL8/Hzmzp1Lr1692n1+o/MiIvWqmhvheEy60BF9joWEjWGq6iac+AWX78BNCQkuPwl8MsSppgC/DKj3JPBkiPZrgHEhyg8BM2Lpu2GkNA0NZ4oluP2G+Py+NDY2niGWADU1NVgwnNFeYtWFRCIixwDBRekC9ATq/YdVtXdwm4QG/cQDVf1FsvtgGClFjx7Oogy2MHv0iMvps7KyyM/PP8vCzMpK+Z8Lw4gaVc0L3BeRdao6IVIbS41nGOlGYSG8+KITSXCvL77oyuNAbm4uc+bMIT8/H3BiOWfOHHJzw3raDKNLkLAxzHTGxjCNlMfng+pq54bt0cOJZUb8nn9Vlbq6OhobG8nKyiI3N5ez84gYxpm0NoaZyojIelU9y10ciPlYDCMdyciIS4BPOETEAnyMTk1gLnKPgsAyVf190HETTMMwDKNLcnHQ/ksBZUJLStVmzCUbAnPJGoZhxE46u2SjwYJ+DMMwjC6HiBSLyF9FpFpEDojISyJSEqmNCaaR1lgKt+Rg37vRCXgKt0LKYGAQLsHNWXP8AzGXbAjMJZse+FO4LVy4kJqamubpD4WFhRbRmUDsezfCkU4uWRHZoKoXtlYWiFmYRtpSV1fX/KMNLhvNwoULsYedxGLfu9FJOCAinxWRLG/7HHAgUgMTTCNtsRRuycG+d6OT8DngBmCvt33MKwuLTSsx0hZL4ZYc7Hs3OgOquhv4RCxtzMI0Ek6iAkQshVtyaO/3bgFDRiogIs97q6QgIo+LyCYRCbUASEsbu1nPxoJ+4keiA0QshVtyaOv3bgFDnZs0C/rZpKrni8glwA9w63D+Q1XHh2tjFqaRUBIdIOJP4danTx969eplP7odRFu/dwsYMlKQ63Hrcu4DIg7E26CDkVAsQMQIxO4HI4X4p4isBAqBi0UkD6iJ1KBNFqaIxDRQanRd/AEigViASNfF7gcjVVDVbwBfACao6mFVPaaq0yO1adMYpohUqmpx27qZ+tgYZvywMSsjELsfOjfJGsMUkUHAncB1wGic9/QI8DbwR1V9OaBu8ColZxBqlZLmtm0UzF2qOqyVOsOAp4GBgALzVXWeiHwXp+r+CaLf9n8YEXkQ96GbgHtV9dUQ570Y+B2QA7wMfE1VVUT6As8CpUA5cIuqHhH3XzgPuBaoBz6rqusi9d0EM75YYE7Hk8rfeWt9S+W+G5FpTTDbogtRXPNmr+1zwBJV3eGV5+JWH5kNXKiqM7zyx7ymRcAkYDFudZKZwGpVvTHstRJlYYrIYGCwqq4Tkd7AWuBG4BagVlV/ElR/LLDQ+wBDgH8C56pqU1C9VcC9wEqcYD6mqq+IyMPAYVV9SEQeAApU9Vsici1wD04wJwPzVHVypL6bYBrpTDpbcencdyMqwYxJF6K8Zi9VrY21joi8Dtykqke9/T7AX1T1ynDnCTuGKSJl3ryU4K0M93QQEVXd57fkVPU48C9gaIQms4FFqtqgqjuB7TjxDOzTYCBPVVeoU/qncV+2v73flP59UPnT6lgB9PHOYxidknSORE3nvhut0wZdiOacEcUyQp2hQEPA/kmcsRaWSCPt17fWiWgRkVLgIpxVeBlwt4jcDqwBvqmqR3CdXxHQbDdnf5FDvfJQdQZ6YcEAVbSI+lBgV4g2+wLKEJG7gLsAunfvHtPnM4xUIp0jUdO57wYAWSKyJmB/vqrOD1UxSl2IGhE5hnOtKtATOOG9B8hV1cygJs8AK0Tkr97+TcAfI10jrIWpqhWBG24A9XjAFu2H6AU8D9ynqseAx4FzgAtxovXTaM8VLZ71GZOvWVXnq+pEVZ1oEXtGOpPOkajp3HcDgEb/76i3hRPLuOuCquapam9VzQM2+d97+xtD1P8+8FngIHAIuENVvxfpGq1OKxGRL4pIFbAJ529ei3sCaBUR6Yb7Uv6oqi94ndyvqk2q6gN+S4vbdQ8QGEhU5JUFsscrD1Vnv9/V6r1Wx3Bew+g0pHPKwHTuuxEdMepCe64T+JTVLVQdVV2vqj9X1cdaCwaFKIJ+ROR9YKqqHoyxs4IbSzysqvcFlA/2u05F5OvAZFW9VUTOA/6blqCf14BRqtokIq8Bt6vqnhBBPz9X1ZdF5MfAoYCgn76qer+IXAfcTUvQz2OqGvGPYUE/RrqTzpGm6dz3rk4UQT8x6UI7+jEPKAFeAT4MdFPVTwXVCXThNheram8ReUtVLz/rvFEI5t9xkUT1MXb4cuBNoAzwecXfBubgzG7FTf/4YsAX9e/A53Hpie7zol8zgJ3AGFU9ISITaZlW8gpwjzetpB8urLgYqMBNKzns/YF+AczCTSv5nKpGtJBNMA3DMGInCsGMWRfa2I8s3BTF8bgA0l+r6sm2nq/5vFEI5kXAUziLrjmiSFXvbe/Fo0FExgGfV5eVoUMwwTQMw4idZCZfF5HuOMOoRlWXRVE/Czgf6B1Q/BDwILDTi905g2hG0n8DvM6ZTwQdhqq+C3SYWBqGYRhpyV+B00CBN8fyUeApVf14mPqv4uJ4AoNYxwDfxA0Ptkkwu3WkdWcYhmEYbaBEVc8TkR7ASlX9rogURajfV1UvCiwQkXWqekO4BtEkX39FRO4SkcEi0te/RfkBDMMwugxqi2Mnk20iMkZVGwBEJBvIjlD/dyHKwuaRhegszDne64MBZQqMiKKtYRhGl8DS+iWdPsB6EVmBi5BdjcsjHo6/icjPcEt6PQKcwrl1wxI26EdEhqjq3rb0Ot2xoB/DMGKltraWBQsWnJGpKD8/n7lz59KrV68k9qzjSHLQzxUBuyeB9yNlCxKRjTgrczAwCJgL/FNVrwjXJpKFucBzvS4F/g68paqWn8owDCMEltYvuajqG8FlIvJDVX0wVH3Ap6qPePU2qeopEekZ6RqRUuNdC0zHCebHcTn3XvDGMzvtWpiGYRhtwdL6JRcR+aGIHBKR4yJyzEtM8G/efijR/LuIfE5EMoEmERnZ6jViGZQWkeHANbi5LoNay5iTrphL1ugI2pPRxrLhpB42hpl0l+wW4PxAT6gX9TohTP1jQC5u/eVTuJVT7vFWtQpJTI8+qrrTW95rHHBzLG0Nw2ihPT+u9sOcmogIhYWFzJ071x5kksPaEMOGm8NV9pKyx0RUFqaX7efTwCdxaeqeV9VfxHqxdMEsTCPRtCdAxIJLjFQlmRZmrIhISaTjMWX6EZFzcVNK5uCWP3kWJ7BhV6M2DCM62hMgYsElhnE2bUim/j8h6je3w+WhPYNILtmtuCS516vqdq9DX4+h/4ZhhMEfIBJsJUYTINKetobRWYnkYg218oiqnh/rNSJl+rkJt5DnEhH5rYjMwKmuYRjtpD3rPtqakYYRGhGZKSI/EZEfi8jMuJ8/itVKcoHZONfsVcDTwF9U9R/x7kyqYGOYnYdkRpO2du32HLcoWSMVSXKU7D24WJsncZnpFuOSF/wkbteIcVpJAS7w51ZVvSpenUg1TDA7B8mMJm3vtS0S1khHkiyYm4Cpqlrnn04iIqviOf0xmuTrzajqEVWd35nF0ug81NXVNQsOuMCYhQsX0hEPQ+29djL73m58PqiqgooK9+rr8FUBjS6Kqvr/QUTck2X3cHVFZLuI/MBbczkqIkXJHqclesj/SKtem+6qahEGRkqTzGjS9l472vaqSlNTE01NTfh8vub3/n1VJSsr66wtYVaqzwdlZTB7thPMkhJ48UUYPx4yYno+N4xYOR6QAz0XeAn4S4T61+Ky2K0RkfeBhcAiVd0RrkFY0VPVwFWoEZFewFeBL7bSCcNICZIZTRqPa/fq1Yva2trm/e7du/Pyyy/T0NBAXV0ddXV1nDx5sk39y8zMJCsri+7du5OTk0PPnj3Jyck5a+vduzd5eXnk5eWRnZ3dutBWV7eIJbjX2bNhxQoYNKhNfTWMKPkMbgFpgB/ixi/fCldZVd8DfiQiXwRuBD4FvCgidTjhfDS4TTRBP32A+4DbcatQP6Kqh1rruYgMwwUIDcRZpvNVdZ6X0P1ZoBQoB25R1SOe+TwPp/r1wGdVdV2I887y6mUCC1T1Ia98OLAI6AesBW7zkun28PpxMXAI+JSqlkfqu41hdg7SYQzz1KlTHDhwgIMHD57xeuTIkbPWUszOzqZ3797k5uY2bzk5OWRmZjZvGRkZZ+yDs1bDbQ0NDZw8eZL6+npOnDjBiRMnqK+vxxfCjZqVldUsnnl5efTu3ZuCggIKCgro27cveXl5ZOzahd56K3Xf+x6NAwaQdeAAuf/xH8iiRc7ajNN3awFPqUlrY5ix6kKM1w55g4VKQBDUboeqjvA06HLg+8AkVc05q26E5b36A9/Eqe6TwM9VtSZk5dDtBwODVXWdiPTGidiNwGeBw6r6kIg8ABSo6rdE5FrgHpxgTgbmqerkoHNmAu8BM4HduPXO5qjqFhF5DnhBVReJyK+Bjar6uIh8BZdf8EsicivwcVX9VKS+m2B2HlItSra+vp7KykoqKiqorKykqqqqWRgzMjLo168fAwYMoH///vTv35/c3Fyys7Pp1asXvXv37pC+qyqnT5+mvr6e2tpajh071rwdP378jP1AYc3IyKBP7970zM6m+sgRTp06Rc+ePfnYtGmMLC0ls7AwLn2zYKjUJQrBjEkXYrz2JloSEfQAhgPbVXVsK+1248T6Jpy+/DdOS46fVTeCYNYBB4CngLMaqurPYvwwLwK/8LbpqrrP+/KWqupoEfmN936hV3+bv17AOaYC31XVj3r7/gz0D3l9HaSqjYH1RORV7/1yEckCqoABGsG0NsE04sWxY8fYuXMnlZWVVFZWcvDgQcBZa0OHDqW4uJjBgwczYMAACgoKmq3CdMDn83H8+HEOHz7MkSNHOHLkCNXV1XzwwQc0NTWdUTcjI4OCggL69+9Pv379mh8MCgsL6dGjR9TXtLSAqU2sUbKt6UI7+zIel0z9rjDHf4Ibw9yLG7/8k6oeiHTOSAMqP6Yl6Kd3hHqtIiKlwEXASmBggAhW4UxzgKHAroBmu72yfQFloepMxrlhjwYk3vW3PaONJ6Y1Xv2DQX28C7gL3FiRYbQFVWX//v1s27aNbdu2sW+fu32zs7MZNmwYF1xwASUlJQwePDjtM/NkZGSQn59Pfn4+w4cPB+Do0aPMm3f2IvcTJkygrq6OgwcPsn379jMEtU+fPgwaNIjCwkIGDhzIwIEDKSgoICNEkJClBUx5skRkTcD+fFWdH6pilLrQZlS1TEQujVBlP06kd0WocwaRgn6+G+6Yl8wgKrxgoeeB+1T1WNDEbBWR6CeCJhDvjzofnIWZ5O4YaURTUxMVFRXNIun/QS8qKmLGjBmMHDmSgQMHdgmXYbhgp2nTpjkL0OfDt2kTNZ/5DAdOn6Zq9GiqP/Yx9h84wLZt25rd0926daOoqIiRI0cycuRIBgwYgIhYWsDUp1FVJ7ZWKRG6ICLfDNjNxMWt7I7Q5DHgKyJyhbf/BvC4qoaNpIt4l4nIUGAw4F+NuhAXAPRZYEgUH6Ab7kv5o6q+4BXvF5HBAaZ3tVe+BxgW0LzIKwskXJ1DQB8RyfKszMC2/ja7PZdsvlffMNqMz+ejvLycTZs2sXXrVhoaGsjKymLEiBFcccUVnHvuuV3SRehP2xc8xtictq+6mowbb6SgooIC4Nz33oNNm2DFCk7368eBAwfYv38/VVVVlJeXs3jxYhYvXkxeXh4jR47knHPO4ROf+ATPP/986PMbKU+MuhALgTdBI/Cid51w/A6oxQWRgouyfQqX1S4kkeZh3gf8O7Ad6CEivwJ+REvEaUS8iKMngH8FjXe+BNyBG3e8A/eh/OV3i8ginJu1xm+ii8hWVR2DC/IZ5UXE7gFuBT7tPZEswa3RuSjEee8AlnvHX480fmkY4fC7Wzdt2kRZWRm1tbX06NGDMWPGMGbMGM455xy6deuW7G4mFWltTciGhpYpJ34qKqChgW7dujFkyBCGDGl5Fj927Bjbt29n+/btbN68mXXr1pGRkcGQIUP40Ic+RElJCX379u0S1nsyiTaN49ixY8+KLA2kDboQSx+/F+J6nyC8aJ6vqucF7C8VkbDrZ0LkoJ8twOWqelhEinHRQ5ep6tpoOi8il+NWOykD/KF038b5q58DioEKXPjwYe+L/AUwCzet5HOqusaL1n3bPwDsRdM+ijO5n1TVH3jlI3Bi2RdYD3xGVRtEJBt4BucrP4xL6xd2YipY0I9xJjU1NZSVlbFp0yYOHDhARkYGo0aNYvz48Zx77rldXiRjoqoKpkw5UzRLSqKap9nU1MTu3bubBbSqqgpwbuCioiJKSkooLS2lqKjIXLRxpLXI5MDjDz/8MHv37g379BKrLsTSTxG5Cef9DIy5mQisAX6nqr8Pqv9nXEDou97+eOA/VfWmsNeIIJjrVHVCwP5GVb0glg8QD0TkemCEqj7WUdc0wTSamprYtm0ba9euZccO93w1bNgwzj//fMaOHUvPnj2T3MME4/O5JAQNDdCjBxQWxidTTxwzAZ04cYLKykrKy8upqKhoDrDKzMxsFtCSkhKKiooskK8dtBaZHHj8N7/5TUTBTCQishX4EnDMK1LcFJHPAHtUtSqo/nKcoG7y6l6AE9eTABpi7edIj2FFIhIoUoMD91X13lg/UFtQ1b91xHUMA+DQoUOsW7eOjRs3UldXR15eHtOmTeOCCy6goKAg2d3rGBKZ3i4jw51nxYp2i3FOTg6jR49m9Gg3+yBYQN98803eeOONZhducXExJSUlFBcXk52d3b7P0YVoLTI51PEkUa+qSwMLROREBK/oV2O9QCTB/Leg/ahcsYaRbjQ2NrJ161bWrl1LeXk5IsLo0aOZMGEC55xzTsjpDZ2aRKe3y8hISJq8YAFtaGhg165dVFRUUFFRwYoVK3jnnXcAGDRoULMFWlxcbEFDEWgtMjnU8SRxjYjkqlutJBMX4Dk1XOVQmeRaI6blvboK5pLtGtTU1LB69WrWr19PfX09ffr0YcKECVx44YX07t2uqcfpTUUFlJaeXV5eHrf0dsng9OnT7Nmzp1lAd+3a1WwlDRgwgOLiYkpLSykpKQn790/ntHxt7Xs8xzATiYisBK4HjgKrgBzgr6r6QJj6x2jJDNQNlx2oToPyqJ/RxgTzbEwwOy+qSkVFBatWrWLr1q0AjB49mokTJzJixIi0+fFLKO0IzEknmpqa2Lt3b7OAVlZWcurUKQD69u3bHERUUlJCfn5+Wqfla2/fo42SnTx5sm7evDkpLhl/nI2Xb/xmVZ0rIpuDImEjtb8a+Iiq3h+2jgnm2Zhgdj5OnTpFWVkZq1atorq6mpycHCZMmMDEiRPp06dPsruXWnTRJbp8Ph9VVVXNAlpRUdG8GkyfPn0YOnQoO3bs4MSJE81t0iUtX0elFIw1NV48EZGNwBXAz4C/qepfRGSDql4YwzneVdWw62Na7LXRqampqWHVqlWsW7eOkydPMmjQID72sY8xbty49k8HSVQkabKJY2BOOuEPDhoyZAhTp07F5/NRXV3dHAdy9U0AACAASURBVET0wQcfnLWcWk1NDYcPH055wewiKQV/jMsbsAn4m4jkAf8IV1lEvhOwmwGMx+WVDUurgikiD+OWOzkB/B04H/i6qv6htbaGkSz27dvH8uXL2bx5M6rK2LFjmTRpEsOGDYuP+6yzW2EJCsxJJzIyMhg0aBCDBg1iypQpHD9+nPnz55+xRinAU089RZ8+fSgtLW3e8vPzk9Tr0HT2lIIi8n+Avwfp0mkgrHuVMxcV6QH0xCXrCX+d1lyyfpNWRD6OG1D9BvBGMuZkdhRdzSWb7ECGeF1fVXn//fdZvnw55eXldO/enQkTJjBlypQ2/YBF7FcXGeczWggeB8zLy+OjH/0ox44da3bh+t21BQUFzeI5fPjwqILIEvl/2FHjr8lyyYrI/cDVuECf14BXgJWqevbCruHPkQG8Fmr+ZXOdKATzXVUdJyILgD+r6t+TlcSgo+hKgpnsQIZ4XL+xsZFNmzaxfPlyDh48SF5eHpMnT2bChAltnm/Xar86aSSpEZlIouZPnVheXt7sxvW7cPv163eGBRrswu2I/8OOeDBO5himd/1c4CO4jHGTcRnq/o6zPqsitBOcS/YFVR0Ztl4UgvkQboHPE8AkoA9uQHVyxIZpTFcSzGSvL9ie6zc1NbF+/XreeOMNjh8/zqBBg5g6dSrnnXdeu9eVbLVfZmEareDz+di/fz87d+5sFlB/FO6AAQPOEFCfz9cp1vlMtmAGIyJjgGuAWeqtoxxwLHBaSSMuJd//jZQsp1UHtqo+4I1j1qhqk7iFpWe34zMYKUSygwHacn2fz0dZWRlLly7l6NGjDBs2jBtvvJHhw4fH7Ym51X4VFroxy+AxzMLCuFzfSH8yMjIYPHgwgwcP5tJLL8Xn87Fv375mAd2wYQOrV68GoH///l0hKCehiEgo184J4AVvOwNVzYv1GtEE/dwe8D7w0NOxXsxIPZIdDBDL9VWVLVu2sHTpUg4ePMigQYP49Kc/zciRI+PuWmq1X100ktRoOxkZGQwdOpShQ4dy+eWX09TUxJ49eygvL+eDDz44q3737t3ZtWsXPXv2tFy40bEF2IGzGIMZiQvqaUZEpoU6iaouE5GLQ6XUi8Yl+/OA3WxgBrBOVW+O3Pf0pSu5ZNNhDNMfzLNkyRKqqqro378/V155JR/60IcS1sdkfy9G10JV2bt3LwsXLqSuro7MzExUFZ/P1yy0/gCioqKilF0hJ8nzMNeqasilJ0MdE5GXQlVV1RtE5FFVvS/UwVg71QdYpKqzYmqYRnQlwYTUjZJVVXbs2MHSpUvZvXs3BQUFTJs2jfHjx3dIftdkfy9G1yL4fuvWrRu7d+9uduHu3bsXVW1ejWX48OHNy5m1d8w+XiRZMLeo6thYj8V0jTYIZjfgXf/6lJ2RriaYqUh5eTlLliyhsrKSvLw8PvzhD3PRRRelzA+DYXQ0DQ0NVFRUUF5ezs6dO89YD9SfB7e0tJQhQ4Yk7f8kyYK5ArfecXlQeSnOyJsSVP4Y8Kiq7hCRbwGXAT9V1WVhrxGFS/Z/aPEJZwBjgefCJbTtDJhgJo9du3axZMkSdu7cSa9evfjwhz/MhAkTOs0Ea8OIFydOnKCioqLZAq2urgbc2GeggA4ePLjDVtxJsmDeBtyHW2lrHU63JuBS5T0aYgHpMlUdLyLnAU8AXwd+paoXhb1GFIIZODDaCFSo6u42fJ60wQSz49mzZw9Lly5l+/bt5ObmctlllzFx4sSUHasxjFSjrq6ueQ5oeXk5Bw8eBKBHjx5nCOigQYMSJqDJnlbiBaneC3zIK9oGzAsWS6+uP1n7/cAJVf25iKxT1Qlhzx9JML01xf4ZKfNBZ8QEs+PYu3cvy5Yt47333iMnJ4fLLruMSy65xKICDaOd1NbWniGghw4dApyA+tcCjbeAJlswY0FE/gh0By7BuWOPAG+GCxyC6CzM14CbVDWm1UFF5ElcKr1qf/Z3Efku8AXggFft26r6snfsQeBOoAm4V1VfDXHOi4Hf4dIfvQx8TVVVRPoCzwKlQDlwi6oe8bI3zAOuBeqBz0azaKgJZuLZs2cPy5Yt4/333ycnJ4cpU6YwefJkevTo0aH9aE9gjwUFGalEa/fj8ePHz8hCFCigxcXFzQLaHhdua4IZqy7EeO0puNyxNcD/xa2L+SFVXR2mfjdcUoOtqvqeZyBmq2rYH/9oBPNF4CJgMdB8IlW9t5V2VwC1wNNBX0ytqv4kqO5YYCEuk9AQ4J/AuaraFFRvFc7cXokTzMdU9RUvscJhVX1IRB4AClT1WyJyLXAPTjAn40zzVjMUmWAmjmChnDp1KpMmTepwoYT2TR2xaSdGKtGW+/H48ePNQUQVFRXNLtzu3bszbNiwZit0yJAhUccQRCGYUetCrIjINuABYChwFXAz8JaqXhqm/u3AUlWtFJEP4yzN/46YQi8KwbwjVHkon3CItqW4NHqtCeaD3jl/6O2/CnxXVZcH1BkMLFHVMd7+HGC6qn7R+6Kmq+o+r95SVR0tIr/x3i/02jTXi9RvE8z4k0pC6ac9afmSnVLQMAKJx/3od+H6F9P2BxFlZWVRVFRESUkJxcXFFBUVhR0yicYlG60uxErgXMuAgJ5IczPLgAuAQuAN4ClcCr2QCQ0gutR4rQpjjNztKfsa4JuqegT3RLAioM5uryyQoV55qDoDA0SwChgY0GZXiDZnCaaI3AXcBdj4WRzZvXs3y5YtY/v27eTk5HDVVVclXSj9tCctYLJTChpGIPG4H3v16sW4ceMYN86tn1xfX09lZWXzSixvvPEGqtqc8q+4uLh569mzOYlOloisCTjtfFWdH8XlQ+lCrPyvJ75PASoiM3Cp8cJxWlV9InId8AdV/aGIfDLSBcIKpog8p6q3eCp8lhmqqudH9RHO5HHgv7zz/RfwU+DzbThPWLwxzdgml7p284H54CzMePapK1JZWcmyZcvYsWMHOTk5zJgxg0suuSQlhNJPe9ICJjuloGEEkoj7sWfPnowZM4YxY8YAcPLkSXbt2kVlZSWVlZWsWrWK5cudE7B///4UFxcDNKrqxBgvFS9d+Iz3ejtwEvgy8LkI9Y+LyFeBucBnvZiXiF9YpINf816vj66vraOq+/3vReS3gD8r/B5gWEDVIq8skD1eeag6+0VkcIBLtjqG8xpxpLy8nGXLllFeXk5ubi4zZ85k4sSJKWm15+bmMmfOnLPGfXJzWw/ya09bw4g3HXE/ZmdnM2rUKEaNGgU4q3bv3r3NLtzNmze36bwRdCHW84yIsYl/3uZPVXWjtzTYlyM1iGYM8xu4LAl7Y+xMKF/1YL/rVES+DkxW1Vu9iaP/TUvQz2vAKHWro7wG3K6qe0IE/fxcVV8WkR8DhwKCfvqq6v2eqX03LUE/j6nqpNb6bWOYZ9Ja9J3P5+Nf//oXy5cvZ8+ePfTq1YvLLruMiy++OOXnUWpTE3U1NTT6fGRlZJCbn49EmSXFomTTj878N0v2Z/P5fGRmZrZlDDOkLkR7XRG5TVWfCSrLwQX9zI00Jhkr0djrvYHFInIYN3XjT4FPBOEQkYXAdKC/iOwGvgNMF5ELcaZ3OfBFAFXdLCLP4bLNNwJf9cQyA5dl/rB32q/QMq3kFW8DeAh4TkTuxK1pdotX/jJOLLfjppVEMs+NEESKvgN4//33ef3119m/390S2dnZ3HrrrQwZMiT1f4h8PuTdd+kVvETX+PFRrToiIhbgk0Z09sjmZN+P0UxFiUUXYuD7IvKaqu4VkQm46Ykzgf/FzZII15cdtKyH2VysqsNF5G+qepZ3NepcsiJyPvAp4BPAblX9SNQfp42IyDjg86r6jURfKxCzMFsIFX2Xl5fHlVdeyerVq9m7d29zonQ/aRMtaotAtx2fD6qr02ppM4tsTjzJSFwgIjcB/w8nfHk4L+Rfg6clhmjXN1S5qh4WkTxVPRZ8LJY7vBoXgXoIF4abcFT13Y4WS+NMQkXfHTt2jBdffJETJ04wY8YMgh+60iZatKHhTLEEt9/QkJz+pAs+H5SVuYeN0lL3WlbmylMYi2zunKjqC950w3uBZcCPcFbnqFbaHQZKcMFCnwFKvTJCiSVEIZgi8hURWYobV+wHfKGNEbJGGpKVlUVe3pkLk2dkZDBr1izuvvtuLrzwQvLz8884njbRoj16OIsykJISV26Ep7oa/G5scK+zZ7vyFMYfSRpI2tyrRquo6mJV/TRwMW464R9F5M1w9UXkPtwUlL7e9pQXsxOWaIJ+fgg8q6obYux/2mIuWYfP56OsrIylS5dy9OhRAHJycrjtttsYNGhQsys2bceF/JZSG8cw43L9NHNrAu67Ki09u7y8/OwHkBQire9Vkh/UEw3JcMlGQkTGq2pZmGNlwCWqetLbzwZWq+r4sOeLdgyzK9HVBbOpqYlNmzbx5ptvcuTIEQoLC5k0aRIjRoygW7duZ/2jpsM/cliSJVrJFuv2kMZjv+l6r6aL2CdpDHMiUKmqYV0cInKdqv5vUNkmYKo/d6w3rWR5JA+qCWYIuqpgNjU1sWHDBt566y2OHj3KoEGDmDZtGqNHj06pf8pOQRqLTlqLfZqSLgFLSRLM83HJDk4CS4D3vPeFOPfsdNw0lu8EtbsHl93tL17RTcBvVXVeuGuZ896gsbGR9evX89Zbb3Hs2DGGDBnCNddcw6hRo0woE0U0AUep6rLNyHDiuGJF6vWtk2IBS+FR1U3ATG8+/3W4aYU5wF7gTeDf/W7XoHY/F5E3gA97Rf9HVTdGupYJZhemsbGRtWvX8vbbb3P8+HGKioq44YYbOOecc0woE40/4CjYwvQHHKW6FZeRkfqWcCfCUjG2jqpuBqJONyQiJbglwP4nsExVK8K2MZfs2XR2l+zp06ebhbK2tpbi4mKmTZvG8OHDTSg7itYEMZ1dtkbcsTHM+OONYfoTF/QAhgPbVXVs2DYmmGeTioIZj2CFU6dOsWbNGt555x3q6uooLS1l2rRplJSUpNQ/XSDpGqQRFZFcru2NRE1Vd26K055Uia2eu533cjr8L6STYAYjIuOBe1T1rnB1zJ5PA9r7dHnq1ClWr17NO++8Q319PcOHD28WylQmXZ6q20wkt2ZrLttIpLo7N0XRpiaqq6pY+Kc/tdxvn/wkhYMGtVs043EvJzv1XWdHVctEJORi037MwgxBqlmYbY2QCxbKESNGMG3aNP8yPClP0iMD22ulnToFlZWwfz8cOwY1Ne41+P2pU+5aqu7V54OmJjh6FFauhBMnIDcXrr4ahgyBnBzIznZ96tkT+vSBgoKW7fRpuOkm2BWwFKy5c1ul9vBhFjz99Nn32+2306tvyCxq0Z872fdyB5FOFqaIfDNgNxMXUZuvqrPCtTELMw2INULu1KlTzWvV1dfXc8455zBt2jSGDRsWsn6qktTIwGittEOH4N13YedOt5WXt7zfs8eJYChEIC/PbT16uP2MjDM3ERgxwgn2qVMtUaknTzoRjeVht6ICvvENOPdcJ7pDh7rXYcOgXz93rS5Oo88X+n6LQ8o/i3JNSQKFvRF4EXg+UgMTzDQg2gi5hoaGZovyxIkTjBw5kmnTplFUVBR8yrQgqZGBodK/XX89/PSnTgzXrHFbeXlLGxEnRMOHw1VXudfSUhg8GPLznTj6X3v1ap9IqUJjI9TVOUv0yJGWrbwcvv99V+6ne3f45z9h0aKzhbZ3b9fXESNaXgPfd5FUgVkZGaHvtzi4sS3KNSX5f6oa0xOLuWRDkGou2dbGP/wWZWcRSj9RjfskKrilvNwJRjhGjICJE912wQVuv7jYCVOyiWQdNzU5F/Hevc4CrqyEHTvctnOnez1xouVcGRnuexgzxm2jR7e879+/U1mmqT6GmUjiFVCUZi7ZN3FJ1ytxS4FdAPxIVR8L28YE82xSTTAh9A19+vTps4Ry+vTpDB06NNndjRsR/5HjHdxSWQmvvw6vveassaqqM4/36QO//jXMnAntHNNKOG19kFB17XbsgA8+gPfeg23bYOtW9/5kwPzvvn3hvPPc9z1uXMtrnz6J+1wJJpWjZBNFPMU8zQRzo6peICLTgLuBTwMbVPW8sG1MMM8mFQUzkOBgns5iUcZMe+cq1tY6cXzlFfe6fbsrHzAArrzSjfc9+aSzxizS1IlwZaUTz23bYMsW2LzZPbQcC1gNadgwJ5znnw8TJrhtxIiu+72lOPEMSEozwdyEC/R5CFinqn8UkQ2qemG4NuZATyNOnz7dLJR1dXVpG8xzFm21htqynuUHH8D//q/bli51wTS9e8P06XD33W7s8bzz3PV9PvjqV9vu7u1scyEzMtyYbGkpzAoIJFR1EbnvvuvEs6zMvf/nP13ELrjv+KKLWgR0wgTn3rUxvKTThQOSngF2AseA/xCRPGBLpAZ2t6YB/sw8b731FnV1dYwYMYLp06env1BC+9yq0cxVbGyEN9+E//kfJ5LvvefKx4xxAnnddXD55aHHHtuT/q0rzYUUceO3xcVw7bUt5adOOQt03bqW7Te/aRkj7dnTjQFPnuy2KVNc0JTRoXTVgCRV/bGI/AaoVVV/KPSnI7VJmEtWRJ4ErgeqVXWcV9YXeBYoBcqBW1T1iDhH+TzgWqAe+KyqrgtxzllevUxggao+5JUPBxbhFrheC9ymqqdEpAfwNM7sPgR8SlXLW+t7OJdsIscgQp27qamJdevW8eabb1JbW8vw4cOZPn16m+ZRpur4SbvcquFE6ZxznHXzl7/A3/4Ghw87EZ0+3Qnkddc5F2EiSfHUdkm7HxobnTt37VoXZbxyJaxf32KJDh3aIqAf/jBccklsVmhns+o7gI4cw4xFF9r4cRJKIgXzCqAWeDrgi3kYOKyqD4nIA0CBqn5LRK4F7sEJ5mRgnqpODjpfJm7ZlpnAbmA1MEdVt4jIc8ALqrpIRH4NbFTVx0XkK8D5qvolEbkV+Liqfqq1vocSzERGuQWfOy8vjwsuuICNGzdy7NgxiouLufLKKykNlSqtDedPqQi9eKWAq6qCN95wY5GLFzsrpqAAbrgBbrzRTfrP7cChlRReZDnl7oeTJ2HDBiee/m3HDncsP9+NJ8+cCR/5CIwaFT4ytytZ9XGmo6JkY9GFNn6UhJLQoB8RKcWtQ+b/YrYB01V1n4gMBpaq6mjPLF6qqguD6wWcayrwXVX9qLf/oHfoIeAAMEhVGwPricir3vvlIpIFVAEDtJUPHUowE5mpI9S5AQYPHsxHPvKRdidFT+ksI+2xxE6dgpdeggULnEXZ1ARFRU4gP/5xZ6F065bY/ocjhS3MlL4f/FRXuzHmxYvd5v8ei4udeM6c6YS0sLClTQp/512FaIJ+otWFhHe2DXS0k3pggAhWAQO990OBgDxe7PbK9gWUhaozGeeGPRowAdXf9ow2npjWePUPBndMRO7CLSZK9xDjWYkcGD916tRZ5wb45Cc/SUFBQbvPn/RB/UhussJCZwUEWwWBP4TBbNkCTzwBTz8NBw86kbz/fpcO7uKLU2NuYFs+VweR9PshGgoL4ZZb3KbqgrX84vnnP7u/P8CHPuRc7dOnOzd7rEFgRrzJEpE1AfvzVXV+K23C6ULKkbRRXVVVEUmZOS3eH3U+OAsz+HgiBsZVlc2bN7NkyZKzjuXn59MtTtZRUgf1W3OTRbsYcW0tPPussyZXrHCW48c+BnPnOmsjTnPl4kYKL7KcdkEeIjBypNu+/GU3Drp2LSxb5qzQZ56Bxx93dbt1axkPhegT1hvxolFVJ7a1carpQjAd/d+73zO58V6rvfI9QGDIZ5FXFki4OoeAPp7LNbhtcxvveL5XP2Zyc3OZM2cO+fn5AM3jPrltGBdTVbZu3cqvf/1rnn/+eTIzM7n66qvJy8tr97kT3feYCZVibvZsV+7HH41aUuJeA0WlrAy+9CWXXm7uXJe0/Cc/gd27naUxa1bqiaWfSJ8riST1fogHWVkuKOj+++Hll106wJUr4aGHYNKkMz0MPh88+CD87nfOSrV556lIOF1IOTp6DPPHwKGAwd2+qnq/iFyHy7TgD/p5TFUneW22quoYT/DeA2bghHA18GlV3SwifwKeDwj62aSqvxKRrwLjA4J+blLVW1rrd6KiZFWV7du3s3TpUvbu3Uu/fv2YNm0a5513HiKS0KjFhEZFxntdx9On4a9/hV/8wgXyZGfDrbfCF74AU6emhss1zUnZqOn24vO5RBNr17ro2/Xr4e23XZJ8cA9ekyY51/1llyV3nDuYdI7w9fres7RU60+ejNjpaHUh4X1uC6qakA1YiBuDPI0bV7wTN374GvA+8E/cFwNu1etfAh8AZcBEr7w/sC3gnNfiRPMD4N8DykcAq4DtwJ+AHl55tre/3Ts+Ipq+9+zZU+PNzp079YknntDvfve7+uijj+r69eu1qakp7tfpcJqaVDdsUC0pUQX3umGDK1dV3bev5Zh/Kylx5cHs26f6n/+pOmSIqzd8uOrDD6sePNhxn8fofDQ1qb77ruovf6l6zTWqmZkt92KvXqqzZqn+4Aeqy5ap1tcnr4+R/o9SmYC+93Re1bjoQipuKZ0aT0Sux4lc2GS4iSCeqfF2797NkiVL2LFjB7179+aKK67goosuIjNV3Yix0lpkYmtjmKqu7s9/7lysp087N+tXvwrXXJO67lYj/Qh1r/bq5eZ+btvm9rt1c8kULr/cbZde6pLMJ6Nv6RLhG9D3XKBOtRO4KkKT0oKZLOIhmFVVVSxZsoT33nuPnj17cvnllzNx4sS4BfLElfa4gqJxuYY6/+nT8Kc/wbx5bgJ7fj58/vMuqGPUqHh9MsNoIdK92ru3c92+9ZbLDLVmTUvw0JgxTjwvu8y9nnNO/IcF4jFvN5Eu3SiHXTq7YKZoWFz6cvDgQZYuXcrmzZvJzs7mqquuYvLkySGnqqQE7Z3sHU16usAUc1VV8L3vuVU/9u93P0a/+hXcdpt72jeMRBHpXu3b1yW5uOEGV37ypBPNt95y2/PPuwhtgIEDnXj6t4suav+ybtH8H0UikUkbWjt3qL53UszCDMG4ceN0xYoVMQVCHD16lGXLlrFx40aysrKYMmUKl156KdnZ2UDigyxaXZYo3BNie11B0f6jrlnjrMlnn3VP7tdeC/fe66aEpEtgg9ExJMpSao+o+Hzwr3+1WKFvv92SjSgnxwUS+QV06lSXZSrefYv0vSTSpRvDsEtuRUWntjBNMEMwZMgQvf/++6NKF1ZbW8sbb7zB2rVrEREuueQSLr/88jNC9DXBqchaXfg20j/jrl2JcwU1NsILLzihfOcdZ0F+7nMu6fm557b7cxudkESnt4unGO/b54TTv61f7+55cCveBFqhI0a07saN1LfWvpdEpmKMYdglmijZdMYEMwRDhgzRL37xixHThZ04cYK3336blStX0tTUxEUXXcS0adOa51IGkuhUZLWHD7Pg6afPPv/tt9Orb9/IT4gQ/yfTw4fht7+FX/7SCfKIEXDPPU4svbl/hhGSdA5+qauDVataBHT5cjdvGNrvxm3te0mmhRlAOq2H2RZsDDMCodKFnTp1ihUrVvDOO+/Q0NDA+PHjmT59On379g17nkSnImv0+UKf3+etWBNp3chhw1pP4RbtU/mWLfDYYy5l3YkTLtfnL37hVgexaFcjGtqyxmmqkJvr7vkrr3T7Pp9b3izQCn3hBXcsO9utxHLppS1u3EjRuK19L4lMxZjCaR47GhPMCASmC2tsbGTNmjW8+eab1NfXM3r0aK688koGDmw97WGiU5FlZWSEPr9f1CIFFLSWwq01V5DPB6++6tyur77q2n/mM2588vzz4/L5jC5Ee4NfUgn//9b48S5bFbS4cd95x20/+xn86Efu2OjRLQJ66aVuP5r/4cBrJSIVYwqneexozCUbgsAxzP79+7Np0yaWLVtGTU0NpaWlzJgxg6KioqjPl9JjmK3d9OHcMa+/3iKU27a5DCpf+Qp88YswYEC7P5PRRelqS3SdOOEC4vwW6DvvuCENcIFDU6c68ZwyxS24PWdOSn8vnd0la4IZgnHjxuny5cuprKxk6dKlHDx4kCFDhjBjxow2L7WVslGyrRFuwD8vD44dc5O8v/51uPnm9ofWJ4J0TjfWVenKfzNVeO+9Fgt0+XLn1gX3HYwd68Y/J02Cj37UJaRPoZSGJphdkPz8fP3xj3/Mvn376N+/P1dddRVjxozpHLk2Y6WqyiW6rqw8s/yGG+CBB1I7t2tXs1aMzok/ubxfRFeudKv3gHuYmDLF/R9OneoeYJOYRL+zC6aNYYYgNzeX+vp6Zs+ezfnnn09GV/1xPX3aLczcu3dLWV4eLFrknm5T/XsJt1JKOkRcGuHpahZoQYFLFzlrlttvanJW5/LlLdtLL7ljmZkudmDqVCekU6aknBWazpiFGYKCggI9cOBA6q4PmGgOHYL5812E6969bs7k5z7n1p/s2zd9fqASOTfNSA7mNQjNoUPuQXD5cve6ahUcP+6O9evnvER+AZ00KWHTuzq7hWmCGYJ4Jl9PK7ZscUE8Tz/tUoPNnAn33eeebNPxxyid5/QZobG/aXQ0NbnMRCtWtGxbtrgxUhGXknLyZCeekye7B4445Lk2weyCdCnB9PngH/+ARx91Ua/Z2S6v6733wrhxye5d+zBrpPNhXoO2U1PjLM+VK1u2AwfcsZwcmDDBiefkyW6OaGlpzK5cE8wuSJcQzPp6eOYZJ5Rbt3beaSFdbbyrs2MWZvxQdQ8agQK6bl1LMoT+/Z0F6t8uuaTVpc5MMLsgnVow9+xxKet+8xs332vCBDct5JZbUnNaiGEEYl6DxHLqFLz7rrNE/ZvflQswfLiLxL3kEvd68cUuENDDBLML0ikFc/VqZ00+95z70bnxRieUl11mEXRGepHKXoNU7lskIvX7+HFnea5c6ZIsrF7tLFM/o0c3C6jcd58JZlej0whmeHDe9wAAD/xJREFUUxP89a/wyCMui0jv3jB3rkuEPnx4sntnGJ2LdLV+29LvgwedePoFdM0a2LsXARPMuF9UpBw4DjQBjao6UUT6As8CpUA5cIuqHgnR9kHgTq/tvar6qlc+C5gHZAILVPUhr3w4sAjoB6wFblPVU5H6l/aCeewYPPGES4ReXu7E8Wtfc1NDQqymYhhGHEjX8dV49XvvXmTo0FYFM9Tvf1u6nQyS+dhzpapeGPBlPQC8pqqjgNe8/TMQkbHArcB5wCzgVyKSKSKZwC+Ba4CxwByvLsCPgEdUdSRwBCe2qYfP527cigr36l9pJBZ27nRu1qIi+MY33EokL7wA77/vBNPE0ujkqCq1tbUcPXqU2tpaOtQgSNeVVuLV7yFDYqkd/PufFqSSn2A28Hvv/e+BG8PUWaSqDaq6E9gOTPK27aq6w7MeFwGzxeWyuwr4cyvnTS5+l8iUKS6Ue8oUtx+NaKo6d+snPuEyevziFy7BwJo18MYb8PGP29JaRpfAv8jBggULmDdvHgsWLKC6urrjRNO/okgg6bDSSrr2OwkkSzAV+IeIrBWRu7yygaq6z3tfBYRaN2sosCtgf7dXFq68H3BUVRuDylOLcCncqqvDtzl9GhYudHOmLr8cliyBb33LuWD/8AcXvdYZiIflbXQJ6urqmlcEArcm7MKFC4nr8Eqk+9G/bqRffILXjUzVe7m1fsefUL//aUGycr9drqp7RKQQWCwiWwMPqqqKSIcOrnp/uLsAunf09IpYXCJHj8Jvf+vGJ3fvdmnrfvUruP32pCZdTgjpGkRhJIVEL9Te6v0Yad3IVL6X47veZZaIrAnYn6+q84PqnPX7r6pvtLn/HUhSBFNV93iv1SLyF5xLdb+IDFbVfSIyGAhlXu0BhgXsF3llhCk/BPQRkSzPygysH9yn+cB8cEE/bf5wbSGaRXM/+MClrXvySairg6uugscfh2uvTf4/XKKw5OlGDCR6ofao7seMjND3Zqrfy+H6HTutBvGE+f1PC8Hs8F9aEckVkd7+98DVwLvAS8AdXrU7gBe9OpNE5Gmv/CXgVhHp4UW/jgJWAauBUSIyXES64wKDXlI3eLEEuDn4vClFOJfIgAFufPKmm2DUKPj1r91Y5YYN8NprcP31nVcsIX2DKIykkJuby5w5c8j3Eov7F2rPjZfnpT33o93LQMTf/7QgGRbmQOAv3tqSWcB/q+rfRWQ18JyI3AlUALd49YuBEwCqullEngO2AI3AV1W1CUBE7gZexU0reVJVvVVX+RawSES+D6wHnuiAzxgbwS6RzEx480246y6XaaOgAB58EO6+26Ww6ypEY3kbhoeIUFhYyNy5cxOzUHt77ke7l/2E/P1PbpeiJ+UTF4jIj4FnVHVTR10zafMwjx2DBQuc67Wy0kW9fv3rcMcdnW98MhpSedzH6Hq0537sIveypcbrgnS4YO7a5UTyt791onnFFW4e5Q03dKp/pjaRrqnGjM5Je+7HLnAvm2B2QTpMMNetg5/+1OV3VYVPftIJ5SWXJP7ahmEYcaazC2ayppV0SlSVurq6yOMnPh/8/e9OKF9/HXr1crldv/a1yOv5dYGnU8MwjFTGBDNO+LOM+CdO+yP0CgsLnWg2NMAf/+iEcssWGDoUHn4YvvAF6NMn8sm7yPiHYRhGKmO/tnEibJaRmhr44Q9dyrs773RrTv7hDy7v67/9W+tiCW3LBGQYhmHEFbMw40TYLCMVFfDtb8PVV8Mzz8CMGbGvP2lzuAzDaC82rNNu7NuKE1mZmc0Tpv3k5+eT5fO5RAOvvgof+UjbFmu25MiGYbSH9izwYDRjghkPVMldvpw5M2eemWXk6qvJHTYMLrigfefv+OTIhmF0JmxYJy6YS7Y9nD4Nf/oTPPwwsnEjhVOmMPd736NxwACyDhwg9847kUWLoH//9l0nvsmRDcPoatiwTlwwwWwLdXXwxBPws5+5m27MGPjZz5BHH6XX1Ve31Iun2zR+yZENw+hqWGq+uGAmSiwcOAD/8R9QXOzmTQ4bBi+9BJs3u/2XXjK3qWEYqYcN68QFy/QTgrMy/ezc6eZPPvkknDjhfP/33w+XXnpmQ4tCMwwjVemA36fOnunHBDMEzYK5caNLLvDss+7Guu02N3dyzJhkd9EwDCPl6OyCaWOYoWhqgmuucSnsevWC++5zq4YMHZrsnhmGYRhJwizMEOSKaF1hoRuX/PKX3XqUhmEYRkTMwuyKdO8O5eWQk5PsnhiGYRgpgglmKDIzoabGDYwHDopbUI9hGEaXxX7tQ3HixNmpoyy1lGEYRpfGxjBDkCuideDmKq1Y4RIGVFU5kQye+Os/bhiG0cXp7GOYXcLCFJFZIrJNRLaLyANRNwxMHdXe1FI+nxPdigr3apapYRhdkDb/HqcAnV4wRSQT+CVwDTAWmCMiY6NqHJg6qj0rhpg71zAMo32/xylApxdMYBKwXf9/e/ceK0dZxnH8+zv2KgWaIjFcjK2gNuV2WkoJhkqtl5KI4ahV2gi0UKiQoH8YiCGNLYkUI2qq8RJsoFS5WKSGWEtpU0u1rdpC7b1aexEqSkubcpFePKX08Y/3XRiWnj1zdmf27M55Pslkd97ZmX2f2Z33nXl39n3N/mlmR4F5wNWdrlXedVQtXUv5SAHOOQfVlscNoifcJXsW8EJi/t/ApeUvkjQVmFqa77937/GjI0YcPX787cvAlpaWlj69e/dR376yvXutfHlH+vft21/t7QI4Rtzpu3djgwfbkfb2I1VH1rP0Iuy+IipqbEWNC4odWy3eK2ltYn62mc1OzKcqjxtVT6gwU4kf6uxOX1gjSWvbzUbm/T5FI2mtFXS/FTW2osYFxY6tFkXfLz2hSfY/wAcS82fHNOecc/XV1OVxT6gwnwU+LGmIpD7ABGBBN+fJOed6oqYujwvfJGtmxyTdBiwB3gPMMbOt3Zil3Jt9C6rI+62osRU1Lih2bLWouF8asDzuEu+4wDnnnEuhJzTJOuecczXzCtM555xLwSvMLpJkkh5OzPeStF/Swhq3+w1Jf5O0SdIySR+M6a2S/iJpa1x2Ta0xdCdJbXEfDs1gW5Mk7YjTpET6TEkvSDpY63t0IS/1iGuipM3xe7BY0vtqfa+U+ckytsWSXi0/XiTNlfScpA1xaq31vVLkJZO4Kh2jkh6J3cBtkTRHUu/ac569vMq1xPaGSFoTu8N7LN7wg6SPS1on6Zik8Vm8V568wuy6Q8D5kkqDZX6aLt4WLelEN1utB0aa2YXAfODemH4YuN7MzgOuBH4oaWBVOW8ME4FV8TG12KVWcn4QMIPwp+dRwAxJpZG+fxfT6inXuOJ35kfAJ+J3ZBNwWxYZTyGT2KLvAdd1sModZtYapw1dzGM1soqr0jH6CDAUuADoD9xUU47zU3O51onvArPM7FzgFWBKTP8XMBl4NMP3yo1XmNVZBHw2Pp8I/Kq0QNKoeLa5XtKfJX00pk+WtEDS08Cy8g2a2XIzOxxnVxP+n4SZbTezHfH5i8A+4PS8AsuTpAHA5YSDZUJMGyNphaQn45n4fZJa4rKDkn4gaSNwWdnmxgFLzexlM3sFWEoorDCz1Wa2p2BxKU4nSRJwCvBik8WGmS0DXs87353JMq5Kx6iZLbIIeIZ4XDeoasq1FcnWAEmrJF2U3Gj8vo4lXAgA/AJoAzCz581sE9AUHWt7hVmdecAESf2AC4E1iWXbgNFmNhyYDtyTWDYCGG9mV3Sy/SnAU+WJkkYBfYBdNeS9O10NLDaz7cABSRfH9FHA1widMZ8DfCGmnwSsMbOLzGxV2bZO1MXWWbnlvLLc4zKzN4Bbgc2EinIY8EAewZTJMrbOzIxNmrMkpRjVoCa5xNXRMRqbYq8DFmcaRbaqKdceIFwhIukjQD8z21i23dOAV82s1JVgdx6rNfEKswrxjGgw4SxsUdniU4HHJW0BZgHnJZYtNbOXK21b0rXASELTVTL9DOAh4AYza4qzsROYSDgoiY+lprBnYmfMbxLOai+P6W8Cv6lvFquSe1yxwL0VGA6cSWiSvbPGfKdRr8/sTkLT5SXAIOCbVec4nczj6uQY/RmwwsxWZpH5PFRZrj0OXBW/nzcCc+uS2W5S+I4LcrQA+D4whnAGVfJtYLmZfV7SYOAPiWWHSk8kzSQ2f5hZa0z7FDANuMLM2hOvPQV4EphmZquzDyV/8be5scAFkozwp2UjxFX+Z+DS/P9iwYWkS4Gfx/TphN9XxiTWOZt37uu6qGNcrQBmtiuu92sg17EEs47NzDrs0SXRhN4u6UHg9myieLc84qp0jEqaQWii/Woe8WSsS+WamR2WtJRwxf5l4GIASUuA9wNrgZuBgZJ6xavMpuoOL8krzOrNITQzbJY0JpF+Km9/GSZ3tLKZTSNUjgBIGk44CK80s32J9D7AE8AvzWz+uzbUPMYDD5nZW4WGpD8Co4FRkoYAu4FrOEFvIWa2hlhpxHUHAfckbvT5DPW54ipXr7j6AcMknW5m+wk3Zfw9n5DekmlslUg6w8z2xN+72oAtGeS/I1l/Zh0eo5JuIvwu/ckmaRmqply7n3Cj3cr4uztmNi75AknLCft9HjAJ+G3mOa8HM/OpCxNw8ARpY4CF8fllwHbCXa93A8/H9MnATyps9/fAS8CGOC2I6dcCbyTSNwCt3b0fqthvywknA8m0rxMK/RWEs/N/APcBLR3t67L1bwR2xumGRPq9hN9JjsfHuwoS1y1xu5sIBdRpTfiZrQT2A0fiZzMupj9N+H12C/AwMKBZ4qp0jBKGANuVSJ+e52dWwz6pqlxLvHZb+T4tW/4hwk1POwnNuH1j+iXxe3AIOABs7e59UWnyrvFct4pnsbeb2VXdnZcsFTUuKG5sRY0rb5LOJDTRDrXmuIqumt/045xzriqSrifcTTut6JUleOfrzjnnXCp+hemcc86l4BWmc845l4JXmM4551wKXmE614Ak3SWpwz/vK4y0MayeeXKup/MK07nm1Ebo79Q5Vyd+l6xzDULSNEIvKPsIHbD/FXgNmEro0HsnoQPvVmBhXPYa8MW4iZ8SumA7DNxsZtvqmX/nis4rTOcaQBwtYy5hHMxewDpCTzMPmtmB+Jq7gZfM7MeS5hJ6YZkfly0DbjGzHbGv0++Y2dj6R+JccXlfss41htHAExbHRJVU6qj8/FhRDgQGAEvKV4xjO36MMJpEKTnv4bGc63G8wnSusc0F2sxso6TJvHMkk5IWQofZqTo6d85Vx2/6ca4xrADaJPWXdDLwuZh+MrAnjjf4lcTrX4/LMLP/As9J+hKEEe7LR713ztXOK0znGoCZrQMeAzYCTwHPxkXfIvTV+SfCiBAl84A7JK2XdA6hMp0iaSOwlTA+oXMuQ37Tj3POOZeCX2E655xzKXiF6ZxzzqXgFaZzzjmXgleYzjnnXApeYTrnnHMpeIXpnHPOpeAVpnPOOZfC/wHOJVLvMrGT/AAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax1 = plt.subplots()\n",
"ax2 = ax1.twinx()\n",
"\n",
"seaborn.scatterplot(df.index, df.viral_rna, ax=ax1, color='red')\n",
"ax1.plot(df.index, df.smoothed_viral_rna, color='red')\n",
"ax1.set_ylim(0, 400_000)\n",
"ax1.set_ylabel('virus RNA/mL')\n",
"ax1.yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}'))\n",
"\n",
"seaborn.scatterplot(df.index, df.hospitalizations, ax=ax2, color='gray')\n",
"ax2.plot(df.index, df.smoothed_hospitalizations, color='gray')\n",
"ax2.set_ylim(0, 40)\n",
"ax2.set_ylabel('Hospital Admissions \\n (Patients/Day)', rotation=270, labelpad=30)\n",
"\n",
"ax1.set_xlim(df.index.min(), df.index.max())\n",
"ax1.xaxis.set_major_formatter(mdates.DateFormatter('%b-%d'))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<style type=\"text/css\" >\n",
"</style><table id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bc\" ><thead> <tr> <th class=\"blank level0\" ></th> <th class=\"col_heading level0 col0\" >model_rsquared</th> <th class=\"col_heading level0 col1\" >smoothed_model_rsquared</th> </tr> <tr> <th class=\"index_name level0\" >offset</th> <th class=\"blank\" ></th> <th class=\"blank\" ></th> </tr></thead><tbody>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row0\" class=\"row_heading level0 row0\" >0</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow0_col0\" class=\"data row0 col0\" >0.101</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow0_col1\" class=\"data row0 col1\" >0.923</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row1\" class=\"row_heading level0 row1\" >1</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow1_col0\" class=\"data row1 col0\" >0.426</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow1_col1\" class=\"data row1 col1\" >0.972</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row2\" class=\"row_heading level0 row2\" >2</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow2_col0\" class=\"data row2 col0\" >0.199</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow2_col1\" class=\"data row2 col1\" >0.998</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row3\" class=\"row_heading level0 row3\" >3</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow3_col0\" class=\"data row3 col0\" >0.185</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow3_col1\" class=\"data row3 col1\" >0.993</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row4\" class=\"row_heading level0 row4\" >4</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow4_col0\" class=\"data row4 col0\" >0.165</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow4_col1\" class=\"data row4 col1\" >0.954</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row5\" class=\"row_heading level0 row5\" >5</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow5_col0\" class=\"data row5 col0\" >0.075</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow5_col1\" class=\"data row5 col1\" >0.878</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row6\" class=\"row_heading level0 row6\" >6</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow6_col0\" class=\"data row6 col0\" >-0.006</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow6_col1\" class=\"data row6 col1\" >0.764</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row7\" class=\"row_heading level0 row7\" >7</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow7_col0\" class=\"data row7 col0\" >-0.026</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow7_col1\" class=\"data row7 col1\" >0.617</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row8\" class=\"row_heading level0 row8\" >8</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow8_col0\" class=\"data row8 col0\" >-0.028</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow8_col1\" class=\"data row8 col1\" >0.447</td>\n",
" </tr>\n",
" <tr>\n",
" <th id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bclevel0_row9\" class=\"row_heading level0 row9\" >9</th>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow9_col0\" class=\"data row9 col0\" >0.026</td>\n",
" <td id=\"T_1afa4f3c_9fd3_11ea_8a2f_e4a4718180bcrow9_col1\" class=\"data row9 col1\" >0.273</td>\n",
" </tr>\n",
" </tbody></table>"
],
"text/plain": [
"<pandas.io.formats.style.Styler at 0x7f4153030390>"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model_summary = pandas.DataFrame(columns=['model_rsquared', \n",
" 'smoothed_model_rsquared', \n",
" 'point_model_summary', \n",
" 'smoothed_model_summary'])\n",
"model_summary.index.name = 'offset'\n",
"\n",
"for offset in range(0, 10):\n",
" point_model = sm.ols(\n",
" formula='hospitalizations ~ viral_rna', \n",
" data={'hospitalizations': df.hospitalizations, \n",
" 'viral_rna': df.viral_rna.shift(offset)})\n",
" \n",
" smoothed_model = sm.ols(\n",
" formula='smoothed_hospitalizations ~ smoothed_viral_rna',\n",
" data={'smoothed_hospitalizations': df.smoothed_hospitalizations, \n",
" 'smoothed_viral_rna': df.smoothed_viral_rna.shift(offset)})\n",
" \n",
" point_fit = point_model.fit()\n",
" smoothed_fit = smoothed_model.fit()\n",
"\n",
" model_summary.loc[offset] = [\n",
" point_fit.rsquared_adj, \n",
" smoothed_fit.rsquared_adj,\n",
" point_fit.summary(),\n",
" smoothed_fit.summary(),\n",
" ]\n",
"\n",
"model_summary[['model_rsquared', 'smoothed_model_rsquared']].style.format('{:.3f}')"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>OLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>hospitalizations</td> <th> R-squared: </th> <td> 0.206</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>OLS</td> <th> Adj. R-squared: </th> <td> 0.185</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 10.10</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 27 May 2020</td> <th> Prob (F-statistic):</th> <td>0.00290</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>00:06:21</td> <th> Log-Likelihood: </th> <td> -126.38</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 41</td> <th> AIC: </th> <td> 256.8</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 39</td> <th> BIC: </th> <td> 260.2</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 1</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 13.3979</td> <td> 1.476</td> <td> 9.078</td> <td> 0.000</td> <td> 10.413</td> <td> 16.383</td>\n",
"</tr>\n",
"<tr>\n",
" <th>viral_rna</th> <td> 3.846e-05</td> <td> 1.21e-05</td> <td> 3.178</td> <td> 0.003</td> <td> 1.4e-05</td> <td> 6.29e-05</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td> 2.062</td> <th> Durbin-Watson: </th> <td> 2.030</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.357</td> <th> Jarque-Bera (JB): </th> <td> 1.138</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td>-0.354</td> <th> Prob(JB): </th> <td> 0.566</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 3.405</td> <th> Cond. No. </th> <td>2.13e+05</td>\n",
"</tr>\n",
"</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 2.13e+05. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: hospitalizations R-squared: 0.206\n",
"Model: OLS Adj. R-squared: 0.185\n",
"Method: Least Squares F-statistic: 10.10\n",
"Date: Wed, 27 May 2020 Prob (F-statistic): 0.00290\n",
"Time: 00:06:21 Log-Likelihood: -126.38\n",
"No. Observations: 41 AIC: 256.8\n",
"Df Residuals: 39 BIC: 260.2\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 13.3979 1.476 9.078 0.000 10.413 16.383\n",
"viral_rna 3.846e-05 1.21e-05 3.178 0.003 1.4e-05 6.29e-05\n",
"==============================================================================\n",
"Omnibus: 2.062 Durbin-Watson: 2.030\n",
"Prob(Omnibus): 0.357 Jarque-Bera (JB): 1.138\n",
"Skew: -0.354 Prob(JB): 0.566\n",
"Kurtosis: 3.405 Cond. No. 2.13e+05\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 2.13e+05. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model_summary.loc[3].point_model_summary"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>OLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>smoothed_hospitalizations</td> <th> R-squared: </th> <td> 0.993</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>OLS</td> <th> Adj. R-squared: </th> <td> 0.993</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 5899.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 27 May 2020</td> <th> Prob (F-statistic):</th> <td>3.51e-44</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>00:06:21</td> <th> Log-Likelihood: </th> <td> -2.2840</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 41</td> <th> AIC: </th> <td> 8.568</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 39</td> <th> BIC: </th> <td> 12.00</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 1</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 7.8792</td> <td> 0.124</td> <td> 63.293</td> <td> 0.000</td> <td> 7.627</td> <td> 8.131</td>\n",
"</tr>\n",
"<tr>\n",
" <th>smoothed_viral_rna</th> <td> 0.0001</td> <td> 1.38e-06</td> <td> 76.805</td> <td> 0.000</td> <td> 0.000</td> <td> 0.000</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td>18.865</td> <th> Durbin-Watson: </th> <td> 0.030</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.000</td> <th> Jarque-Bera (JB): </th> <td> 3.384</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td> 0.075</td> <th> Prob(JB): </th> <td> 0.184</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 1.601</td> <th> Cond. No. </th> <td>2.73e+05</td>\n",
"</tr>\n",
"</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 2.73e+05. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" OLS Regression Results \n",
"=====================================================================================\n",
"Dep. Variable: smoothed_hospitalizations R-squared: 0.993\n",
"Model: OLS Adj. R-squared: 0.993\n",
"Method: Least Squares F-statistic: 5899.\n",
"Date: Wed, 27 May 2020 Prob (F-statistic): 3.51e-44\n",
"Time: 00:06:21 Log-Likelihood: -2.2840\n",
"No. Observations: 41 AIC: 8.568\n",
"Df Residuals: 39 BIC: 12.00\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"======================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"--------------------------------------------------------------------------------------\n",
"Intercept 7.8792 0.124 63.293 0.000 7.627 8.131\n",
"smoothed_viral_rna 0.0001 1.38e-06 76.805 0.000 0.000 0.000\n",
"==============================================================================\n",
"Omnibus: 18.865 Durbin-Watson: 0.030\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 3.384\n",
"Skew: 0.075 Prob(JB): 0.184\n",
"Kurtosis: 1.601 Cond. No. 2.73e+05\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 2.73e+05. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model_summary.loc[3].smoothed_model_summary"
]
}
],
"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.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@AlexRiina
Copy link
Author

trying to understand analysis in https://www.medrxiv.org/content/10.1101/2020.05.19.20105999v1.full.pdf

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment