カプラン・マイヤー推定量で2群の生存曲線を描いたとき、「2つの曲線は本当に異なるのか、それとも偶然の変動で離れて見えているだけか」を統計的に判定したい場面が生じます。たとえば、新薬群と対照群の生存曲線の差が統計的に有意かどうかを判定する場面です。
この問いに答える最も標準的な手法がログランク検定(log-rank test)です。ログランク検定は、各イベント時刻において「帰無仮説(2群の生存関数が同一)のもとで期待されるイベント数」と「実際に観測されたイベント数」を比較し、全時点にわたって集計します。
ログランク検定を理解すると、以下のような応用が可能です。
- 臨床試験: 治療群と対照群の生存期間の比較に最も広く使われる検定です
- 信頼性工学: 異なる製造条件での製品寿命の比較に利用されます
- 疫学: 曝露群と非曝露群の疾病発症リスクの比較に使われます
本記事の内容
- ログランク検定の直感的な理解
- 検定統計量の導出
- 重み付きログランク検定の変種
- 多群への拡張
- Pythonによる実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ログランク検定の直感
各時点での「期待値との比較」
ログランク検定のアイデアは、各イベント時刻での超幾何分布に基づいています。
時刻 $t_{(j)}$ でイベントが発生したとき、リスク集合($t_{(j)}$ の直前に生存している全個体)から $d_j$ 人がイベントを経験します。もし帰無仮説(2群の生存関数が同一)が正しければ、このイベントは両群から均等に発生するはずです。
群1のリスク集合のサイズを $n_{1j}$、全体のリスク集合のサイズを $n_j$、全体のイベント数を $d_j$ とすると、帰無仮説のもとで群1の期待イベント数は次のようになります。
$$ e_{1j} = d_j \cdot \frac{n_{1j}}{n_j} $$
これは「リスクにある人数の割合に比例してイベントが配分される」という自然な期待値です。
群1の観測イベント数を $d_{1j}$ とすると、各時点での「観測 – 期待」の差 $d_{1j} – e_{1j}$ が帰無仮説からのずれを測ります。ログランク検定統計量は、この差を全時点にわたって集計したものです。
検定統計量
ログランク検定統計量は次のように定義されます。
$$ \begin{equation} \chi^2_{LR} = \frac{\left(\sum_{j=1}^{m}(d_{1j} – e_{1j})\right)^2}{\sum_{j=1}^{m} v_{1j}} \end{equation} $$
ここで分子は各時点の「観測 – 期待」の総和の2乗、分母は各時点の分散の総和です。
各時点の分散は超幾何分布の分散から次のように計算されます。
$$ \begin{equation} v_{1j} = \frac{n_{1j} \cdot n_{2j} \cdot d_j \cdot (n_j – d_j)}{n_j^2(n_j – 1)} \end{equation} $$
帰無仮説のもとで $\chi^2_{LR} \xrightarrow{d} \chi^2(1)$ です。$k$ 群の場合は $\chi^2(k-1)$ です。
超幾何分布に基づく正確な導出
ログランク検定統計量の導出をより詳しく見ていきましょう。各イベント時刻 $t_{(j)}$ において、帰無仮説(2群のハザード関数が同一)のもとで、群1のイベント数 $d_{1j}$ がどのような分布に従うかを考えます。
時刻 $t_{(j)}$ の直前にリスク集合にいる人数は、群1が $n_{1j}$、群2が $n_{2j}$、合計 $n_j = n_{1j} + n_{2j}$ です。この時刻に $d_j$ 人がイベントを経験します。帰無仮説のもとでは、両群のハザード関数は同じなので、イベントは両群から無差別に発生します。
これは「$n_j$ 個の中から $d_j$ 個を非復元抽出するとき、$n_{1j}$ 個ある群1のものがいくつ含まれるか」という問題と同じ構造です。したがって、$d_{1j}$ は超幾何分布に従います。
$$ d_{1j} \mid n_{1j}, n_{2j}, d_j \sim \text{Hypergeometric}(n_j, n_{1j}, d_j) $$
超幾何分布の期待値と分散は以下の通りです。
$$ E[d_{1j}] = e_{1j} = d_j \cdot \frac{n_{1j}}{n_j} $$
$$ \text{Var}(d_{1j}) = v_{1j} = \frac{n_{1j} \cdot n_{2j} \cdot d_j \cdot (n_j – d_j)}{n_j^2 \cdot (n_j – 1)} $$
分散の式を直感的に理解しましょう。分子の $n_{1j} \cdot n_{2j}$ は2群のリスク集合サイズの積で、2群のサイズが均等に近いほど分散は大きくなります。$d_j \cdot (n_j – d_j)$ はイベント数と非イベント数の積で、イベントがリスク集合に対して中程度の割合のときに分散が最大になります。分母の $(n_j – 1)$ は有限母集団補正に対応します。
検定統計量の漸近的な正当化
「観測 – 期待」の累積和 $U = \sum_{j=1}^m (d_{1j} – e_{1j})$ は、帰無仮説のもとで期待値0です。各時点のイベントが(条件付きで)独立なので、分散は各時点の分散の和 $V = \sum_{j=1}^m v_{1j}$ です。
中心極限定理の一般化により、$m$(イベント時刻の数)が十分大きいとき、
$$ \frac{U}{\sqrt{V}} \xrightarrow{d} N(0, 1) \quad \text{under } H_0 $$
したがって $U^2/V \xrightarrow{d} \chi^2(1)$ であり、これがログランク検定統計量 $\chi^2_{LR}$ の漸近分布です。
ここで注意が必要なのは、各時点の超幾何分布は条件付き独立である(リスク集合に条件づけた上での独立性)という点です。厳密にはマルチンゲール中心極限定理を使って正当化します。
ログランク検定とCox回帰の関係
ログランク検定統計量は、実はCox比例ハザードモデルのスコア検定統計量と一致します。
Cox回帰モデルで群の指示変数 $Z_i$($Z_i = 1$ なら群1、$Z_i = 0$ なら群2)を共変量として、$h(t \mid Z_i) = h_0(t) \exp(\beta Z_i)$ とモデル化します。帰無仮説 $H_0: \beta = 0$ に対するスコア検定統計量は、$\beta = 0$ における偏対数尤度のスコアの二乗をフィッシャー情報量で割ったものです。
Cox偏尤度のスコアは、
$$ S(\beta)\big|_{\beta=0} = \sum_{j=1}^m (d_{1j} – e_{1j}) $$
これはまさにログランク検定統計量の分子の $U$ です。フィッシャー情報量は $V = \sum v_{1j}$ に等しく、したがってスコア検定統計量 $U^2/V$ はログランク検定統計量と一致します。
この関係は、ログランク検定が比例ハザードの対立仮説に対して最も検出力が高いことの理論的な裏付けを与えます。スコア検定が局所的に最良な性質を持つことから、ログランク検定は比例ハザード型の対立仮説に対してローカルに最強力な検定です。
ログランク検定の理論的背景がわかったところで、重み付き検定への拡張を見ていきましょう。
重み付きログランク検定
標準的なログランク検定はすべての時点に等しい重みを置きますが、異なる重みを使う変種があります。
| 検定 | 重み $w_j$ | 特徴 |
|---|---|---|
| ログランク | 1 | 比例ハザードの対立仮説に最適 |
| Gehan-Wilcoxon | $n_j$ | 早期の差に敏感 |
| Peto-Prentice | $\tilde{S}(t_j)$ | 中間的 |
| Fleming-Harrington | $[\hat{S}(t)]^p[1-\hat{S}(t)]^q$ | 柔軟な重み |
重み付き検定統計量は次のようになります。
$$ \chi^2_W = \frac{\left(\sum_{j} w_j(d_{1j} – e_{1j})\right)^2}{\sum_{j} w_j^2 v_{1j}} $$
Gehan-Wilcoxon検定はリスク集合のサイズ $n_j$ で重み付けするため、観測者数が多い早期の時点に大きな重みを置きます。早期にリスク集合が大きく、時間が経つにつれて打ち切りやイベントによりリスク集合が小さくなるのが典型的なので、Gehan-Wilcoxon検定は実質的に早期の差を重視します。
一方、ログランク検定は後期の差にも等しい重みを置くため、比例ハザードの対立仮説のもとで最も検出力が高い検定です。
Fleming-Harrington検定は $[\hat{S}(t)]^p[1-\hat{S}(t)]^q$ という形の重みを使い、$p, q$ の値によって早期・後期のどちらを重視するかを調整できます。$p = 0, q = 0$ のときログランク検定、$p = 1, q = 0$ のときは早期重視、$p = 0, q = 1$ のときは後期重視に対応します。
検定の選択指針
どの重み付き検定を使うべきかは、対立仮説の形に依存します。
- 2群のハザード比が時間を通じて一定であると想定される場合(比例ハザード)→ ログランク検定が最適
- 2群の差が早期に顕著であると想定される場合 → Gehan-Wilcoxon検定が適切
- 2群の差が後期に顕著であると想定される場合(例: 治療効果の遅延発現)→ Fleming-Harrington ($p = 0, q = 1$) が適切
- 対立仮説の形が不明な場合 → 複数の検定を事前に指定して実施するか、ログランク検定をデフォルトとする
重要なのは、検定の選択はデータを見る前に行うべきだという点です。データを見てから最も有意な結果を出す検定を選ぶと、検定の有意水準が名目値より大きくなってしまいます(多重比較の問題)。
ここまでで2群の比較を扱ってきましたが、3群以上の比較に拡張する方法を見ていきましょう。
多群への拡張
$k$ 群のログランク検定
$k$ 群の比較の場合、ログランク検定は自然に拡張されます。帰無仮説は $H_0: S_1(t) = S_2(t) = \cdots = S_k(t)$ です。
各イベント時刻 $t_{(j)}$ において、群 $g$ ($g = 1, \dots, k$) の観測イベント数 $d_{gj}$ と期待イベント数 $e_{gj} = d_j \cdot n_{gj} / n_j$ を計算します。
「観測 – 期待」のベクトル $\bm{U} = (U_1, \dots, U_{k-1})^\top$ を構成します($U_g = \sum_j (d_{gj} – e_{gj})$)。最後の群は $\sum U_g = 0$ の制約から冗長なので、$k-1$ 個だけ取ります。
$\bm{U}$ の共分散行列 $\bm{V}$ の $(g, h)$ 成分は次のように計算されます。
$$ V_{gh} = \sum_{j=1}^m v_{ghj} $$
$g = h$ のとき: $v_{ggj} = \frac{n_{gj}(n_j – n_{gj})d_j(n_j – d_j)}{n_j^2(n_j – 1)}$
$g \neq h$ のとき: $v_{ghj} = -\frac{n_{gj} \cdot n_{hj} \cdot d_j(n_j – d_j)}{n_j^2(n_j – 1)}$
多群のログランク検定統計量は、
$$ \begin{equation} \chi^2_{LR} = \bm{U}^\top \bm{V}^{-1} \bm{U} \xrightarrow{d} \chi^2(k-1) \quad \text{under } H_0 \end{equation} $$
層化ログランク検定
交絡因子(たとえば年齢層、施設など)の影響を制御するために、層化ログランク検定(stratified log-rank test)が使われます。
データを $L$ 個の層に分け、各層内でログランク統計量を計算し、それらを合計します。
$$ U_{\text{strat}} = \sum_{\ell=1}^L U_\ell, \quad V_{\text{strat}} = \sum_{\ell=1}^L V_\ell $$
$$ \chi^2_{\text{strat}} = \frac{U_{\text{strat}}^2}{V_{\text{strat}}} $$
層化検定は、各層内では帰無仮説が成立する(つまり層内で群間に差がない)ことを検定します。交絡因子の影響を除去した上で処理効果を評価できるため、臨床試験のプライマリーエンドポイント解析でしばしば使用されます。
検出力の解析
サンプルサイズの計算
臨床試験の設計においては、必要なサンプルサイズを事前に計算することが不可欠です。比例ハザードの対立仮説 $H_1: \text{ハザード比} = \Delta$ のもとで、ログランク検定の検出力を確保するために必要なイベント数 $D$ は、
$$ \begin{equation} D = \frac{(z_{1-\alpha/2} + z_{1-\beta})^2}{(\ln \Delta)^2} \cdot \frac{(1 + r)^2}{r} \end{equation} $$
ここで $z_p$ は標準正規分布の上側 $p$ 分位点、$\alpha$ は有意水準、$\beta$ は第二種の過誤確率($1 – \beta$ が検出力)、$r = n_1 / n_2$ はサンプルサイズ比です。
たとえば、ハザード比 $\Delta = 0.7$(治療群のハザードが対照群の70%)を検出したい場合、$\alpha = 0.05$(両側)、検出力 $1 – \beta = 0.80$、$r = 1$(等サイズ)として、
$$ D = \frac{(1.96 + 0.84)^2}{(\ln 0.7)^2} \cdot \frac{4}{1} = \frac{7.84}{0.1278} \cdot 4 \approx 245 $$
245件のイベントが必要です。総サンプルサイズは、予想イベント率に基づいて計算されます。
この検出力の計算は、ログランク検定がスコア検定として比例ハザードの対立仮説に対して最適であることを前提としています。
Pythonによる実装
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
def log_rank_test(times1, events1, times2, events2, weight="logrank"):
"""
ログランク検定(2群)
Parameters
----------
times1, events1 : 群1のデータ
times2, events2 : 群2のデータ
weight : str 重み ('logrank', 'wilcoxon', 'peto')
"""
# 全イベント時刻
all_event_times = np.unique(np.concatenate([
times1[events1 == 1], times2[events2 == 1]
]))
all_event_times.sort()
O_E_sum = 0.0
var_sum = 0.0
S_km = 1.0 # Peto用
for t in all_event_times:
# リスク集合
n1 = np.sum(times1 >= t)
n2 = np.sum(times2 >= t)
n = n1 + n2
# イベント数
d1 = np.sum((times1 == t) & (events1 == 1))
d2 = np.sum((times2 == t) & (events2 == 1))
d = d1 + d2
if n <= 1:
continue
# 期待値
e1 = d * n1 / n
# 分散
v1 = n1 * n2 * d * (n - d) / (n**2 * (n - 1))
# 重み
if weight == "logrank":
w = 1.0
elif weight == "wilcoxon":
w = n
elif weight == "peto":
w = S_km
else:
w = 1.0
O_E_sum += w * (d1 - e1)
var_sum += w**2 * v1
# KM推定量の更新(Peto用)
S_km *= (1 - d / n)
if var_sum <= 0:
return 0, 1.0
chi2_stat = O_E_sum**2 / var_sum
p_value = 1 - stats.chi2.cdf(chi2_stat, df=1)
return chi2_stat, p_value
# --- テストデータ ---
np.random.seed(42)
n1, n2 = 60, 60
t1 = np.random.weibull(1.5, n1) * 30
c1 = np.random.uniform(15, 50, n1)
obs1 = np.minimum(t1, c1)
ev1 = (t1 <= c1).astype(int)
t2 = np.random.weibull(1.5, n2) * 20
c2 = np.random.uniform(15, 50, n2)
obs2 = np.minimum(t2, c2)
ev2 = (t2 <= c2).astype(int)
# 3つの重みでの検定
for w_name in ["logrank", "wilcoxon", "peto"]:
chi2, p = log_rank_test(obs1, ev1, obs2, ev2, weight=w_name)
print(f"{w_name:12}: χ² = {chi2:.3f}, p = {p:.6f}")
# --- カプラン・マイヤー曲線と検定結果の可視化 ---
from itertools import accumulate
def km_curve(times, events):
t_sorted = np.sort(np.unique(times[events == 1]))
km_t = [0.0]
km_s = [1.0]
s = 1.0
for t in t_sorted:
n_risk = np.sum(times >= t)
d = np.sum((times == t) & (events == 1))
if n_risk > 0:
s *= (1 - d / n_risk)
km_t.append(t)
km_s.append(s)
return np.array(km_t), np.array(km_s)
chi2_lr, p_lr = log_rank_test(obs1, ev1, obs2, ev2)
fig, ax = plt.subplots(figsize=(9, 6))
t1_km, s1_km = km_curve(obs1, ev1)
t2_km, s2_km = km_curve(obs2, ev2)
ax.step(t1_km, s1_km, where="post", linewidth=2, color="blue", label="Group 1")
ax.step(t2_km, s2_km, where="post", linewidth=2, color="red", label="Group 2")
ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Survival probability", fontsize=12)
ax.set_title(f"Log-rank test: χ² = {chi2_lr:.2f}, p = {p_lr:.4f}", fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)
plt.tight_layout()
plt.savefig("log_rank_test.png", dpi=150, bbox_inches="tight")
plt.show()
この結果から、ログランク検定が2群の生存曲線の差を統計的に評価できることが確認できます。p値が小さければ2群の生存関数が異なることを示し、p値が大きければ差があるとは言えません。
3つの重み付き検定の結果を比べると、重みの選択によってp値が異なることがわかります。データの構造(差が早期に現れるか後期に現れるか)によって、各検定の感度が変わるためです。カプラン・マイヤー曲線のグラフから視覚的に2群の差が確認でき、その差の統計的有意性がログランク検定のp値によって定量的に評価されています。
検出力のシミュレーション
異なるハザード比に対するログランク検定の検出力を、シミュレーションにより確認します。
import numpy as np
from scipy import stats
def log_rank_statistic(t1, e1, t2, e2):
"""ログランク検定統計量を計算(高速版)"""
event_times = np.unique(np.concatenate([t1[e1==1], t2[e2==1]]))
event_times.sort()
U, V = 0.0, 0.0
for t in event_times:
n1 = np.sum(t1 >= t)
n2 = np.sum(t2 >= t)
n = n1 + n2
d1 = np.sum((t1 == t) & (e1 == 1))
d2 = np.sum((t2 == t) & (e2 == 1))
d = d1 + d2
if n <= 1:
continue
e = d * n1 / n
v = n1 * n2 * d * (n - d) / (n**2 * (n - 1))
U += d1 - e
V += v
return U**2 / V if V > 0 else 0
np.random.seed(42)
n_sim = 1000
alpha = 0.05
sample_sizes = [30, 60, 100, 200]
hazard_ratios = np.linspace(0.3, 1.0, 15)
print("=== ログランク検定の検出力(シミュレーション) ===")
print(f"{'HR':>6} | " + " | ".join([f"n={n:>3}" for n in sample_sizes]))
print("-" * 50)
for hr in [0.4, 0.6, 0.8, 1.0]:
powers = []
for n in sample_sizes:
rejections = 0
for _ in range(n_sim):
t1 = np.random.exponential(1.0, n)
c1 = np.random.uniform(0, 5, n)
obs1 = np.minimum(t1, c1)
ev1 = (t1 <= c1).astype(int)
t2 = np.random.exponential(1.0 / hr, n)
c2 = np.random.uniform(0, 5, n)
obs2 = np.minimum(t2, c2)
ev2 = (t2 <= c2).astype(int)
stat = log_rank_statistic(obs1, ev1, obs2, ev2)
if stat > stats.chi2.ppf(1 - alpha, 1):
rejections += 1
powers.append(rejections / n_sim)
print(f"{hr:>6.2f} | " + " | ".join([f"{p:>5.3f}" for p in powers]))
この結果から、ログランク検定の検出力に関する以下の知見が得られます。
-
ハザード比が1に近いほど検出力が低い: $HR = 1.0$(帰無仮説が正しい場合)では、棄却率は名目有意水準 $\alpha = 0.05$ に近くなるはずです。$HR$ が1から離れるほど検出力が上がります
-
サンプルサイズが大きいほど検出力が高い: 同じハザード比に対して、$n$ が大きいほど差を検出しやすくなります
-
ハザード比0.6を検出するには $n \approx 100$ が必要: 実用的な検出力(80%以上)を達成するために必要なサンプルサイズが確認できます
まとめ
本記事では、ログランク検定の理論を基礎から体系的に解説しました。
- ログランク検定は各イベント時刻での超幾何分布に基づいて「観測イベント数 – 期待イベント数」を全時点にわたって集計し、$\chi^2$ 検定を行います
- 帰無仮説は2群(または多群)の生存関数(ハザード関数)が同一であることです
- Cox回帰との関係: ログランク検定統計量はCox比例ハザードモデルのスコア検定統計量と一致し、比例ハザードの対立仮説に対してローカルに最強力です
- 重み付き変種: Gehan-Wilcoxon(早期重視)、Peto-Prentice(中間的)、Fleming-Harrington(柔軟な重み)があり、対立仮説の形に応じて選択します
- 多群への拡張: $k$ 群の場合は $(k-1)$ 次元のスコアベクトルの二次形式として構成し、$\chi^2(k-1)$ 分布に従います
- 層化検定: 交絡因子の影響を制御するために層化ログランク検定が使われます
- サンプルサイズ計算: 必要なイベント数はハザード比の対数の二乗の逆数に比例します
ログランク検定は生存時間分析において最も広く使われる検定であり、臨床試験のプライマリーエンドポイント解析の標準的な手法です。比例ハザードの仮定が成り立つ場合に最も強力ですが、仮定が崩れる場合は重み付き変種を検討する必要があります。
次のステップとして、以下の記事も参考にしてください。
- Cox比例ハザードモデル — 共変量を含む生存時間回帰
- 生存時間分析入門 — カプラン・マイヤー推定量の詳細
- 偽発見率の制御 — 多重検定問題への対処法
- スコア検定 — ログランク検定の理論的背景