syadlowsky/confidence-intervals-for-sars-cov-2-antibody-prevalence-using-the-parametric-bootstrap.ipynb Secret

## confidence-intervals-for-sars-cov-2-antibody-prevalence-using-the-parametric-bootstrap.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Confidence Intervals for SARS-CoV-2 Antibody Prevalence Using the Parametric Bootstrap ",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyOvp9xLSzNWEYZ46gjqS2OG",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/syadlowsky/a23a2f67358ef64a7d80d8c5cc0ad307/confidence-intervals-for-sars-cov-2-antibody-prevalence-using-the-parametric-bootstrap.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fDlzRj8g8s3d",
        "colab_type": "text"
      },
      "source": [
        "# Confidence Intervals for SARS-CoV-2 Antibody Prevalence Using the Parametric Bootstrap "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R6P878uLC6k0",
        "colab_type": "text"
      },
      "source": [
        "Here are the values reported in the study."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TaSXJMnloknz",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# manufacturer validation data\n",
        "manufacturer_neg_tests = 371\n",
        "manufacturer_false_positives = 2\n",
        "manufacturer_pos_tests = 85\n",
        "manufacturer_true_pos = 78\n",
        "\n",
        "# currently not used\n",
        "stanford_val_neg_tests = 30\n",
        "stanford_val_false_positives = 0\n",
        "\n",
        "# community test results\n",
        "study_tests = 3439\n",
        "study_positives = 50"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0IMIlvJp8-ct",
        "colab_type": "text"
      },
      "source": [
        "We can calculate the false positive rate (FPR), the true positive rate (TPR), and an estimate of the prevalence using the estimating equation, $$ (1-\\text{Prevalence})\\cdot \\text{FPR} + \\text{Prevalence} \\cdot \\text{TPR} = \\text{Observed Positive Rate}.$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "szE0nH2s3XYz",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "def calc_parameters(val_neg_tests, val_false_positives, val_pos_tests,\n",
        "                    val_true_positives, study_tests, study_positives):\n",
        "  fpr = val_false_positives / val_neg_tests\n",
        "  tpr = val_true_positives / val_pos_tests\n",
        "  observed_positive_fraction = study_positives / study_tests\n",
        "\n",
        "  prev = (observed_positive_fraction - fpr) / (tpr - fpr)\n",
        "  return fpr, tpr, prev"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SdeNu68v5H5-",
        "colab_type": "code",
        "outputId": "375a027c-9b48-4cba-e4a6-92f97042ee8a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "fpr, tpr, prev = calc_parameters(manufacturer_neg_tests, manufacturer_false_positives,\n",
        "                                 manufacturer_pos_tests, manufacturer_true_pos,\n",
        "                                 study_tests, study_positives)\n",
        "print(fpr, tpr, prev)"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0.005390835579514825 0.9176470588235294 0.010028185496404683\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nznlCkRn9zHl",
        "colab_type": "text"
      },
      "source": [
        "The estimated FPR is 0.5%, the estimated TPR is 91.7%, and the estimated prevalence is 1.00%.\n",
        "\n",
        "Now, using these estimates and pretending like they are the true parameters, let's simulate the data again, using the same number of true negatives, true positives, and community members tested. Then, we aggregate these simulations and compute the mean and quantiles to get the parametric bootstrap estimate of the confidence intervals."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0AeKfSOa6CpZ",
        "colab_type": "code",
        "outputId": "a8e12d25-87fd-40a2-d437-18b7e2dcf2a1",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "boot_prev = list()\n",
        "for b in range(50000):\n",
        "  boot_false_positives = np.random.binomial(manufacturer_neg_tests, fpr)\n",
        "  boot_true_pos = np.random.binomial(manufacturer_pos_tests, tpr)\n",
        "  boot_study_has_antibody = np.random.binomial(study_tests, prev)\n",
        "  boot_study_pos = np.random.binomial(boot_study_has_antibody, tpr) + \\\n",
        "                   np.random.binomial(study_tests - boot_study_has_antibody, fpr)\n",
        "  _, _, boot_sample = calc_parameters(manufacturer_neg_tests, boot_false_positives,\n",
        "                                      manufacturer_pos_tests, boot_true_pos,\n",
        "                                      study_tests, boot_study_pos)\n",
        "  boot_prev.append(boot_sample)\n",
        "print(max(np.quantile(boot_prev, 0.025), 0), np.mean(boot_prev), min(np.quantile(boot_prev, 0.975), 1))"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0 0.01002662169206058 0.018212148574401982\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PWKJZXwn8gVT",
        "colab_type": "text"
      },
      "source": [
        "Therefore, we estimate that the 95% confidence interval for the prevalence is (0, 1.82%)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cqOPCkxyT9pU",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 491
        },
        "outputId": "30fa9aab-1478-4e39-90e3-3153eb1a6542"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "plt.hist(boot_prev, bins = 30)"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([1.000e+00, 0.000e+00, 4.000e+00, 5.000e+00, 1.500e+01, 3.100e+01,\n",
              "        4.200e+01, 8.600e+01, 1.380e+02, 3.020e+02, 4.610e+02, 7.120e+02,\n",
              "        1.150e+03, 1.676e+03, 2.410e+03, 3.417e+03, 4.097e+03, 5.050e+03,\n",
              "        5.819e+03, 5.971e+03, 5.817e+03, 4.686e+03, 3.687e+03, 2.397e+03,\n",
              "        1.240e+03, 5.430e+02, 1.840e+02, 4.300e+01, 1.300e+01, 3.000e+00]),\n",
              " array([-0.01589747, -0.01450886, -0.01312025, -0.01173164, -0.01034303,\n",
              "        -0.00895442, -0.0075658 , -0.00617719, -0.00478858, -0.00339997,\n",
              "        -0.00201136, -0.00062275,  0.00076586,  0.00215447,  0.00354308,\n",
              "         0.00493169,  0.0063203 ,  0.00770891,  0.00909752,  0.01048613,\n",
              "         0.01187474,  0.01326335,  0.01465196,  0.01604057,  0.01742918,\n",
              "         0.01881779,  0.0202064 ,  0.02159501,  0.02298363,  0.02437224,\n",
              "         0.02576085]),\n",
              " <a list of 30 Patch objects>)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAARoUlEQVR4nO3df4xl5V3H8fen0FJTtbvAuuIucWm6saEmbXECNBpTi10WaLokthVjZIIkGyP+Skx0sSYolWTpH9Y2tTWbsro0KkW0YSNNcbpAjCb8GAqlAuJO6ZLdFdiRpShtSkP9+sc8U2/pDHNnd+bO7D7vV3Jzz/me55x5zsnu556c89xzU1VIkvrwmpXugCRpdAx9SeqIoS9JHTH0Jakjhr4kdeTUle7AqznzzDNr06ZNK90NSTqhPPjgg/9VVevmWraqQ3/Tpk1MTk6udDck6YSS5Kn5lnl5R5I6YuhLUkcMfUnqiKEvSR0x9CWpI0OFfpI1SW5L8u9JHk/yziSnJ5lIsr+9r21tk+TjSaaSPJLkvIHtjLf2+5OML9dOSZLmNuyZ/seAL1TVW4C3AY8DO4B9VbUZ2NfmAS4BNrfXduBTAElOB64DLgDOB66b/aCQJI3GgqGf5I3AzwI3AVTVt6vq68A2YE9rtge4vE1vA26uGfcCa5KcBVwMTFTV0ap6HpgAti7p3kiSXtUwZ/rnANPAXyZ5KMmnk7wBWF9VT7c2zwDr2/QG4ODA+odabb7690iyPclkksnp6enF7Y0k6VUN843cU4HzgN+sqvuSfIz/v5QDQFVVkiX5NZaq2gXsAhgbG/MXXqRjtGnHHUO1O7DzsmXuiVaTYc70DwGHquq+Nn8bMx8Cz7bLNrT3I235YeDsgfU3ttp8dUnSiCwY+lX1DHAwyU+00kXAY8BeYHYEzjhwe5veC1zZRvFcCLzQLgPdCWxJsrbdwN3SapKkERn2gWu/Cfx1ktcBTwJXMfOBcWuSq4GngA+2tp8HLgWmgG+2tlTV0SQfBh5o7a6vqqNLsheSpKEMFfpV9TAwNseii+ZoW8A182xnN7B7MR2UJC2dVf1oZUnfb9gbtMuxPW/6nvh8DIMkdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0JakjQ4V+kgNJvpLk4SSTrXZ6kokk+9v72lZPko8nmUrySJLzBrYz3trvTzK+PLskSZrPYs70f66q3l5VY21+B7CvqjYD+9o8wCXA5vbaDnwKZj4kgOuAC4DzgetmPygkSaNx6nGsuw14V5veA9wD/H6r31xVBdybZE2Ss1rbiao6CpBkAtgK/O1x9EE6aWzaccdKd0EdGPZMv4B/SvJgku2ttr6qnm7TzwDr2/QG4ODAuodabb7690iyPclkksnp6ekhuydJGsawZ/o/U1WHk/wIMJHk3wcXVlUlqaXoUFXtAnYBjI2NLck2JUkzhjrTr6rD7f0I8Dlmrsk/2y7b0N6PtOaHgbMHVt/YavPVJUkjsmDoJ3lDkh+anQa2AP8G7AVmR+CMA7e36b3AlW0Uz4XAC+0y0J3AliRr2w3cLa0mSRqRYS7vrAc+l2S2/d9U1ReSPADcmuRq4Cngg63954FLgSngm8BVAFV1NMmHgQdau+tnb+pKkkZjwdCvqieBt81Rfw64aI56AdfMs63dwO7Fd1OStBT8Rq4kdcTQl6SOGPqS1BFDX5I6YuhLUkeO59k7kjoz7POBDuy8bJl7omPlmb4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjoy9A+jJzkFmAQOV9V7k5wD3AKcATwI/EpVfTvJacDNwE8BzwG/WFUH2jauBa4GvgP8VlXduZQ7I61Gw/6YuDQKiznT/23g8YH5G4GPVtWbgeeZCXPa+/Ot/tHWjiTnAlcAbwW2Ap9sHySSpBEZKvSTbAQuAz7d5gO8G7itNdkDXN6mt7V52vKLWvttwC1V9VJVfQ2YAs5fip2QJA1n2DP9PwN+D/jfNn8G8PWqernNHwI2tOkNwEGAtvyF1v679TnW+a4k25NMJpmcnp5exK5IkhayYOgneS9wpKoeHEF/qKpdVTVWVWPr1q0bxZ+UpG4McyP3p4H3JbkUeD3ww8DHgDVJTm1n8xuBw639YeBs4FCSU4E3MnNDd7Y+a3AdSdIILHimX1XXVtXGqtrEzI3Yu6rql4G7gfe3ZuPA7W16b5unLb+rqqrVr0hyWhv5sxm4f8n2RJK0oKGHbM7h94FbkvwJ8BBwU6vfBHwmyRRwlJkPCqrq0SS3Ao8BLwPXVNV3juPvS5IWaVGhX1X3APe06SeZY/RNVX0L+MA8698A3LDYTkqSlobfyJWkjhj6ktQRQ1+SOmLoS1JHDH1J6sjxDNmUpDkN+2TRAzsvW+ae6JU805ekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xC9nScdo2C8gSauJZ/qS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHFgz9JK9Pcn+SLyd5NMkft/o5Se5LMpXks0le1+qntfmptnzTwLaubfUnkly8XDslSZrbMGf6LwHvrqq3AW8Htia5ELgR+GhVvRl4Hri6tb8aeL7VP9rakeRc4ArgrcBW4JNJTlnKnZEkvboFQ79mvNhmX9teBbwbuK3V9wCXt+ltbZ62/KIkafVbquqlqvoaMAWcvyR7IUkaylDX9JOckuRh4AgwAXwV+HpVvdyaHAI2tOkNwEGAtvwF4IzB+hzrDP6t7Ukmk0xOT08vfo8kSfMaKvSr6jtV9XZgIzNn529Zrg5V1a6qGquqsXXr1i3Xn5GkLi1q9E5VfR24G3gnsCbJ7I+wbAQOt+nDwNkAbfkbgecG63OsI0kagWFG76xLsqZN/wDwHuBxZsL//a3ZOHB7m97b5mnL76qqavUr2uiec4DNwP1LtSOSpIUN83OJZwF72kib1wC3VtU/JnkMuCXJnwAPATe19jcBn0kyBRxlZsQOVfVokluBx4CXgWuq6jtLuzuSpFezYOhX1SPAO+aoP8kco2+q6lvAB+bZ1g3ADYvvpiRpKfiNXEnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUkWGepy91ZdOOO1a6C9KyMfQlrZhhP2AP7LxsmXvSDy/vSFJHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRxYM/SRnJ7k7yWNJHk3y261+epKJJPvb+9pWT5KPJ5lK8kiS8wa2Nd7a708yvny7JUmayzBn+i8Dv1tV5wIXAtckORfYAeyrqs3AvjYPcAmwub22A5+CmQ8J4DrgAuB84LrZDwpJ0mgsGPpV9XRVfalN/w/wOLAB2Absac32AJe36W3AzTXjXmBNkrOAi4GJqjpaVc8DE8DWJd0bSdKrWtQ1/SSbgHcA9wHrq+rptugZYH2b3gAcHFjtUKvNV3/l39ieZDLJ5PT09GK6J0lawNChn+QHgb8Hfqeq/ntwWVUVUEvRoaraVVVjVTW2bt26pdikJKkZKvSTvJaZwP/rqvqHVn62XbahvR9p9cPA2QOrb2y1+eqSpBEZZvROgJuAx6vqTwcW7QVmR+CMA7cP1K9so3guBF5ol4HuBLYkWdtu4G5pNUnSiAzzIyo/DfwK8JUkD7faHwA7gVuTXA08BXywLfs8cCkwBXwTuAqgqo4m+TDwQGt3fVUdXZK9kCQNZcHQr6p/ATLP4ovmaF/ANfNsazewezEdlCQtHb+RK0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjoyzDh96YS3accdK90FaVXwTF+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcRHK0ta9YZ9NPaBnZctc09OfJ7pS1JHDH1J6oiXd3RC8xexpMXxTF+SOrJg6CfZneRIkn8bqJ2eZCLJ/va+ttWT5ONJppI8kuS8gXXGW/v9ScaXZ3ckSa9mmDP9vwK2vqK2A9hXVZuBfW0e4BJgc3ttBz4FMx8SwHXABcD5wHWzHxSSpNFZMPSr6p+Bo68obwP2tOk9wOUD9Ztrxr3AmiRnARcDE1V1tKqeByb4/g8SSdIyO9Zr+uur6uk2/Qywvk1vAA4OtDvUavPVv0+S7Ukmk0xOT08fY/ckSXM57hu5VVVALUFfZre3q6rGqmps3bp1S7VZSRLHHvrPtss2tPcjrX4YOHug3cZWm68uSRqhYw39vcDsCJxx4PaB+pVtFM+FwAvtMtCdwJYka9sN3C2tJkkaoQW/nJXkb4F3AWcmOcTMKJydwK1JrgaeAj7Ymn8euBSYAr4JXAVQVUeTfBh4oLW7vqpeeXNYkrTMFgz9qvqleRZdNEfbAq6ZZzu7gd2L6p0kaUn5jVxJ6oihL0kdMfQlqSOGviR1xEcra1XykcnS8vBMX5I6YuhLUkcMfUnqiKEvSR0x9CWpI47ekXTSGHbU14Gdly1zT1Yvz/QlqSOGviR1xMs7Gim/dCWtLM/0Jakjhr4kdcTQl6SOGPqS1BFv5GpJeINWOjF4pi9JHTH0Jakjhr4kdcRr+pK60/MzejzTl6SOeKavV+WoHOnk4pm+JHVk5Gf6SbYCHwNOAT5dVTtH3Qd5Bi/1aqShn+QU4M+B9wCHgAeS7K2qx0bZj5OVQS4trcX8nzpRbvqO+kz/fGCqqp4ESHILsA3oMvQNaenkcaKMCBp16G8ADg7MHwIuGGyQZDuwvc2+mOSJEfVtpZ0J/NdKd2KV8xgtzGO0sBU9RrlxJH/mx+dbsOpG71TVLmDXSvdj1JJMVtXYSvdjNfMYLcxjtLDej9GoR+8cBs4emN/YapKkERh16D8AbE5yTpLXAVcAe0fcB0nq1kgv71TVy0l+A7iTmSGbu6vq0VH2YRXr7pLWMfAYLcxjtLCuj1GqaqX7IEkaEb+RK0kdMfQlqSOG/ogkOT3JRJL97X3tPO3GW5v9ScYH6jckOZjkxdH1ejSSbE3yRJKpJDvmWH5aks+25fcl2TSw7NpWfyLJxaPs9ygd6zFKckaSu5O8mOQTo+73KB3HMXpPkgeTfKW9v3vUfR+pqvI1ghfwEWBHm94B3DhHm9OBJ9v72ja9ti27EDgLeHGl92WJj8spwFeBNwGvA74MnPuKNr8O/EWbvgL4bJs+t7U/DTinbeeUld6nVXaM3gD8DPBrwCdWel9W6TF6B/BjbfongcMrvT/L+fJMf3S2AXva9B7g8jnaXAxMVNXRqnoemAC2AlTVvVX19Eh6OlrffTRHVX0bmH00x6DBY3cbcFGStPotVfVSVX0NmGrbO9kc8zGqqm9U1b8A3xpdd1fE8Ryjh6rqP1v9UeAHkpw2kl6vAEN/dNYPhPYzwPo52sz1mIoNy92xFTbMPn+3TVW9DLwAnDHkuieD4zlGvViqY/QLwJeq6qVl6ueKW3WPYTiRJfki8KNzLPrQ4ExVVRLHykqrSJK3AjcCW1a6L8vJ0F9CVfXz8y1L8mySs6rq6SRnAUfmaHYYeNfA/EbgniXt5OozzKM5ZtscSnIq8EbguSHXPRkczzHqxXEdoyQbgc8BV1bVV5e/uyvHyzujsxeYHY0zDtw+R5s7gS1J1rbRPVta7WQ2zKM5Bo/d+4G7auau217gijYq4xxgM3D/iPo9SsdzjHpxzMcoyRrgDmYGWvzryHq8Ulb6TnIvL2auHe4D9gNfBE5v9TFmfkFstt2vMnNDcgq4aqD+EWauU/5ve/+jld6nJTw2lwL/wczoiw+12vXA+9r064G/a8fkfuBNA+t+qK33BHDJSu/LKj1GB4CjwIvt3865o+7/aj5GwB8C3wAeHnj9yErvz3K9fAyDJHXEyzuS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXk/wAJOIK3s2NRowAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}
	{
	"nbformat": 4,
	"nbformat_minor": 0,
	"metadata": {
	"colab": {
	"name": "Confidence Intervals for SARS-CoV-2 Antibody Prevalence Using the Parametric Bootstrap ",
	"provenance": [],
	"collapsed_sections": [],
	"authorship_tag": "ABX9TyOvp9xLSzNWEYZ46gjqS2OG",
	"include_colab_link": true
	},
	"kernelspec": {
	"name": "python3",
	"display_name": "Python 3"
	}
	},
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "view-in-github",
	"colab_type": "text"
	},
	"source": [
	"<a href=\"https://colab.research.google.com/gist/syadlowsky/a23a2f67358ef64a7d80d8c5cc0ad307/confidence-intervals-for-sars-cov-2-antibody-prevalence-using-the-parametric-bootstrap.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "fDlzRj8g8s3d",
	"colab_type": "text"
	},
	"source": [
	"# Confidence Intervals for SARS-CoV-2 Antibody Prevalence Using the Parametric Bootstrap "
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "R6P878uLC6k0",
	"colab_type": "text"
	},
	"source": [
	"Here are the values reported in the study."
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "TaSXJMnloknz",
	"colab_type": "code",
	"colab": {}
	},
	"source": [
	"# manufacturer validation data\n",
	"manufacturer_neg_tests = 371\n",
	"manufacturer_false_positives = 2\n",
	"manufacturer_pos_tests = 85\n",
	"manufacturer_true_pos = 78\n",
	"\n",
	"# currently not used\n",
	"stanford_val_neg_tests = 30\n",
	"stanford_val_false_positives = 0\n",
	"\n",
	"# community test results\n",
	"study_tests = 3439\n",
	"study_positives = 50"
	],
	"execution_count": 0,
	"outputs": []
	},
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "0IMIlvJp8-ct",
	"colab_type": "text"
	},
	"source": [
	"We can calculate the false positive rate (FPR), the true positive rate (TPR), and an estimate of the prevalence using the estimating equation, $$ (1-\\text{Prevalence})\\cdot \\text{FPR} + \\text{Prevalence} \\cdot \\text{TPR} = \\text{Observed Positive Rate}.$$"
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "szE0nH2s3XYz",
	"colab_type": "code",
	"colab": {}
	},
	"source": [
	"import numpy as np\n",
	"\n",
	"def calc_parameters(val_neg_tests, val_false_positives, val_pos_tests,\n",
	" val_true_positives, study_tests, study_positives):\n",
	" fpr = val_false_positives / val_neg_tests\n",
	" tpr = val_true_positives / val_pos_tests\n",
	" observed_positive_fraction = study_positives / study_tests\n",
	"\n",
	" prev = (observed_positive_fraction - fpr) / (tpr - fpr)\n",
	" return fpr, tpr, prev"
	],
	"execution_count": 0,
	"outputs": []
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "SdeNu68v5H5-",
	"colab_type": "code",
	"outputId": "375a027c-9b48-4cba-e4a6-92f97042ee8a",
	"colab": {
	"base_uri": "https://localhost:8080/",
	"height": 34
	}
	},
	"source": [
	"fpr, tpr, prev = calc_parameters(manufacturer_neg_tests, manufacturer_false_positives,\n",
	" manufacturer_pos_tests, manufacturer_true_pos,\n",
	" study_tests, study_positives)\n",
	"print(fpr, tpr, prev)"
	],
	"execution_count": 3,
	"outputs": [
	{
	"output_type": "stream",
	"text": [
	"0.005390835579514825 0.9176470588235294 0.010028185496404683\n"
	],
	"name": "stdout"
	}
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "nznlCkRn9zHl",
	"colab_type": "text"
	},
	"source": [
	"The estimated FPR is 0.5%, the estimated TPR is 91.7%, and the estimated prevalence is 1.00%.\n",
	"\n",
	"Now, using these estimates and pretending like they are the true parameters, let's simulate the data again, using the same number of true negatives, true positives, and community members tested. Then, we aggregate these simulations and compute the mean and quantiles to get the parametric bootstrap estimate of the confidence intervals."
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "0AeKfSOa6CpZ",
	"colab_type": "code",
	"outputId": "a8e12d25-87fd-40a2-d437-18b7e2dcf2a1",
	"colab": {
	"base_uri": "https://localhost:8080/",
	"height": 34
	}
	},
	"source": [
	"boot_prev = list()\n",
	"for b in range(50000):\n",
	" boot_false_positives = np.random.binomial(manufacturer_neg_tests, fpr)\n",
	" boot_true_pos = np.random.binomial(manufacturer_pos_tests, tpr)\n",
	" boot_study_has_antibody = np.random.binomial(study_tests, prev)\n",
	" boot_study_pos = np.random.binomial(boot_study_has_antibody, tpr) + \\\n",
	" np.random.binomial(study_tests - boot_study_has_antibody, fpr)\n",
	" _, _, boot_sample = calc_parameters(manufacturer_neg_tests, boot_false_positives,\n",
	" manufacturer_pos_tests, boot_true_pos,\n",
	" study_tests, boot_study_pos)\n",
	" boot_prev.append(boot_sample)\n",
	"print(max(np.quantile(boot_prev, 0.025), 0), np.mean(boot_prev), min(np.quantile(boot_prev, 0.975), 1))"
	],
	"execution_count": 4,
	"outputs": [
	{
	"output_type": "stream",
	"text": [
	"0 0.01002662169206058 0.018212148574401982\n"
	],
	"name": "stdout"
	}
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "PWKJZXwn8gVT",
	"colab_type": "text"
	},
	"source": [
	"Therefore, we estimate that the 95% confidence interval for the prevalence is (0, 1.82%)."
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "cqOPCkxyT9pU",
	"colab_type": "code",
	"colab": {
	"base_uri": "https://localhost:8080/",
	"height": 491
	},
	"outputId": "30fa9aab-1478-4e39-90e3-3153eb1a6542"
	},
	"source": [
	"import matplotlib.pyplot as plt\n",
	"plt.hist(boot_prev, bins = 30)"
	],
	"execution_count": 5,
	"outputs": [
	{
	"output_type": "execute_result",
	"data": {
	"text/plain": [
	"(array([1.000e+00, 0.000e+00, 4.000e+00, 5.000e+00, 1.500e+01, 3.100e+01,\n",
	" 4.200e+01, 8.600e+01, 1.380e+02, 3.020e+02, 4.610e+02, 7.120e+02,\n",
	" 1.150e+03, 1.676e+03, 2.410e+03, 3.417e+03, 4.097e+03, 5.050e+03,\n",
	" 5.819e+03, 5.971e+03, 5.817e+03, 4.686e+03, 3.687e+03, 2.397e+03,\n",
	" 1.240e+03, 5.430e+02, 1.840e+02, 4.300e+01, 1.300e+01, 3.000e+00]),\n",
	" array([-0.01589747, -0.01450886, -0.01312025, -0.01173164, -0.01034303,\n",
	" -0.00895442, -0.0075658 , -0.00617719, -0.00478858, -0.00339997,\n",
	" -0.00201136, -0.00062275, 0.00076586, 0.00215447, 0.00354308,\n",
	" 0.00493169, 0.0063203 , 0.00770891, 0.00909752, 0.01048613,\n",
	" 0.01187474, 0.01326335, 0.01465196, 0.01604057, 0.01742918,\n",
	" 0.01881779, 0.0202064 , 0.02159501, 0.02298363, 0.02437224,\n",
	" 0.02576085]),\n",
	" <a list of 30 Patch objects>)"
	]
	},
	"metadata": {
	"tags": []
	},
	"execution_count": 5
	},
	{
	"output_type": "display_data",
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAARoUlEQVR4nO3df4xl5V3H8fen0FJTtbvAuuIucWm6saEmbXECNBpTi10WaLokthVjZIIkGyP+Skx0sSYolWTpH9Y2tTWbsro0KkW0YSNNcbpAjCb8GAqlAuJO6ZLdFdiRpShtSkP9+sc8U2/pDHNnd+bO7D7vV3Jzz/me55x5zsnu556c89xzU1VIkvrwmpXugCRpdAx9SeqIoS9JHTH0Jakjhr4kdeTUle7AqznzzDNr06ZNK90NSTqhPPjgg/9VVevmWraqQ3/Tpk1MTk6udDck6YSS5Kn5lnl5R5I6YuhLUkcMfUnqiKEvSR0x9CWpI0OFfpI1SW5L8u9JHk/yziSnJ5lIsr+9r21tk+TjSaaSPJLkvIHtjLf2+5OML9dOSZLmNuyZ/seAL1TVW4C3AY8DO4B9VbUZ2NfmAS4BNrfXduBTAElOB64DLgDOB66b/aCQJI3GgqGf5I3AzwI3AVTVt6vq68A2YE9rtge4vE1vA26uGfcCa5KcBVwMTFTV0ap6HpgAti7p3kiSXtUwZ/rnANPAXyZ5KMmnk7wBWF9VT7c2zwDr2/QG4ODA+odabb7690iyPclkksnp6enF7Y0k6VUN843cU4HzgN+sqvuSfIz/v5QDQFVVkiX5NZaq2gXsAhgbG/MXXqRjtGnHHUO1O7DzsmXuiVaTYc70DwGHquq+Nn8bMx8Cz7bLNrT3I235YeDsgfU3ttp8dUnSiCwY+lX1DHAwyU+00kXAY8BeYHYEzjhwe5veC1zZRvFcCLzQLgPdCWxJsrbdwN3SapKkERn2gWu/Cfx1ktcBTwJXMfOBcWuSq4GngA+2tp8HLgWmgG+2tlTV0SQfBh5o7a6vqqNLsheSpKEMFfpV9TAwNseii+ZoW8A182xnN7B7MR2UJC2dVf1oZUnfb9gbtMuxPW/6nvh8DIMkdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0JakjQ4V+kgNJvpLk4SSTrXZ6kokk+9v72lZPko8nmUrySJLzBrYz3trvTzK+PLskSZrPYs70f66q3l5VY21+B7CvqjYD+9o8wCXA5vbaDnwKZj4kgOuAC4DzgetmPygkSaNx6nGsuw14V5veA9wD/H6r31xVBdybZE2Ss1rbiao6CpBkAtgK/O1x9EE6aWzaccdKd0EdGPZMv4B/SvJgku2ttr6qnm7TzwDr2/QG4ODAuodabb7690iyPclkksnp6ekhuydJGsawZ/o/U1WHk/wIMJHk3wcXVlUlqaXoUFXtAnYBjI2NLck2JUkzhjrTr6rD7f0I8Dlmrsk/2y7b0N6PtOaHgbMHVt/YavPVJUkjsmDoJ3lDkh+anQa2AP8G7AVmR+CMA7e36b3AlW0Uz4XAC+0y0J3AliRr2w3cLa0mSRqRYS7vrAc+l2S2/d9U1ReSPADcmuRq4Cngg63954FLgSngm8BVAFV1NMmHgQdau+tnb+pKkkZjwdCvqieBt81Rfw64aI56AdfMs63dwO7Fd1OStBT8Rq4kdcTQl6SOGPqS1BFDX5I6YuhLUkeO59k7kjoz7POBDuy8bJl7omPlmb4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjoy9A+jJzkFmAQOV9V7k5wD3AKcATwI/EpVfTvJacDNwE8BzwG/WFUH2jauBa4GvgP8VlXduZQ7I61Gw/6YuDQKiznT/23g8YH5G4GPVtWbgeeZCXPa+/Ot/tHWjiTnAlcAbwW2Ap9sHySSpBEZKvSTbAQuAz7d5gO8G7itNdkDXN6mt7V52vKLWvttwC1V9VJVfQ2YAs5fip2QJA1n2DP9PwN+D/jfNn8G8PWqernNHwI2tOkNwEGAtvyF1v679TnW+a4k25NMJpmcnp5exK5IkhayYOgneS9wpKoeHEF/qKpdVTVWVWPr1q0bxZ+UpG4McyP3p4H3JbkUeD3ww8DHgDVJTm1n8xuBw639YeBs4FCSU4E3MnNDd7Y+a3AdSdIILHimX1XXVtXGqtrEzI3Yu6rql4G7gfe3ZuPA7W16b5unLb+rqqrVr0hyWhv5sxm4f8n2RJK0oKGHbM7h94FbkvwJ8BBwU6vfBHwmyRRwlJkPCqrq0SS3Ao8BLwPXVNV3juPvS5IWaVGhX1X3APe06SeZY/RNVX0L+MA8698A3LDYTkqSlobfyJWkjhj6ktQRQ1+SOmLoS1JHDH1J6sjxDNmUpDkN+2TRAzsvW+ae6JU805ekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xC9nScdo2C8gSauJZ/qS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHFgz9JK9Pcn+SLyd5NMkft/o5Se5LMpXks0le1+qntfmptnzTwLaubfUnkly8XDslSZrbMGf6LwHvrqq3AW8Htia5ELgR+GhVvRl4Hri6tb8aeL7VP9rakeRc4ArgrcBW4JNJTlnKnZEkvboFQ79mvNhmX9teBbwbuK3V9wCXt+ltbZ62/KIkafVbquqlqvoaMAWcvyR7IUkaylDX9JOckuRh4AgwAXwV+HpVvdyaHAI2tOkNwEGAtvwF4IzB+hzrDP6t7Ukmk0xOT08vfo8kSfMaKvSr6jtV9XZgIzNn529Zrg5V1a6qGquqsXXr1i3Xn5GkLi1q9E5VfR24G3gnsCbJ7I+wbAQOt+nDwNkAbfkbgecG63OsI0kagWFG76xLsqZN/wDwHuBxZsL//a3ZOHB7m97b5mnL76qqavUr2uiec4DNwP1LtSOSpIUN83OJZwF72kib1wC3VtU/JnkMuCXJnwAPATe19jcBn0kyBRxlZsQOVfVokluBx4CXgWuq6jtLuzuSpFezYOhX1SPAO+aoP8kco2+q6lvAB+bZ1g3ADYvvpiRpKfiNXEnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUkWGepy91ZdOOO1a6C9KyMfQlrZhhP2AP7LxsmXvSDy/vSFJHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRxYM/SRnJ7k7yWNJHk3y261+epKJJPvb+9pWT5KPJ5lK8kiS8wa2Nd7a708yvny7JUmayzBn+i8Dv1tV5wIXAtckORfYAeyrqs3AvjYPcAmwub22A5+CmQ8J4DrgAuB84LrZDwpJ0mgsGPpV9XRVfalN/w/wOLAB2Absac32AJe36W3AzTXjXmBNkrOAi4GJqjpaVc8DE8DWJd0bSdKrWtQ1/SSbgHcA9wHrq+rptugZYH2b3gAcHFjtUKvNV3/l39ieZDLJ5PT09GK6J0lawNChn+QHgb8Hfqeq/ntwWVUVUEvRoaraVVVjVTW2bt26pdikJKkZKvSTvJaZwP/rqvqHVn62XbahvR9p9cPA2QOrb2y1+eqSpBEZZvROgJuAx6vqTwcW7QVmR+CMA7cP1K9so3guBF5ol4HuBLYkWdtu4G5pNUnSiAzzIyo/DfwK8JUkD7faHwA7gVuTXA08BXywLfs8cCkwBXwTuAqgqo4m+TDwQGt3fVUdXZK9kCQNZcHQr6p/ATLP4ovmaF/ANfNsazewezEdlCQtHb+RK0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjoyzDh96YS3accdK90FaVXwTF+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcRHK0ta9YZ9NPaBnZctc09OfJ7pS1JHDH1J6oiXd3RC8xexpMXxTF+SOrJg6CfZneRIkn8bqJ2eZCLJ/va+ttWT5ONJppI8kuS8gXXGW/v9ScaXZ3ckSa9mmDP9vwK2vqK2A9hXVZuBfW0e4BJgc3ttBz4FMx8SwHXABcD5wHWzHxSSpNFZMPSr6p+Bo68obwP2tOk9wOUD9Ztrxr3AmiRnARcDE1V1tKqeByb4/g8SSdIyO9Zr+uur6uk2/Qywvk1vAA4OtDvUavPVv0+S7Ukmk0xOT08fY/ckSXM57hu5VVVALUFfZre3q6rGqmps3bp1S7VZSRLHHvrPtss2tPcjrX4YOHug3cZWm68uSRqhYw39vcDsCJxx4PaB+pVtFM+FwAvtMtCdwJYka9sN3C2tJkkaoQW/nJXkb4F3AWcmOcTMKJydwK1JrgaeAj7Ymn8euBSYAr4JXAVQVUeTfBh4oLW7vqpeeXNYkrTMFgz9qvqleRZdNEfbAq6ZZzu7gd2L6p0kaUn5jVxJ6oihL0kdMfQlqSOGviR1xEcra1XykcnS8vBMX5I6YuhLUkcMfUnqiKEvSR0x9CWpI47ekXTSGHbU14Gdly1zT1Yvz/QlqSOGviR1xMs7Gim/dCWtLM/0Jakjhr4kdcTQl6SOGPqS1BFv5GpJeINWOjF4pi9JHTH0Jakjhr4kdcRr+pK60/MzejzTl6SOeKavV+WoHOnk4pm+JHVk5Gf6SbYCHwNOAT5dVTtH3Qd5Bi/1aqShn+QU4M+B9wCHgAeS7K2qx0bZj5OVQS4trcX8nzpRbvqO+kz/fGCqqp4ESHILsA3oMvQNaenkcaKMCBp16G8ADg7MHwIuGGyQZDuwvc2+mOSJEfVtpZ0J/NdKd2KV8xgtzGO0sBU9RrlxJH/mx+dbsOpG71TVLmDXSvdj1JJMVtXYSvdjNfMYLcxjtLDej9GoR+8cBs4emN/YapKkERh16D8AbE5yTpLXAVcAe0fcB0nq1kgv71TVy0l+A7iTmSGbu6vq0VH2YRXr7pLWMfAYLcxjtLCuj1GqaqX7IEkaEb+RK0kdMfQlqSOG/ogkOT3JRJL97X3tPO3GW5v9ScYH6jckOZjkxdH1ejSSbE3yRJKpJDvmWH5aks+25fcl2TSw7NpWfyLJxaPs9ygd6zFKckaSu5O8mOQTo+73KB3HMXpPkgeTfKW9v3vUfR+pqvI1ghfwEWBHm94B3DhHm9OBJ9v72ja9ti27EDgLeHGl92WJj8spwFeBNwGvA74MnPuKNr8O/EWbvgL4bJs+t7U/DTinbeeUld6nVXaM3gD8DPBrwCdWel9W6TF6B/BjbfongcMrvT/L+fJMf3S2AXva9B7g8jnaXAxMVNXRqnoemAC2AlTVvVX19Eh6OlrffTRHVX0bmH00x6DBY3cbcFGStPotVfVSVX0NmGrbO9kc8zGqqm9U1b8A3xpdd1fE8Ryjh6rqP1v9UeAHkpw2kl6vAEN/dNYPhPYzwPo52sz1mIoNy92xFTbMPn+3TVW9DLwAnDHkuieD4zlGvViqY/QLwJeq6qVl6ueKW3WPYTiRJfki8KNzLPrQ4ExVVRLHykqrSJK3AjcCW1a6L8vJ0F9CVfXz8y1L8mySs6rq6SRnAUfmaHYYeNfA/EbgniXt5OozzKM5ZtscSnIq8EbguSHXPRkczzHqxXEdoyQbgc8BV1bVV5e/uyvHyzujsxeYHY0zDtw+R5s7gS1J1rbRPVta7WQ2zKM5Bo/d+4G7auau217gijYq4xxgM3D/iPo9SsdzjHpxzMcoyRrgDmYGWvzryHq8Ulb6TnIvL2auHe4D9gNfBE5v9TFmfkFstt2vMnNDcgq4aqD+EWauU/5ve/+jld6nJTw2lwL/wczoiw+12vXA+9r064G/a8fkfuBNA+t+qK33BHDJSu/LKj1GB4CjwIvt3865o+7/aj5GwB8C3wAeHnj9yErvz3K9fAyDJHXEyzuS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXk/wAJOIK3s2NRowAAAABJRU5ErkJggg==\n",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"tags": [],
	"needs_background": "light"
	}
	}
	]
	}
	]
	}