临床 Python 进阶路线图
⌕ /
路线图 › 工具包 clinic/ · 数据与报表

clinic/derive.py

derive.py · 250 行 · 8 KB
clinic.derive —— 派生逻辑(与 SAS 结果对齐)
Python
"""
clinic.derive —— 派生逻辑(与 SAS 结果对齐)

本模块集中实现"跨语言最容易出现差异"的派生逻辑,
所有函数都以**与 SAS 结果一致**为第一优先级。

包含
----
- ``sas_round`` / ``sas_round_series`` : 复刻 SAS ``ROUND`` 的舍入语义
- ``derive_agegr1``      : 年龄分组(含边界归属说明)
- ``parse_partial_date`` : 处理 ISO 8601 部分日期(2014 / 2014-01 / 2014-01-05)
- ``sas_date_to_datetime`` / ``datetime_to_sas_date`` : 与 SAS 日期数字互换
- ``derive_baseline_flag`` : 派生基线记录标记
- ``to_categorical``     : 用 Categorical 固定分组顺序(替代 SAS 的 CLASS 顺序)
"""

from __future__ import annotations

import numpy as np
import pandas as pd

# SAS 日期基准日(SAS 里 '01JAN1960'd = 0)
SAS_EPOCH = pd.Timestamp("1960-01-01")


# --------------------------------------------------------------------------
# 舍入:Python 与 SAS 的头号差异来源
# --------------------------------------------------------------------------
def sas_round(x, unit: float = 1):
    """复刻 SAS 的 ``ROUND(x, unit)``。

    SAS 语义:
    - 第二个参数是**舍入单位**(不是小数位数)
    - 采用**四舍五入**(round half away from zero)

    Python 的 ``round()`` 采用**银行家舍入**(round half to even),
    且第二个参数是小数位数 —— 两者对临床数值都不可直接替换。

    Examples
    --------
    >>> sas_round(2.5, 1)
    3.0
    >>> sas_round(2.675, 0.01)      # Python round(2.675, 2) 会给 2.67
    2.68
    >>> sas_round(1234, 10)
    1230.0
    >>> sas_round(-2.5, 1)
    -3.0
    """
    import math
    from decimal import Decimal, ROUND_HALF_UP

    if x is None:
        return None
    try:
        if math.isnan(float(x)):
            return np.nan
    except (TypeError, ValueError):
        return None

    unit_d = Decimal(str(unit))
    if unit_d == 0:
        return float(x)
    # 用 Decimal 避免二进制浮点误差(2.675 实际存储为 2.67499999...)
    q = (Decimal(str(x)) / unit_d).quantize(Decimal("1"), rounding=ROUND_HALF_UP)
    return float(q * unit_d)


def sas_round_series(s: pd.Series, decimals: int = 0) -> pd.Series:
    """向量化版本的 SAS 风格舍入(四舍五入)。

    比 ``s.apply(sas_round)`` 快得多,适合整列处理。
    ``decimals`` 是**小数位数**(与 SAS 的 unit 语义不同,这里是常用写法)。
    """
    factor = 10 ** decimals
    s = pd.to_numeric(s, errors="coerce")
    return np.sign(s) * np.floor(np.abs(s) * factor + 0.5) / factor


def sas_round_to(s: pd.Series, unit: float) -> pd.Series:
    """按 SAS 的"舍入单位"语义整列舍入(如把 1234 舍入到最近的 10)。"""
    return sas_round_series(s / unit, 0) * unit


# --------------------------------------------------------------------------
# 分组派生:边界必须显式
# --------------------------------------------------------------------------
def derive_agegr1(df: pd.DataFrame, age_var: str = "AGE",
                  out_var: str = "AGEGR1",
                  num_var: str = "AGEGR1N") -> pd.DataFrame:
    """派生年龄分组。

    ★ 边界规则(与 CDISC 试点项目 ADSL 一致,也是最常见的惯例)★

    ==========  =====================
    区间        归属
    ==========  =====================
    AGE < 65    ``'<65'``    (不含 65)
    65 ≤ AGE ≤ 80  ``'65-80'`` (**含 65、含 80**)
    AGE > 80    ``'>80'``    (不含 80)
    ==========  =====================

    .. warning::
       不要用 ``pd.cut(bins=[0, 65, 80, 200], right=False)`` 来实现!
       那样产生的区间是 ``[0,65) / [65,80) / [80,200)``,
       **AGE = 80 会落到 ``>80`` 组**,与临床惯例冲突。
       这是分组派生最隐蔽的错误之一。

    >>> import pandas as pd
    >>> derive_agegr1(pd.DataFrame({"AGE": [64, 65, 80, 81]}))["AGEGR1"].tolist()
    ['<65', '65-80', '65-80', '>80']
    """
    df = df.copy()
    a = pd.to_numeric(df[age_var], errors="coerce")
    df[out_var] = np.select([a < 65, a <= 80], ["<65", "65-80"], default=">80")
    df.loc[a.isna(), out_var] = np.nan
    if num_var:
        df[num_var] = df[out_var].map({"<65": 1, "65-80": 2, ">80": 3})
    return df


def to_categorical(s: pd.Series, levels: list, ordered: bool = True) -> pd.Series:
    """把列转成有序 Categorical —— 用来**固定分组输出顺序**。

    这是复刻 SAS ``CLASS`` / 自定义 format 输出顺序的关键:
    SAS 报表里治疗组按 SAP 指定的顺序出现,而 pandas 默认按字母序。
    不做这一步,报表的行序就会错。

    >>> to_categorical(pd.Series(["B", "A"]), ["A", "B"]).tolist()
    ['B', 'A']           # categories 顺序为 A, B,导致 groupby 按 A, B 输出
    """
    return pd.Categorical(s, categories=levels, ordered=ordered)


# --------------------------------------------------------------------------
# 日期
# --------------------------------------------------------------------------
def sas_date_to_datetime(value):
    """SAS 日期数字 → pandas Timestamp(SAS 以 1960-01-01 为 0)。"""
    return SAS_EPOCH + pd.to_timedelta(value, unit="D")


def datetime_to_sas_date(value):
    """pandas Timestamp / datetime → SAS 日期数字。"""
    return (pd.Timestamp(value) - SAS_EPOCH).days


def parse_partial_date(series: pd.Series) -> pd.DataFrame:
    """解析 ISO 8601 日期并保留"精度"信息。

    临床数据里 ``--DTC`` 常有部分日期(只知道年月)。

    ==============  ===============  ========
    输入            解析结果         精度
    ==============  ===============  ========
    ``2014-01-05``  2014-01-05       DAY
    ``2014-01``     2014-01-01       MONTH
    ``2014``        2014-01-01       YEAR
    ==============  ===============  ========

    调用方应根据精度决定是否参与计算 —— **不要直接用补全后的日期做运算**,
    否则会引入虚假精度。
    """
    s = series.astype("string").str.strip()
    precision = np.select(
        [s.str.len() >= 10, s.str.len() == 7], ["DAY", "MONTH"], default="YEAR"
    )
    filled = s.where(
        s.str.len() >= 10,
        s.str.pad(10, side="right", fillchar="0")
         .str.replace(r"^(\d{4})$", r"\1-01-01", regex=True)
         .str.replace(r"^(\d{4}-\d{2})$", r"\1-01", regex=True),
    )
    return pd.DataFrame(
        {"date": pd.to_datetime(filled, errors="coerce"), "precision": precision}
    )


def add_datetime(df: pd.DataFrame, dtc_var: str,
                 out_var: str | None = None) -> pd.DataFrame:
    """把 ISO 8601 字符日期列转成 datetime,**保留原始列**。

    保留原始值是 CDISC 的要求(原始值必须可追溯),
    同时也是 QC 时对照"原始 vs 派生"的依据。
    """
    out_var = out_var or f"{dtc_var}_DT"
    df = df.copy()
    df[out_var] = pd.to_datetime(df[dtc_var], errors="coerce")
    return df


def derive_baseline_flag(df: pd.DataFrame, date_var: str, ref_var: str,
                         by: list[str], out_var: str = "_BLFL_NEW") -> pd.DataFrame:
    """派生基线标记:**首次给药日当天或之前的最后一次测量**。

    用于数据集没有原生 ``--BLFL`` 的情况。
    ``ref_var`` 通常是 ``TRTSDT``;若为空则视为""受试者尚未给药",
    其所有记录都参与基线候选。
    """
    df = df.copy()
    d = pd.to_datetime(df[date_var], errors="coerce")
    ref = pd.to_datetime(df[ref_var], errors="coerce") if ref_var in df.columns else pd.NaT

    cand = df[(ref.isna()) | (d <= ref)].copy()
    cand["_d"] = d[cand.index]
    idx = cand.dropna(subset=["_d"]).sort_values("_d").groupby(by)["_d"].idxmax()

    df[out_var] = "N"
    df.loc[idx.dropna(), out_var] = "Y"
    return df


# --------------------------------------------------------------------------
# 受试者标识
# --------------------------------------------------------------------------
def split_usubjid(series: pd.Series, names: list[str] | None = None) -> pd.DataFrame:
    """拆解 USUBJID(如 ``01-701-1015`` → study/site/subjid)。

    .. note::
       USUBJID 的拼接规则**因公司/项目而异**(分隔符可能是 ``-`` / ``_``,
       段数也可能是 2~4 段)。生产代码里应先做正则校验,不要盲目 split。
    """
    import re
    names = names or ["STUDY_PART", "SITE_PART", "SUBJ_PART"]
    ext = series.astype(str).str.extract(r"^(?P<a>[^-_]+)[-_](?P<b>[^-_]+)[-_](?P<c>.+)$")
    ext.columns = names
    ok = series.astype(str).str.match(r"^[^-_]+[-_][^-_]+[-_].+$")
    return ext.where(ok, np.nan)


# --------------------------------------------------------------------------
# 人群标记
# --------------------------------------------------------------------------
def derive_saffl(df: pd.DataFrame, trtsdt_var: str = "TRTSDT_DT",
                 arm_var: str = "ARM",
                 excluded_arms: tuple[str, ...] = (
                     "Screen Failure", "Not Assigned", "Not Randomized"),
                 out_var: str = "SAFFL") -> pd.DataFrame:
    """派生安全性人群标记:接受了至少一次研究药物,且不属于排除组。

    这对应 SAS 里最常见的 ``if not missing(trtsdt) and arm not in (...) then SAFFL='Y';``
    """
    df = df.copy()
    has_trt = pd.to_datetime(df[trtsdt_var], errors="coerce").notna() \
        if trtsdt_var in df.columns else pd.Series(False, index=df.index)
    bad_arm = df[arm_var].isin(excluded_arms) if arm_var in df.columns \
        else pd.Series(False, index=df.index)
    df[out_var] = np.where(has_trt & ~bad_arm, "Y", "N")
    return df