第 13 章 · ADaM 衍生与 TFL 报表生成

import decimal
from decimal import Decimal, ROUND_HALF_UP

def sas_round(x, unit=1):
    """复刻 SAS 的 ROUND(x, unit)：四舍五入（远离零）。

    sas_round(2.5, 1)    -> 3.0
    sas_round(2.675, 0.01) -> 2.68
    sas_round(1234, 10)  -> 1230
    """
    if x is None or (isinstance(x, float) and pd.isna(x)):
        return x
    unit = Decimal(str(unit))
    d = Decimal(str(x)) / unit
    # SAS 用"远离零"的舍入
    rounded = d.quantize(Decimal("1"), rounding=ROUND_HALF_UP)
    return float(rounded * unit)


def sas_round_series(s: pd.Series, decimals: int = 2) -> pd.Series:
    """整列按 SAS 规则舍入到 N 位小数。向量化实现（比 apply 快）。"""
    factor = 10 ** decimals
    return (np.sign(s) * np.floor(np.abs(s) * factor + 0.5) / factor)


# 验证
import numpy as np, pandas as pd
print(sas_round(2.5, 1))          # 3.0      （Python 内置 round(2.5) 给 2）
print(sas_round(2.675, 0.01))     # 2.68
print(sas_round(1234, 10))        # 1230.0
print(sas_round(-2.5, 1))         # -3.0     （SAS 也是 -3）

s = pd.Series([2.5, 1.5, 0.5, -2.5, 2.675])
print(sas_round_series(s, 0).tolist())   # [3.0, 2.0, 1.0, -3.0, 3.0]

"""clinic/report.py —— TFL 报表构造工具"""

import numpy as np
import pandas as pd


def pct_format(n: int, denom: int, decimals: int = 1,
               zero_as: str = "0", min_n: int = 0) -> str:
    """生成 "n (x.x%)" 格式的单元格。

    Parameters
    ----------
    n : 分子（受试者数）
    denom : 分母（该治疗组的安全性人群人数）
    decimals : 百分比小数位
    zero_as : 分子为 0 时显示什么（CSR 惯例通常显示 "0"）
    min_n : 只显示 ">0" 的单元格（用于"按频率降序只显示前 N 个"场景）
    """
    if denom == 0 or pd.isna(denom):
        return "-"
    if n == 0:
        return zero_as
    p = n / denom * 100
    # 用 SAS 规则舍入
    factor = 10 ** decimals
    p = np.sign(p) * np.floor(np.abs(p) * factor + 0.5) / factor
    return f"{int(n)} ({p:.{decimals}f}%)"


def summarize_continuous(s: pd.Series, decimals: int = 1) -> dict:
    """连续变量的标准描述统计（复刻 PROC MEANS 的常用输出）。

    返回 n / mean / sd / median / min / max，均按 SAS 规则舍入。
    """
    s = s.dropna()
    n = len(s)
    if n == 0:
        return {k: None for k in ("n", "mean", "sd", "median", "min", "max")}
    r = lambda v: float(np.sign(v) * np.floor(abs(v) * 10**decimals + 0.5) / 10**decimals)
    return {
        "n": n,
        "mean": r(s.mean()),
        "sd": r(s.std(ddof=1)) if n > 1 else None,   # ★ ddof=1 = SAS 的 STD
        "median": r(s.median()),
        "min": r(s.min()),
        "max": r(s.max()),
    }

# 演示这个坑
import numpy as np, pandas as pd
s = pd.Series([63, 64, 71, 74, 55])
s.std()          # 7.60    ← 与 SAS 的 STD 一致 ✅
np.std(s)        # 6.80    ← 与 SAS 不一致 ❌
np.std(s, ddof=1)  # 7.60  ← 必须显式指定

"""case02_人口学表.py —— 生成 CSR 标准的人口学与基线特征表（Table 1）"""
from pathlib import Path
import numpy as np
import pandas as pd

BASE = Path(__file__).resolve().parent.parent
SAMPLES = BASE / "data" / "samples"

# ★ 治疗组顺序必须显式声明（绝不能依赖字母序！）
TRT_LEVELS = ["Placebo", "Xanomeline Low Dose", "Xanomeline High Dose"]

adsl = pd.read_csv(SAMPLES / "adsl.csv")
adsl["TRT01P"] = pd.Categorical(adsl["TRT01P"], categories=TRT_LEVELS, ordered=True)

# 分析人群：安全性人群
adsl = adsl.query("SAFFL == 'Y'").copy()

# ---- 分母：各治疗组人数 ----
denom = adsl.groupby("TRT01P", observed=False)["USUBJID"].nunique()
# Placebo 86, Xanomeline Low Dose 84, Xanomeline High Dose 84  (+ 合计 254)


def cont_row(label, var, df=adsl, decimals=1):
    """连续变量行：Mean (SD) / Median / Min-Max"""
    out = {"变量": label, "类别": ""}
    for trt in TRT_LEVELS:
        s = df.loc[df["TRT01P"] == trt, var].dropna()
        if len(s) == 0:
            out[trt] = "-"
            continue
        mean = s.mean(); sd = s.std(ddof=1); med = s.median()
        r = lambda v, d=decimals: f"{v:.{d}f}"
        out[trt] = f"{r(mean)} ({r(sd)})"
    out["合计"] = _cont_total(df[var].dropna(), decimals)
    return out


def _cont_total(s, decimals=1):
    if len(s) == 0:
        return "-"
    return f"{s.mean():.{decimals}f} ({s.std(ddof=1):.{decimals}f})"


def cat_row(label, var, level, df=adsl, decimals=1, indent="  "):
    """分类变量行：n (%)"""
    out = {"变量": label, "类别": indent + str(level)}
    total_n = int((df[var] == level).sum())
    for trt in TRT_LEVELS:
        sub = df[df["TRT01P"] == trt]
        n = int((sub[var] == level).sum())
        out[trt] = _fmt_pct(n, int(denom[trt]), decimals)
    out["合计"] = _fmt_pct(total_n, int(denom.sum()), decimals)
    return out


def _fmt_pct(n, d, decimals):
    if d == 0:
        return "-"
    if n == 0:
        return "0"
    return f"{n} ({n / d * 100:.{decimals}f}%)"


def n_row(label, var, df=adsl):
    """计数行：n"""
    out = {"变量": label, "类别": ""}
    for trt in TRT_LEVELS:
        out[trt] = str(int(df.loc[df["TRT01P"] == trt, var].nunique()))
    out["合计"] = str(int(df[var].nunique()))
    return out


# ---- 组装表格 ----
rows = []
rows.append({"变量": "受试者数", "类别": "", **{t: str(int(denom[t])) for t in TRT_LEVELS},
             "合计": str(int(denom.sum()))})
rows.append(cont_row("年龄 (岁)", "AGE"))
rows.append({"变量": "", "类别": "", **{t: "" for t in TRT_LEVELS}, "合计": ""})

rows.append({"变量": "年龄分组, n (%)", "类别": "", **{t: "" for t in TRT_LEVELS}, "合计": ""})
for lvl in ["<65", "65-80", ">80"]:
    rows.append(cat_row("", "AGEGR1", lvl))

rows.append({"变量": "", "类别": "", **{t: "" for t in TRT_LEVELS}, "合计": ""})
rows.append({"变量": "性别, n (%)", "类别": "", **{t: "" for t in TRT_LEVELS}, "合计": ""})
for lvl, lab in [("F", "女"), ("M", "男")]:
    rows.append(cat_row("", "SEX", lvl))

rows.append({"变量": "", "类别": "", **{t: "" for t in TRT_LEVELS}, "合计": ""})
rows.append({"变量": "种族, n (%)", "类别": "", **{t: "" for t in TRT_LEVELS}, "合计": ""})
for lvl in adsl["RACE"].value_counts().index:
    rows.append(cat_row("", "RACE", lvl))

table = pd.DataFrame(rows)

# ---- 输出 ----
pd.set_option("display.width", 220)
pd.set_option("display.max_colwidth", 45)
print("表 1. 人口学与基线特征（安全性人群）")
print(table.to_string(index=False))

table.to_csv(BASE / "outputs" / "table01_demographics.csv", index=False, encoding="utf-8-sig")

"""case03_不良事件汇总表.py —— TEAE 按 SOC / PT 的受试者数与百分比"""

adae = pd.read_csv(SAMPLES / "adae.csv")
adsl = pd.read_csv(SAMPLES / "adsl.csv")

TRT_LEVELS = ["Placebo", "Xanomeline Low Dose", "Xanomeline High Dose"]

# ---- Step 1: 确定分母（安全性人群）----
saf = adsl.query("SAFFL == 'Y'")[["USUBJID", "TRT01P"]].copy()
saf["TRT01P"] = pd.Categorical(saf["TRT01P"], categories=TRT_LEVELS, ordered=True)
denom = saf.groupby("TRT01P", observed=False)["USUBJID"].nunique()

# ---- Step 2: 限定分析事件（TEAE）----
# TRTEMFL='Y' 表示"治疗中出现的不良事件"
teae = adae.query("TRTEMFL == 'Y'").copy()
# 只保留安全性人群中的记录
teae = teae.merge(saf[["USUBJID", "TRT01P"]], on="USUBJID", how="inner")
# 有没有治疗组信息？ADAE 自带 TRTA，这里用 ADSL 的 TRT01P 更严谨
print(f"TEAE 记录数: {len(teae)}，涉及受试者: {teae['USUBJID'].nunique()}")

# ---- Step 3: ★核心★ 受试者层级去重 ----
# 同一受试者同一 SOC 发生多次 → 只计 1 人
subj_soc = teae.drop_duplicates(["USUBJID", "TRT01P", "AEBODSYS"])
subj_pt = teae.drop_duplicates(["USUBJID", "TRT01P", "AEBODSYS", "AEDECOD"])

n_any = teae.drop_duplicates(["USUBJID", "TRT01P"]).groupby(
    "TRT01P", observed=False)["USUBJID"].nunique()      # 有任一 TEAE 的人数


def count_table(df_uniq, keys):
    """按 keys 分组统计各治疗组的受试者数，返回长表。"""
    return (df_uniq.groupby(keys + ["TRT01P"], observed=False)["USUBJID"]
                   .nunique().rename("N").reset_index())


def pivot_counts(df_uniq, keys):
    c = count_table(df_uniq, keys)
    w = c.pivot_table(index=keys, columns="TRT01P", values="N",
                      aggfunc="sum", fill_value=0, observed=False)
    # ★ 确保所有治疗组列都存在（即使某组为 0）
    for t in TRT_LEVELS:
        if t not in w.columns:
            w[t] = 0
    return w[TRT_LEVELS].reset_index()


def fmt_cells(w):
    """把计数转成 n (%) 字符串。"""
    out = pd.DataFrame(w[["AEBODSYS", "AEDECOD"]])
    for t in TRT_LEVELS:
        out[t] = [
            "0" if n == 0 else f"{n} ({n / denom[t] * 100:.1f}%)"
            for n in w[t]
        ]
    return out


# ---- Step 4: SOC 层级汇总 ----
soc_w = pivot_counts(subj_soc, ["AEBODSYS"])
soc_pct = fmt_cells(soc_w)
# 按总计数降序（CSR 惯例：按发生频率降序）
soc_pct["_tot"] = soc_w[TRT_LEVELS].sum(axis=1)
soc_pct = soc_pct.sort_values("_tot", ascending=False).drop(columns="_tot")

# ---- Step 5: PT 层级汇总（按 SOC 分组）----
pt_w = pivot_counts(subj_pt, ["AEBODSYS", "AEDECOD"])
pt_pct = pd.DataFrame(pt_w[["AEBODSYS", "AEDECOD"]])
for t in TRT_LEVELS:
    pt_pct[t] = ["0" if n == 0 else f"{n} ({n / denom[t] * 100:.1f}%)" for n in pt_w[t]]

# ---- Step 6: 拼装成 CSR 形态 ----
print("表 2. 治疗中出现的不良事件（按 MedDRA SOC / PT）")
print(f"{'':52s}" + "".join(f"{t[:16]:>20s}" for t in TRT_LEVELS))
print("-" * 120)

# 总行
print(f"{'有任一 TEAE 的受试者':52s}" +
      "".join(f"{int(n_any[t])} ({n_any[t]/denom[t]*100:.1f}%)".rjust(20) for t in TRT_LEVELS))
print()

for _, soc in soc_pct.iterrows():
    name = soc["AEBODSYS"]
    cells = "".join(str(soc[t]).rjust(20) for t in TRT_LEVELS)
    print(f"{name[:50]:52s}{cells}")
    pts = pt_pct[pt_pct["AEBODSYS"] == name]
    for _, pt in pts.iterrows():
        cells = "".join(str(pt[t]).rjust(20) for t in TRT_LEVELS)
        print(f"{'  ' + str(pt['AEDECOD'])[:48]:52s}{cells}")

def export_table(df, path, sheet="Table", col_widths=None):
    """导出到 Excel，做基本的排版美化。"""
    from openpyxl.styles import Alignment, Font, PatternFill
    from openpyxl.utils import get_column_letter

    with pd.ExcelWriter(path, engine="openpyxl") as writer:
        df.to_excel(writer, sheet_name=sheet, index=False)
        ws = writer.sheets[sheet]

        # 表头加粗 + 灰底 + 自动换行
        for c in range(1, len(df.columns) + 1):
            cell = ws.cell(row=1, column=c)
            cell.font = Font(bold=True)
            cell.fill = PatternFill("solid", fgColor="D9D9D9")
            cell.alignment = Alignment(horizontal="center", vertical="center",
                                      wrap_text=True)
        # 冻结表头 + 列宽
        ws.freeze_panes = "A2"
        for i, col in enumerate(df.columns, start=1):
            width = col_widths.get(col, 18) if col_widths else max(
                12, min(45, int(df[col].astype(str).str.len().max() or 12) + 2))
            ws.column_dimensions[get_column_letter(i)].width = width

def to_html_table(df, path, title="", footnote=""):
    """输出带样式的 HTML 表格，可直接贴进报告或邮件。"""
    html = f"""<!DOCTYPE html>
<html><head><meta charset="utf-8"><title>{title}</title>
<style>
  body {{ font-family: "Segoe UI", "Microsoft YaHei", sans-serif; margin: 24px; }}
  table {{ border-collapse: collapse; font-size: 13px; }}
  th, td {{ border: 1px solid #bbb; padding: 4px 10px; }}
  th {{ background: #f0f2f5; text-align: center; }}
  td.num {{ text-align: right; }}
  caption {{ text-align: left; font-weight: bold; margin-bottom: 8px; }}
  tfoot td {{ font-size: 12px; color: #666; border: none; padding-top: 10px; }}
</style></head><body>
<table>
<caption>{title}</caption>
{df.to_html(index=False, escape=False, border=0)}
<tfoot><tr><td colspan="{len(df.columns)}">{footnote}</td></tr></tfoot>
</table></body></html>"""
    Path(path).write_text(html, encoding="utf-8")

# pip install great-tables
from great_tables import GT

gt_tbl = (GT(df)
          .tab_header(title="表 1. 人口学与基线特征",
                      subtitle="安全性人群")
          .cols_align(align="center", columns=[c for c in df.columns if c not in ("变量", "类别")])
          .tab_source_note("数据来源：CDISC 试点项目 CDISCPILOT01（公开数据）")
          .tab_options(table_font_size="12px"))
gt_tbl.save("outputs/table01.html")

def make_listing(adae, out_csv, vars_=None, sort_by=("USUBJID", "ASTDT")):
    """生成受试者级的不良事件列表（对应 CSR 的 Listing）。"""
    vars_ = vars_ or ["USUBJID", "TRTA", "AEDECOD", "AESEV", "AESER",
                      "ASTDT", "AENDT", "AEREL", "AEOUT"]
    lis = adae[vars_].sort_values(list(sort_by))
    # 日期格式化（对应 SAS 的 yymmdd10.）
    for c in [c for c in vars_ if c.endswith("DT")]:
        lis[c] = pd.to_datetime(lis[c], errors="coerce").dt.strftime("%Y-%m-%d")
    lis.to_csv(out_csv, index=False, encoding="utf-8-sig")
    return lis

"""一个典型的 TFL 生产脚本骨架（对应公司里的 'TLF 驱动程序'）"""

def main():
    out_dir = BASE / "outputs"
    out_dir.mkdir(exist_ok=True)

    # 1) 载入
    adsl = pd.read_csv(SAMPLES / "adsl.csv")
    adae = pd.read_csv(SAMPLES / "adae.csv")

    # 2) 数据准备（人群、治疗组顺序、派生变量）
    adsl, ae_saf = prepare(adsl, adae)

    # 3) 生产表
    t1 = make_table01(adsl)                  # 人口学
    t2 = make_table02(ae_saf, adsl)          # TEAE SOC/PT
    l1 = make_listing(adae)                  # AE 列表

    # 4) 输出
    t1.to_csv(out_dir / "t01_demographics.csv", index=False, encoding="utf-8-sig")
    t2.to_csv(out_dir / "t02_teae.csv", index=False, encoding="utf-8-sig")
    l1.to_csv(out_dir / "l01_ae_listing.csv", index=False, encoding="utf-8-sig")

    # 5) 日志（对应 SAS 日志，必须留痕）
    logging.info("TFL 生产完成，共 3 个输出")


if __name__ == "__main__":
    main()