カプラン・マイヤー推定量では生存曲線を群別に推定できましたが、「年齢が10歳上がるとリスクは何倍になるか」や「喫煙者と非喫煙者で、他の条件を揃えた上で生存時間に差があるか」といった問いには答えられません。共変量（年齢、性別、治療法など）が生存時間に与える影響を定量的に評価するには、回帰モデルが必要です。
Cox比例ハザードモデル（Cox proportional hazards model）は、1972年にデイビッド・コックス（David Cox）が提案した半パラメトリックな生存時間回帰モデルです。このモデルの革新的なアイデアは、ベースラインハザード関数を指定しないまま共変量の効果を推定できることです。
Cox回帰モデルを理解すると、以下のような応用が可能です。
- 臨床試験: 交絡因子を調整した上での治療効果の評価に最も広く使われています
- 疫学: リスク因子と疾病発症の関連をハザード比として定量化します
- 信頼性工学: 使用条件（温度、負荷など）が故障リスクに与える影響を評価します
- 経済学: 失業期間に影響する要因の分析に利用されます
本記事の内容
- 比例ハザードの仮定の直感的な意味
- 部分尤度（partial likelihood）の導出と推定手順
- ハザード比の解釈
- 比例ハザードの仮定の検証
- Pythonによる実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
比例ハザードの仮定
モデルの定義
Cox比例ハザードモデルは、共変量ベクトル $\bm{x} = (x_1, x_2, \dots, x_p)^T$ を持つ個体のハザード関数を次のように表します。
$$ \begin{equation} h(t | \bm{x}) = h_0(t) \exp(\bm{\beta}^T \bm{x}) = h_0(t) \exp(\beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p) \end{equation} $$
- $h_0(t)$: ベースラインハザード関数。全共変量がゼロの個体のハザード。形状は未指定（ノンパラメトリック）
- $\exp(\bm{\beta}^T \bm{x})$: リスクスコア。共変量による倍率。常に正
- $\bm{\beta}$: 回帰係数ベクトル。推定の対象
比例ハザードの意味
モデルの名前「比例ハザード」は、2人の個体のハザード関数の比が時間に依存しないことに由来します。
共変量が $\bm{x}_1$ と $\bm{x}_2$ の2人を比較すると、
$$ \frac{h(t | \bm{x}_1)}{h(t | \bm{x}_2)} = \frac{h_0(t)\exp(\bm{\beta}^T \bm{x}_1)}{h_0(t)\exp(\bm{\beta}^T \bm{x}_2)} = \exp[\bm{\beta}^T(\bm{x}_1 – \bm{x}_2)] $$
ベースラインハザード $h_0(t)$ が約分されるため、ハザード比は時間 $t$ に依存しません。つまり、2人のリスクの比は時間を通じて一定です。
たとえば、年齢が10歳高い人のハザード比が1.3であれば、どの時点でもリスクが1.3倍であることを意味します。
ハザード比の解釈
共変量 $x_k$ が1単位増加したときのハザード比は $\exp(\beta_k)$ です。
$$ \begin{equation} \text{HR}_k = \exp(\beta_k) \end{equation} $$
- $\beta_k > 0$（HR $> 1$）: $x_k$ の増加がリスクを高める
- $\beta_k < 0$（HR $< 1$）: $x_k$ の増加がリスクを下げる（保護的）
- $\beta_k = 0$（HR $= 1$）: $x_k$ はリスクに影響しない
この解釈はロジスティック回帰のオッズ比の解釈と類似しています。
部分尤度の導出
コックスの革新的アイデア
Cox回帰の最も革新的な点は、ベースラインハザード $h_0(t)$ を推定せずに $\bm{\beta}$ を推定する方法を見出したことです。この方法が部分尤度（partial likelihood）です。
アイデアは次の通りです。イベントが時刻 $t_{(j)}$ で発生したとき、「リスク集合（$t_{(j)}$ の直前に生存している全個体の集合）$\mathcal{R}_j$ の中で、実際にイベントを経験したのがどの個体か」という条件付き確率を考えます。
時刻 $t_{(j)}$ でイベントが起きた個体の共変量を $\bm{x}_{(j)}$ とすると、リスク集合の中からこの個体が選ばれる条件付き確率は次のようになります。
$$ \frac{h(t_{(j)} | \bm{x}_{(j)})}{\sum_{l \in \mathcal{R}_j} h(t_{(j)} | \bm{x}_l)} = \frac{h_0(t_{(j)})\exp(\bm{\beta}^T \bm{x}_{(j)})}{\sum_{l \in \mathcal{R}_j} h_0(t_{(j)})\exp(\bm{\beta}^T \bm{x}_l)} $$
ベースラインハザード $h_0(t_{(j)})$ が分子・分母で約分されます。
$$ \begin{equation} = \frac{\exp(\bm{\beta}^T \bm{x}_{(j)})}{\sum_{l \in \mathcal{R}_j} \exp(\bm{\beta}^T \bm{x}_l)} \end{equation} $$
部分尤度関数
全イベント時刻について上の条件付き確率を掛け合わせたものが部分尤度関数です。
$$ \begin{equation} L_P(\bm{\beta}) = \prod_{j=1}^{m} \frac{\exp(\bm{\beta}^T \bm{x}_{(j)})}{\sum_{l \in \mathcal{R}_j} \exp(\bm{\beta}^T \bm{x}_l)} \end{equation} $$
対数部分尤度は次のようになります。
$$ \begin{equation} \ell_P(\bm{\beta}) = \sum_{j=1}^{m} \left[\bm{\beta}^T \bm{x}_{(j)} – \ln\left(\sum_{l \in \mathcal{R}_j} \exp(\bm{\beta}^T \bm{x}_l)\right)\right] \end{equation} $$
この対数部分尤度を $\bm{\beta}$ について最大化することで、$\hat{\bm{\beta}}$ が得られます。
推定の性質
部分尤度推定量は完全な尤度の最尤推定量ではありませんが、次の望ましい性質を持つことがコックスによって示されています。
- 一致性: $n \to \infty$ で $\hat{\bm{\beta}} \to \bm{\beta}_0$（真のパラメータ値）
- 漸近正規性: $\sqrt{n}(\hat{\bm{\beta}} – \bm{\beta}_0) \xrightarrow{d} N(\bm{0}, \bm{I}^{-1})$
- 漸近効率性: セミパラメトリックな効率の意味で最適
スコア方程式と情報行列
推定量 $\hat{\bm{\beta}}$ は対数部分尤度のスコア方程式を解いて得られます。
$$ \bm{U}(\bm{\beta}) = \frac{\partial \ell_P}{\partial \bm{\beta}} = \sum_{j=1}^m \left[\bm{x}_{(j)} – \bar{\bm{x}}_j(\bm{\beta})\right] = \bm{0} $$
ここで $\bar{\bm{x}}_j(\bm{\beta})$ はリスク集合における共変量の重み付き平均です。
$$ \bar{\bm{x}}_j(\bm{\beta}) = \frac{\sum_{l \in \mathcal{R}_j} \bm{x}_l \exp(\bm{\beta}^T \bm{x}_l)}{\sum_{l \in \mathcal{R}_j} \exp(\bm{\beta}^T \bm{x}_l)} $$
スコア方程式の意味は明確です。各イベント時刻 $t_{(j)}$ において、イベントを経験した個体の共変量 $\bm{x}_{(j)}$ と、リスク集合の「リスクスコアで重み付けした平均的な共変量」$\bar{\bm{x}}_j$ の差を全時点にわたって合計し、それがゼロになるように $\bm{\beta}$ を選びます。
情報行列は次のように計算されます。
$$ \bm{I}(\bm{\beta}) = -\frac{\partial^2 \ell_P}{\partial \bm{\beta} \partial \bm{\beta}^T} = \sum_{j=1}^m \bm{V}_j(\bm{\beta}) $$
ここで $\bm{V}_j(\bm{\beta})$ はリスク集合における共変量の重み付き分散行列です。
$$ \bm{V}_j(\bm{\beta}) = \frac{\sum_{l \in \mathcal{R}_j} (\bm{x}_l – \bar{\bm{x}}_j)(\bm{x}_l – \bar{\bm{x}}_j)^T \exp(\bm{\beta}^T \bm{x}_l)}{\sum_{l \in \mathcal{R}_j} \exp(\bm{\beta}^T \bm{x}_l)} $$
$\hat{\bm{\beta}}$ の標準誤差は $\bm{I}(\hat{\bm{\beta}})$ の逆行列の対角要素の平方根です。
同順イベントの処理（タイ処理）
実データでは、2人以上の個体が同じ時刻にイベントを経験する（同順、tie）場合があります。正確な部分尤度の計算は同順イベントが多いと組合せ的に計算量が爆発するため、近似法が使われます。
Breslow近似（最も広く使われる）:
$$ L_P(\bm{\beta}) \approx \prod_{j=1}^m \frac{\exp(\bm{\beta}^T \bm{s}_j)}{\left[\sum_{l \in \mathcal{R}_j} \exp(\bm{\beta}^T \bm{x}_l)\right]^{d_j}} $$
ここで $\bm{s}_j = \sum_{k: \text{イベントが} t_{(j)} \text{で発生}} \bm{x}_k$ であり、$d_j$ は時刻 $t_{(j)}$ でのイベント数です。
Efron近似（Breslowより正確）:
$$ L_P(\bm{\beta}) \approx \prod_{j=1}^m \frac{\exp(\bm{\beta}^T \bm{s}_j)}{\prod_{r=1}^{d_j}\left[\sum_{l \in \mathcal{R}_j}\exp(\bm{\beta}^T\bm{x}_l) – \frac{r-1}{d_j}\sum_{k \in \mathcal{D}_j}\exp(\bm{\beta}^T\bm{x}_k)\right]} $$
ここで $\mathcal{D}_j$ は時刻 $t_{(j)}$ でイベントが起きた個体の集合です。
同順イベントが少ない場合は両近似の差は小さいですが、同順が多い場合はEfron近似の方が正確です。
ベースラインハザード関数の推定
ブレスロー推定量
回帰係数 $\hat{\bm{\beta}}$ が得られた後、ベースラインハザード $h_0(t)$ を推定することもできます。ブレスロー（Breslow）推定量は、各イベント時刻 $t_{(j)}$ でのベースラインハザードの累積を推定します。
$$ \begin{equation} \hat{H}_0(t) = \sum_{j: t_{(j)} \leq t} \frac{d_j}{\sum_{l \in \mathcal{R}_j} \exp(\hat{\bm{\beta}}^T \bm{x}_l)} \end{equation} $$
これはカプラン・マイヤー推定量の回帰モデル版と考えることができます。ベースラインの生存関数は、
$$ \hat{S}_0(t) = \exp(-\hat{H}_0(t)) $$
共変量 $\bm{x}$ を持つ個体の生存関数の推定は、
$$ \hat{S}(t | \bm{x}) = [\hat{S}_0(t)]^{\exp(\hat{\bm{\beta}}^T \bm{x})} $$
これにより、個体ごとの生存確率の予測が可能になります。
比例ハザードの仮定の検証
なぜ検証が必要か
Cox回帰の推定結果の妥当性は、比例ハザードの仮定（ハザード比が時間に依存しない）に依存します。この仮定が破れていると、推定されたハザード比の解釈が曖昧になります。仮定の検証は、モデル診断の重要なステップです。
ショーンフェルド残差による検定
ショーンフェルド残差（Schoenfeld residuals）は、各イベント時刻 $t_{(j)}$ に対して次のように定義されます。
$$ \hat{r}_{jk} = x_{(j)k} – \bar{x}_{jk}(\hat{\bm{\beta}}) $$
ここで $x_{(j)k}$ はイベントを経験した個体の第 $k$ 共変量、$\bar{x}_{jk}$ はリスク集合の重み付き平均の第 $k$ 共変量です。
比例ハザードの仮定が正しければ、ショーンフェルド残差は時間に対して無相関のはずです。残差を時間に対してプロットし、系統的なパターンがないかを確認します。
グレンビル・テルノー検定（Grambsch-Therneau test）は、スケーリングしたショーンフェルド残差と時間（または時間の関数）の相関を検定します。帰無仮説は「ショーンフェルド残差と時間が無相関」（＝比例ハザードが成立）です。
比例ハザードの仮定が破れた場合
仮定が破れた場合の対処法には以下があります。
- 時間依存共変量の導入: $\bm{\beta}(t) = \bm{\beta} + \bm{\gamma} \cdot g(t)$ のように、係数に時間の関数を加える
- 層化Cox回帰: 比例ハザードを満たさない共変量（例: 施設）で層化し、層内でのみ比例ハザードを仮定する
- 加速故障時間モデル: 比例ハザードではなく、生存時間のスケールパラメータが共変量に依存するモデルへの切り替え
これらの理論的背景を踏まえて、Pythonでの実装を行いましょう。
Pythonによる実装
Cox回帰のスクラッチ実装
import numpy as np
from scipy.optimize import minimize
from scipy import stats
import matplotlib.pyplot as plt
def cox_regression(times, events, X):
"""
Cox比例ハザードモデルの部分尤度による推定
Parameters
----------
times : array (n,) 観測時間
events : array (n,) イベント指示変数
X : array (n, p) 共変量行列
Returns
-------
result : dict 推定結果
"""
n, p = X.shape
def neg_partial_loglik(beta):
"""負の対数部分尤度"""
risk_scores = X @ beta
loglik = 0.0
for i in range(n):
if events[i] == 1:
risk_set = times >= times[i]
loglik += risk_scores[i] - np.log(
np.sum(np.exp(risk_scores[risk_set]))
)
return -loglik
def gradient(beta):
"""対数部分尤度の勾配"""
risk_scores = X @ beta
grad = np.zeros(p)
for i in range(n):
if events[i] == 1:
risk_set = times >= times[i]
exp_scores = np.exp(risk_scores[risk_set])
X_risk = X[risk_set]
weighted_mean = (X_risk.T @ exp_scores) / np.sum(exp_scores)
grad += X[i] - weighted_mean
return -grad
# 最適化
result = minimize(neg_partial_loglik, x0=np.zeros(p),
jac=gradient, method="BFGS")
beta_hat = result.x
# 標準誤差（ヘッセ行列の逆行列の対角要素の平方根）
hessian = result.hess_inv
se = np.sqrt(np.diag(hessian))
# ハザード比と信頼区間
hr = np.exp(beta_hat)
hr_lower = np.exp(beta_hat - 1.96 * se)
hr_upper = np.exp(beta_hat + 1.96 * se)
# ワルド検定のp値
z = beta_hat / se
p_values = 2 * (1 - stats.norm.cdf(np.abs(z)))
return {
"beta": beta_hat, "se": se,
"hr": hr, "hr_lower": hr_lower, "hr_upper": hr_upper,
"z": z, "p_values": p_values,
"loglik": -result.fun
}
# --- テストデータの生成 ---
np.random.seed(42)
n = 200
# 共変量
age = np.random.normal(60, 10, n)
treatment = np.random.binomial(1, 0.5, n)
biomarker = np.random.normal(0, 1, n)
# 真のパラメータ
beta_true = np.array([0.03, -0.5, 0.3])
# ワイブル分布でイベント時間を生成
linear_pred = beta_true[0] * age + beta_true[1] * treatment + beta_true[2] * biomarker
scale = np.exp(-linear_pred / 1.5) * 20
event_times = np.random.weibull(1.5, n) * scale
# 打ち切り
censor_times = np.random.uniform(10, 60, n)
observed_times = np.minimum(event_times, censor_times)
events = (event_times <= censor_times).astype(int)
# 共変量行列
X = np.column_stack([age, treatment, biomarker])
var_names = ["Age", "Treatment", "Biomarker"]
# Cox回帰
result = cox_regression(observed_times, events, X)
print("=" * 65)
print("Cox比例ハザードモデルの推定結果")
print("=" * 65)
print(f"\nイベント数: {np.sum(events)}, 打ち切り数: {n - np.sum(events)}")
print(f"\n{'変数':<12} {'β':>8} {'SE':>8} {'HR':>8} {'95% CI':>18} {'p値':>10}")
print("-" * 65)
for i, name in enumerate(var_names):
print(f"{name:<12} {result['beta'][i]:>8.4f} {result['se'][i]:>8.4f} "
f"{result['hr'][i]:>8.3f} "
f"({result['hr_lower'][i]:.3f}, {result['hr_upper'][i]:.3f}) "
f"{result['p_values'][i]:>10.4f}")
# --- ハザード比のフォレストプロット ---
fig, ax = plt.subplots(figsize=(8, 4))
y_pos = np.arange(len(var_names))
ax.errorbar(result["hr"], y_pos,
xerr=[result["hr"] - result["hr_lower"],
result["hr_upper"] - result["hr"]],
fmt="o", color="steelblue", markersize=8, capsize=5,
linewidth=2, elinewidth=2)
ax.axvline(1.0, color="red", linestyle="--", linewidth=1)
ax.set_yticks(y_pos)
ax.set_yticklabels(var_names, fontsize=11)
ax.set_xlabel("Hazard Ratio (95% CI)", fontsize=12)
ax.set_title("Cox Proportional Hazards Model", fontsize=13)
ax.grid(True, alpha=0.3, axis="x")
plt.tight_layout()
plt.savefig("cox_forest_plot.png", dpi=150, bbox_inches="tight")
plt.show()
この結果から、以下のことが確認できます。
-
推定されたハザード比が真のパラメータと整合している。年齢のHR $> 1$（加齢がリスクを増加）、治療のHR $< 1$（治療が保護的）、バイオマーカーのHR $> 1$（高値がリスクを増加）という方向性が正しく推定されています
-
フォレストプロットはハザード比と95%信頼区間を視覚的に示している。赤い破線（HR = 1）を信頼区間がまたぐかどうかで、統計的有意性が判断できます
まとめ
本記事では、Cox比例ハザードモデルの理論を部分尤度の導出から実装まで解説しました。
- 比例ハザードの仮定は、共変量の効果がハザードの乗法的な倍率として表され、その倍率が時間に依存しないことです
- 部分尤度はベースラインハザードを指定せずに回帰係数を推定する革新的な方法であり、条件付き確率の積として構成されます
- ハザード比 $\exp(\beta_k)$ は共変量 $x_k$ が1単位増加したときのリスクの倍率を表し、直感的に解釈可能です
- 比例ハザードの仮定の検証にはショーンフェルド残差に基づく検定が使われます
次のステップとして、以下の記事も参考にしてください。