notebooks/supplementary.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Supplementary notebook\n",
    " \n",
    "This notebook contains the estimation of the dilution parameter d and the plot showing the effect of dilution on a pooled test's sensitivity as a function of the concentration of viral load."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir(\"../\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import glob\n",
    "import json\n",
    "from src import utils\n",
    "from src.utils import get_fig_dim\n",
    "from IPython.display import clear_output\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression, Ridge\n",
    "from matplotlib.offsetbox import OffsetImage, AnnotationBbox\n",
    "\n",
    "sns.set_theme(context='paper', style='ticks', font_scale=1)\n",
    "width_pt = 469\n",
    "palette = sns.color_palette('husl', 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation of dilution parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated dilution parameter: 0.04549745359504318\n"
     ]
    }
   ],
   "source": [
    "se=0.99 # Analytic sensitivity\n",
    "sp=0.99 # Analytic specificity\n",
    "S_val = [0.81, 0.91, 0.93, 0.99]\n",
    "c_val = [0.02, 0.1, 0.2, 1.0]\n",
    "xs = np.array([np.log(c) for c in c_val])\n",
    "ys = np.array([np.log((S-1+sp)/(sp+se-1)) for S in S_val])\n",
    "xs = xs.reshape(-1,1)\n",
    "reg = Ridge(alpha=1.0).fit(xs, ys)\n",
    "d = reg.coef_[0]\n",
    "print('Estimated dilution parameter: {d}'.format(d=d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization of dilution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.linspace(0, 1.0, 100000)\n",
    "S = 1 - sp + (sp+se-1)*np.power(c, d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARYAAACmCAYAAAAf8KCOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2UUlEQVR4nO2dd3wVVfr/33NreiNAEkIJCSShigklK7CiiLuooLiwoljAL7IW1nUVxFVsi2XBZaXYsCzNtiBW8CdSJIIQUigJSYD0ACEhBNJz6/z+uMmQS9oN3CQEzvv1uq+ZOW2eucl87pkz5zyPJMuyjEAgEDgRVUcbIBAIrj6EsAgEAqcjhEUgEDgdISwCgcDpCGERCARORwiLQCBwOkJYBAKB0xHCIrgqaWp6ltVqxWg0UlFRgdVqbWerrh00HW2AQHCp5OfnU1hYiMlkYtCgQXh6egI2UZEkiV27dlFUVERJSQk333wzYWFhbN68maqqKk6cOMH48eMZOnRoB1/F1YnosQiapKVJ2WazmaqqKkpLSzGZTACUlJSQlZVFWloalZWVl9V+S+zatYvz58+Tl5dHYmKiki5JEkePHiUnJ4dx48bh4+NDUlISO3bs4Ny5c9x1111ERkayd+9eSkpKLssGQeOIHosAsL/JJUlSfvUPHjxIdnY2ZrOZUaNG0bNnT6xWKyqVip9//pnS0lIKCgqIiYlh1KhR7Ny5E1dXV44fP86AAQO45ZZbmjynJEnN2mQymaisrMRsNmMymTAYDGi1Wnr06EFqaipnzpxhwoQJlJWVsW3bNnr27EloaCgA6enplJaW4u7uTkREBHFxcWRkZODt7U1VVRVhYWGkpKRw4sQJ/Pz8nPMlChSEsAiAhje5JEnU1NSQmJjIsGHDKCsrY//+/fj6+uLh4UFCQgKnTp1ixowZ7N69m7S0NPLz86msrGTy5Mn4+vqyb98+wsPD6dWrl13bJpOJo0ePsnPnTvz8/Jg4cSK+vr4NbMrLyyM1NRW1Wo1OpwOga9eu9OjRg7y8PGRZRqPR4OPjQ3l5OSdOnFCEpXv37uzevZusrCw8PT0xmUz069eP1NRUUlJScHd3R6VSKY9PAufSqYUlPDyco0ePdrQZVzRWq5XKykosFguenp6o1Wq7/LNnz5Kbm4vRaMRsNlNTU4OXlxcjRowgNTWV0tJSAgMDCQwMZOPGjaSmpjJixAiKioooKiqisrKSkJAQMjIyOH78OFqtlqKiIkJCQti9ezeFhYX06tVL6QGBTTD27NlDTEwMPXv2xGAwNGp73759CQ4OxmQyYTabMRqNSs9Kq9VitVqVdlUqlfI4ZrVaGTVqFAUFBWzbto2ioiL8/PyYOnUqu3fv5siRIxw/fpw+ffoQEhLSht/+tUunFhZB81RUVBAfH09iYiJdu3bllltuISgoyK6MJEnodDrc3d3R6XRoNBpcXV0BOHfuHCaTCZ1Oh1qtxmQyUVxcDNhu+t9++42kpCT8/f2prq5mzJgxxMfHs2vXLnx8fLBarQQHBzewq2/fvvj7+7N48WImT57Mfffd16j9ubm57Nu3D0mS0Gq1yLJM165d8fPzo0+fPqSlpVFSUoKnpydarZY+ffoAoFLZhg7vuusuSktL+eGHH/D39wdg9OjRXH/99ezYsYOAgAAAO9ETOIc2FZb4+HhWrlzJmjVr7NKTkpJISkrCarVy991306VLl0bTBJfH2bNnSU1NJTo6mmHDhjX6erVOTGpqaigrK6O6uhoXFxf8/f3x9PREkiSMRiN6vR6NRqM8kkRERHDnnXdy4MABtm/fjp+fH2PGjCE4OJgjR44QFxeHv78/gYGBgP2jVnV1NePGjSMsLIxvvvmGHTt2cNNNNzWwrU+fPopYXExoaCj9+/fn119/xd/fn+DgYMLCwti0aRMWi4WpU6eSlZVFRkYG4eHhREdHYzQayc/P5+jRo/Tt25cBAwYo40UC59KmwjJ8+PBGu7lLly5l3bp1FBQUsHz5cl555ZVG0wSXR+/evRk5ciSLFy8mIiKCZ555pkGZmpoaiouLUavV6PV6fHx88PDwQJZl+vbtS0pKCvn5+QQHB6NSqRg4cCAA5eXlREdHM2DAAGJjY5Ufgh49euDm5oa7uzuhoaFYrVa+/PJLNBoNU6dOBeD06dNs3boVlUrFgAEDKCsr4+OPP6a6uhpJknj88ccb2Fn3CFR/YHnkyJH4+vri6uqq9D6io6MVoQgODiYwMBC9Xg/YRDQ0NFQZh5FlWYhKPeRqA5aUY8i5BUi9A1EP6o/kqr+kttr8UUir1dodZ2dno1arkSSJoKAgEhMTG027mBUrVrBy5cq2Nveq4syZM3Tp0oX333+fNWvWsG3bNqZMmUJ2djbe3t74+fnh7+/P6NGjG63v7+/P9ddfrwhHv379CAwM5PPPP+fs2bPMmDGDPXv20Lt3bwYNGoTRaKSoqIjc3FxCQkLo1auXMu5x+vRppd0ePXrwxz/+EZVKRXBwMGlpaSQmJjJq1CiSk5MpLCyke/fudrbU7/HU7Xt7ezN8+HC7cvUHiut6V01xtT3+yLKM2WrCLJswWU22fasJs9WM2WrCZDUq+XVp9dMNWTmY9ycz6kxP3H/TwXTQDB98Sba0+xjLuXPn8Pb2Vo5LS0sbTbuYuXPnMnfuXLu08PDwtjP0KsBkMvH9999jMpkYOnQowcHBbNmyBS8vLyorK5XXrBe/aq7P4MGD6dOnD1qtFhcXFwCmT5+u5N92221KGzqdjuDgYGVcxWq1otVqGT9+vN2cFr1eT+/evZVjT09PKioqSE9Px83NrYGoNEdHjo/IsozJarR9LLVb5Ua9sG/LN9Wm1UtX6hhrheDivIZ1zFYTRouhybyWUElqtCotGpUWjVS7VWnRqDRoymrQBproX+qPu1mHnHcaOouw+Pn5UV1drRzrdLpG0wSXT1BQEH/961+V44yMDFJSUnjwwQfRaC786Zu7MdVqNT4+Pi2eq7E2VCoVZrMZtVrdYNAYLoiCxWKhoKCA2267rcGraUfPa5EtGC01GC1G29ZquGi/9lO7b7IaMFgu7BstDY/rl2+sPZPV2KxtGkmDRqVDq9ahVWnRqnRoare2T+2+Wqek1eW76t1qj+uVq923b1Nv145G0ijt2cRDo4iHWqVBLambtNccn4z58x8vfLe9Alr1t7C79kuu2UosFgvV1dX06dNHGXfJz89nxIgRjaYJnEfdDRwcHEx4eDiffvopkydPbpcBco1Gg5+fH1bZSrW5khpzNQZLDQZLde1+NdW6Kn7/0HCOyUkkZ/1GjaW6Xr59WdtxjX2epbrJX2uVpEan0qNT69CpXWr3a49VLujUerQqHXq1bd9N645W5YtO7YK+Nu9CvYvb0NuLhFpvJyAqqXON36gH9YfpIOedRuoVYDu+RKS2dKZ99OhRHnnkET788ENOnjxJQkIC8+bNY+/evaSkpGAwGJg2bRrdunVrNK0lxDyWtsVkNVFtqqDKXEm1uYrq2m1N7bbKXHFRum3fVt4+rcZS3aB9FSr0Gldc1K7o1S64aFzRqV1sx5raNLUrerUrLhoX9OoLabra8nX5eo0LerWLIgR1oqBRaRu5MkFb06bC0tYIYWkeq2yl0lRe71NGhbGcCmMZVZaKi/JsnwpTmW3fXIHRUmPXnkalxVXjjpvGHVeNGy61W9farZvGvUGafXk39GqbQLioXdGotFfdAKrAhpgg18E4+opPlmVqLFWUGs5RZjxPufE8ZbWfcqMtrcxwnjLjhfwKYxlWLsxdcdN44K71wF3ribvWq3briZfehyCPXrhpPXHX2NI8dJ64aTxx13rgUisQWvHrL3AQISwdjDn5KOe++pZzumrOp1VTVtWPsiAXztUUc85wlnM1Z2r3izHU9iBUqPDQeeOl88FL54OnzgdvvS/BniF46YbVpvnipfPBQ+eFh9YLN61HswN3AoEzcUhYPv30U+65554G60yuRVo7icgiWzhbXciZ6gLOVBVwpvq0/X75SYw32AYe3U06fM+n4aftg4++C11dA+jvOxg/vT8+Ll3w1fvjrffDXevZ6QYGBdcWDglL165d+fLLL5FlmZiYGPr27dvWdl2xWFKOXXgl99tBZRJRubGUE+XZnKrM5VRFHqcqczlZkcvpynxMViM6lR5/1wC6ugbQ1S2QSL9hjHULxC+3Gt8tyfgYXdFZ1Wim//GSJyUJBFcKDgnLhAkTAJtjnZkzZ3LHHXcQGRnJH/7wh2uuF2PNPUWhSwU5nufI8ThH3slD5J0/R3H1aTQqLYHuPQly702QRy+iuo0myKMXAe498dX7NzpQKfsbsOgjnPKK71qlsUlyYmFhx+LQW6G///3vFBQUMHDgQB588EF69uxJUVERr7/+Om+//XY7mNk47fFWyCJbyDyfRlrJQdJLDpJWmEiptRRfgyu9K3wI6RVF34gb6OPVn0CPXmIco40oLy+ntLQUs9lMUFCQMomyTkDqlgKUl5cTExNDQEAAW7ZswWQyUVhYyKhRoxgyZEgHX8W1g0M9Fp1Ox3vvvWc3A9PX1xez2dxWdnUoFcYyEgp/5UDRbxw4s5cqUwV9vSOI8BvC6EHz6XfaHd+TBqRBAZe1UOtKp6Vf/To/KQaDATc3t2ZnTFutViRJuuRexM8//4wsy5SXlzNo0CCio6MB28zbwsJCDh8+zLBhw/jtt9+Ii4sjPDycM2fOcNNNN1FQUMDhw4cJDAyka9eul3R+Qetw+FGovqhs2LCBqVOnXlWLAi1WM/GFsew6sYWEwl/x1vlyffcbeGzoCwzxH4Gb1uNC4T4dZmabcSmuKb/77jsAiouLCQ0N5eabb25SOFqzithisWA2m7Farbi6upKTk0NeXh5//vOfKSkpYd++fQQGBtKjRw8AEhMTKS0tJTw8nOrqanbt2kX//v3Jzs4mLi6OgIAA+vbtK0SlHWlWWAoKCli5ciVJSUmsXr3a7p+vbgl8Z6fGXM3W3K/4IetzKk3ljO4xgVdi3iPC77pr6s1La11TZmRkUFRUxO23305VVRUJCQkUFxc3uHllWebkyZPEx8fj5ubGDTfcgIeHBxdz7NgxDhw4gCzLqNVqLBYLgYGB/P73vycjIwOr1YpOp6Nr166cP3+enJwcRVjc3NzIz88nPz+fHj16YDKZ0Gq13H///SxbtgwfHx/mzJnTdl+eoAHNCktgYCDz5s0jMzOTqKgoJf1qiMciyzKxJ39kXepy1CoNk0Pv56Zek3DVuHW0aU7F2a4p09PTiY6OJiwsjK+++oovv/ySAQMGMHDgQLp27drg8enEiRPs3LmT7t27ExgYSGlpaaPC0r9/f/r3vzBwbbValUdtq9WKxWJBkiQ0Gg2yLNstWh01ahTHjh1jy5YtGAwG1Go1wcHB/Prrr9xzzz3k5eXxwQcfMG/ePOHjtp1oUljWrl3LAw88gI+PDydOnCA/P1/JO3LkCM8//3y7GNgWlBtLeefgqxwu3s+0/rO5LeQetOqrb0V1W7qmBLj33nv55JNP2LVrF/fcc0+jNtS5qNyxYwe33XYbQ4YMQZblBk6WcnJy2LNnD2azWVkV3a1bN2655RaCg4NJTU2lvLwcV1dXtFqt3VoyFxcX7r//fk6ePMmPP/5ISEgIZrOZ3NxcbrjhBoYMGUJlZSWpqamMHDnSmV+xoAmaFJaePXsq+2lpaYSGhtotte+slNSc4aXf/oKLxpWlv/+cAPeGPlmvFlrjmrK6utoh15R1/wNGo5GEhAQmT56M2Wxm/fr1+Pv72/3fgM3Xyh//+Ed0Oh3/7//9P7y8vLjuuusaPHr16tWLnj17Njp9oV+/fhw7doz4+Hi6d+9Oly5dlLhABoOBG2+8kbKyMs6dO0d0dDTR0dFotVrCwsLYtWsXvXr1IiAgQIhKO9KkUowbN07Zv/vuu+nXr59yPHny5La1qo2oKi/hpV9n0sXowbygv+GmuroH89rCNeWAAQMAm3uL6upqysvLiYyMZMSIEZw6dYrS0lIqKioYNWoUYOuJJCUl0aVLF/r378+xY8c4deoU586dQ61Wc8899zTrIlKWZbRaLSNHjiQuLg6DwUBMTIzimbDO7aRWq6V79+74+fkpeaNHjyY6OhqdTtfAk6GgbXFoHsuSJUsICwujT58+DBs2rD3scojWzmNZ+vMT5J5J5ZWk8bhYNVf9LNczZ85QUVGBt7c3a9asoXfv3g1cU7bEwYMHFdeU/v7+3HrrrXz66aeYTCYmTJjADz/8QEhICC4uLowZM4YlS5ZQWVnJyy+/DNi8AR48eBCAgQMHkp+fT3x8PIMHDyYtLY077rhDvK25CmmV24TMzEyWL1/OkCFDePjhh9vSLodojbCknT3IC7v/jzcSJtCr0gcA9Q3D0N7ddKS+zs6pU6fYuHGjnWvKrKwsvLy88PHxYdCgQUDzriktFgvl5eWKa8qWZlqXlJRQWFhIZGRko/knTpzg559/pnv37lRWVrbq7WJ9h9r1067VGbayLGORwWwFs1Wu3YJZlrFYL0qXbVuLXVr9erYtwLjeOjx1l/dG1KFBk9WrV3PixAkSEhIYP3684ue0M/Ft5jrGuv1OERW4PNd7nYH2dE0JtnGXmpqaBu4lZVnGarWiVqsxGo0cO3aMqKioVq85a8xOZ4uKVb5wgxottn2TVcZUe/OZLLZjJd3uuLZMbZ7ZcqHuxcd19c1WMFnqlal3rNzwtQJgqXcOsxUsrfCkpJJAI4FaBRqVhKZuWy9NqwKdWmJQV81lC4tDPZb77ruPv/zlL4wePfqK+nVwtMdSZjjHzK0TeH34B4Tma+3W5Vyts2brU/erXlNTw08//URmZiaTJ09WwmC0JzU1NZw4cYKQkJAGvR+zVabGbLupDBYZo8V2c5ssYLTKGOulKVtr3bGtXl2ewe7YvrzJrp0LaSYHb1a1BFoVaNS2m1Fbe6Ne2IJWfSHt4mOljEpCo7avr1VfqKOuvfE1dmJQb1+ybe3EQrpIOGrrqNr5vnVIWEpLS+286BcXFyuR5ToSR4Vlz6mf+Sh5MZ9M2HpFCWNnweZkCqpNMtVmmSqzjMEsY7BAjVmmxmJ/bKhLq3dssKCUq7n42AIGs9zkTS0BWrXt11Rf+6tad6xTS+hUdfu1ebXH+nr7OrXt5tarQaeS7NuoLVP/5q77Ba8vBHVC0t43aWekyUehWbNm8fHHHyNJEi+88IISkkOWZYqKivjpp59abPyTTz7Bz8+PiooKZsyYoaQ//PDDZGVlKWtHtm/fDsDjjz/OwYMHGTduHIsWLbrca1NIPZvEoC5R15SoGC0yFUaZCpNMhdFKpUmm0iRTZbRSabJSY1FRba4VilrBUITDZBMPW75NUC6+5zUqcKm9UfUaCRdN7b66bl/CRWM79nZR2cpqbHV0anDVSOiVchfa0dfd/Bfd8NfS3+5qoElh+eSTT5T9e++9lxEjRihd15SUlBYbTkhI4Pz588yaNYt33nmHQ4cOMXToUCoqKliwYAH9+vXDYDAo640OHz7MPffcwzvvvHO519SAk+U5DOhyvdPbbWtqzDKlBiulBvttuVGmUhGNettaAakw2rr3dUiAu07CTSPhprV9XDUSbhqUfS83Fa5aWxnXuvy6Y01dHZQ8jUrc6IKmcWjwtm7OQR2//vqr8kahKWJjY5XBubqJSkOHDsXDw0OZE7N7924lCl9cXBzr169n1KhRvPzyy8rszzouJxJiQVU+N/e+85LqOhOLVea8QeZstVX5FFdbOVcjc76moYAYLBfqeugkvHUSXnoVnjoJD52Eh1Yi0EOFh07CXWtLd9deyKvbumol0X0XtCstLkJcsmQJaWlpfPHFF8rrPkfGV+pHN9Tr9XZTwetITEzk6aefBmD27NnMnDmTt956i1WrVvHkk0/alb3USIiyLFNSc4YuLi2HE7kcZFnmXI1MQaWF05VWTldYOV1pobDSSnGVTUTOGWSstc8UPnqJLq4quriq8HWR8HdTEearwksv4aOX8NbX7du2oocg6Ey0uAjxzTffJCsri4iICCW9sUDvF1M/umFlZWWDV5Z1a0Lq94Q0Gg3PPPMMzz33XGuuoVmMVgNmqwkPrZdT2qswWskts5BXaiG3zEJuqYW8MisFlRaMFtsbg65uKgLcVQR6qOnvp+F3PVS1ImITEz8XFVq1EArB1UuTwvKvf/2L+fPno9PpiI2NZdu2bUpeRkZGi57jxo4dS2xsLBMnTiQjI4PRo0dTXl6urC6Ni4uzW7tR90q0srLSbiX15VJpKgfAXdv6Va1lBivpZ82kl5g5etZCeomZwkoragl6eKro7aWml7ea0T11BHmoCXRX4e+mEr0LwTVPk8IyefJkZSRekiSioqKU3oWbW8uuBaKiooiLi2Pjxo14eXnh5eXFSy+9xNKlSwGbsNR/tJk+fTpDhgyhX79+TvX10hphMVpkEk+bSDxtIum0iePnLLhqJfr7qgnvouH3vdzo56sm2FMtehwCQTM4NI/FaDQqzndiY2O5/vrrHVpn0tY4Mo8lqzSdp3fdy8Y74hv1R2uVZeJOmfg528Duk7YwHNd313B9dy3DArSE+qjFwKdA0Eoceis0e/ZsXnrpJVauXIler+fw4cP8/e9/b2vbnIqEvTiYLDLfZRj4X1o1Z6utjOut5+UbPIgO1KITvRGB4LJwSFimTp1KSUkJeXl5bNy4kfXr17e1XU6jsQ5Z0mkTi+MqqDbJTB/gym1h+steGyEQCC7g0N3k6urKd999x5IlS0hJSVGWwXcm6nos649U89T2Msb21PH5ZF/uGeAqREUgcDIO9Vhuvvlmbr75ZuV45syZbWZQW7IupZo1yVUsHufJyKCrzxWlQHCl4JCwPPfcc/z222+KI+PKykri4uLa2jYnYXsUOlxk4qNDVSwZ58kIISoCQZvikLAYjUZ27dqlHOfl5bWZQW2BLMM7B6q5s7+LEBWBoB1waHAhPDyc7du3Ex8fT3x8PFu3bm1ru5yGDJit/Th61sz9A11bLC8QCC4fh3os2dnZVFdXKxPk0tPT29QoZ2M03cDQbhr83cQgrUDQHjgkLH/729/Iz88nKiqK5ORkHnnkkba2y6mYLOEM7Sa8tAsE7YVDP+Gvvvoq+/btQ5IkevTowZIlS9raLqchyzIWawhhvs07gRYIBM7DIWEZP3481113na2CSsWPP/7YljY5FYsVrLIv3d2FsAgE7YVDj0J6vZ5ffvmF5ORkvvrqKyZNmtTWdjmNCpMEqPDRi2n6AkF74ZCwTJw4kdDQUCWuUF00vM5ApckmKF56MXArELQXDt1ty5YtQ5Ik3N3dWbhwIZ999llb2+U0jBabsLh0/rDTAkGnwSFh6dKlC7169eKf//wnS5Ysccgfi0AguHZxSFiqqqqYP38+f/rTn/D392fz5s1tbZfTESMsAkH74dADwsyZMykuLiYwMJD8/HzFAXZnwOp4aGqBQOAkHBIWrVZLYGAgAD179mxTg9oKEfBKIGg/rv5XJaLDIhC0Ow6vbr4Un7dNhViFhuFUs7Oz2bJlC66urowbN46QkJBLu6KLaBgcVCAQtDUO9Vhmz55Nbm4uCxYsYPv27axevbrFOnUhVu+8805KS0s5dOiQklcXTnXPnj1KjObXXnuNhx56iBkzZvDvf//70q6mEYSsCATtj0PCUt/n7RtvvEG3bi1HFWwsxGodcXFxvPDCCzz77LNUV1dTU1NDfn4+7u7u6HQ6Tpw4gdlstmtvxYoVhIeH230EAsGVSZv5vG0uxOrs2bPZvn07vr6+rFq1itLSUtzd3ZV8tVpNSUmJXXtz587l6NGjdh9HsPVYrC2UEggEzqTZMZa9e/ei1+sZMWKEnc/bt956q8WGWwqxWj+cqq+vL0ajUcmrqanBy8s5IVEFAkH702yPZc+ePcTGxhIbG4vFYmlVw2PHjlV6FfVDrMKFkBx14VR1Oh1BQUFUV1djMBgICAjAxcXlUq6nAbZTiZEWgaA9aVZYnnnmGf72t78hSZJd8Pb33nuvxYajoqLQ6/UNQqyCLZzq66+/ztatW5Vwqs888wwfffQRq1evdmpQeIFA0P40G2K1oKCAJUuWkJaWRteuXZWehr+/P//5z3/azcimcCTE6g+ZR3hzb1d2z2h5wFkgEDiHZsdYAgMDefPNN8nKyiIiIqK9bHIu4ilIIGh3WnwrpNPp+Pzzzzl27Bhff/01EyZMYNWqVe1hmxMR6iIQtCcOvW4eMmQI3bt3Z8WKFXz22WfKuqHOgJh5KxC0Pw4JS2lpKXPnzuXxxx+ntLSUtWvXtrVdTkaIi0DQnjQ7eFsfk8mEVqulpqYGk8mEp6dnW9vWIo4M3n6XkcLiff7snhHQTlYJBAKHeiz/+c9/ePHFFwHYt28fSUlJbWqUQCDo3DgkLFqtlokTJwJw44038q9//atNjXIm4iFIIGh/HHKb0KNHD8xmM4WFhXz44YdXxGOQo4iZtwJB++NQj+X3v/89+/fv54UXXgBg5cqVbWqUQCDo3DjUY0lOTubZZ58FbOt8Nm7cqEzFv/KRkUSPRSBoV5oVloKCAlauXElSUhIff/wxcGEBYWcRFuFLWyBof1qc0j9v3jwyMzOJiopqL5sEAkEnp8UxFh8fn04tKqLDIhC0P1e/l35AyItA0L5ckrDs37/f2Xa0GSJgmUDQ/jQ5xnLnnXfi4uKCVqslOzuboKAgXFxckGUZi8XSqQLDix6LQNC+NCksK1asUKIefvHFF9xzzz1K3gcffND2lgmuKqxWK6dOncJkMnW0KYJWotVqCQoKQqVy/AGnSWGpH0r1/Pnz5OTk4O/vT1JSElu3bmXOnDmXZ63gmuLUqVOKi1JB56KsrIxTp04RHBzscB2HJshNmzaNt99+mwMHDtCtWzf++c9/OtR4U5EQf/jhB9asWUNFRQWLFy9m8ODBACxatIgff/yRiIgIZd7M5SKGWK4MTCaTEJVOipeXF2fPnm1VHYf6NiqVikmTJvHdd9/x5JNPMmDAgBbrNBUJUZZlXFxc2LBhA7NmzWLFihUAFBYWEhkZyZ49e5wmKhcQ6iIQtCcOCcvzzz/Pvn37kCSJHj168Nprr7VYp6lIiJIkMX78eAAGDx5M165dAVsMo3fffZc5c+Y0CFYGlx4JUUiKQND+OCQs48eP57rrrrNVUKn48ccfW6zTXCTEOvbu3cvMmTMB21uobdu2ERMTw5tvvtmg7KVGQgRAEvIi6FhmzpyJ1XppETmNRiOLFi3im2++ca5RbYhDwqLX6/nll1947733mDp1KpMmTWqxTkuREHNzcwkKCiIsLExJkySJhx56yKlvDhx0kCe4yklISGDo0KGsWLGCTz75hL/+9a9s3LixyfIVFRWUlZU57fwfffRRk29VrFYrp0+fbrKuTqejX79+lyxMHYFDg7cTJ04kNDSUzMxMli1bRmRkZIt1xo4dS2xsLBMnTrSLhOjp6UlxcTHp6enceuutVFZWIkkSrq6uSJKE0WhUBnOdheTU1gSdkejoaLp06cIjjzyCXq+nqqqKKVOm4O3tzS233NKg/LJly3jwwQedNuBcP+Dfxaxdu5bIyEgCApp2n6rVap1iR3vRpLB8++23jf7aHz9+nK+//lrxzdIUUVFRxMXFNYiEuHDhQh5++GHUarUyH+arr77iySefxNvbm8jISLs5M5eL6K8IGsPNzY0HHniAzz77jJiYGF577TWio6PJy8tjxowZpKSk8NNPPzFx4kSWL1+u5D311FNKG5s3b2bHjh14e3uzb98+li9fTlhYGP/973/x8fEhISGBp59+mqKiIl555RU++ugjFi1aRJ8+fTh69Ch6vZ4XXniB+Ph4iouL6d27N5s3b8bLy4vk5GReffVVO5tTUlLYvHkzsizz/vvvs3DhQgwGA2+//TY//fQTRUVF3H///e39VTZKk8KSlpZGWFhYs0rbEo899pjd8dKlSwGbaF3M8uXLL/k8LSPkRdCQHj16UFhYSFlZGSNGjGDChAnMnDmTp556it69e3Prrbciy3KDvDqGDBnC999/z7///W++/PJLVq1axaRJk6isrGTmzJlotVpWrlyp+It2d3ene/fuhIaGMmfOHMaMGYO7uzsRERGMGDECf39/fvnlF1asWNForz0kJIQXX3yRxx9/nB9//JEFCxbw6KOPAlBcXOzUH+TLpUlhmTdvnp2oFBcXk5OTQ1hYGHfddVe7GOcMhKQImqKgoICePXsSGBhIdXU1+/fvb/DI0VyeJEnK2GFUVBRbtmwhNTUVd3d3ACIiIvjyyy8B0Ghst5parVZcu17cnkajYfLkydx111088cQTDaKP1rU7fPhwzp49i6+vLwMHDmTXrl1YrdYr6nGpycHb+qKyYcMGZs6cyfr163nxxRf54Ycf2sU45yAj5EVwMQaDgU8//ZQHH3yQ2NhYMjMzGTdunPL4L0kSsiw3mlcfi8UC2GanDhw4kH79+pGSkgJAVVUVAwcObNGWunNVVVUxbNgw/ve///Hee+81Wb6wsFB5S/vQQw+xcOFCfve737X2K2hTHBq8PXv2LN9//71yvGHDhjYzyNnY/hfE8O21zqFDhygpKWHNmjW4uLhw5MgRHn30UX73u9+RkZHBzp076dq1K1arlYSEBCIiIli+fDlTp05tkBcdHa20m56ezubNm8nMzGT27Nn4+Piwe/duPv74YwwGA3PmzCEvL4+8vDyOHz9Oeno6er2egIAAysvLSU5OJjw8nNWrV/PEE0/w8ssvc//99/OnP/3Jzv6IiAjWrVuHxWIhODiY66+/HrAtvbn99tsJDQ1t1++zRWQH+O9//ytbrVZZlmU5JydHnj9/viPV2pz+/fu3WObTIwfkMeuz294YQbNkZWV1tAlOJz8/X3722Wc77PwGg0E+ffq0/N1337X5uVr793OoxxITE8MDDzxAdnY2AQEBvPzyy20sd85DPAR1TuRqA5aUY8i5BUi9A1EP6o/kqu9os+w4ePAgJ06coKSkBD8/v3Y//1NPPYVarebf//53u5+7JRwOsXrmzBlycnLo169fg8luHYUjIVY/TT3IBwe8ib0vpJ2sEjRGdnY2ISGO/w3M8cmYP78ww1sz/Y9ohjt3fpPAcVr793Oox7JhwwbWrl1LWFgYFouFCRMmcPvtt1+yke2JmHjbOZFzC+yP806DEBaFb7/9lkOHDimvsq80HJrSXzd4+5///Ifly5crU/UFgrZC6h1of9yr6VmpjlA3pX/dunWsWbOGuXPnkpub22T5U6dOXdb5LgdHzj1s2DCqqqqazF+wYAErV64kMTFRSYuLi+Mf//iHXbn6aVu2bGHevHls2rTpEi2/gEPCUueSEmxrfBISEi77xO2F6LB0TtSD+qOZ/kfUNwyzbQf1v6z26qb0T5s2jQcffJAJEyYoEzYv5tChQ3ZvQdsTg8HAkiVLWizniDe3O++8U4mwce7cOfbt22e33ujitIkTJxITE3OJltvTqsHbrKwsAgMDeeWVV5xy8vZBzGPpjEiuetuYShs9/hQVFSkuO95//328vb05fPgwzz33HHv27CElJYXk5GQkSSI3N5fdu3dz88034+7uzpo1axgxYgSbN29mzpw5HDp0iIMHD/LRRx+RnZ3Npk2b8PHxITw8nL59+7Jw4UImTpzI999/z6OPPsqNN96o2FFSUsIvv/zCyZMnMZvN/OEPfyAtLY3Nmzdz2223KeVee+01Bg4cSFxcnMOO1urz3XffMWXKFN55551m05xFi7KXkZGBn58f69atY8+ePWzcuNGhST9XCrIsZrEILrBhwwaWLFnC+vXrlbkiZrOZ6dOno9frOXbsGFFRUURGRjJ48GD++9//4unpyciRIzl+/DjDhw+nqKiImTNnMnfuXBITE5k3bx7XXXcdaWlpnDlzhkmTJjFmzBh27txJaGgoNTU1TJs2jccee4zt27fb2bNq1Sq8vb0ZOnQohYWFhIeH07VrVztRAfD29mbSpEmUlJRw5syZVl3zli1bmDhxIpIkNZvmTJrtsdQtMddqtSxbtsxp3aT2R/RYBDamTp2KXq9n5MiRzJw5k127dhETE8Pnn38OXJhJW8fRo0cZO3YsYHNvoFKpcHd3V1bke3h4ALbhApPJxPXXX8///vc/AgMDlUcMd3d3dDqdUqY+GRkZPPjggwQGBjJ69OgmH3EGDhzIN998g06na2BjS3z22Wd8/vnnGAwGTp48yYYNG/j2228bpDkzbHKzPZbMzEzi4uLYtm0bBw4ccNpJ2xMhKYLGGDRoEOXl5RiNRp599lmmT5+Oj48PsiwrU+ytViuSJLF161YAxQtic7z33nsEBwczaNAgh+zo3bu34op1165dyLKsfOqoqKhg2bJlTJkyxW68s47y8vJm/Q6tX7+edevWsXTpUsaMGcPUqVMbTXMmzQpL//79kSQJFxcXQkNDMRqNGI3GTrZWSCCwn9L/xRdf8Morr/Dyyy/j4eFBr169ePrppzl//jw7d+6kZ8+e/Prrr+zfv5+FCxfyxhtvcO+99xIWFkZycjKFhYWcPn2aQ4cOkZGRQWFhIcePHyclJYX+/fuzcuVKfvnlF7Kzs9m/fz+5ubnk5ORw5MgR8vLyKC0tVex65JFHOHz4MH/4wx+QJAlJkujSpQvvvvuuUkav12OxWHj11VcxmUxs376dlJQUcnJyOH/+PC+++GLrPCq2A81OkBs2bBh+fn4NFmadP3+epKSkdjOyKRyZILc6OYn/HvZm131X2FqKa4zWTrASXB4LFizgiSeeaFXIDkB51TxlyhS7dKdOkPvoo48aDQh/+PBhh09wRSB83gquQb755htiYmIavYcbY8uWLezdu9cpY6nNCktTBg0ZMuSyT9xeCEkRXIs05pC+JSZOnMjEiROdcv5LCgovEAgEzeHQBLlLpalIiNnZ2WzZsgVXV1fGjRtHSEhIo2nOoIurETdtKhDWYlmBQOAc2qzH0lQkRLDNInzooYeYMWOGsuS7sTRnMLirH6N7bW+5oKBN0Wq1Tg2nIWg/ysrKWu32ss16LI1FQhw6dCg1NTXk5+cr/jtPnDhBRUVFgzSz2az4CQXbZL2VK1e22o4Q7/48P3KZE65IcDkEBQVx6tSpVscAFnQ8Wq2WoKCgVtVpM2FpKhJiaWmpIiBg861bWVnZIK2kpIRu3bopaXPnzmXu3Ll253A0zKqg41GpVK1+9SnovLTZo1BTkRB9fX0xGo1KuZqaGtzc3BqkOStQlEAgaH/aTFjGjh2rTF6rHwlRp9MRFBREdXU1BoOBgIAAPD09G6S5uLi0lWkCgaCNcdg15aXw7rvv0q1bN8rLy4mJiWHVqlUsXbqUY8eO8dNPP6HT6bj55psJCwtrNK0lHJl5KxAI2p82FZa2RoyxCAQdS1M/7J1aWBylM/RshI2Xz5VuH1w7NoqZtwKBwOkIYREIBE5HCItAIHA614SwPPHEEx1tQosIGy+fK90+uHZsvCYGbwUCQftyTfRYBAJB+yKERSAQOB0hLAKBwOm0qaOnjqA1zqWuNBt/+OEH1qxZQ0VFBYsXL2bw4I4Lgt6UjXX84x//YPLkyYwcObIDrLPRnI2ZmZkkJibSr18/hg0bdkXZ9/PPP3P+/HnAFo/ojjvu6BD7AOLj41m5ciVr1qyxS09KSiIpKQmr1crdd99Nly5dWtXuVdVjaa1zqSvJRlmWcXFxYcOGDcyaNYsVK1ZccTbWsWPHjmYDkrcHzdmYlZXFhg0bmDZtWoeJSnP2rV27lqlTpzJ16lS++uqrDrGvjuHDh2MwGBqkL126lIcffpjbb7+d5cuXt7rdq0pYGnMuBdg5l9LpdIojqSvJRkmSGD9+PACDBw9W4gpfSTYC5OfnY7FYCA3t2HAqzdn42muvERQUxKJFi0hISLji7Bs4cCDLli0jOTmZe++9t0Psq8/F3uGys7NRq9VIkkRQUBCJiYmtbvOqEpbWOJcqKSm5omysz969e5k5c2Z7m6bQlI1ms5nY2FhuueWWDrOtjqZsrKqq4sSJE8yYMYNZs2bxt7/9zc7XT0fbB/Dkk0+Sl5fH4sWLiY6ObnfbWqK+7YBdgDVHuaqEpTXOpTrKkVRTNtaRm5tLUFCQQ24j2oqmbIyPj+f777/n/vvv5+uvv+b111+nsLDwirLRZDLh4uKCSqUiKCiIbt26NSreHWUfwJIlS3jllVd4+OGH+fvf/97utrVEfdsBdDpdq9u4qoSlNc6lOsqRVFM2AhQXF5Oens6tt95KZWVlh41jNGVjTEwMX3zxBevWreOuu+7iH//4B927d7+ibPT29kan01FZWQnYflQ6wsbm/s7Hjx/Hw8ODG2+8scMeyRvDYrFQUVFBnz59lHGX/Px8RowY0eq2riphiYqKQq/Xs3HjRry8vPDy8uKll14C4JlnnuGjjz5i9erVPPfcc1ecjefOnePhhx/mgw8+YMqUKdx///24urpeUTZeSTRn48KFC1m1ahWbN29m9uzZqNXqK8q+WbNmsW7dOrZu3cqf//zndretPkePHiUvL49jx44RGxvLe++9B8Cjjz7Khx9+yLfffstTTz3V6nbFlH6BQOB0rqoei0AguDIQwiIQCJyOEBaBQOB0hLAIBAKnI4RFIBA4HSEsAqeTnp7OG2+84dQ2jx07xoYNG7jvvvuand+zY8cOVq9e3WJ7e/fubRCyNz09nfvuu6/VtplMJhYtWsQ333zT6rpXK1fd6uarlbqFdX379kWSJOLj45k7d26bx0O2Wq0UFRUREBDQYtlTp04RFBREREQE8+fPd6odq1ev5vXXX6dPnz5oNE3/2950001YrdYW2xsxYgTvvvuuXVpERASXMvtCq9USGhrq0HmvFYSwdAJKSkp48skn+eyzz/D09ASwW8vRlqxdu5bIyMgWhSU/P5/169crkw+dPSnt5MmTgG01bkuoVC13xJuy71Lt1uv1l1TvakUISyfg22+/ZfDgwYqoAIwfPx6j0UhBQQHfffcdrq6ulJeX89BDD7Fo0SL69OnD0aNH0ev1vPHGG5w8eZLY2FiOHTvGoEGDuPvuu9myZQs1NTVs27aNF198kWXLltnVe+GFF4iPj6e4uBi1Ws2rr77KpEmTyMzMZPjw4RQXF3Py5En+7//+j6SkJA4fPkxcXByHDh3CarXyl7/8hS1btlBVVUVaWhqTJ0/G3d2dhQsXMnHiRL7//nseffRRbrzxRuW6ysvLWb9+PV26dCEjI4MFCxYQGxvL6dOn2bBhA1OmTEGtVnP+/Hkee+wxxo8fz6xZs1i8eDGTJ09m06ZNREZGEhERwSuvvMLAgQPx8vJCp9Ph7e3N4cOHee6551pcK3bx9/r444+zadMmu2vu2bMnq1evxt/fn82bN3Pbbbe11b9Ap0OMsXQCsrOzG7hRkCQJvV7PW2+9xe23384DDzxAUlISWVlZdO/endDQUJYuXcru3bsBWLlyJVOnTmXu3LnodDrS09NJTEzE39+fyMhI8vPzG9Rzd3cnIiKCMWPGEB0djYeHB3fffTeLFi2irKyM++67jx49epCYmEhUVBS9e/dm5MiRDBo0CJPJRHl5OV9//TV/+tOfeOSRR3j++ecJDQ2lpqaGadOm8dhjj7F9+3a76/rwww+Jiopi2rRp1NTUsHXrVm688Ua6devG1KlTlR6Fj48Pf/nLX5SeTGhoKOHh4YSFhWG1WomIiKC6upr58+fz+OOPYzabmT59Onq9nmPHjrX4nV/8vSYnJze45gMHDnD27Fluv/12brrpJmf8qa8ahLB0AgICAjh79myjeampqYpLiPDwcDIyMlCr1Urvps7XRnZ2NhqNBj8/P+644w6ysrLw8/Nj7NixzJ07l6ioqEbr1UetVuPr64tarWb06NFs3LgRo9HYYGyh7ubPyspS9rt37865c+eQZVnxi+Pi4oLJZGrxeppi9OjRJCQkkJGRQZ8+fezODeDp6YlOp0OtVhMTE8Pnn38O2BbbtURjdlx8zUlJSYpntY5a13WlIoSlE3DXXXexZ88eO3HJzc2lqKiIfv36kZKSAtjcQURGRjbaho+PDzt27ADgl19+oVevXmzatImSkhKKi4uVNi5GkiRkWW4gHgsWLGDSpEnKyuHGyvXq1YvMzEysVitWq5WQkBAkSWr2WutfT3V1NQMGDABodFBVpVJxyy23sHz5cqKioppt99lnn2X69On4+Pg4NEDb2Pd68TV369aNAwcOKHXE4O0FhLB0AgIDA3nzzTd5/vnnefvtt/nyyy85c+YM3bp1Y968eWzatIn//e9/hIaGEhISQnp6OocOHSI3N5fy8nKSk5N5+umneeutt3jggQcICAhg0KBBTJgwgTvuuIP333+f8PDwRuuFh4ezevVqEhISOH36tPJoFRkZyfz588nPz+e3335TxkS2bdtGSkoKx48fx9XVlTlz5vDGG2+wdu1a5s+fT05ODrm5ueTk5HDkyBHy8vLsHAnNnj2b/fv388UXXyDLMuPGjSMtLY2MjIxGPZndfffdim9gi8VCamoqqampZGVlkZmZSXp6OmATuaeffprz58+zc+dOkpOTKSgosPMnk5+fT05ODtnZ2Q2+14iIiAbXPG7cOEwmEwsWLODgwYNkZmY26IFdq4jVzQKBwOmIHotAIHA6QlgEAoHTEcIiEAicjhAWgUDgdISwCAQCpyOERSAQOB0hLAKBwOn8f6Ep3w+Gprv5AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 280.349x173.265 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "utils.latexify(font_serif='Times New Roman', font_size=10, usetex=False)\n",
    "\n",
    "fig_height, fig_aspect = get_fig_dim(width_pt, fraction=0.6)\n",
    "fig, ax = plt.subplots(figsize=(fig_height*fig_aspect,fig_height))\n",
    "\n",
    "sns.scatterplot(x=c_val, y=S_val, ax=ax, color=palette[0], label='Data points by\\nBateman et al. [44]')\n",
    "ax.legend()\n",
    "\n",
    "se, sp=0.99, 0.99\n",
    "S = 1 - sp + (sp+se-1)*np.power(c, d) \n",
    "sns.lineplot(x=c, y=S, ax=ax, color=palette[2], legend=None)\n",
    "ax.text(x=0.3, y=0.99, s=r'$\\mathregular{s_e=}$' + r'{se}, '.format(se=se) + r'$\\mathregular{s_p=}$' +r'{sp}'.format(sp=sp), rotation=2.5)\n",
    "\n",
    "se, sp=0.8, 0.98\n",
    "S = 1 - sp + (sp+se-1)*np.power(c, d) \n",
    "sns.lineplot(x=c, y=S, ax=ax, color=palette[4], legend=None)\n",
    "ax.text(x=0.3, y=0.81, s=r'$\\mathregular{s_e=}$' + r'{se}, '.format(se=se) + r'$\\mathregular{s_p=}$' +r'{sp}'.format(sp=sp), rotation=2.5)\n",
    "# ax.text(0.3,0.8,r'$s_e={se}, s_p={sp}$'.format(se=se, sp=sp),rotation=2.5)\n",
    "\n",
    "sns.despine()\n",
    "\n",
    "ax.set_ylabel('Pooled test\\'s sensitivity')\n",
    "ax.set_xlabel('Concentration of viral load')\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "fig.savefig('figures/figs1.tif', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e312f55c3c97e43ef328a86651ae8cf61ffd8a19c7888f1a2843b6f17cf5ecdd"
  },
  "kernelspec": {
   "display_name": "Python 3.9.1 64-bit ('env')",
   "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"
  },
  "metadata": {
   "interpreter": {
    "hash": "c425b3d18afb7ef01beae38734bbc4c3f66a0926b3e7f0a803887fb6dc8ff96a"
   }
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}