scikit-survival-analysis について
scikit-survivalは、右打ち切りデータを扱うscikit-learnベースのPython生存分析ライブラリです。複数のモデル、ノンパラメトリック推定、打ち切りを考慮した評価指標を、scikit-learnのfit/predict APIで利用できます。
主な機能
- 01Cox比例ハザードモデルの適合
- 02生存フォレスト・ブースティング・SVMの学習
- 03Kaplan-Meier生存曲線の推定
- 04競合リスクの分析
- 05打ち切り対応指標による予測評価
活用シーン
- 01臨床試験のイベント時間の分析
- 02信頼性データの生存モデル構築
- 03複数の生存モデルの予測性能比較
- 04高次元生存データでの特徴選択
インストール
コマンドをコピーしてターミナルで実行します。
npx skills add https:/ / github.com/ jaechang-hits/ sciagent-skills --skill scikit-survival-analysisインストール後、Claude Code で/scikit-survival-analysisを入力して読み込まれたか確認します。
同梱ファイルreferences/data_competing_risks.mdreferences/models_evaluation.md
詳細
scikit-survival — 生存時間解析
概要
scikit-survival は、scikit-learn を基盤とするイベント発生までの時間解析用 Python ライブラリです。Cox モデル、アンサンブル手法、生存時間 SVM、ノンパラメトリック推定量を用いて、右打ち切りデータ(イベントがまだ発生していない観測データ)を扱います。すべてのモデルが scikit-learn の fit/predict API に準拠し、Pipeline、交差検証、GridSearchCV と統合できます。
使用する場面
- 右打ち切りデータを用いてイベント発生までの時間をモデル化する場合(臨床試験、信頼性解析)
- Cox 比例ハザードモデル(標準または Elastic Net ペナルティ付き)を当てはめる場合
- アンサンブル生存時間モデル(Random Survival Forest、Gradient Boosting)を構築する場合
- 中規模データセットでマージンベース学習を行う生存時間 SVM を訓練する場合
- 打ち切りを考慮した指標(C-index、Brier スコア、AUC)で生存予測を評価する場合
- ノンパラメトリックな生存曲線(Kaplan-Meier、Nelson-Aalen)を推定する場合
- 累積発生関数を用いて競合リスクを分析する場合
- 特徴量の自動選択(CoxNet L1/L2)を用いて高次元生存データを扱う場合
- より単純なパラメトリックモデル(Weibull、対数正規 AFT)や統計検定(ログランク検定)には、
lifelinesを使用してください - 深層学習による生存時間モデルには、
pycoxまたはtorchlifeを使用してください
前提条件
pip install scikit-survival scikit-learn pandas numpy matplotlib
Python:>= 3.9。依存関係:scikit-learn、numpy、scipy、pandas、joblib、osqp(一部の SVM ソルバーで使用)。
データ形式:生存アウトカムは (event, time) フィールドを持つ NumPy 構造化配列です。イベントはブール値です(True = イベント発生、False = 打ち切り)。時間は正の浮動小数点数です。
事前確認
分析コードを書く前に、ユーザーと以下の点を確認してください。
decisions:
- id: D1
param: timeAndEventColumns
kind: required
source: data
ask: "Which column holds follow-up time, which holds the event indicator, and does the indicator mark the event or censoring?"
default: null
- id: D2
param: eventDefinition
kind: required
source: user
depends_on: [D1]
ask: "Which event is being modelled, and are competing events treated as censoring?"
default: null
- id: D3
param: covariates
kind: required
source: data
ask: "Which variables enter the model, and are any measured after baseline?"
default: null
- id: D4
param: modelFamily
kind: required
source: user
ask: "A proportional-hazards model whose coefficients are interpretable, or an ensemble that predicts better but does not yield hazard ratios?"
default: "Cox proportional hazards"
- id: D5
param: regularization
kind: required
source: user
depends_on: [D3, D4]
ask: "With more covariates than events, should coefficients be penalized - and toward selection or toward shrinkage?"
default: "none"
- id: D6
param: validationScheme
kind: required
source: user
ask: "How should performance be estimated - cross-validation, a held-out split, or an external cohort?"
default: "cross-validation"
- id: D7
param: ensembleHyperparameters
kind: optional_conditional
source: user
depends_on: [D4]
ask: "How many trees, how deep, and at what learning rate?"
default: "100 estimators, depth 3, learning rate 0.1"
skip_if: "a proportional-hazards model was chosen"
D1 と D2 は、何よりも先に確認します。イベント指標を反転させると、イベントが発生しない時間をモデル化することになり、警告なしにすべてのハザード比が逆転するためです。D3 にも別の見落とされやすい誤りがあります。ベースライン測定後の共変量を使うと、アウトカムが予測変数に漏れ込み、見栄えのよい一致指数が得られても再現できなくなります。
クイックスタート
from sksurv.datasets import load_breast_cancer
from sksurv.ensemble import RandomSurvivalForest
from sksurv.metrics import concordance_index_ipcw
from sklearn.model_selection import train_test_split
X, y = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
rsf = RandomSurvivalForest(n_estimators=100, random_state=42)
rsf.fit(X_train, y_train)
risk_scores = rsf.predict(X_test)
c_index = concordance_index_ipcw(y_train, y_test, risk_scores)[0]
print(f"C-index: {c_index:.3f}") # e.g., 0.68
# Individual survival curves
surv_fns = rsf.predict_survival_function(X_test[:2])
for fn in surv_fns:
print(f"5-year survival: {fn(365 * 5):.3f}")
コア API
モジュール 1:データの準備
構造化生存時間配列を作成し、特徴量を前処理します。
import numpy as np
import pandas as pd
from sksurv.util import Surv
from sksurv.preprocessing import OneHotEncoder, encode_categorical
from sksurv.datasets import load_gbsg2, load_breast_cancer
from sklearn.preprocessing import StandardScaler
# Create survival outcome from arrays
event = np.array([True, False, True, True, False])
time = np.array([120.0, 365.0, 200.0, 90.0, 400.0])
y = Surv.from_arrays(event=event, time=time)
print(y.dtype) # [('event', '?'), ('time', '<f8')]
# From DataFrame columns
# y = Surv.from_dataframe("event_col", "time_col", df)
# Load built-in datasets
# Available: load_gbsg2, load_breast_cancer, load_veterans_lung_cancer,
# load_whas500, load_aids, load_flchain
X, y = load_gbsg2()
print(f"Shape: {X.shape}, Events: {y['event'].sum()}, "
f"Censoring rate: {1 - y['event'].mean():.1%}")
# Encode categoricals (survival-aware one-hot)
X_encoded = encode_categorical(X) # auto-detect and encode all categorical cols
# Standardize (critical for Cox and SVM models)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_encoded)
from sksurv.io import loadarff
# Load ARFF format (Weka format)
data = loadarff("survival_data.arff")
X_arff, y_arff = data[0], data[1] # DataFrame, structured array
モジュール 2:Cox 比例ハザード
セミパラメトリックモデル:h(t|x) = h_0(t) * exp(beta^T x)。係数は対数ハザード比として解釈できます。
from sksurv.linear_model import CoxPHSurvivalAnalysis, CoxnetSurvivalAnalysis, IPCRidge
# Standard Cox PH model
cox = CoxPHSurvivalAnalysis(alpha=0.0, ties="breslow")
cox.fit(X_train, y_train)
print(f"Coefficients: {cox.coef_}") # log hazard ratios
# Hazard ratio interpretation: exp(coef) = HR for 1-unit increase
risk_scores = cox.predict(X_test) # Higher = higher risk
# Survival function for individual patients
surv_funcs = cox.predict_survival_function(X_test[:3])
for fn in surv_funcs:
print(f"5-year survival: {fn(365 * 5):.3f}")
# Penalized Cox (elastic net) -- for high-dimensional data (p > n)
coxnet = CoxnetSurvivalAnalysis(
l1_ratio=0.9, # 0=Ridge, 1=Lasso, between=Elastic Net
alpha_min_ratio=0.01, # smallest alpha / largest alpha ratio
n_alphas=100, # steps in regularization path
)
coxnet.fit(X_train, y_train)
# Feature selection: non-zero coefficients
selected = np.where(coxnet.coef_ != 0)[0]
print(f"Selected {len(selected)} / {X_train.shape[1]} features")
# IPCRidge: accelerated failure time model (predicts log survival time)
ipcridge = IPCRidge(alpha=1.0)
ipcridge.fit(X_train, y_train)
log_survival_time = ipcridge.predict(X_test)
モジュール 3:アンサンブル手法
複雑な非線形関係に対応する、ノンパラメトリックな木ベースのモデルです。
from sksurv.ensemble import (
RandomSurvivalForest,
GradientBoostingSurvivalAnalysis,
ComponentwiseGradientBoostingSurvivalAnalysis,
ExtraSurvivalTrees,
)
# Random Survival Forest -- robust, minimal tuning
rsf = RandomSurvivalForest(
n_estimators=200, min_samples_split=10, min_samples_leaf=15,
max_features="sqrt", random_state=42, n_jobs=-1,
)
rsf.fit(X_train, y_train)
risk = rsf.predict(X_test)
# Gradient Boosting -- best performance, needs tuning
gbs = GradientBoostingSurvivalAnalysis(
loss="coxph", # "coxph" or "ipcwls" (AFT)
n_estimators=300, learning_rate=0.05, max_depth=3,
subsample=0.8, dropout_rate=0.1, random_state=42,
)
gbs.fit(X_train, y_train)
# ComponentwiseGB -- linear model with automatic feature selection
cgbs = ComponentwiseGradientBoostingSurvivalAnalysis(
n_estimators=100, learning_rate=0.1,
)
cgbs.fit(X_train, y_train)
print(f"Non-zero coefficients: {np.sum(cgbs.coef_ != 0)}")
# ExtraSurvivalTrees -- more regularized than RSF, faster training
est = ExtraSurvivalTrees(n_estimators=100, random_state=42)
est.fit(X_train, y_train)
# Survival curves from any ensemble model
surv_funcs = rsf.predict_survival_function(X_test[:1])
chf_funcs = rsf.predict_cumulative_hazard_function(X_test[:1])
モジュール 4:生存時間 SVM
生存順位付けのためのマージンベース学習です。特徴量は必ず標準化してください。
from sksurv.svm import FastSurvivalSVM, FastKernelSurvivalSVM, HingeLossSurvivalSVM
# Linear SVM -- fast, for linear relationships
lsvm = FastSurvivalSVM(alpha=1.0, rank_ratio=1.0, max_iter=100, random_state=42)
lsvm.fit(X_train_scaled, y_train)
risk = lsvm.predict(X_test_scaled)
# Kernel SVM -- for non-linear relationships (rbf, poly, sigmoid)
ksvm = FastKernelSurvivalSVM(
alpha=1.0, kernel="rbf", gamma="scale",
max_iter=50, random_state=42,
)
ksvm.fit(X_train_scaled, y_train)
# Hinge loss variant
hsvm = HingeLossSurvivalSVM(alpha=1.0, random_state=42)
hsvm.fit(X_train_scaled, y_train)
# NaiveSurvivalSVM also available but slower (O(n^3))
from sksurv.kernels import ClinicalKernelTransform
# Clinical kernel: combines clinical + molecular features
# Weighs clinical variables separately from high-dimensional molecular data
transform = ClinicalKernelTransform(fit_once=True)
transform.prepare(X_train) # auto-detect clinical features
X_kern = transform.fit_transform(X_train)
モジュール 5:ノンパラメトリック推定
モデルの仮定を置かずに、生存曲線とハザード曲線を推定します。
from sksurv.nonparametric import kaplan_meier_estimator, nelson_aalen_estimator
import matplotlib.pyplot as plt
# Kaplan-Meier survival curve
time_km, surv_prob = kaplan_meier_estimator(y["event"], y["time"])
plt.step(time_km, surv_prob, where="post")
plt.xlabel("Time (days)")
plt.ylabel("Survival probability")
plt.title("Kaplan-Meier Estimate")
plt.savefig("km_curve.png", dpi=150)
# With confidence intervals
time_km, surv_prob, conf_int = kaplan_meier_estimator(
y["event"], y["time"], conf_type="log-log",
)
# Nelson-Aalen cumulative hazard
time_na, cum_hazard = nelson_aalen_estimator(y["event"], y["time"])
# Stratified KM by group
for group_name, mask in [("Treated", treated_mask), ("Control", control_mask)]:
t, s = kaplan_meier_estimator(y[mask]["event"], y[mask]["time"])
plt.step(t, s, where="post", label=group_name)
plt.legend()
モジュール 6:評価指標
識別と較正のための、打ち切りを考慮した指標です。
from sksurv.metrics import (
concordance_index_censored,
concordance_index_ipcw,
cumulative_dynamic_auc,
integrated_brier_score,
brier_score,
as_concordance_index_ipcw_scorer,
as_integrated_brier_score_scorer,
)
import numpy as np
risk_scores = model.predict(X_test)
# Harrell's C-index (simple, biased with high censoring)
c_harrell = concordance_index_censored(
y_test["event"], y_test["time"], risk_scores
)[0]
# Uno's C-index (recommended -- robust to censoring)
c_uno = concordance_index_ipcw(y_train, y_test, risk_scores)[0]
print(f"C-index (Harrell): {c_harrell:.3f}, (Uno): {c_uno:.3f}")
# Time-dependent AUC at clinically relevant timepoints
times = np.array([365, 730, 1095]) # 1, 2, 3 years
auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk_scores, times)
print(f"AUC at 1/2/3yr: {auc}, Mean: {mean_auc:.3f}")
# Integrated Brier Score (measures discrimination + calibration)
surv_funcs = model.predict_survival_function(X_test)
preds = np.row_stack([fn(times) for fn in surv_funcs])
ibs = integrated_brier_score(y_train, y_test, preds, times)
print(f"IBS: {ibs:.3f}") # Lower is better; compare vs KM baseline
モジュール 7:競合リスク
複数の相互に排他的なイベント種別(例:がんによる死亡と心血管疾患による死亡)を扱います。
from sksurv.nonparametric import cumulative_incidence_competing_risks
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.util import Surv
import numpy as np
# Cumulative Incidence Function (CIF) estimation
# y_competing: structured array with integer event codes
# (0=censored, 1=relapse, 2=death_in_remission)
time_pts, cif_relapse, cif_death = cumulative_incidence_competing_risks(y_competing)
import matplotlib.pyplot as plt
plt.step(time_pts, cif_relapse, where="post", label="Relapse")
plt.step(time_pts, cif_death, where="post", label="Death in remission")
plt.xlabel("Time")
plt.ylabel("Cumulative Incidence")
plt.legend()
# Cause-specific hazard models: fit separate Cox per event type
event_types = np.array([0, 1, 2, 1, 0, 2, 1])
times = np.array([10.2, 5.3, 8.1, 3.7, 12.5, 6.8, 4.2])
# Model for relapse (event type 1); treat type 2 as censored
y_relapse = Surv.from_arrays(event=(event_types == 1), time=times)
cox_relapse = CoxPHSurvivalAnalysis()
cox_relapse.fit(X, y_relapse)
# Model for death (event type 2); treat type 1 as censored
y_death = Surv.from_arrays(event=(event_types == 2), time=times)
cox_death = CoxPHSurvivalAnalysis()
cox_death.fit(X, y_death)
print(f"Relapse HR (age): {np.exp(cox_relapse.coef_[0]):.3f}")
print(f"Death HR (age): {np.exp(cox_death.coef_[0]):.3f}")
主要な概念
モデル選択ガイド
High-dimensional (p > n)?
├── Yes → CoxnetSurvivalAnalysis (elastic net)
└── No
├── Need interpretable coefficients?
│ ├── Yes → CoxPHSurvivalAnalysis or ComponentwiseGB
│ └── No
│ ├── Large data (n > 1000) → GradientBoostingSurvivalAnalysis
│ ├── Medium data → RandomSurvivalForest or FastKernelSurvivalSVM
│ └── Small data → RandomSurvivalForest
└── For max performance → compare multiple models
打ち切りの種類
scikit-survival は主に右打ち切り(研究終了時までにイベントが観測されない場合)を扱います。左打ち切りや区間打ち切りには、別のツール(例:lifelines)が必要です。打ち切り率は指標の選択に影響します。打ち切り率が 40% を超える場合、Harrell の C-index より Uno の C-index が適しています。
一致指数の解釈
- 0.5 = ランダムな予測(識別能力なし)
- 0.6-0.7 = 中程度の識別能力
- 0.7-0.8 = 良好な識別能力(臨床モデルでは一般的)
- 0.8+ = 非常に優れた識別能力(臨床の生存時間データではまれ)
- C-index が測定するのは順位付け能力であり、較正ではありません。必ず Brier スコアも併用してください。
評価指標の選択表
| 指標 | 測定対象 | 最適な用途 |
|---|---|---|
| Harrell's C-index | 順位付けのみ | 迅速な開発、打ち切りが少ない場合 |
| Uno's C-index | 順位付け(打ち切りを補正) | 論文での報告、打ち切り率を問わない場合 |
| Time-dependent AUC | 固定時点 t における識別能 | 特定の予測期間における臨床上の判断 |
| Brier Score (single t) | 固定時点 t における較正 | ある時点での確率の精度 |
| Integrated Brier Score | 全体的な較正 | 包括的なモデル評価 |
一般的なワークフロー
ワークフロー 1:標準的な生存分析パイプライン
目的:前処理、モデルの学習、評価を含むエンドツーエンドの分析です。
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV, train_test_split
from sksurv.datasets import load_breast_cancer
from sksurv.ensemble import GradientBoostingSurvivalAnalysis
from sksurv.metrics import (
concordance_index_ipcw, cumulative_dynamic_auc,
integrated_brier_score, as_concordance_index_ipcw_scorer,
)
import numpy as np
X, y = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Build pipeline
pipe = Pipeline([
("scaler", StandardScaler()),
("model", GradientBoostingSurvivalAnalysis(random_state=42)),
])
# Hyperparameter search
param_grid = {
"model__learning_rate": [0.01, 0.05, 0.1],
"model__n_estimators": [100, 200],
"model__max_depth": [3, 5],
}
cv = GridSearchCV(pipe, param_grid, cv=5,
scoring=as_concordance_index_ipcw_scorer(), n_jobs=-1)
cv.fit(X_train, y_train)
print(f"Best CV C-index: {cv.best_score_:.3f}, Params: {cv.best_params_}")
# Final evaluation on test set
best = cv.best_estimator_
risk = best.predict(X_test)
c_uno = concordance_index_ipcw(y_train, y_test, risk)[0]
times = np.percentile(y_test["time"][y_test["event"]], [25, 50, 75])
auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk, times)
print(f"Test C-index: {c_uno:.3f}, Mean AUC: {mean_auc:.3f}")
ワークフロー 2:モデルの比較
目的:複数のモデルファミリーを比較し、最適なものを選択します。
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
from sklearn.preprocessing import StandardScaler
import numpy as np
scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s = scaler.transform(X_test)
models = {
"CoxPH": CoxPHSurvivalAnalysis(),
"RSF": RandomSurvivalForest(n_estimators=200, random_state=42),
"GBS": GradientBoostingSurvivalAnalysis(n_estimators=200, random_state=42),
"LinearSVM": FastSurvivalSVM(alpha=1.0, random_state=42),
}
times = np.percentile(y_test["time"][y_test["event"]], [25, 50, 75])
for name, model in models.items():
X_tr = X_train_s if name in ("CoxPH", "LinearSVM") else X_train
X_te = X_test_s if name in ("CoxPH", "LinearSVM") else X_test
model.fit(X_tr, y_train)
risk = model.predict(X_te)
c = concordance_index_ipcw(y_train, y_test, risk)[0]
print(f"{name:12s}: C-index = {c:.3f}")
ワークフロー 3:高次元の特徴量選択
- l1_ratio=0.9(Lasso 優位)を指定して CoxnetSurvivalAnalysis を学習します
- C-index スコアラーを使用し、alpha_min_ratio の値を GridSearchCV で探索します
- 最適なモデルからゼロ以外の係数を抽出します
- 選択した特徴量に対して標準的な CoxPH または RSF を再学習し、最終モデルを作成します
- 保留したテストセットで Uno's C-index を使用して評価します
主要なパラメーター
| パラメーター | モジュール | デフォルト | 範囲 | 効果 |
|---|---|---|---|---|
alpha | CoxPH | 0.0 | >= 0 | 正則化(0 = なし) |
l1_ratio | CoxNet | 0.5 | 0-1 | L1 と L2 のペナルティ(1=Lasso、0=Ridge) |
n_estimators | RSF, GBS | 100 | 50-1000 | 木/反復の数 |
max_depth | RSF, GBS | None/3 | 1-20 | 木の深さ(GBS:3-5 を推奨) |
learning_rate | GBS | 0.1 | 0.001-0.3 | 各反復でのステップサイズの縮小率 |
subsample | GBS | 1.0 | 0.5-1.0 | 各木で使用するサンプルの割合 |
dropout_rate | GBS | 0.0 | 0-0.3 | 学習中に前の学習器を除外する割合 |
max_features | RSF | None | "sqrt","log2",int | 分割ごとの特徴量数 |
alpha | FastSurvivalSVM | 1.0 | 0.01-100 | SVM の正則化 |
kernel | FastKernelSVM | "rbf" | "linear","rbf","poly","sigmoid" | カーネル関数 |
gamma | FastKernelSVM | "scale" | "scale","auto",float | カーネル係数 |
ベストプラクティス
-
Cox モデルと SVM モデルでは、必ず特徴量を標準化します:係数とマージンはスケールに依存します。データリークを防ぐため、Pipeline 内で
StandardScalerを使用します。 -
Harrell's C-index より Uno's C-index を使用します:打ち切り分布の違いに対してより頑健で、データセット間でモデルを比較する際に必要です。最初の引数として
y_trainを渡します。 -
アンチパターン -- 組み込みの RSF 特徴量重要度を使用する:連続特徴量やカテゴリ数の多い特徴量を優遇するバイアスがあります。代わりに
sklearn.inspection.permutation_importance()を使用します。 -
複数の指標を報告します:C-index だけでは識別能しか測定できません。較正には Brier score、時点ごとの性能には time-dependent AUC を追加します。
-
アンチパターン -- 比例ハザードの仮定を無視する:Cox モデルでは、ハザード比が時間を通じて一定であると仮定します。Schoenfeld 残差(lifelines で利用可能)を使って検証するか、ノンパラメトリックな代替手法(RSF、GBS)を使用します。
-
特徴量あたりのイベント数を十分に確保します:目安は Cox で特徴量あたり 10 件以上のイベントです。イベント数が少ない場合は、CoxNet(l1_ratio=0.9)またはアンサンブル手法を使用します。
-
scikit-learn の Pipeline を使用します:前処理とモデルをまとめることで、データリークを防ぎ、交差検証やデプロイを簡素化できます。
-
アンチパターン -- KM で競合イベントを打ち切りとして扱う:競合リスクが存在する場合、Kaplan-Meier はイベント確率を過大評価します。代わりに
cumulative_incidence_competing_risks()を使用します。 -
GBS の早期停止:GradientBoosting では、検証性能が頭打ちになったときに自動停止するよう、
n_estimators=1000とn_iter_no_change=10を設定します。
よく使うレシピ
レシピ:交差検証によるハイパーパラメーター調整
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sksurv.ensemble import RandomSurvivalForest
from sksurv.metrics import as_concordance_index_ipcw_scorer
pipe = Pipeline([
("scaler", StandardScaler()),
("model", RandomSurvivalForest(random_state=42)),
])
param_grid = {
"model__n_estimators": [100, 300, 500],
"model__min_samples_split": [6, 10, 20],
"model__max_depth": [None, 10, 20],
}
cv = GridSearchCV(pipe, param_grid, cv=5,
scoring=as_concordance_index_ipcw_scorer(), n_jobs=-1)
cv.fit(X, y)
print(f"Best: {cv.best_score_:.3f}, Params: {cv.best_params_}")
レシピ:置換による特徴量重要度
from sklearn.inspection import permutation_importance
from sksurv.metrics import as_concordance_index_ipcw_scorer
result = permutation_importance(
rsf, X_test, y_test, n_repeats=15,
scoring=as_concordance_index_ipcw_scorer(), random_state=42, n_jobs=-1,
)
sorted_idx = result.importances_mean.argsort()[::-1]
for i in sorted_idx[:10]:
print(f"{X.columns[i]:20s}: {result.importances_mean[i]:.4f} "
f"+/- {result.importances_std[i]:.4f}")
レシピ:信頼帯付き Kaplan-Meier
from sksurv.nonparametric import kaplan_meier_estimator
import matplotlib.pyplot as plt
time, surv, conf_int = kaplan_meier_estimator(
y["event"], y["time"], conf_type="log-log",
)
plt.step(time, surv, where="post", label="KM estimate")
plt.fill_between(time, conf_int[0], conf_int[1], alpha=0.25, step="post")
plt.xlabel("Time")
plt.ylabel("Survival probability")
plt.legend()
plt.savefig("km_with_ci.png", dpi=150)
トラブルシューティング
| 問題 | 原因 | 解決策 |
|---|---|---|
ValueError: y must be a structured array | アウトカム形式が正しくありません | Surv.from_arrays(event, time) を使用してください。event は bool、time は float である必要があります |
concordance_index_ipcw でエラーが発生する | y_train 引数がありません | IPCW では打ち切り分布の算出に学習セットが必要です。最初の引数に y_train を渡してください |
| Cox モデルが収束しない | 特徴量間の相関が高い、または標準化されていません | StandardScaler で標準化してください。CoxnetSurvivalAnalysis への切り替えも検討してください |
| C-index が非常に低い(< 0.5) | 符号が誤っている、または予測性のない特徴量です | event/time のコーディングを確認してください。特徴量エンジニアリングも確認してください |
CV で as_concordance_index_ipcw_scorer がエラーになる | フォールド間で時間範囲が一致していません | 評価対象の時点におけるイベントが、すべてのフォールドに含まれることを確認してください |
| SVM の処理が極端に遅い | カーネル SVM に大規模データセットを使用しています | 大規模な n には FastSurvivalSVM(線形、O(np))を使用してください。カーネルを使用する場合は、サブサンプリングしてください |
| RSF の特徴量重要度が誤解を招く | 組み込みの不純度重要度には偏りがあります | scikit-learn の permutation_importance() を使用してください |
| Kaplan-Meier が過大評価する | 競合イベントが打ち切りとして扱われています | 代わりに cumulative_incidence_competing_risks() を使用してください |
同梱リソース
references/models_evaluation.md
内容:Cox モデルファミリーの詳細(CoxPH、CoxNet、IPCRidge のパラメーター設定の指針と解釈)、アンサンブル手法(RSF、GBS、ComponentwiseGB、ExtraSurvivalTrees のハイパーパラメーター調整、早期停止、損失関数)、SVM モデル(FastSurvivalSVM、FastKernelSurvivalSVM、HingeLossSurvivalSVM、NaiveSurvivalSVM、MinlipSurvivalAnalysis のカーネル選択)、評価指標(Harrell と Uno の C-index の選択、時間依存 AUC のプロット、ヌルモデルとの比較を含む Brier スコア、包括的な評価パイプライン)、モデル比較表。
インラインに再配置:すべてのモデルファミリーの基本的な使用パターンを、Core API のモジュール 2〜4 およびモジュール 6 に集約しました。モデル選択の判断フローチャートを Key Concepts に再配置しました。
省略:スコアリング文字列による CoxNet の交差検証時の alpha 選択(有効な sklearn スコアラーの使用方法ではありません。Core API に示すように as_concordance_index_ipcw_scorer() を使用してください)。NaiveSurvivalSVM の詳細な使用方法(低速で、主にベンチマーク用)。MinlipSurvivalAnalysis の詳細な使用方法(研究用途のみ)。元のリファレンスファイル間で重複していた「使用すべき状況」の説明。
元の 4 ファイルを統合:cox-models.md(182 行)、ensemble-models.md(327 行)、svm-models.md(411 行)、evaluation-metrics.md(378 行)。合計:1,298 行。
references/data_competing_risks.md
内容:データの前処理ワークフロー全般(Surv 配列、カスタム CSV の読み込み、ARFF の入出力、カテゴリーデータのエンコード方法、標準化、欠損データの補完戦略、特徴量選択)、データ検証チェック(負の時間、打ち切り率、特徴量あたりのイベント数)、学習・テスト分割(ランダム分割とイベント/時間で層化した分割)、scikit-learn を使った完全な前処理パイプライン、競合リスク分析(CIF 推定、イベントタイプを含むデータ形式、層化グループ比較、原因別ハザードモデリング、シミュレーションデータを使った完全な実用例)。
インラインに再配置:Surv 配列の基本的な作成方法、encode_categorical、StandardScaler の使用方法を Core API のモジュール 1 に集約しました。CIF 推定と原因別 Cox モデルを Core API のモジュール 7 に再配置しました。
省略:docstring 付きのカスタム前処理関数(ワークフロー 1 で示した sklearn Pipeline のパターンと重複します)。時間依存共変量に関する注記(scikit-survival ではサポートされていないため、元の資料では 1 行の注記としてのみ記載)。Fine-Gray のサブ分布モデル(scikit-survival では利用できず、lifelines の代替手段として記載)。
元の 2 ファイルを統合:data-handling.md(494 行)、competing-risks.md(397 行)。合計:891 行。
関連スキル
- scikit-learn-machine-learning -- 同じ Pipeline/GridSearchCV API を使用する一般的な ML。生存時間解析以外のタスク向け
- statsmodels-statistical-modeling -- 頻度論的回帰(OLS、GLM)。打ち切りのないアウトカム向け
- pymc-bayesian-modeling -- 事前分布を指定するベイズ生存モデル
- matplotlib-scientific-plotting -- Kaplan-Meier 曲線と生存時間プロットのカスタマイズ
- statistical-analysis -- 検定の選択フレームワークと報告ガイドライン
参考資料
- scikit-survival のドキュメント -- 公式 API リファレンス
- scikit-survival の GitHub -- ソースコードと例
- API リファレンス -- クラス/関数の完全な一覧
- Polsterl (2020)「scikit-survival: scikit-learn を基盤とするイベント発生時間解析ライブラリ」 -- JMLR 21(212):1-6