Updated with Jupyter notebook with methods from Blood manuscript

IWP_from_OperationalData.ipynb

@@ -0,0 +1,175 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Estimating the infectious window period from operational data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The simple algorithm in this notebook is based on the idea that you can estimate the duration of the infectious window period (or window of residual risk) from lookback investigation data.\n",
+    "\n",
+    "Lookback investigations are triggered when repeat donors seroconvert and trigger an investigation into transfused material collected at the donor's previous donation. The prior donation could have been collected while the donor was already infected (and infectious) but not yet detectable.\n",
+    "\n",
+    "Imagine the following situation: \n",
+    "* 100 donors seroconvert and become HIV-positive\n",
+    "* all 100 had prior donations exactly 105 days before the donation at which HIV was detected\n",
+    "* 1 recipient of material from prior donations became infected and it is confirmed that transfusion transmission occurred\n",
+    "* the least sensitive positive test at the seroconversion donation has a diagnostic delay of 10 days\n",
+    "* the most sensitive negative test at the prior donation has a diagnostic delay of 5 days\n",
+    "\n",
+    "Each of these IDIs of 105 days should be adjusted by the diagnostic delays relevant to the time of last negative and first positive donations to derive the interval during which infection could have occurred. In this case to 105+5-10 = 100 days.\n",
+    "\n",
+    "In this simple example, the infectious window period associated with the screening strategy in use would be 1% of the duration of the shared infection interval or 1 day. The 1% comes from the 1/100 donors who transmitted the infection. This can be generalized to a situation where each donor has a different IDI by considering the number of transmissions as a Poisson process and treating the adjusted IDIs as \"inverse exposure time\".\n",
+    "\n",
+    "$$IWP = \\frac{n}{\\sum_{i = 1}^{N} \\frac{1}{IDI_{i}} }$$\n",
+    "\n",
+    "with $n$ the number of transmissions, $N$ the number of seroconverting donors and $IDI_i$ the i<sup>th</sup> donor's adjusted interdonation interval.\n",
+    "\n",
+    "The properties of the Poisson and Chi-square distributions then allow us to obtain confidence intervals on the IWP estimate:\n",
+    "\n",
+    "$$IWP_{lb} = \\frac{ \\left. \\chi_{2n,\\alpha/2}^{2} \\middle/ 2 \\right. }{ \\sum_{i = 1}^{N} \\frac{1}{IDI_{i}} }$$\n",
+    "\n",
+    "$$IWP_{ub} = \\frac{ \\left. \\chi_{2n,1-\\alpha/2}^{2} \\middle/ 2 \\right. }{ \\sum_{i = 1}^{N} \\frac{1}{IDI_{i}} }$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import scipy.stats as stats"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def residual_risk_operational(n_transmissions,\n",
+    "                              intervals, \n",
+    "                              negative_diagnositc_day,\n",
+    "                              positive_diagnostic_delay,\n",
+    "                              alpha = 0.05\n",
+    "                             ):\n",
+    "    inverse_intervals = [1/(x + negative_diagnositc_day - positive_diagnostic_delay) for x in intervals]\n",
+    "    total_exposure = sum(inverse_intervals)\n",
+    "    iwp = n_transmissions / total_exposure\n",
+    "    if n_transmissions > 0:\n",
+    "        iwp_lb = stats.chi2.ppf(alpha/2, df=2*n_transmissions)/2/total_exposure\n",
+    "    else:\n",
+    "        iwp_lb = 0.0\n",
+    "    iwp_ub = stats.chi2.ppf(1.0-alpha/2, df=2*(n_transmissions+1))/2/total_exposure\n",
+    "    return (iwp, (iwp_lb, iwp_ub))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example with equal IDIs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(0.9999999999999993, (0.02531780798428986, 5.5716433909388945))"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "IDIs = np.repeat(105,100)\n",
+    "residual_risk_operational(1, IDIs, 5, 10)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example with more realistic IDIs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 295,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def truncated_normal(mu=0, sigma=1, lower=-3, upper=3):\n",
+    "    return stats.truncnorm(a = (lower - mu)/sigma,\n",
+    "                           b = (upper - mu)/sigma,\n",
+    "                           loc = mu,\n",
+    "                           scale = sigma\n",
+    "                          )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 298,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(123)\n",
+    "X = truncated_normal(mu = 0, sigma = 600, lower = 56, upper = 2000)\n",
+    "IDIs = X.rvs(500)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 299,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(2.733677254711141, (0.887617563599965, 6.379490801132148))"
+      ]
+     },
+     "execution_count": 299,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "residual_risk_operational(5, IDIs, 5, 10)"
+   ]
+  }
+ ],
+ "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.8.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}