Trajectory Analysis Skill (轨迹分析)
Overview
本技能封装了基于组的轨迹模型（GBTM）的完整分析流程，从MIMIC-IV等数据库的纵向数据出发，
经过数据预处理、GBTM轨迹建模、描述性统计、Logistic回归、生存分析，最终产出含三线表和
图表的药学期刊风格Word报告。核心脚本使用Python实现，无需依赖R/lcmm。
Workflow
Step 0: 环境准备
所需Python包（venv内安装，不可全局安装）：
pip install pandas numpy scipy scikit-learn matplotlib python-docx
关键依赖说明：
- scikit-learn.GaussianMixture: 用于GBTM的核心GMM聚类
- scipy.optimize.minimize: 自定义Logistic回归MLE（替代statsmodels）
- scipy.stats: 统计检验（Kruskal-Wallis、卡方、Wald z检验）
- python-docx: 三线表和Word报告生成
- KM曲线和log-rank检验：手动实现，无需lifelines
已知避坑：
- statsmodels 在Python 3.13下安装可能失败（patsy/numpy版本冲突），本技能使用自定义MLE替代
- lifelines 依赖链过长（autograd/formulaic/autograd_gamma），本技能手动实现KM和log-rank
- matplotlib 中文字体需设置 plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
Step 1: 数据加载与预处理
执行 scripts/trajectory_analysis.py，核心逻辑：
- 时间列解析: CSV中的日期格式为
dd/mm/YYYY HH:MM:SS，需用pd.to_datetime(format='%d/%m/%Y %H:%M:%S')解析VCI_starttime；TDM/Scr时间戳用format='mixed', dayfirst=True - 相对天数计算:
(charttime - VCI_starttime).dt.total_seconds() / 86400.0 - 纵向AKI定义: AKI分期是纵向变量（每行代表一个时间点），需取每个患者所有行的
max(aki_stage, aki_stage_creat)，而非首行值。正确AKI率约67.6% - 日级数据构建: 按
np.ceil(day).astype(int).clip(upper=14)分箱到Day 1-14，每日取均值 - 宽格式轨迹矩阵:
pivot_table(index='stay_id', columns='day')，保留 ≥2 个非空时间点的患者 - 线性插补:
df.interpolate(method='linear', axis=1, limit_direction='both')填充Day 1-7缺失值 - 联合队列: TDM和Scr轨迹取交集，确保同一批患者有双变量轨迹
关键参数:
- 分析时间窗: Day 1-7（按研究方案），数据提取可扩展到Day 1-14
- 最小数据点: ≥2（TDM非每日检测，要求≥3会丢失过多患者）
- 离群值过滤: TDM 1.0-100.0 mg/L, Scr 0.1-30.0 mg/dL
Step 2: GBTM轨迹模型构建
使用 scripts/trajectory_analysis.py 中的 gbtm_model() 函数：
- GMM聚类初始化: 使用
sklearn.mixture.GaussianMixture(n_components=n_g, covariance_type='diag', max_iter=500, n_init=10)对Day 1-7的宽格式矩阵直接聚类 - 多项式拟合: 对每个GMM组的成员，使用
np.polyfit(days, weighted_mean, poly_deg=2)拟合二次多项式轨迹曲线 - 模型选择: 依次拟合2-5组模型，计算以下指标：
- BIC:gmm.bic(data)— 越小越好
- AvePP（平均后验概率）: 每组成员的gmm.predict_proba()对角线均值 — 要求 >0.70
- OCC（正确分类优势比）:(AvePP/prop) / ((1-AvePP)/(1-prop))— 要求 >5.0
- 最小组样本量: ≥总样本5%或≥20例 - 最优选择策略: 优先选3组（临床可解释性最优），其次4组；需满足AvePP>0.60且最小组≥15
轨迹命名规则（基于组内均值轨迹的起点值和变化量）:
- TDM: "低浓度-稳定组" / "中浓度-稳定组" / "高浓度-平台组" / "中低浓度-上升组"
- Scr: "肌酐-低位稳定组" / "肌酐-高位稳定组" / "肌酐-持续升高组" / "肌酐-下降组"
Step 3: 描述性统计（表1）
按轨迹组比较基线特征：
- 连续变量: Shapiro-Wilk正态性检验 → 正态用均值±标准差，非正态用中位数(IQR)
- 组间比较: 2组用Mann-Whitney U, ≥3组用Kruskal-Wallis
- 分类变量: 频数(百分比), 卡方检验
- P值格式: P≥0.001保留3位小数, P<0.001显示"<0.001"
Step 4: AKI关联分析（表3 + 表4）
- 表3: 轨迹组 × AKI分期(0/1/2/3)交叉表, 卡方检验比较组间AKI率
- 表4: 三层Logistic回归（自定义MLE实现，见
scripts/logistic_mle.py）: - 模型1: 未校正（仅轨迹组虚拟变量）
- 模型2: 校正年龄、性别
- 模型3: 进一步校正LODS评分、基线eGFR、日剂量(mg/kg)
- 结果: OR(95%CI) + Wald z检验P值
- 参照组: 组1（最低浓度组）
Logistic回归MLE实现要点:
- 目标函数: 负对数似然 -[y*log(p) + (1-y)*log(1-p)]
- 优化器: scipy.optimize.minimize(method='BFGS'), 失败则回退 Nelder-Mead
- 数值稳定性: np.clip(z, -30, 30) 防止exp溢出
- 协方差矩阵: inv(X^T @ W @ X), W = diag(p(1-p))，使用pinv防止奇异
- 标准误: sqrt(diag(cov)), Wald z = beta/se, P = 2(1-Phi(|z|))
Step 5: 生存分析（图7）
手动实现KM曲线和log-rank检验（见 scripts/km_survival.py）:
- KM曲线: 按时间排序，逐步计算 S(t) = prod(1 - d_i/n_i)
- 多组log-rank: 计算每组的观测死亡数O和期望死亡数E，chi2 = sum((O-E)^2/E), df=K-1
- 生存时间: 死亡患者用 (dod - VCI_start)/86400, 截尾患者用28天
Step 6: 图表生成
全部图表使用matplotlib, 300dpi PNG输出:
- 图1: 研究流程图（FancyBboxPatch + 箭头）
- 图2: TDM轨迹图（组均值 + 95%CI填充, 15/20mg/L参考线）
- 图3: Scr轨迹图（组均值 + 95%CI, 1.5x/2x/3x基线参考线）
- 图4: 双变量联合轨迹图（上下双面板, TDM组对应Scr）
- 图5: AKI柱状图（堆叠柱: 0/1/2/3期, 标注AKI率和P值）
- 图6: 森林图（OR对数尺度, 模型3, 红色=显著/蓝色=非显著）
- 图7: KM生存曲线（step函数, log-rank P值）
- 补充图S1: BIC-组数诊断图
- 补充图S2: 后验概率直方图
- 补充图S3: 个体轨迹叠加图
配色: ['#2166AC', '#D6604D', '#4DAF4A', '#FF7F00', '#984EA3']
Step 7: Word报告生成
执行 scripts/generate_report.py, 使用python-docx:
- 三线表: 顶线1.5pt(sz=12), 栏目线0.75pt(sz=6), 底线1.5pt; 无竖线; 表头SimHei加粗
- 图表: 图片居中插入(5.5英寸宽), 图题置于图下方
- 正文字体: SimSun 10.5pt, 标题SimHei加粗
- 结构: 标题页 → 摘要 → 研究背景 → 研究方法 → 研究结果(含全部表图) → 讨论 → 结论 → 图表清单 → 参考文献
Resources
scripts/
trajectory_analysis.py— 主分析脚本: 数据预处理 + GBTM建模 + 描述统计 + AKI关联 + 生存分析 + 全部图表生成logistic_mle.py— 自定义Logistic回归MLE(Wald检验), 替代statsmodelskm_survival.py— 手动实现KM曲线和多组log-rank检验, 替代lifelinesgenerate_report.py— Word报告生成脚本(三线表 + 图表 + 药学期刊风格)
references/
gbtm_methodology.md— GBTM方法论参考: Nagin(2005)模型原理, lcmm包对应关系, Python实现要点kdigo_aki.md— KDIGO 2012 AKI定义和分期标准, 基线肌酐确定优先级strobe_checklist.md— STROBE声明观察性研究报告条目清单
Quick Start
典型使用流程:
# 1. 安装依赖（在venv内）
pip install pandas numpy scipy scikit-learn matplotlib python-docx
# 2. 运行主分析（修改DATA_PATH为实际CSV路径）
python scripts/trajectory_analysis.py
# 3. 生成Word报告
python scripts/generate_report.py
输入要求:
- CSV文件: long-format, 每行一个时间点, 含 stay_id (患者ID), 时间戳列, 实验室值列, AKI分期列, 死亡列
- 研究方案: docx文件, 含纳排标准、变量定义、预期图表清单
输出产物:
- analysis_output/: 中间数据(CSV) + 图表(PNG) + Word报告(.docx)
