[[R - Cox比例ハザードモデルの基礎]] の続きです。Cox回帰を実行して結果を読めるようになったあと、結果を報告する前に確認しておきたいことを扱います。 ## 使うデータ 基礎編と同じ、サブスクリプションサービスの契約者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")) fit2 <- coxph(Surv(months, churn) ~ plan + usage, data = churn) ``` `plan`は料金プラン（月払い/年払い）、`months`は解約までの月数、`churn`は1が解約で0が打ち切り、`usage`は契約初月の週あたり利用時間です。24人のうち11人が解約しています。 ## 比例ハザードの仮定を確かめる Cox回帰は「ハザード比が追跡期間を通じて一定」という仮定に立っています。この仮定が崩れる典型例は、手術と薬物治療の比較のように、初期はある群が不利で後半は有利になる場合です。 `cox.zph`はこの仮定を検定します。p値が小さいと、ハザード比が時間とともに変化している疑いが強くなります。 ```r cox.zph(fit2) ``` ``` chisq df p plan 0.0848 1 0.77 usage 0.4845 1 0.49 GLOBAL 0.5180 2 0.77 ``` いずれもp値が大きく、比例ハザードを疑う根拠はありません。ただしイベントが11個では検出力がほとんどないため、この結果は「確かめた」というより「確かめようがなかった」に近いものです。 仮定が破れている例も見ておきます。`survival`パッケージに含まれる`veteran`（肺がん治療の臨床試験データ、137人）で、Karnofskyスコア（全身状態の指標）を共変量にした場合です。 ```r fit_v <- coxph(Surv(time, status) ~ trt + karno, data = veteran) zph_v <- cox.zph(fit_v) zph_v plot(zph_v[2], main = "Scaled Schoenfeld residuals: karno (veteran)") abline(h = coef(fit_v)["karno"], col = "red", lty = 2) ``` ``` chisq df p trt 0.307 1 0.57942 karno 11.959 1 0.00054 GLOBAL 13.577 2 0.00113 ``` ![[cox_zph_veteran.png|500]] 図の実線は各時点での係数の推定値、赤い破線はCox回帰が推定した一定の係数です。実線が右上がりで、追跡初期は係数が赤線より小さく（スコアの効果が強く）、後期には0に近づいています。この場合に報告される単一のHRは、追跡期間全体の平均のようなものであり、時期による違いを覆い隠します。 ### log-logプロットで見る もうひとつの定番の確認方法が log-log プロットです。群ごとの生存曲線から $\log(-\log S(t))$ を計算し、時間の対数に対して描きます。**比例ハザードが成り立っていれば、2本の線は上下に平行にずれるだけ**になります。交差したり、間隔が広がったり狭まったりしていれば、ハザード比が時間とともに変わっていることを意味します。Stataの`stphplot`に相当する図です。 Rでは`plot`に`fun = "cloglog"`を指定します。比較のため、同じ`veteran`データで、仮定に問題のない変数（組織型）と問題のある変数（Karnofskyスコア）を並べます。 ```r v <- veteran v$cell_g <- factor(ifelse(v$celltype == "squamous", "squamous", "other"), levels = c("other", "squamous")) v$karno_g <- factor(ifelse(v$karno >= 60, "high", "low"), levels = c("low", "high")) par(mfrow = c(1, 2)) plot(survfit(Surv(time, status) ~ cell_g, data = v), fun = "cloglog", col = c("steelblue", "tomato"), lwd = 2, xlab = "Days (log scale)", ylab = "log(-log(S(t)))", main = "Parallel: cell type") plot(survfit(Surv(time, status) ~ karno_g, data = v), fun = "cloglog", col = c("steelblue", "tomato"), lwd = 2, xlab = "Days (log scale)", ylab = "log(-log(S(t)))", main = "Not parallel: Karnofsky score") ``` ![[cox_loglog_veteran.png|600]] 左の組織型は、序盤を除けば2本の間隔がおおむね一定に保たれています。`cox.zph`のp値も0.23で、仮定に問題は見当たりません。右のKarnofskyスコアは、序盤に大きく開いていた差が時間とともに縮まり、終盤ではほとんど重なっています。`cox.zph`のp値は0.00005で、先ほどのSchoenfeld残差の図と同じ現象を別の角度から見ていることになります。 読むときの注意点があります。 - 序盤はイベントが少なく推定が不安定なので、左端の乱れは気にしません - カテゴリ変数にしか使えません。連続変数は群に分ける必要があり、分け方によって見え方が変わります - 群ごとのイベント数が少ないと階段が数段しか描かれず、平行かどうか判断できません。今回の24人のデータがこれにあたります > [!note]- なぜ平行になるのか > 比例ハザードモデルのもとで、生存関数は次のように書けます。 > > $ > S(t \mid x) = \exp\left(-H_0(t)\, e^{\beta x}\right) > $ > > 両辺の対数をとって符号を変え、もう一度対数をとると、 > > $ > \log(-\log S(t \mid x)) = \log H_0(t) + \beta x > $ > > となります。右辺の第1項は共変量によらず、第2項は時間によりません。したがって群による違いは $\beta x$ の分だけの縦方向の平行移動になり、線は平行に並びます。横軸を対数にとるのは初期の挙動を見やすくするための慣例で、平行性の判断そのものには影響しません。 仮定が満たされないときの対処には次のような選択肢があります。 | 対処 | 書き方 | 使いどころ | | --- | --- | --- | | 層別 | `coxph(Surv(t, d) ~ x + strata(z))` | `z`の効果自体には関心がなく、調整だけしたい場合 | | 時間依存係数 | `coxph(Surv(t, d) ~ x + tt(x), tt = function(x, t, ...) x * log(t))` | 係数の時間変化そのものを表現したい場合 | | 期間を区切る | `survSplit`で期間分割し、期間ごとにHRを推定 | 前半と後半でHRを分けて報告したい場合 | | 別の指標を使う | RMST（制限付き平均生存時間）など | 比例ハザードに依存しない要約が欲しい場合 | なお、`cox.zph`の検定はサンプルサイズが大きいと実質的に問題のない逸脱でも有意になり、小さいと検出力が足りません。p値だけで判断せず、図を見て逸脱の向きと大きさを確認します。 ## 交互作用を入れる 「年払いの効果は、よく使う人と使わない人とで違うのではないか」という問いは交互作用（interaction、効果修飾）の問題です。2つの変数の掛け算にあたる項を`plan:usage`と書いて、式に足します。 ```r fit_int <- coxph(Surv(months, churn) ~ plan + usage + plan:usage, data = churn) summary(fit_int)$coefficients ``` ``` coef exp(coef) se(coef) z Pr(>|z|) planyearly 0.669 1.953 3.331 0.201 0.8407 usage -0.258 0.772 0.148 -1.743 0.0814 planyearly:usage -0.187 0.829 0.350 -0.535 0.5929 ``` 3行目が交互作用項です。読み方は次のようになります。 | 行 | 意味 | | --- | --- | | `planyearly` | 利用時間が0の人での、年払いのハザード比。単独では意味を持ちません | | `usage` | 月払いの人での、利用時間1時間あたりのハザード比 | | `planyearly:usage` | ハザード比の比。利用時間が1時間増えるごとに、年払いのハザード比が0.83倍ずつ変わります | なお、`plan + usage + plan:usage`は`plan * usage`と略記できます。結果は同じで、他の資料ではこちらの書き方をよく見かけます。 交互作用の有無を判定するときは、交互作用を入れたモデルと入れないモデルを尤度比検定で比べます。 ```r anova(fit2, fit_int) ``` ``` loglik Chisq Df Pr(>|Chi|) 1 -27.1 2 -27.0 0.31 1 0.58 ``` 交互作用項の行のp値は0.59で、こちらも似た値です。今回のように交互作用項が1行しかない場合は、どちらを見ても結論は変わりません。ただし次の場合には、尤度比検定でまとめて判定する必要があります。 **交互作用項が複数行に分かれるとき**：3水準以上のカテゴリ変数との交互作用では、交互作用項が複数行になります。3か月ごとの支払いプラン（quarterly）の契約者12人を加えて、プランを3水準にしてみます。 ```r quarterly <- data.frame( id = 25:36, plan = "quarterly", months = c( 5, 10, 14, 21, 24, 24, 24, 24, 24, 24, 24, 24), churn = c( 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0), usage = c( 8, 6, 11, 13, 10, 12, 9, 15, 7, 14, 11, 9) ) churn3 <- rbind(churn[, c("id", "plan", "months", "churn", "usage")], quarterly) churn3$plan <- factor(churn3$plan, levels = c("monthly", "quarterly", "yearly")) f0 <- coxph(Surv(months, churn) ~ plan + usage, data = churn3) f1 <- coxph(Surv(months, churn) ~ plan + usage + plan:usage, data = churn3) summary(f1)$coefficients ``` ``` coef exp(coef) se(coef) z Pr(>|z|) planquarterly -1.2782 0.279 2.460 -0.520 0.6034 planyearly 0.7925 2.209 3.323 0.238 0.8115 usage -0.2517 0.777 0.145 -1.734 0.0829 planquarterly:usage 0.0322 1.033 0.253 0.127 0.8990 planyearly:usage -0.1987 0.820 0.350 -0.569 0.5697 ``` 交互作用項が2行あるので、どちらか一方のp値だけでは「プランによって利用時間の効果が違うか」を判定できません。しかもこの2つは月払いを基準にした値なので、基準を3か月プランに変えると数字が変わります。 ```r churn3$plan_q <- relevel(churn3$plan, ref = "quarterly") summary(coxph(Surv(months, churn) ~ plan_q + usage + plan_q:usage, data = churn3))$coefficients ``` ``` coef exp(coef) se(coef) z Pr(>|z|) plan_qmonthly:usage -0.0322 0.968 0.253 -0.127 0.899 plan_qyearly:usage -0.2309 0.794 0.382 -0.605 0.545 ``` （プランの主効果の行は省略しています）。基準の取り方で値が動くものを個別に見ても、「どこかに交互作用があるか」の答えにはなりません。2行をまとめて検定します。 ```r anova(f0, f1) ``` ``` loglik Chisq Df Pr(>|Chi|) 1 -44.4 2 -44.2 0.44 2 0.8 ``` 自由度が2になっているのが、2行分をまとめて検定したことを表しています。 **推定が不安定なとき**：p値の計算に使われるWald検定は、係数の推定が不安定な状況では尤度比検定より信頼できなくなります。 実際、`planyearly`の標準誤差が3.33に膨らんでいることからわかるように、24人・11イベントでは交互作用の推定はまったく安定していません。交互作用の検出には主効果の場合よりはるかに多くのイベントが必要で、このp値は「交互作用がない」ではなく「何も言えない」と読むべきものです。交互作用を調べる予定があるなら、その検出力を見込んで研究を設計しておく必要があります。 なお、比例ハザードが破れている状態は「共変量と時間の交互作用がある」と言い換えられます。前の節の表に挙げた`tt()`は、まさに時間との交互作用項をモデルに入れる書き方です。 ## 実務上の注意点 **イベント数が実質的なサンプルサイズです**：24人いても解約が11人なら、投入できる変数はごく少数です。目安として共変量1つあたり10イベント程度は必要とされるので、この例の多変量モデル（2変数で11イベント）は本来なら変数が多すぎます。ここでは書き方を示すために2変数を入れていますが、実際の解析では、投入する変数の数をイベント数から逆算して決めます。 **打ち切りが無情報である必要があります**：「解約しそうな人ほど追跡から脱落する」といった状況では推定が偏ります。追跡期間の終了による打ち切り（今回の24か月）は通常この条件を満たします。 **時点0の定義を先に決めます**：曝露の判定より後に開始時点を置くと、その間はイベントが起こりえない期間が生まれ、曝露群が有利に見えます（immortal time bias）。開始時点、曝露の定義、追跡終了の定義を先に固めます。 **同時点データ（tie）の扱い**：`coxph`の既定は`ties = "efron"`で、通常はこのままで問題ありません。今回のように月単位で丸めた時間は同着が生じやすく、同着が非常に多い場合は`ties = "exact"`も検討します。 ## まとめ Cox回帰の結果を報告する前に確認することを3つ扱いました。 | 確認すること | 方法 | | --- | --- | | ハザード比が期間を通じて一定か | `cox.zph()`の検定と、Schoenfeld残差または log-log プロット | | 効果が群によって違わないか | 交互作用項を入れ、`anova()`で尤度比検定 | | そもそも推定に足る情報があるか | イベント数を数え、共変量の数を決める | いずれも、p値ひとつで機械的に判定できるものではありません。イベント数が少なければ検定の力は弱く、多ければ小さな逸脱でも有意になります。図を見て、逸脱の向きと大きさを自分で判断します。 ## 関連項目 - [[R - Cox比例ハザードモデルの基礎]] - [[R - カプランマイヤー法をマスターしよう]] - [[R - ハザード関数を手動計算で理解しよう]] - [[R - 交互作用項を理解しよう]] - [[生物統計における仮定とは]] - [[Stata - 生存時間解析]]