Merge pull request #20 from DeepPSP/fix-make-risk-report

wenh06 · web-flow · commit 9fb2f5cdbe8e · 2026-03-10T00:39:28.000+08:00
Fix errors in make risk report function
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -18,7 +18,7 @@ repos:
     hooks:
       - id: flake8
         additional_dependencies: [pycodestyle>=2.11.0]
-        args: [--max-line-length=128, '--exclude=./.*,build,dist', '--ignore=E501,W503,E231,E203,E251,E202', --count, --statistics, --show-source]
+        args: [--max-line-length=128, '--exclude=./.*,build,dist', '--ignore=E501,W503,E231,E203,E251,E202,E226', --count, --statistics, --show-source]
   - repo: https://github.com/pycqa/isort
     rev: 7.0.0
     hooks:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -26,9 +26,12 @@ Removed
 Fixed
 ~~~~~
 
-- Fix errors in edge cases (mainly when `n_total` equals `ref_total`)
+- (NOT fully resolved yet) Fix errors in edge cases (mainly when `n_total` equals `ref_total`)
   in the computation of difference of two proportions confidence intervals
   using `wang` method.
+- Fix errors in generating LaTeX table for risk report (function ``make_risk_report``):
+  - Escape special characters in LaTeX.
+  - Fix the header for the generated LaTeX table.
 
 Security
 ~~~~~~~~
diff --git a/diff_binom_confint/_applications.py b/diff_binom_confint/_applications.py
@@ -1,3 +1,5 @@
+import re
+import textwrap
 import warnings
 from pathlib import Path
 from typing import Dict, Optional, Sequence, Tuple, Union
@@ -13,6 +15,27 @@
 ]
 
 
+def _latex_escape(line: str) -> str:
+    """Escape LaTeX special characters in a string."""
+    # LaTeX special characters that need escaping:
+    #    \  &  %  $  #  _  {  }  ~  ^
+    # We replace each occurrence with its escaped form.
+    specials = {
+        "\\": r"\textbackslash{}",
+        "&": r"\&",
+        "%": r"\%",
+        "$": r"\$",
+        "#": r"\#",
+        "_": r"\_",
+        "{": r"\{",
+        "}": r"\}",
+        "~": r"\textasciitilde{}",
+        "^": r"\textasciicircum{}",
+    }
+    pattern = re.compile(r"([\\&%$#_{}~^])")
+    return pattern.sub(lambda m: specials[m.group(1)], line)
+
+
 def make_risk_report(
     data_source: Union[pd.DataFrame, Tuple[pd.DataFrame, pd.DataFrame]],
     target: str,
@@ -220,7 +243,7 @@ def make_risk_report(
                     n_positive[ref_item],
                     n_affected[ref_item],
                     conf_level,
-                    method,
+                    diff_method,
                     **kwargs,
                 ).astuple(),
             }
@@ -283,17 +306,37 @@ def make_risk_report(
     if return_type.lower() == "pd":
         return df_risk_table
     elif return_type.lower() == "latex":
-        rows = [line.replace("%", r"\%") for line in df_risk_table.to_latex(header=False, index=False).splitlines()]
-        rows[0] = r"\begin{tabular}{@{\extracolsep{6pt}}lllllll@{}}"
-        rows[2] = (
-            r"\multicolumn{2}{l}{Feature} & \multicolumn{affected_cols}{l}{Affected} & \multicolumn{2}{l}{risk_name Risk ($95\%$ CI)} & risk_name Risk Difference  ($95\%$ CI) \\ \cline{1-2}\cline{3-4}\cline{5-6}\cline{7-7}"
-        )
-        rows[2].replace("risk_name", risk_name).replace("95", str(int(conf_level * 100)))
+        latex_body = df_risk_table.to_latex(header=False, index=False)
+        body_lines = latex_body.splitlines()
+
         if is_split:
-            rows[2].replace("affected_cols", "3")
+            header = rf"""
+                \begin{{tabular}}{{@{{\extracolsep{{6pt}}}}llllllll@{{}}}}
+                \toprule
+                \multicolumn{{2}}{{l}}{{Feature}} &
+                \multicolumn{{3}}{{l}}{{Affected}} &
+                \multicolumn{{2}}{{l}}{{{risk_name} Risk (${int(conf_level*100)}\%$ CI)}} &
+                {risk_name} Risk Difference (${int(conf_level*100)}\%$ CI) \\
+                \cmidrule(lr){{1-2}}\cmidrule(lr){{3-5}}\cmidrule(lr){{6-7}}\cmidrule(lr){{8-8}}
+                & & n & \% & t/v & n & \% & \\ \midrule
+                """
         else:
-            rows[2].replace("affected_cols", "2")
-        ret_lines = "\n".join(rows)
+            header = rf"""
+                \begin{{tabular}}{{@{{\extracolsep{{6pt}}}}lllllll@{{}}}}
+                \toprule
+                \multicolumn{{2}}{{l}}{{Feature}} &
+                \multicolumn{{2}}{{l}}{{Affected}} &
+                \multicolumn{{2}}{{l}}{{{risk_name} Risk (${int(conf_level*100)}\%$ CI)}} &
+                {risk_name} Risk Difference (${int(conf_level*100)}\%$ CI) \\
+                \cmidrule(lr){{1-2}}\cmidrule(lr){{3-4}}\cmidrule(lr){{5-6}}\cmidrule(lr){{7-7}}
+                & & n & \% & n & \% & \\ \midrule
+                """
+        # remove extra leading spaces
+        header = textwrap.dedent(header).strip()
+
+        body = "\n".join(_latex_escape(line) for line in body_lines[5:-1])
+
+        ret_lines = header + "\n" + body + "\n\\end{tabular}"
         if save_path is not None:
             save_path.with_suffix(".tex").write_text(ret_lines)
         return ret_lines
@@ -303,3 +346,5 @@ def make_risk_report(
         return df_risk_table.to_html(index=False)
     elif return_type.lower() == "dict":
         return ret_dict
+    else:
+        raise ValueError(f"Unsupported return_type {repr(return_type)}")
diff --git a/diff_binom_confint/_specials/_wang.py b/diff_binom_confint/_specials/_wang.py
@@ -43,16 +43,17 @@ def wang_binomial_ci(
     sides : Union[str, int], optional
         sides: str or int, default "two-sided",
         the sides of the confidence interval, should be one of
-        "two-sided" (aliases "2-sided", "two_sided", "2_sided", "2-sides", "two_sides", "two-sides", "2_sides", "ts", "t", "two", "2", 2),
+        "two-sided" (aliases "2-sided", "two_sided", "2_sided", "2-sides",
+        "two_sides", "two-sides", "2_sides", "ts", "t", "two", "2", 2),
         "left-sided" (aliases "left_sided", "left", "ls", "l"),
         "right-sided" (aliases "right_sided", "right", "rs", "r"),
         case insensitive.
     precision : float, optional
-        Precision for the search algorithm, by default 1e-5
+        Precision for the search algorithm, by default 1e-5.
     grid_one : int, optional
-        Number of grid points in first step, by default 30
+        Number of grid points in first step, by default 30.
     grid_two : int, optional
-        Number of grid points in second step, by default 20
+        Number of grid points in second step, by default 20.
     verbose : bool, optional
         Verbosity for debug message.
 
@@ -101,7 +102,6 @@ def wang_binomial_ci(
             precision,
             grid_one,
             grid_two,
-            verbose,
         )
         if verbose:
             print(f"Left CI: {ci_l}")
@@ -115,7 +115,6 @@ def wang_binomial_ci(
             precision,
             grid_one,
             grid_two,
-            verbose,
         )
         if verbose:
             print(f"Right CI: {ci_u}")
@@ -124,7 +123,7 @@ def wang_binomial_ci(
         return ConfidenceInterval(lower, upper, estimate, conf_level, "wang", sides_val)
     else:
         ci = binomial_ci_one_sided(
-            n_positive, n_total, ref_positive, ref_total, conf_level, sides_val, precision, grid_one, grid_two, verbose
+            n_positive, n_total, ref_positive, ref_total, conf_level, sides_val, precision, grid_one, grid_two
         )
         return ConfidenceInterval(ci[1], ci[2], ci[0], conf_level, "wang", sides_val)
 
@@ -139,7 +138,6 @@ def binomial_ci_one_sided(
     precision: float,
     grid_one: int,
     grid_two: int,
-    verbose: bool = False,
 ) -> List[float]:
     """Helper function that calculates one-sided confidence interval.
 
@@ -154,11 +152,12 @@ def binomial_ci_one_sided(
     ref_total : int
         total number of samples of the reference.
     conf_level : float, optional
-        Confidence level, by default 0.95
+        Confidence level, by default 0.95.
     sides : Union[str, int], optional
         sides: str or int, default "two-sided",
         the sides of the confidence interval, should be one of
-        "two-sided" (aliases "2-sided", "two_sided", "2_sided", "2-sides", "two_sides", "two-sides", "2_sides", "ts", "t", "two", "2", 2),
+        "two-sided" (aliases "2-sided", "two_sided", "2_sided", "2-sides",
+        "two_sides", "two-sides", "2_sides", "ts", "t", "two", "2", 2),
         "left-sided" (aliases "left_sided", "left", "ls", "l"),
         "right-sided" (aliases "right_sided", "right", "rs", "r"),
         case insensitive.
@@ -168,8 +167,6 @@ def binomial_ci_one_sided(
         Number of grid points in first step, by default 30.
     grid_two : int, optional
         Number of grid points in second step, by default 20.
-    verbose : bool, optional
-        Verbosity for debug message, by default False.
 
     Returns
     -------
@@ -209,7 +206,7 @@ def binomial_ci_one_sided(
     f[:, 2] = (p1hat - p0hat) / np.sqrt(denom)
 
     # Sort f by the third column in descending order
-    f = f[(-f[:, 2]).argsort(kind="stable"), :]
+    f = f[(-f[:, 2]).argsort(), :]
 
     allvector = np.round(f[:, 0] * (m + 2) + f[:, 1]).astype(int)
     allvectormove = np.round((f[:, 0] + 1) * (m + 3) + (f[:, 1] + 1)).astype(int)
@@ -273,7 +270,7 @@ def binomial_ci_one_sided(
 
         # Generate N
         n_arr = np.unique(np.vstack((a, b)), axis=0)
-        nvector = ((n_arr[:, 0] + 1) * (m + 3) + n_arr[:, 1] + 1).astype(int)  # type: ignore
+        nvector = ((n_arr[:, 0] + 1) * (m + 3) + n_arr[:, 1] + 1).astype(int)
         nvector = nvector[np.isin(nvector, allvectormove)]
 
         skvector = ((s[:kk, 0] + 1) * (m + 3) + s[:kk, 1] + 1).astype(int)
@@ -308,7 +305,6 @@ def binomial_ci_one_sided(
         else:
             length_nc = nc_arr.shape[0]
 
-        ncmax = 0  # avoid pylance warning
         for ci in range(length_nc):
             ls_arr[kk, 0:2] = nc_arr[ci, 0:2]
             i1_vec = ls_arr[: (kk + 1), 0]
@@ -369,7 +365,7 @@ def binomial_ci_one_sided(
         if length_nc >= 2:
             valid = ~np.isnan(nc_arr[:, 0])
             ncnomiss = nc_arr[valid]
-            ncnomiss = ncnomiss[(-ncnomiss[:, 2]).argsort(kind="stable"), :]
+            ncnomiss = ncnomiss[(-ncnomiss[:, 2]).argsort(), :]
             morepoint = np.sum(ncnomiss[:, 2] >= ncnomiss[0, 2] - delta)
             if morepoint >= 2:
                 ls_arr[kk : kk + morepoint, 0:2] = ncnomiss[:morepoint, 0:2]
@@ -435,8 +431,7 @@ def binomial_ci_one_sided(
 
         kk1 = kk
 
-    # output = [val.item() if isinstance(val, np.generic) else val for val in output]
-    output = np.array(output).tolist()
+    output = [val.item() if isinstance(val, np.generic) else val for val in output]
 
     return output
 
@@ -463,32 +458,27 @@ def _prob2step(delv, delta, n, m, i1, i2, grid_one, grid_two):
         p0 = np.linspace(-delv + delta, 1 - delta, grid_one)
     else:
         p0 = np.linspace(delta, 1 - delv - delta, grid_one)
-
     i1 = np.atleast_1d(i1)
     i2 = np.atleast_1d(i2)
     part1 = np.log(comb(n, i1))[:, None] + np.outer(i1, np.log(p0 + delv)) + np.outer(n - i1, np.log(1 - p0 - delv))
     part2 = np.log(comb(m, i2))[:, None] + np.outer(i2, np.log(p0)) + np.outer(m - i2, np.log(1 - p0))
     sumofprob = np.exp(part1 + part2).sum(axis=0)
 
-    # mansum = sumofprob.max()
-    # atol = 1e-14 * (mansum if mansum > 0 else 1.0)
-    # plateau_idx = np.where(np.isclose(sumofprob, mansum, rtol=0.0, atol=atol))[0]
-    plateau_idx = np.where(sumofprob == sumofprob.max())[0]
+    # plateau-aware refinement (R: which(sumofprob == max(sumofprob)))
+    mansum = sumofprob.max()
+    atol = 1e-14 * (mansum if mansum > 0 else 1.0)
+    plateau_idx = np.where(np.isclose(sumofprob, mansum, rtol=0.0, atol=atol))[0]
     leftmost = plateau_idx.min()
     rightmost = plateau_idx.max()
 
-    denom = (grid_one - 1) if grid_one > 1 else 1
-    stepv = (p0[-1] - p0[0]) / denom
-    # lowerb = max(p0[0], p0[rightmost] - stepv) + delta
-    # upperb = min(p0[-1], p0[leftmost] + stepv) - delta
+    stepv = (p0[-1] - p0[0]) / grid_one
+    lowerb = max(p0[0], p0[rightmost] - stepv) + delta
+    upperb = min(p0[-1], p0[leftmost] + stepv) - delta
 
-    raw_lowerb = max(p0[0], p0[rightmost] - stepv) + delta
-    raw_upperb = min(p0[-1], p0[leftmost] + stepv) - delta
-    if raw_lowerb <= raw_upperb:
-        lowerb, upperb = raw_lowerb, raw_upperb
-    else:
-        # Ensure bounds are ordered for linspace; swap if necessary
-        lowerb, upperb = raw_upperb, raw_lowerb
+    # stepv = (p0[-1] - p0[0]) / grid_one
+    # maxloc = np.argmax(sumofprob)
+    # lowerb = max(p0[0], p0[maxloc] - stepv) + delta
+    # upperb = min(p0[-1], p0[maxloc] + stepv) - delta
 
     p0 = np.linspace(lowerb, upperb, grid_two)
     part1 = np.log(comb(n, i1))[:, None] + np.outer(i1, np.log(p0 + delv)) + np.outer(n - i1, np.log(1 - p0 - delv))
@@ -502,25 +492,28 @@ def _prob2steplmin(delv, delta, n, m, i1, i2, grid_one, grid_two):
         p0 = np.linspace(-delv + delta, 1 - delta, grid_one)
     else:
         p0 = np.linspace(delta, 1 - delv - delta, grid_one)
-
     i1 = np.atleast_1d(i1)
     i2 = np.atleast_1d(i2)
     part1 = np.log(comb(n, i1))[:, None] + np.outer(i1, np.log(p0 + delv)) + np.outer(n - i1, np.log(1 - p0 - delv))
     part2 = np.log(comb(m, i2))[:, None] + np.outer(i2, np.log(p0)) + np.outer(m - i2, np.log(1 - p0))
     sumofprob = np.exp(part1 + part2).sum(axis=0)
 
-    # mansum = sumofprob.min()
-    # atol = 1e-14 * (abs(mansum) if mansum != 0 else 1.0)
-    # plateau_idx = np.where(np.isclose(sumofprob, mansum, rtol=0.0, atol=atol))[0]
-    plateau_idx = np.where(sumofprob == sumofprob.min())[0]
+    # plateau-aware refinement for minima (R: which(sumofprob == min(sumofprob)))
+    mansum = sumofprob.min()
+    atol = 1e-14 * (abs(mansum) if mansum != 0 else 1.0)
+    plateau_idx = np.where(np.isclose(sumofprob, mansum, rtol=0.0, atol=atol))[0]
     leftmost = plateau_idx.min()
     rightmost = plateau_idx.max()
 
-    denom = (grid_one - 1) if grid_one > 1 else 1
-    stepv = (p0[-1] - p0[0]) / denom
+    stepv = (p0[-1] - p0[0]) / grid_one
     lowerb = max(p0[0], p0[rightmost] - stepv) + delta
     upperb = min(p0[-1], p0[leftmost] + stepv) - delta
 
+    # stepv = (p0[-1] - p0[0]) / grid_one
+    # minloc = np.argmin(sumofprob)
+    # lowerb = max(p0[0], p0[minloc] - stepv) + delta
+    # upperb = min(p0[-1], p0[minloc] + stepv) - delta
+
     p0 = np.linspace(lowerb, upperb, grid_two)
     part1 = np.log(comb(n, i1))[:, None] + np.outer(i1, np.log(p0 + delv)) + np.outer(n - i1, np.log(1 - p0 - delv))
     part2 = np.log(comb(m, i2))[:, None] + np.outer(i2, np.log(p0)) + np.outer(m - i2, np.log(1 - p0))
diff --git a/test/test_applications.py b/test/test_applications.py
@@ -31,10 +31,10 @@ def test_make_risk_report():
         [str(_TMP_DIR / "risk-report"), None],  # save_path
     )
 
-    for data_source, ref_classes, risk_name, return_type, save_path in grid:
+    for data_source, rc, risk_name, return_type, save_path in grid:
         report = make_risk_report(
             data_source=data_source,
-            ref_classes=ref_classes,
+            ref_classes=rc,
             risk_name=risk_name,
             return_type=return_type,
             save_path=save_path,
@@ -49,6 +49,16 @@ def test_make_risk_report():
         elif return_type in ("latex", "md", "markdown", "html"):
             assert isinstance(report, str)
 
+    with pytest.raises(ValueError, match="Unsupported return_type"):
+        make_risk_report(
+            data_source=df_test,
+            ref_classes=ref_classes,
+            risk_name="Seizure",
+            target="HasSeizure",
+            positive_class="Yes",
+            return_type="xxx",
+        )
+
     with pytest.raises(ValueError, match=f"target {repr('xxx')} not in the columns"):
         make_risk_report(
             data_source=df_test,
