# -*- coding: utf-8 -*- """ ObesityRecruitFeas-Kor — 핵심 로직 모듈 ======================================= 항비만 임상연구 모집 feasibility 계산기의 순수 계산 로직. pandas / numpy 만 사용하며 Streamlit·외부 네트워크에 의존하지 않는다. 도메인 : Obesity (비만) 카테고리 : 인체실험 도구 (RCT/관찰연구 설계·계산기) 주요 함수 --------- load_microdata(path) 합성 KNHANES microdata 로드 build_eligibility_mask(df, crit) 적격성 기준 → boolean mask estimate_population(df, crit) 가중 외삽 적격 모집단 + 95% CI sensitivity_analysis(df, crit) 기준 민감도 (tornado 데이터) screen_fail_analysis(df, crit) 제외기준 적용 전후 + screen-fail율 recruitment_timeline(...) Monte-Carlo 모집 timeline generate_report(...) feasibility 리포트(국문 텍스트) 면책: 본 도구는 참고용·연구용이다. 실제 임상시험 feasibility는 CRO/생물통계 전문가 검토가 필수다. """ import os import math import numpy as np import pandas as pd # -------------------------------------------------------------------------- # 면책 문구 (앱·CLI·리포트 공통 사용) # -------------------------------------------------------------------------- DISCLAIMER = ( "[면책] 본 도구는 참고용·연구용 feasibility 설계 보조 도구입니다. " "합성 데이터(실제 KNHANES 아님) 기반 추정치이며, 실제 임상시험 " "feasibility 산정은 CRO 및 생물통계 전문가의 검토가 필수입니다." ) # 적격성 기준 기본값 — 항비만제 RCT 전형적 설계 예시 DEFAULT_CRITERIA = { # --- 포함(inclusion) 기준 --- "bmi_min": 27.0, # BMI 하한 (>=) "bmi_max": 45.0, # BMI 상한 (<=) "age_min": 19, # 연령 하한 (>=) "age_max": 70, # 연령 상한 (<=) # 동반질환 요구: 아래 목록 중 'comorbidity_required_count' 개 이상 보유 "comorbidities": [], # 후보 ["t2dm","htn","dyslipidemia","masld","osa"] "comorbidity_required_count": 0, "waist_min": 0.0, # 허리둘레 하한(cm). 0 이면 미적용 # --- 제외(exclusion) 기준 --- "exclude_t2dm": False, # T2DM 보유자 제외 "exclude_recent_glp1ra": True, # 최근 GLP-1RA 사용자 제외 "exclude_current_smoker": False, # 현재 흡연자 제외 "exclude_severe_htn": True, # 중증 고혈압(SBP>=160 또는 DBP>=100) 제외 "exclude_high_alt": True, # ALT > 3xULN (>120 U/L) 제외 } # 동반질환 코드 → microdata 컬럼 / 라벨 _COMORBIDITY_MAP = { "t2dm": ("t2dm", "제2형 당뇨병"), "htn": ("htn", "고혈압"), "dyslipidemia": ("dyslipidemia", "이상지질혈증"), "masld": ("masld_surrogate", "MASLD(대용지표)"), "osa": ("osa_high_risk", "OSA 고위험"), } # -------------------------------------------------------------------------- # 데이터 로드 # -------------------------------------------------------------------------- def default_data_path(): """패키지에 동봉된 합성 microdata CSV의 기본 경로.""" return os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "knhanes_synthetic.csv") def load_microdata(path=None): """합성 KNHANES microdata CSV를 DataFrame으로 로드한다. 파일이 없으면 안내 메시지와 함께 FileNotFoundError를 발생시킨다. """ path = path or default_data_path() if not os.path.exists(path): raise FileNotFoundError( f"합성 microdata를 찾을 수 없습니다: {path}\n" "먼저 data/generate_knhanes_synthetic.py 를 실행해 CSV를 생성하세요." ) df = pd.read_csv(path) required = {"age", "sex", "bmi", "weight"} missing = required - set(df.columns) if missing: raise ValueError(f"microdata에 필수 컬럼 누락: {sorted(missing)}") return df def merge_criteria(criteria=None): """사용자 입력 기준을 기본값과 병합해 완전한 기준 dict를 반환한다.""" merged = dict(DEFAULT_CRITERIA) if criteria: merged.update({k: v for k, v in criteria.items() if v is not None}) return merged # -------------------------------------------------------------------------- # 적격성 마스크 # -------------------------------------------------------------------------- def _inclusion_mask(df, crit): """포함 기준만 적용한 boolean Series를 반환한다.""" m = pd.Series(True, index=df.index) m &= df["bmi"] >= crit["bmi_min"] m &= df["bmi"] <= crit["bmi_max"] m &= df["age"] >= crit["age_min"] m &= df["age"] <= crit["age_max"] if crit.get("waist_min", 0) and crit["waist_min"] > 0: m &= df["waist_cm"] >= crit["waist_min"] # 동반질환 요구: 지정 목록 중 N개 이상 보유 comorbs = crit.get("comorbidities") or [] need = int(crit.get("comorbidity_required_count", 0) or 0) if comorbs and need > 0: count = pd.Series(0, index=df.index) for code in comorbs: col = _COMORBIDITY_MAP.get(code, (None,))[0] if col and col in df.columns: count = count + (df[col] == 1).astype(int) m &= count >= need return m def _exclusion_mask(df, crit): """제외 기준에 *해당되는*(즉 제외 대상) boolean Series를 반환한다.""" excl = pd.Series(False, index=df.index) if crit.get("exclude_t2dm"): excl |= df["t2dm"] == 1 if crit.get("exclude_recent_glp1ra"): excl |= df["recent_glp1ra"] == 1 if crit.get("exclude_current_smoker"): excl |= df["current_smoker"] == 1 if crit.get("exclude_severe_htn"): excl |= (df["sbp_mmhg"] >= 160) | (df["dbp_mmhg"] >= 100) if crit.get("exclude_high_alt"): excl |= df["alt_ul"] > 120 return excl def build_eligibility_mask(df, crit): """포함 충족 AND 제외 미해당 인 최종 적격 boolean Series.""" inc = _inclusion_mask(df, crit) excl = _exclusion_mask(df, crit) return inc & (~excl) # -------------------------------------------------------------------------- # 가중 모집단 외삽 + 95% 신뢰구간 # -------------------------------------------------------------------------- def _weighted_total(df, mask): """mask가 참인 행의 가중치 합.""" return float(df.loc[mask, "weight"].sum()) def estimate_population(df, crit): """적격성 기준을 적용해 한국 성인 적격 모집단 규모와 95% CI를 추정한다. 복합표본 분산을 단순 모사하기 위해, 가중 비율 p 의 표준오차를 이항 근사 + 설계효과(deff) 보정으로 계산한다. 반환 dict: n_sample_eligible : 적격 표본 수(비가중) weighted_eligible : 가중 외삽 적격 모집단 weighted_total : 가중 전체 성인 인구 prevalence : 적격 비율(0~1) ci_low / ci_high : 적격 모집단 95% CI prevalence_ci : 적격 비율 95% CI (low, high) """ crit = merge_criteria(crit) mask = build_eligibility_mask(df, crit) w = df["weight"] w_total = float(w.sum()) w_elig = _weighted_total(df, mask) n_sample = int(mask.sum()) n_total = len(df) p = w_elig / w_total if w_total > 0 else 0.0 # 설계효과(deff): 가중치 변동 기반 Kish 근사 deff = (n_total * (w ** 2).sum()) / (w.sum() ** 2) if w.sum() > 0 else 1.0 n_eff = n_total / deff if deff > 0 else n_total # 비율 p 의 표준오차 (유효표본 기준) if 0 < p < 1 and n_eff > 1: se = math.sqrt(p * (1 - p) / n_eff) else: se = 0.0 z = 1.96 p_low = max(0.0, p - z * se) p_high = min(1.0, p + z * se) return { "n_sample_eligible": n_sample, "n_sample_total": n_total, "weighted_eligible": w_elig, "weighted_total": w_total, "prevalence": p, "prevalence_ci": (p_low, p_high), "ci_low": p_low * w_total, "ci_high": p_high * w_total, "design_effect": deff, "effective_n": n_eff, } # -------------------------------------------------------------------------- # 기준 민감도 분석 (tornado) # -------------------------------------------------------------------------- def sensitivity_analysis(df, crit): """각 기준을 단독으로 완화했을 때 적격 모집단 증가량을 계산한다. 반환: 리스트(dict) — 항목별 label : 기준 설명 baseline : 현재 기준의 적격 모집단 relaxed : 해당 기준만 완화 시 적격 모집단 delta : 증가량(절대) delta_pct : 증가율(%) delta 내림차순 정렬 → tornado 차트 입력으로 사용. """ crit = merge_criteria(crit) base = estimate_population(df, crit)["weighted_eligible"] scenarios = [] def add(label, modified): m = merge_criteria({**crit, **modified}) relaxed = estimate_population(df, m)["weighted_eligible"] delta = relaxed - base scenarios.append({ "label": label, "baseline": base, "relaxed": relaxed, "delta": delta, "delta_pct": (100.0 * delta / base) if base > 0 else 0.0, }) # BMI 하한 1단위 완화 add(f"BMI 하한 {crit['bmi_min']:.0f}→{crit['bmi_min']-1:.0f}", {"bmi_min": crit["bmi_min"] - 1}) # BMI 상한 1단위 완화 add(f"BMI 상한 {crit['bmi_max']:.0f}→{crit['bmi_max']+1:.0f}", {"bmi_max": crit["bmi_max"] + 1}) # 연령 상한 5세 확대 add(f"연령 상한 {crit['age_max']}→{crit['age_max']+5}", {"age_max": crit["age_max"] + 5}) # 연령 하한 5세 확대 add(f"연령 하한 {crit['age_min']}→{max(19, crit['age_min']-5)}", {"age_min": max(19, crit["age_min"] - 5)}) # 동반질환 요구 개수 1 감소 need = int(crit.get("comorbidity_required_count", 0) or 0) if need > 0: add(f"동반질환 요구 {need}→{need-1}개", {"comorbidity_required_count": need - 1}) # 허리둘레 하한 해제 if crit.get("waist_min", 0) and crit["waist_min"] > 0: add(f"허리둘레 하한({crit['waist_min']:.0f}cm) 해제", {"waist_min": 0.0}) # 제외 기준 단독 해제 for key, label in [ ("exclude_t2dm", "제외기준 해제: T2DM"), ("exclude_recent_glp1ra", "제외기준 해제: 최근 GLP-1RA"), ("exclude_current_smoker", "제외기준 해제: 현재 흡연"), ("exclude_severe_htn", "제외기준 해제: 중증 고혈압"), ("exclude_high_alt", "제외기준 해제: ALT>3xULN"), ]: if crit.get(key): add(label, {key: False}) scenarios.sort(key=lambda s: s["delta"], reverse=True) return scenarios # -------------------------------------------------------------------------- # screen-fail 분석 # -------------------------------------------------------------------------- def screen_fail_analysis(df, crit): """제외 기준 적용 전후 적격자 차이로 예상 screen-fail율을 계산한다. screen-fail 정의(본 도구): 포함 기준은 충족했으나 제외 기준에 걸려 등록에 실패하는 비율. (운영상 발생하는 동의철회·기타 사유는 별도 attrition 파라미터로 timeline에서 다룬다.) 반환 dict: weighted_inclusion : 포함기준만 통과한 가중 인구 weighted_eligible : 포함 AND 제외통과 가중 인구 screen_fail_rate : (inclusion - eligible)/inclusion per_exclusion : 제외기준별 단독 영향(가중 인구) """ crit = merge_criteria(crit) inc = _inclusion_mask(df, crit) elig = build_eligibility_mask(df, crit) w_inc = _weighted_total(df, inc) w_elig = _weighted_total(df, elig) sf_rate = (w_inc - w_elig) / w_inc if w_inc > 0 else 0.0 # 제외기준별 단독 기여 (포함기준 통과자 중 해당 제외기준에 걸리는 가중 비율) per_exclusion = [] excl_specs = [ ("exclude_t2dm", "T2DM 보유", lambda d: d["t2dm"] == 1), ("exclude_recent_glp1ra", "최근 GLP-1RA 사용", lambda d: d["recent_glp1ra"] == 1), ("exclude_current_smoker", "현재 흡연", lambda d: d["current_smoker"] == 1), ("exclude_severe_htn", "중증 고혈압(SBP>=160/DBP>=100)", lambda d: (d["sbp_mmhg"] >= 160) | (d["dbp_mmhg"] >= 100)), ("exclude_high_alt", "ALT>3xULN(>120 U/L)", lambda d: d["alt_ul"] > 120), ] for key, label, fn in excl_specs: if crit.get(key): hit = inc & fn(df) w_hit = _weighted_total(df, hit) per_exclusion.append({ "label": label, "weighted_excluded": w_hit, "pct_of_inclusion": (100.0 * w_hit / w_inc) if w_inc > 0 else 0.0, }) per_exclusion.sort(key=lambda x: x["weighted_excluded"], reverse=True) return { "weighted_inclusion": w_inc, "weighted_eligible": w_elig, "screen_fail_rate": sf_rate, "per_exclusion": per_exclusion, } def screened_to_enrolled(target_n, screen_fail_rate, extra_attrition=0.10): """목표 등록 N 달성에 필요한 스크리닝 환자 수를 환산한다. extra_attrition: 동의철회/기타 운영 손실 추가 가정(기본 10%). """ keep = (1.0 - screen_fail_rate) * (1.0 - extra_attrition) if keep <= 0: return float("inf") return target_n / keep # -------------------------------------------------------------------------- # 모집 timeline Monte-Carlo # -------------------------------------------------------------------------- def recruitment_timeline(target_n, n_sites, screening_per_site_month, screen_fail_rate, extra_attrition=0.10, site_activation_months=2.0, n_sim=2000, seed=20260517, screening_cv=0.35, max_months=120): """사이트 기반 모집 timeline을 Monte-Carlo로 시뮬레이션한다. 매개변수 --------- target_n 목표 등록 환자 수 n_sites 참여 사이트 수 screening_per_site_month 사이트당 월간 스크리닝 capacity(평균) screen_fail_rate screen-fail율(0~1) extra_attrition 동의철회 등 추가 손실(0~1) site_activation_months 사이트 가동(IRB·계약)까지 평균 소요 개월 n_sim Monte-Carlo 반복수 screening_cv 월간 스크리닝 수의 변동계수(불확실성) max_months 시뮬레이션 상한(개월) 반환 dict: p10 / p50 / p90 목표 N 도달 개월 분위수 mean 평균 개월 enroll_yield 스크리닝→등록 환산율 feasible_ratio max_months 내 목표 도달한 시뮬레이션 비율 monthly_expected_enroll 결정론적 월간 기대 등록 수 """ rng = np.random.default_rng(seed) yield_rate = (1.0 - screen_fail_rate) * (1.0 - extra_attrition) yield_rate = max(yield_rate, 1e-6) # 결정론적 월간 기대 등록 수 (전 사이트 가동 시) monthly_expected = n_sites * screening_per_site_month * yield_rate months_to_target = np.full(n_sim, float(max_months)) for s in range(n_sim): # 사이트별 가동 시점(개월): 평균 site_activation_months 의 지수분포 근사 if site_activation_months > 0: activation = rng.exponential(site_activation_months, size=n_sites) else: activation = np.zeros(n_sites) enrolled = 0.0 month = 0 while enrolled < target_n and month < max_months: month += 1 # 이번 달 가동된 사이트 수 active = int(np.sum(activation <= month)) if active <= 0: continue # 사이트별 월간 스크리닝 수: 감마분포(평균=capacity, CV=screening_cv) mean_s = screening_per_site_month if mean_s > 0 and screening_cv > 0: shape = 1.0 / (screening_cv ** 2) scale = mean_s / shape screened = rng.gamma(shape, scale, size=active).sum() else: screened = mean_s * active # 스크리닝 → 등록 (이항 변동) new_enroll = rng.binomial(int(round(screened)), min(yield_rate, 1.0)) enrolled += new_enroll months_to_target[s] = month if enrolled >= target_n else float(max_months) feasible = months_to_target < max_months return { "p10": float(np.percentile(months_to_target, 10)), "p50": float(np.percentile(months_to_target, 50)), "p90": float(np.percentile(months_to_target, 90)), "mean": float(np.mean(months_to_target)), "enroll_yield": yield_rate, "monthly_expected_enroll": monthly_expected, "feasible_ratio": float(np.mean(feasible)), "max_months": max_months, "n_sim": n_sim, } # -------------------------------------------------------------------------- # feasibility 리포트 (국문 텍스트) # -------------------------------------------------------------------------- def _fmt(n): """정수형 천단위 구분 포맷.""" try: return f"{n:,.0f}" except (ValueError, TypeError): return str(n) def criteria_summary_text(crit): """적격성 기준을 사람이 읽는 국문 요약으로 변환한다.""" crit = merge_criteria(crit) lines = [] lines.append(f" - BMI: {crit['bmi_min']:.0f} ~ {crit['bmi_max']:.0f} kg/m^2") lines.append(f" - 연령: {crit['age_min']} ~ {crit['age_max']}세") if crit.get("waist_min", 0) and crit["waist_min"] > 0: lines.append(f" - 허리둘레: >= {crit['waist_min']:.0f} cm") comorbs = crit.get("comorbidities") or [] need = int(crit.get("comorbidity_required_count", 0) or 0) if comorbs and need > 0: labels = [_COMORBIDITY_MAP[c][1] for c in comorbs if c in _COMORBIDITY_MAP] lines.append(f" - 동반질환: [{', '.join(labels)}] 중 {need}개 이상") else: lines.append(" - 동반질환 요구: 없음") excl = [] if crit.get("exclude_t2dm"): excl.append("T2DM 보유자") if crit.get("exclude_recent_glp1ra"): excl.append("최근 GLP-1RA 사용자") if crit.get("exclude_current_smoker"): excl.append("현재 흡연자") if crit.get("exclude_severe_htn"): excl.append("중증 고혈압") if crit.get("exclude_high_alt"): excl.append("ALT>3xULN") lines.append(f" - 제외기준: {', '.join(excl) if excl else '없음'}") return "\n".join(lines) def generate_report(df, crit, target_n=200, n_sites=10, screening_per_site_month=8.0, extra_attrition=0.10, site_activation_months=2.0): """적격 모집단·민감도·screen-fail·timeline을 종합한 국문 feasibility 리포트.""" crit = merge_criteria(crit) pop = estimate_population(df, crit) sens = sensitivity_analysis(df, crit) sf = screen_fail_analysis(df, crit) tl = recruitment_timeline( target_n=target_n, n_sites=n_sites, screening_per_site_month=screening_per_site_month, screen_fail_rate=sf["screen_fail_rate"], extra_attrition=extra_attrition, site_activation_months=site_activation_months, ) need_screen = screened_to_enrolled(target_n, sf["screen_fail_rate"], extra_attrition) L = [] L.append("=" * 66) L.append(" 항비만 임상연구 모집 FEASIBILITY 리포트") L.append(" (ObesityRecruitFeas-Kor)") L.append("=" * 66) L.append("") L.append("1. 적격성 기준") L.append(criteria_summary_text(crit)) L.append("") L.append("2. 적격 모집단 추정 (합성 KNHANES 가중 외삽)") L.append(f" - 적격 표본 수: {pop['n_sample_eligible']:,} / " f"{pop['n_sample_total']:,} 명") L.append(f" - 적격 비율: {pop['prevalence']*100:.2f}% " f"(95% CI {pop['prevalence_ci'][0]*100:.2f}~" f"{pop['prevalence_ci'][1]*100:.2f}%)") L.append(f" - 한국 성인 적격 모집단(외삽): 약 {_fmt(pop['weighted_eligible'])} 명") L.append(f" 95% CI: {_fmt(pop['ci_low'])} ~ {_fmt(pop['ci_high'])} 명") L.append(f" - 설계효과(deff) {pop['design_effect']:.2f}, " f"유효표본 {pop['effective_n']:.0f}") L.append("") L.append("3. 기준 민감도 분석 (단독 완화 시 모집단 증가 — 상위 5개)") if sens: for s in sens[:5]: L.append(f" - {s['label']}: +{_fmt(s['delta'])} 명 " f"({s['delta_pct']:+.1f}%)") top = sens[0] L.append(f" => 가장 '비싼' 기준: {top['label']} " f"(완화 시 +{_fmt(top['delta'])} 명)") else: L.append(" - (완화 가능한 기준 없음)") L.append("") L.append("4. screen-fail 예측") L.append(f" - 포함기준 통과 인구: 약 {_fmt(sf['weighted_inclusion'])} 명") L.append(f" - 최종 적격 인구: 약 {_fmt(sf['weighted_eligible'])} 명") L.append(f" - 예상 screen-fail율(제외기준 기준): " f"{sf['screen_fail_rate']*100:.1f}%") if sf["per_exclusion"]: L.append(" - 제외기준별 단독 영향(포함기준 통과자 대비):") for e in sf["per_exclusion"]: L.append(f" · {e['label']}: {e['pct_of_inclusion']:.1f}%") L.append(f" - 추가 운영 손실(동의철회 등) 가정: " f"{extra_attrition*100:.0f}%") L.append(f" - 목표 등록 {target_n}명 달성에 필요한 스크리닝 수: " f"약 {_fmt(need_screen)} 명") L.append("") L.append("5. 모집 timeline 시뮬레이션 (Monte-Carlo, " f"n_sim={tl['n_sim']})") L.append(f" - 설정: 사이트 {n_sites}개, 사이트당 월 스크리닝 " f"{screening_per_site_month:.0f}명, " f"사이트 가동 평균 {site_activation_months:.0f}개월") L.append(f" - 스크리닝→등록 환산율: {tl['enroll_yield']*100:.1f}%") L.append(f" - 전 사이트 가동 시 월 기대 등록: " f"{tl['monthly_expected_enroll']:.1f} 명") L.append(f" - 목표 N={target_n} 도달 기간: " f"P10 {tl['p10']:.0f} / P50 {tl['p50']:.0f} / " f"P90 {tl['p90']:.0f} 개월") L.append(f" - {tl['max_months']}개월 내 목표 도달 확률: " f"{tl['feasible_ratio']*100:.1f}%") L.append("") L.append("6. 종합 판단") L.append(_feasibility_verdict(pop, tl, target_n)) L.append("") L.append("-" * 66) L.append(DISCLAIMER) L.append("주: 적격 모집단은 합성 KNHANES microdata 기반 외삽치이며 " "실제 KNHANES 데이터가 아닙니다.") L.append("=" * 66) return "\n".join(L) def _feasibility_verdict(pop, tl, target_n): """모집단·timeline 기반 간이 종합 판정 문구.""" msgs = [] elig = pop["weighted_eligible"] # 모집단 충분성: 적격 모집단이 목표 N의 최소 1000배 이상이면 여유 ratio = elig / target_n if target_n > 0 else 0 if ratio >= 2000: msgs.append(" - 적격 모집단 규모는 목표 대비 충분합니다.") elif ratio >= 200: msgs.append(" - 적격 모집단은 확보 가능하나 사이트 접근성이 " "관건입니다.") else: msgs.append(" - 적격 모집단이 제한적입니다. 기준 완화 검토를 " "권고합니다.") # timeline 판정 if tl["feasible_ratio"] >= 0.9 and tl["p90"] <= 24: msgs.append(" - 통상적 임상시험 모집 기간(<=24개월) 내 달성이 " "유력합니다.") elif tl["feasible_ratio"] >= 0.7: msgs.append(" - 모집 기간이 다소 길어질 수 있어 사이트 추가 또는 " "기준 완화를 고려하십시오.") else: msgs.append(" - 현재 설정으로는 모집 지연 위험이 큽니다. " "사이트 수·capacity·기준 재검토가 필요합니다.") return "\n".join(msgs) # -------------------------------------------------------------------------- # 자기 점검용 (모듈 직접 실행 시) # -------------------------------------------------------------------------- if __name__ == "__main__": _df = load_microdata() print(generate_report(_df, DEFAULT_CRITERIA))