rwd-cohort-analysis

Compare original and translation side by side

🇺🇸

Original

English

🇨🇳

Translation

Chinese

Real-World Data Cohort Analysis Pipeline

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

Overview

概述

This skill encodes a reproducible pipeline for analyzing administrative claims or EHR data to estimate comparative treatment effects. It covers the full workflow from raw claims tables to adjusted survival estimates:

Cohort identification using diagnosis (ICD-10-CM), procedure (CPT/HCPCS), and drug (NDC) code lists
Medication adherence measurement (PDC and MPR)
Propensity score estimation and balance diagnostics
Kaplan-Meier curves and Cox proportional hazards regression

Code snippets use Python (pandas, scikit-learn, lifelines) and R (MatchIt, cobalt, survival, survminer). The user is responsible for data access, IRB/privacy compliance, and clinical code-list validation.

Out of scope: data extraction from source systems, natural language processing of clinical notes, interrupted time series, instrumental variable analyses.

本 skill 封装了一套可复现的分析流水线，用于分析行政理赔数据或 EHR 数据并估计比较性治疗效果。它涵盖从原始理赔表到经调整生存估计的完整工作流：

使用诊断（ICD-10-CM）、操作（CPT/HCPCS）和药品（NDC）代码列表识别队列
测量用药依从性（PDC 和 MPR）
估计倾向评分并进行平衡诊断
绘制 Kaplan-Meier 曲线并执行 Cox 比例风险回归

代码示例使用 Python（pandas、scikit-learn、lifelines）和 R（MatchIt、cobalt、survival、survminer）。用户负责数据访问、IRB/隐私合规以及临床代码列表验证。

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

Usage

使用方法

Invoke when the user asks to:

Build a cohort from claims or EHR data using diagnosis/procedure/drug codes
Calculate PDC or MPR for medication adherence
Estimate propensity scores and check covariate balance
Run Kaplan-Meier or Cox regression on a claims-derived cohort
Assess the proportional hazards assumption

The skill emits exact code. The user supplies the data and code lists.

在用户提出以下要求时调用：

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

该 skill 输出准确可用的代码。用户需提供数据和代码列表。

Response Format

响应格式

Lead with the command or code the user needs — explain after
Structure as: confirm inputs → working code → key parameters explained → gotchas
One complete working example per task; do not show every alternative
Keep code comments minimal and functional (what, not why-it-exists)
Target: 50-100 lines of code with brief surrounding explanation

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

Core Concepts

核心概念

Data model assumptions

数据模型假设

The pipeline expects these tables (column names are illustrative):

Table	Key columns
`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`

该流水线需要以下表（列名仅为示例）：

表	关键列
`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. Cohort Identification

1. 队列识别

Define code lists

定义代码列表

python

undefined

python

undefined

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"]

undefined

OUTCOME_CPT = ["27590", "27591", "27592"]

undefined

Identify new users with washout

使用洗脱期识别新使用者

python

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

python

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

Apply inclusion/exclusion criteria

应用纳入/排除标准

python

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

python

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. Medication Adherence

2. 用药依从性

PDC (Proportion of Days Covered)

PDC（Proportion of Days Covered）

PDC is the preferred measure (endorsed by PQA). It counts the number of days in the measurement period covered by at least one fill, capped at 1.0.

python

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)

PDC 是首选指标（由 PQA 认可）。它统计测量周期内至少有一次配药覆盖的天数，并将结果上限设为 1.0。

python

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 (Medication Possession Ratio)

MPR（Medication Possession Ratio）

python

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

Metric	Formula	Range	Notes
PDC	Days covered / days in period	0–1.0	Preferred; no double-counting overlapping fills
MPR	Total days supply / days in period	0–∞	Can exceed 1.0 with early refills; less preferred

Adherence threshold: PDC ≥ 0.80 is the standard cutoff for "adherent."

python

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. Propensity Score Estimation

3. 倾向评分估计

Python (scikit-learn)

Python（scikit-learn）

python

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

python

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)

R（MatchIt + cobalt）

library(MatchIt)
library(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)

undefined

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

undefined

4. Balance Diagnostics

4. 平衡诊断

SMD calculation (Python)

SMD 计算（Python）

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)

Diagnostic	Target	Action if failed
SMD	< 0.1 for all covariates	Add interactions/non-linear terms to PS model
Variance ratio	0.5–2.0	Re-specify PS model or trim extremes
PS overlap	Substantial overlap in both groups	Trim non-overlapping regions; report trimmed N
Weight distribution (IPTW)	Mean ≈ 1.0; max < 10	Truncate at 1st/99th percentile

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 (Inverse Probability of Treatment Weighting)

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

python

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

python

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 Survival Analysis

6. Kaplan-Meier 生存分析

Python (lifelines)

Python（lifelines）

python

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

python

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)

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"))

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 Proportional Hazards Model

7. Cox 比例风险模型

Python (lifelines)

Python（lifelines）

python

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

python

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)

R（survival）

undefined

undefined

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)

undefined

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

undefined

8. Proportional Hazards Assumption Check

8. 比例风险假设检查

Schoenfeld residuals (Python)

Schoenfeld residuals（Python）

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)

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 test result	Action
All p > 0.05	PH assumption holds; proceed
Covariate p < 0.05	Stratify on that covariate or add time-interaction term
Global p < 0.05	Consider restricted mean survival time (RMST) or parametric AFT model

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 模型

Pipeline summary

流水线摘要

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

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

Common Mistakes

常见错误

Wrong: Using MPR instead of PDC for medication adherence Right: Use PDC (Proportion of Days Covered) as the adherence metric Why: MPR double-counts overlapping fills, can exceed 1.0, and is not endorsed by PQA
Wrong: Not requiring continuous enrollment during the study period Right: Verify continuous enrollment for both washout and follow-up windows Why: Gaps in enrollment create unobservable periods where events and fills are missed
Wrong: Omitting a washout period before the index date Right: Require a clean washout (typically 180 days) with no prior fills of the study drug Why: Without washout, prevalent users contaminate the new-user cohort and bias results
Wrong: Including post-index covariates in the propensity score model Right: Only use pre-index (baseline) variables in the PS model Why: Post-index variables may be affected by treatment, introducing collider bias
Wrong: Using p-values to assess covariate balance after matching/weighting Right: Use standardized mean differences (SMD < 0.1) for all covariates Why: P-values conflate balance with sample size and are uninformative for this purpose
Wrong: Not checking propensity score overlap between treatment groups Right: Plot PS distributions for both groups and trim non-overlapping regions Why: Estimates in non-overlapping regions are extrapolated and unreliable
Wrong: Reporting Cox model results without testing the proportional hazards assumption Right: Always run Schoenfeld residual tests and inspect plots for time trends Why: Violated PH produces misleading hazard ratios; consider RMST or stratification instead
Wrong: Using naive standard errors with IPTW-weighted analyses Right: Use robust (sandwich) variance estimators when fitting weighted models Why: Naive SEs underestimate uncertainty because they ignore the weight estimation step
Wrong: Truncating IPTW weights too aggressively (e.g., at 5th/95th percentile) Right: Truncate conservatively (1st/99th) and report sensitivity across thresholds Why: Over-truncation reintroduces confounding by effectively unweighting key observations
Wrong: Defining follow-up end as a fixed calendar date for all patients Right: End follow-up at the earliest of: event, disenrollment, end of study, or death Why: Incorrect censoring biases survival estimates and violates non-informative censoring assumptions

错误： 在用药依从性分析中使用 MPR 而不是 PDC 正确： 使用 PDC（Proportion of Days Covered）作为依从性指标 原因： MPR 会重复计算重叠配药，结果可能超过 1.0，且未获 PQA 认可
错误： 未要求研究期间持续参保 正确： 验证洗脱期和随访期内的持续参保情况 原因： 参保中断会造成不可观察时段，导致事件和配药记录遗漏
错误： 省略 index 日期之前的洗脱期 正确： 要求存在干净的洗脱期（通常为 180 天），且研究药物此前无任何配药记录 原因： 没有洗脱期时，存量使用者会污染新使用者队列并使结果产生偏倚
错误： 将 index 后协变量纳入倾向评分模型 正确： PS 模型仅使用 index 前（基线）变量 原因： index 后变量可能受治疗影响，从而引入碰撞因子偏倚
错误： 匹配或加权后使用 p 值评估协变量平衡 正确： 对所有协变量使用标准化均数差（SMD < 0.1） 原因： p 值会将平衡性与样本量混为一谈，不适合用于此目的
错误： 未检查治疗组之间的倾向评分重叠 正确： 绘制两组的 PS 分布，并裁剪不重叠区域 原因： 不重叠区域中的估计属于外推且不可靠
错误： 报告 Cox 模型结果但未检验比例风险假设 正确： 始终运行 Schoenfeld residuals 检验并检查图中的时间趋势 原因： PH 假设不成立时，hazard ratio 可能产生误导；应改用 RMST 或分层分析
错误： 在 IPTW 加权分析中使用朴素标准误 正确： 拟合加权模型时使用 robust（sandwich）方差估计量 原因： 朴素标准误会忽略权重估计步骤，从而低估不确定性
错误： 过度截尾 IPTW 权重（例如截尾至第 5/95 百分位） 正确： 保守截尾至第 1/99 百分位，并报告不同阈值下的敏感性分析 原因： 过度截尾会通过实质性取消关键观测的权重而重新引入混杂
错误： 将所有患者的随访结束定义为同一个固定日历日期 正确： 在事件、退出研究、研究结束或死亡中最早发生时结束随访 原因： 不正确的删失会偏倚生存估计，并违反非信息性删失假设

References

参考资料

Hernán MA, Robins JM. Causal Inference: What If. Chapman & Hall/CRC 2020, https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/
Austin PC. An introduction to propensity score methods. Multivariate Behav Res 2011, https://doi.org/10.1080/00273171.2011.568786
Pharmacy Quality Alliance. PDC measure specifications, https://www.pqaalliance.org/
Davidson-Pilon C. lifelines documentation, https://lifelines.readthedocs.io/
Ho DE et al. MatchIt: nonparametric preprocessing for parametric causal inference. J Stat Softw 2011, https://doi.org/10.18637/jss.v042.i08
Therneau TM, Grambsch PM. Modeling Survival Data. Springer 2000

Hernán MA, Robins JM. Causal Inference: What If. Chapman & Hall/CRC 2020, https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/
Austin PC. An introduction to propensity score methods. Multivariate Behav Res 2011, https://doi.org/10.1080/00273171.2011.568786
Pharmacy Quality Alliance. PDC measure specifications, https://www.pqaalliance.org/
Davidson-Pilon C. lifelines 文档, https://lifelines.readthedocs.io/
Ho DE et al. MatchIt: nonparametric preprocessing for parametric causal inference. J Stat Softw 2011, https://doi.org/10.18637/jss.v042.i08
Therneau TM, Grambsch PM. Modeling Survival Data. Springer 2000