diff --git "a/src/streamlit_app.py" "b/src/streamlit_app.py" --- "a/src/streamlit_app.py" +++ "b/src/streamlit_app.py" @@ -1,4285 +1,344 @@ import streamlit as st +import plotly.graph_objects as go import numpy as np import pandas as pd -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, +# ============================================================================= +# 1. CONFIGURATION & DATA +# ============================================================================= + +st.set_page_config( + page_title="Biomaterials Stress-Strain Explorer", + page_icon="🦷", + layout="wide", + initial_sidebar_state="expanded" +) + +# Material Database: Dental & Biomaterials +MATERIALS = { + "🦷 Dental Amalgam": { + "E": 80000, "yield": 200, "uts": 260, "strain": 0.006, "type": "brittle", + "desc": "Traditional restorative alloy. High stiffness, low ductility.", + "color": "#4A90D9" + }, + "🦷 Composite Resin": { + "E": 12000, "yield": 90, "uts": 120, "strain": 0.035, "type": "semi-ductile", + "desc": "Tooth-colored restorative. Moderate stiffness and ductility.", + "color": "#E8A838" + }, + "🦷 Zirconia (Y-TZP)": { + "E": 200000, "yield": 900, "uts": 1100, "strain": 0.003, "type": "brittle", + "desc": "High-strength ceramic. Extremely stiff, very brittle failure.", + "color": "#9B59B6" + }, + "🦷 Gold Alloy (Type IV)": { + "E": 85000, "yield": 350, "uts": 500, "strain": 0.12, "type": "ductile", + "desc": "Cast noble metal. Excellent ductility and fatigue resistance.", + "color": "#F39C12" + }, + "🦷 Titanium (Grade 5)": { + "E": 110000, "yield": 850, "uts": 950, "strain": 0.10, "type": "ductile", + "desc": "Ti-6Al-4V alloy. Standard for implants. High strength-to-weight.", + "color": "#3498DB" + }, + "🦷 Glass Ionomer": { + "E": 8000, "yield": 40, "uts": 55, "strain": 0.008, "type": "brittle", + "desc": "Adhesive cement. Low strength, fluoride release.", + "color": "#1ABC9C" + }, + "🦴 Cortical Bone": { + "E": 17000, "yield": 120, "uts": 150, "strain": 0.018, "type": "semi-ductile", + "desc": "Dense outer bone. Anisotropic, moderate stiffness.", + "color": "#ECF0F1" + }, + "🧬 PEEK": { + "E": 3600, "yield": 95, "uts": 100, "strain": 0.06, "type": "ductile", + "desc": "Polyetheretherketone. Used in spinal implants. Low modulus matches bone.", + "color": "#E74C3C" + }, + "🧬 UHMWPE": { + "E": 800, "yield": 20, "uts": 30, "strain": 0.35, "type": "ductile", + "desc": "Ultra-high molecular weight polyethylene. Joint bearing surface.", + "color": "#F1C40F" + }, + "⚙️ Stainless Steel 316L": { + "E": 193000, "yield": 205, "uts": 505, "strain": 0.40, "type": "ductile", + "desc": "Austenitic steel. Surgical instruments and temporary implants.", + "color": "#7F8C8D" } - 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}} -""" - - 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}}}}" +# ============================================================================= +# 2. PHYSICS ENGINE (Curve Generation) +# ============================================================================= - 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) +def generate_stress_strain(E, sigma_y, sigma_uts, epsilon_break, material_type): + """ + Generates a realistic stress-strain curve using a piecewise model. + 1. Elastic: Linear (Hooke's Law) + 2. Plastic: Cosine transition to UTS + 3. Necking: Power-law decay to fracture + """ + # Calculate yield strain + epsilon_y = sigma_y / E + + # Determine strain at UTS based on material type + if material_type == "brittle": + epsilon_uts = epsilon_y * 1.1 + elif material_type == "semi-ductile": + epsilon_uts = epsilon_y + (epsilon_break - epsilon_y) * 0.5 + else: # ductile + epsilon_uts = epsilon_y + (epsilon_break - epsilon_y) * 0.4 + + # Generate strain points + n_points = 1000 + epsilon = np.linspace(0, epsilon_break, n_points) + sigma = np.zeros_like(epsilon) + + for i, eps in enumerate(epsilon): + if eps <= epsilon_y: + # Elastic Region + sigma[i] = E * eps + elif eps <= epsilon_uts: + # Strain Hardening Region (Smooth transition) + # Normalize position between yield and UTS + t = (eps - epsilon_y) / (epsilon_uts - epsilon_y) + # Use cosine for smooth hardening + sigma[i] = sigma_y + (sigma_uts - sigma_y) * (0.5 * (1 - np.cos(np.pi * t))) 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", + # Necking / Softening Region + t = (eps - epsilon_uts) / (epsilon_break - epsilon_uts) + # Stress drops from UTS to near 0 + sigma[i] = sigma_uts * (1 - 0.95 * (t ** 0.7)) + + return epsilon, np.maximum(sigma, 0) + +# ============================================================================= +# 3. ACCESSIBILITY HELPERS +# ============================================================================= + +def inject_accessibility_css(): + """Injects CSS for high contrast and large text modes.""" + st.markdown(""" + + """, unsafe_allow_html=True) - 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) +# ============================================================================= +# 4. MAIN APP LOGIC +# ============================================================================= - elif analysis_mode == "Compromise": - col1, col2, col3 = st.columns(3) +def main(): + inject_accessibility_css() + + # --- Sidebar: Controls & Accessibility --- + with st.sidebar: + st.header("⚙️ Controls") + + # Material Selection + selected_material = st.selectbox( + "Select Material Preset", + options=list(MATERIALS.keys()), + index=0 + ) + + st.markdown("---") + st.header("📊 Parameters") + + # Get current preset values + preset = MATERIALS[selected_material] + + # Widgets + col1, col2 = st.columns(2) with col1: - st.metric("Total Sample Size (N)", n_total) + E = st.slider("Young's Modulus (MPa)", 100, 300000, int(preset["E"]), 100) + sigma_y = st.slider("Yield Strength (MPa)", 10, 1500, int(preset["yield"]), 10) 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" - + sigma_uts = st.slider("UTS (MPa)", 10, 1500, int(preset["uts"]), 10) + epsilon_break = st.slider("Strain at Break", 0.001, 0.50, float(preset["strain"]), 0.001) + + st.markdown("---") + st.header("♿ Accessibility") + + # Accessibility Toggles + high_contrast = st.toggle("High Contrast Mode", value=False) + large_text = st.toggle("Large Text Mode", value=False) + screen_reader_mode = st.toggle("Verbose Descriptions", value=False) + + # Apply CSS classes based on toggles + css_classes = [] + if high_contrast: css_classes.append("high-contrast") + if large_text: css_classes.append("large-text") + + if css_classes: + st.markdown(f'
', unsafe_allow_html=True) + + # --- Main Content --- + + # Title + st.title("Biomaterials Stress-Strain Explorer") + + # Description Area + st.subheader(f"Material: {selected_material}") + st.info(preset["desc"]) + + # Screen Reader Description (Hidden visually but readable) + if screen_reader_mode: + st.markdown(f""" +

+ Chart displays stress-strain curve for {selected_material}. + Young's Modulus is {E} MPa. + Yield Strength is {sigma_y} MPa. + Ultimate Tensile Strength is {sigma_uts} MPa. + Fracture occurs at strain {epsilon_break}. +

+ """, unsafe_allow_html=True) + + # --- Calculation --- + strain, stress = generate_stress_strain(E, sigma_y, sigma_uts, epsilon_break, preset["type"]) + + # --- Plotly Chart --- 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 - - 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}", - ) + + # 1. The Stress-Strain Curve + fig.add_trace(go.Scatter( + x=strain, + y=stress, + mode='lines', + name='Stress-Strain', + line=dict(color=preset["color"], width=4), + hovertemplate="Strain: %{x:.3f}
Stress: %{y:.1f} MPa" + )) + + # 2. Elastic Modulus Line (Dashed, up to yield) + epsilon_y = sigma_y / E + fig.add_trace(go.Scatter( + x=[0, epsilon_y], + y=[0, sigma_y], + mode='lines', + name='Elastic Modulus', + line=dict(color='red', width=2, dash='dash'), + hoverinfo="skip" + )) + + # 3. Annotations for Key Points + annotations = [] + + # Yield Point + annotations.append(dict( + x=epsilon_y, y=sigma_y, + xref="x", yref="y", + text=f"Yield Point
σ={sigma_y} MPa", + showarrow=True, arrowhead=2, arrowsize=1, arrowwidth=2, arrowcolor="orange", + bordercolor="orange", borderwidth=1, borderpad=4, bgcolor="white", + font=dict(family="Arial", size=12, color="orange") + )) + + # UTS Point + # Find index closest to UTS for accurate annotation placement + uts_idx = np.argmax(stress) + uts_strain = strain[uts_idx] + annotations.append(dict( + x=uts_strain, y=sigma_uts, + xref="x", yref="y", + text=f"UTS
σ={sigma_uts} MPa", + showarrow=True, arrowhead=2, arrowsize=1, arrowwidth=2, arrowcolor="purple", + bordercolor="purple", borderwidth=1, borderpad=4, bgcolor="white", + font=dict(family="Arial", size=12, color="purple") + )) + + # Fracture Point + annotations.append(dict( + x=epsilon_break, y=0, + xref="x", yref="y", + text=f"Fracture
ε={epsilon_break}", + showarrow=True, arrowhead=2, arrowsize=1, arrowwidth=2, arrowcolor="black", + bordercolor="black", borderwidth=1, borderpad=4, bgcolor="white", + font=dict(family="Arial", size=12, color="black") + )) + 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." + annotations=annotations, + xaxis_title="Strain (ε) - Dimensionless", + yaxis_title="Stress (σ) - MPa", + title="Stress-Strain Curve", + template="plotly_white" if not high_contrast else "plotly_dark", + hovermode="x unified", + legend=dict(x=0.01, y=0.99, bgcolor="rgba(255,255,255,0.8)"), + margin=dict(l=50, r=50, t=50, b=50) ) + + # Accessibility: If high contrast, force dark background on plot + if high_contrast: + fig.update_layout(paper_bgcolor="black", plot_bgcolor="black") + fig.update_xaxes(gridcolor="gray", zerolinecolor="white", tickfont=dict(color="white")) + fig.update_yaxes(gridcolor="gray", zerolinecolor="white", tickfont=dict(color="white")) - 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} -""" - - 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. - """) - - -st.set_page_config(page_title="Power Analysis & Sample Size", page_icon="⚡", layout="wide") - -def main(): - _render_power_analysis() - render_footer() - + st.plotly_chart(fig, use_container_width=True) + + # --- Data Table (For Accessibility/Reference) --- + st.subheader("📋 Key Mechanical Properties") + + df_props = pd.DataFrame({ + "Property": ["Young's Modulus", "Yield Strength", "Ultimate Tensile Strength", "Strain at Break", "Yield Strain", "Toughness (Approx)"], + "Value": [ + f"{E} MPa", + f"{sigma_y} MPa", + f"{sigma_uts} MPa", + f"{epsilon_break}", + f"{epsilon_y:.4f}", + f"{np.trapz(stress, strain):.1f} MJ/m³" + ], + "Description": [ + "Stiffness in elastic region", + "Onset of plastic deformation", + "Maximum stress before necking", + "Strain at fracture", + "Calculated as σ_y / E", + "Area under curve (Energy absorbed)" + ] + }) + + st.dataframe(df_props, use_container_width=True, hide_index=True) + + # Close the div if we opened it for CSS classes + if css_classes: + st.markdown("
", unsafe_allow_html=True) if __name__ == "__main__": main()