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回帰を実行する
• 比例ハザード仮定を評価する

このスキルは、そのまま実行できるコードを出力します。データとコードリストはユーザーが用意します。

応答形式

• 必要なコマンドまたはコードを最初に提示し、その後に説明する
• 構成は「入力の確認 → 動作するコード → 主要パラメーターの説明 → 注意点」とする
• 各タスクにつき、動作する例を1つ提示する。代替案をすべて示さない
• コードコメントは最小限で実用的にする（存在理由ではなく、何をするかを記述）
• 目安：コード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回の処方で カバーされた日数を数え、上限を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より推奨度は低い

アドヒアランスのしきい値：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	傾向スコアモデルに交互作用項／非線形項を追加する
分散比	0.5～2.0	傾向スコアモデルを再指定するか、極端な値を除外する
傾向スコアの重なり	両群で十分な重なりがある	重なりのない領域を除外し、除外後の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

よくある間違い

• **誤り：**服薬アドヒアランスの指標に PDC ではなく MPR を使用する **正解：**アドヒアランスの指標として 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. 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