1. #!/usr/bin/env python3
2. """
3. Case 20 arithmetic.
4. ​
5. Source: Pain Resilience Therapy vs. Pain Neuroscience Education, Healthcare
6. 2026, 14(13):1940, doi 10.3390/healthcare14131940, PMID 42450950, PMC13362344.
7. n = 60 completers, 30 per arm. Counts read from their Table 4; RR/ARR/NNT
8. compared against their Table 5.
9. ​
10. Nothing here comes from prose: every count is a cell of Table 4 and every
11. printed value is a cell of Table 5. Run this file to regenerate all numbers.
12. """
13. import numpy as np
14. from scipy.stats import norm, hypergeom
15. ​
16. n1 = n2 = 30 # completers per arm (PNE, PRT); Table 4 rows sum to 30 each
17. ​
18. # name: (PNE yes, PRT yes, printed RR, printed ARR, printed NNT)
19. T = {
20. "STAI": (8, 13, 1.62, 0.17, 6),
21. "FABQ-PA": (8, 13, 1.62, 0.17, 6),
22. "FABQ-W": (4, 9, 2.25, 0.17, 6),
23. "PCS": (6, 17, 2.83, 0.37, 3),
24. "TSK-11": (4, 16, 4.00, 0.40, 3),
25. "PASS-20": (5, 9, 1.80, 0.13, 8),
26. "SER": (15, 24, 1.60, 0.30, 2),
27. "CPAQ": (2, 13, 6.50, 0.37, 3),
28. }
29. ​
30. ​
31. def wilson(x, n, alpha=0.05):
32. z = norm.ppf(1 - alpha / 2)
33. p = x / n
34. d = 1 + z * z / n
35. c = (p + z * z / (2 * n)) / d
36. h = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
37. return p, max(0.0, c - h), min(1.0, c + h)
38. ​
39. ​
40. def newcombe_arr(x1, n1, x2, n2, alpha=0.05):
41. """Newcombe hybrid score interval for p2 - p1 (method 10, 1998)."""
42. p1, l1, u1 = wilson(x1, n1, alpha)
43. p2, l2, u2 = wilson(x2, n2, alpha)
44. d = p2 - p1
45. lo = d - np.sqrt((p2 - l2) ** 2 + (u1 - p1) ** 2)
46. hi = d + np.sqrt((u2 - p2) ** 2 + (p1 - l1) ** 2)
47. return d, lo, hi
48. ​
49. ​
50. def wald_arr(x1, n1, x2, n2, alpha=0.05):
51. p1, p2 = x1 / n1, x2 / n2
52. se = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
53. z = norm.ppf(1 - alpha / 2)
54. return p2 - p1, (p2 - p1) - z * se, (p2 - p1) + z * se
55. ​
56. ​
57. def koopman_rr(x1, n1, x2, n2, alpha=0.05):
58. """Koopman (1984) score interval for RR = p2/p1. Good coverage at small counts."""
59. from scipy.stats import chi2
60. from scipy.optimize import brentq, minimize_scalar
61. crit = chi2.ppf(1 - alpha, 1)
62. ​
63. def T(r):
64. def f(p1):
65. p2 = r * p1
66. if p1 <= 0 or p1 >= 1 or p2 >= 1:
67. return 1e9
68. return ((x1 / n1 - p1) ** 2 / (p1 * (1 - p1) / n1) +
69. (x2 / n2 - p2) ** 2 / (p2 * (1 - p2) / n2))
70. hi = min(1.0, 1 / r - 1e-9)
71. return minimize_scalar(f, bounds=(1e-9, hi), method="bounded").fun
72. ​
73. rhat = (x2 / n2) / (x1 / n1)
74. return (brentq(lambda r: T(r) - crit, 1e-6, rhat),
75. brentq(lambda r: T(r) - crit, rhat, 1e5))
76. ​
77. ​
78. def reciprocal_ci(lo, hi):
79. """NNT interval from the ARR interval; flags sign change (no finite bound)."""
80. if hi <= 0 or lo >= 0:
81. vals = sorted([1 / lo, 1 / hi]) if lo * hi > 0 else (float("-inf"), float("inf"))
82. return vals
83. return (float("-inf"), float("inf"))
84. ​
85. ​
86. print(f"{'outcome':9} {'p1':>6} {'p2':>6} {'ARR':>6} {'ARR 95% CI (Newcombe)':>26} "
87. f"{'half-w':>7} {'ARR/hw':>7} {'NNT':>6} {'NNT CI':>16} {'printed':>14}")
88. rows = []
89. for k, (x1, x2, rr_p, arr_p, nnt_p) in T.items():
90. p1, p2 = x1 / n1, x2 / n2
91. arr, lo, hi = newcombe_arr(x1, n1, x2, n2)
92. hw = (hi - lo) / 2
93. dist = arr / hw if hw > 0 else float("inf")
94. nnt = 1 / arr
95. nlo, nhi = reciprocal_ci(lo, hi)
96. rows.append((k, arr, lo, hi, hw, dist, nnt, nlo, nhi, nnt_p))
97. print(f"{k:9} {p1:6.3f} {p2:6.3f} {arr:6.3f} {'[' + format(lo, '.3f') + ', ' + format(hi, '.3f') + ']':>26} "
98. f"{hw:7.3f} {dist:7.2f} {nnt:6.2f} "
99. f"{('[' + format(nlo, '.2f') + ', ' + format(nhi, '.2f') + ']'):>16} "
100. f"{f'RR {rr_p} ARR {arr_p} NNT {nnt_p}':>14}")
101. ​
102. print("\n-- arithmetic check of the printed table --")
103. for k, (x1, x2, rr_p, arr_p, nnt_p) in T.items():
104. arr = x2 / n2 - x1 / n1
105. rr = (x2 / n2) / (x1 / n1)
106. print(f"{k:9} printed RR {rr_p:5.2f} vs {rr:5.2f} | printed ARR {arr_p:5.2f} vs {arr:5.3f} "
107. f"| 1/ARR = {1 / arr:5.2f} -> printed NNT {nnt_p}")
108. ​
109. print("\n-- the two headline numbers --")
110. for k in ("CPAQ", "SER"):
111. x1, x2 = T[k][0], T[k][1]
112. arr, lo, hi = newcombe_arr(x1, n1, x2, n2)
113. w_lo, w_hi = wald_arr(x1, n1, x2, n2)[1:]
114. rlo, rhi = koopman_rr(x1, n1, x2, n2)
115. print(f"{k}: RR={(x2/n2)/(x1/n1):.2f} Koopman 95% CI [{rlo:.2f}, {rhi:.2f}]")
116. print(f" ARR={arr:.3f} Newcombe [{lo:.3f}, {hi:.3f}] Wald [{w_lo:.3f}, {w_hi:.3f}]")
117. print(f" half-width Newcombe {(hi-lo)/2:.3f} -> ARR/half-width = {arr/((hi-lo)/2):.2f}")
118. nlo, nhi = reciprocal_ci(lo, hi)
119. print(f" NNT = {1/arr:.2f}, 95% CI [{nlo:.2f}, {nhi:.2f}] "
120. f"(a factor of {nhi/nlo:.1f} between the edges)")
121. if w_lo <= 0:
122. print(" under the Wald reading the ARR interval touches 0 -> no finite NNT bound")
123. else:
124. print(f" under the Wald reading NNT CI [{1/w_hi:.2f}, {1/w_lo:.2f}]")
125. ​
126. print("\n-- what the rounding hides (CPAQ) --")
127. x1, x2 = T["CPAQ"][0], T["CPAQ"][1]
128. for d in range(-1, 3):
129. xx1 = x1 + d
130. arr = x2 / n2 - xx1 / n1
131. print(f" if {xx1} of 30 PNE patients had improved instead of {x1}: "
132. f"ARR {arr:.3f}, NNT {1/arr:6.2f}, RR {(x2/n2)/(xx1/n1):6.2f}")