""" power.py — 검정력 시뮬레이션. - z-test 기반 두 비율 비교(예: NASH 해소율 placebo vs treatment). - 컴플라이언스(입력률) 보정: effective N = nominal N × compliance. - monte carlo 옵션으로 power 추정. CLI 진입: python3 power.py --target-effect 0.15 --p-control 0.10 --alpha 0.05 \ --power 0.80 --compliance 0.85 """ from __future__ import annotations import argparse import math import random def _norm_ppf(p: float) -> float: """역표준정규 CDF 근사 (Beasley-Springer-Moro 단순판).""" if not (0.0 < p < 1.0): raise ValueError("p must be in (0,1)") a = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00] b = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, 6.680131188771972e+01, -1.328068155288572e+01] c = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00] d = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00] plow = 0.02425 phigh = 1 - plow if p < plow: q = math.sqrt(-2 * math.log(p)) return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \ ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) if p <= phigh: q = p - 0.5 r = q*q return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / \ (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1) q = math.sqrt(-2 * math.log(1 - p)) return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \ ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) def _norm_cdf(x: float) -> float: return 0.5 * (1.0 + math.erf(x / math.sqrt(2))) def sample_size_two_proportions(p1: float, p2: float, alpha: float = 0.05, power: float = 0.80, two_sided: bool = True) -> int: """두 비율 비교에 필요한 군당 표본수 (간단 z-test).""" if p1 == p2: return float("inf") z_alpha = _norm_ppf(1 - alpha / (2 if two_sided else 1)) z_beta = _norm_ppf(power) p_bar = (p1 + p2) / 2 num = (z_alpha * math.sqrt(2 * p_bar * (1 - p_bar)) + z_beta * math.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) ** 2 n = num / ((p1 - p2) ** 2) return math.ceil(n) def adjust_for_compliance(n_per_arm: int, compliance: float, dropout: float = 0.0) -> int: """컴플라이언스+dropout 보정. effective n = nominal × c × (1-d).""" if compliance <= 0: return float("inf") eff_factor = compliance * (1 - dropout) if eff_factor <= 0: return float("inf") return math.ceil(n_per_arm / eff_factor) def compliance_scenarios(p_control: float, p_treatment: float, alpha: float = 0.05, power: float = 0.80, dropout: float = 0.10) -> list[dict]: base_n = sample_size_two_proportions(p_control, p_treatment, alpha, power) out = [] for c in [0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60]: adj = adjust_for_compliance(base_n, c, dropout) out.append({ "compliance": c, "dropout": dropout, "n_per_arm_nominal": base_n, "n_per_arm_required": adj, "total_required": adj * 2, }) return out def power_at_n(n_per_arm: int, p1: float, p2: float, alpha: float = 0.05, two_sided: bool = True) -> float: if n_per_arm <= 0: return 0.0 z_alpha = _norm_ppf(1 - alpha / (2 if two_sided else 1)) p_bar = (p1 + p2) / 2 se_null = math.sqrt(2 * p_bar * (1 - p_bar) / n_per_arm) se_alt = math.sqrt((p1 * (1 - p1) + p2 * (1 - p2)) / n_per_arm) if se_alt == 0: return 1.0 if p1 != p2 else 0.0 z = (abs(p2 - p1) - z_alpha * se_null) / se_alt return _norm_cdf(z) def monte_carlo_power(n_per_arm: int, p1: float, p2: float, alpha: float = 0.05, n_sims: int = 2000, seed: int = 42) -> float: rng = random.Random(seed) z_crit = _norm_ppf(1 - alpha / 2) sig = 0 for _ in range(n_sims): x1 = sum(1 for _ in range(n_per_arm) if rng.random() < p1) x2 = sum(1 for _ in range(n_per_arm) if rng.random() < p2) ph1 = x1 / n_per_arm ph2 = x2 / n_per_arm p_pool = (x1 + x2) / (2 * n_per_arm) if p_pool in (0, 1): continue se = math.sqrt(2 * p_pool * (1 - p_pool) / n_per_arm) if se == 0: continue z = (ph2 - ph1) / se if abs(z) >= z_crit: sig += 1 return sig / n_sims def cli(): ap = argparse.ArgumentParser(description="MASHePRO-Kor power simulation") ap.add_argument("--p-control", type=float, default=0.10, help="대조군 NASH 해소율 (default 0.10)") ap.add_argument("--p-treatment", type=float, default=0.25, help="치료군 NASH 해소율 (default 0.25)") ap.add_argument("--alpha", type=float, default=0.05) ap.add_argument("--power", type=float, default=0.80) ap.add_argument("--dropout", type=float, default=0.10) ap.add_argument("--compliance", type=float, default=None, help="단일 시나리오 입력률; 미지정 시 다단 시나리오 표") ap.add_argument("--n-current", type=int, default=None) ap.add_argument("--target-n", type=int, default=None) ap.add_argument("--monte-carlo", action="store_true") ap.add_argument("--n-sims", type=int, default=2000) args = ap.parse_args() print(f"두 비율 비교: control={args.p_control}, treatment={args.p_treatment}") base_n = sample_size_two_proportions( args.p_control, args.p_treatment, args.alpha, args.power) print(f"군당 nominal sample size (alpha={args.alpha}, power={args.power}): {base_n}") print(f"총 nominal: {base_n * 2}") if args.compliance is not None: adj = adjust_for_compliance(base_n, args.compliance, args.dropout) print(f"\n컴플라이언스 {args.compliance:.0%}, dropout {args.dropout:.0%} 보정:") print(f" → 군당 필요 N: {adj}, 총 {adj*2}") else: print("\n컴플라이언스 시나리오 표:") print(f" {'comp':>6} | {'n/arm nominal':>14} | {'n/arm req':>10} | {'total':>7}") for row in compliance_scenarios(args.p_control, args.p_treatment, args.alpha, args.power, args.dropout): print(f" {row['compliance']:>6.0%} | {row['n_per_arm_nominal']:>14} | " f"{row['n_per_arm_required']:>10} | {row['total_required']:>7}") if args.n_current: pw = power_at_n(args.n_current, args.p_control, args.p_treatment, args.alpha) print(f"\n현재 군당 N={args.n_current} → 추정 power = {pw:.3f}") if args.monte_carlo: mc = monte_carlo_power(args.n_current, args.p_control, args.p_treatment, args.alpha, args.n_sims) print(f" Monte Carlo ({args.n_sims} sims): power = {mc:.3f}") if args.target_n: pw = power_at_n(args.target_n, args.p_control, args.p_treatment, args.alpha) print(f"\n목표 군당 N={args.target_n} → 추정 power = {pw:.3f}") if __name__ == "__main__": cli()