Cox比例ハザードモデルは、生存時間解析でもっともよく使われる回帰モデルです。まず24人分の小さなデータで実行して結果を読むところまでを済ませ、そのあとで「ハザードとは何か」「モデルは何をしているのか」を解説します。数式は折りたたみの中に入れてあるので、必要になったときに開いてください。なお、コードの中の`options(digits = 3)`は出力の桁数をそろえるための設定で、解析結果そのものには影響しません。
## まずは動かしてみる
サブスクリプションサービスの契約者24人を24か月追跡し、解約までの月数を記録したデータとします。全員分を書き出します。
```r
library(survival)
options(digits = 3) # 出力の桁数をそろえる
churn <- data.frame(
id = 1:24,
plan = c(rep("monthly", 12), rep("yearly", 12)),
months = c( 2, 4, 7, 9, 13, 15, 19, 23, 24, 24, 24, 24,
8, 15, 20, 24, 24, 24, 24, 24, 24, 24, 24, 24),
churn = c( 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,
1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
usage = c( 6, 11, 5, 9, 13, 7, 12, 10, 8, 14, 10, 15,
7, 12, 9, 13, 9, 16, 11, 14, 8, 12, 15, 10)
)
churn$plan <- factor(churn$plan, levels = c("monthly", "yearly"))
churn
```
```
id plan months churn usage
1 monthly 2 1 6
2 monthly 4 1 11
3 monthly 7 1 5
4 monthly 9 1 9
5 monthly 13 1 13
6 monthly 15 1 7
7 monthly 19 1 12
8 monthly 23 1 10
9 monthly 24 0 8
10 monthly 24 0 14
11 monthly 24 0 10
12 monthly 24 0 15
13 yearly 8 1 7
14 yearly 15 1 12
15 yearly 20 1 9
16 yearly 24 0 13
17 yearly 24 0 9
18 yearly 24 0 16
19 yearly 24 0 11
20 yearly 24 0 14
21 yearly 24 0 8
22 yearly 24 0 12
23 yearly 24 0 15
24 yearly 24 0 10
```
変数は4つです。
| 変数 | 内容 |
| --- | --- |
| `plan` | 料金プラン。月払い（monthly）12人、年払い（yearly）12人 |
| `months` | 契約から解約までの月数。解約しなかった人は追跡終了の24 |
| `churn` | 1が解約、0が打ち切り（24か月時点で契約継続中） |
| `usage` | 契約初月の週あたり利用時間 |
解約は11人、残り13人は24か月の追跡が終わった時点でまだ契約中です。この13人は「解約しなかった」のではなく「少なくとも24か月は解約しなかった」ことしかわかりません。これを打ち切り（censoring）と呼びます。
Cox回帰は2ステップです。`Surv()`で「時間」と「イベントが起きたかどうか」を1つにまとめ、それを`coxph()`の左辺に置きます。
```r
fit1 <- coxph(Surv(months, churn) ~ plan, data = churn)
summary(fit1)
```
```
n= 24, number of events= 11
coef exp(coef) se(coef) z Pr(>|z|)
planyearly -1.336 0.263 0.680 -1.97 0.049 *
exp(coef) exp(-coef) lower .95 upper .95
planyearly 0.263 3.8 0.0694 0.997
Concordance= 0.67 (se = 0.068 )
Likelihood ratio test= 4.43 on 1 df, p=0.04
Wald test = 3.86 on 1 df, p=0.05
Score (logrank) test = 4.45 on 1 df, p=0.03
```
見るべきは`exp(coef)`の列で、これがハザード比（hazard ratio, HR）です。ここでは0.26なので、
**年払いの契約者は、月払いの契約者に比べて、解約が起こるスピードが0.26倍**
と読みます。逆向きに言えば、月払いは年払いの3.8倍（`exp(-coef)`の値）のスピードで解約していきます。95%信頼区間は0.07から1.00で、24人しかいないため非常に広くなっています。
ここまでがCox回帰の最短の流れです。以降で、ハザード比の読み方から順に見ていきます。
## ハザード比の読み方
ハザード比は「基準の群を1としたときの、イベントの起こりやすさの比」です。したがって1が「差がない」を意味します。
```
0.25 0.5 1.0 2.0 4.0
───────┼─────────┼─────────┼─────────┼─────────┼───────
←──── 起こりにくい ──┼── 起こりやすい ────→
基準（差なし）
▲
年払い 0.26
```
| ハザード比 | 意味 |
| --- | --- |
| 1.0 | 基準の群とまったく同じ。イベントの起こりやすさに差がありません |
| 2.0 | 基準の2倍のスピードでイベントが起こります |
| 0.5 | 基準の半分のスピード。50%低いと表現します |
| 0.26 | 74%低い（今回の年払いの結果） |
パーセントに直すときは、1との差をとります。1より小さければ $(1 - \text{HR}) \times 100$ が減少率、1より大きければ $(\text{HR} - 1) \times 100$ が増加率です。今回のHRは0.26なので $(1 - 0.26) \times 100 = 74$ となり、「年払いでは解約が74%少ない」と読めます。
**信頼区間が1をまたぐかどうか**が判断の目安になります。95%信頼区間に1が含まれていれば、「基準と差がある」とは言い切れません。今回の上限は0.997で、かろうじて1を下回っています。
### 基準（reference）はどう決まるか
Rでは、因子の最初の水準が基準になります。今回は`levels = c("monthly", "yearly")`と指定したので`monthly`が基準です。基準の群そのものは出力に現れません。「`planyearly`の行しかない」のはそのためです。
基準を入れ替えるには`relevel`を使います。
```r
churn$plan2 <- relevel(churn$plan, ref = "yearly")
fit1b <- coxph(Surv(months, churn) ~ plan2, data = churn)
print(summary(fit1b)$conf.int, digits = 4)
```
```
exp(coef) exp(-coef) lower .95 upper .95
plan2monthly 3.804 0.2629 1.003 14.42
```
年払いを基準にすると、月払いのHRは3.8になりました。単変量の0.26の逆数にあたります。同じ結果を裏返しに表示しているだけで、どちらを基準にしても解析の中身は変わりません。読み手に伝わりやすいほうを選びます。
### 連続変数の場合
連続変数には群がないので、「1単位低い人」が基準になります。`usage`のHR 0.74は「利用時間が1時間長い人は、1時間短い人に比べてハザードが0.74倍」という意味です。
## 出力の各項目
| 項目 | 意味 |
| --- | --- |
| `coef` | 回帰係数。ハザード比の対数です |
| `exp(coef)` | ハザード比。ふつう報告するのはこの値です |
| `se(coef)` | 係数の標準誤差 |
| `z` | 係数を標準誤差で割った値 |
| `Pr(>\|z\|)` | Wald検定のp値 |
| `exp(-coef)` | 逆向きのハザード比（1÷HR）。参照カテゴリを入れ替えた場合の値です |
| `lower .95` `upper .95` | ハザード比の95%信頼区間 |
| `Concordance` | 予測の順位一致度（c統計量）。0.5が偶然と同じ、1が完全 |
| 3つの検定 | モデル全体の検定。スコア検定はログランク検定と一致します |
行名が`planyearly`となっているのは、`plan`が因子で`monthly`が参照カテゴリだからです。「monthlyに比べてyearlyは」という読み方になります。
`Surv()`が何を作っているかも見ておきます。
```r
Surv(churn$months, churn$churn)
```
```
[1] 2 4 7 9 13 15 19 23 24+ 24+ 24+ 24+ 8 15 20 24+ 24+ 24+ 24+
[20] 24+ 24+ 24+ 24+ 24+
```
数値の後ろの`+`が打ち切りです。`24+`は「24か月より後に解約する」ことだけがわかっている状態を表し、24という数値そのものとしては扱われません。第2引数は1がイベント発生、0が打ち切りです（`TRUE`/`FALSE`でも構いません）。この0と1を逆に指定するのはよくある間違いで、係数の符号が反転します。
## 絵で確認する
回帰の前後で、群ごとの生存曲線も見ておきます。
```r
km <- survfit(Surv(months, churn) ~ plan, data = churn)
km
plot(km, col = c("steelblue", "tomato"), lwd = 2, ylim = c(0, 1),
xlab = "Months since contract", ylab = "Probability of staying",
main = "Kaplan-Meier curves by plan")
legend("bottomleft", legend = c("monthly", "yearly"),
col = c("steelblue", "tomato"), lwd = 2, bty = "n")
```
```
n events median 0.95LCL 0.95UCL
plan=monthly 12 8 17 9 NA
plan=yearly 12 3 NA NA NA
```
![[churn_km_plan.png|500]]
月払いは半数が解約するまで17か月、年払いは24か月時点でも半数に達しないため中央値が`NA`になっています。階段が数えられるほど粗いのは、イベントが11個しかないからです。
ログランク検定でも差を確認できます。
```r
survdiff(Surv(months, churn) ~ plan, data = churn)
```
```
N Observed Expected (O-E)^2/E (O-E)^2/V
plan=monthly 12 8 4.58 2.56 4.48
plan=yearly 12 3 6.42 1.83 4.48
Chisq= 4.5 on 1 degrees of freedom, p= 0.03
```
ログランク検定でわかるのは差の有無だけです。差の大きさを数値にし、他の変数で調整できることがCox回帰の利点になります。
## 変数を増やす
式の右辺に`+`でつなぐだけです。利用時間を加えます。
```r
fit2 <- coxph(Surv(months, churn) ~ plan + usage, data = churn)
summary(fit2)
```
```
n= 24, number of events= 11
coef exp(coef) se(coef) z Pr(>|z|)
planyearly -1.111 0.329 0.693 -1.6 0.109
usage -0.297 0.743 0.135 -2.2 0.027 *
exp(coef) exp(-coef) lower .95 upper .95
planyearly 0.329 3.04 0.0848 1.280
usage 0.743 1.35 0.5709 0.968
Concordance= 0.774 (se = 0.068 )
Likelihood ratio test= 10.2 on 2 df, p=0.006
```
ハザード比と信頼区間だけを取り出したいときは、`confint`を指数変換します。
```r
exp(cbind(HR = coef(fit2), confint(fit2)))
```
```
HR 2.5 % 97.5 %
planyearly 0.329 0.0848 1.280
usage 0.743 0.5709 0.968
```
## 係数を解釈する
**2値変数の場合**：`plan`のHRは0.33です。「利用時間が同じなら、年払いの解約スピードは月払いの0.33倍」と読みます。
**多変量にすると値が動きます**：`plan`のHRは単変量の0.26から0.33へ1に近づき、p値は0.049から0.109になりました。このデータでは年払いの人のほうが利用時間がやや長く（平均11.3時間 対 10.0時間）、単変量の推定にはその分が混ざっていたためです。
**連続変数は単位を変えられます**：`usage`のHR 0.74は1時間あたりの値ですが、1時間の差では実感しにくいところです。係数を5倍してから指数変換すれば、5時間あたりのハザード比になります。
```r
b <- coef(fit2)["usage"]
exp(5 * b) # 5時間あたりのHR
exp(5 * confint(fit2)["usage", ]) # その95%信頼区間
```
```
usage
0.227
2.5 % 97.5 %
0.0606 0.8481
```
このとき、モデル上は「5時間の差」がどの水準から測っても同じHRになる前提が置かれています。連続変数をそのまま入れるのは、対数ハザードが直線的に変化するという仮定であり、値の範囲が広い変数ではスプラインなどで確かめます。
**ハザード比はリスク比ではありません**：HRは各時点での「まだ解約していない人」の中での比です。24か月時点の解約割合の比とは別の量であり、追跡期間が長いほど、累積の差はHRから受ける印象より小さく見えることがあります。絶対的な差を知りたいときは、次の生存曲線を使います。
## 調整済み生存曲線を描く
Cox回帰の結果から、共変量を特定の値に固定した場合の生存曲線を描けます。`survfit`に`newdata`を渡します。
```r
nd <- data.frame(plan = factor(c("monthly", "yearly"), levels = levels(churn$plan)),
usage = c(10, 10))
sf <- survfit(fit2, newdata = nd)
summary(sf, times = c(6, 12, 24))
```
```
time n.risk n.event survival1 survival2
6 22 2 0.902 0.966
12 19 3 0.722 0.898
24 13 6 0.345 0.704
```
```r
plot(sf, col = c("steelblue", "tomato"), lwd = 2, ylim = c(0, 1),
xlab = "Months since contract", ylab = "Probability of staying",
main = "Adjusted survival curves (usage = 10)")
legend("bottomleft", legend = c("monthly", "yearly"),
col = c("steelblue", "tomato"), lwd = 2, bty = "n")
```
![[churn_adjusted_survival.png|500]]
利用時間が週10時間の人どうしで比べると、24か月継続している割合は月払いで34.5%、年払いで70.4%と推定されます。ハザード比0.33という相対的な指標が、絶対的な差としてはこの程度になる、という対応を確認できます。
`newdata`を省略すると共変量を平均値に固定した曲線が描かれますが、因子変数の「平均」には意味がないため、`newdata`を明示的に指定するほうが安全です。
## 解説：ハザードとは何か
ハザードは、日本語では瞬間死亡率、瞬間発生率などと訳されます。工学分野では故障率（failure rate）と呼びます。ひとことで言えば、**その瞬間にイベントが起こるスピード**です。
車のスピードメーターを思い浮かべると理解しやすくなります。針が指しているのは「今この瞬間の速さ」であって、走った距離ではありません。ハザードも同じで、時点ごとに値が変わります。
| 性質 | 内容 |
| --- | --- |
| 確率ではありません | 1を超えることがあります。1か月あたり2.0は「このペースが続けば平均0.5か月で1件起こる」という意味です |
| 単位を持ちます | 1か月あたり0.01は、1年あたりでは0.12になります |
| 瞬間の値です | 累積した量ではありません。走行距離にあたるのは累積ハザードという別の量です |
| 条件付きです | 分母は「その時点まで残っている人」です。時間が進むにつれて、比べている集団の中身が入れ替わっていきます |
生存曲線とは表裏の関係にあります。ハザードが速度、それを時間で積み上げた累積ハザードが走行距離、生存率が残っている割合です。ハザードが大きい時期ほど生存曲線は急に下がります。
> [!note]- 数式で確認する
> ハザード関数の定義は次のとおりです。
>
> $
> h(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}
> $
>
> 分子は「時点 $t$ までイベントが起きていない対象が、直後の短い時間 $\Delta t$ の間にイベントを起こす確率」です。これを $\Delta t$ で割っているため、全体としては単位時間あたりの量、つまり率（rate）になります。確率を時間で割っているので1を超えうるわけです。
>
> 累積ハザード $H(t)$ と生存関数 $S(t)$ は次の関係にあります。
>
> $
> H(t) = \int_0^t h(u)\,du, \qquad S(t) = \exp(-H(t))
> $
ハザードを実際に手で計算して確かめる手順は [[R - ハザード関数を手動計算で理解しよう]] にあります。
## 解説：Cox回帰は何をしているのか
Cox回帰の考え方は次の1行に集約されます。
**ハザード = ベースラインのハザード × 共変量による倍率**
ベースラインのハザードとは、基準となる人（すべての共変量が0の人）のハザードです。これは時間とともに自由に上下してよく、形を一切仮定しません。解約が契約直後に多かろうと1年後に多かろうと、モデルは気にしません。
一方、共変量による倍率のほうには時間が入りません。年払いの人のハザードは、契約直後だろうと20か月後だろうと、月払いの人の0.26倍で一定、という構造です。この「倍率が時間によらず一定」という部分が比例ハザードという名前の由来であり、`coxph`が推定しているのはこの倍率だけです。
倍率が時間によらないので、2人のハザードを割り算するとベースラインの部分が消えてしまいます。消えるからこそ、ベースラインの形を知らないままハザード比を推定できます。この性質のおかげで、Cox回帰は「分布を仮定しない」回帰として使えます（セミパラメトリック）。
推定にあたって使っているのは、イベントが起きた時点の情報だけです。「この時点で解約したのが、まだ残っていた人たちの中でこの人だったのはどのくらいもっともらしいか」を全イベントについて掛け合わせ、それを最大にする係数を探します。ここから次の性質が出ます。
- 実質的なサンプルサイズはイベント数で決まります。24人いても解約が11人なら情報は11個分です
- 時間の単位を変えても係数は変わりません。使っているのはイベントが起きた順序だけだからです
- 打ち切りの人も「その時点まで残っていた人」として分母に貢献し、情報を捨てずに済みます
> [!note]- 数式で確認する
> Cox比例ハザードモデルは、ハザード関数を次の形で表します。
>
> $
> h(t \mid x) = h_0(t) \exp(\beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p)
> $
>
> $h_0(t)$ がベースラインハザード、$\exp(\beta_1 x_1 + \cdots)$ が共変量による倍率です。両辺の対数をとると線形モデルの形になります。
>
> $
> \log h(t \mid x) = \log h_0(t) + \beta_1 x_1 + \cdots + \beta_p x_p
> $
>
> **係数がハザード比になる理由**：$x_1$ だけが1違う2人のハザードを比べると、
>
> $
> \frac{h(t \mid x_1 = a+1)}{h(t \mid x_1 = a)}
> = \frac{h_0(t)\exp(\beta_1 (a+1) + \cdots)}{h_0(t)\exp(\beta_1 a + \cdots)}
> = \exp(\beta_1)
> $
>
> となり、$h_0(t)$ が約分されて $t$ が残りません。したがって $\text{HR} = \exp(\beta_1)$ であり、この値は時点によらず一定になります。
>
> **部分尤度**：Coxは、イベントが起きた各時点について「リスク集合の中で実際にイベントを起こしたのがその人である確率」を掛け合わせた部分尤度を最大化します。
>
> $
> L(\beta) = \prod_{i:\, \delta_i = 1} \frac{\exp(\beta^{\top} x_i)}{\sum_{j \in R(t_i)} \exp(\beta^{\top} x_j)}
> $
>
> $R(t_i)$ は時点 $t_i$ の直前まで観察されていた人の集合（リスク集合）、$\delta_i$ はイベント発生を表す指標です。分子も分母も $h_0(t_i)$ を含むため約分され、ベースラインハザードを推定せずに $\beta$ が求まります。
## まとめ
Cox回帰の手順は、`Surv()`でアウトカムを作り、`coxph()`で推定し、`summary()`の`exp(coef)`とその信頼区間を読む、という3段階です。得られるハザード比は「イベントが起こるスピードの比」で、1が「差がない」を意味します。絶対的な差を示したいときは`survfit()`で調整済み生存曲線を描きます。
結果を報告する前の確認（比例ハザードの仮定、交互作用、イベント数の目安）は [[R - Cox比例ハザードモデルの発展]] で扱います。
## 関連項目
- [[R - Cox比例ハザードモデルの発展]]
- [[R - カプランマイヤー法をマスターしよう]]
- [[R - ハザード関数を手動計算で理解しよう]]
- [[R - 生存時間解析の基礎]]
- [[R - 回帰係数を理解しよう]]
- [[Stata - 生存時間解析]]