diff --git "a/src/streamlit_app.py" "b/src/streamlit_app.py"
--- "a/src/streamlit_app.py"
+++ "b/src/streamlit_app.py"
@@ -1,40 +1,4274 @@
-import altair as alt
+import streamlit as st
import numpy as np
import pandas as pd
-import streamlit as st
+import plotly.graph_objects as go
+from scipy.stats import norm
+
+
+def _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=None,
+ nobs=None,
+ alpha=None,
+ power=None,
+ nobs_name="nobs",
+ **extra_kw,
+):
+ """Call solver.solve_power with the parameter to solve for set to None."""
+ effective = {}
+ for k, v in [("effect_size", effect_size), ("alpha", alpha), ("power", power)]:
+ if v is not None:
+ effective[k] = v
+ effective[nobs_name] = nobs
+ effective.update(extra_kw)
+
+ if analysis_mode == "Post Hoc":
+ effective["power"] = None
+ elif analysis_mode == "Sensitivity":
+ effective["effect_size"] = None
+ elif analysis_mode == "Criterion":
+ effective["alpha"] = None
+ else: # A Priori (or Compromise, handled separately)
+ effective[nobs_name] = None
+
+ return solver.solve_power(**effective)
+
+
+def _solve_n_binary(func_n, target_n, lo=2, hi=1000000):
+ """Binary search an integer N to make func_n(N) >= target_n (max 100 iterations)."""
+ for _ in range(100):
+ mid = (lo + hi) // 2
+ if func_n(mid):
+ hi = mid
+ else:
+ lo = mid + 1
+ if lo >= hi:
+ break
+ return hi
+
+
+def _solve_compromise(
+ solver,
+ effect_size,
+ nobs,
+ cost_ratio,
+ alternative="two-sided",
+ nobs_name="nobs",
+ **extra_kw,
+):
+ """Brent's method to find α such that β/α = cost_ratio for given N and effect size.
+ Returns dict {'alpha': adjusted_alpha, 'power': achieved_power}.
+ """
+ from scipy.optimize import brentq
+
+ def f(alpha_candidate):
+ kw = {
+ nobs_name: nobs,
+ "effect_size": effect_size,
+ "alpha": alpha_candidate,
+ "power": None,
+ }
+ kw.update(extra_kw)
+ try:
+ achieved = solver.solve_power(**kw)
+ except Exception:
+ return 1e6
+ beta = max(1e-10, 1 - achieved)
+ return beta / alpha_candidate - cost_ratio
+
+ lo, hi = 1e-8, 0.5
+ f_lo, f_hi = f(lo), f(hi)
+ if f_lo * f_hi > 0:
+ return None
+ alpha_solved = brentq(f, lo, hi)
+ kw = {
+ nobs_name: nobs,
+ "effect_size": effect_size,
+ "alpha": alpha_solved,
+ "power": None,
+ }
+ kw.update(extra_kw)
+ power_solved = solver.solve_power(**kw)
+ return {"alpha": alpha_solved, "power": power_solved}
+
+
+def render_power_calculator(params, analysis_mode="A Priori"):
+ """Render power analysis results (a_priori / post_hoc / sensitivity / criterion / compromise)."""
+
+ atype = params["type"]
+ alpha = params["alpha"]
+ power = params["power"]
+ tails = params["tails"]
+ alternative = "two-sided" if tails == "Two-tailed" else "larger"
+
+ # Determine z critical values
+ if alternative == "two-sided":
+ z_alpha = norm.ppf(1 - alpha / 2)
+ else:
+ z_alpha = norm.ppf(1 - alpha)
+ z_beta = norm.ppf(power)
+ _z_sub = r"z_{1-\alpha/2}" if alternative == "two-sided" else r"z_{1-\alpha}"
+
+ # Extract global adjustment parameters (with defaults for backward compat)
+ dropout_rate = params.get("dropout_rate", 0.0)
+ num_tests = params.get("num_tests", 1)
+ mc_method = params.get("mc_method", "None")
+ cost_per = params.get("cost_per", 0.0)
+ recruitment_rate = params.get("recruitment_rate", 0.0)
+
+ # Apply multiple testing correction
+ alpha_raw = alpha
+ if num_tests > 1 and mc_method != "None":
+ if mc_method == "Bonferroni" or mc_method == "Holm-Bonferroni":
+ alpha = alpha / num_tests
+ elif mc_method == "Benjamini-Hochberg (FDR)":
+ alpha = alpha * (num_tests + 1) / (2 * num_tests)
+ if alternative == "two-sided":
+ z_alpha = norm.ppf(1 - alpha / 2)
+ else:
+ z_alpha = norm.ppf(1 - alpha)
+
+ # --- Computation ---
+ result = None
+ n_total = None
+ n_total_raw = None
+ n_per_group = None
+ explanation = ""
+ formula_latex = ""
+ computed_value = None
+ computed_label = ""
+ _skip_computation = False
+
+ # Compromise mode: generic solver using z-approximation for all test types
+ if analysis_mode == "Compromise":
+ from scipy.optimize import brentq
+
+ n_total = params.get("n_total", 0)
+ if n_total > 0:
+ es_val = params.get("effect_size")
+ n_eff = n_total
+ if atype == "one_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ es_val = abs(
+ proportion_effectsize(params["prop_alt"], params["prop_null"])
+ )
+ elif atype == "two_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ es_val = abs(proportion_effectsize(params["p2"], params["p1"]))
+ ratio = params.get("ratio", 1.0)
+ n_eff = int(n_total / (1 + ratio)) if ratio > 0 else int(n_total / 2)
+ elif atype == "mannwhitney":
+ es_val = np.sqrt(3) * (params["effect_size"] - 0.5)
+ ratio = params.get("ratio", 1.0)
+ n_eff = int(n_total / (1 + ratio)) if ratio > 0 else int(n_total / 2)
+ elif atype == "regression":
+ es_val = params.get("effect_size", 0.15)
+ elif atype == "logistic":
+ ev_cr = params.get("event_rate", 0.3)
+ or_cr = params.get("or", 2.0)
+ p1_cr = (or_cr * ev_cr) / (1 - ev_cr + or_cr * ev_cr)
+ d_cr = abs(p1_cr - ev_cr)
+ p_bar_cr = (ev_cr + p1_cr) / 2
+ se_cr = np.sqrt(2 * p_bar_cr * (1 - p_bar_cr))
+ es_val = d_cr / se_cr if se_cr > 0 else 0
+
+ if es_val is not None and es_val > 0:
+
+ def f(alpha_candidate):
+ z_a = (
+ norm.ppf(1 - alpha_candidate / 2)
+ if alternative == "two-sided"
+ else norm.ppf(1 - alpha_candidate)
+ )
+ z_b = es_val * np.sqrt(n_eff) - z_a
+ achieved = norm.cdf(z_b)
+ beta = max(1e-10, 1 - achieved)
+ return beta / alpha_candidate - params.get("cost_ratio", 1.0)
+
+ lo, hi = 1e-8, 0.5
+ try:
+ if f(lo) * f(hi) < 0:
+ alpha_solved = brentq(f, lo, hi)
+ z_a_solved = (
+ norm.ppf(1 - alpha_solved / 2)
+ if alternative == "two-sided"
+ else norm.ppf(1 - alpha_solved)
+ )
+ z_b_solved = es_val * np.sqrt(n_eff) - z_a_solved
+ power_solved = norm.cdf(z_b_solved)
+ computed_value = {"alpha": alpha_solved, "power": power_solved}
+ n_per_group = n_total
+ _skip_computation = True
+ except Exception:
+ pass
+
+ if computed_value is None:
+ explanation = "Compromise analysis could not converge with current parameters. Try adjusting N or cost ratio."
+ else:
+ explanation = (
+ f"Compromise power analysis with N = {n_total}, "
+ f"effect size = {es_val:.4f}, cost ratio β/α = {params.get('cost_ratio', 1.0):.2f}. "
+ f"Adjusted α = {computed_value['alpha']:.4f}, achieved power = {computed_value['power']:.1%}."
+ )
+ formula_latex = r"\text{Find } \alpha, \beta \text{ s.t. } \beta/\alpha = q \text{ via } \Phi^{-1}"
+
+ if _skip_computation:
+ pass
+ elif atype == "one_mean":
+ d = params["effect_size"]
+ from statsmodels.stats.power import TTestPower
+ solver = TTestPower()
+ if d > 0:
+ if analysis_mode == "Compromise":
+ n_total = params.get("n_total", 0)
+ n_per_group = n_total
+ comp = _solve_compromise(
+ solver,
+ d,
+ n_total,
+ params.get("cost_ratio", 1.0),
+ alternative=alternative,
+ )
+ computed_value = comp
+ else:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=d,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n_total = int(np.ceil(raw))
+ n_per_group = n_total
+ else:
+ n_total = params.get("n_total", 0)
+ n_per_group = n_total
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required total sample size for a one-sample {tails.lower()} t-test "
+ f"to detect Cohen's d = {d:.3f} with α = {alpha} and power = {power}."
+ )
+ formula_latex = (
+ rf"n = \left( \frac{{{_z_sub} + z_{{1-\beta}}}}{{d}} \right)^2"
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for a one-sample {tails.lower()} t-test "
+ f"with N = {n_total}, Cohen's d = {d:.3f}, α = {alpha}."
+ )
+ formula_latex = rf"\text{{Power}} = \int_{{t_\text{{crit}}}}^{{\infty}} f(t; n-1, \delta) dt"
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable effect size for a one-sample {tails.lower()} t-test "
+ f"with N = {n_total}, α = {alpha}, power = {power}."
+ )
+ formula_latex = (
+ rf"d \approx \frac{{{_z_sub} + z_{{1-\beta}}}}{{\sqrt{{n}}}}"
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for a one-sample {tails.lower()} t-test "
+ f"with N = {n_total}, Cohen's d = {d:.3f}, power = {power}."
+ )
+ formula_latex = rf"\alpha = \int_{{t_\text{{crit}}}}^{{\infty}} f(t; n-1, \delta = 0) dt \quad \text{{where }} t_\text{{crit}} = t_{{1-\alpha, n-1}}"
+ elif analysis_mode == "Compromise" and isinstance(computed_value, dict):
+ explanation = (
+ f"Compromise power analysis for a one-sample {tails.lower()} t-test "
+ f"with N = {n_total}, Cohen's d = {d:.3f}, "
+ f"cost ratio β/α = {params.get('cost_ratio', 1.0):.2f}. "
+ f"Adjusted α = {computed_value['alpha']:.4f}, achieved power = {computed_value['power']:.1%}."
+ )
+ formula_latex = rf"""
+ \beta = q \cdot \alpha \quad \text{{Find }} \alpha, \beta \text{ via } \Phi^{-1}
+ \alpha = 2 \left(1 - \Phi\left(c)\right)
+ \beta = \Phi \left(c - \delta \right) -\Phi\left(-c - \delta \right)
+ \delta = d \sqrt{{n}} \quad c = z_{{1-\alpha/2}} \text{{ or }} z_{{1-\alpha}}
"""
-# Welcome to Streamlit!
-Edit `/streamlit_app.py` to customize this app to your heart's desire :heart:.
-If you have any questions, checkout our [documentation](https://docs.streamlit.io) and [community
-forums](https://discuss.streamlit.io).
+ elif atype == "two_means":
+ d = params["effect_size"]
+ ratio = params["ratio"]
+ from statsmodels.stats.power import TTestIndPower
+
+ solver = TTestIndPower()
+ if d > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=d,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ nobs_name="nobs1",
+ ratio=ratio,
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n1 = int(np.ceil(raw))
+ n2 = int(np.ceil(n1 * ratio))
+ n_total = n1 + n2
+ n_per_group = (n1, n2)
+ else:
+ n_total = params.get("n_total", 0)
+ n1 = int(np.ceil(n_total / (1 + ratio))) if ratio > 0 else 0
+ n2 = n_total - n1
+ n_per_group = (n1, n2)
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required sample size per group for an independent {tails.lower()} t-test "
+ f"to detect Cohen's d = {d:.3f} with α = {alpha} and power = {power}, "
+ f"allocation ratio n₂/n₁ = {ratio:.2f}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for an independent {tails.lower()} t-test "
+ f"with N = {n_total}, Cohen's d = {d:.3f}, α = {alpha}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable effect size for an independent {tails.lower()} t-test "
+ f"with N = {n_total}, α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for an independent {tails.lower()} t-test "
+ f"with N = {n_total}, Cohen's d = {d:.3f}, power = {power}, ratio = {ratio:.2f}."
+ )
+ formula_latex = rf"n_1 = 2 \left( \frac{{{_z_sub} + z_{{1-\beta}}}}{{d}} \right)^2 \quad n_2 = r \cdot n_1"
+
+ elif atype == "paired":
+ d = params["effect_size"]
+ from statsmodels.stats.power import TTestPower
+
+ solver = TTestPower()
+ if d > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=d,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n_total = int(np.ceil(raw))
+ n_per_group = n_total
+ else:
+ n_total = params.get("n_total", 0)
+ n_per_group = n_total
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required number of pairs for a paired {tails.lower()} t-test "
+ f"to detect Cohen's d_z = {d:.3f} with α = {alpha} and power = {power}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for a paired {tails.lower()} t-test "
+ f"with {n_total} pairs, Cohen's d_z = {d:.3f}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable effect size for a paired {tails.lower()} t-test "
+ f"with {n_total} pairs, α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for a paired {tails.lower()} t-test "
+ f"with {n_total} pairs, Cohen's d_z = {d:.3f}, power = {power}."
+ )
+ formula_latex = (
+ rf"n = \left( \frac{{{_z_sub} + z_{{1-\beta}}}}{{d_z}} \right)^2"
+ )
+
+ elif atype == "one_prop":
+ p0 = params["prop_null"]
+ p1 = params["prop_alt"]
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ d_eff = proportion_effectsize(p1, p0)
+ solver = NormalIndPower()
+ if abs(d_eff) > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=abs(d_eff),
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ nobs_name="nobs1",
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n_total = int(np.ceil(raw))
+ n_per_group = n_total
+ else:
+ n_total = params.get("n_total", 0)
+ n_per_group = n_total
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required sample size for a one-sample proportion test "
+ f"to detect a difference from {p0} to {p1} "
+ f"with α = {alpha} and power = {power}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for a one-sample proportion test "
+ f"with N = {n_total}, difference {p0} → {p1}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable proportion difference "
+ f"for a one-sample proportion test with N = {n_total}, α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for a one-sample proportion test "
+ f"with N = {n_total}, difference {p0} → {p1}, power = {power}."
+ )
+ formula_latex = rf"n = \left( \frac{{{_z_sub} \sqrt{{p_0(1-p_0)}} + z_{{1-\beta}} \sqrt{{p_1(1-p_1)}}}}{{{{p_1 - p_0}}}} \right)^2"
+
+ elif atype == "two_prop":
+ p1 = params["p1"]
+ p2 = params["p2"]
+ ratio = params["ratio"]
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ d_eff = proportion_effectsize(p2, p1)
+ solver = NormalIndPower()
+ if abs(d_eff) > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=abs(d_eff),
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ nobs_name="nobs1",
+ ratio=ratio,
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n1 = int(np.ceil(raw))
+ n2 = int(np.ceil(n1 * ratio))
+ n_total = n1 + n2
+ n_per_group = (n1, n2)
+ else:
+ n_total = params.get("n_total", 0)
+ n1 = int(np.ceil(n_total / (1 + ratio))) if ratio > 0 else 0
+ n2 = n_total - n1
+ n_per_group = (n1, n2)
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required sample size per group for a two-proportion z-test "
+ f"to detect a difference between {p1} and {p2} "
+ f"with α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for a two-proportion z-test "
+ f"with N = {n_total}, proportions {p1} vs {p2}, α = {alpha}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable proportion difference "
+ f"for a two-proportion z-test with N = {n_total}, α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for a two-proportion z-test "
+ f"with N = {n_total}, proportions {p1} vs {p2}, power = {power}, ratio = {ratio:.2f}."
+ )
+ formula_latex = rf"n_1 = \left( \frac{{{_z_sub} \sqrt{{2\bar{{p}}(1-\bar{{p}})}} + z_{{1-\beta}} \sqrt{{p_1(1-p_1) + p_2(1-p_2)}}}}{{{{p_1 - p_2}}}} \right)^2"
+
+ elif atype == "anova":
+ f_eff = params["effect_size"]
+ k = params["k"]
+ from statsmodels.stats.power import FTestAnovaPower
+
+ solver = FTestAnovaPower()
+ if f_eff > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=f_eff,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ k_groups=k,
+ )
+ if analysis_mode == "A Priori":
+ n_per_g = int(np.ceil(raw))
+ n_total = n_per_g * k
+ n_per_group = n_per_g
+ else:
+ n_total = params.get("n_total", 0)
+ n_per_g = max(1, n_total // k)
+ n_per_group = n_per_g
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required sample size per group for a one-way ANOVA with {k} groups "
+ f"to detect Cohen's f = {f_eff:.3f} with α = {alpha} and power = {power}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for a one-way ANOVA with {k} groups "
+ f"with N = {n_total}, Cohen's f = {f_eff:.3f}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable effect size for a one-way ANOVA with {k} groups "
+ f"with N = {n_total}, α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for a one-way ANOVA with {k} groups "
+ f"with N = {n_total}, Cohen's f = {f_eff:.3f}, power = {power}."
+ )
+ formula_latex = r"n = \frac{\text{from non-central }F\text{ distribution}}{k} \quad f = \frac{\sigma_{\text{between}}}{\sigma_{\text{within}}}"
+
+ elif atype == "correlation":
+ r_val = params["effect_size"]
+ import math as cmath
+
+ fisher_z = cmath.atanh(r_val)
+ if analysis_mode == "A Priori":
+ n_total = int(np.ceil(3 + ((z_alpha + z_beta) / fisher_z) ** 2))
+ elif analysis_mode == "Post Hoc":
+ n_total = params.get("n_total", 0)
+ z_beta_c = (cmath.sqrt(max(0, n_total - 3)) * fisher_z) - z_alpha
+ computed_value = norm.cdf(z_beta_c)
+ elif analysis_mode == "Sensitivity":
+ n_total = params.get("n_total", 0)
+ z_beta_c = norm.ppf(power)
+ fisher_z_c = (z_alpha + z_beta_c) / cmath.sqrt(max(1, n_total - 3))
+ computed_value = cmath.tanh(fisher_z_c)
+ elif analysis_mode == "Criterion":
+ n_total = params.get("n_total", 0)
+ z_beta_c = norm.ppf(power)
+ z_alpha_c = cmath.sqrt(max(0, n_total - 3)) * fisher_z - z_beta_c
+ computed_value = 2 * (1 - norm.cdf(z_alpha_c))
+ n_per_group = n_total
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required sample size to detect a Pearson correlation of r = {r_val:.3f} "
+ f"with α = {alpha} and power = {power} ({tails.lower()}), "
+ f"based on Fisher's z-transformation."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power to detect Pearson correlation r = {r_val:.3f} "
+ f"with N = {n_total}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable Pearson correlation "
+ f"with N = {n_total}, α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level to detect Pearson correlation r = {r_val:.3f} "
+ f"with N = {n_total}, power = {power}."
+ )
+ formula_latex = rf"n = 3 + \left( \frac{{{_z_sub} + z_{{1-\beta}}}}{{\text{{arctanh}}(r)}} \right)^2"
+
+ elif atype == "regression":
+ f2 = params["effect_size"]
+ k = params["k"]
+ if f2 > 0:
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ if analysis_mode == "A Priori":
+ for n_candidate in range(k + 2, 10000):
+ dfd = n_candidate - k - 1
+ ncp = f2 * n_candidate
+ f_crit = f_dist.ppf(1 - alpha, k, dfd)
+ pwr_cur = 1 - noncentral_f.cdf(f_crit, k, dfd, ncp)
+ if pwr_cur >= power:
+ n_total = n_candidate
+ break
+ elif analysis_mode == "Post Hoc":
+ n_total = params.get("n_total", 0)
+ dfd = n_total - k - 1
+ ncp = f2 * n_total
+ f_crit = f_dist.ppf(1 - alpha, k, dfd)
+ computed_value = 1 - noncentral_f.cdf(f_crit, k, dfd, ncp)
+ elif analysis_mode == "Sensitivity":
+ n_total = params.get("n_total", 0)
+ from scipy.optimize import brentq
+
+ def power_for_f2(f2_try):
+ dfd = n_total - k - 1
+ ncp = f2_try * n_total
+ f_crit = f_dist.ppf(1 - alpha, k, dfd)
+ return 1 - noncentral_f.cdf(f_crit, k, dfd, ncp) - power
+
+ computed_value = brentq(power_for_f2, 1e-6, 5.0)
+ elif analysis_mode == "Criterion":
+ n_total = params.get("n_total", 0)
+ from scipy.optimize import brentq
+
+ def power_for_alpha(a_try):
+ dfd = n_total - k - 1
+ ncp = f2 * n_total
+ f_crit = f_dist.ppf(1 - a_try, k, dfd)
+ return 1 - noncentral_f.cdf(f_crit, k, dfd, ncp) - power
+
+ computed_value = brentq(power_for_alpha, 1e-8, 0.5)
+ n_per_group = n_total
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required total sample size for multiple linear regression "
+ f"with {k} predictor(s) to detect Cohen's f² = {f2:.3f} "
+ f"(R² = {f2 / (1 + f2):.3f}) with α = {alpha} and power = {power}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for multiple linear regression "
+ f"with {k} predictor(s), N = {n_total}, Cohen's f² = {f2:.3f}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable Cohen's f² for multiple linear regression "
+ f"with {k} predictor(s), N = {n_total}, α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for multiple linear regression "
+ f"with {k} predictor(s), N = {n_total}, Cohen's f² = {f2:.3f}, power = {power}."
+ )
+ formula_latex = r"n = \text{from non-central }F\text{ distribution} \quad f^2 = \frac{R^2}{1-R^2}"
+
+ elif atype == "logistic":
+ k = params["k"]
+ ev_rate = params["event_rate"]
+ or_val = params["or"]
+ p1_log = (or_val * ev_rate) / (1 - ev_rate + or_val * ev_rate)
+ d_log = abs(p1_log - ev_rate)
+ p_bar = (ev_rate + p1_log) / 2
+ from statsmodels.stats.power import NormalIndPower
+
+ solver = NormalIndPower()
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ d_eff_log = d_log / se if se > 0 else 0
+ if d_eff_log > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=d_eff_log,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ nobs_name="nobs1",
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n_base = int(np.ceil(raw))
+ n_total = max(n_base, 10 * k)
+ else:
+ n_total = params.get("n_total", 0)
+ n_per_group = n_total
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required total sample size for logistic regression "
+ f"with {k} predictor(s) to detect OR = {or_val:.2f} "
+ f"with baseline event rate = {ev_rate:.2f}, α = {alpha}, power = {power}. "
+ f"Lower bound of 10 × {k} = {10 * k} events per predictor applied."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for logistic regression "
+ f"with {k} predictor(s), N = {n_total}, OR = {or_val:.2f}, "
+ f"event rate = {ev_rate:.2f}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable odds ratio for logistic regression "
+ f"with {k} predictor(s), N = {n_total}, event rate = {ev_rate:.2f}, "
+ f"α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for logistic regression "
+ f"with {k} predictor(s), N = {n_total}, OR = {or_val:.2f}, "
+ f"event rate = {ev_rate:.2f}, power = {power}."
+ )
+ formula_latex = rf"n = \frac{{({_z_sub} + z_{{1-\beta}})^2 \bar{{p}}(1-\bar{{p}})}}{{{{(p_1 - p_0)^2}}}} \quad \text{{min }} 10k"
+
+ elif atype == "chisq":
+ w = params["effect_size"]
+ df = params["df"]
+ from statsmodels.stats.power import GofChisquarePower
+
+ solver = GofChisquarePower()
+ if w > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=w,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ nobs_name="nobs",
+ n_bins=df + 1,
+ )
+ if analysis_mode == "A Priori":
+ n_total = int(np.ceil(raw))
+ else:
+ n_total = params.get("n_total", 0)
+ n_per_group = n_total
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required total sample size for a chi-square test "
+ f"with {df} degree(s) of freedom to detect Cohen's w = {w:.3f} "
+ f"with α = {alpha} and power = {power}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for a chi-square test "
+ f"with N = {n_total}, {df} DF, Cohen's w = {w:.3f}, α = {alpha}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable Cohen's w for a chi-square test "
+ f"with N = {n_total}, {df} DF, α = {alpha}, power = {power}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for a chi-square test "
+ f"with N = {n_total}, {df} DF, Cohen's w = {w:.3f}, power = {power}."
+ )
+ formula_latex = r"n = \text{from non-central }\chi^2\text{ distribution} \quad w = \sqrt{\sum \frac{(p_{0i} - p_{1i})^2}{p_{0i}}}"
+
+ elif atype == "mannwhitney":
+ p_val = params["effect_size"]
+ ratio = params["ratio"]
+ are = params["are"]
+ from statsmodels.stats.power import NormalIndPower
+
+ solver = NormalIndPower()
+ d_mw = np.sqrt(3) * (p_val - 0.5)
+ if d_mw > 0:
+ raw = _solve_with_mode(
+ solver,
+ analysis_mode,
+ effect_size=d_mw,
+ nobs=params.get("n_total"),
+ alpha=alpha,
+ power=power,
+ nobs_name="nobs1",
+ ratio=ratio,
+ alternative=alternative,
+ )
+ if analysis_mode == "A Priori":
+ n1 = int(np.ceil(raw / are))
+ n2 = int(np.ceil(n1 * ratio))
+ n_total = n1 + n2
+ n_per_group = (n1, n2)
+ else:
+ n_total = params.get("n_total", 0)
+ n1 = int(np.ceil(n_total / (1 + ratio))) if ratio > 0 else 0
+ n2 = n_total - n1
+ n_per_group = (n1, n2)
+ computed_value = raw
+ if analysis_mode == "A Priori":
+ explanation = (
+ f"Required sample size for Mann-Whitney / Wilcoxon test "
+ f"to detect P(X>Y) = {p_val:.3f} (d ≈ {d_mw:.3f}) with ARE = {are:.3f}, "
+ f"α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Post Hoc":
+ explanation = (
+ f"Achieved power for Mann-Whitney / Wilcoxon test "
+ f"with N = {n_total}, P(X>Y) = {p_val:.3f}, α = {alpha}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Sensitivity":
+ explanation = (
+ f"Minimum detectable P(X>Y) for Mann-Whitney / Wilcoxon test "
+ f"with N = {n_total}, α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ elif analysis_mode == "Criterion":
+ explanation = (
+ f"Required significance level for Mann-Whitney / Wilcoxon test "
+ f"with N = {n_total}, P(X>Y) = {p_val:.3f}, power = {power}, ratio = {ratio:.2f}."
+ )
+ formula_latex = r"n_{\text{nonparam}} = \frac{n_{\text{param}}}{\text{ARE}} \quad \text{ARE} \approx 0.955"
+
+ elif atype == "logrank":
+ hr = params["hr"]
+ ratio = params["ratio"]
+ med_ctrl = params["median_survival"]
+ study_dur = params["study_duration"]
+ log_hr = np.log(hr)
+ num_events = (
+ ((z_alpha + z_beta) ** 2) * ((ratio + 1) ** 2) / (ratio * (log_hr**2))
+ )
+ lambda_ctrl = np.log(2) / med_ctrl
+ p_event_ctrl = 1 - np.exp(-lambda_ctrl * study_dur)
+ p_event_trt = 1 - np.exp(-lambda_ctrl / hr * study_dur)
+ p_event = (p_event_ctrl + ratio * p_event_trt) / (1 + ratio)
+ if p_event > 0:
+ n_total = int(np.ceil(num_events / p_event))
+ n1 = int(np.ceil(n_total / (1 + ratio)))
+ n2 = n_total - n1
+ n_per_group = (n1, n2)
+ explanation = (
+ f"Required sample size for log-rank test to detect HR = {hr:.2f} "
+ f"with median survival {med_ctrl:.0f} months, study duration {study_dur:.0f} months, "
+ f"α = {alpha}, power = {power}, ratio = {ratio:.2f}. "
+ f"Events needed: {int(np.ceil(num_events))}."
+ )
+ formula_latex = rf"E = \frac{{({_z_sub}+z_{{1-\beta}})^2 (r+1)^2}}{{r (\log HR)^2}} \quad N = \frac{{E}}{{P(\text{{event}})}}"
+
+ elif atype == "cox":
+ hr = params["hr"]
+ k = params["k"]
+ sd_x = params["sd_x"]
+ r2_x = params["r2_x"]
+ ev_rate = params["event_rate"]
+ log_hr = np.log(hr)
+ var_denom = (sd_x**2) * (log_hr**2) * (1 - r2_x)
+ if var_denom > 0:
+ num_events = ((z_alpha + z_beta) ** 2) / var_denom
+ num_events = max(num_events, 10 * k)
+ n_total = int(np.ceil(num_events / ev_rate))
+ n_per_group = n_total
+ explanation = (
+ f"Required total sample size for Cox regression with {k} predictor(s) "
+ f"to detect HR = {hr:.2f} (SD = {sd_x:.1f}, R² = {r2_x:.2f}) "
+ f"with α = {alpha}, power = {power}, event rate = {ev_rate:.2f}. "
+ f"Events needed: {int(np.ceil(num_events))}."
+ )
+ formula_latex = rf"E = \frac{{({_z_sub}+z_{{1-\beta}})^2}}{{\sigma_x^2 \beta^2 (1-R^2)}} \quad N = \frac{{E}}{{P(\text{{event}})}}"
+
+ elif atype == "equiv":
+ margin = params["margin"]
+ exp_diff = params["expected_diff"]
+ sd = params["sd"]
+ ratio = params["ratio"]
+ equiv_type = params.get("equiv_param_type", "Mean")
+ if equiv_type == "Proportion":
+ p1_eq = params.get("p1_eq", 0.2)
+ p2_eq = params.get("p2_eq", 0.2)
+ d_prop = abs(p1_eq - p2_eq)
+ d_e = margin - d_prop
+ if d_e > 0:
+ p_bar = (p1_eq + p2_eq) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ es_e = d_e / se if se > 0 else 0
+ from statsmodels.stats.power import NormalIndPower
+
+ n1 = NormalIndPower().solve_power(
+ effect_size=es_e,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative="larger",
+ )
+ n1 = int(np.ceil(n1))
+ n2 = int(np.ceil(n1 * ratio))
+ n_total = n1 + n2
+ n_per_group = (n1, n2)
+ explanation = (
+ f"Required sample size for non-inferiority test of proportions "
+ f"with margin = {margin:.3f}, p₁ = {p1_eq:.3f}, p₂ = {p2_eq:.3f}, "
+ f"expected difference = {d_prop:.3f}, α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ formula_latex = r"n_1 = \frac{(z_{1-\alpha}+z_{1-\beta})^2 \, 2\bar{p}(1-\bar{p})}{(\delta - |p_1-p_2|)^2}"
+
+ else:
+ d_e = margin - abs(exp_diff)
+ if d_e > 0:
+ from statsmodels.stats.power import NormalIndPower
+
+ es_e = d_e / sd
+ n1 = NormalIndPower().solve_power(
+ effect_size=es_e,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative="larger",
+ )
+ n1 = int(np.ceil(n1))
+ n2 = int(np.ceil(n1 * ratio))
+ n_total = n1 + n2
+ n_per_group = (n1, n2)
+ explanation = (
+ f"Required sample size for non-inferiority test of means "
+ f"with margin = {margin:.2f}, expected difference = {exp_diff:.2f}, "
+ f"SD = {sd:.1f}, α = {alpha}, power = {power}, ratio = {ratio:.2f}."
+ )
+ formula_latex = r"n_1 = \frac{(z_{1-\alpha}+z_{1-\beta})^2 \sigma^2 (1+1/r)}{(\delta - |d|)^2}"
+
+ elif atype == "rm_anova":
+ f_eff = params["effect_size"]
+ k = params["k"]
+ m = params["m"]
+ rho = params["rho"]
+ epsilon = params["epsilon"]
+ from statsmodels.stats.power import FTestAnovaPower
+
+ solver = FTestAnovaPower()
+ if f_eff > 0:
+ n_per_g = solver.solve_power(
+ effect_size=f_eff,
+ alpha=alpha,
+ power=power,
+ k_groups=k,
+ )
+ design_effect = (1 + (m - 1) * rho) / m
+ df_adj = (m - 1) * epsilon
+ n_per_g_adj = int(
+ np.ceil(n_per_g * design_effect * k / (k * df_adj / (k - 1)))
+ )
+ n_per_g_adj = max(n_per_g_adj, int(np.ceil(n_per_g)))
+ n_total = n_per_g_adj * k
+ n_per_group = n_per_g_adj
+ explanation = (
+ f"Required sample size per group for repeated measures ANOVA "
+ f"with {k} group(s), {m} measurement(s), ρ = {rho:.2f}, ε = {epsilon:.2f}, "
+ f"Cohen's f = {f_eff:.3f}, α = {alpha}, power = {power}."
+ )
+ formula_latex = r"\text{Adjustment: } n_{\text{adj}} = n \times \frac{1+(m-1)\rho}{m} \quad \text{GG: } \varepsilon"
+
+ elif atype == "twoway_anova":
+ f_a = params["f_a"]
+ f_b = params["f_b"]
+ f_ab = params["f_ab"]
+ rows = params["rows"]
+ cols = params["cols"]
+ focus = params["focus"]
+ from statsmodels.stats.power import FTestAnovaPower
+
+ solver = FTestAnovaPower()
+ if focus == "Main Effect A":
+ f_use = f_a
+ k_use = rows
+ elif focus == "Main Effect B":
+ f_use = f_b
+ k_use = cols
+ else:
+ f_use = f_ab
+ k_use = rows * cols
+ if f_use > 0:
+ n_per_cell = solver.solve_power(
+ effect_size=f_use,
+ alpha=alpha,
+ power=power,
+ k_groups=k_use,
+ )
+ n_per_cell = int(np.ceil(n_per_cell))
+ n_total = n_per_cell * rows * cols
+ n_per_group = n_per_cell
+ explanation = (
+ f"Required sample size per cell for two-way ANOVA "
+ f"({rows} × {cols} design), focus on {focus}, "
+ f"Cohen's f = {f_use:.3f}, α = {alpha}, power = {power}. "
+ f"Total N = {n_total} across {rows * cols} cells."
+ )
+ formula_latex = r"n_{\text{per cell}} = \text{from non-central }F \quad N = n \times r \times c"
+
+ elif atype == "roc_auc":
+ auc = params["auc"]
+ null_auc = params["null_auc"]
+ ratio = params["ratio"]
+ v_auc = auc * (1 - auc) + (ratio - 1) * (auc / (2 - auc) - auc**2) / (1 + ratio)
+ delta_auc = auc - null_auc
+ if delta_auc > 0:
+ n_cases = ((z_alpha + z_beta) ** 2 * v_auc) / (delta_auc**2)
+ n_cases = int(np.ceil(n_cases))
+ n_controls = int(np.ceil(n_cases * ratio))
+ n_total = n_cases + n_controls
+ n_per_group = (n_cases, n_controls)
+ explanation = (
+ f"Required sample size for ROC/AUC analysis to detect AUC = {auc:.3f} "
+ f"(null = {null_auc:.1f}) with case:control ratio {ratio:.2f}, "
+ f"α = {alpha}, power = {power}."
+ )
+ formula_latex = rf"n_{{\text{{cases}}}} = \frac{{({_z_sub}+z_{{1-\beta}})^2 V(AUC)}}{{{{(AUC - 0.5)^2}}}}"
+
+ elif atype == "kappa":
+ kappa = params["kappa"]
+ null_kappa = params["null_kappa"]
+ raters = params["raters"]
+ cats = params["categories"]
+ delta_k = kappa - null_kappa
+ if delta_k > 0:
+ n_total = int(
+ np.ceil(
+ ((z_alpha + z_beta) ** 2 * null_kappa * (1 - null_kappa))
+ / (delta_k**2)
+ )
+ )
+ n_total = max(n_total, raters * cats * 5)
+ n_per_group = n_total
+ explanation = (
+ f"Required sample size for Cohen's Kappa with {raters} rater(s), {cats} category(ies), "
+ f"expected κ = {kappa:.3f}, null κ = {null_kappa:.2f}, "
+ f"α = {alpha}, power = {power}."
+ )
+ formula_latex = rf"n = \frac{{({_z_sub}+z_{{1-\beta}})^2 \kappa_0(1-\kappa_0)}}{{(\kappa - \kappa_0)^2}}"
+
+ elif atype == "cluster_rct":
+ d = params["effect_size"]
+ icc = params["icc"]
+ cluster_m = params["cluster_size"]
+ ratio = params["ratio"]
+ deff = 1 + (cluster_m - 1) * icc
+ from statsmodels.stats.power import TTestIndPower
+
+ solver = TTestIndPower()
+ if d > 0:
+ n1_ind = solver.solve_power(
+ effect_size=d,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1_ind = int(np.ceil(n1_ind))
+ n1_clust = int(np.ceil(n1_ind * deff))
+ n1_clust = int(np.ceil(n1_clust / cluster_m)) * cluster_m
+ n2_clust = int(np.ceil(n1_clust * ratio))
+ n_total = n1_clust + n2_clust
+ num_clusters_1 = int(np.ceil(n1_clust / cluster_m))
+ num_clusters_2 = int(np.ceil(n2_clust / cluster_m))
+ n_per_group = (n1_clust, n2_clust)
+ explanation = (
+ f"Required sample size for cluster RCT with ICC = {icc:.3f}, "
+ f"cluster size m = {cluster_m}, DEFF = {deff:.2f}, "
+ f"Cohen's d = {d:.3f}, α = {alpha}, power = {power}, ratio = {ratio:.2f}. "
+ f"Clusters needed: {num_clusters_1} in group 1, {num_clusters_2} in group 2."
+ )
+ formula_latex = (
+ r"N_{\text{clust}} = N_{\text{ind}} \times [1+(m-1)\rho] \quad \text{DEFF}"
+ )
+
+ elif atype == "precision":
+ half_width = params["half_width"]
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and half_width > 0:
+ from scipy.stats import t as t_dist
+
+ n_total = int(
+ np.ceil((norm.ppf(1 - conf_alpha / 2) * sd / half_width) ** 2)
+ )
+ n_total = max(n_total, 3)
+ for _ in range(20):
+ t_val = t_dist.ppf(1 - conf_alpha / 2, df=n_total - 1)
+ n_next = int(np.ceil((t_val * sd / half_width) ** 2))
+ n_next = max(n_next, 3)
+ if n_next == n_total:
+ break
+ n_total = n_next
+ else:
+ n_total = 3
+ else:
+ prop = params["prop"]
+ if prop > 0 and half_width > 0:
+ n_total = int(
+ np.ceil(
+ (norm.ppf(1 - conf_alpha / 2) ** 2 * prop * (1 - prop))
+ / (half_width**2)
+ )
+ )
+ else:
+ n_total = 3
+ n_total = max(n_total, 3)
+ n_per_group = n_total
+ explanation = (
+ f"Required sample size for precision-based estimation of a {param_type.lower()} "
+ f"with {conf_level:.0f}% CI half-width = {half_width:.2f}. "
+ f"σ = {sd:.2f}"
+ if param_type == "Mean"
+ else f"p = {prop:.3f}."
+ )
+ formula_latex = (
+ r"n = \left(\frac{t_{1-\alpha/2, n-1} \sigma}{w}\right)^2"
+ if param_type == "Mean"
+ else r"n = \frac{z_{1-\alpha/2}^2 p(1-p)}{w^2}"
+ )
+
+ elif atype == "pilot":
+ method = params["method"]
+ if method == "Rule of thumb":
+ n_total = params["n_per_group"]
+ n_per_group = n_total
+ elif method == "Precision-based":
+ half_width = params["half_width"]
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and half_width > 0:
+ from scipy.stats import t as t_dist
+
+ n_total = int(
+ np.ceil((norm.ppf(1 - conf_alpha / 2) * sd / half_width) ** 2)
+ )
+ n_total = max(n_total, 3)
+ for _ in range(20):
+ t_val = t_dist.ppf(1 - conf_alpha / 2, df=n_total - 1)
+ n_next = int(np.ceil((t_val * sd / half_width) ** 2))
+ n_next = max(n_next, 3)
+ if n_next == n_total:
+ break
+ n_total = n_next
+ else:
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ prop = params["prop"]
+ if prop > 0 and half_width > 0:
+ n_total = int(
+ np.ceil((z_hw**2 * prop * (1 - prop)) / (half_width**2))
+ )
+ n_total = max(n_total, 5)
+ n_per_group = n_total
+ else:
+ main_n = params["main_n"]
+ fraction = params["fraction"]
+ n_total = max(int(np.ceil(main_n * fraction)), 5)
+ n_per_group = n_total
+ explanation = (
+ f"Pilot/feasibility study sample size using '{method}' method. "
+ f"N = {n_total} per group. Recommended for estimating parameters "
+ f"for a future definitive trial."
+ )
+ formula_latex = r"n_{\text{pilot}} = \text{rule of thumb: } 12\text{/group}"
+
+ elif atype == "wilcoxon_sr":
+ p_val = params["effect_size"]
+ are = params["are"]
+ d_wsr = np.sqrt(3) * (p_val - 0.5) * 2
+ if d_wsr > 0:
+ from statsmodels.stats.power import NormalIndPower
+
+ n_total = int(
+ np.ceil(
+ NormalIndPower().solve_power(
+ effect_size=d_wsr,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ / are
+ )
+ )
+ n_per_group = n_total
+ explanation = (
+ f"Required sample size for a paired Wilcoxon signed-rank test ({tails.lower()}) "
+ f"to detect Pr(positive diff) = {p_val:.3f} with α = {alpha}, power = {power}, ARE = {are:.3f}."
+ )
+ formula_latex = rf"n = \frac{{1}}{{\text{{ARE}}}} \left( \frac{{{_z_sub} + z_{{1-\beta}}}}{{\sqrt{{3}} \cdot (P - 0.5)}} \right)^2"
+
+ elif atype == "kruskal":
+ f_eff = params["effect_size"]
+ k = params["k"]
+ are_kw = params.get("are", 0.955)
+ if f_eff > 0:
+ from statsmodels.stats.power import FTestAnovaPower
+
+ n_per_g = FTestAnovaPower().solve_power(
+ effect_size=f_eff, alpha=alpha, power=power, k_groups=k
+ )
+ n_per_g = int(np.ceil(n_per_g / are_kw))
+ n_total = n_per_g * k
+ n_per_group = n_per_g
+ explanation = (
+ f"Required sample size per group for Kruskal-Wallis test with {k} groups "
+ f"to detect Cohen's f = {f_eff:.3f} with α = {alpha}, power = {power}. "
+ f"Inflated by {1/are_kw:.2f}× vs ANOVA (ARE = {are_kw:.3f}) to account for efficiency loss."
+ )
+ formula_latex = r"n_{\text{KW}} = \frac{n_{\text{ANOVA}}}{\text{ARE}} \quad f = \sqrt{\frac{\sum(\bar{R}_i - \bar{R})^2}{(N(N+1)/12)}}"
+
+ elif atype == "friedman":
+ k = params["k"]
+ m = params["m"]
+ w = params["w"]
+ are = params["are"]
+ if w > 0:
+ from statsmodels.stats.power import FTestAnovaPower
+
+ f_fr = np.sqrt(w / (1 - w))
+ n_per_g = FTestAnovaPower().solve_power(
+ effect_size=f_fr, alpha=alpha, power=power, k_groups=k
+ )
+ n_per_g = int(np.ceil(n_per_g * m / are))
+ n_total = n_per_g * k
+ n_per_group = n_per_g
+ explanation = (
+ f"Required sample size per group for Friedman test with {k} groups, {m} measurements, "
+ f"Kendall's W = {w:.3f}, α = {alpha}, power = {power}."
+ )
+ formula_latex = (
+ r"W = \frac{12 \sum R_i^2}{m^2 k(k^2-1)} \quad f = \sqrt{\frac{W}{1-W}}"
+ )
+
+ elif atype == "mcnemar":
+ p_b = params["p_b"]
+ p_c = params["p_c"]
+ d_mc = abs(p_b - p_c)
+ p_discordant = p_b + p_c
+ if d_mc > 0 and p_discordant > 0:
+ n_total = int(np.ceil((z_alpha + z_beta) ** 2 * p_discordant / (d_mc**2)))
+ n_per_group = n_total
+ explanation = (
+ f"Required number of pairs for McNemar's test ({tails.lower()}) "
+ f"to detect discordant proportions p_b = {p_b:.3f}, p_c = {p_c:.3f} "
+ f"with α = {alpha}, power = {power}."
+ )
+ formula_latex = rf"n = \frac{{({_z_sub} + z_{{1-\beta}})^2 (p_b + p_c)}}{{{{(p_b - p_c)^2}}}}"
+
+ elif atype == "fisher":
+ p1 = params["p1"]
+ p2 = params["p2"]
+ ratio = params["ratio"]
+ are_f = params.get("are", 0.833)
+ p_bar = (p1 + ratio * p2) / (1 + ratio)
+ d_f = abs(p1 - p2)
+ if d_f > 0 and p_bar > 0 and p_bar < 1:
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ d_eff = proportion_effectsize(p2, p1)
+ n1 = NormalIndPower().solve_power(
+ effect_size=abs(d_eff),
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1 = int(np.ceil(n1 / are_f))
+ n2 = int(np.ceil(n1 * ratio))
+ n_total = n1 + n2
+ n_per_group = (n1, n2)
+ explanation = (
+ f"Required sample size per group for Fisher's exact test ({tails.lower()}) "
+ f"to detect difference between {p1:.3f} and {p2:.3f} "
+ f"with α = {alpha}, power = {power}, ratio = {ratio:.2f}. "
+ f"Inflated by {1/are_f:.2f}× vs z-test (ARE = {are_f:.3f}) for exact-method efficiency."
+ )
+ formula_latex = r"n_{\text{Fisher}} \approx \frac{n_{\text{two-prop}}}{\text{ARE}} \quad \text{(exact conditional inference)}"
+
+ elif atype == "manova":
+ k = params["k"]
+ dv = params["dv"]
+ f2 = params["f2"]
+ rho = params["rho"]
+ manova_test = params.get("manova_test", "Pillai's Trace")
+ if f2 > 0:
+ u = dv
+ v_num = k - 1
+ v_den = 1e9
+ n_total = None
+ for n_try in range(k * dv + 2, 5000):
+ v_den = n_try - k - dv
+ if v_den <= 0:
+ continue
+ s_val = min(u, v_num)
+ df1 = u * v_num
+ if manova_test == "Pillai's Trace":
+ df2 = s_val * (v_den - dv + 1) + 4
+ elif manova_test == "Wilks' Lambda":
+ t_val = max(
+ np.sqrt((u**2 * v_num**2 - 4) / max(u**2 + v_num**2 - 5, 1)), 1
+ )
+ df2 = (v_den - (u - v_num + 1) / 2) * t_val - (u * v_num - 2) / 2
+ elif manova_test == "Hotelling-Lawley Trace":
+ df2 = s_val * (v_den - dv - 1) + 4
+ else:
+ df2 = s_val * (v_den - dv + 1) + 4
+ ncp = f2 * n_try * (1 - rho)
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ f_crit = f_dist.ppf(1 - alpha, df1, df2)
+ p_cur = 1 - noncentral_f.cdf(f_crit, df1, df2, ncp)
+ if p_cur >= power:
+ n_total = n_try
+ break
+ if n_total is not None:
+ n_per_group = int(np.ceil(n_total / k))
+ n_total = n_per_group * k
+ explanation = (
+ f"Required total sample size for MANOVA with {k} groups, {dv} dependent variables, "
+ f"effect size f²(V) = {f2:.4f}, α = {alpha}, power = {power}, "
+ f"DV correlation ρ = {rho:.2f}. Test: {manova_test}."
+ )
+ formula_latex = r"n \text{ from non-central } F \text{ with } df_1 = u \cdot v_{\text{num}}, \quad f^2(V) = \frac{V}{s - V}"
+
+ elif atype == "binomial":
+ p0 = params["p0"]
+ p1 = params["p1"]
+ if p0 != p1:
+ from scipy.stats import binom
+
+ n_total = None
+ for n_try in range(3, 10000):
+ if alternative == "two-sided":
+ alpha_lo = binom.ppf(alpha / 2, n_try, p0)
+ alpha_hi = binom.ppf(1 - alpha / 2, n_try, p0)
+ p_pow = binom.cdf(alpha_hi, n_try, p1) - binom.cdf(
+ alpha_lo - 1, n_try, p1
+ )
+ elif p1 > p0:
+ crit = binom.ppf(1 - alpha, n_try, p0)
+ p_pow = 1 - binom.cdf(crit - 1, n_try, p1)
+ else:
+ crit = binom.ppf(alpha, n_try, p0)
+ p_pow = binom.cdf(crit, n_try, p1)
+ if p_pow >= power:
+ n_total = n_try
+ break
+ n_per_group = n_total
+ explanation = (
+ f"Required sample size for binomial exact test ({tails.lower()}) "
+ f"to detect π₁ = {p1:.3f} vs π₀ = {p0:.3f} "
+ f"with α = {alpha}, power = {power}. Uses exact binomial CDF."
+ )
+ _z_sub_phi = r"1-\alpha/2" if alternative == "two-sided" else r"1-\alpha"
+ formula_latex = rf"n = \min\{{n: 1 - \Phi(\Phi^{{-1}}({_z_sub_phi}) - \sqrt{{n}} d) \geq \text{{power}}\}}"
+
+ elif atype == "simulation":
+ sim_test = params["sim_test"]
+ n_sim = params["n_sim"]
+ n_per = params["n_per"]
+ n_total = n_per * 2
+ n_per_group = (n_per, n_per)
+ alpha_sim = alpha
+ st.info(
+ f"Running {n_sim} Monte Carlo simulations with N = {n_per} per group..."
+ )
+ progress_bar = st.progress(0)
+ rejects = 0
+ np.random.seed(42)
+ if sim_test == "Two-proportion z-test":
+ p1_s = params["p1_s"]
+ p2_s = params["p2_s"]
+ for i in range(n_sim):
+ x1 = np.random.binomial(1, p1_s, n_per)
+ x2 = np.random.binomial(1, p2_s, n_per)
+ from scipy.stats import chi2_contingency
+
+ _, p_val, _, _ = chi2_contingency(
+ pd.crosstab(
+ pd.Series(np.concatenate([x1, x2])),
+ pd.Series(["G1"] * n_per + ["G2"] * n_per),
+ )
+ )
+ if p_val < alpha_sim:
+ rejects += 1
+ if (i + 1) % max(1, n_sim // 20) == 0:
+ progress_bar.progress((i + 1) / n_sim)
+ else:
+ mu1 = params["mu1"]
+ mu2 = params["mu2"]
+ sd = params["sd"]
+ dist = params.get("dist", "Normal")
+ for i in range(n_sim):
+ if dist == "Normal":
+ g1 = np.random.normal(mu1, sd, n_per)
+ g2 = np.random.normal(mu2, sd, n_per)
+ elif dist == "Skewed (Exponential)":
+ g1 = np.random.exponential(sd, n_per) + mu1
+ g2 = np.random.exponential(sd, n_per) + mu2
+ else:
+ g1 = np.random.uniform(mu1 - sd * 1.732, mu1 + sd * 1.732, n_per)
+ g2 = np.random.uniform(mu2 - sd * 1.732, mu2 + sd * 1.732, n_per)
+ from scipy.stats import ttest_ind, mannwhitneyu
+
+ if sim_test == "Independent t-test (pooled)":
+ _, p_val = ttest_ind(g1, g2, equal_var=True)
+ elif sim_test == "Welch's t-test":
+ _, p_val = ttest_ind(g1, g2, equal_var=False)
+ else:
+ _, p_val = mannwhitneyu(g1, g2, alternative="two-sided")
+ if p_val < alpha_sim:
+ rejects += 1
+ if (i + 1) % max(1, n_sim // 20) == 0:
+ progress_bar.progress((i + 1) / n_sim)
+ progress_bar.empty()
+ empirical_power = rejects / n_sim
+ st.success(f"Monte Carlo simulation complete — {rejects}/{n_sim} rejections.")
+ params["_empirical_power"] = empirical_power
+ explanation = (
+ f"Monte Carlo simulation ({n_sim} iterations) using {sim_test} "
+ f"with N = {n_per} per group. "
+ f"Estimated power: {empirical_power:.1%} "
+ f"(95% CI: {max(0, empirical_power - 1.96*np.sqrt(empirical_power*(1-empirical_power)/n_sim)):.1%} – "
+ f"{min(1, empirical_power + 1.96*np.sqrt(empirical_power*(1-empirical_power)/n_sim)):.1%})."
+ )
+ formula_latex = r"\text{Power} = \frac{1}{N_{\text{sim}}} \sum_{i=1}^{N_{\text{sim}}} I(p_i < \alpha)"
+
+ # --- Compromise fallback for test types not explicitly handled ---
+ if (
+ analysis_mode == "Compromise"
+ and computed_value is None
+ and n_total is not None
+ and n_total > 0
+ ):
+ from scipy.optimize import brentq
+
+ es_val = params.get("effect_size", None)
+ if es_val is not None and es_val > 0:
+ from statsmodels.stats.power import TTestPower, NormalIndPower
+
+ if atype in ("one_mean", "paired", "correlation"):
+ solver = TTestPower()
+ elif atype in ("one_prop", "mannwhitney", "wilcoxon_sr"):
+ solver = NormalIndPower()
+ else:
+ solver = None
+ if solver is not None:
+ cost_ratio = params.get("cost_ratio", 1.0)
+ try:
+ comp = _solve_compromise(
+ solver, es_val, n_total, cost_ratio, alternative=alternative
+ )
+ if comp is not None:
+ computed_value = comp
+ explanation = (
+ f"Compromise power analysis for {atype} with N = {n_total}, "
+ f"effect size = {es_val:.4f}, cost ratio β/α = {cost_ratio:.2f}. "
+ f"Adjusted α = {computed_value['alpha']:.4f}, achieved power = {computed_value['power']:.1%}."
+ )
+ except Exception:
+ pass
+
+ # --- Apply Attrition Adjustment ---
+ n_total_raw = n_total
+ if analysis_mode == "A Priori" and dropout_rate > 0 and n_total is not None:
+ n_total = int(np.ceil(n_total / (1 - dropout_rate)))
+ if isinstance(n_per_group, tuple):
+ n1_adj = int(np.ceil(n_per_group[0] / (1 - dropout_rate)))
+ n2_adj = int(np.ceil(n_per_group[1] / (1 - dropout_rate)))
+ n_per_group = (n1_adj, n2_adj)
+ elif n_per_group is not None:
+ n_per_group = int(np.ceil(n_per_group / (1 - dropout_rate)))
+
+ # --- Display Results ---
+
+ if n_total is None and computed_value is None:
+ st.error(
+ "Effect size is too small — consider a larger effect or different design."
+ )
+ return
+
+ mode_titles = {
+ "A Priori": "Sample Size Results",
+ "Post Hoc": "Post-Hoc Power Analysis Results",
+ "Sensitivity": "Sensitivity Analysis Results",
+ "Criterion": "Criterion Analysis Results",
+ "Compromise": "Compromise Power Analysis Results",
+ }
+ st.subheader(mode_titles.get(analysis_mode, "Power Analysis Results"))
+
+ if formula_latex:
+ st.latex(formula_latex)
+
+ if analysis_mode == "A Priori":
+ col1, col2, col3 = st.columns(3)
+ with col1:
+ st.metric("Total Sample Size (N)", n_total)
+ with col2:
+ if isinstance(n_per_group, tuple):
+ st.metric("Group 1 (n₁)", n_per_group[0])
+ st.metric("Group 2 (n₂)", n_per_group[1])
+ else:
+ st.metric("Per Group (n)", n_per_group)
+ with col3:
+ if atype == "simulation":
+ emp_pwr = params.get("_empirical_power", 0)
+ st.metric("Empirical Power", f"{emp_pwr:.1%}")
+ else:
+ st.metric("Power", f"{power:.0%}")
+ st.metric("Alpha (α)", f"{alpha:.3f}")
+
+ adjustments = []
+ if dropout_rate > 0 and n_total_raw is not None and n_total_raw != n_total:
+ adjustments.append(
+ f"Raw N = {n_total_raw} (before {dropout_rate:.0%} dropout adjustment)"
+ )
+ if num_tests > 1:
+ adjustments.append(
+ f"{mc_method} α = {alpha:.4f} (original: {alpha_raw:.4f}) for {num_tests} comparisons"
+ )
+ if adjustments:
+ st.caption(" | ".join(adjustments))
+
+ elif analysis_mode == "Post Hoc":
+ col1, col2, col3 = st.columns(3)
+ with col1:
+ st.metric("Total Sample Size (N)", n_total)
+ with col2:
+ if isinstance(n_per_group, tuple):
+ st.metric("Group 1 (n₁)", n_per_group[0])
+ st.metric("Group 2 (n₂)", n_per_group[1])
+ else:
+ st.metric("Per Group (n)", n_per_group)
+ with col3:
+ pwr_display = f"{computed_value:.1%}" if computed_value is not None else "—"
+ st.metric("Achieved Power", pwr_display)
+ st.metric("Alpha (α)", f"{alpha:.3f}")
+
+ elif analysis_mode == "Sensitivity":
+ col1, col2, col3 = st.columns(3)
+ with col1:
+ st.metric("Total Sample Size (N)", n_total)
+ with col2:
+ st.metric("Power", f"{power:.0%}")
+ st.metric("Alpha (α)", f"{alpha:.3f}")
+ with col3:
+ es_display = f"{computed_value:.4f}" if computed_value is not None else "—"
+ st.metric("Min Detectable Effect", es_display)
+
+ elif analysis_mode == "Criterion":
+ col1, col2, col3 = st.columns(3)
+ with col1:
+ st.metric("Total Sample Size (N)", n_total)
+ with col2:
+ if isinstance(n_per_group, tuple):
+ st.metric("Group 1 (n₁)", n_per_group[0])
+ st.metric("Group 2 (n₂)", n_per_group[1])
+ else:
+ st.metric("Per Group (n)", n_per_group)
+ with col3:
+ st.metric("Power", f"{power:.0%}")
+ a_display = f"{computed_value:.4f}" if computed_value is not None else "—"
+ st.metric("Required α", a_display)
+
+ elif analysis_mode == "Compromise":
+ col1, col2, col3 = st.columns(3)
+ with col1:
+ st.metric("Total Sample Size (N)", n_total)
+ with col2:
+ st.metric(
+ "Adjusted α",
+ (
+ f"{computed_value['alpha']:.4f}"
+ if isinstance(computed_value, dict)
+ else "—"
+ ),
+ )
+ with col3:
+ st.metric(
+ "Achieved Power",
+ (
+ f"{computed_value['power']:.1%}"
+ if isinstance(computed_value, dict)
+ else "—"
+ ),
+ )
+
+ st.info(explanation)
+
+ if isinstance(n_per_group, tuple) and n_total is not None and n_total > 0:
+ labels = [
+ f"Group 1 (n₁ = {n_per_group[0]})",
+ f"Group 2 (n₂ = {n_per_group[1]})",
+ ]
+ values = [n_per_group[0], n_per_group[1]]
+ fig_pie = go.Figure(
+ data=[
+ go.Pie(labels=labels, values=values, hole=0.4, textinfo="label+percent")
+ ]
+ )
+ fig_pie.update_layout(
+ template="plotly_dark",
+ height=280,
+ title="Group-Size Distribution",
+ margin=dict(t=40, b=10, l=10, r=10),
+ )
+ st.plotly_chart(fig_pie, use_container_width=True)
+
+ if atype == "simulation":
+ st.info(
+ "Power curve and sensitivity plots use parametric formulas and are not available for Monte Carlo simulation. Change N above and re-run to test different sample sizes."
+ )
+ return
+
+ # --- Power Curve ---
+ st.subheader("Power Analysis Curve")
+ max_n = max(n_total * 3, 30)
+ n_range = np.arange(5, max_n + 1, 2)
+
+ from statsmodels.stats.power import (
+ TTestIndPower,
+ TTestPower,
+ NormalIndPower,
+ FTestAnovaPower,
+ GofChisquarePower,
+ )
+
+ alpha_curve_values = sorted(set([0.01, 0.05, 0.10, alpha]))
+ alpha_curve_colors = {0.01: "#FF6B6B", 0.05: "#00BFFF", 0.10: "#51CF66"}
+ if alpha not in alpha_curve_values:
+ alpha_curve_colors[alpha] = "#FFD43B"
+
+ fig = go.Figure()
+ for a in alpha_curve_values:
+ saved_alpha = alpha
+ saved_z_alpha = z_alpha
+ alpha = a
+ if alternative == "two-sided":
+ z_alpha = norm.ppf(1 - alpha / 2)
+ else:
+ z_alpha = norm.ppf(1 - alpha)
+
+ power_vals = []
+ for n in n_range:
+ try:
+ if atype == "one_mean":
+ solver = TTestPower()
+ pv = solver.solve_power(
+ effect_size=params["effect_size"],
+ nobs=n,
+ alpha=alpha,
+ alternative=alternative,
+ )
+ elif atype == "two_means":
+ solver = TTestIndPower()
+ pv = solver.solve_power(
+ effect_size=params["effect_size"],
+ nobs1=n,
+ alpha=alpha,
+ ratio=params["ratio"],
+ alternative=alternative,
+ )
+ elif atype == "paired":
+ solver = TTestPower()
+ pv = solver.solve_power(
+ effect_size=params["effect_size"],
+ nobs=n,
+ alpha=alpha,
+ alternative=alternative,
+ )
+ elif atype == "one_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ solver = NormalIndPower()
+ d_eff = proportion_effectsize(
+ params["prop_alt"], params["prop_null"]
+ )
+ pv = solver.solve_power(
+ effect_size=abs(d_eff),
+ nobs1=n,
+ alpha=alpha,
+ alternative=alternative,
+ )
+ elif atype == "two_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ solver = NormalIndPower()
+ d_eff = proportion_effectsize(params["p2"], params["p1"])
+ pv = solver.solve_power(
+ effect_size=abs(d_eff),
+ nobs1=n,
+ alpha=alpha,
+ ratio=params["ratio"],
+ alternative=alternative,
+ )
+ elif atype == "anova":
+ solver = FTestAnovaPower()
+ pv = solver.solve_power(
+ effect_size=params["effect_size"],
+ nobs=n,
+ alpha=alpha,
+ k_groups=params["k"],
+ )
+ elif atype == "correlation":
+ import math
-In the meantime, below is an example of what you can do with just a few lines of code:
+ fisher_z = math.atanh(params["effect_size"])
+ z_alpha_c = (
+ norm.ppf(1 - alpha / 2)
+ if alternative == "two-sided"
+ else norm.ppf(1 - alpha)
+ )
+ z_beta_c = (np.sqrt(n - 3) * fisher_z) - z_alpha_c
+ pv = norm.cdf(z_beta_c)
+ elif atype == "regression":
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ k_r = params["k"]
+ dfd = n - k_r - 1
+ if dfd > 0:
+ ncp = params["effect_size"] * n
+ f_crit = f_dist.ppf(1 - alpha, k_r, dfd)
+ pv = 1 - noncentral_f.cdf(f_crit, k_r, dfd, ncp)
+ else:
+ pv = 0
+ elif atype == "logistic":
+ k = params["k"]
+ ev_rate = params["event_rate"]
+ or_val = params["or"]
+ p1_log = (or_val * ev_rate) / (1 - ev_rate + or_val * ev_rate)
+ d_log = abs(p1_log - ev_rate)
+ p_bar = (ev_rate + p1_log) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ d_eff_log = d_log / se if se > 0 else 0
+ solver = NormalIndPower()
+ pv = solver.solve_power(
+ effect_size=d_eff_log,
+ nobs1=n,
+ alpha=alpha,
+ alternative=alternative,
+ )
+ elif atype == "chisq":
+ solver = GofChisquarePower()
+ pv = solver.solve_power(
+ effect_size=params["effect_size"],
+ nobs=n,
+ alpha=alpha,
+ n_bins=params["df"] + 1,
+ )
+ elif atype == "mannwhitney":
+ p_val = params["effect_size"]
+ ratio = params["ratio"]
+ are = params["are"]
+ d_mw = np.sqrt(3) * (p_val - 0.5)
+ if d_mw > 0:
+ solver = NormalIndPower()
+ pv = solver.solve_power(
+ effect_size=d_mw,
+ nobs1=n * are,
+ alpha=alpha,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ else:
+ pv = 0
+ elif atype == "logrank":
+ hr = params["hr"]
+ ratio = params["ratio"]
+ med_ctrl = params["median_survival"]
+ study_dur = params["study_duration"]
+ log_hr = np.log(hr)
+ lambda_ctrl = np.log(2) / med_ctrl
+ p_event_ctrl = 1 - np.exp(-lambda_ctrl * study_dur)
+ p_event_trt = 1 - np.exp(-lambda_ctrl / hr * study_dur)
+ p_event = (p_event_ctrl + ratio * p_event_trt) / (1 + ratio)
+ num_events = n * p_event
+ if log_hr != 0 and num_events > 0:
+ z_val = (
+ abs(log_hr)
+ * np.sqrt(num_events * ratio / ((ratio + 1) ** 2))
+ - z_alpha
+ )
+ pv = norm.cdf(z_val)
+ else:
+ pv = 0
+ elif atype == "cox":
+ hr = params["hr"]
+ sd_x = params["sd_x"]
+ r2_x = params["r2_x"]
+ ev_rate = params["event_rate"]
+ log_hr = np.log(hr)
+ num_events = n * ev_rate
+ if log_hr != 0 and num_events > 0:
+ z_val = (
+ abs(log_hr) * sd_x * np.sqrt(num_events * (1 - r2_x))
+ - z_alpha
+ )
+ pv = norm.cdf(z_val)
+ else:
+ pv = 0
+ elif atype == "equiv":
+ margin = params["margin"]
+ exp_diff = params["expected_diff"]
+ sd = params["sd"]
+ ratio = params["ratio"]
+ equiv_type = params.get("equiv_param_type", "Mean")
+ if equiv_type == "Proportion":
+ p1_eq = params.get("p1_eq", 0.2)
+ p2_eq = params.get("p2_eq", 0.2)
+ d_prop = abs(p1_eq - p2_eq)
+ d_e = margin - d_prop
+ if d_e > 0:
+ p_bar = (p1_eq + p2_eq) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ es_e = d_e / se if se > 0 else 0
+ solver = NormalIndPower()
+ pv = solver.solve_power(
+ effect_size=es_e,
+ nobs1=n,
+ alpha=alpha,
+ ratio=ratio,
+ alternative="larger",
+ )
+ else:
+ pv = 0
+ else:
+ d_e = margin - abs(exp_diff)
+ if d_e > 0:
+ es_e = d_e / sd
+ solver = NormalIndPower()
+ pv = solver.solve_power(
+ effect_size=es_e,
+ nobs1=n,
+ alpha=alpha,
+ ratio=ratio,
+ alternative="larger",
+ )
+ else:
+ pv = 0
+ elif atype == "rm_anova":
+ f_eff = params["effect_size"]
+ k = params["k"]
+ m = params["m"]
+ rho = params["rho"]
+ epsilon = params["epsilon"]
+ if f_eff > 0:
+ solver = FTestAnovaPower()
+ design_effect = (1 + (m - 1) * rho) / m
+ df_adj = (m - 1) * epsilon
+ n_indep = max(n * df_adj / (design_effect * (k - 1)), k + 1)
+ pv = solver.solve_power(
+ effect_size=f_eff,
+ nobs=n_indep,
+ alpha=alpha,
+ k_groups=k,
+ )
+ else:
+ pv = 0
+ elif atype == "twoway_anova":
+ f_a = params["f_a"]
+ f_b = params["f_b"]
+ f_ab = params["f_ab"]
+ rows = params["rows"]
+ cols = params["cols"]
+ focus = params["focus"]
+ if focus == "Main Effect A":
+ f_use = f_a
+ k_use = rows
+ elif focus == "Main Effect B":
+ f_use = f_b
+ k_use = cols
+ else:
+ f_use = f_ab
+ k_use = rows * cols
+ if f_use > 0:
+ solver = FTestAnovaPower()
+ pv = solver.solve_power(
+ effect_size=f_use,
+ nobs=n,
+ alpha=alpha,
+ k_groups=k_use,
+ )
+ else:
+ pv = 0
+ elif atype == "roc_auc":
+ auc = params["auc"]
+ null_auc = params["null_auc"]
+ ratio = params["ratio"]
+ v_auc = auc * (1 - auc) + (ratio - 1) * (
+ auc / (2 - auc) - auc**2
+ ) / (1 + ratio)
+ delta_auc = auc - null_auc
+ if delta_auc > 0:
+ n_cases_eff = n / (1 + ratio)
+ z_val = delta_auc * np.sqrt(n_cases_eff / v_auc) - z_alpha
+ pv = norm.cdf(z_val)
+ else:
+ pv = 0
+ elif atype == "kappa":
+ kappa = params["kappa"]
+ null_kappa = params["null_kappa"]
+ raters = params["raters"]
+ cats = params["categories"]
+ delta_k = kappa - null_kappa
+ if delta_k > 0:
+ z_val = (
+ delta_k * np.sqrt(n / (null_kappa * (1 - null_kappa)))
+ - z_alpha
+ )
+ pv = norm.cdf(z_val)
+ else:
+ pv = 0
+ elif atype == "cluster_rct":
+ d = params["effect_size"]
+ icc = params["icc"]
+ cluster_m = params["cluster_size"]
+ ratio = params["ratio"]
+ deff = 1 + (cluster_m - 1) * icc
+ if d > 0:
+ n1_equiv = (n / (1 + ratio)) / deff
+ solver = TTestIndPower()
+ if n1_equiv > 1:
+ pv = solver.solve_power(
+ effect_size=d,
+ nobs1=n1_equiv,
+ alpha=alpha,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ else:
+ pv = 0
+ else:
+ pv = 0
+ elif atype == "precision":
+ half_width = params["half_width"]
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and half_width > 0:
+ from scipy.stats import t as t_dist
+
+ n_req = max(
+ (norm.ppf(1 - conf_alpha / 2) * sd / half_width) ** 2, 3
+ )
+ for _ in range(20):
+ t_val = t_dist.ppf(
+ 1 - conf_alpha / 2, df=int(n_req) - 1
+ )
+ n_next = (t_val * sd / half_width) ** 2
+ if abs(n_next - n_req) < 0.5:
+ break
+ n_req = max(n_next, 3)
+ else:
+ n_req = 0
+ else:
+ prop = params["prop"]
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ if prop > 0 and half_width > 0:
+ n_req = (z_hw**2 * prop * (1 - prop)) / (half_width**2)
+ else:
+ n_req = 0
+ if n_req > 0:
+ pv = min(n / n_req, 1.0)
+ else:
+ pv = 0
+ elif atype == "pilot":
+ method = params["method"]
+ if method == "Rule of thumb":
+ n_req = params["n_per_group"]
+ pv = min(n / n_req, 1.0) if n_req > 0 else 0
+ elif method == "Precision-based":
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ if param_type == "Mean":
+ sd = params["sd"]
+ half_width = params["half_width"]
+ if sd > 0 and half_width > 0:
+ from scipy.stats import t as t_dist
+
+ n_req = max(
+ (norm.ppf(1 - conf_alpha / 2) * sd / half_width)
+ ** 2,
+ 3,
+ )
+ for _ in range(20):
+ t_val = t_dist.ppf(
+ 1 - conf_alpha / 2, df=int(n_req) - 1
+ )
+ n_next = (t_val * sd / half_width) ** 2
+ if abs(n_next - n_req) < 0.5:
+ break
+ n_req = max(n_next, 3)
+ else:
+ n_req = 0
+ else:
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ prop = params["prop"]
+ half_width = params["half_width"]
+ n_req = (
+ (z_hw**2 * prop * (1 - prop)) / (half_width**2)
+ if prop > 0 and half_width > 0
+ else 0
+ )
+ pv = min(n / n_req, 1.0) if n_req > 0 else 0
+ else:
+ main_n = params["main_n"]
+ fraction = params["fraction"]
+ n_req = max(int(np.ceil(main_n * fraction)), 5)
+ pv = min(n / n_req, 1.0)
+ elif atype == "wilcoxon_sr":
+ p_val = params["effect_size"]
+ are = params["are"]
+ d_wsr = np.sqrt(3) * (p_val - 0.5) * 2
+ if d_wsr > 0:
+ solver = NormalIndPower()
+ pv = solver.solve_power(
+ effect_size=d_wsr,
+ nobs1=n * are,
+ alpha=alpha,
+ alternative=alternative,
+ )
+ else:
+ pv = 0
+ elif atype == "kruskal":
+ f_eff = params["effect_size"]
+ k = params["k"]
+ are_kw = params.get("are", 0.955)
+ if f_eff > 0:
+ solver = FTestAnovaPower()
+ pv = solver.solve_power(
+ effect_size=f_eff, nobs=n * are_kw, alpha=alpha, k_groups=k
+ )
+ else:
+ pv = 0
+ elif atype == "friedman":
+ w = params["w"]
+ k = params["k"]
+ m = params["m"]
+ are = params["are"]
+ f_fr = np.sqrt(w / (1 - w)) if w < 1 else 1
+ if f_fr > 0:
+ solver = FTestAnovaPower()
+ n_indep = n * are / m
+ pv = solver.solve_power(
+ effect_size=f_fr, nobs=n_indep, alpha=alpha, k_groups=k
+ )
+ else:
+ pv = 0
+ elif atype == "mcnemar":
+ p_b = params["p_b"]
+ p_c = params["p_c"]
+ d_mc = abs(p_b - p_c)
+ p_discordant = p_b + p_c
+ if d_mc > 0 and p_discordant > 0:
+ z_val = d_mc * np.sqrt(n / p_discordant) - z_alpha
+ pv = norm.cdf(z_val)
+ else:
+ pv = 0
+ elif atype == "fisher":
+ p1 = params["p1"]
+ p2 = params["p2"]
+ ratio = params["ratio"]
+ are_f = params.get("are", 0.833)
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ d_eff = proportion_effectsize(p2, p1)
+ if abs(d_eff) > 0:
+ solver = NormalIndPower()
+ pv = solver.solve_power(
+ effect_size=abs(d_eff),
+ nobs1=n / (1 + ratio) * are_f,
+ alpha=alpha,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ else:
+ pv = 0
+ elif atype == "manova":
+ k = params["k"]
+ dv = params["dv"]
+ f2 = params["f2"]
+ rho = params["rho"]
+ manova_test = params.get("manova_test", "Pillai's Trace")
+ if f2 > 0:
+ u = dv
+ v_num = k - 1
+ v_den = n - k - dv
+ if v_den > 0:
+ s_val = min(u, v_num)
+ df1 = u * v_num
+ if manova_test == "Pillai's Trace":
+ df2 = s_val * (v_den - dv + 1) + 4
+ elif manova_test == "Wilks' Lambda":
+ t_val = max(
+ np.sqrt(
+ (u**2 * v_num**2 - 4)
+ / max(u**2 + v_num**2 - 5, 1)
+ ),
+ 1,
+ )
+ df2 = (v_den - (u - v_num + 1) / 2) * t_val - (
+ u * v_num - 2
+ ) / 2
+ elif manova_test == "Hotelling-Lawley Trace":
+ df2 = s_val * (v_den - dv - 1) + 4
+ else:
+ df2 = s_val * (v_den - dv + 1) + 4
+ ncp = f2 * n * (1 - rho)
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ f_crit = f_dist.ppf(1 - alpha, df1, df2)
+ pv = 1 - noncentral_f.cdf(f_crit, df1, df2, ncp)
+ else:
+ pv = 0
+ else:
+ pv = 0
+ elif atype == "binomial":
+ p0 = params["p0"]
+ p1 = params["p1"]
+ if p0 != p1:
+ from scipy.stats import binom
+
+ if alternative == "two-sided":
+ alpha_lo = binom.ppf(alpha / 2, n, p0)
+ alpha_hi = binom.ppf(1 - alpha / 2, n, p0)
+ pv = binom.cdf(alpha_hi, n, p1) - binom.cdf(
+ alpha_lo - 1, n, p1
+ )
+ elif p1 > p0:
+ crit = binom.ppf(1 - alpha, n, p0)
+ pv = 1 - binom.cdf(crit - 1, n, p1)
+ else:
+ crit = binom.ppf(alpha, n, p0)
+ pv = binom.cdf(crit, n, p1)
+ else:
+ pv = 0
+ else:
+ pv = 0
+ except Exception:
+ pv = 0
+ power_vals.append(pv)
+
+ alpha = saved_alpha
+ z_alpha = saved_z_alpha
+
+ is_selected = abs(a - saved_alpha) < 1e-10
+ line_w = 4 if is_selected else 2
+ line_d = "solid" if is_selected else "dash"
+ fig.add_trace(
+ go.Scatter(
+ x=n_range,
+ y=power_vals,
+ mode="lines",
+ name=f"α = {a:.3f}",
+ line=dict(
+ color=alpha_curve_colors.get(a, "#FFD43B"),
+ width=line_w,
+ dash=line_d,
+ ),
+ )
+ )
+
+ target_power_val = params["power"]
+ fig.add_hline(
+ y=target_power_val,
+ line_dash="dash",
+ line_color="orange",
+ annotation_text=f"Target Power ({target_power_val:.0%})",
+ )
+ fig.add_vline(
+ x=n_total,
+ line_dash="dot",
+ line_color="green",
+ annotation_text=f"N = {n_total}",
+ )
+ fig.update_layout(
+ template="plotly_dark",
+ height=450,
+ xaxis_title="Sample Size (N)",
+ yaxis_title="Power (1 − β)",
+ )
+ st.plotly_chart(fig, use_container_width=True)
+
+ # --- Distribution Plot ---
+ st.subheader(":orange[Distribution Plot (Null & Alternative)]")
+ st.caption(
+ "Shows the sampling distributions under H₀ (null) and H₁ (alternative) with critical regions."
+ )
+ x_min = -4.5
+ x_max = 4.5
+ if n_total is not None and n_total > 0:
+ sqrt_n = np.sqrt(n_total)
+ if atype in ("one_mean", "paired", "wilcoxon_sr"):
+ d = (
+ params.get("effect_size", 0.5)
+ if atype != "wilcoxon_sr"
+ else np.sqrt(3) * (params["effect_size"] - 0.5) * 2
+ )
+ ncp = d * sqrt_n
+ elif atype == "two_means":
+ d = params.get("effect_size", 0.5)
+ ncp = d * sqrt_n / np.sqrt(1 + 1 / params.get("ratio", 1))
+ elif atype in ("one_prop", "binomial"):
+ p0 = params.get("prop_null", params.get("p0", 0.5))
+ p1 = params.get("prop_alt", params.get("p1", 0.7))
+ d_eff = abs(p1 - p0) / np.sqrt(p0 * (1 - p0)) if p0 > 0 and p0 < 1 else 0
+ ncp = d_eff * sqrt_n
+ elif atype in ("two_prop", "fisher"):
+ p1 = params["p1"]
+ ratio = params.get("ratio", 1)
+ d_eff = abs(p1 - params.get("p2", 0.5)) / np.sqrt(p1 * (1 - p1))
+ ncp = d_eff * sqrt_n / np.sqrt(1 + ratio)
+ elif atype in ("mcnemar",):
+ p_b = params.get("p_b", 0.2)
+ p_c = params.get("p_c", 0.4)
+ d_mc = abs(p_b - p_c)
+ p_disc = p_b + p_c
+ ncp = d_mc * sqrt_n / np.sqrt(p_disc) if p_disc > 0 else 0
+ elif atype in ("correlation",):
+ r = params.get("effect_size", 0.3)
+ ncp = np.arctanh(r) * sqrt_n
+ elif atype in ("mannwhitney",):
+ p = params.get("effect_size", 0.65)
+ d_mw = np.sqrt(3) * (p - 0.5)
+ ncp = d_mw * sqrt_n * np.sqrt(params.get("are", 1))
+ elif atype in ("logistic",):
+ or_val = params.get("or", 2)
+ ev_rate = params.get("event_rate", 0.3)
+ p1_log = (or_val * ev_rate) / (1 - ev_rate + or_val * ev_rate)
+ d_log = abs(p1_log - ev_rate)
+ p_bar = (ev_rate + p1_log) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ d_eff = d_log / se if se > 0 else 0
+ ncp = d_eff * sqrt_n
+ elif atype in ("logrank",):
+ hr = params.get("hr", 2)
+ ncp = abs(np.log(hr)) * sqrt_n / 2
+ elif atype in ("cox",):
+ hr = params.get("hr", 2)
+ sd_x = params.get("sd_x", 1)
+ r2_x = params.get("r2_x", 0)
+ ev_rate = params.get("event_rate", 0.5)
+ ncp = abs(np.log(hr)) * sd_x * np.sqrt(ev_rate * (1 - r2_x)) * sqrt_n
+ else:
+ ncp = 2.5
+
+ if ncp > 0:
+ x_range = np.linspace(max(-4.5, -ncp - 4), max(4.5, ncp + 4), 400)
+ null_pdf = norm.pdf(x_range, 0, 1)
+ alt_pdf = norm.pdf(x_range, ncp, 1)
+ z_crit = z_alpha
+ fig_dist = go.Figure()
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_range,
+ y=null_pdf,
+ mode="lines",
+ name="H₀ (null)",
+ line=dict(color="#00BFFF", width=2),
+ )
+ )
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_range,
+ y=alt_pdf,
+ mode="lines",
+ name="H₁ (alternative)",
+ line=dict(color="#FF6B6B", width=2),
+ )
+ )
+ if alternative == "two-sided":
+ x_shade_left = x_range[x_range <= -z_crit]
+ x_shade_right = x_range[x_range >= z_crit]
+ if len(x_shade_left) > 0:
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_shade_left,
+ y=null_pdf[: len(x_shade_left)],
+ mode="lines",
+ fill="tozeroy",
+ name=f"α/2 ({alpha/2:.4f})",
+ line=dict(color="rgba(255,0,0,0.3)"),
+ fillcolor="rgba(255,0,0,0.2)",
+ )
+ )
+ if len(x_shade_right) > 0:
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_shade_right,
+ y=null_pdf[len(null_pdf) - len(x_shade_right) :],
+ mode="lines",
+ fill="tozeroy",
+ name=f"α/2 ({alpha/2:.4f})",
+ line=dict(color="rgba(255,0,0,0.3)"),
+ fillcolor="rgba(255,0,0,0.2)",
+ )
+ )
+ x_beta = x_range[x_range <= z_crit]
+ if len(x_beta) > 0:
+ beta_pdf = alt_pdf[: len(x_beta)]
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_beta,
+ y=beta_pdf,
+ mode="lines",
+ fill="tozeroy",
+ name=f"β ({1-power:.2f})",
+ line=dict(color="rgba(255,165,0,0.3)"),
+ fillcolor="rgba(255,165,0,0.2)",
+ )
+ )
+ else:
+ ncp_pos = ncp > 0
+ if ncp_pos:
+ x_shade = x_range[x_range >= z_crit]
+ x_beta = x_range[x_range <= z_crit]
+ else:
+ x_shade = x_range[x_range <= -z_crit]
+ x_beta = x_range[x_range >= -z_crit]
+ if len(x_shade) > 0:
+ null_shade = (
+ null_pdf[-len(x_shade) :]
+ if ncp_pos
+ else null_pdf[: len(x_shade)]
+ )
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_shade,
+ y=null_shade,
+ mode="lines",
+ fill="tozeroy",
+ name=f"α ({alpha:.4f})",
+ line=dict(color="rgba(255,0,0,0.3)"),
+ fillcolor="rgba(255,0,0,0.2)",
+ )
+ )
+ if len(x_beta) > 0:
+ alt_shade = (
+ alt_pdf[-len(x_beta) :]
+ if not ncp_pos
+ else alt_pdf[: len(x_beta)]
+ )
+ fig_dist.add_trace(
+ go.Scatter(
+ x=x_beta,
+ y=alt_shade,
+ mode="lines",
+ fill="tozeroy",
+ name=f"β ({1-power:.2f})",
+ line=dict(color="rgba(255,165,0,0.3)"),
+ fillcolor="rgba(255,165,0,0.2)",
+ )
+ )
+
+ fig_dist.add_vline(
+ x=0,
+ line_dash="dot",
+ line_color="#00BFFF",
+ annotation_text="H₀ center",
+ )
+ fig_dist.add_vline(
+ x=ncp,
+ line_dash="dot",
+ line_color="#FF6B6B",
+ annotation_text="H₁ center",
+ )
+ fig_dist.update_layout(
+ template="plotly_dark",
+ height=350,
+ xaxis_title="Test Statistic (z)",
+ yaxis_title="Density",
+ legend=dict(orientation="h", y=1.1),
+ )
+ st.plotly_chart(fig_dist, use_container_width=True)
+ else:
+ st.info(
+ "Distribution plot not available for this test type with current parameters."
+ )
+
+ # --- Sensitivity Table ---
+ st.subheader("Sensitivity Analysis")
+ st.caption("How required sample size changes with different effect sizes.")
+
+ effect_multipliers = [0.5, 0.67, 0.8, 1.0, 1.25, 1.5, 2.0]
+ sens_data = []
+ for mult in effect_multipliers:
+ try:
+ if atype in ("one_mean", "paired", "correlation"):
+ adj_es = params["effect_size"] * mult
+ if atype == "correlation":
+ import math
+
+ adj_es = max(min(adj_es, 0.99), 0.01)
+ fisher_z_s = math.atanh(adj_es)
+ n_s = int(np.ceil(3 + ((z_alpha + z_beta) / fisher_z_s) ** 2))
+ elif atype in ("one_mean", "paired"):
+ solver = TTestPower()
+ n_s = int(
+ np.ceil(
+ solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ )
+ )
+ elif atype == "two_means":
+ adj_es = params["effect_size"] * mult
+ solver = TTestIndPower()
+ n1_s = solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ ratio=params["ratio"],
+ alternative=alternative,
+ )
+ n1_s = int(np.ceil(n1_s))
+ n_s = n1_s + int(np.ceil(n1_s * params["ratio"]))
+ elif atype == "one_prop":
+ adj_es = (
+ params["prop_null"]
+ + (params["prop_alt"] - params["prop_null"]) * mult
+ )
+ adj_es = max(min(adj_es, 0.99), 0.01)
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ d_eff_s = proportion_effectsize(adj_es, params["prop_null"])
+ solver = NormalIndPower()
+ n_s = int(
+ np.ceil(
+ solver.solve_power(
+ effect_size=abs(d_eff_s),
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ )
+ )
+ elif atype == "two_prop":
+ adj_p2 = params["p1"] + (params["p2"] - params["p1"]) * mult
+ adj_p2 = max(min(adj_p2, 0.99), 0.01)
+ adj_es = abs(adj_p2 - params["p1"])
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ d_eff_s = proportion_effectsize(adj_p2, params["p1"])
+ solver = NormalIndPower()
+ n1_s = solver.solve_power(
+ effect_size=abs(d_eff_s),
+ alpha=alpha,
+ power=power,
+ ratio=params["ratio"],
+ alternative=alternative,
+ )
+ n1_s = int(np.ceil(n1_s))
+ n_s = n1_s + int(np.ceil(n1_s * params["ratio"]))
+ elif atype == "logistic":
+ adj_es = max(params["or"] ** mult, 1.01)
+ ev_rate = params["event_rate"]
+ k_log = params["k"]
+ p1_log = (adj_es * ev_rate) / (1 - ev_rate + adj_es * ev_rate)
+ d_log = abs(p1_log - ev_rate)
+ p_bar = (ev_rate + p1_log) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ d_eff_log = d_log / se if se > 0 else 0
+ solver = NormalIndPower()
+ n_base = int(
+ np.ceil(
+ solver.solve_power(
+ effect_size=d_eff_log,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ )
+ )
+ n_s = max(n_base, 10 * k_log)
+ elif atype == "anova":
+ adj_es = params["effect_size"] * mult
+ solver = FTestAnovaPower()
+ n_s = (
+ int(
+ np.ceil(
+ solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ k_groups=params["k"],
+ )
+ )
+ )
+ * params["k"]
+ )
+ elif atype == "regression":
+ adj_es = params["effect_size"] * mult
+ k_r = params["k"]
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ n_s = None
+ for nc in range(k_r + 2, 10000):
+ dfd = nc - k_r - 1
+ ncp = adj_es * nc
+ f_crit = f_dist.ppf(1 - alpha, k_r, dfd)
+ p_cur = 1 - noncentral_f.cdf(f_crit, k_r, dfd, ncp)
+ if p_cur >= power:
+ n_s = nc
+ break
+ elif atype == "chisq":
+ adj_es = params["effect_size"] * mult
+ solver = GofChisquarePower()
+ n_s = int(
+ np.ceil(
+ solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ n_bins=params["df"] + 1,
+ )
+ )
+ )
+ elif atype == "mannwhitney":
+ adj_p = params["effect_size"]
+ adj_p = max(min(adj_p * mult, 0.99), 0.01)
+ adj_es = np.sqrt(3) * (adj_p - 0.5)
+ are = params["are"]
+ ratio = params["ratio"]
+ solver = NormalIndPower()
+ n1_s = solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1_s = int(np.ceil(n1_s / are))
+ n_s = n1_s + int(np.ceil(n1_s * ratio))
+ elif atype == "logrank":
+ adj_es = max(params["hr"] ** mult, 1.001)
+ ratio = params["ratio"]
+ med_ctrl = params["median_survival"]
+ study_dur = params["study_duration"]
+ log_hr = np.log(adj_es)
+ num_events = (
+ ((z_alpha + z_beta) ** 2)
+ * ((ratio + 1) ** 2)
+ / (ratio * (log_hr**2))
+ )
+ lambda_ctrl = np.log(2) / med_ctrl
+ p_event_ctrl = 1 - np.exp(-lambda_ctrl * study_dur)
+ p_event_trt = 1 - np.exp(-lambda_ctrl / adj_es * study_dur)
+ p_event = (p_event_ctrl + ratio * p_event_trt) / (1 + ratio)
+ if p_event > 0:
+ n_s = int(np.ceil(num_events / p_event))
+ else:
+ n_s = None
+ elif atype == "cox":
+ adj_es = max(params["hr"] ** mult, 1.001)
+ k = params["k"]
+ sd_x = params["sd_x"]
+ r2_x = params["r2_x"]
+ ev_rate = params["event_rate"]
+ log_hr = np.log(adj_es)
+ var_denom = (sd_x**2) * (log_hr**2) * (1 - r2_x)
+ if var_denom > 0:
+ num_events = ((z_alpha + z_beta) ** 2) / var_denom
+ num_events = max(num_events, 10 * k)
+ n_s = int(np.ceil(num_events / ev_rate))
+ else:
+ n_s = None
+ elif atype == "equiv":
+ adj_margin = params["margin"] * mult
+ exp_diff = params["expected_diff"]
+ sd = params["sd"]
+ ratio = params["ratio"]
+ equiv_type = params.get("equiv_param_type", "Mean")
+ if equiv_type == "Proportion":
+ p1_eq = params.get("p1_eq", 0.2)
+ p2_eq = params.get("p2_eq", 0.2)
+ d_prop = abs(p1_eq - p2_eq)
+ d_e = adj_margin - d_prop
+ if d_e > 0:
+ p_bar = (p1_eq + p2_eq) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ es_e = d_e / se if se > 0 else 0
+ adj_es = es_e
+ solver = NormalIndPower()
+ n1_s = solver.solve_power(
+ effect_size=es_e,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative="larger",
+ )
+ n1_s = int(np.ceil(n1_s))
+ n_s = n1_s + int(np.ceil(n1_s * ratio))
+ else:
+ n_s = None
+ else:
+ d_e = adj_margin - abs(exp_diff)
+ adj_es = d_e / sd if d_e > 0 else 0
+ if d_e > 0:
+ es_e = d_e / sd
+ solver = NormalIndPower()
+ n1_s = solver.solve_power(
+ effect_size=es_e,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative="larger",
+ )
+ n1_s = int(np.ceil(n1_s))
+ n_s = n1_s + int(np.ceil(n1_s * ratio))
+ else:
+ n_s = None
+ elif atype == "rm_anova":
+ adj_es = params["effect_size"] * mult
+ k = params["k"]
+ m = params["m"]
+ rho = params["rho"]
+ epsilon = params["epsilon"]
+ solver = FTestAnovaPower()
+ n_per_g = solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ k_groups=k,
+ )
+ design_effect = (1 + (m - 1) * rho) / m
+ df_adj = (m - 1) * epsilon
+ n_per_g_adj = int(
+ np.ceil(n_per_g * design_effect * k / (k * df_adj / (k - 1)))
+ )
+ n_per_g_adj = max(n_per_g_adj, int(np.ceil(n_per_g)))
+ n_s = n_per_g_adj * k
+ elif atype == "twoway_anova":
+ adj_f_a = params["f_a"] * mult
+ adj_f_b = params["f_b"] * mult
+ adj_f_ab = params["f_ab"] * mult
+ rows = params["rows"]
+ cols = params["cols"]
+ focus = params["focus"]
+ if focus == "Main Effect A":
+ adj_es = adj_f_a
+ k_use = rows
+ elif focus == "Main Effect B":
+ adj_es = adj_f_b
+ k_use = cols
+ else:
+ adj_es = adj_f_ab
+ k_use = rows * cols
+ solver = FTestAnovaPower()
+ n_per_cell = solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ k_groups=k_use,
+ )
+ n_per_cell = int(np.ceil(n_per_cell))
+ n_s = n_per_cell * rows * cols
+ elif atype == "roc_auc":
+ auc = params["auc"]
+ adj_auc = 0.5 + (auc - 0.5) * mult
+ adj_auc = max(min(adj_auc, 0.99), 0.51)
+ null_auc = params["null_auc"]
+ ratio = params["ratio"]
+ v_auc = adj_auc * (1 - adj_auc) + (ratio - 1) * (
+ adj_auc / (2 - adj_auc) - adj_auc**2
+ ) / (1 + ratio)
+ delta_auc = adj_auc - null_auc
+ adj_es = delta_auc
+ if delta_auc > 0:
+ n_cases = ((z_alpha + z_beta) ** 2 * v_auc) / (delta_auc**2)
+ n_cases = int(np.ceil(n_cases))
+ n_controls = int(np.ceil(n_cases * ratio))
+ n_s = n_cases + n_controls
+ else:
+ n_s = None
+ elif atype == "kappa":
+ adj_kappa = params["kappa"]
+ adj_kappa = 0.5 + (adj_kappa - 0.5) * mult
+ adj_kappa = max(min(adj_kappa, 0.99), 0.01)
+ null_kappa = params["null_kappa"]
+ adj_kappa = max(adj_kappa, null_kappa + 0.01)
+ delta_k = adj_kappa - null_kappa
+ raters = params["raters"]
+ cats = params["categories"]
+ adj_es = delta_k
+ if delta_k > 0:
+ n_s = int(
+ np.ceil(
+ ((z_alpha + z_beta) ** 2 * null_kappa * (1 - null_kappa))
+ / (delta_k**2)
+ )
+ )
+ n_s = max(n_s, raters * cats * 5)
+ else:
+ n_s = None
+ elif atype == "cluster_rct":
+ adj_es = params["effect_size"] * mult
+ icc = params["icc"]
+ cluster_m = params["cluster_size"]
+ ratio = params["ratio"]
+ deff = 1 + (cluster_m - 1) * icc
+ solver = TTestIndPower()
+ n1_ind = solver.solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1_ind = int(np.ceil(n1_ind))
+ n1_clust = int(np.ceil(n1_ind * deff))
+ n1_clust = int(np.ceil(n1_clust / cluster_m)) * cluster_m
+ n2_clust = int(np.ceil(n1_clust * ratio))
+ n_s = n1_clust + n2_clust
+ elif atype == "precision":
+ adj_hw = params["half_width"] * mult
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ adj_es = adj_hw
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and adj_hw > 0:
+ from scipy.stats import t as t_dist
+
+ n_s = int(
+ np.ceil((norm.ppf(1 - conf_alpha / 2) * sd / adj_hw) ** 2)
+ )
+ n_s = max(n_s, 3)
+ for _ in range(20):
+ t_val = t_dist.ppf(1 - conf_alpha / 2, df=n_s - 1)
+ n_next = int(np.ceil((t_val * sd / adj_hw) ** 2))
+ n_next = max(n_next, 3)
+ if n_next == n_s:
+ break
+ n_s = n_next
+ else:
+ n_s = None
+ else:
+ prop = params["prop"]
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ if prop > 0 and adj_hw > 0:
+ n_s = int(np.ceil((z_hw**2 * prop * (1 - prop)) / (adj_hw**2)))
+ else:
+ n_s = None
+ if n_s is not None:
+ n_s = max(n_s, 3)
+ elif atype == "pilot":
+ method = params["method"]
+ adj_es = mult
+ if method == "Rule of thumb":
+ n_s = int(np.ceil(params["n_per_group"] / max(mult, 0.1)))
+ elif method == "Precision-based":
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ adj_hw = params["half_width"] * mult
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and adj_hw > 0:
+ from scipy.stats import t as t_dist
+
+ n_s = int(
+ np.ceil(
+ (norm.ppf(1 - conf_alpha / 2) * sd / adj_hw) ** 2
+ )
+ )
+ n_s = max(n_s, 3)
+ for _ in range(20):
+ t_val = t_dist.ppf(1 - conf_alpha / 2, df=n_s - 1)
+ n_next = int(np.ceil((t_val * sd / adj_hw) ** 2))
+ n_next = max(n_next, 3)
+ if n_next == n_s:
+ break
+ n_s = n_next
+ else:
+ n_s = None
+ else:
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ prop = params["prop"]
+ if prop > 0 and adj_hw > 0:
+ n_s = int(
+ np.ceil((z_hw**2 * prop * (1 - prop)) / (adj_hw**2))
+ )
+ else:
+ n_s = None
+ if n_s is not None:
+ n_s = max(n_s, 5)
+ else:
+ main_n = params["main_n"]
+ fraction = params["fraction"]
+ n_s = max(int(np.ceil(main_n * fraction / max(mult, 0.1))), 5)
+ elif atype == "wilcoxon_sr":
+ p_val = params["effect_size"]
+ are = params["are"]
+ d_wsr = np.sqrt(3) * (p_val - 0.5) * 2
+ adj_es = d_wsr * max(mult, 0.1)
+ if d_wsr > 0:
+ solver = NormalIndPower()
+ n_s = int(
+ np.ceil(
+ solver.solve_power(
+ effect_size=adj_es / are,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ / are
+ )
+ )
+ else:
+ n_s = None
+ elif atype == "kruskal":
+ adj_es = params["effect_size"] * mult
+ k = params["k"]
+ are_kw = params.get("are", 0.955)
+ if adj_es > 0:
+ solver = FTestAnovaPower()
+ n_per_g = solver.solve_power(
+ effect_size=adj_es, alpha=alpha, power=power, k_groups=k
+ )
+ n_per_g = int(np.ceil(n_per_g / are_kw))
+ n_s = n_per_g * k
+ else:
+ n_s = None
+ elif atype == "friedman":
+ w = params["w"]
+ k = params["k"]
+ m = params["m"]
+ are = params["are"]
+ adj_w = max(min(w * mult, 0.99), 0.01)
+ adj_es = np.sqrt(adj_w / (1 - adj_w))
+ if adj_es > 0:
+ solver = FTestAnovaPower()
+ n_per_g = solver.solve_power(
+ effect_size=adj_es, alpha=alpha, power=power, k_groups=k
+ )
+ n_per_g = int(np.ceil(n_per_g * m / are))
+ n_s = n_per_g * k
+ else:
+ n_s = None
+ elif atype == "mcnemar":
+ p_b = params["p_b"]
+ p_c = params["p_c"]
+ p_b_adj = max(min(p_b * mult, 0.99), 0.01)
+ p_c_adj = p_c * mult
+ p_c_adj = max(min(p_c_adj, 0.99), 0.01)
+ d_mc = abs(p_b_adj - p_c_adj)
+ p_discordant = p_b_adj + p_c_adj
+ adj_es = d_mc
+ if d_mc > 0 and p_discordant > 0:
+ n_s = int(
+ np.ceil((z_alpha + z_beta) ** 2 * p_discordant / (d_mc**2))
+ )
+ n_s = max(n_s, 5)
+ else:
+ n_s = None
+ elif atype == "fisher":
+ p1 = params["p1"]
+ p2 = params["p2"]
+ ratio = params["ratio"]
+ are_f = params.get("are", 0.833)
+ adj_p2 = p1 + (p2 - p1) * mult
+ adj_p2 = max(min(adj_p2, 0.99), 0.01)
+ adj_es = abs(adj_p2 - p1)
+ from statsmodels.stats.proportion import proportion_effectsize
+
+ d_eff = proportion_effectsize(adj_p2, p1)
+ if abs(d_eff) > 0:
+ solver = NormalIndPower()
+ n1 = solver.solve_power(
+ effect_size=abs(d_eff),
+ alpha=alpha,
+ power=power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1 = int(np.ceil(n1 / are_f))
+ n_s = n1 + int(np.ceil(n1 * ratio))
+ else:
+ n_s = None
+ elif atype == "manova":
+ k = params["k"]
+ dv = params["dv"]
+ f2 = params["f2"]
+ rho = params["rho"]
+ adj_f2 = f2 * mult
+ adj_es = adj_f2
+ manova_test = params.get("manova_test", "Pillai's Trace")
+ if adj_f2 > 0:
+ u = dv
+ v_num = k - 1
+ n_s = None
+ for n_try in range(k * dv + 2, 5000):
+ v_den = n_try - k - dv
+ if v_den <= 0:
+ continue
+ s_val = min(u, v_num)
+ df1 = u * v_num
+ if manova_test == "Pillai's Trace":
+ df2 = s_val * (v_den - dv + 1) + 4
+ elif manova_test == "Wilks' Lambda":
+ t_val = max(
+ np.sqrt(
+ (u**2 * v_num**2 - 4) / max(u**2 + v_num**2 - 5, 1)
+ ),
+ 1,
+ )
+ df2 = (v_den - (u - v_num + 1) / 2) * t_val - (
+ u * v_num - 2
+ ) / 2
+ elif manova_test == "Hotelling-Lawley Trace":
+ df2 = s_val * (v_den - dv - 1) + 4
+ else:
+ df2 = s_val * (v_den - dv + 1) + 4
+ ncp = adj_f2 * n_try * (1 - rho)
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ f_crit = f_dist.ppf(1 - alpha, df1, df2)
+ p_cur = 1 - noncentral_f.cdf(f_crit, df1, df2, ncp)
+ if p_cur >= power:
+ n_s = n_try
+ break
+ else:
+ n_s = None
+ elif atype == "binomial":
+ p0 = params["p0"]
+ p1 = params["p1"]
+ p1_adj = p0 + (p1 - p0) * mult
+ p1_adj = max(min(p1_adj, 0.99), 0.01)
+ adj_es = abs(p1_adj - p0)
+ if p0 != p1_adj:
+ from scipy.stats import binom
+
+ n_s = None
+ for n_try in range(3, 10000):
+ if alternative == "two-sided":
+ alpha_lo = binom.ppf(alpha / 2, n_try, p0)
+ alpha_hi = binom.ppf(1 - alpha / 2, n_try, p0)
+ p_pow = binom.cdf(alpha_hi, n_try, p1_adj) - binom.cdf(
+ alpha_lo - 1, n_try, p1_adj
+ )
+ elif p1_adj > p0:
+ crit = binom.ppf(1 - alpha, n_try, p0)
+ p_pow = 1 - binom.cdf(crit - 1, n_try, p1_adj)
+ else:
+ crit = binom.ppf(alpha, n_try, p0)
+ p_pow = binom.cdf(crit, n_try, p1_adj)
+ if p_pow >= power:
+ n_s = n_try
+ break
+ else:
+ n_s = None
+ else:
+ n_s = None
+ except Exception:
+ n_s = None
+ if n_s is not None and n_s > 0 and n_s < 100000:
+ sens_data.append(
+ {
+ "Effect Size Multiplier": f"{mult:.1f}×",
+ "Adjusted Effect": (
+ f"{adj_es:.3f}" if mult != 1.0 else f"{adj_es:.3f} (baseline)"
+ ),
+ "Required N": n_s,
+ }
+ )
+
+ if sens_data:
+ st.dataframe(pd.DataFrame(sens_data), use_container_width=True, hide_index=True)
+
+ # --- Interpretation Guide ---
+ with st.expander("📖 How to Interpret These Results"):
+ if n_per_group is not None:
+ if isinstance(n_per_group, tuple):
+ n1, n2 = n_per_group
+ st.markdown(f"""
+ - You need **at least {n1} participants in Group 1** and **{n2} in Group 2**.
+ - **Total**: {n_total} participants.
+ - This assumes α = {alpha}, power = {power:.0%}, and your estimated effect size.
+ - Use the **Power Curve** above to see how N affects your study's power.
+ - The **Sensitivity Table** shows how N changes with different effect sizes.
+ """)
+ else:
+ st.markdown(f"""
+ - You need **at least {n_per_group} participants per group** (total N = {n_total}).
+ - This assumes α = {alpha}, power = {power:.0%}, and your estimated effect size.
+ - Use the **Power Curve** above to see how N affects your study's power.
+ - The **Sensitivity Table** shows how N changes with different effect sizes.
+ """)
+ st.warning(
+ "**Disclaimer**: Sample size estimation is based on statistical assumptions. "
+ "Always consult a biostatistician and consider practical constraints (budget, "
+ "dropout rates, feasibility) when finalizing your study size."
+ )
+
+ # --- Budget & Feasibility ---
+ if cost_per > 0 and n_total is not None:
+ st.subheader("💰 Budget & Feasibility")
+ total_cost = n_total * cost_per
+ col_b1, col_b2, col_b3 = st.columns(3)
+ with col_b1:
+ st.metric("Total Study Cost", f"${total_cost:,.0f}")
+ with col_b2:
+ st.metric("Cost per Participant", f"${cost_per:,.0f}")
+ with col_b3:
+ if recruitment_rate > 0:
+ months = int(np.ceil(n_total / recruitment_rate))
+ st.metric("Est. Recruitment Duration", f"{months} months")
+ else:
+ st.metric("Recruitment Rate", "Not specified")
+
+ if dropout_rate > 0 and n_total_raw is not None:
+ st.caption(
+ f"Base N = {n_total_raw}, adjusted for {dropout_rate:.0%} dropout → {n_total}. "
+ f"Extra cost due to dropout: ${(n_total - n_total_raw) * cost_per:,.0f}."
+ )
+
+ # --- What-If Scenario Explorer ---
+ st.subheader("🔍 What-If Scenario Explorer")
+ st.caption(
+ "Explore how sample size changes across different combinations of power and effect size."
+ )
+
+ whatif_powers = [0.70, 0.80, 0.90]
+ whatif_mults = [0.5, 0.75, 1.0, 1.25, 1.5]
+
+ # Compute a heatmap: rows = power, cols = effect multiplier
+ heatmap_data = []
+ for w_power in whatif_powers:
+ w_z_beta = norm.ppf(w_power)
+ row = []
+ for w_mult in whatif_mults:
+ try:
+ if atype == "one_mean":
+ from statsmodels.stats.power import TTestPower
+
+ d = params["effect_size"] * w_mult
+ w_n = int(
+ np.ceil(
+ TTestPower().solve_power(
+ effect_size=d,
+ alpha=alpha,
+ power=w_power,
+ alternative=alternative,
+ )
+ )
+ )
+ elif atype == "two_means":
+ from statsmodels.stats.power import TTestIndPower
+
+ d = params["effect_size"] * w_mult
+ n1 = TTestIndPower().solve_power(
+ effect_size=d,
+ alpha=alpha,
+ power=w_power,
+ ratio=params["ratio"],
+ alternative=alternative,
+ )
+ n1 = int(np.ceil(n1))
+ w_n = n1 + int(np.ceil(n1 * params["ratio"]))
+ elif atype == "paired":
+ from statsmodels.stats.power import TTestPower
+
+ d = params["effect_size"] * w_mult
+ w_n = int(
+ np.ceil(
+ TTestPower().solve_power(
+ effect_size=d,
+ alpha=alpha,
+ power=w_power,
+ alternative=alternative,
+ )
+ )
+ )
+ elif atype == "one_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ adj_p1 = (
+ params["prop_null"]
+ + (params["prop_alt"] - params["prop_null"]) * w_mult
+ )
+ adj_p1 = max(min(adj_p1, 0.99), 0.01)
+ d_eff = proportion_effectsize(adj_p1, params["prop_null"])
+ w_n = int(
+ np.ceil(
+ NormalIndPower().solve_power(
+ effect_size=abs(d_eff),
+ alpha=alpha,
+ power=w_power,
+ alternative=alternative,
+ )
+ )
+ )
+ elif atype == "two_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ adj_p2 = params["p1"] + (params["p2"] - params["p1"]) * w_mult
+ adj_p2 = max(min(adj_p2, 0.99), 0.01)
+ d_eff = proportion_effectsize(adj_p2, params["p1"])
+ n1 = NormalIndPower().solve_power(
+ effect_size=abs(d_eff),
+ alpha=alpha,
+ power=w_power,
+ ratio=params["ratio"],
+ alternative=alternative,
+ )
+ n1 = int(np.ceil(n1))
+ w_n = n1 + int(np.ceil(n1 * params["ratio"]))
+ elif atype == "anova":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ f_eff = params["effect_size"] * w_mult
+ n_per_g = FTestAnovaPower().solve_power(
+ effect_size=f_eff,
+ alpha=alpha,
+ power=w_power,
+ k_groups=params["k"],
+ )
+ n_per_g = int(np.ceil(n_per_g))
+ w_n = n_per_g * params["k"]
+ elif atype == "correlation":
+ import math
+
+ r_val = max(min(params["effect_size"] * w_mult, 0.99), 0.01)
+ fisher_z = math.atanh(r_val)
+ w_n = int(np.ceil(3 + ((z_alpha + w_z_beta) / fisher_z) ** 2))
+ elif atype == "regression":
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ adj_f2 = params["effect_size"] * w_mult
+ k_r = params["k"]
+ w_n = None
+ for nc in range(k_r + 2, 10000):
+ dfd = nc - k_r - 1
+ ncp = adj_f2 * nc
+ f_crit = f_dist.ppf(1 - alpha, k_r, dfd)
+ p_cur = 1 - noncentral_f.cdf(f_crit, k_r, dfd, ncp)
+ if p_cur >= w_power:
+ w_n = nc
+ break
+ elif atype == "logistic":
+ from statsmodels.stats.power import NormalIndPower
+
+ adj_or = max(params["or"] ** w_mult, 1.01)
+ ev_rate = params["event_rate"]
+ k_log = params["k"]
+ p1_log = (adj_or * ev_rate) / (1 - ev_rate + adj_or * ev_rate)
+ d_log = abs(p1_log - ev_rate)
+ p_bar = (ev_rate + p1_log) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ d_eff = d_log / se if se > 0 else 0
+ w_n = int(
+ np.ceil(
+ max(
+ NormalIndPower().solve_power(
+ effect_size=d_eff,
+ alpha=alpha,
+ power=w_power,
+ alternative=alternative,
+ ),
+ 10 * k_log,
+ )
+ )
+ )
+ elif atype == "chisq":
+ from statsmodels.stats.power import GofChisquarePower
+
+ adj_w = params["effect_size"] * w_mult
+ w_n = int(
+ np.ceil(
+ GofChisquarePower().solve_power(
+ effect_size=adj_w,
+ alpha=alpha,
+ power=w_power,
+ n_bins=params["df"] + 1,
+ )
+ )
+ )
+ elif atype == "mannwhitney":
+ from statsmodels.stats.power import NormalIndPower
+
+ adj_p = max(min(params["effect_size"] * w_mult, 0.99), 0.01)
+ d_mw = np.sqrt(3) * (adj_p - 0.5)
+ are = params["are"]
+ ratio = params["ratio"]
+ n1 = NormalIndPower().solve_power(
+ effect_size=d_mw,
+ alpha=alpha,
+ power=w_power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1 = int(np.ceil(n1 / are))
+ w_n = n1 + int(np.ceil(n1 * ratio))
+ elif atype == "logrank":
+ adj_hr = max(params["hr"] ** w_mult, 1.001)
+ ratio = params["ratio"]
+ med_ctrl = params["median_survival"]
+ study_dur = params["study_duration"]
+ log_hr = np.log(adj_hr)
+ num_events = (
+ ((z_alpha + w_z_beta) ** 2)
+ * ((ratio + 1) ** 2)
+ / (ratio * (log_hr**2))
+ )
+ lambda_ctrl = np.log(2) / med_ctrl
+ p_event_ctrl = 1 - np.exp(-lambda_ctrl * study_dur)
+ p_event_trt = 1 - np.exp(-lambda_ctrl / adj_hr * study_dur)
+ p_event = (p_event_ctrl + ratio * p_event_trt) / (1 + ratio)
+ w_n = int(np.ceil(num_events / p_event)) if p_event > 0 else None
+ elif atype == "cox":
+ adj_hr = max(params["hr"] ** w_mult, 1.001)
+ k = params["k"]
+ sd_x = params["sd_x"]
+ r2_x = params["r2_x"]
+ ev_rate = params["event_rate"]
+ log_hr = np.log(adj_hr)
+ var_denom = (sd_x**2) * (log_hr**2) * (1 - r2_x)
+ if var_denom > 0:
+ num_events = ((z_alpha + w_z_beta) ** 2) / var_denom
+ num_events = max(num_events, 10 * k)
+ w_n = int(np.ceil(num_events / ev_rate))
+ else:
+ w_n = None
+ elif atype == "equiv":
+ adj_margin = params["margin"] * w_mult
+ exp_diff = params["expected_diff"]
+ sd = params["sd"]
+ ratio = params["ratio"]
+ equiv_type = params.get("equiv_param_type", "Mean")
+ if equiv_type == "Proportion":
+ p1_eq = params.get("p1_eq", 0.2)
+ p2_eq = params.get("p2_eq", 0.2)
+ d_prop = abs(p1_eq - p2_eq)
+ d_e = adj_margin - d_prop
+ if d_e > 0:
+ p_bar = (p1_eq + p2_eq) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ es_e = d_e / se if se > 0 else 0
+ from statsmodels.stats.power import NormalIndPower
+
+ n1 = NormalIndPower().solve_power(
+ effect_size=es_e,
+ alpha=alpha,
+ power=w_power,
+ ratio=ratio,
+ alternative="larger",
+ )
+ n1 = int(np.ceil(n1))
+ w_n = n1 + int(np.ceil(n1 * ratio))
+ else:
+ w_n = None
+ else:
+ d_e = adj_margin - abs(exp_diff)
+ if d_e > 0:
+ es_e = d_e / sd
+ from statsmodels.stats.power import NormalIndPower
+
+ n1 = NormalIndPower().solve_power(
+ effect_size=es_e,
+ alpha=alpha,
+ power=w_power,
+ ratio=ratio,
+ alternative="larger",
+ )
+ n1 = int(np.ceil(n1))
+ w_n = n1 + int(np.ceil(n1 * ratio))
+ else:
+ w_n = None
+ elif atype == "rm_anova":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ adj_es = params["effect_size"] * w_mult
+ k = params["k"]
+ m = params["m"]
+ rho = params["rho"]
+ epsilon = params["epsilon"]
+ n_per_g = FTestAnovaPower().solve_power(
+ effect_size=adj_es, alpha=alpha, power=w_power, k_groups=k
+ )
+ design_effect = (1 + (m - 1) * rho) / m
+ df_adj = (m - 1) * epsilon
+ n_per_g_adj = max(
+ int(
+ np.ceil(
+ n_per_g * design_effect * k / (k * df_adj / (k - 1))
+ )
+ ),
+ int(np.ceil(n_per_g)),
+ )
+ w_n = n_per_g_adj * k
+ elif atype == "twoway_anova":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ rows = params["rows"]
+ cols = params["cols"]
+ focus = params["focus"]
+ if focus == "Main Effect A":
+ adj_es = params["f_a"] * w_mult
+ k_use = rows
+ elif focus == "Main Effect B":
+ adj_es = params["f_b"] * w_mult
+ k_use = cols
+ else:
+ adj_es = params["f_ab"] * w_mult
+ k_use = rows * cols
+ n_per_cell = int(
+ np.ceil(
+ FTestAnovaPower().solve_power(
+ effect_size=adj_es,
+ alpha=alpha,
+ power=w_power,
+ k_groups=k_use,
+ )
+ )
+ )
+ w_n = n_per_cell * rows * cols
+ elif atype == "roc_auc":
+ auc = params["auc"]
+ adj_auc = 0.5 + (auc - 0.5) * w_mult
+ adj_auc = max(min(adj_auc, 0.99), 0.51)
+ null_auc = params["null_auc"]
+ ratio = params["ratio"]
+ v_auc = adj_auc * (1 - adj_auc) + (ratio - 1) * (
+ adj_auc / (2 - adj_auc) - adj_auc**2
+ ) / (1 + ratio)
+ delta_auc = adj_auc - null_auc
+ if delta_auc > 0:
+ n_cases = int(
+ np.ceil(
+ ((z_alpha + w_z_beta) ** 2 * v_auc) / (delta_auc**2)
+ )
+ )
+ n_controls = int(np.ceil(n_cases * ratio))
+ w_n = n_cases + n_controls
+ else:
+ w_n = None
+ elif atype == "kappa":
+ adj_kappa = params["kappa"]
+ adj_kappa = 0.5 + (adj_kappa - 0.5) * w_mult
+ adj_kappa = max(min(adj_kappa, 0.99), 0.01)
+ null_kappa = params["null_kappa"]
+ adj_kappa = max(adj_kappa, null_kappa + 0.01)
+ delta_k = adj_kappa - null_kappa
+ if delta_k > 0:
+ w_n = int(
+ np.ceil(
+ (
+ (z_alpha + w_z_beta) ** 2
+ * null_kappa
+ * (1 - null_kappa)
+ )
+ / (delta_k**2)
+ )
+ )
+ w_n = max(w_n, params["raters"] * params["categories"] * 5)
+ else:
+ w_n = None
+ elif atype == "cluster_rct":
+ from statsmodels.stats.power import TTestIndPower
+
+ adj_d = params["effect_size"] * w_mult
+ icc = params["icc"]
+ cluster_m = params["cluster_size"]
+ ratio = params["ratio"]
+ deff = 1 + (cluster_m - 1) * icc
+ n1_ind = TTestIndPower().solve_power(
+ effect_size=adj_d,
+ alpha=alpha,
+ power=w_power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1_ind = int(np.ceil(n1_ind))
+ n1_clust = int(np.ceil(n1_ind * deff))
+ n1_clust = int(np.ceil(n1_clust / cluster_m)) * cluster_m
+ n2_clust = int(np.ceil(n1_clust * ratio))
+ w_n = n1_clust + n2_clust
+ elif atype == "precision":
+ adj_hw = params["half_width"] * w_mult
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and adj_hw > 0:
+ from scipy.stats import t as t_dist
+
+ w_n = int(
+ np.ceil(
+ (norm.ppf(1 - conf_alpha / 2) * sd / adj_hw) ** 2
+ )
+ )
+ w_n = max(w_n, 3)
+ for _ in range(20):
+ t_val = t_dist.ppf(1 - conf_alpha / 2, df=w_n - 1)
+ n_next = int(np.ceil((t_val * sd / adj_hw) ** 2))
+ n_next = max(n_next, 3)
+ if n_next == w_n:
+ break
+ w_n = n_next
+ else:
+ w_n = None
+ else:
+ prop = params["prop"]
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ w_n = (
+ int(np.ceil((z_hw**2 * prop * (1 - prop)) / (adj_hw**2)))
+ if prop > 0 and adj_hw > 0
+ else None
+ )
+ if w_n is not None:
+ w_n = max(w_n, 3)
+ elif atype == "pilot":
+ method = params["method"]
+ if method == "Rule of thumb":
+ w_n = int(np.ceil(params["n_per_group"] / max(w_mult, 0.1)))
+ elif method == "Precision-based":
+ conf_level = params["conf_level"]
+ conf_alpha = 1 - conf_level / 100
+ param_type = params["param_type"]
+ adj_hw = params["half_width"] * w_mult
+ if param_type == "Mean":
+ sd = params["sd"]
+ if sd > 0 and adj_hw > 0:
+ from scipy.stats import t as t_dist
+
+ w_n = int(
+ np.ceil(
+ (norm.ppf(1 - conf_alpha / 2) * sd / adj_hw)
+ ** 2
+ )
+ )
+ w_n = max(w_n, 3)
+ for _ in range(20):
+ t_val = t_dist.ppf(1 - conf_alpha / 2, df=w_n - 1)
+ n_next = int(np.ceil((t_val * sd / adj_hw) ** 2))
+ n_next = max(n_next, 3)
+ if n_next == w_n:
+ break
+ w_n = n_next
+ else:
+ w_n = None
+ else:
+ z_hw = norm.ppf(1 - conf_alpha / 2)
+ prop = params["prop"]
+ w_n = (
+ int(
+ np.ceil((z_hw**2 * prop * (1 - prop)) / (adj_hw**2))
+ )
+ if prop > 0 and adj_hw > 0
+ else None
+ )
+ if w_n is not None:
+ w_n = max(w_n, 5)
+ else:
+ main_n = params["main_n"]
+ fraction = params["fraction"]
+ w_n = max(int(np.ceil(main_n * fraction / max(w_mult, 0.1))), 5)
+ elif atype == "wilcoxon_sr":
+ p_val = params["effect_size"]
+ are = params["are"]
+ d_wsr = np.sqrt(3) * (p_val - 0.5) * 2
+ if d_wsr > 0:
+ from statsmodels.stats.power import NormalIndPower
+
+ w_n = int(
+ np.ceil(
+ NormalIndPower().solve_power(
+ effect_size=d_wsr,
+ alpha=alpha,
+ power=w_power,
+ alternative=alternative,
+ )
+ / are
+ )
+ )
+ else:
+ w_n = None
+ elif atype == "kruskal":
+ f_eff = params["effect_size"] * w_mult
+ k = params["k"]
+ are_kw = params.get("are", 0.955)
+ if f_eff > 0:
+ from statsmodels.stats.power import FTestAnovaPower
+
+ n_per_g = FTestAnovaPower().solve_power(
+ effect_size=f_eff, alpha=alpha, power=w_power, k_groups=k
+ )
+ n_per_g = int(np.ceil(n_per_g / are_kw))
+ w_n = n_per_g * k
+ else:
+ w_n = None
+ elif atype == "friedman":
+ w = params["w"]
+ k = params["k"]
+ m = params["m"]
+ are = params["are"]
+ adj_w = max(min(w * w_mult, 0.99), 0.01)
+ f_fr = np.sqrt(adj_w / (1 - adj_w))
+ if f_fr > 0:
+ from statsmodels.stats.power import FTestAnovaPower
+
+ n_per_g = FTestAnovaPower().solve_power(
+ effect_size=f_fr, alpha=alpha, power=w_power, k_groups=k
+ )
+ n_per_g = int(np.ceil(n_per_g * m / are))
+ w_n = n_per_g * k
+ else:
+ w_n = None
+ elif atype == "mcnemar":
+ p_b = params["p_b"]
+ p_c = params["p_c"]
+ p_b_adj = max(min(p_b * w_mult, 0.99), 0.01)
+ p_c_adj = max(min(p_c * w_mult, 0.99), 0.01)
+ d_mc = abs(p_b_adj - p_c_adj)
+ p_discordant = p_b_adj + p_c_adj
+ if d_mc > 0 and p_discordant > 0:
+ w_n = max(
+ int(
+ np.ceil(
+ ((z_alpha + w_z_beta) ** 2 * p_discordant)
+ / (d_mc**2)
+ )
+ ),
+ 5,
+ )
+ else:
+ w_n = None
+ elif atype == "fisher":
+ p1 = params["p1"]
+ p2 = params["p2"]
+ ratio = params["ratio"]
+ are_f = params.get("are", 0.833)
+ adj_p2 = p1 + (p2 - p1) * w_mult
+ adj_p2 = max(min(adj_p2, 0.99), 0.01)
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ d_eff = proportion_effectsize(adj_p2, p1)
+ if abs(d_eff) > 0:
+ n1 = NormalIndPower().solve_power(
+ effect_size=abs(d_eff),
+ alpha=alpha,
+ power=w_power,
+ ratio=ratio,
+ alternative=alternative,
+ )
+ n1 = int(np.ceil(n1 / are_f))
+ w_n = n1 + int(np.ceil(n1 * ratio))
+ else:
+ w_n = None
+ elif atype == "manova":
+ k = params["k"]
+ dv = params["dv"]
+ f2 = params["f2"]
+ rho = params["rho"]
+ adj_f2 = f2 * w_mult
+ manova_test = params.get("manova_test", "Pillai's Trace")
+ if adj_f2 > 0:
+ u = dv
+ v_num = k - 1
+ w_n = None
+ for n_try in range(k * dv + 2, 5000):
+ v_den = n_try - k - dv
+ if v_den <= 0:
+ continue
+ s_val = min(u, v_num)
+ df1 = u * v_num
+ if manova_test == "Pillai's Trace":
+ df2 = s_val * (v_den - dv + 1) + 4
+ elif manova_test == "Wilks' Lambda":
+ t_val = max(
+ np.sqrt(
+ (u**2 * v_num**2 - 4)
+ / max(u**2 + v_num**2 - 5, 1)
+ ),
+ 1,
+ )
+ df2 = (v_den - (u - v_num + 1) / 2) * t_val - (
+ u * v_num - 2
+ ) / 2
+ elif manova_test == "Hotelling-Lawley Trace":
+ df2 = s_val * (v_den - dv - 1) + 4
+ else:
+ df2 = s_val * (v_den - dv + 1) + 4
+ ncp = adj_f2 * n_try * (1 - rho)
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ f_crit = f_dist.ppf(1 - alpha, df1, df2)
+ p_cur = 1 - noncentral_f.cdf(f_crit, df1, df2, ncp)
+ if p_cur >= w_power:
+ w_n = n_try
+ break
+ else:
+ w_n = None
+ elif atype == "binomial":
+ p0 = params["p0"]
+ p1 = params["p1"]
+ p1_adj = p0 + (p1 - p0) * w_mult
+ p1_adj = max(min(p1_adj, 0.99), 0.01)
+ if p0 != p1_adj:
+ from scipy.stats import binom
+
+ w_n = None
+ for n_try in range(3, 10000):
+ if alternative == "two-sided":
+ alpha_lo = binom.ppf(alpha / 2, n_try, p0)
+ alpha_hi = binom.ppf(1 - alpha / 2, n_try, p0)
+ p_pow = binom.cdf(alpha_hi, n_try, p1_adj) - binom.cdf(
+ alpha_lo - 1, n_try, p1_adj
+ )
+ elif p1_adj > p0:
+ crit = binom.ppf(1 - alpha, n_try, p0)
+ p_pow = 1 - binom.cdf(crit - 1, n_try, p1_adj)
+ else:
+ crit = binom.ppf(alpha, n_try, p0)
+ p_pow = binom.cdf(crit, n_try, p1_adj)
+ if p_pow >= w_power:
+ w_n = n_try
+ break
+ else:
+ w_n = None
+ else:
+ w_n = None
+ except Exception:
+ w_n = None
+ row.append(w_n if w_n and w_n < 1000000 else None)
+ heatmap_data.append(row)
+
+ if any(any(r is not None for r in row) for row in heatmap_data):
+ fig_heat = go.Figure(
+ data=go.Heatmap(
+ z=heatmap_data,
+ x=[f"{m:.2f}×" for m in whatif_mults],
+ y=[f"Power = {p:.0%}" for p in whatif_powers],
+ text=[[str(v) if v else "—" for v in row] for row in heatmap_data],
+ texttemplate="%{text}",
+ colorscale="Blues",
+ hovertemplate="Power: %{y}
Effect Size: %{x}
N: %{text}",
+ )
+ )
+ fig_heat.update_layout(
+ template="plotly_dark",
+ height=250,
+ xaxis_title="Effect Size Multiplier",
+ yaxis_title="",
+ )
+ st.plotly_chart(fig_heat, use_container_width=True)
+ st.caption(
+ "Cells show required N at each power × effect size combination. Adjust your design assumptions accordingly."
+ )
+
+ # --- Inverse Power Analysis (Minimum Detectable Effect) ---
+ st.subheader("🎯 Sensitivity Analysis: Minimum Detectable Effect")
+ st.caption(
+ "Given a fixed sample size, what is the smallest effect size your study can detect?"
+ )
+
+ with st.expander(
+ "Enter a candidate sample size to compute the minimum detectable effect"
+ ):
+ c1, c2 = st.columns([1, 3])
+ with c1:
+ candidate_n = st.number_input(
+ "Candidate N (total)",
+ min_value=5,
+ max_value=100000,
+ value=n_total if n_total else 100,
+ step=10,
+ )
+ with c2:
+ st.caption(" ")
+ if st.button("Compute Minimum Detectable Effect", type="secondary"):
+ try:
+ from scipy.stats import t as t_dist_inv
+
+ mde = None
+ mde_label = ""
+ mde_note = ""
+
+ if atype == "one_mean":
+ from statsmodels.stats.power import TTestPower
+
+ mde = TTestPower().solve_power(
+ effect_size=None,
+ nobs=candidate_n,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ mde_label = "Cohen's d"
+ elif atype == "two_means":
+ from statsmodels.stats.power import TTestIndPower
+
+ ratio_val = params.get("ratio", 1)
+ n1 = candidate_n / (1 + ratio_val)
+ mde = TTestIndPower().solve_power(
+ effect_size=None,
+ nobs1=n1,
+ alpha=alpha,
+ power=power,
+ ratio=ratio_val,
+ alternative=alternative,
+ )
+ mde_label = "Cohen's d"
+ elif atype == "paired":
+ from statsmodels.stats.power import TTestPower
+
+ mde = TTestPower().solve_power(
+ effect_size=None,
+ nobs=candidate_n,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ mde_label = "Cohen's d_z"
+ elif atype == "one_prop":
+ from statsmodels.stats.power import NormalIndPower
+
+ d_eff = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=candidate_n,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ p0 = params.get("prop_null", 0.5)
+ mde = p0 + d_eff * np.sqrt(p0 * (1 - p0))
+ mde = max(0.01, min(0.99, mde))
+ mde_label = "Detectable proportion p₁"
+ mde_note = f"(null p₀ = {p0})"
+ elif atype == "two_prop":
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ ratio_val = params.get("ratio", 1)
+ n1 = candidate_n / (1 + ratio_val)
+ d_eff = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=n1,
+ alpha=alpha,
+ power=power,
+ ratio=ratio_val,
+ alternative=alternative,
+ )
+ p1_base = params.get("p1", 0.3)
+ mde = p1_base + d_eff * np.sqrt(2 * p1_base * (1 - p1_base))
+ mde = max(0.01, min(0.99, mde))
+ mde_label = "Detectable proportion p₂"
+ mde_note = f"(group 1 p₁ = {p1_base})"
+ elif atype == "anova":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ k_val = params.get("k", 3)
+ mde = FTestAnovaPower().solve_power(
+ effect_size=None,
+ nobs=candidate_n,
+ alpha=alpha,
+ power=power,
+ k_groups=k_val,
+ )
+ mde_label = "Cohen's f"
+ elif atype == "correlation":
+ import math
+
+ fisher_z = (
+ norm.ppf(
+ 1 - alpha / 2 if alternative == "two-sided" else 1 - alpha
+ )
+ + z_beta
+ ) / np.sqrt(candidate_n - 3)
+ mde = min(math.tanh(fisher_z), 0.99)
+ mde_label = "Pearson r"
+ elif atype == "chisq":
+ from statsmodels.stats.power import GofChisquarePower
+
+ df_val = params.get("df", 2)
+ mde = GofChisquarePower().solve_power(
+ effect_size=None,
+ nobs=candidate_n,
+ alpha=alpha,
+ power=power,
+ n_bins=df_val + 1,
+ )
+ mde_label = "Cohen's w"
+ elif atype == "mannwhitney":
+ from statsmodels.stats.power import NormalIndPower
+
+ are_val = params.get("are", 0.955)
+ ratio_val = params.get("ratio", 1)
+ n1_eff = candidate_n / (1 + ratio_val) * are_val
+ d_mw = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=n1_eff,
+ alpha=alpha,
+ power=power,
+ ratio=ratio_val,
+ alternative=alternative,
+ )
+ mde = 0.5 + d_mw / np.sqrt(3)
+ mde = max(0.51, min(0.99, mde))
+ mde_label = "P(X>Y)"
+ elif atype == "wilcoxon_sr":
+ from statsmodels.stats.power import NormalIndPower
+
+ are_val = params.get("are", 0.955)
+ d_z = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=candidate_n * are_val,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ mde = 0.5 + d_z / (2 * np.sqrt(3))
+ mde = max(0.51, min(0.99, mde))
+ mde_label = "Pr(positive diff)"
+ elif atype == "kruskal":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ are_val = params.get("are", 0.955)
+ k_val = params.get("k", 3)
+ mde = FTestAnovaPower().solve_power(
+ effect_size=None,
+ nobs=candidate_n * are_val,
+ alpha=alpha,
+ power=power,
+ k_groups=k_val,
+ )
+ mde_label = "Cohen's f"
+ elif atype == "mcnemar":
+ p_b = params.get("p_b", 0.2)
+ p_c = params.get("p_c", 0.4)
+ p_disc = p_b + p_c
+ if p_disc > 0:
+ delta = (
+ norm.ppf(
+ 1 - alpha / 2
+ if alternative == "two-sided"
+ else 1 - alpha
+ )
+ + z_beta
+ ) * np.sqrt(p_disc / candidate_n)
+ mde = max(delta, 0.01)
+ mde_label = "|p_b − p_c|"
+ elif atype == "fisher":
+ from statsmodels.stats.proportion import proportion_effectsize
+ from statsmodels.stats.power import NormalIndPower
+
+ are_val = params.get("are", 0.833)
+ ratio_val = params.get("ratio", 1)
+ p1_base = params.get("p1", 0.3)
+ n1 = candidate_n / (1 + ratio_val) * are_val
+ d_eff = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=n1,
+ alpha=alpha,
+ power=power,
+ ratio=ratio_val,
+ alternative=alternative,
+ )
+ mde = p1_base + d_eff * np.sqrt(2 * p1_base * (1 - p1_base))
+ mde = max(0.01, min(0.99, mde))
+ mde_label = "Detectable proportion p₂"
+ elif atype == "logrank":
+ ratio_val = params.get("ratio", 1)
+ hr_guess = params.get("hr", 2)
+ num_events_est = candidate_n * 0.5
+ log_hr_min = (
+ norm.ppf(
+ 1 - alpha / 2 if alternative == "two-sided" else 1 - alpha
+ )
+ + z_beta
+ ) * np.sqrt((ratio_val + 1) ** 2 / (ratio_val * num_events_est))
+ if log_hr_min > 0:
+ mde = np.exp(log_hr_min)
+ mde = max(1.001, min(10.0, mde))
+ mde_label = "Hazard Ratio (HR)"
+ mde_note = "(approximate, depends on event probability)"
+ elif atype == "cox":
+ sd_x = params.get("sd_x", 1)
+ r2_x = params.get("r2_x", 0)
+ ev_rate = params.get("event_rate", 0.5)
+ num_events_est = candidate_n * ev_rate
+ var_denom_inv = (
+ norm.ppf(
+ 1 - alpha / 2 if alternative == "two-sided" else 1 - alpha
+ )
+ + z_beta
+ ) ** 2 / num_events_est
+ if var_denom_inv > 0 and sd_x > 0:
+ log_hr_min = np.sqrt(var_denom_inv / (sd_x**2 * (1 - r2_x)))
+ mde = np.exp(log_hr_min)
+ mde = max(1.001, min(10.0, mde))
+ mde_label = "Hazard Ratio (HR)"
+ mde_note = "(approximate)"
+ elif atype == "logistic":
+ import math
+
+ ev_rate = params.get("event_rate", 0.3)
+ or_val = params.get("or", 2)
+ p1_log = (or_val * ev_rate) / (1 - ev_rate + or_val * ev_rate)
+ d_log = abs(p1_log - ev_rate)
+ p_bar = (ev_rate + p1_log) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ if se > 0:
+ d_eff = d_log / se
+ from statsmodels.stats.power import NormalIndPower
+
+ mde_d = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=candidate_n,
+ alpha=alpha,
+ power=power,
+ alternative=alternative,
+ )
+ p1_mde = ev_rate + mde_d * se
+ p1_mde = max(0.01, min(0.99, p1_mde))
+ mde = (p1_mde * (1 - ev_rate)) / (ev_rate * (1 - p1_mde))
+ mde = max(1.001, min(100.0, mde))
+ mde_label = "Odds Ratio (OR)"
+ mde_note = "(approximate)"
+ elif atype == "regression":
+ k_r = params.get("k", 3)
+ from scipy.stats import ncf as noncentral_f, f as f_dist
+
+ dfd = candidate_n - k_r - 1
+ if dfd > 0:
+ mde = None
+ for f2_try in np.linspace(0.001, 2.0, 2000):
+ ncp = f2_try * candidate_n
+ f_crit = f_dist.ppf(1 - alpha, k_r, dfd)
+ p_cur = 1 - noncentral_f.cdf(f_crit, k_r, dfd, ncp)
+ if p_cur >= power:
+ mde = f2_try
+ break
+ mde_label = "Cohen's f²"
+ else:
+ mde = None
+ elif atype == "binomial":
+ p0 = params.get("p0", 0.5)
+ from scipy.stats import binom
+
+ mde = None
+ for p1_try in np.linspace(p0 + 0.001, 0.99, 990):
+ if alternative == "two-sided":
+ alpha_lo = binom.ppf(alpha / 2, candidate_n, p0)
+ alpha_hi = binom.ppf(1 - alpha / 2, candidate_n, p0)
+ p_pow = binom.cdf(
+ alpha_hi, candidate_n, p1_try
+ ) - binom.cdf(alpha_lo - 1, candidate_n, p1_try)
+ elif p1_try > p0:
+ crit = binom.ppf(1 - alpha, candidate_n, p0)
+ p_pow = 1 - binom.cdf(crit - 1, candidate_n, p1_try)
+ else:
+ crit = binom.ppf(alpha, candidate_n, p0)
+ p_pow = binom.cdf(crit, candidate_n, p1_try)
+ if p_pow >= power:
+ mde = p1_try
+ break
+ mde_label = "Detectable proportion p₁"
+ elif atype == "equiv":
+ margin = params.get("margin", 1.0)
+ sd = params.get("sd", 1.0)
+ ratio_v = params.get("ratio", 1)
+ equiv_type = params.get("equiv_param_type", "Mean")
+ n1_eff = candidate_n / (1 + ratio_v)
+ if n1_eff > 1:
+ from statsmodels.stats.power import NormalIndPower
+
+ es_mde = NormalIndPower().solve_power(
+ effect_size=None,
+ nobs1=n1_eff,
+ alpha=alpha,
+ power=power,
+ ratio=ratio_v,
+ alternative="larger",
+ )
+ if equiv_type == "Proportion":
+ p_bar = (
+ params.get("p1_eq", 0.2) + params.get("p2_eq", 0.2)
+ ) / 2
+ se = np.sqrt(2 * p_bar * (1 - p_bar))
+ mde = es_mde * se
+ mde_label = "Detectable margin remaining (δ - |p₁−p₂|)"
+ else:
+ mde = es_mde * sd
+ mde_label = "Detectable margin remaining (δ - |d|)"
+ elif atype == "rm_anova":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ f_eff = params.get("effect_size", 0.25)
+ k_v = params.get("k", 2)
+ m_v = params.get("m", 3)
+ rho_v = params.get("rho", 0.5)
+ eps_v = params.get("epsilon", 0.75)
+ design_effect = (1 + (m_v - 1) * rho_v) / m_v
+ df_adj = (m_v - 1) * eps_v
+ n_eff = (
+ candidate_n
+ * design_effect
+ * (k_v - 1)
+ / (k_v * df_adj / (k_v - 1))
+ )
+ if n_eff > k_v:
+ mde = FTestAnovaPower().solve_power(
+ effect_size=None,
+ nobs=n_eff,
+ alpha=alpha,
+ power=power,
+ k_groups=k_v,
+ )
+ mde_label = "Cohen's f"
+ elif atype == "twoway_anova":
+ from statsmodels.stats.power import FTestAnovaPower
+
+ rows_v = params.get("rows", 2)
+ cols_v = params.get("cols", 2)
+ focus_v = params.get("focus", "Main Effect A")
+ k_use = (
+ rows_v
+ if focus_v == "Main Effect A"
+ else (cols_v if focus_v == "Main Effect B" else rows_v * cols_v)
+ )
+ n_per_cell_eff = candidate_n / (rows_v * cols_v)
+ if n_per_cell_eff > k_use:
+ mde = FTestAnovaPower().solve_power(
+ effect_size=None,
+ nobs=n_per_cell_eff,
+ alpha=alpha,
+ power=power,
+ k_groups=k_use,
+ )
+ mde_label = "Cohen's f"
+
+ if mde is not None and mde > 0:
+ st.success(
+ f"**Minimum detectable {mde_label}** with N = {candidate_n}: **{mde:.4f}** {mde_note}"
+ )
+ else:
+ st.warning(
+ "Could not compute minimum detectable effect for this analysis type."
+ )
+ except Exception as e:
+ st.error(f"Could not compute minimum detectable effect: {e}")
+
+ # --- Sample Size Justification ---
+ st.subheader("📝 Sample Size Justification")
+ st.caption(
+ "Copy the full protocol below for your grant application, IRB submission, or research protocol."
+ )
+
+ from datetime import datetime
+
+ now_str = datetime.now().strftime("%Y-%m-%d %H:%M")
+
+ atype_label = {
+ "one_mean": "One-sample Mean t/z-test",
+ "two_means": "Two Independent Means t-test",
+ "paired": "Paired Means t-test",
+ "one_prop": "One-sample Proportion",
+ "two_prop": "Two Proportions",
+ "anova": "One-way ANOVA",
+ "correlation": "Pearson Correlation",
+ "regression": "Multiple Linear Regression",
+ "logistic": "Logistic Regression",
+ "chisq": "Chi-Square Test",
+ "mannwhitney": "Mann-Whitney / Wilcoxon (Non-parametric)",
+ "logrank": "Log-Rank Test (Survival)",
+ "cox": "Cox Regression",
+ "equiv": "Equivalence / Non-Inferiority",
+ "rm_anova": "Repeated Measures ANOVA",
+ "twoway_anova": "Two-way / Factorial ANOVA",
+ "roc_auc": "ROC / AUC Analysis",
+ "kappa": "Cohen's Kappa / ICC Agreement",
+ "cluster_rct": "Cluster-RCT / Multilevel",
+ "precision": "Precision-based (CI Width)",
+ "pilot": "Pilot / Feasibility Study",
+ "wilcoxon_sr": "Wilcoxon Signed-Rank (paired)",
+ "kruskal": "Kruskal-Wallis Test",
+ "friedman": "Friedman Test",
+ "mcnemar": "McNemar's Test",
+ "fisher": "Fisher's Exact Test",
+ "manova": "MANOVA (Multivariate ANOVA)",
+ "binomial": "Binomial Exact Test",
+ "simulation": "Simulation-based Power (Monte Carlo)",
+ }
+
+ es_val = params.get(
+ "effect_size",
+ params.get("or", params.get("w", params.get("f2", params.get("f", "N/A")))),
+ )
+ es_str = (
+ f"{es_val:.4f}"
+ if isinstance(es_val, (int, float, np.integer, np.floating))
+ else str(es_val)
+ )
+ if n_per_group is not None:
+ n_desc = (
+ f"{n_per_group} per group"
+ if not isinstance(n_per_group, tuple)
+ else f"{n_per_group[0]} (Group 1) and {n_per_group[1]} (Group 2)"
+ )
+ protocol = f"""SAMPLE SIZE ESTIMATION PROTOCOL
+Generated: {now_str}
+Application: Statistical Test Finder (opencode.ai)
+
+Analysis: {atype_label.get(atype, atype)}
+Direction: {tails.lower()}
+Significance Level (α): {alpha * num_tests:.4f}"""
+ if num_tests > 1:
+ protocol += (
+ f" ({mc_method}-adjusted from {alpha_raw:.4f}, {num_tests} comparisons)"
+ )
+ protocol += f"""
+Statistical Power (1−β): {power:.0%}
+Effect Size: {es_str}
+"""
+ if isinstance(n_per_group, tuple):
+ protocol += f"""Allocation Ratio (n₂/n₁): {n_per_group[1] / n_per_group[0]:.2f}
+"""
+ protocol += f"""
+Required Sample Size: {n_desc}
+Total N: {n_total}"""
+ if dropout_rate > 0 and n_total_raw is not None and n_total_raw != n_total:
+ protocol += f"""
+Dropout Rate: {dropout_rate:.0%}
+Raw N (pre-dropout): {n_total_raw}
+Adjusted N (post-dropout): {n_total}
+"""
+ else:
+ protocol += "\n"
+ if cost_per > 0:
+ protocol += f"\nEstimated Study Cost: ${n_total * cost_per:,.0f}"
+ if recruitment_rate > 0:
+ protocol += f"\nEst. Recruitment Duration: {int(np.ceil(n_total / recruitment_rate))} months"
+ protocol += "\n"
+ else:
+ protocol = f"""SAMPLE SIZE ESTIMATION PROTOCOL
+Generated: {now_str}
+Application: Statistical Test Finder (opencode.ai)
+
+Analysis: {atype_label.get(atype, atype)}
+Direction: {tails.lower()}
+Significance Level (α): {alpha * num_tests:.4f}
+Statistical Power (1−β): {power:.0%}
+Effect Size: {es_str}
+Required Total N: {n_total}
"""
-num_points = st.slider("Number of points in spiral", 1, 10000, 1100)
-num_turns = st.slider("Number of turns in spiral", 1, 300, 31)
-
-indices = np.linspace(0, 1, num_points)
-theta = 2 * np.pi * num_turns * indices
-radius = indices
-
-x = radius * np.cos(theta)
-y = radius * np.sin(theta)
-
-df = pd.DataFrame({
- "x": x,
- "y": y,
- "idx": indices,
- "rand": np.random.randn(num_points),
-})
-
-st.altair_chart(alt.Chart(df, height=700, width=700)
- .mark_point(filled=True)
- .encode(
- x=alt.X("x", axis=None),
- y=alt.Y("y", axis=None),
- color=alt.Color("idx", legend=None, scale=alt.Scale()),
- size=alt.Size("rand", legend=None, scale=alt.Scale(range=[1, 150])),
- ))
\ No newline at end of file
+ fields_str = []
+ for k, v in params.items():
+ if k in (
+ "type",
+ "alpha",
+ "power",
+ "tails",
+ "dropout_rate",
+ "num_tests",
+ "cost_per",
+ "recruitment_rate",
+ ):
+ continue
+ fields_str.append(f" {k}: {v}")
+ if fields_str:
+ protocol += "\nFull Parameters:\n" + "\n".join(fields_str) + "\n"
+
+ protocol += "\n--- Generated by Statistical Test Finder ---"
+
+ if n_per_group is not None:
+ if isinstance(n_per_group, tuple):
+ n1, n2 = n_per_group
+ justification = (
+ f"A sample size of {n1} in Group 1 and {n2} in Group 2 "
+ f"(total N = {n_total}) was determined to provide {power:.0%} power "
+ f"at a significance level of α = {alpha * num_tests:.3f} "
+ )
+ if num_tests > 1:
+ justification += f"({mc_method}-adjusted for {num_tests} comparisons, per-test α = {alpha:.4f}) "
+ justification += (
+ f"to detect the anticipated effect size. "
+ f"Sample size estimation was performed using a {tails.lower()} {explanation.split('to detect')[0].strip().lower()}. "
+ )
+ if dropout_rate > 0:
+ justification += (
+ f"To account for an anticipated dropout rate of {dropout_rate:.0%}, "
+ f"the required sample was inflated from {n_total_raw} to {n_total} participants. "
+ )
+ justification += "All calculations were performed using the Statistical Test Finder application."
+ else:
+ justification = (
+ f"A sample size of {n_per_group} per group (total N = {n_total}) "
+ f"was determined to provide {power:.0%} power at a significance level of α = {alpha * num_tests:.3f} "
+ )
+ if num_tests > 1:
+ justification += f"({mc_method}-adjusted for {num_tests} comparisons, per-test α = {alpha:.4f}) "
+ justification += (
+ f"to detect the anticipated effect size. "
+ f"Sample size estimation was performed using a {tails.lower()} {explanation.split('to detect')[0].strip().lower()}. "
+ )
+ if dropout_rate > 0:
+ justification += (
+ f"To account for an anticipated dropout rate of {dropout_rate:.0%}, "
+ f"the required sample was inflated from {n_total_raw} to {n_total} participants. "
+ )
+ justification += "All calculations were performed using the Statistical Test Finder application."
+ else:
+ justification = f"A sample size of {n_total} participants was determined to provide {power:.0%} power at α = {alpha * num_tests:.3f}. All calculations were performed using the Statistical Test Finder application."
+
+ st.text(justification)
+
+ st.divider()
+ with st.expander("📄 Full Protocol Text"):
+ st.text_area(
+ ":orange[Full Protocol (select all, Ctrl+C / Cmd+C to copy)]",
+ protocol,
+ height=350,
+ )
+
+ # --- Key References ---
+ with st.expander("📚 Key References"):
+ st.markdown("""
+ **General Sample Size & Power:**
+ - Cohen, J. (1988). *Statistical Power Analysis for the Behavioral Sciences* (2nd ed.). Lawrence Erlbaum Associates.
+ - Faul, F., Erdfelder, E., Lang, A.-G., & Buchner, A. (2007). G*Power 3: A flexible statistical power analysis program. *Behavior Research Methods*, 39(2), 175–191.
+
+ **Means & t-tests:**
+ - Julious, S. A. (2004). Sample sizes for clinical trials with normal data. *Statistics in Medicine*, 23(12), 1921–1986.
+
+ **Proportions:**
+ - Fleiss, J. L., Levin, B., & Paik, M. C. (2003). *Statistical Methods for Rates and Proportions* (3rd ed.). Wiley.
+
+ **ANOVA & F-tests:**
+ - Cohen, J. (1988). *Statistical Power Analysis for the Behavioral Sciences* (2nd ed.). Lawrence Erlbaum Associates.
+
+ **Regression:**
+ - Green, S. B. (1991). How many subjects does it take to do a regression analysis? *Multivariate Behavioral Research*, 26(3), 499–510.
+ - Hsieh, F. Y., Bloch, D. A., & Larsen, M. D. (1998). A simple method of sample size calculation for linear and logistic regression. *Statistics in Medicine*, 17(14), 1623–1634.
+
+ **Correlation:**
+ - Fisher, R. A. (1915). Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population. *Biometrika*, 10(4), 507–521.
+
+ **Survival Analysis:**
+ - Schoenfeld, D. (1983). Sample-size formula for the proportional-hazards regression model. *Biometrics*, 39(2), 499–503.
+ - Freedman, L. S. (1982). Tables of the number of patients required in clinical trials using the logrank test. *Statistics in Medicine*, 1(2), 121–129.
+
+ **Non-parametric Tests:**
+ - Lehmann, E. L. (2006). *Nonparametrics: Statistical Methods Based on Ranks*. Springer.
+
+ **Equivalence / Non-Inferiority:**
+ - Blackwelder, W. C. (1982). "Proving the null hypothesis" in clinical trials. *Controlled Clinical Trials*, 3(4), 345–353.
+
+ **Cluster-RCT:**
+ - Donner, A., & Klar, N. (2000). *Design and Analysis of Cluster Randomization Trials in Health Research*. Arnold.
+
+ **Repeated Measures:**
+ - Greenhouse, S. W., & Geisser, S. (1959). On methods in the analysis of profile data. *Psychometrika*, 24(2), 95–112.
+
+ **ROC / AUC:**
+ - Obuchowski, N. A. (1994). Sample size calculations in studies of test accuracy. *Statistical Methods in Medical Research*, 7(4), 371–392.
+
+ **Kappa / Agreement:**
+ - Cantor, A. B. (1996). Sample-size calculations for Cohen's kappa. *Psychological Methods*, 1(2), 150–153.
+
+ **Pilot Studies:**
+ - Julious, S. A. (2005). Sample size of 12 per group rule of thumb for a pilot study. *Pharmaceutical Statistics*, 4(4), 287–291.
+ - Whitehead, A. L., et al. (2016). Estimating the sample size for a pilot randomised trial to minimise the overall trial sample size. *Journal of Clinical Epidemiology*, 71, 23–29.
+
+ **Multiple Testing:**
+ - Bonferroni, C. E. (1936). Teoria statistica delle classi e calcolo delle probabilità. *Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze*, 8, 3–62.
+ """)