Updating logistic reg boostrapping plots

python/logistic/src/logit.py

@@ -49,7 +49,7 @@ def logit_poly_model(x, b):
     X = np.column_stack([x**i for i in range(len(b))])
     F = X @ b
 
-    return 1/(1 + np.exp(-F))
+    return scipy.special.expit(F)
 
 
 """
@@ -230,7 +230,7 @@ def logit_poly_model_quantiles_normal(x, probs, mean_b, cov):
     Q = locs + np.outer(scipy.stats.norm.ppf(probs),scales)
     
     # convert logit to expit
-    return 1/(1 + np.exp(-Q))
+    return scipy.special.expit(Q)
 
 """
     Calculating quantiles using delta method of the model values 
@@ -267,8 +267,9 @@ def logit_poly_model_quantiles_delta(x, probs, mean_b, cov):
     X = np.column_stack([x**i for i in range(len(mean_b))])
     F = X@mean_b
 
-    locs = 1/(1 + np.exp(-F))
-    scales = np.sqrt(np.diag(X@cov@X.T))/(4*np.cosh(F/2)**2)
+    locs = scipy.special.expit(F)
+    deriv = locs * (1.0 - locs)
+    scales = np.sqrt(np.diag(X@cov@X.T))* deriv
 
     # computing quantiles of logit
     Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)

python/logistic/src/report_utils.py

@@ -3,6 +3,286 @@ import pickle
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
+import seaborn as sns
+import scipy
+
+def plot_pdf(
+    data1,
+    data2,
+    results_path,
+    file1,
+    file2,
+    pars,
+    cov_pars,
+    fontsize = {"ticks":16, "axes_labels":16, "plot_labels": 18},
+    width = 6,
+    height = 4
+):
+    """
+    Plot joint KDEs for two 2D datasets.
+
+    Parameters
+    ----------
+    data1 : ndarray of shape (n, 2)
+        First dataset.
+    data2 : ndarray of shape (m, 2)
+        Second dataset, typically a zoomed/subselected version of data1.
+    file1 : str
+        Output filename for the first figure.
+    file2 : str
+        Output filename for the second figure.
+    pars : array-like of length 2
+        Mean/MLE parameter estimates [beta0, beta1].
+    cov_pars : ndarray of shape (2, 2)
+        Covariance matrix of the parameter estimates.
+    fontsize : int, optional
+        Base font size.
+    fontsize_label : int, optional
+        Font size for subplot labels.
+    """
+    data1 = np.asarray(data1)
+    data2 = np.asarray(data2)
+    pars = np.asarray(pars)
+    cov_pars = np.asarray(cov_pars)
+
+    if data1.ndim != 2 or data1.shape[1] != 2:
+        raise ValueError("data1 must have shape (n, 2)")
+    if data2.ndim != 2 or data2.shape[1] != 2:
+        raise ValueError("data2 must have shape (m, 2)")
+    if pars.shape != (2,):
+        raise ValueError("pars must have shape (2,)")
+    if cov_pars.shape != (2, 2):
+        raise ValueError("cov_pars must have shape (2, 2)")
+
+    sns.set_theme(rc={"figure.figsize": (width, height)})
+
+    # FIRST
+    g1 = sns.jointplot(x=data1[:, 0], y=data1[:, 1], kind="kde", fill=True)
+    g1.ax_joint.tick_params(labelsize=fontsize["ticks"])
+    g1.set_axis_labels(xlabel=r"$\beta_0$", ylabel=r"$\beta_1$", fontsize=fontsize["axes_labels"])
+    g1.figure.text(0.04, 0.9, "(a)", fontsize=fontsize["plot_labels"], ha="center", va="center")
+    
+    g1.figure.savefig(os.path.join(results_path,  file1), bbox_inches="tight")
+    plt.show()
+
+    # SECOND
+    g2 = sns.jointplot(x=data2[:, 0], y=data2[:, 1], kind="kde", fill=True)
+    g2.ax_joint.tick_params(labelsize=fontsize["ticks"])
+    g2.set_axis_labels(xlabel=r"$\beta_0$", ylabel=r"$\beta_1$", fontsize=fontsize["axes_labels"])
+    g2.figure.text(0.04, 0.9, "(b)", fontsize=fontsize["plot_labels"], ha="center", va="center")
+
+    # MLE point on the joint axes
+    g2.ax_joint.plot(
+        pars[0], pars[1],
+        marker="o",
+        markersize=10,
+        markeredgecolor="red",
+        markerfacecolor="red",
+        linestyle="None",
+    )
+
+    # Normal curves on the margins
+    x_vals = np.linspace(np.min(data2[:, 0]), np.max(data2[:, 0]), 100)
+    g2.ax_marg_x.plot(
+        x_vals,
+        scipy.stats.norm.pdf(x_vals, loc=pars[0], scale=np.sqrt(cov_pars[0, 0])),
+        "r",
+    )
+
+    y_vals = np.linspace(np.min(data2[:, 1]), np.max(data2[:, 1]), 100)
+    g2.ax_marg_y.plot(
+        scipy.stats.norm.pdf(y_vals, loc=pars[1], scale=np.sqrt(cov_pars[1, 1])),
+        y_vals,
+        "r",
+    )
+    
+    g2.figure.savefig(os.path.join(results_path, file2), bbox_inches="tight")
+    plt.show()
+    
+
+def plot_logit_ci_comparison(
+    x,
+    y,
+    pars,
+    cov_pars,
+    probs,
+    results_path,
+    logit,
+    bpars_nonpar=None,
+    bpars_nonpar_strat=None,
+    bpars_param=None,
+    filename="logit_cmp_CI.pdf",
+    title="cond. prob. for AE",
+    figsize=(8, 5),
+    alpha_fill=0.12,
+    boundary_linewidth=1.4,
+):
+    """
+    Plot observed binary data, fitted logistic model, and confidence bands
+    with styled boundaries for better visibility in overlapping regions.
+
+    Parameters
+    ----------
+    x : array-like
+        Predictor values.
+    y : array-like
+        Binary response values, expected to contain 0/1.
+    pars : array-like
+        Fitted model parameters.
+    cov_pars : array-like
+        Covariance matrix of fitted parameters.
+    probs : sequence of float
+        Quantiles for CI bands, e.g. [0.025, 0.975].
+    results_path : str
+        Directory where the figure will be saved.
+    logit : module or object
+        Object providing:
+            - logit_poly_model(x, pars)
+            - logit_poly_model_quantiles_normal(x, probs, pars, cov_pars)
+            - logit_poly_model_quantiles_delta(x, probs, pars, cov_pars)
+    bpars_nonpar : sequence, optional
+        Bootstrapped parameter samples for nonparametric bootstrap.
+    bpars_nonpar_strat : sequence, optional
+        Bootstrapped parameter samples for stratified nonparametric bootstrap.
+    bpars_param : sequence, optional
+        Bootstrapped parameter samples for parametric bootstrap.
+    filename : str, optional
+        Output figure filename.
+    title : str, optional
+        Plot title.
+    figsize : tuple, optional
+        Figure size.
+    alpha_fill : float, optional
+        Transparency of CI fills.
+    boundary_linewidth : float, optional
+        Line width of CI boundaries.
+
+    Returns
+    -------
+    fig, ax
+        Matplotlib figure and axes.
+    """
+    x = np.asarray(x)
+    y = np.asarray(y)
+    pars = np.asarray(pars)
+    cov_pars = np.asarray(cov_pars)
+
+    if x.ndim != 1:
+        raise ValueError("x must be a 1D array")
+    if y.ndim != 1:
+        raise ValueError("y must be a 1D array")
+    if len(x) != len(y):
+        raise ValueError("x and y must have the same length")
+    if pars.ndim != 1:
+        raise ValueError("pars must be a 1D array")
+    if cov_pars.ndim != 2:
+        raise ValueError("cov_pars must be a 2D array")
+
+    fig, ax = plt.subplots(figsize=figsize)
+
+    ax.set_title(title)
+    ax.set_xlabel(r"$x = \max_{visit}\ \mathrm{SUV}(visit, p)$")
+    ax.set_ylabel(r"$P(\mathrm{AE}\mid X=x)$")
+
+    # observed data
+    class_labels = {0: "NC", 1: "AE"}
+    for cls, lab in class_labels.items():
+        mask = (y == cls)
+        ax.scatter(x[mask], y[mask], label=lab, zorder=3)
+
+    # fitted model
+    xp = np.linspace(np.min(x), np.max(x), 200)
+    yp = logit.logit_poly_model(xp, pars)
+    ax.plot(xp, yp, label="logistic reg", linewidth=2, zorder=4)
+
+    # style map for all CI types
+    ci_styles = {
+        "normal": {
+            "color": "red",
+            "ls_low": "--",
+            "ls_up": "--",
+        },
+        "delta": {
+            "color": "green",
+            "ls_low": "-.",
+            "ls_up": "-.",
+        },
+        "nonpar": {
+            "color": "blue",
+            "ls_low": ":",
+            "ls_up": ":",
+        },
+        "nonpar-strat": {
+            "color": "orange",
+            "ls_low": (0, (5, 2)),
+            "ls_up": (0, (5, 2)),
+        },
+        "param": {
+            "color": "purple",
+            "ls_low": "-",
+            "ls_up": "-",
+        },
+    }
+
+    def add_ci_band(xgrid, bounds, label, style):
+        """Add a filled confidence band and styled boundary lines."""
+        lower, upper = bounds
+
+        ax.fill_between(
+            xgrid,
+            lower,
+            upper,
+            color=style["color"],
+            alpha=alpha_fill,
+            label=f"CI:{label}",
+            zorder=1,
+        )
+
+        ax.plot(
+            xgrid,
+            lower,
+            color=style["color"],
+            linestyle=style["ls_low"],
+            linewidth=boundary_linewidth,
+            zorder=2,
+        )
+        ax.plot(
+            xgrid,
+            upper,
+            color=style["color"],
+            linestyle=style["ls_up"],
+            linewidth=boundary_linewidth,
+            zorder=2,
+        )
+
+    # asymptotic CIs
+    asymptotic_methods = ["normal", "delta"]
+    for lab in asymptotic_methods:
+        qm = getattr(logit, f"logit_poly_model_quantiles_{lab}")
+        bounds = qm(xp, probs, pars, cov_pars)
+        add_ci_band(xp, bounds, lab, ci_styles[lab])
+
+    # bootstrap CIs
+    bootstrap_sets = [
+        ("nonpar", bpars_nonpar),
+        ("nonpar-strat", bpars_nonpar_strat),
+        ("param", bpars_param),
+    ]
+
+    for lab, bpars in bootstrap_sets:
+        if bpars is None:
+            continue
+
+        model_values = np.array([logit.logit_poly_model(xp, p) for p in bpars])
+        bounds = np.quantile(model_values, probs, axis=0)
+        add_ci_band(xp, bounds, lab, ci_styles[lab])
+
+    ax.legend(loc="lower right")
+
+    os.makedirs(results_path, exist_ok=True)
+    outpath = os.path.join(results_path, filename)
+    fig.savefig(outpath, bbox_inches="tight")
+    plt.show()
 
 
 # =========================