alterlab-scikit-survival

scikit-survival：Python での生存時間分析

概要

scikit-survival は、scikit-learn を基盤とする生存時間分析用の Python ライブラリです。イベント発生までの時間を分析するための専用ツールを提供し、一部の観測値が部分的にしか分からない打ち切りデータ特有の課題に対応します。

生存時間分析は、共変量とイベント発生までの時間の関連を明らかにすることを目的とし、打ち切り記録（特に、観察期間中に参加者がイベントを経験しなかった研究から得られる右打ち切りデータ）を考慮します。

このスキルを使用する場面

以下の場合にこのスキルを使用します：

• 生存時間分析またはイベント発生までの時間のモデリングを行う場合
• 右打ち切りデータを扱う場合（scikit-survival の対象範囲です。左打ち切りデータや区間打ち切りデータには別のツールが必要です）
• Cox 比例ハザードモデル（標準またはペナルティ付き）を当てはめる場合
• アンサンブル生存モデル（ランダム生存フォレスト、勾配ブースティング）を構築する場合
• 生存サポートベクターマシンを学習させる場合
• 生存モデルの性能を評価する場合（一致指数、ブライアスコア、時間依存 AUC）
• Kaplan-Meier 曲線または Nelson-Aalen 曲線を推定する場合
• 競合リスクを分析する場合
• 生存データを前処理する場合、または生存データセットの欠損値を扱う場合
• scikit-survival ライブラリを使用する分析を行う場合

該当しない場合

シナリオ	代わりに使用
打ち切りのない、固定された追跡時点での二値アウトカム（ロジスティック回帰、オッズ比）	alterlab-statsmodels
イベント発生までの時間をアウトカムとしない標準的な回帰または分類	alterlab-scikit-learn
生存曲線とハザード比を含む、臨床意思決定支援レポートまたはバイオマーカーで層別化したコホート文書の作成	alterlab-clinical-decision

scikit-survival は予測を重視しています。CoxPHSurvivalAnalysis は係数を返しますが、標準誤差、信頼区間、p 値は返しません。ハザード比の推論に基づく報告が必要な場合は、lifelines の CoxPHFitter または statsmodels の PHReg（statsmodels.duration.hazard_regression）を使用してください。

主な機能

1. モデルの種類と選択

scikit-survival は複数のモデルファミリーを提供しており、それぞれ異なるシナリオに適しています。

Cox 比例ハザードモデル

用途：係数の解釈が可能な標準的な生存時間分析

• CoxPHSurvivalAnalysis：基本的な Cox モデル
• CoxnetSurvivalAnalysis：高次元データ向けの Elastic Net を用いたペナルティ付き Cox モデル
• IPCRidge：加速故障時間モデル用のリッジ回帰

参照：Cox モデル、正則化、解釈に関する詳しいガイダンスは references/cox-models.md を参照してください。

アンサンブル手法

用途：複雑な非線形関係に対する高い予測性能

• RandomSurvivalForest：堅牢なノンパラメトリック・アンサンブル手法
• GradientBoostingSurvivalAnalysis：最高の性能を目指す、木ベースのブースティング
• ComponentwiseGradientBoostingSurvivalAnalysis：特徴量選択を伴う線形ブースティング
• ExtraSurvivalTrees：追加の正則化を行う、極端にランダム化された決定木

参照：アンサンブル手法、ハイパーパラメーター調整、各モデルの使い分けに関する包括的なガイダンスは references/ensemble-models.md を参照してください。

生存サポートベクターマシン

用途：マージンベースの学習を用いる中規模データセット

• FastSurvivalSVM：速度を重視して最適化された線形 SVM
• FastKernelSurvivalSVM：非線形関係に対応するカーネル SVM
• HingeLossSurvivalSVM：ヒンジ損失を用いる SVM
• ClinicalKernelTransform：臨床データと分子データを組み合わせた専用カーネル

参照：SVM、カーネルの選択、ハイパーパラメーター調整に関する詳しいガイダンスは references/svm-models.md を参照してください。

モデル選択の決定木

Start
├─ High-dimensional data (p > n)?
│  ├─ Yes → CoxnetSurvivalAnalysis (elastic net)
│  └─ No → Continue
│
├─ Need interpretable coefficients?
│  ├─ Yes → CoxPHSurvivalAnalysis or ComponentwiseGradientBoostingSurvivalAnalysis
│  └─ No → Continue
│
├─ Complex non-linear relationships expected?
│  ├─ Yes
│  │  ├─ Large dataset (n > 1000) → GradientBoostingSurvivalAnalysis
│  │  ├─ Medium dataset → RandomSurvivalForest or FastKernelSurvivalSVM
│  │  └─ Small dataset → RandomSurvivalForest
│  └─ No → CoxPHSurvivalAnalysis or FastSurvivalSVM
│
└─ For maximum performance → Try multiple models and compare

2. データの準備と前処理

モデリングの前に、生存データを適切に準備してください。

生存アウトカムの作成

from sksurv.util import Surv

# From separate arrays
y = Surv.from_arrays(event=event_array, time=time_array)

# From DataFrame
y = Surv.from_dataframe('event', 'time', df)

組み込みデータセットでは、それぞれ独自のフィールド名とカテゴリ特徴量が使われています。 load_breast_cancer() は ('e.tdm', 't.tdm')、load_gbsg2() は ('cens', 'time')、load_whas500() は ('fstat', 'lenfol')、load_veterans_lung_cancer() は ('Status', 'Survival_in_days') のフィールドを返します。また、これらの X には、スケーラーやほとんどの推定器で扱えない pandas のカテゴリカル列が含まれています。y['event'] / y['time'] をインデックス参照するコードを再利用する前に、両方のフィールド名を標準化してください。

from sksurv.preprocessing import encode_categorical

X = encode_categorical(X)                          # one-hot encode categorical columns
event_field, time_field = y.dtype.names
y = Surv.from_arrays(event=y[event_field], time=y[time_field])   # fields 'event', 'time'

重要な前処理手順

1. 欠損値の処理：特徴量に対する補完方法
2. カテゴリ変数のエンコード：ワンホットエンコーディングまたはラベルエンコーディング
3. 特徴量の標準化：SVM と正則化 Cox モデルでは特に重要です
4. データ品質の検証：負の時間や、特徴量あたりのイベント数が十分かを確認します
5. 訓練データとテストデータの分割：各分割で打ち切り率が同程度になるようにします

参照：references/data-handling.md では、前処理のワークフロー、データ検証、ベストプラクティスについて詳しく説明しています

3. モデル評価

生存モデルでは適切な評価が重要です。打ち切りを考慮した適切な指標を使用してください。

一致指数（C-index）

順位付け／識別の主要な指標です。

• Harrell の C-index：打ち切りが少ない場合（40% 未満）に使用します
• Uno の C-index：打ち切りが中程度から多い場合（40% 超）に使用します。より頑健です

from sksurv.metrics import concordance_index_censored, concordance_index_ipcw

# Harrell's C-index
c_harrell = concordance_index_censored(y_test['event'], y_test['time'], risk_scores)[0]

# Uno's C-index (recommended)
c_uno = concordance_index_ipcw(y_train, y_test, risk_scores)[0]

時間依存 AUC

特定の時点における識別性能を評価します。

from sksurv.metrics import cumulative_dynamic_auc

times = [365, 730, 1095]  # 1, 2, 3 years
auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk_scores, times)

ブライアスコア

識別性能とキャリブレーションの両方を評価します。

from sksurv.metrics import integrated_brier_score

ibs = integrated_brier_score(y_train, y_test, survival_functions, times)

参照：references/evaluation-metrics.md では、包括的な評価方法、指標の選び方、交差検証でのスコアラーの使用方法について説明しています

4. 競合リスク分析

複数の相互排他的なイベントタイプが存在する状況を扱います。

from sksurv.nonparametric import cumulative_incidence_competing_risks

# Pass SEPARATE arrays: integer-coded event status (0=censored, 1, 2, ...)
# and the observed time. Do NOT collapse the status to a boolean.
times, cif = cumulative_incidence_competing_risks(event_status, time)
# cif[0] = total risk (any event); cif[1:] = CIF for each event type k
cif_event1, cif_event2 = cif[1], cif[2]

競合リスクを使用するケース：

• 複数の相互排他的なイベントタイプが存在する場合（例：死因が複数ある場合）
• あるイベントが発生すると、ほかのイベントが発生しなくなる場合
• 特定のイベントタイプの確率推定が必要な場合

参照：references/competing-risks.md では、競合リスクの手法、原因別ハザードモデル、解釈について詳しく説明しています

5. ノンパラメトリック推定

パラメトリックな仮定を置かずに生存関数を推定します。

Kaplan-Meier 推定量

from sksurv.nonparametric import kaplan_meier_estimator

time, survival_prob = kaplan_meier_estimator(y['event'], y['time'])

Nelson-Aalen 推定量

from sksurv.nonparametric import nelson_aalen_estimator

time, cumulative_hazard = nelson_aalen_estimator(y['event'], y['time'])

一般的なワークフロー

ワークフロー 1：標準的な生存分析

from sksurv.datasets import load_breast_cancer
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.metrics import concordance_index_ipcw
from sksurv.preprocessing import encode_categorical
from sksurv.util import Surv
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# 1. Load and prepare data (encode categoricals; rename outcome fields to event/time)
X, y = load_breast_cancer()
X = encode_categorical(X)
event_field, time_field = y.dtype.names          # ('e.tdm', 't.tdm') for this dataset
y = Surv.from_arrays(event=y[event_field], time=y[time_field])
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y['event'])

# 2. Preprocess
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 3. Fit model (alpha > 0 adds a ridge penalty; unpenalized Cox fails when
#    features are many relative to events, as in this 198-patient, 80-feature dataset)
estimator = CoxPHSurvivalAnalysis(alpha=0.1)
estimator.fit(X_train_scaled, y_train)

# 4. Predict
risk_scores = estimator.predict(X_test_scaled)

# 5. Evaluate
c_index = concordance_index_ipcw(y_train, y_test, risk_scores)[0]
print(f"C-index: {c_index:.3f}")

ワークフロー 2：特徴選択を伴う高次元データ

IPCW スコアラーラッパー（as_concordance_index_ipcw_scorer、as_integrated_brier_score_scorer、as_cumulative_dynamic_auc_scorer）は推定器をラップし、その .score() メソッドをオーバーライドします。これらを scoring= に渡すことはありません。ラップしたオブジェクトを GridSearchCV の推定器として渡し、チューニングするパラメーターには estimator__ を接頭辞として付けてください。有効な scoring='concordance_index_ipcw' という文字列はありません。

CoxnetSurvivalAnalysis は正則化パス全体に適合します。coef_ の形状は (n_features, n_alphas) で、alphas_ にあるペナルティごとに列が 1 つずつあります。そのパスから単一のペナルティをチューニングし、選択した列の非ゼロ係数を確認してください。

import numpy as np
from sksurv.linear_model import CoxnetSurvivalAnalysis
from sklearn.model_selection import GridSearchCV
from sksurv.metrics import as_concordance_index_ipcw_scorer

# 1. Fit once to get the penalty path (l1_ratio near 1 = lasso-like sparsity)
path = CoxnetSurvivalAnalysis(l1_ratio=0.9, alpha_min_ratio=0.01).fit(X_train_scaled, y_train)

# 2. Wrap the estimator so .score() uses Uno's C-index, then tune one alpha per fit.
#    tau caps the evaluation horizon to avoid unstable IPCW weights in the tail.
wrapped = as_concordance_index_ipcw_scorer(
    CoxnetSurvivalAnalysis(l1_ratio=0.9, fit_baseline_model=True), tau=y_train['time'].max())
param_grid = {'estimator__alphas': [[a] for a in path.alphas_]}
# With p >> n the smallest penalties can fail to converge; GridSearchCV reports them as
# FitFailedWarning / NaN scores and still selects among the penalties that fit.
cv = GridSearchCV(wrapped, param_grid, cv=5)
cv.fit(X_train_scaled, y_train)

# 3. Identify selected features (unwrap to reach the Coxnet estimator; coef_ is (n_features, 1))
best_model = cv.best_estimator_.estimator_
selected_features = X.columns[np.flatnonzero(best_model.coef_[:, 0])]

ワークフロー 3：最高の性能を目指すアンサンブル手法

from sksurv.ensemble import GradientBoostingSurvivalAnalysis
from sklearn.model_selection import GridSearchCV
from sksurv.metrics import as_concordance_index_ipcw_scorer, concordance_index_ipcw

# 1. Define parameter grid
param_grid = {
    'learning_rate': [0.01, 0.05, 0.1],
    'n_estimators': [100, 200, 300],
    'max_depth': [3, 5, 7]
}

# 2. Grid search (wrap estimator so .score() is Uno's C-index; prefix params)
gbs = GradientBoostingSurvivalAnalysis(random_state=42)
wrapped = as_concordance_index_ipcw_scorer(gbs, tau=y_train['time'].max())
param_grid = {f'estimator__{k}': v for k, v in param_grid.items()}
cv = GridSearchCV(wrapped, param_grid, cv=5, n_jobs=-1)
cv.fit(X_train, y_train)

# 3. Evaluate best model on held-out test set
best_model = cv.best_estimator_.estimator_
risk_scores = best_model.predict(X_test)
c_index = concordance_index_ipcw(y_train, y_test, risk_scores)[0]

ワークフロー 4：包括的なモデル比較

from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.ensemble import RandomSurvivalForest, GradientBoostingSurvivalAnalysis
from sksurv.svm import FastSurvivalSVM
from sksurv.metrics import concordance_index_ipcw, integrated_brier_score

# Define models
models = {
    'Cox': CoxPHSurvivalAnalysis(alpha=0.1),
    'RSF': RandomSurvivalForest(n_estimators=100, random_state=42),
    'GBS': GradientBoostingSurvivalAnalysis(random_state=42),
    'SVM': FastSurvivalSVM(random_state=42)
}

# Evaluate each model
results = {}
for name, model in models.items():
    model.fit(X_train_scaled, y_train)
    risk_scores = model.predict(X_test_scaled)
    c_index = concordance_index_ipcw(y_train, y_test, risk_scores)[0]
    results[name] = c_index
    print(f"{name}: C-index = {c_index:.3f}")

# Select best model
best_model_name = max(results, key=results.get)
print(f"\nBest model: {best_model_name}")

scikit-learn との統合

scikit-survival は scikit-learn のエコシステムに完全に統合されています。

from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.metrics import as_concordance_index_ipcw_scorer

# Build the pipeline, then wrap it so .score() is Uno's C-index.
# Wrapping the whole pipeline keeps scaling inside each CV fold (no leakage).
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('model', CoxPHSurvivalAnalysis())
])
wrapped = as_concordance_index_ipcw_scorer(pipeline, tau=y['time'].max())

# Cross-validation uses the wrapped estimator's .score(); leave scoring=None
scores = cross_val_score(wrapped, X, y, cv=5)

# Grid search: params live under estimator__ (wrapper) then the pipeline step
param_grid = {'estimator__model__alpha': [0.1, 1.0, 10.0]}
cv = GridSearchCV(wrapped, param_grid, cv=5)
cv.fit(X, y)

ベストプラクティス

1. SVM と正則化 Cox モデルでは必ず特徴量を標準化します
2. 打ち切り率が 40％を超える場合は、Harrell の C-index ではなく Uno の C-index を使用します
3. 複数の評価指標（C-index、統合 Brier スコア、時間依存 AUC）を報告します
4. Cox モデルでは比例ハザードの仮定を確認します
5. 適切なスコアラーを使用し、ハイパーパラメーターの調整には交差検証を使用します
6. モデリングの前にデータ品質を検証します（負の時間がないか、特徴量あたりのイベント数が十分かを確認します）
7. 最適な性能を見つけるために複数のモデルタイプを比較します
8. Random Survival Forests では組み込みの重要度ではなく、置換重要度を使用します
9. 複数のイベントタイプがある場合は競合リスクを考慮します
10. 分析では打ち切りのメカニズムと割合を記録します

よくある落とし穴

1. 打ち切り率が高い場合に Harrell の C-index を使用する → Uno の C-index を使用します
2. SVM の特徴量を標準化しない → 必ず標準化します
3. concordance_index_ipcw に y_train を渡し忘れる → IPCW の計算に必要です
4. 競合イベントを打ち切りとして扱う → 競合リスクの手法を使用します
5. 特徴量あたりのイベント数が十分か確認しない → 目安は特徴量あたり 10 件以上です
6. RSF で組み込みの特徴量重要度を使用する → 置換重要度を使用します
7. 比例ハザードの仮定を無視する → 検証するか、別のモデルを使用します
8. as_concordance_index_ipcw_scorer() を scoring= に渡す → これは推定器をラップし（.score() をオーバーライドします）、ラップ済みオブジェクトを推定器として渡し、パラメーターには estimator__ を付けます

バージョンに関する注意事項

• scikit-survival 0.28（2026 年 7 月）では scikit-learn ≥ 1.9 が必要です。すべての推定器で polars DataFrame が使えるようになり、GradientBoostingSurvivalAnalysis から criterion パラメーターが削除されました。0.27 では pandas 3 に対応しました。
• CoxnetSurvivalAnalysis では 0 < l1_ratio ≤ 1 が必要です。純粋なリッジペナルティを使用するには、CoxPHSurvivalAnalysis(alpha=...) を使用します。

参考ファイル

このスキルには、特定のトピックに関する詳細な参考ファイルが含まれています。

• references/cox-models.md：Cox 比例ハザードモデル、ペナルティ付き Cox（CoxNet）、IPCRidge、正則化戦略、解釈に関する完全なガイド
• references/ensemble-models.md：Random Survival Forests、Gradient Boosting、ハイパーパラメーター調整、特徴量重要度、モデル選択
• references/evaluation-metrics.md：一致指数（Harrell と Uno の比較）、時間依存 AUC、Brier スコア、包括的な評価パイプライン
• references/data-handling.md：データの読み込み、前処理のワークフロー、欠損データの処理、特徴量のエンコード、検証チェック
• references/svm-models.md：生存時間サポートベクターマシン、カーネルの選択、臨床カーネル変換、ハイパーパラメーター調整
• references/competing-risks.md：競合リスク分析、累積発生関数、原因別ハザードモデル

特定のタスクで詳細な情報が必要な場合は、これらの参考ファイルを読み込んでください。

その他のリソース

• 公式ドキュメント：https://scikit-survival.readthedocs.io/
• GitHub リポジトリ：https://github.com/sebp/scikit-survival
• 組み込みデータセット：練習用データセット（GBSG2、WHAS500、veterans lung cancer など）には sksurv.datasets を使用します
• API リファレンス：クラスと関数の完全な一覧は https://scikit-survival.readthedocs.io/en/stable/api/index.html にあります

クイックリファレンス：主なインポート

# Models
from sksurv.linear_model import CoxPHSurvivalAnalysis, CoxnetSurvivalAnalysis, IPCRidge
from sksurv.ensemble import RandomSurvivalForest, GradientBoostingSurvivalAnalysis
from sksurv.svm import FastSurvivalSVM, FastKernelSurvivalSVM
from sksurv.tree import SurvivalTree

# Evaluation metrics
from sksurv.metrics import (
    concordance_index_censored,
    concordance_index_ipcw,
    cumulative_dynamic_auc,
    brier_score,
    integrated_brier_score,
    as_concordance_index_ipcw_scorer,
    as_integrated_brier_score_scorer
)

# Non-parametric estimation
from sksurv.nonparametric import (
    kaplan_meier_estimator,
    nelson_aalen_estimator,
    cumulative_incidence_competing_risks
)

# Data handling
from sksurv.util import Surv
from sksurv.preprocessing import OneHotEncoder, encode_categorical
from sksurv.datasets import load_gbsg2, load_breast_cancer, load_veterans_lung_cancer

# Kernels
from sksurv.kernels import ClinicalKernelTransform