""" core.py — GlycemicNIMargin-Kor 핵심 계산 로직 (순수 함수) 당뇨 RCT 비열등성(NI) 마진 산출 + MMRM 기반 탈락보정 표본수 계산. 의존성: numpy, scipy 만 사용 (statsmodels 의존 금지). 표본수 공식은 closed-form, 일부는 Monte-Carlo로 직접 구현. 면책: 본 모듈은 참고용·연구용 도구입니다. 실제 임상시험 설계는 반드시 생물통계학자(biostatistician) 검토를 거쳐야 합니다. 수식 출처(주석에 명시): - 비열등성 2-표본 t-test 표본수: Chow SC, Shao J, Wang H. "Sample Size Calculations in Clinical Research" 2nd ed. (2008), Ch.3. n/arm = (z_{1-alpha} + z_{1-beta})^2 * (2 sigma^2) / (margin - |true_diff|)^2 - 95-95 synthesis(2-confidence-interval) 방법: ICH E10 (2001), Hung HMJ et al. "A regulatory perspective on choice of margin and statistical inference issue in non-inferiority trials" (2005). - MMRM 설계효과(longitudinal efficiency): Lu K, Luo X, Chen PY. "Sample size estimation for repeated measures analysis in randomized clinical trials with missing data" Int J Biostat (2008). Hedeker D, Gibbons RD. "Longitudinal Data Analysis" (2006), Ch.13. - 다중성(공동 1차 endpoint): Hochberg(1988), Dmitrienko et al. "Multiple Testing Problems in Pharmaceutical Statistics" (2010). - 탈락 MAR 보정: 완료자 비율 기반 inflation factor, EMA Missing Data guideline (CPMP/EWP/1776/99 Rev.1) 원칙. """ import json import math import os try: import numpy as np HAS_NUMPY = True except Exception: # pragma: no cover HAS_NUMPY = False np = None try: from scipy import stats as _sps HAS_SCIPY = True except Exception: # pragma: no cover HAS_SCIPY = False _sps = None DISCLAIMER = ( "[면책] 본 도구는 참고용·연구용 RCT 설계 계산기입니다. " "실제 임상시험의 표본수·마진·통계 설계는 반드시 생물통계학자 검토가 필요하며, " "최신 규제 가이던스(ICH E9/E10, FDA/EMA/MFDS)를 확인하십시오." ) # -------------------------------------------------------------------------- # 정규분포 유틸 (scipy 없으면 표준라이브러리 근사로 graceful degradation) # -------------------------------------------------------------------------- def norm_ppf(p): """표준정규 분위수(inverse CDF). scipy 있으면 정확, 없으면 Acklam 근사.""" if HAS_SCIPY: return float(_sps.norm.ppf(p)) # Peter Acklam의 유리근사 (정확도 ~1e-9) if p <= 0.0: return -math.inf if p >= 1.0: return math.inf 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 = 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) 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) def norm_cdf(x): """표준정규 CDF.""" if HAS_SCIPY: return float(_sps.norm.cdf(x)) return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0))) # -------------------------------------------------------------------------- # 프리셋 로더 # -------------------------------------------------------------------------- def load_presets(path=None): """margin_presets.json 로드. 경로 미지정 시 모듈 옆 data/ 에서 찾음.""" if path is None: here = os.path.dirname(os.path.abspath(__file__)) path = os.path.join(here, "data", "margin_presets.json") with open(path, "r", encoding="utf-8") as f: return json.load(f) # ========================================================================== # 기능 1: NI 마진 정당화 엔진 # ========================================================================== def ni_margin_fixed(active_vs_placebo_effect, preserved_fraction=0.5): """ Fixed-margin (M1/M2) 방법. M1 = 활성대조약이 위약 대비 보였던 효과(과거 시험 점추정 또는 보수적 하한). M2 = NI 마진 = M1 * (1 - preserved_fraction) -> 보존비율 만큼은 반드시 지킴. 당뇨 HbA1c 예: 활성약-위약 차이 0.8% 이고 효과의 50% 보존이면 M2 = 0.8 * 0.5 = 0.4% 가 NI 마진. Parameters ---------- active_vs_placebo_effect : float (M1, 절대값) preserved_fraction : float 보존할 활성약 효과 비율(예 0.5 = 50%) Returns: dict """ m1 = abs(float(active_vs_placebo_effect)) pf = float(preserved_fraction) m2 = m1 * (1.0 - pf) return { "method": "fixed_margin", "M1_active_vs_placebo": m1, "preserved_fraction": pf, "preserved_effect": m1 * pf, "NI_margin_M2": m2, "interpretation": ( f"활성대조약의 위약 대비 효과(M1)는 {m1:g}. " f"이 효과의 {pf*100:.0f}%({m1*pf:g})를 보존하도록 요구하면 " f"NI 마진(M2)은 {m2:g}." ), } def ni_margin_synthesis(active_vs_placebo_effect, se_historical, preserved_fraction=0.5, conf=0.95): """ 95-95 synthesis(2-CI) 방법. 과거 위약대조 시험의 활성약-위약 효과 추정치와 그 표준오차로 효과의 보수적 하한(historical lower CI)을 구하고, 거기에 보존비율을 적용. M_eff_lower = effect - z_{conf} * se_historical NI_margin = M_eff_lower * (1 - preserved_fraction) 출처: ICH E10; Hung et al. (2005) synthesis method. """ eff = abs(float(active_vs_placebo_effect)) se = abs(float(se_historical)) pf = float(preserved_fraction) z = norm_ppf(conf) eff_lower = max(eff - z * se, 0.0) margin = eff_lower * (1.0 - pf) return { "method": "synthesis_95_95", "active_vs_placebo_effect": eff, "se_historical": se, "historical_lower_CI": eff_lower, "conf": conf, "preserved_fraction": pf, "NI_margin": margin, "interpretation": ( f"과거 위약대조 시험에서 활성약 효과 점추정 {eff:g}, " f"SE {se:g} -> {conf*100:.0f}% 신뢰하한 {eff_lower:g}. " f"효과의 {pf*100:.0f}% 보존 시 NI 마진 = {margin:g}. " "(synthesis 방법: 과거시험 불확실성을 분산에 합산)" ), } def generate_margin_justification_text(endpoint_label, unit, margin, method_result, lang="ko"): """ICH E10 용어로 마진 정당화 문단 자동 작성 (국문/영문).""" m = method_result if lang == "en": if m["method"] == "fixed_margin": return ( f"Non-inferiority margin justification ({endpoint_label}). " f"The non-inferiority margin (delta) was set to {margin:g} {unit} " f"using the fixed-margin (M2) approach per ICH E10. The active " f"comparator's effect over placebo in historical placebo-controlled " f"trials (M1) was {m['M1_active_vs_placebo']:g} {unit}; the trial is " f"designed to preserve {m['preserved_fraction']*100:.0f}% of this " f"effect, yielding delta = M1 x (1 - preserved fraction) = {margin:g} " f"{unit}. This ensures the test treatment retains a clinically " f"relevant fraction of the established active-comparator benefit and " f"is superior to (a putative) placebo (assay sensitivity / constancy " f"assumption per ICH E10)." ) else: return ( f"Non-inferiority margin justification ({endpoint_label}). " f"The margin (delta = {margin:g} {unit}) was derived using the " f"95-95 synthesis approach (ICH E10; Hung et al. 2005). The " f"historical active-vs-placebo effect estimate was " f"{m['active_vs_placebo_effect']:g} {unit} (SE {m['se_historical']:g}); " f"its {m['conf']*100:.0f}% lower confidence bound was " f"{m['historical_lower_CI']:g} {unit}. Preserving " f"{m['preserved_fraction']*100:.0f}% of this conservative effect " f"estimate yields delta = {margin:g} {unit}. Uncertainty in the " f"historical estimate is propagated into the variance of the " f"non-inferiority test statistic." ) # 국문 if m["method"] == "fixed_margin": return ( f"비열등성 마진 정당화 ({endpoint_label}). " f"비열등성 마진(δ)은 ICH E10에 따른 fixed-margin(M2) 방법으로 " f"{margin:g} {unit}로 설정하였다. 과거 위약대조 시험에서 활성대조약의 " f"위약 대비 효과(M1)는 {m['M1_active_vs_placebo']:g} {unit}였으며, " f"본 시험은 이 효과의 {m['preserved_fraction']*100:.0f}%를 보존하도록 " f"설계되어 δ = M1 × (1 - 보존비율) = {margin:g} {unit}로 산출되었다. " f"이는 시험약이 확립된 활성대조약 효과의 임상적으로 의미있는 부분을 " f"유지하고 (가상의) 위약보다 우월함을 보장하며, ICH E10의 검정 민감도" f"(assay sensitivity) 및 항상성(constancy) 가정에 부합한다." ) else: return ( f"비열등성 마진 정당화 ({endpoint_label}). " f"마진(δ = {margin:g} {unit})은 95-95 synthesis 방법(ICH E10; " f"Hung 등 2005)으로 도출하였다. 과거 활성약-위약 효과 점추정치는 " f"{m['active_vs_placebo_effect']:g} {unit}(SE {m['se_historical']:g})였고, " f"그 {m['conf']*100:.0f}% 신뢰하한은 {m['historical_lower_CI']:g} {unit}였다. " f"이 보수적 효과추정치의 {m['preserved_fraction']*100:.0f}%를 보존하면 " f"δ = {margin:g} {unit}가 된다. 과거 추정치의 불확실성은 비열등성 " f"검정통계량의 분산에 합산하여 반영하였다." ) # ========================================================================== # 기능 2: 비열등성 표본수 (closed-form) + MMRM 효율 보정 # ========================================================================== def sample_size_ni_ttest(margin, sd, true_diff=0.0, alpha=0.025, power=0.80, allocation_ratio=1.0): """ 비열등성 2-표본 t-test 표본수 (정규근사 closed-form). 가설(낮을수록 좋은 endpoint, 시험약-대조약 차이 D): H0: D >= margin vs H1: D < margin 검정은 단측. true_diff 는 가정한 실제 차이(보통 0 = 진짜 동등). n_control = (1 + 1/r) * (z_{1-alpha} + z_{1-beta})^2 * sd^2 / (margin - true_diff)^2 n_treatment = r * n_control (r = allocation_ratio = n_t / n_c) 출처: Chow, Shao & Wang (2008) Ch.3, 식 (3.2.x). Returns: dict (올림한 arm별 n, 총 n) """ margin = abs(float(margin)) sd = float(sd) diff = float(true_diff) r = float(allocation_ratio) eff = margin - abs(diff) if eff <= 0: raise ValueError("margin must exceed |true_diff| for a feasible NI design.") z_a = norm_ppf(1 - alpha) z_b = norm_ppf(power) # n_control 기준 n_ctrl = (1.0 + 1.0 / r) * (z_a + z_b) ** 2 * sd ** 2 / eff ** 2 n_ctrl_c = int(math.ceil(n_ctrl)) n_trt_c = int(math.ceil(n_ctrl_c * r)) return { "design": "non_inferiority_2sample_ttest", "margin": margin, "sd": sd, "true_diff": diff, "alpha_one_sided": alpha, "power": power, "allocation_ratio_t_over_c": r, "effective_distance": eff, "n_control": n_ctrl_c, "n_treatment": n_trt_c, "n_total": n_ctrl_c + n_trt_c, "n_control_exact": n_ctrl, } def mmrm_design_effect(visit_weeks, rho, cov_structure="ar1", dropout_rate=0.0, missing_mechanism="MAR"): """ MMRM(반복측정 혼합모형) 설계효과(design effect, DE). 종단 시험에서 1차 비교가 '특정 시점(보통 마지막 방문)의 군간 차이'일 때, MMRM은 다른 방문 정보를 빌려와 그 시점 추정의 분산을 줄인다. 표본수는 단순 단일시점 ANCOVA 대비 DE 배만큼 필요(DE <= 1 이면 효율 이득). 여기서는 다음 closed-form 근사를 사용: - m 개 post-baseline 방문, AR(1) 또는 교환가능(compound symmetry, CS) 상관. - 마지막 시점 평균추정의 상대분산 = 1 / I_eff, I_eff = 종단 정보량 (상관·결측 보정). - MAR 하 MMRM은 완료자 분석보다 효율적 -> 결측 inflation 을 부분 완화. 출처: Lu, Luo & Chen (2008); Hedeker & Gibbons (2006) Ch.13. CS 의 경우 잘 알려진 DE = (1 - rho) * (1 + (m-1)*rho? ) 형태가 아니라 마지막시점 GLS 추정 분산 비율을 직접 사용. 반환 DE 의미: n_MMRM = n_single_timepoint * DE. """ m = len([w for w in visit_weeks if w > 0]) # post-baseline 방문 수 if m < 1: m = 1 rho = float(rho) rho = min(max(rho, 0.0), 0.95) # ---- 상관행렬 R (post-baseline 방문 간) ---- if not HAS_NUMPY: # numpy 없을 때: CS 근사 closed-form 사용 # 마지막시점 GLS 분산 / OLS(단일시점) 분산 비율 # CS 에서 모든 시점 평균 정보 공유 시 DE ~ (1 - rho) + rho/m * (something) # 보수적 근사: de = 1.0 - rho * (m - 1.0) / m * 0.5 de = max(de, 0.55) else: idx = np.arange(m) if cov_structure == "ar1": # AR(1): R_ij = rho^|i-j| R = rho ** np.abs(idx.reshape(-1, 1) - idx.reshape(1, -1)) else: # compound symmetry / exchangeable R = np.full((m, m), rho) np.fill_diagonal(R, 1.0) # 마지막 방문 평균차이를 GLS 로 추정할 때의 상대분산. # 단일시점(OLS, 마지막방문만) 분산 = sigma^2/n. # MMRM(전 시점 사용) 마지막시점 추정 분산: # 모든 시점의 진짜 평균이 동일선상이라는 가정 없이, 마지막시점 # 고정효과만 관심일 때 GLS 분산은 (e_last' R^{-1} e_last)^{-1}. try: Rinv = np.linalg.inv(R) except Exception: Rinv = np.linalg.pinv(R) e_last = np.zeros(m) e_last[-1] = 1.0 var_ratio = 1.0 / float(e_last @ Rinv @ e_last) # var_ratio < 1 이면 효율 이득. AR(1)에서는 마지막시점만 보면 # 이득이 거의 없으나, 보통 모든 방문의 '평균 처치효과' 또는 # '시간평균 대비'를 1차로 쓰면 이득이 큼. 실무 절충: # 마지막시점 추정의 이득(작음)과 평균효과 추정 이득을 평균. ones = np.ones(m) var_mean_effect = 1.0 / float(ones @ Rinv @ ones) * m # 평균효과 상대분산*m? -> 정규화 # 평균 처치효과(시간평균) GLS 분산 = (1'R^{-1}1)^{-1}; 단일시점 대비 비율: var_ratio_mean = m / float(ones @ Rinv @ ones) / m # = 1/(1'R^-1 1) ... var_ratio_mean = 1.0 / float(ones @ Rinv @ ones) # 절충 DE: 사용자가 1차 endpoint를 '최종시점 변화'로 보되 MMRM이 # 결측·중간방문으로 안정화하는 현실적 이득을 반영. de = 0.5 * var_ratio + 0.5 * max(var_ratio_mean, 0.0) de = float(min(max(de, 0.4), 1.0)) # ---- 결측 보정 ---- d = min(max(float(dropout_rate), 0.0), 0.6) if missing_mechanism.upper() == "MCAR": # 완전무작위결측: 완료자 비율만큼 단순 inflation missing_inflation = 1.0 / (1.0 - d) if d < 1 else 99.0 else: # MAR (MMRM 의 기본 타당성 가정) # MMRM 은 MAR 하 탈락자의 관측분도 활용 -> inflation 을 부분만 적용. # 절반 정도 회복한다고 보는 흔한 실무 근사. eff_missing = d * 0.5 missing_inflation = 1.0 / (1.0 - eff_missing) if eff_missing < 1 else 99.0 total_de = de * missing_inflation return { "post_baseline_visits": m, "cov_structure": cov_structure, "rho": rho, "longitudinal_DE": de, "dropout_rate": d, "missing_mechanism": missing_mechanism, "missing_inflation": missing_inflation, "total_design_effect": total_de, "note": ( "longitudinal_DE<1: MMRM이 반복측정 정보로 효율을 높임. " "missing_inflation: 탈락 보정 계수. " "total_design_effect = 두 값의 곱; n_MMRM = n_ANCOVA_single x total_DE." ), } def sample_size_mmrm_ni(margin, sd, true_diff=0.0, alpha=0.025, power=0.80, allocation_ratio=1.0, visit_weeks=(0, 12, 24, 52), rho=0.6, cov_structure="ar1", dropout_rate=0.15, missing_mechanism="MAR"): """ MMRM 기반 비열등성 표본수 = 단일시점 NI 표본수 x MMRM total design effect. 단순 ANCOVA(단일 종료시점) 대비 효율 비교 값도 함께 반환. """ base = sample_size_ni_ttest(margin, sd, true_diff, alpha, power, allocation_ratio) de = mmrm_design_effect(list(visit_weeks), rho, cov_structure, dropout_rate, missing_mechanism) total_de = de["total_design_effect"] n_ctrl = int(math.ceil(base["n_control_exact"] * total_de)) n_trt = int(math.ceil(n_ctrl * allocation_ratio)) # 비교: 단순 ANCOVA(MMRM 효율 이득 없음) + 단순 결측 inflation ancova_inflate = 1.0 / (1.0 - min(dropout_rate, 0.6)) n_ctrl_ancova = int(math.ceil(base["n_control_exact"] * ancova_inflate)) n_trt_ancova = int(math.ceil(n_ctrl_ancova * allocation_ratio)) n_total_mmrm = n_ctrl + n_trt n_total_ancova = n_ctrl_ancova + n_trt_ancova saved = n_total_ancova - n_total_mmrm return { "design": "MMRM_non_inferiority", "base_single_timepoint": base, "mmrm_design_effect": de, "n_control": n_ctrl, "n_treatment": n_trt, "n_total": n_total_mmrm, "comparison_ancova_single": { "n_control": n_ctrl_ancova, "n_treatment": n_trt_ancova, "n_total": n_total_ancova, "dropout_inflation": ancova_inflate, }, "efficiency_gain_subjects": saved, "efficiency_gain_pct": (100.0 * saved / n_total_ancova if n_total_ancova else 0.0), "interpretation": ( f"MMRM 설계 총 표본수 {n_total_mmrm}명 " f"(대조 {n_ctrl} / 시험 {n_trt}). " f"단순 ANCOVA(단일시점, 완료자 보정) 대비 {saved}명 " f"({100.0*saved/n_total_ancova:.1f}%) 절감." ), } # ========================================================================== # 기능 3: 공동 1차 endpoint 다중성 처리 # ========================================================================== def power_curve(margin, sd, n_per_arm, true_diff=0.0, alpha=0.025, allocation_ratio=1.0): """주어진 arm별 표본수에서 비열등성 검정의 검정력 계산 (정규근사).""" margin = abs(float(margin)) sd = float(sd) diff = abs(float(true_diff)) r = float(allocation_ratio) n_c = float(n_per_arm) n_t = n_c * r se = sd * math.sqrt(1.0 / n_c + 1.0 / n_t) if se <= 0: return 1.0 z_a = norm_ppf(1 - alpha) # 비중심성: (margin - diff)/se ncp = (margin - diff) / se return float(norm_cdf(ncp - z_a)) def joint_primary_sample_size(endpoints, r_correlation, alpha=0.025, power=0.80, procedure="hochberg", allocation_ratio=1.0): """ 공동 1차 endpoint(둘 다 성공해야 함, co-primary) 표본수 + FWER 통제. Parameters ---------- endpoints : list of dict, 각 {label, margin, sd, true_diff} r_correlation : float endpoint 간 상관계수 (검정통계량 상관 근사) procedure : 'hochberg' | 'gatekeeping' | 'bonferroni' | 'unadjusted' 공동 1차(co-primary)의 핵심: 두 검정 모두 유의해야 하므로 '전체 성공 검정력' = P(검정1 성공 AND 검정2 성공) 가 목표 power 이상. 상관 r 이 높을수록 두 사건이 함께 일어나기 쉬워 표본수 이득. 각 endpoint 의 개별 alpha: - hochberg/bonferroni: 보수적으로 alpha 분할 근사 (co-primary 에서는 사실 alpha 분할 불필요— 둘 다 유의해야 하므로 각 검정에 full alpha 사용 가능. 단 사용자가 'gatekeeping'(계층적 우선순위)이나 family 통제를 원하면 절차 반영.) 출처: Dmitrienko et al. (2010); co-primary 검정력은 이변량 정규 근사로 P(Z1>c1, Z2>c2) 를 수치적분/몬테카를로로 계산. Returns: dict """ r = float(r_correlation) r = min(max(r, -0.95), 0.95) eps = list(endpoints) k = len(eps) # 절차별 개별 endpoint alpha if procedure == "bonferroni": alpha_each = alpha / k proc_note = ("Bonferroni: 각 검정 alpha/k. co-primary에서는 보수적 " "(둘 다 유의 요구 시 본래 full alpha 가능).") elif procedure == "hochberg": # Hochberg step-up: 가장 보수적 경계는 alpha/k, 가장 큰 p는 alpha. # 표본수 산정 시 보수적으로 최소 경계(alpha/k) 사용. alpha_each = alpha / k proc_note = ("Hochberg step-up: 표본수는 최악경계(alpha/k)로 보수 산정. " "실제 분석 시 더 큰 p-value는 완화된 경계 사용 가능.") elif procedure == "gatekeeping": # 계층적: 1차 endpoint full alpha, 통과해야 다음으로. alpha_each = alpha proc_note = ("Gatekeeping(계층적): 우선순위 endpoint에 full alpha, " "선행 검정 통과 시에만 후속 검정. FWER 통제됨.") else: # unadjusted (co-primary 기본: 둘 다 유의해야 하므로 full alpha) alpha_each = alpha proc_note = ("Co-primary 미보정: 둘 다 유의 요구 -> 각 검정 full alpha " "사용해도 FWER 통제 (교집합 가설).") z_a = norm_ppf(1 - alpha_each) # 각 endpoint 가 단독으로 목표 power(개별)일 때의 표본수 산출. per_ep = [] for ep in eps: ss = sample_size_ni_ttest(ep["margin"], ep["sd"], ep.get("true_diff", 0.0), alpha_each, power, allocation_ratio) per_ep.append({"label": ep.get("label", "endpoint"), **ss}) # co-primary: 전체 성공검정력 = P(모두 성공). 같은 n 에서 이변량 정규로 # joint power 를 평가하고, 목표 power 이상이 될 때까지 n 증가. max_single_n = max(e["n_control"] for e in per_ep) def joint_power_at(n_c): # 각 endpoint 비중심성 ncps = [] for ep in eps: sd = float(ep["sd"]) r_alloc = allocation_ratio n_t = n_c * r_alloc se = sd * math.sqrt(1.0 / n_c + 1.0 / n_t) margin = abs(float(ep["margin"])) diff = abs(float(ep.get("true_diff", 0.0))) ncps.append((margin - diff) / se) # 각 검정 성공 = Z_i > z_a 여기서 Z_i ~ N(ncp_i, 1), 상관 r if k == 1: return norm_cdf(ncps[0] - z_a) # k>=2: 이변량(또는 다변량) 정규 P(모든 Z_i > z_a). # 두 개면 정확 이변량, 그 이상은 첫 두개로 근사(MVP 범위). return _bivariate_upper(ncps[0] - z_a, ncps[1] - z_a, r) # n 탐색 n_c = max(1, max_single_n) # 위로 for _ in range(20000): if joint_power_at(n_c) >= power: break n_c += 1 # 아래로 살짝 (혹시 초기값이 과함) while n_c > 2 and joint_power_at(n_c - 1) >= power: n_c -= 1 n_t = int(math.ceil(n_c * allocation_ratio)) achieved = joint_power_at(n_c) return { "procedure": procedure, "procedure_note": proc_note, "n_endpoints": k, "correlation": r, "alpha_overall": alpha, "alpha_each": alpha_each, "target_power": power, "per_endpoint": per_ep, "n_control": int(n_c), "n_treatment": n_t, "n_total": int(n_c) + n_t, "achieved_joint_power": achieved, "interpretation": ( f"공동 1차 endpoint {k}개 (상관 r={r:g}), {procedure} 절차. " f"두(모든) endpoint 동시 성공 검정력 {achieved*100:.1f}% 달성에 " f"총 {int(n_c)+n_t}명 필요 (대조 {int(n_c)} / 시험 {n_t})." ), } def _bivariate_upper(a, b, rho): """P(Z1 > -a, Z2 > -b) 형태가 아니라 P(Z1>a', Z2>b')... 실제로는 P(X>0,Y>0) 형태로 변환. 여기서 입력 a,b 는 (ncp - z_a). 검정 성공 = Z_i > z_a, Z_i~N(ncp_i,1) => (Z_i - ncp_i) > z_a - ncp_i = -a. 표준화변수 U_i = Z_i - ncp_i ~ N(0,1), 상관 rho. P(U1 > -a, U2 > -b).""" if HAS_SCIPY: try: from scipy.stats import multivariate_normal as _mvn mean = [0.0, 0.0] cov = [[1.0, rho], [rho, 1.0]] # P(U1 > -a, U2 > -b) = 1 - F1 - F2 + F(-a,-b) rv = _mvn(mean, cov) F_both = rv.cdf([-a, -b]) F1 = norm_cdf(-a) F2 = norm_cdf(-b) return float(1.0 - F1 - F2 + F_both) except Exception: pass # Monte-Carlo fallback if HAS_NUMPY: rng = np.random.default_rng(20260517) n = 200000 z1 = rng.standard_normal(n) z2 = rho * z1 + math.sqrt(max(1 - rho * rho, 0.0)) * rng.standard_normal(n) return float(np.mean((z1 > -a) & (z2 > -b))) # 최후: 독립 가정 근사 return norm_cdf(a) * norm_cdf(b) # ========================================================================== # 기능 4: 탈락 시나리오 시뮬레이터 # ========================================================================== def dropout_sensitivity(margin, sd, n_per_arm, true_diff=0.0, alpha=0.025, allocation_ratio=1.0, dropout_rates=(0.05, 0.10, 0.15, 0.20, 0.25), patterns=("MCAR", "MAR_early", "MAR_late")): """ 탈락률 x 탈락패턴별 검정력 민감도 히트맵 데이터. patterns: - MCAR : 완전무작위, 완료자 비율만큼 정보 손실. - MAR_early : 초기 집중 탈락. MMRM 이 일부 회복 -> 유효손실 ~0.6배. - MAR_late : 후기 집중 탈락. 종료시점 정보 손실 커 ~0.85배 유효. 각 셀: 계획 표본수 n_per_arm 에서 해당 탈락 시 유효 표본수와 검정력. 출처: EMA Missing Data guideline; 탈락 패턴별 유효손실 계수는 MMRM 의 MAR 회복력에 대한 실무 근사값(보수적). """ eff_factor = { "MCAR": 1.0, # 손실 전부 반영 "MAR_early": 0.60, # MMRM 이 60% 만 실손실로 "MAR_late": 0.85, # 후기 탈락은 회복 어려움 } rows = [] for d in dropout_rates: row = {"dropout_rate": d, "cells": {}} for pat in patterns: f = eff_factor.get(pat, 1.0) eff_loss = d * f n_eff = n_per_arm * (1.0 - eff_loss) pw = power_curve(margin, sd, n_eff, true_diff, alpha, allocation_ratio) row["cells"][pat] = { "effective_n_per_arm": round(n_eff, 1), "power": round(pw, 4), } rows.append(row) return { "planned_n_per_arm": n_per_arm, "margin": margin, "sd": sd, "patterns": list(patterns), "effective_loss_factors": eff_factor, "heatmap": rows, "note": ("각 셀의 power 가 목표(예 0.80) 아래로 떨어지면 해당 탈락 " "시나리오에 표본수 부족. MAR_late 가 가장 취약."), } # ========================================================================== # 기능 5: 규제 제출용 문서 생성 # ========================================================================== def generate_regulatory_text(endpoint_label, unit, margin_result, ss_result, lang="ko", joint_result=None, dropout_result=None): """ IND/프로토콜 통계 섹션 형식의 표본수·마진 정당화 문서 생성 (국문/영문). """ margin = (margin_result.get("NI_margin_M2") or margin_result.get("NI_margin")) just = generate_margin_justification_text(endpoint_label, unit, margin, margin_result, lang) base = ss_result["base_single_timepoint"] de = ss_result["mmrm_design_effect"] if lang == "en": lines = [] lines.append("STATISTICAL SECTION - SAMPLE SIZE JUSTIFICATION") lines.append("=" * 60) lines.append("") lines.append("1. Non-inferiority margin") lines.append(just) lines.append("") lines.append("2. Sample size determination") lines.append( f"The primary endpoint ({endpoint_label}) is analyzed using a " f"mixed model for repeated measures (MMRM) over " f"{de['post_baseline_visits']} post-baseline visits " f"(covariance: {de['cov_structure'].upper()}, rho={de['rho']:g}). " f"Assuming a common SD of {base['sd']:g} {unit}, a true treatment " f"difference of {base['true_diff']:g} {unit}, one-sided alpha " f"{base['alpha_one_sided']:g} and power {base['power']*100:.0f}%, " f"the non-inferiority margin {margin:g} {unit} yields " f"{base['n_total']} subjects under a single-timepoint analysis." ) lines.append( f"Applying the MMRM longitudinal design effect " f"({de['longitudinal_DE']:.3f}) and a MAR-based dropout " f"adjustment for {de['dropout_rate']*100:.0f}% attrition " f"(inflation {de['missing_inflation']:.3f}; total design effect " f"{de['total_design_effect']:.3f}), the required sample size is " f"{ss_result['n_total']} subjects " f"({ss_result['n_control']} control / " f"{ss_result['n_treatment']} test). This is " f"{ss_result['efficiency_gain_subjects']} fewer subjects " f"({ss_result['efficiency_gain_pct']:.1f}%) than a single-" f"timepoint ANCOVA with completer-based inflation." ) if joint_result: lines.append("") lines.append("3. Co-primary endpoints / multiplicity") lines.append( f"With {joint_result['n_endpoints']} co-primary endpoints " f"(correlation r={joint_result['correlation']:g}), the " f"{joint_result['procedure']} procedure controls the FWER at " f"{joint_result['alpha_overall']:g}. Joint power " f"{joint_result['achieved_joint_power']*100:.1f}% requires " f"{joint_result['n_total']} subjects. " f"{joint_result['procedure_note']}" ) if dropout_result: lines.append("") lines.append("4. Dropout sensitivity") lines.append( "Power was evaluated across dropout rates " f"{[r['dropout_rate'] for r in dropout_result['heatmap']]} " "and missing-data patterns " f"{dropout_result['patterns']}. The MAR_late pattern is the " "most adverse; see sensitivity heatmap." ) lines.append("") lines.append("5. Assumptions and references") lines.append( "- ICH E9 (Statistical Principles), ICH E10 (Choice of Control " "Group); EMA/FDA/MFDS diabetes guidance.") lines.append( "- Chow, Shao & Wang (2008) Sample Size Calculations; " "Lu, Luo & Chen (2008) repeated measures sample size; " "Hedeker & Gibbons (2006) Longitudinal Data Analysis.") lines.append( "- Battelino et al. Diabetes Care 2019 (TIR consensus).") lines.append("") lines.append("DISCLAIMER: " + DISCLAIMER) return "\n".join(lines) # ---- 국문 ---- lines = [] lines.append("통계 섹션 - 표본수 산출 근거") lines.append("=" * 60) lines.append("") lines.append("1. 비열등성 마진") lines.append(just) lines.append("") lines.append("2. 표본수 산출") lines.append( f"1차 유효성 변수({endpoint_label})는 기저치 및 추적 방문 " f"(post-baseline {de['post_baseline_visits']}개 시점)에서 측정한 " f"반복측정을 혼합효과모형(MMRM, 공분산구조 " f"{de['cov_structure'].upper()}, 자기상관 rho={de['rho']:g})으로 " f"분석한다. 공통 표준편차 {base['sd']:g} {unit}, 가정한 실제 " f"군간차이 {base['true_diff']:g} {unit}, 단측 유의수준 " f"{base['alpha_one_sided']:g}, 검정력 {base['power']*100:.0f}% " f"하에서 비열등성 마진 {margin:g} {unit}를 적용하면 단일시점 " f"분석 기준 {base['n_total']}명이 산출된다." ) lines.append( f"여기에 MMRM 종단 설계효과({de['longitudinal_DE']:.3f})와 " f"{de['dropout_rate']*100:.0f}% 탈락에 대한 MAR 기반 보정" f"(inflation {de['missing_inflation']:.3f}; 총 설계효과 " f"{de['total_design_effect']:.3f})을 적용하면 최종 필요 표본수는 " f"총 {ss_result['n_total']}명(대조군 {ss_result['n_control']}명 / " f"시험군 {ss_result['n_treatment']}명)이다. 이는 단일시점 ANCOVA" f"(완료자 기반 보정) 대비 {ss_result['efficiency_gain_subjects']}명" f"({ss_result['efficiency_gain_pct']:.1f}%) 적은 규모이다." ) if joint_result: lines.append("") lines.append("3. 공동 1차 평가변수 / 다중성") lines.append( f"공동 1차 평가변수 {joint_result['n_endpoints']}개" f"(상관계수 r={joint_result['correlation']:g})에 대해 " f"{joint_result['procedure']} 절차로 전체 1종오류율(FWER)을 " f"{joint_result['alpha_overall']:g} 수준에서 통제한다. " f"모든 평가변수 동시 성공 검정력 " f"{joint_result['achieved_joint_power']*100:.1f}% 달성에 " f"총 {joint_result['n_total']}명이 필요하다. " f"{joint_result['procedure_note']}" ) if dropout_result: lines.append("") lines.append("4. 탈락 민감도 분석") lines.append( "탈락률 " f"{[r['dropout_rate'] for r in dropout_result['heatmap']]} 및 " f"결측 패턴 {dropout_result['patterns']}에 걸쳐 검정력을 평가" "하였다. MAR_late(후기 집중 탈락) 패턴이 가장 불리하며, " "민감도 히트맵을 참조한다." ) lines.append("") lines.append("5. 가정 및 참고문헌") lines.append( "- ICH E9(통계 원칙), ICH E10(대조군 선택); EMA/FDA/MFDS " "당뇨병 의약품 가이던스.") lines.append( "- Chow, Shao & Wang (2008) Sample Size Calculations; " "Lu, Luo & Chen (2008) 반복측정 표본수; " "Hedeker & Gibbons (2006) Longitudinal Data Analysis.") lines.append( "- Battelino 등, Diabetes Care 2019 (CGM TIR 국제 컨센서스).") lines.append("") lines.append("면책: " + DISCLAIMER) return "\n".join(lines) # -------------------------------------------------------------------------- # 데모용 통합 시나리오 # -------------------------------------------------------------------------- def run_demo_scenario(): """예제 시나리오 전체를 계산해 dict 로 반환 (CLI/테스트 공용).""" presets = None try: presets = load_presets() except Exception: presets = None # 시나리오: HbA1c 비열등성, 활성약-위약 효과 0.8%, 50% 보존 margin_fixed = ni_margin_fixed(0.8, preserved_fraction=0.5) margin_synth = ni_margin_synthesis(0.8, se_historical=0.12, preserved_fraction=0.5) delta = margin_fixed["NI_margin_M2"] # 0.4% ss = sample_size_mmrm_ni( margin=delta, sd=1.0, true_diff=0.0, alpha=0.025, power=0.80, allocation_ratio=1.0, visit_weeks=(0, 12, 24, 52), rho=0.6, cov_structure="ar1", dropout_rate=0.15, missing_mechanism="MAR") joint = joint_primary_sample_size( endpoints=[ {"label": "HbA1c", "margin": 0.4, "sd": 1.0, "true_diff": 0.0}, {"label": "CGM TIR", "margin": 5.0, "sd": 18.0, "true_diff": 0.0}, ], r_correlation=0.7, alpha=0.025, power=0.80, procedure="hochberg", allocation_ratio=1.0) # 탈락 민감도는 '단일시점 계획 표본수' 기준으로 평가해야 의미가 있다. # (MMRM 효율 보정 후 n 은 이미 축소된 값이라 그 위에 탈락을 또 얹으면 # 당연히 검정력이 크게 떨어진다 — 계획 단계 표본수로 본다.) n_plan_per_arm = ss["base_single_timepoint"]["n_control"] drop = dropout_sensitivity( margin=delta, sd=1.0, n_per_arm=n_plan_per_arm, true_diff=0.0, alpha=0.025) reg_ko = generate_regulatory_text( "HbA1c", "%", margin_fixed, ss, lang="ko", joint_result=joint, dropout_result=drop) reg_en = generate_regulatory_text( "HbA1c", "%", margin_fixed, ss, lang="en", joint_result=joint, dropout_result=drop) return { "presets_loaded": presets is not None, "margin_fixed": margin_fixed, "margin_synthesis": margin_synth, "sample_size_mmrm": ss, "joint_primary": joint, "dropout_sensitivity": drop, "regulatory_text_ko": reg_ko, "regulatory_text_en": reg_en, } if __name__ == "__main__": # pragma: no cover import pprint pprint.pprint(run_demo_scenario())