Adding section on CI of betas

README.md

@@ -1,24 +1,28 @@
-# Uncertainty study
+# Uncertainty-aware quantitative risk assessment (QRA)
 
+This project studies uncertainty in estimates of the probability of adverse events based on the SUVmax biomarker. The main focus is polynomial logistic regression, including monotonicity constraints, regularization, goodness-of-fit assessment, and asymptotic and bootstrap uncertainty estimates.
 
-Studing modeling of uncertainty in estimating probability of adverse events based on SUV max biomarker
+## Repository structure
 
-Team:
+- `src/irae_risk/`: core implementation of the logistic-regression models and supporting numerical methods.
+- `notebooks/`: analyses used to generate and review the results.
+  - `logit_all_models.ipynb`: comparison of the fitted logistic-regression models.
+  - `logit_mono-cubic4paper.ipynb`: monotonic cubic analysis used for the manuscript.
+  - `logit_review_boots.ipynb`: review of bootstrap-based uncertainty estimates.
+- `tests/`: automated tests for the core implementation.
+- `data/`: raw, interim, and processed data.
+- `results/`: generated tables, figures, and cached calculations.
+- `manuscript/`: manuscript sources and related material.
 
-  * Marija Delić <makidelic@gmail.com>  PI
-  * Strašek, Katja <Katja.Strasek@fmf.uni-lj.si>
-  * Martin Horvat <martin.horvat@fmf.uni-lj.si>
-  * Robert Jeraj <robert.jeraj@fmf.uni-lj.si>
+## Team
 
-Other online resources:
+- Marija Delić <makidelic@gmail.com>
+- Zahra Alirezaei <Zahra.Alirezaei@fmf.uni-lj.si>
+- Katja Strašek <Katja.Strasek@fmf.uni-lj.si>
+- Martin Horvat <martin.horvat@fmf.uni-lj.si>
+- Robert Jeraj <robert.jeraj@fmf.uni-lj.si>
 
-  A. Minutes of the meetings:
+## Online resources
 
-    https://docs.google.com/document/d/10vQ1hQK9TbRgyw8s8O_KXbIIqJ8lp8AsJKMbpaviFuA/edit?usp=sharing
-
-
-  B. Repository
-
-    https://med1.fmf.uni-lj.si/owncloud/index.php/apps/files/?dir=/Optimisation%20Group&fileid=689823
-
-  prepared by Strašek, Katja <Katja.Strasek@fmf.uni-lj.si>. Currently this repo is not available :( .
+- [Meeting minutes](https://docs.google.com/document/d/10vQ1hQK9TbRgyw8s8O_KXbIIqJ8lp8AsJKMbpaviFuA/edit?usp=sharing)
+- [Original project repository](https://med1.fmf.uni-lj.si/owncloud/index.php/apps/files/?dir=/Optimisation%20Group&fileid=689823), prepared by Katja Strašek. This repository is currently unavailable.

notebooks/logit_review_boots.ipynb

@@ -13,8 +13,16 @@
    ]
   },
   {
+   "cell_type": "markdown",
+   "id": "7b2cffa0",
+   "metadata": {},
+   "source": [
+    "## Common"
+   ]
+  },
+  {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -25,8 +33,16 @@
     "project_root = Path.cwd().parent\n",
     "src_dir = project_root / 'src'\n",
     "if str(src_dir) not in sys.path:\n",
-    "    sys.path.insert(0, str(src_dir))\n",
-    "\n",
+    "    sys.path.insert(0, str(src_dir))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "063c7f0a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "import pandas as pd\n",
@@ -910,7 +926,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1026,23 +1042,23 @@
     "if not res_fit['success']:\n",
     "    raise RuntimeError(f\"Logistic fit failed: {res_fit}\")\n",
     "\n",
-    "pars = res_fit['pars']\n",
-    "cov_pars = logit.get_cov(x, y, pars)\n",
+    "theta = res_fit['theta']\n",
+    "cov_theta = logit.get_cov(x, y, theta)\n",
     "alpha = 0.05\n",
     "probs = [alpha / 2, 1 - alpha / 2]\n",
     "\n",
     "print(res_fit)\n",
-    "print('covariance matrix:\\n', cov_pars)\n",
+    "print('covariance matrix:\\n', cov_theta)\n",
     "pd.DataFrame(\n",
-    "    {'estimate': pars, 'LCL': logit.get_pars_quantiles_normal(probs, pars, cov_pars)[0],\n",
-    "     'UCL': logit.get_pars_quantiles_normal(probs, pars, cov_pars)[1]},\n",
-    "    index=[f'b{i}' for i in range(len(pars))],\n",
+    "    {'estimate': theta, 'LCL': logit.get_theta_quantiles_normal(probs, theta, cov_theta)[0],\n",
+    "     'UCL': logit.get_theta_quantiles_normal(probs, theta, cov_theta)[1]},\n",
+    "    index=[f'b{i}' for i in range(len(theta))],\n",
     ")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1181,12 +1197,12 @@
     }
    ],
    "source": [
-    "pd.DataFrame([logit.goodness_of_fit(x, y, pars)])"
+    "pd.DataFrame([logit.goodness_of_fit(x, y, theta)])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1207,11 +1223,11 @@
     "for cls, label in [(0, 'NC'), (1, 'AE')]:\n",
     "    ax.scatter(x[y == cls], y[y == cls], label=label, zorder=3)\n",
     "\n",
-    "ax.plot(xp, logit.model(xp, pars), color='black', label='logistic fit', zorder=4)\n",
+    "ax.plot(xp, logit.model(xp, theta), color='black', label='logistic fit', zorder=4)\n",
     "\n",
     "for label, method, color in [('normal', logit.get_model_quantiles_normal, 'tab:red'),\n",
     "                             ('delta', logit.get_model_quantiles_delta, 'tab:green')]:\n",
-    "    lo, hi = method(xp, probs, pars, cov_pars)\n",
+    "    lo, hi = method(xp, probs, theta, cov_theta)\n",
     "    ax.fill_between(xp, lo, hi, color=color, alpha=0.18, label=f'95% CI ({label})')\n",
     "\n",
     "ax.set(xlabel='maximum SUV percentile across visits', \n",
@@ -1233,7 +1249,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1257,7 +1273,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1288,9 +1304,9 @@
     "    return samples.mean(axis=0), np.cov(samples, rowvar=False, ddof=1)\n",
     "\n",
     "bootstrap_results = {\n",
-    "    \"nonparametric\": bpars_nonpar,\n",
-    "    \"stratified\": bpars_nonpar_strat,\n",
-    "    \"parametric\": bpars_param,\n",
+    "    \"nonparametric\": btheta_nonpar,\n",
+    "    \"stratified\": btheta_nonpar_strat,\n",
+    "    \"parametric\": btheta_param,\n",
     "}\n",
     "\n",
     "for name, samples in bootstrap_results.items():\n",
@@ -1302,7 +1318,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1339,10 +1355,10 @@
    "source": [
     "def plot_bootstrap_parameters(\n",
     "    samples,\n",
-    "    pars,\n",
-    "    cov_pars,\n",
+    "    theta,\n",
+    "    cov_theta,\n",
     "    title = None,\n",
-    "    pars_labels=(r\"$\\beta_0$\", r\"$\\beta_1$\"),\n",
+    "    theta_labels=(r\"$\\beta_0$\", r\"$\\beta_1$\"),\n",
     "    kind=\"kde\",\n",
     "    fill=True,\n",
     "    fontsize=None,\n",
@@ -1355,17 +1371,17 @@
     "    fontsize = fontsize or {\"ticks\": 16, \"axes_labels\": 16, \"plot_label\": 18, \"plot_title\": 20}\n",
     "\n",
     "    samples = np.asarray(samples, dtype=float)\n",
-    "    pars = np.asarray(pars, dtype=float)\n",
-    "    cov_pars = np.asarray(cov_pars, dtype=float)\n",
+    "    theta = np.asarray(theta, dtype=float)\n",
+    "    cov_theta = np.asarray(cov_theta, dtype=float)\n",
     "\n",
     "    if samples.ndim != 2 or samples.shape[1] != 2:\n",
     "        raise ValueError(\"samples must have shape (n_samples, 2)\")\n",
-    "    if pars.shape != (2,):\n",
-    "        raise ValueError(\"pars must have shape (2,)\")\n",
-    "    if cov_pars.shape != (2, 2):\n",
-    "        raise ValueError(\"cov_pars must have shape (2, 2)\")\n",
-    "    if np.any(np.diag(cov_pars) <= 0):\n",
-    "        raise ValueError(\"cov_pars must have positive diagonal elements\")\n",
+    "    if theta.shape != (2,):\n",
+    "        raise ValueError(\"theta must have shape (2,)\")\n",
+    "    if cov_theta.shape != (2, 2):\n",
+    "        raise ValueError(\"cov_theta must have shape (2, 2)\")\n",
+    "    if np.any(np.diag(cov_theta) <= 0):\n",
+    "        raise ValueError(\"cov_theta must have positive diagonal elements\")\n",
     "    if kind not in {\"kde\", \"scatter\"}:\n",
     "        raise ValueError(\"kind must be 'kde' or 'scatter'\")\n",
     "\n",
@@ -1388,20 +1404,20 @@
     "                      kind=kind, height=figsize[1], **plot_options,)\n",
     "\n",
     "    g.figure.set_size_inches(*figsize)\n",
-    "    g.set_axis_labels(*pars_labels, fontsize=fontsize[\"axes_labels\"],)\n",
+    "    g.set_axis_labels(*theta_labels, fontsize=fontsize[\"axes_labels\"],)\n",
     "    g.ax_joint.tick_params(labelsize=fontsize[\"ticks\"])\n",
     "\n",
     "    # MLE point\n",
-    "    g.ax_joint.plot(*pars, \"ro\", markersize=8, label=\"MLE\", zorder=10,)\n",
+    "    g.ax_joint.plot(*theta, \"ro\", markersize=8, label=\"MLE\", zorder=10,)\n",
     "    g.ax_joint.legend(fontsize=fontsize[\"ticks\"])\n",
     "\n",
     "    # Marginal normal approximations\n",
     "    x = np.linspace(*g.ax_joint.get_xlim(), 300)\n",
     "    y = np.linspace(*g.ax_joint.get_ylim(), 300)\n",
     "\n",
-    "    g.ax_marg_x.plot(x, norm.pdf(x, pars[0], np.sqrt(cov_pars[0, 0])), \n",
+    "    g.ax_marg_x.plot(x, norm.pdf(x, theta[0], np.sqrt(cov_theta[0, 0])), \n",
     "                     color=\"red\", linewidth=2,)\n",
-    "    g.ax_marg_y.plot(norm.pdf(y, pars[1], np.sqrt(cov_pars[1, 1])), y,\n",
+    "    g.ax_marg_y.plot(norm.pdf(y, theta[1], np.sqrt(cov_theta[1, 1])), y,\n",
     "                     color=\"red\", linewidth=2,)\n",
     "\n",
     "    if plot_label is not None:\n",
@@ -1435,8 +1451,8 @@
     "    if len(samples_sel) < 2:\n",
     "        raise ValueError( f\"Too few samples remain inside the selected window for {name!r}.\")\n",
     "\n",
-    "    g = plot_bootstrap_parameters(samples_sel, pars, cov_pars,\n",
-    "                                  pars_labels = [r\"$\\theta_{0}$\", r\"$\\theta_{1}$\"], \n",
+    "    g = plot_bootstrap_parameters(samples_sel, theta, cov_theta,\n",
+    "                                  theta_labels = [r\"$\\theta_{0}$\", r\"$\\theta_{1}$\"], \n",
     "                                  plot_label=plot_labels.get(name))\n",
     "\n",
     "    safe_name = name.lower().replace(\" \", \"_\")\n",
@@ -1459,7 +1475,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1474,22 +1490,22 @@
     }
    ],
    "source": [
-    "def bootstrap_band(xgrid, bootstrap_pars, quantiles):\n",
-    "    values = np.vstack([logit.model(xgrid, sample) for sample in bootstrap_pars])\n",
+    "def bootstrap_band(xgrid, bootstrap_theta, quantiles):\n",
+    "    values = np.vstack([logit.model(xgrid, sample) for sample in bootstrap_theta])\n",
     "    return np.quantile(values, quantiles, axis=0)\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(8, 5))\n",
     "for cls, label in [(0, 'NC'), (1, 'AE')]:\n",
     "    ax.scatter(x[y == cls], y[y == cls], label=label, zorder=4)\n",
     "\n",
-    "ax.plot(xp, logit.model(xp, pars), color='black', lw=2, label='logistic fit', zorder=5)\n",
+    "ax.plot(xp, logit.model(xp, theta), color='black', lw=2, label='logistic fit', zorder=5)\n",
     "\n",
     "bands = [\n",
-    "    ('normal', logit.get_model_quantiles_normal(xp, probs, pars, cov_pars), 'tab:red'),\n",
-    "    ('delta', logit.get_model_quantiles_delta(xp, probs, pars, cov_pars), 'tab:green'),\n",
-    "    ('non-parametric bootstrap', bootstrap_band(xp, bpars_nonpar, probs), 'tab:blue'),\n",
-    "    ('stratified bootstrap', bootstrap_band(xp, bpars_nonpar_strat, probs), 'tab:purple'),\n",
-    "    ('parametric bootstrap', bootstrap_band(xp, bpars_param, probs), 'tab:orange'),\n",
+    "    ('normal', logit.get_model_quantiles_normal(xp, probs, theta, cov_theta), 'tab:red'),\n",
+    "    ('delta', logit.get_model_quantiles_delta(xp, probs, theta, cov_theta), 'tab:green'),\n",
+    "    ('non-parametric bootstrap', bootstrap_band(xp, btheta_nonpar, probs), 'tab:blue'),\n",
+    "    ('stratified bootstrap', bootstrap_band(xp, btheta_nonpar_strat, probs), 'tab:purple'),\n",
+    "    ('parametric bootstrap', bootstrap_band(xp, btheta_param, probs), 'tab:orange'),\n",
     "]\n",
     "\n",
     "for label, (lo, hi), color in bands:\n",