多変量Coxモデルを組んで比例ハザード性を確かめたところ、治療変数は問題ないのに、調整のために入れた全身状態の指標が仮定を破っているケースがあります。関心があるのは治療効果のほうで、全身状態のハザード比を報告する予定はありません。
この場合、その変数を層別因子に回すという手があります。層ごとに別々の基準ハザードを許し、係数は層をまたいで共通と置く。仮定を破っていた変数がモデルの係数から消えるので、比例ハザード性の問題も一緒に消えます。
代わりに失うものがあります。その変数のハザード比は推定されなくなり、層を細かくしすぎると治療効果の精度も落ちます。何を捨てて何を守るのかを、実データで確認します。
層ごとに基準ハザードを分ける
層別Coxモデルは、層$s$ごとに違う基準ハザードを持たせます。
$$h(t \mid \mathbf{x}, s) = h_{0s}(t)\exp(\boldsymbol{\beta}^\top \mathbf{x}) \tag{1}$$
式(1)で層によって変わるのは$h_{0s}(t)$だけで、係数$\boldsymbol{\beta}$は共通です。層ごとに時間的な形が違ってよい代わりに、共変量の効果は層をまたいで同じと仮定します。
推定も層の中で閉じます。部分尤度のリスク集合を、同じ層の人だけで作ります。
$$L(\boldsymbol{\beta}) = \prod_{s} \prod_{j \in s} \frac{\exp(\boldsymbol{\beta}^\top \mathbf{x}_{(j)})}{\sum_{i \in R_s(t_j)} \exp(\boldsymbol{\beta}^\top \mathbf{x}_i)} \tag{2}$$
式(2)の分母$R_s(t_j)$は、時刻$t_j$に層$s$で追跡中だった人の集合です。層をまたいだ比較は一切行われません。そのため、層の効果そのものは推定されず、$h_{0s}(t)$は係数として現れないまま約分されます。
veteranデータで確かめます。Karnofsky指標（70以上か未満か）は比例ハザード性を明確に破っており、関心のある治療変数trtは破っていません。
何が消えて何が残るか
library(survival)
veteran$kg <- factor(ifelse(veteran$karno >= 70, "high", "low"))
adjusted <- coxph(Surv(time, status) ~ trt + kg, data = veteran)
stratified <- coxph(Surv(time, status) ~ trt + strata(kg), data = veteran)
cox.zph(adjusted)
cox.zph(stratified)#> [共変量として投入] chisq df p
#> trt 1.68 1 0.19528
#> kg 12.46 1 0.00042
#> GLOBAL 15.94 2 0.00035
#>
#> [層別] chisq df p
#> trt 3.19 1 0.074
#> GLOBAL 3.19 1 0.074層別にした途端、比例ハザード性の検定に引っかかる項目がなくなります。kgは係数として推定されていないので、検定の対象にもなりません。仮定違反が解決したのではなく、仮定を課す対象から外しただけです。
ここが選択の分かれ目になります。層別因子にできるのは、調整のためだけに入れている変数です。その変数自体の効果量を報告したいなら、共変量として入れたうえで、比例ハザード性の破れを別の方法で扱うことになります。時変係数を入れる、期間を分けて推定する、あるいはRMSTのように比に依存しない指標に切り替える。
実務では、無作為化のときに層別した因子をそのまま解析でも層別する、という使い方が多いと思います。施設、病期、年齢層。これらは割付のバランスを取るために使った変数で、効果量の報告対象ではありません。層別解析にすれば、割付の設計と解析の構造が揃います。
もうひとつの用途は、比例ハザード性を破る交絡因子です。今回のKarnofsky指標がこれにあたります。ただし、破っているという事実自体が臨床的に重要な情報である場合もあります。「全身状態の影響は診断直後に強く、時間が経つと薄れる」というのは、それ自体が知見です。層別するとこの情報は本文に書けなくなるので、記述として残したいなら時変係数のほうを選びます。
層を細かくすると精度が落ちる
式(2)のリスク集合は層の中だけで作られます。層が小さいと、比較できる相手がいなくなります。極端な話、1つの層に1人しかいなければ、その人はどのイベント時点でも分母を独占するので、部分尤度に何も情報を与えません。
どのくらいで問題になるのかを確かめます。400例、治療効果を真のハザード比0.7に固定し、層の数だけを2から200まで変えました。層ごとの予後差は乱数で与えています。以下は説明のための架空データです。
推定値の偏りはほとんど生じません。落ちるのは精度だけです。それでも、検出力の観点では無視できない損失になります。
多施設研究で施設ごとに層別したくなる場面がありますが、1施設あたり数例という研究では勧めません。この場合の選択肢は2つあります。施設をいくつかのグループにまとめて層の数を減らすか、施設を層別せずにロバスト分散で施設内の相関を扱うか。施設の予後差そのものを推定したいなら、frailtyモデルという選択肢もあります。
目安として、私は1層あたりのイベント数が10を下回るなら層の統合を勧めます。図3の横軸は人数ですが、実際に効くのはイベント数です。イベントが1件も起きていない層は、部分尤度に一切寄与しません。
層をまたいで効果が同じという仮定
式(1)は係数$\boldsymbol{\beta}$を層で共通と置いています。これは検証すべき仮定です。層と共変量の交互作用項を入れて、係数が層によって違わないかを確かめられます。
interaction <- coxph(Surv(time, status) ~ trt * kg + strata(kg), data = veteran)
anova(stratified, interaction)交互作用が有意なら、治療効果が層によって違うことになります。全体で1つのハザード比を報告する意味が薄れるので、層ごとの推定値を並べる形に変えることになります。ただし、この検定は検出力が低いので、有意でなかったことを効果修飾がない証拠として扱わないほうがよいでしょう。Cox比例ハザードモデルで見たとおり、有意でないことは仮定の成立を意味しません。
層別ログランク検定との対応
層別は検定の側にも同じ形で入ります。層ごとに観測数と期待数の差を計算してから、層をまたいで足し上げる。ログランク検定で見た式の期待値を、層の中で作り直したものです。
survdiff(Surv(time, status) ~ trt + strata(kg), data = veteran)#> N Observed Expected (O-E)^2/E (O-E)^2/V
#> trt=1 69 64 65.5 0.0332 0.0751
#> trt=2 68 64 62.5 0.0347 0.0751
#>
#> Chisq= 0.1 on 1 degrees of freedom, p= 0.8層別Coxモデルのスコア検定は0.08でp=0.778、層別ログランク検定は0.1でp=0.8。両者が一致するのは部分尤度のときと同じ関係で、層別しても崩れません。
層別しないログランク検定ではカイ二乗が0.008、p=0.928でした。この例では治療効果がもともとないので結論は変わりませんが、層ごとの予後差が大きい研究では、層別したほうが検出力が上がることがあります。層内で比べることで、層間のばらつきが誤差から取り除かれるためです。
割付時に層別したなら検定でも層別する、という原則はここから来ています。Kahan and Morris（Stat Med 2012）は、層別割付を使いながら解析で層を無視した論文が多いことと、その場合に検定が保守的になることを報告しています。
連続変数で層別するとき
層別因子はカテゴリでなければなりません。年齢や検査値のような連続変数を層別因子にするには、カテゴリに切る必要があります。
切り方で結果が変わります。層を細かくすれば各層内の均質性は上がりますが、図3のとおり精度が落ちます。粗くすれば精度は保たれますが、層内に予後の違う患者が混ざり、層別した意味が薄れます。
目安として、4分位や5分位で切る形がよく使われます。1層あたりのイベント数が十分に残る範囲で、できるだけ細かく、という判断です。等間隔で切ると層ごとの人数が偏るので、分位点で切るほうが安定します。
そもそも連続変数を層別するのは、その変数が比例ハザード性を破っていて、かつ効果量を報告しなくてよい場合に限られます。仮定が保たれているなら、共変量としてスプラインなどで柔軟に入れるほうが情報を失いません。
層別と交互作用は別のこと
層別Coxモデルは、層ごとに基準ハザードを変えます。共変量の係数は共通です。一方、交互作用項を入れたモデルは、基準ハザードを共通にしたまま係数を層ごとに変えます。
この2つは直交する操作です。「施設によって全体の予後水準が違う」なら層別、「施設によって治療の効き方が違う」なら交互作用。両方あるなら両方入れます。
混同されやすいのは、どちらも「施設を考えに入れた」と表現されるためだと思います。論文で「施設で調整した」とだけ書かれていると、どちらの構造を仮定したのか読み取れません。層別したのか、共変量として入れたのか、交互作用を許したのかを書き分ける必要があります。
報告の仕方
層別した変数を明記します。「Karnofsky指標で層別したCox回帰」と書けば、その変数のハザード比が表にない理由が読者に伝わります。書かないと、単に入れ忘れたように見えます。
層の数と、1層あたりの症例数・イベント数も添えます。層が細かい解析では、これがないと精度の妥当性を判断できません。
そして、なぜ層別したのかを1文書きます。事前規定の層別因子だったのか、比例ハザード性の破れへの対処だったのか。後者なら、破れをどう確認したかも本文か補足資料に置きます。
比例ハザード性の確認そのものは比例ハザード仮定の検証、破れているときの効果指標の選択はRMST、群間比較の重み付けはログランク検定で扱っています。
層別因子をどこまで解析に持ち込むか、施設をどう扱うか。このあたりは試験の設計と施設あたりの症例数を見ないと決まりません。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
Kalbfleisch JD, Prentice RL. The Statistical Analysis of Failure Time Data. 2nd ed. Wiley; 2002.
Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000.
Kahan BC, Morris TP. Improper analysis of trials randomised using stratified blocks or minimisation. Statistics in Medicine. 2012;31(4):328-340.
Glidden DV, Vittinghoff E. Modelling clustered survival data from multicentre clinical trials. Statistics in Medicine. 2004;23(3):369-388.
Klein JP, Moeschberger ML. Survival Analysis: Techniques for Censored and Truncated Data. 2nd ed. Springer; 2003.
