第 14 章 · 统计分析与可视化

本章目标：用 Python 做描述统计、常见假设检验、生存分析与临床常用图形。注意：本章不追求"覆盖所有统计方法"，而是建立"知道有什么、知道去哪查"的能力，并强调跨软件结果一致性这一临床特有要求。

14.1 从 PROC 到库：该用什么工具

任务	SAS	Python	备注
描述统计	`PROC MEANS` / `UNIVARIATE`	`pandas` + `scipy.stats`	pandas 覆盖 80%
频数 / 交叉表	`PROC FREQ`（含 Fisher）	`pandas.crosstab` + `scipy.stats.fisher_exact`
t 检验	`PROC TTEST`	`scipy.stats.ttest_ind` / `ttest_rel`
方差分析	`PROC ANOVA` / `GLM`	`statsmodels` / `pingouin`
非参数检验	`PROC NPAR1WAY`	`scipy.stats.mannwhitneyu` 等
线性/逻辑回归	`PROC REG` / `LOGISTIC`	`statsmodels` / `sklearn`	statsmodels 输出更接近 SAS
混合模型 / MMRM	`PROC MIXED`	`statsmodels.MixedLM` / `pymer4`	⚠️ 需仔细核对
生存分析	`PROC LIFETEST` / `PHREG`	`lifelines`	见 14.4
重复测量	`PROC MIXED`	`statsmodels`
样本量/检验效能	`PROC POWER`	`statsmodels.stats.power`

Shell
pip install scipy statsmodels lifelines pingouin

📌 你的已有资产：你在 jinbeiwang/mmrm-guide 和 lmm-notes 里已经研究过 MMRM / LMM。Python 侧的对应实现是 statsmodels.MixedLM（通用 LMM）与 pymer4（更接近 SAS 的 PROC MIXED 语法）。 跨语言实现一致性可参考 PHUSE 的 CAMIS 项目（https://psiaims.github.io/CAMIS/）——它专门比对 SAS/R/Python 的同一方法实现。

14.2 描述统计与假设检验

Python
import numpy as np
import pandas as pd
from scipy import stats

adsl = pd.read_csv("data/samples/adsl.csv")

# ---- 描述统计（对应 PROC MEANS）----
age = adsl["AGE"].dropna()
print(f"N={len(age)}  Mean={age.mean():.2f}  SD={age.std(ddof=1):.2f}  "
      f"Median={age.median():.1f}  Min={age.min():.0f}  Max={age.max():.0f}")

# ---- 分组描述统计（对应 PROC MEANS + CLASS）----
summary = adsl.groupby("TRT01P")["AGE"].agg(
    N="count", Mean="mean", SD=lambda s: s.std(ddof=1),
    Median="median", Min="min", Max="max",
).round(2)

# ---- 正态性检验（对应 PROC UNIVARIATE 的 NORMAL 选项）----
# Shapiro-Wilk（SAS 里是 W 检验）
stat, p = stats.shapiro(adsl["AGE"].dropna())
print(f"Shapiro-Wilk W={stat:.4f}, p={p:.4f}  "
      f"-> {'不符合正态分布' if p < 0.05 else '未发现偏离正态'}")

# ---- 两独立样本 t 检验（对应 PROC TTEST）----
g1 = adsl.loc[adsl["TRT01P"] == "Placebo", "AGE"].dropna()
g2 = adsl.loc[adsl["TRT01P"] == "Xanomeline High Dose", "AGE"].dropna()

# SAS 默认输出的是"合并方差 t"，R/Python 默认 Welch t —— 必须显式指定！
t_eq, p_eq = stats.ttest_ind(g1, g2, equal_var=True)     # ≈ SAS 默认
t_w, p_w = stats.ttest_ind(g1, g2, equal_var=False)      # Welch
print(f"合并方差 t={t_eq:.4f}, p={p_eq:.4f}")
print(f"Welch    t={t_w:.4f}, p={p_w:.4f}")

# ---- 非参数检验（对应 PROC NPAR1WAY WILCOXON）----
u, p = stats.mannwhitneyu(g1, g2, alternative="two-sided")
print(f"Mann-Whitney U={u:.1f}, p={p:.4f}")

# ---- 卡方 / Fisher（对应 PROC FREQ 的 CHISQ / EXACT）----
ct = pd.crosstab(adsl["TRT01P"], adsl["SEX"])
chi2, p, dof, expected = stats.chi2_contingency(ct)
print(f"Chi-square={chi2:.4f}, df={dof}, p={p:.4f}")

# 2x2 表用 Fisher 精确检验（对应 / exact 选项）
ct2 = pd.crosstab(adsl["SEX"], adsl["AGEGR1"] == ">80")
odds, p_fisher = stats.fisher_exact(ct2)
print(f"Fisher exact p={p_fisher:.4f}, OR={odds:.4f}")

⚠️ 跨软件差异的三个经典陷阱（PHUSE CAMIS 专门在研究这些）： 1. t 检验的方差假设：SAS PROC TTEST 默认合并方差； R 的 t.test 和 scipy 默认 Welch。不指定就会得到不同结果。 2. 列联表的卡方：SAS 默认连续性校正（Yates），scipy 默认不校正（correction=True 才是校正版）。 3. 分位数算法：SAS 与 numpy 的 percentile 在插值方式上有默认差异（numpy 有 method="linear"/"lower"/"higher" 等多种，SAS 用其中一种）。

临床编程的应对：这些差异必须在 ADRG / ADRG-like 文档里 明确定义并说明理由，不能"默认"。

14.3 回归分析：`statsmodels` ≈ SAS 输出

临床场景常用的是"输出一张系数表"，statsmodels 的 output 最接近 SAS。

Python
import statsmodels.api as sm
import statsmodels.formula.api as smf

adsl = pd.read_csv("data/samples/adsl.csv")

# ---- 线性回归（对应 PROC REG / PROC GLM）----
model = smf.ols("AGE ~ C(TRT01P) + C(SEX) + BMIBL", data=adsl).fit()
print(model.summary())

# ---- 逻辑回归（对应 PROC LOGISTIC）----
# 例：预测"是否出现严重不良事件"
adae = pd.read_csv("data/samples/adae.csv")
adsl["HAS_SAE"] = adsl["USUBJID"].isin(
    adae.loc[adae["AESER"] == "Y", "USUBJID"]
).astype(int)

logit = smf.logit("HAS_SAE ~ AGE + C(SEX) + C(TRT01P)", data=adsl).fit(disp=0)

# ---- 输出 OR 与 95% CI（这是 CSR 里真正要的东西）----
import numpy as np
res = pd.DataFrame({
    "coef": logit.params,
    "se": logit.bse,
    "p": logit.pvalues,
})
res["OR"] = np.exp(res["coef"])
res["CI_low"] = np.exp(res["coef"] - 1.96 * res["se"])
res["CI_high"] = np.exp(res["coef"] + 1.96 * res["se"])
print(res.round(4).to_string())

📌 两个必须知道的差异： 1. 置信区间：SAS 的 PROC LOGISTIC 默认用 Wald CI； R 的 glm 默认 profile likelihood CI。 statsmodels 用 Wald，与 SAS 一致。 2. 参考水平：SAS 默认用字母序最大的类别作参考（或 REF= 指定）；pandas 的 C() 默认也是最后一个。 但显式指定永远更安全： C(TRT01P, Treatment(reference='Placebo'))

14.4 生存分析：`lifelines`

ADTTE 数据集就是为这个准备的。

Python
# pip install lifelines
import pandas as pd
from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

adtte = pd.read_csv("data/samples/adtte.csv")
print(adtte["PARAMCD"].value_counts())
print(f"事件数: {(adtte['CNSR'] == 0).sum()}, 删失数: {(adtte['CNSR'] == 1).sum()}")

# 只分析 OS（总生存）
os_df = adtte.query("PARAMCD == 'OS'")

# ---- Kaplan-Meier（对应 PROC LIFETEST）----
kmf = KaplanMeierFitter()
for trt in ["Placebo", "Xanomeline Low Dose", "Xanomeline High Dose"]:
    sub = os_df[os_df["TRT01P"] == trt]
    kmf.fit(sub["AVAL"] / 30.44,          # 天 → 月
            event_observed=1 - sub["CNSR"],  # ★ CNSR=1 是删失，所以事件 = 1-CNSR
            label=trt)
    print(f"{trt:24s} 中位生存 = {kmf.median_survival_time_:.1f} 月")

# ---- log-rank 检验（对应 PROC LIFETEST 的 STRATA）----
a = os_df.query("TRT01P == 'Placebo'")
b = os_df.query("TRT01P == 'Xanomeline High Dose'")
lr = logrank_test(a["AVAL"], b["AVAL"],
                  event_observed_A=1 - a["CNSR"],
                  event_observed_B=1 - b["CNSR"])
print(f"log-rank p = {lr.p_value:.4f}")

# ---- Cox 比例风险（对应 PROC PHREG）----
cox_df = os_df[["AVAL", "CNSR", "AGE", "SEX", "TRT01P"]].copy()
cox_df["EVENT"] = 1 - cox_df["CNSR"]           # ★ 关键：CNSR 语义转换
cox_df = pd.get_dummies(cox_df, columns=["TRT01P"], drop_first=True)
cph = CoxPHFitter()
cph.fit(cox_df[["AVAL", "EVENT", "AGE", "TRT01P_Xanomeline High Dose"]],
        duration_col="AVAL", event_col="EVENT")
cph.print_summary()

🔥 CNSR 的语义转换是生存分析最容易出错的地方： ADaM 里 CNSR=1 表示删失（censored），CNSR=0 表示发生事件。 而 lifelines 的参数是 event_observed（True = 发生事件）。所以必须是 event_observed = 1 - CNSR。

如果忘了取反，你会得到一个"完全颠倒"的结果—— 而且它看起来完全正常（有中位生存时间、有 p 值），不会报错。这类"静默错误"是临床编程中最危险的。

防范方法：拟合后打印事件数，与数据核对： ```python print(f"传给模型的事件数: {cox_df['EVENT'].sum()}，" f"数据中的实际事件数: {(os_df['CNSR'] == 0).sum()}")

两个数字必须相等

```

生存曲线（KM Plot）

Python
import matplotlib.pyplot as plt

# 中文字体（Windows）
plt.rcParams["font.sans-serif"] = ["Microsoft YaHei", "SimHei", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False

fig, ax = plt.subplots(figsize=(8, 5.5))
for trt, color in zip(["Placebo", "Xanomeline Low Dose", "Xanomeline High Dose"],
                      ["#2E5C8A", "#C0392B", "#E67E22"]):
    sub = os_df[os_df["TRT01P"] == trt]
    KaplanMeierFitter().fit(sub["AVAL"] / 30.44,
                            event_observed=1 - sub["CNSR"],
                            label=f"{trt} (n={len(sub)})").plot_survival_function(
        ax=ax, ci_show=False, color=color, linewidth=2)

ax.set_xlabel("时间（月）")
ax.set_ylabel("生存概率")
ax.set_title("图 1. 总生存期 Kaplan-Meier 曲线（安全性人群）")
ax.grid(alpha=0.3, linestyle="--")
ax.legend(loc="upper right", frameon=False)
fig.tight_layout()
fig.savefig("outputs/figure01_km_os.png", dpi=150)

14.5 临床常用图形的 Python 实现

图形对照表

CSR 图形	SAS	Python
KM 生存曲线	`PROC LIFETEST` + `SGPLOT`	`lifelines` + matplotlib
森林图	`SGPLOT` / 自绘	matplotlib `errorbar` / `forestplot`
瀑布图（肿瘤）	`SGPLOT highlow`	matplotlib `barh`
泳道图 / Swimmer plot	`SGPLOT`	matplotlib `barh`
箱线图	`SGPLOT vbox`	`seaborn.boxplot`
均值变化图（Mean±SE）	`SGPLOT series`	matplotlib `errorbar`
移位表热图	`PROC REPORT` + heatmap	`seaborn.heatmap`
AE 频率条形图	`SGPLOT hbar`	matplotlib `barh`

例：均值随访视变化的图（Mean ± SE，MMRM 类分析的标配图）

Python
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

vs = pd.read_csv("data/samples/vs_bp.csv")
sysbp = vs.query("VSTESTCD == 'SYSBP'").copy()
sysbp["VSSTRESN"] = pd.to_numeric(sysbp["VSSTRESN"], errors="coerce")
sysbp["VISITNUM"] = pd.to_numeric(sysbp["VISITNUM"], errors="coerce")

agg = (sysbp.groupby(["TRT_UNKNOWN", "VISITNUM"]) if False else
       sysbp.groupby(["VISITNUM"])["VSSTRESN"]
            .agg(mean="mean", sd=lambda s: s.std(ddof=1), n="count").reset_index())
agg["se"] = agg["sd"] / np.sqrt(agg["n"])

fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(agg["VISITNUM"], agg["mean"], yerr=1.96 * agg["se"],
            marker="o", capsize=4, linewidth=1.8, color="#2E5C8A")
ax.set_xlabel("访视")
ax.set_ylabel("收缩压 (mmHg)，Mean ± 95% CI")
ax.set_title("图 2. 收缩压随访视的变化")
ax.grid(alpha=0.3, linestyle="--")
fig.tight_layout()
fig.savefig("outputs/figure02_sysbp_trend.png", dpi=150)

例：AE 发生频率的水平条形图

Python
top_soc = (adae.drop_duplicates(["USUBJID", "AEBODSYS"])
                .groupby("AEBODSYS")["USUBJID"].nunique()
                .sort_values(ascending=True).tail(10))

fig, ax = plt.subplots(figsize=(9, 5.5))
ax.barh([s.title()[:42] for s in top_soc.index], top_soc.values, color="#C0392B", alpha=0.85)
ax.set_xlabel("受试者数")
ax.set_title("图 3. 治疗中出现的不良事件（Top 10 SOC）")
for i, v in enumerate(top_soc.values):
    ax.text(v + 0.3, i, str(v), va="center", fontsize=9)
ax.grid(axis="x", alpha=0.3, linestyle="--")
fig.tight_layout()
fig.savefig("outputs/figure03_ae_soc.png", dpi=150)

💡 Python 绘图相对 SAS 的真实优势： - 版式控制精细（tight_layout、GridSpec 多面板） - 主题统一容易（定义一次 rcParams） - 交互式输出（plotly 可生成可缩放 HTML，给 DMC 会议用极方便） - 批量出图简单（一个 for 循环出 50 张图，SAS 要写宏）

14.6 移位表（Shift Table）—— 实验室/Vitals 的标准输出

移位表用于展示"基线 → 基线后"的状态迁移（如 Lab 的 LOW/NORMAL/HIGH）。

Python
"""case04_实验室移位表.py —— 生命体征移位表（同一代码可用于 ADLBC）"""
import numpy as np
import pandas as pd

vs = pd.read_csv("data/samples/vs_bp.csv")
vs["VSSTRESN"] = pd.to_numeric(vs["VSSTRESN"], errors="coerce")
vs["VISITNUM"] = pd.to_numeric(vs["VISITNUM"], errors="coerce")

# 只分析收缩压
bp = vs.query("VSTESTCD == 'SYSBP'").copy()

# ---- Step 1: 按临床惯例定义分级（对应 SAS 自定义 format）----
def grade(x):
    """SYSBP 分级：<90 低 / 90-140 正常 / >140 高（示例阈值，实际以 SAP 为准）"""
    if pd.isna(x):
        return np.nan
    if x < 90:
        return "LOW"
    if x <= 140:
        return "NORMAL"
    return "HIGH"

bp["GRADE"] = bp["VSSTRESN"].apply(grade)
GRADE_ORDER = ["LOW", "NORMAL", "HIGH"]

# ---- Step 2: 确定基线与基线后 ----
# 数据集自带 VSBLFL='Y' 标记基线（优先用它）
baseline = bp.loc[bp["VSBLFL"] == "Y", ["USUBJID", "GRADE"]].rename(
    columns={"GRADE": "GRADE_BL"})
post = bp.loc[bp["VISITNUM"] > 1].copy()
post = post.merge(baseline, on="USUBJID", how="inner")

# ---- Step 3: 交叉制表（基线 × 基线后），按受试者去重 ----
# 每位受试者取最后一次基线后访视（惯例之一：END OF TREATMENT）
last_visit = post.groupby("USUBJID")["VISITNUM"].transform("max")
eot = post[post["VISITNUM"] == last_visit].drop_duplicates("USUBJID")

shift = pd.crosstab(
    pd.Categorical(eot["GRADE_BL"], categories=GRADE_ORDER),
    pd.Categorical(eot["GRADE"], categories=GRADE_ORDER),
    dropna=False,
)

# ---- Step 4: 每个基线行内算百分比（行百分比）----
row_tot = shift.sum(axis=1)
shift_pct = (shift.div(row_tot.replace(0, np.nan), axis=0) * 100).round(1)

# ---- Step 5: 组装成 "n (%)" 格式 ----
out = pd.DataFrame(index=shift.index)
for g in GRADE_ORDER:
    if g in shift.columns:
        out[g] = [f"{int(shift.loc[i, g])} ({shift_pct.loc[i, g]:.1f}%)"
                  if row_tot[i] else "-" for i in shift.index]
out.index.name = "基线"
print("表 3. 收缩压分级移位表（基线 vs 治疗结束）")
print(out.to_string())

🧠 移位表最容易搞错的三点： 1. 百分比的分母：是"该基线分类的行合计"还是"该治疗组总人数"？ 两者都有，取决于 SAP，且必须写清楚。 2. 用哪一次"基线后"访视：末次访视？END OF TREATMENT？固定访视？ 3. 同一受试者同一访视有多条记录（重复测量）→ 必须先定规则（取均值？取最后一次？）

14.7 动手练习

描述统计：对 ADSL 的 AGE / BMIBL 分别按治疗组计算 N / Mean / SD / Median / Min / Max，输出成表。
检验差异：用 t 检验比较 Placebo 与 Xanomeline High Dose 的基线 BMI，分别报告合并方差与 Welch 的结果，说明差异。
回归：用逻辑回归分析 AESER='Y'（严重 AE）与 AGE、SEX、治疗组的关系，输出 OR 与 95% CI（注意 SAE 只有 3 例，样本量不足——这本身就是个教学点：何时应该放弃建模）。
生存分析：用 ADTTE 的 PARAMCD='OS' 分别画出三个治疗组的 KM 曲线，做 log-rank 检验。务必验证事件数是否与数据一致。
移位表：把 14.6 节的代码改成"按治疗组分别输出三张移位表"。
（进阶） 把 14.6 节的代码改写成适用于实验室数据的版本（提示：需要先下载 ADLBC，按 PARAMCD 分组，用 ANRLO/ANRHI 参考范围来定义 LOW/NORMAL/HIGH）。

上一章 ← 第 13 章 · ADaM 衍生与 TFL 报表生成 下一章 → 第 15 章 · AI 辅助编程与 SAS→Python 迁移

第 14 章 · 统计分析与可视化

14.1 从 PROC 到库：该用什么工具

14.2 描述统计与假设检验

14.3 回归分析：statsmodels ≈ SAS 输出

14.4 生存分析：lifelines

两个数字必须相等

生存曲线（KM Plot）

14.5 临床常用图形的 Python 实现

图形对照表

例：均值随访视变化的图（Mean ± SE，MMRM 类分析的标配图）

例：AE 发生频率的水平条形图

14.6 移位表（Shift Table）—— 实验室/Vitals 的标准输出

14.7 动手练习

14.3 回归分析：`statsmodels` ≈ SAS 输出

14.4 生存分析：`lifelines`