#!/usr/bin/env python3
"""
Case 20 arithmetic.

Source: Pain Resilience Therapy vs. Pain Neuroscience Education, Healthcare
2026, 14(13):1940, doi 10.3390/healthcare14131940, PMID 42450950, PMC13362344.
n = 60 completers, 30 per arm. Counts read from their Table 4; RR/ARR/NNT
compared against their Table 5.

Nothing here comes from prose: every count is a cell of Table 4 and every
printed value is a cell of Table 5. Run this file to regenerate all numbers.
"""
import numpy as np
from scipy.stats import norm, hypergeom

n1 = n2 = 30  # completers per arm (PNE, PRT); Table 4 rows sum to 30 each

# name: (PNE yes, PRT yes, printed RR, printed ARR, printed NNT)
T = {
    "STAI":     (8, 13, 1.62, 0.17, 6),
    "FABQ-PA":  (8, 13, 1.62, 0.17, 6),
    "FABQ-W":   (4,  9, 2.25, 0.17, 6),
    "PCS":      (6, 17, 2.83, 0.37, 3),
    "TSK-11":   (4, 16, 4.00, 0.40, 3),
    "PASS-20":  (5,  9, 1.80, 0.13, 8),
    "SER":      (15, 24, 1.60, 0.30, 2),
    "CPAQ":     (2, 13, 6.50, 0.37, 3),
}


def wilson(x, n, alpha=0.05):
    z = norm.ppf(1 - alpha / 2)
    p = x / n
    d = 1 + z * z / n
    c = (p + z * z / (2 * n)) / d
    h = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
    return p, max(0.0, c - h), min(1.0, c + h)


def newcombe_arr(x1, n1, x2, n2, alpha=0.05):
    """Newcombe hybrid score interval for p2 - p1 (method 10, 1998)."""
    p1, l1, u1 = wilson(x1, n1, alpha)
    p2, l2, u2 = wilson(x2, n2, alpha)
    d = p2 - p1
    lo = d - np.sqrt((p2 - l2) ** 2 + (u1 - p1) ** 2)
    hi = d + np.sqrt((u2 - p2) ** 2 + (p1 - l1) ** 2)
    return d, lo, hi


def wald_arr(x1, n1, x2, n2, alpha=0.05):
    p1, p2 = x1 / n1, x2 / n2
    se = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    z = norm.ppf(1 - alpha / 2)
    return p2 - p1, (p2 - p1) - z * se, (p2 - p1) + z * se


def koopman_rr(x1, n1, x2, n2, alpha=0.05):
    """Koopman (1984) score interval for RR = p2/p1. Good coverage at small counts."""
    from scipy.stats import chi2
    from scipy.optimize import brentq, minimize_scalar
    crit = chi2.ppf(1 - alpha, 1)

    def T(r):
        def f(p1):
            p2 = r * p1
            if p1 <= 0 or p1 >= 1 or p2 >= 1:
                return 1e9
            return ((x1 / n1 - p1) ** 2 / (p1 * (1 - p1) / n1) +
                    (x2 / n2 - p2) ** 2 / (p2 * (1 - p2) / n2))
        hi = min(1.0, 1 / r - 1e-9)
        return minimize_scalar(f, bounds=(1e-9, hi), method="bounded").fun

    rhat = (x2 / n2) / (x1 / n1)
    return (brentq(lambda r: T(r) - crit, 1e-6, rhat),
            brentq(lambda r: T(r) - crit, rhat, 1e5))


def reciprocal_ci(lo, hi):
    """NNT interval from the ARR interval; flags sign change (no finite bound)."""
    if hi <= 0 or lo >= 0:
        vals = sorted([1 / lo, 1 / hi]) if lo * hi > 0 else (float("-inf"), float("inf"))
        return vals
    return (float("-inf"), float("inf"))


print(f"{'outcome':9} {'p1':>6} {'p2':>6} {'ARR':>6} {'ARR 95% CI (Newcombe)':>26} "
      f"{'half-w':>7} {'ARR/hw':>7} {'NNT':>6} {'NNT CI':>16} {'printed':>14}")
rows = []
for k, (x1, x2, rr_p, arr_p, nnt_p) in T.items():
    p1, p2 = x1 / n1, x2 / n2
    arr, lo, hi = newcombe_arr(x1, n1, x2, n2)
    hw = (hi - lo) / 2
    dist = arr / hw if hw > 0 else float("inf")
    nnt = 1 / arr
    nlo, nhi = reciprocal_ci(lo, hi)
    rows.append((k, arr, lo, hi, hw, dist, nnt, nlo, nhi, nnt_p))
    print(f"{k:9} {p1:6.3f} {p2:6.3f} {arr:6.3f} {'[' + format(lo, '.3f') + ', ' + format(hi, '.3f') + ']':>26} "
          f"{hw:7.3f} {dist:7.2f} {nnt:6.2f} "
          f"{('[' + format(nlo, '.2f') + ', ' + format(nhi, '.2f') + ']'):>16} "
          f"{f'RR {rr_p} ARR {arr_p} NNT {nnt_p}':>14}")

print("\n-- arithmetic check of the printed table --")
for k, (x1, x2, rr_p, arr_p, nnt_p) in T.items():
    arr = x2 / n2 - x1 / n1
    rr = (x2 / n2) / (x1 / n1)
    print(f"{k:9} printed RR {rr_p:5.2f} vs {rr:5.2f} | printed ARR {arr_p:5.2f} vs {arr:5.3f} "
          f"| 1/ARR = {1 / arr:5.2f} -> printed NNT {nnt_p}")

print("\n-- the two headline numbers --")
for k in ("CPAQ", "SER"):
    x1, x2 = T[k][0], T[k][1]
    arr, lo, hi = newcombe_arr(x1, n1, x2, n2)
    w_lo, w_hi = wald_arr(x1, n1, x2, n2)[1:]
    rlo, rhi = koopman_rr(x1, n1, x2, n2)
    print(f"{k}: RR={(x2/n2)/(x1/n1):.2f} Koopman 95% CI [{rlo:.2f}, {rhi:.2f}]")
    print(f"      ARR={arr:.3f} Newcombe [{lo:.3f}, {hi:.3f}] Wald [{w_lo:.3f}, {w_hi:.3f}]")
    print(f"      half-width Newcombe {(hi-lo)/2:.3f} -> ARR/half-width = {arr/((hi-lo)/2):.2f}")
    nlo, nhi = reciprocal_ci(lo, hi)
    print(f"      NNT = {1/arr:.2f}, 95% CI [{nlo:.2f}, {nhi:.2f}] "
          f"(a factor of {nhi/nlo:.1f} between the edges)")
    if w_lo <= 0:
        print("      under the Wald reading the ARR interval touches 0 -> no finite NNT bound")
    else:
        print(f"      under the Wald reading NNT CI [{1/w_hi:.2f}, {1/w_lo:.2f}]")

print("\n-- what the rounding hides (CPAQ) --")
x1, x2 = T["CPAQ"][0], T["CPAQ"][1]
for d in range(-1, 3):
    xx1 = x1 + d
    arr = x2 / n2 - xx1 / n1
    print(f"  if {xx1} of 30 PNE patients had improved instead of {x1}: "
          f"ARR {arr:.3f}, NNT {1/arr:6.2f}, RR {(x2/n2)/(xx1/n1):6.2f}")