長期追跡の生存曲線を見ていると、途中からほとんど下がらなくなることがあります。がんの長期生存データでは古くから扱われてきた形で、Boag（1949）やFarewell（1982）の治癒モデルも、このような長期生存者を統計モデルの中でどう扱うかという問題から発展してきました。
ここで通常のCox回帰を当てはめること自体が間違いというわけではありません。観測期間内の群間差をハザード比として評価することはできます。ただし、Cox回帰には「イベントを起こさない集団」という潜在的な群は入っていないので、曲線の平坦部を治癒割合として分けて推定することはできません。
ワイブル分布など通常のパラメトリックモデルでは、さらに長期の外挿が問題になります。一般的なワイブル生存関数は時間とともに0へ近づくため、明確な平坦部が残るデータにそのまま当てはめると、長期生存の形が合わないことがあります。
そこで使われるのが治癒モデルです。混合治癒モデルでは、イベントを長期的に経験しない群と、イベントを経験し得る群を分け、前者の割合と後者のイベント時間を別々に推定します。
治癒割合と非治癒者を分ける
混合治癒モデルの生存関数は、次のように書けます。
$$S(t) = \pi + (1-\pi)\,S_u(t) \tag{1}$$
$\pi$は治癒割合、$S_u(t)$は非治癒者の生存関数です。非治癒者の生存確率が時間とともに0へ近づけば、集団全体の生存関数$S(t)$は$\pi$に近づきます。生存曲線の平坦部をモデルの中に直接持たせているわけです。
共変量を入れる場所も2か所あります。治癒割合にはロジスティック回帰を使い、非治癒者のイベント時間にはAFTモデルやCox型のモデルを組み合わせます。
たとえば非治癒者側をAFTモデルで書くと、次の形になります。
$$\pi(\mathbf{x}) = \frac{\exp(\boldsymbol{\gamma}^\top\mathbf{x})}{1+\exp(\boldsymbol{\gamma}^\top\mathbf{x})}, \qquad S_u(t \mid \mathbf{x}) = S_0\!\left(\frac{t}{\exp(\boldsymbol{\beta}^\top\mathbf{x})}\right) \tag{2}$$
$\boldsymbol{\gamma}$は治癒割合との関連を表す係数、$\boldsymbol{\beta}$は非治癒者のイベント時間との関連を表す係数です。式(2)では、1つの共変量について、治癒割合側の関連と非治癒者側の時間の関連を分けて推定できます。
ここは解釈を少し慎重にした方がよいところです。治療群を共変量に入れたからといって、観察研究で「治療が治癒者を増やした」とそのまま因果的に読めるわけではありません。モデルから分かるのは、治癒割合側と非治癒者側のどちらに群差が現れているかです。
通常のCox回帰では、この2つを別々のパラメータとしては表示しません。治癒モデルを使う理由は、ハザード比より高度なモデルだからではなく、研究で知りたい対象がこの2つに分かれている場合にあります。
追跡30か月では平坦部を読めない
Kaplan-Meier曲線の右端が水平になったのを見て、「治癒者がいる」と判断したくなることがありますが、ここは注意が必要です。
追跡後半ではリスク集合が小さくなります。イベントが起きていないから曲線が平らなのか、それとも観察できる患者がほとんど残っていないから動かないのかは、曲線だけでは区別できません。Kaplan-Meier推定量で扱った右端の不安定さが、そのまま治癒割合の推定にも効いてきます。
どの程度影響するのかを、治癒割合を0.25と0.45に設定してランダムに生成した600例のデータセットで確認しました。データの生成条件は同じにして、観察できる追跡期間の上限だけを30か月、60か月、120か月、240か月と変えています。
#> 追跡上限 イベント数 治癒割合の推定（真値 0.25 / 0.45）
#> 30か月 280 0.370 / 0.446
#> 60か月 368 0.292 / 0.362
#> 120か月 397 0.265 / 0.403
#> 240か月 400 0.260 / 0.407数字の動きが大きいのは対照群です。30か月では真値0.25に対して0.370と推定されましたが、追跡を120か月まで延ばすと0.265まで下がりました。
このシミュレーションでは、追跡が短いと治癒割合が高めに推定されました。ただし、「短い追跡なら必ず過大評価になる」という一般則として読むのは正確ではなく、どの程度、どちらの方向へずれるかは、非治癒者のイベント時間分布や打ち切りの仕方にも左右されます。
私なら、治癒モデルを当てはめる前に平坦部のリスクテーブルを確認します。曲線が水平になったあとにも相当数が追跡され、その状態が十分な期間続いているかを見るためです。右端に数例しか残っていなければ、曲線の形だけを根拠に治癒割合を推定するのは適切ではありません。
Cox回帰では有意でも内訳は不確か
同じ600例について通常のCox回帰を当てはめると、ハザード比は0.765、95%信頼区間は0.623〜0.939、ログランク検定はp=0.010でした。
治癒モデルで群差を2つに分けると、治癒割合側のオッズ比は1.38（95%信頼区間0.88〜2.15）、非治癒者側の時間比は1.21（0.89〜1.66）となり、どちらの95%信頼区間も1を含みました。
ここだけを見ると、Cox回帰では差が出たのに、治癒モデルでは差がなくなったように見えます。
ただし、検定している帰無仮説が違います。Cox回帰とログランク検定は、生存時間全体の群差を1つのパラメータまたは1自由度の検定として扱います。治癒モデルでは、その群差を治癒割合側と非治癒者側に分けて推定するため、それぞれの推定値には別々の不確実性が付きます。
今回のデータで、治癒割合側と非治癒者側の2係数が同時に0という帰無仮説を尤度比検定で調べると、統計量は6.81、自由度2でp=0.033でした。つまり、このデータではモデル全体として群差を検出している一方、どちらの成分にどの程度の差があるのかは、個別の信頼区間ではかなり不確かです。
「Cox回帰は2つの効果をまとめるので、治癒モデルより常に検出力が高い」と一般化するのは避けた方がよいでしょう。治癒割合だけが変わる場合、非治癒者の分布だけが変わる場合、2つが逆方向に動く場合では、検定の挙動も違います。
研究目的が単純な群間の生存時間差なら、まずCox回帰やログランク検定で答えられることがあります。治癒モデルを主要な解析に置く理由があるのは、「長期的にイベントを経験しない割合」と「非治癒者のイベント時間」を分けて推定すること自体が研究上の問いになっている場合です。
パラメータを2つに分けて個別に推定したいのであれば、その不確実性まで見込んだデータ量が必要です。通常のハザード比だけを前提にした症例数の考え方を、そのまま治癒割合の比較に持ち込むことはできません。
分布だけで治癒割合が13.5ポイント動いた
60か月まで追跡した同じデータに、非治癒者の分布だけを変えて混合治癒モデルを当てはめました。
ワイブル分布では対照群の治癒割合が0.292、対数ロジスティック分布では0.157です。その差は13.5ポイントありました。
なぜここまで動くのでしょうか。
追跡終了時点でイベントを経験していない患者について、モデルは「治癒者なのか、非治癒者だがイベントがまだ起きていないのか」をデータから分けなければなりません。非治癒者側に裾の長い分布を置けば、遅い時点までイベントを起こさない患者を非治癒者側で説明しやすくなり、その分だけ治癒割合の推定値が下がることがあります。
今回のAICはワイブル3421.9、対数ロジスティック3422.8で、差は0.9でした。少なくともこの結果では、AICだけから一方を明確に選び、そのモデルの治癒割合だけを確定値のように報告するのは難しい状態です。
対数正規分布では、さらに極端な結果になりました。治癒割合の推定値は0.018まで下がり、オッズ比の信頼区間上限は2,400万を超えています。
この0.018を「治癒割合は約2%」と解釈するのは適切ではありません。信頼区間がここまで広いときは、データから治癒割合を十分に識別できていない可能性を先に疑います。推定値そのものより、尤度の形や収束、標準誤差を確認すべき場面です。
治癒割合を主要な結果として扱うなら、私は1つの分布だけで結論を出す方法は適切ではありません。ワイブル、対数正規、対数ロジスティックなど、臨床的・統計的に候補となる分布を変えた感度分析を併記した方が、結果がどの程度安定しているかを読者が判断できます。
パラメトリック生存モデルでは分布選択そのものを扱いましたが、治癒モデルではその選択が非治癒者の曲線だけでなく、治癒割合の推定にも直接影響します。
がん登録では相対生存も選択肢になる
がん登録データで「治癒」を扱う場合には、混合治癒モデルとは別に相対生存を使う方法があります。Lambertらは2007年のBiostatisticsで、人口ベースのがん生存データから治癒割合を推定する方法を整理しています。
相対生存では、対象患者の死亡と、同じ年齢・性別・暦年などに対応する一般人口の期待死亡を比較します。そして、がんに関連すると考える超過死亡が長期的に認められなくなる状態を、統計的な治癒として扱います。
この定義は、混合治癒モデルで単純に「イベントを経験しない潜在群」を置く方法とは違います。
特に全死亡を扱うがん登録では、患者はがん以外の原因でも死亡します。「死亡しない人が一定割合いる」という通常の混合モデルを、そのまま全死亡に当てはめるのは解釈が難しくなります。一般人口にも存在する死亡を分けて考えたいのであれば、相対生存率の枠組みの方が研究目的に合うことがあります。
一方、再発をイベントとする研究では話が違います。再発前の死亡は、その後の再発を観察できなくするため、競合リスクとして扱うべきかを考える必要があります。生存曲線が平坦になったという形だけから、治癒と競合リスクのどちらで説明すべきかは決まりません。
つまり、最初に決めるべきなのはモデル名ではなく、「何を治癒と定義するのか」です。再発しないことなのか、疾患による超過死亡がなくなることなのか、全死亡を含む生存なのかで、選ぶ解析も変わります。
Rでは収束と感度分析まで確認する
混合治癒モデルでは、打ち切り例の尤度に治癒者と非治癒者の両方が入ります。式(1)の構造をそのまま書くと、次の形です。
negll <- function(par, time, event, group) {
cure <- plogis(par[1] + par[2] * group)
shape <- exp(par[3])
scale <- exp(par[4] + par[5] * group)
surv_u <- exp(-(time / scale)^shape)
dens_u <- (shape / scale) * (time / scale)^(shape - 1) * surv_u
-sum(ifelse(event == 1, log((1 - cure) * dens_u),
log(cure + (1 - cure) * surv_u)))
}イベントが起きた患者は非治癒者なので、尤度への寄与は$(1-\pi)f_u(t)$です。打ち切り患者は治癒者である可能性と、非治癒者のままイベント前である可能性の両方が残るため、$\pi+(1-\pi)S_u(t)$を使います。
この足し算が、通常のパラメトリック生存時間解析との大きな違いです。
ただし、上のコードは尤度の構造を確認するための最小例です。論文用の解析で独自に尤度関数を実装するより、推定法と分散推定が確認されている専用パッケージを使う方を勧めます。特に治癒モデルは尤度面が平坦になりやすく、初期値によって別の解に収束することがあります。
今回のように対数正規モデルで治癒割合が0付近に張り付き、信頼区間が極端に広くなった場合は、数値が返っただけで収束したと判断しません。初期値を変える、尤度を確認する、別の分布でも同じ推定になるかを見る、といった確認が必要です。
論文では治癒割合だけを報告せず、非治癒者側のパラメータと95%信頼区間も併記します。また、追跡期間の中央値だけでは足りません。最長追跡期間と、平坦部が始まった時点以降に何人がリスク集合に残っているかも、治癒割合を読むうえで必要な情報です。
追跡期間の中央値を出す場合は、逆転打ち切りKaplan-Meier法を使うと、単純な観察期間の中央値より追跡状況を表しやすくなります。推定した治癒モデルの曲線をKaplan-Meier曲線に重ねておけば、平坦部だけでなく追跡途中の当てはまりも確認できます。
「治癒」という言葉にも注意が要ります。統計モデル上の治癒割合は、研究で定義したイベントを長期的に経験しない潜在集団の割合であり、そのまま生物学的な治癒を意味するわけではありません。論文では、何をイベントとし、どの意味でcure fractionを定義したのかを明記した方が誤解がありません。
実際のデータでは、平坦部の長さ、残っている患者数、分布を変えた感度分析、競合イベントの有無まで見てから、治癒割合を主要な結果として扱えるかを判断することになります。Dr.データサイエンスでは、こうした治癒モデルを含む生存時間解析の解析方針についてご相談を承っています。
参考文献
1. Boag JW. Maximum likelihood estimates of the proportion of patients cured by cancer therapy. Journal of the Royal Statistical Society: Series B. 1949;11(1):15-53.
2. Farewell VT. The use of mixture models for the analysis of survival data with long-term survivors. Biometrics. 1982;38(4):1041-1046.
3. Sy JP, Taylor JMG. Estimation in a Cox proportional hazards cure model. Biometrics. 2000;56(1):227-236.
4. Lambert PC, Thompson JR, Weston CL, Dickman PW. Estimating and modeling the cure fraction in population-based cancer survival analysis. Biostatistics. 2007;8(3):576-594.
5. Amico M, Van Keilegom I. Cure models in survival analysis. Annual Review of Statistics and Its Application. 2018;5:311-342.
