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
You can’t perform that action at this time.