Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Confidence Intervals for SARS-CoV-2 Antibody Prevalence Using the Parametric Bootstrap
{
"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"
}
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment