rwd-cohort-analysis

真实世界数据队列分析流水线

概述

此技能提供一套可复现的流水线，用于分析行政理赔数据或 EHR 数据，以估计治疗比较效果。涵盖从原始理赔表到校正后的生存估计的完整工作流：

1. 使用诊断（ICD-10-CM）、操作（CPT/HCPCS）和药物（NDC）代码列表识别队列
2. 测量用药依从性（PDC 和 MPR）
3. 估计倾向评分并诊断协变量平衡
4. Kaplan-Meier 曲线和 Cox 比例风险回归

代码片段使用 Python（pandas、scikit-learn、lifelines）和 R （MatchIt、cobalt、survival、survminer）。用户负责获取数据、遵守 IRB/隐私要求，以及验证临床代码列表。

不在范围内：从源系统提取数据、对临床记录进行自然语言处理、中断时间序列分析、工具变量分析。

用法

在用户要求以下操作时调用：

• 使用诊断/操作/药物代码从理赔或 EHR 数据构建队列
• 计算用药依从性的 PDC 或 MPR
• 估计倾向评分并检查协变量平衡
• 对基于理赔数据构建的队列运行 Kaplan-Meier 分析或 Cox 回归
• 评估比例风险假设

此技能会输出准确代码。用户需提供数据和代码列表。

响应格式

• 先给出用户需要的命令或代码，再作解释
• 按以下顺序组织：确认输入 → 可运行代码 → 解释关键参数 → 注意事项
• 每项任务提供一个完整可运行的示例；不要展示所有替代方案
• 代码注释应精简且实用（说明做什么，而非为何存在）
• 目标：50–100 行代码，并附简短说明

核心概念

数据模型假设

此流水线预期使用以下表格（列名仅供示例）：

表格	关键列
enrollment	patient_id、enroll_start、enroll_end
diagnoses	patient_id、dx_date、icd10_code、dx_position
procedures	patient_id、proc_date、cpt_code
pharmacy	patient_id、fill_date、ndc_code、days_supply、quantity
demographics	patient_id、birth_date、sex、race

1. 队列识别

定义代码列表

# ICD-10-CM codes for type 2 diabetes (example)
T2D_ICD10 = ["E11.0", "E11.1", "E11.2", "E11.3", "E11.4",
             "E11.5", "E11.6", "E11.8", "E11.9"]

# NDC codes for Drug A (user supplies actual NDCs)
DRUG_A_NDC = ["12345678901", "12345678902"]
DRUG_B_NDC = ["98765432101", "98765432102"]

# CPT codes for outcome procedure (e.g., amputation)
OUTCOME_CPT = ["27590", "27591", "27592"]

筛选有洗脱期的新用药者

import pandas as pd

def identify_new_users(pharmacy, ndc_list, washout_days=180):
    """Identify incident users: first fill after washout with no prior fills."""
    drug_fills = pharmacy[pharmacy["ndc_code"].isin(ndc_list)].copy()
    drug_fills = drug_fills.sort_values(["patient_id", "fill_date"])

    # First-ever fill
    first_fill = drug_fills.groupby("patient_id")["fill_date"].min().reset_index()
    first_fill.columns = ["patient_id", "index_date"]

    # Require continuous enrollment for washout_days before index
    cohort = first_fill.merge(enrollment, on="patient_id")
    cohort["washout_start"] = cohort["index_date"] - pd.Timedelta(days=washout_days)
    cohort = cohort[cohort["enroll_start"] <= cohort["washout_start"]]

    # Exclude patients with any fill of the drug during washout
    cohort = cohort.merge(drug_fills, on="patient_id", suffixes=("", "_rx"))
    prior_fills = cohort[
        (cohort["fill_date"] >= cohort["washout_start"]) &
        (cohort["fill_date"] < cohort["index_date"])
    ]
    prevalent_users = prior_fills["patient_id"].unique()
    cohort = first_fill[~first_fill["patient_id"].isin(prevalent_users)]
    return cohort  # columns: patient_id, index_date

应用纳入/排除标准

def apply_criteria(cohort, diagnoses, enrollment, min_age=18,
                   required_dx_codes=None, excluded_dx_codes=None,
                   pre_index_window=365, min_followup=30):
    """Filter cohort by age, required/excluded diagnoses, enrollment."""
    # Age at index
    cohort = cohort.merge(demographics, on="patient_id")
    cohort["age"] = (cohort["index_date"] - cohort["birth_date"]).dt.days / 365.25
    cohort = cohort[cohort["age"] >= min_age]

    # Required diagnosis in pre-index window
    if required_dx_codes:
        pre_dx = diagnoses[diagnoses["icd10_code"].isin(required_dx_codes)]
        pre_dx = pre_dx.merge(cohort[["patient_id", "index_date"]], on="patient_id")
        pre_dx = pre_dx[
            (pre_dx["dx_date"] >= pre_dx["index_date"] - pd.Timedelta(days=pre_index_window)) &
            (pre_dx["dx_date"] < pre_dx["index_date"])
        ]
        cohort = cohort[cohort["patient_id"].isin(pre_dx["patient_id"].unique())]

    # Exclude patients with certain diagnoses
    if excluded_dx_codes:
        excl = diagnoses[diagnoses["icd10_code"].isin(excluded_dx_codes)]
        excl = excl.merge(cohort[["patient_id", "index_date"]], on="patient_id")
        excl = excl[excl["dx_date"] < excl["index_date"]]
        cohort = cohort[~cohort["patient_id"].isin(excl["patient_id"].unique())]

    # Minimum follow-up enrollment
    cohort = cohort.merge(enrollment, on="patient_id")
    cohort["followup_end"] = cohort[["enroll_end", cohort.columns[-1]]].min(axis=1)
    cohort = cohort[
        (cohort["enroll_end"] - cohort["index_date"]).dt.days >= min_followup
    ]
    return cohort

2. 用药依从性

PDC（覆盖天数比例）

PDC 是首选指标（由 PQA 认可）。它统计测量期内至少有一次处方覆盖的天数，最高按 1.0 计。

def calculate_pdc(pharmacy, cohort, ndc_list, period_days=365):
    """Calculate PDC for each patient over a fixed measurement period."""
    fills = pharmacy[pharmacy["ndc_code"].isin(ndc_list)].copy()
    fills = fills.merge(cohort[["patient_id", "index_date"]], on="patient_id")
    fills["period_end"] = fills["index_date"] + pd.Timedelta(days=period_days)

    # Keep fills within measurement period
    fills = fills[fills["fill_date"] < fills["period_end"]]
    fills["cover_start"] = fills["fill_date"]
    fills["cover_end"] = fills["fill_date"] + pd.to_timedelta(fills["days_supply"], unit="D")
    fills["cover_end"] = fills[["cover_end", "period_end"]].min(axis=1)
    fills["cover_start"] = fills[["cover_start", "index_date"]].max(axis=1)

    # Build day-level coverage (vectorized per patient)
    results = []
    for pid, grp in fills.groupby("patient_id"):
        days_covered = set()
        for _, row in grp.iterrows():
            days_covered.update(pd.date_range(row["cover_start"], row["cover_end"] - pd.Timedelta(days=1)))
        results.append({"patient_id": pid, "pdc": min(len(days_covered) / period_days, 1.0)})
    return pd.DataFrame(results)

MPR（药物持有比例）

def calculate_mpr(pharmacy, cohort, ndc_list, period_days=365):
    """MPR = total days supply dispensed / days in period. Can exceed 1.0."""
    fills = pharmacy[pharmacy["ndc_code"].isin(ndc_list)].copy()
    fills = fills.merge(cohort[["patient_id", "index_date"]], on="patient_id")
    fills["period_end"] = fills["index_date"] + pd.Timedelta(days=period_days)
    fills = fills[(fills["fill_date"] >= fills["index_date"]) &
                  (fills["fill_date"] < fills["period_end"])]
    mpr = fills.groupby("patient_id")["days_supply"].sum().reset_index()
    mpr.columns = ["patient_id", "total_supply"]
    mpr["mpr"] = mpr["total_supply"] / period_days
    return mpr

指标	公式	范围	备注
PDC	覆盖天数 / 时段天数	0–1.0	首选；重叠处方不会重复计数
MPR	总供药天数 / 时段天数	0–∞	提前续药时可能超过 1.0；优先级较低

依从性阈值：PDC ≥ 0.80 是“依从”的标准界值。

3. 倾向评分估计

Python（scikit-learn）

from sklearn.linear_model import LogisticRegression
import numpy as np

def estimate_ps(cohort, covariate_cols, treatment_col="treatment"):
    """Fit logistic regression PS model. Returns cohort with ps column."""
    X = cohort[covariate_cols].fillna(0)
    y = cohort[treatment_col]
    model = LogisticRegression(max_iter=1000, solver="lbfgs", C=1e6)
    model.fit(X, y)
    cohort = cohort.copy()
    cohort["ps"] = model.predict_proba(X)[:, 1]
    return cohort, model

R（MatchIt + cobalt）

library(MatchIt)
library(cobalt)

# Estimate PS and perform 1:1 nearest-neighbor matching
m <- matchit(treatment ~ age + sex + charlson + prior_hosp + prior_rx_count,
             data = cohort, method = "nearest", distance = "glm",
             caliper = 0.2, ratio = 1)
matched_data <- match.data(m)

# Balance diagnostics
bal.tab(m, thresholds = c(m = 0.1))  # SMD threshold
love.plot(m, threshold = 0.1)

4. 平衡性诊断

SMD 计算（Python）

def compute_smd(cohort, covariate_cols, treatment_col="treatment"):
    """Compute standardized mean differences for all covariates."""
    treated = cohort[cohort[treatment_col] == 1]
    comparator = cohort[cohort[treatment_col] == 0]
    results = []
    for col in covariate_cols:
        mean_t = treated[col].mean()
        mean_c = comparator[col].mean()
        var_t = treated[col].var()
        var_c = comparator[col].var()
        pooled_sd = np.sqrt((var_t + var_c) / 2)
        smd = abs(mean_t - mean_c) / pooled_sd if pooled_sd > 0 else 0
        results.append({"covariate": col, "smd": round(smd, 4),
                        "balanced": smd < 0.1})
    return pd.DataFrame(results)

诊断项	目标	未达标时的处理
SMD	所有协变量均 < 0.1	在 PS 模型中加入交互项/非线性项
方差比	0.5–2.0	重新设定 PS 模型或截除极端值
PS 重叠	两组之间有较大重叠	截除不重叠区域；报告截除后的 N
权重分布（IPTW）	均值 ≈ 1.0；最大值 < 10	在第 1/99 百分位处截断

5. IPTW（逆概率治疗加权）

def compute_iptw(cohort, treatment_col="treatment", ps_col="ps",
                 stabilized=True, truncate_pct=(1, 99)):
    """Compute stabilized IPTW weights with truncation."""
    cohort = cohort.copy()
    t = cohort[treatment_col]
    ps = cohort[ps_col]

    if stabilized:
        p_treat = t.mean()
        cohort["iptw"] = np.where(t == 1, p_treat / ps, (1 - p_treat) / (1 - ps))
    else:
        cohort["iptw"] = np.where(t == 1, 1 / ps, 1 / (1 - ps))

    # Truncate extreme weights
    lo, hi = np.percentile(cohort["iptw"], truncate_pct)
    cohort["iptw"] = cohort["iptw"].clip(lo, hi)
    return cohort

6. Kaplan-Meier 生存分析

Python（lifelines）

from lifelines import KaplanMeierFitter
import matplotlib.pyplot as plt

def plot_km(cohort, duration_col="followup_days", event_col="event",
            group_col="treatment", labels=("Drug B", "Drug A")):
    """Plot Kaplan-Meier curves by treatment group."""
    fig, ax = plt.subplots(figsize=(8, 6))
    kmf = KaplanMeierFitter()
    for grp, label in zip([0, 1], labels):
        mask = cohort[group_col] == grp
        kmf.fit(cohort.loc[mask, duration_col], cohort.loc[mask, event_col],
                label=label)
        kmf.plot_survival_function(ax=ax, ci_show=True)
    ax.set_xlabel("Days from index date")
    ax.set_ylabel("Survival probability")
    ax.set_title("Kaplan-Meier Survival Curves")
    plt.tight_layout()
    return fig

R（survival + survminer）

library(survival)
library(survminer)

fit <- survfit(Surv(followup_days, event) ~ treatment, data = cohort)
ggsurvplot(fit, data = cohort, pval = TRUE, risk.table = TRUE,
           xlab = "Days from index date", ylab = "Survival probability",
           legend.labs = c("Drug B", "Drug A"))

7. Cox 比例风险模型

Python（lifelines）

from lifelines import CoxPHFitter

def fit_cox(cohort, duration_col="followup_days", event_col="event",
            covariates=None, weight_col=None):
    """Fit Cox PH model, optionally weighted (IPTW)."""
    cols = [duration_col, event_col] + (covariates or [])
    if weight_col:
        cols.append(weight_col)
    df = cohort[cols].dropna()
    cph = CoxPHFitter()
    cph.fit(df, duration_col=duration_col, event_col=event_col,
            weights_col=weight_col, robust=True if weight_col else False)
    cph.print_summary()
    return cph

R（survival）

# Unweighted (matched cohort)
cox_fit <- coxph(Surv(followup_days, event) ~ treatment + age + sex + charlson,
                 data = matched_data)
summary(cox_fit)

# IPTW-weighted
cox_iptw <- coxph(Surv(followup_days, event) ~ treatment,
                  data = cohort, weights = iptw, robust = TRUE)
summary(cox_iptw)

8. 比例风险假设检验

Schoenfeld 残差（Python）

def check_ph_assumption(cph):
    """Test and plot Schoenfeld residuals for PH assumption."""
    cph.check_assumptions(cohort, p_value_threshold=0.05, show_plots=True)

R

ph_test <- cox.zph(cox_fit)
print(ph_test)       # p < 0.05 → PH assumption violated for that covariate
plot(ph_test)         # visual: residuals should show no trend over time

PH 检验结果	操作
所有 p > 0.05	PH 假设成立；继续
协变量 p < 0.05	对该协变量进行分层，或添加时间交互项
全局 p < 0.05	考虑限制平均生存时间（RMST）或参数 AFT 模型

流程总结

1. Define code lists (ICD-10, NDC, CPT)
2. Identify new users with washout          → cohort table
3. Apply inclusion/exclusion criteria       → filtered cohort
4. Build baseline covariate matrix          → pre-index features
5. Estimate propensity scores               → ps column
6. Match or weight (IPTW)                   → balanced cohort
7. Check balance (SMD < 0.1)               → Love plot / SMD table
8. Calculate adherence (PDC/MPR)            → adherence metrics
9. Define outcome + follow-up               → duration + event columns
10. Kaplan-Meier curves                     → survival plot
11. Cox PH model (adjusted or weighted)     → hazard ratios + CIs
12. Check PH assumption (Schoenfeld)        → assumption diagnostics

常见错误

• 错误： 使用 MPR 而非 PDC 衡量药物依从性 正确： 使用 PDC（覆盖天数比例）作为依从性指标 原因： MPR 会对重叠配药重复计数，数值可能超过 1.0，且未获 PQA 认可

• 错误： 未要求研究期间连续参保 正确： 确认洗脱期和随访期内均连续参保 原因： 参保中断会产生无法观测的时段，导致遗漏事件和配药记录

• 错误： 未在索引日期前设置洗脱期 正确： 要求有一段干净的洗脱期（通常为 180 天），期间未配过研究药物 原因： 没有洗脱期，既往用药者会混入新用药者队列，导致结果偏倚

• 错误： 在倾向评分模型中纳入索引日期后的协变量 正确： PS 模型中仅使用索引日期前的（基线）变量 原因： 索引日期后的变量可能受治疗影响，从而引入碰撞偏倚

• 错误： 匹配或加权后用 p 值评估协变量平衡 正确： 对所有协变量使用标准化均值差（SMD < 0.1） 原因： p 值混淆了平衡程度与样本量，因此不适用于此目的

• 错误： 未检查治疗组之间的倾向评分重叠情况 正确： 绘制两组的 PS 分布图，并截去不重叠区域 原因： 对不重叠区域的估计依赖外推，可靠性较低

• 错误： 未检验比例风险假设就报告 Cox 模型结果 正确： 始终运行 Schoenfeld 残差检验，并检查图中是否存在时间趋势 原因： PH 假设违反会导致风险比产生误导；应考虑 RMST 或分层

• 错误： 在 IPTW 加权分析中使用朴素标准误 正确： 拟合加权模型时使用稳健（三明治）方差估计量 原因： 朴素 SE 会忽略权重估计步骤，从而低估不确定性

• 错误： 过度截尾 IPTW 权重（例如截在第 5/95 百分位） 正确： 谨慎截尾（第 1/99 百分位），并报告不同阈值下的敏感性分析 原因： 过度截尾实际上会取消关键观测的加权，从而重新引入混杂

• 错误： 将所有患者的随访结束时间定义为固定的日历日期 正确： 在以下最早发生的时间结束随访：事件、退保、研究结束或死亡 原因： 不正确的删失会使生存估计产生偏倚，并违反非信息性删失假设

参考文献

• Hernán MA、Robins JM。《因果推断：假如》。Chapman & Hall/CRC，2020，https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/
• Austin PC. 倾向评分方法简介。《Multivariate Behav Res》，2011，https://doi.org/10.1080/00273171.2011.568786
• Pharmacy Quality Alliance。PDC 指标规范，https://www.pqaalliance.org/
• Davidson-Pilon C. lifelines 文档，https://lifelines.readthedocs.io/
• Ho DE 等。MatchIt：用于参数因果推断的非参数预处理。《J Stat Softw》，2011，https://doi.org/10.18637/jss.v042.i08
• Therneau TM、Grambsch PM。《生存数据建模》。Springer，2000