競合リスクデータに対する解析

率と割合の1対1対応関係
競合リスクがあると困ること
- 率と割合の1対1対応関係が崩れる
- CIFと1-KMの違い
競合リスクデータに対する解析方針
おわりに
参考資料

以前、競合リスクが存在する生存時間を解析する方法について勉強したのですが、すっかり忘れている（というかあまり深く理解できていなかった）ので、頭の引き出しを整理したいと思います。

ちなみに以前書いた記事はこちら：

necostat.hatenablog.jp

まずは登場する関数と表記方法について。

生存関数（survival function）, S(t)：ある時点tまでにイベントを発生しなかった人の割合。
累積発生関数（cumulative incidence function, CIF）, F(t)：時点tまでにイベントを発生した人の割合。
ハザード関数（hazard function）, h(t)：リスクセットに含まれる対象が、次の瞬間にイベントを発生する条件付き確率密度（瞬間発生率）。
累積ハザード関数（cumulative hazard function）, H(t)：ある時点tまで累積したハザード。ハザード関数を区間0〜tまで積分したもの。

特定のある時点を示したいときは、 $S(i)$ のようにtに代わってiを使います。また、複数あるイベントのうち特定のイベントを示したいときは、 $h_j(t)$ のように添字を用いることにします。

率と割合の1対1対応関係

競合リスクがない場合、ハザード関数の定義

$\begin{aligned} \frac{dF(t)}{dt} &= S(t)h(t) \\ h(t) &= \frac{dF(t)/df}{S(t)}\\ \end{aligned}$

から、以下の関係式が成り立ちます。

$\begin{aligned} F(t) &= \int_0^t S(s) h(s)ds\\ \end{aligned}$

これを（不正確を許容して）日本語で噛み砕くと、

h(s)ds：微小時間 ds において対象集団がイベントを発生する確率
S(s)：時点sまでイベントを発生せず、リスクセットに含まれる確率

なので、その積は微小時間dsの中でイベントが発生する確率になります。この確率を、観察開始（s=0）から時点tまで連続的に足し合わせていけば累積発生割合になる、というわけです。

ここから以下の関係式が導けますし、

$\begin{aligned} S(t) &= e^{-H(t)}\\ \end{aligned}$

S(t) = 1 -F(t)を使って書き直すと、

$\begin{aligned} F(t) &= 1 - e^{-H(t)}\\ \end{aligned}$

となります。累積ハザード関数はハザード関数を積分して得られますので、 発生率（rate）を表すハザード関数と、発生割合（risk）を表す累積発生関数が1対1で対応していることが分かります。

Coxハザードモデルでは、共変量の係数からハザード比（hazard ratio）が推定されますが、ハザード（= 率）と発生割合は1対1対応しているので、ハザードに対する共変量の効果と、発生割合に対する共変量の効果の方向は一致します*1。そのため著者・読者は「率なのか割合なのか」を厳密に区別せず、「危険度が増加・減少する」などのようなやや曖昧な表現を使うことが許容されてい（ると思い）ます。

生存時間解析にまつわる関数同士の関係性はこちらで復習してください。

necostat.hatenablog.jp

競合リスクがあると困ること

2つ以上のイベントを観察していて、あるタイプのイベントが観察されると他のタイプのイベントが観察できなくなるとき、イベントのタイプ同士が発生リスクを競合することになります。例えば、疾患の再発というイベントが、死亡という別のイベントによって観察を妨げられてしまう、といった状況です。

このとき、興味あるイベントの観察を妨げるイベントのことを、競合リスク（competing risk）あるいは競合イベント（competing event）と呼びます。

生存時間解析では「観察打ち切り（censor）」が不可避です。疾患と関係ない理由での脱落や、事前に設定された観察期間終了が理由であれば、打ち切りとイベント観察は独立で無情報であると考えられます（noninformative censoring）。しかし、競合リスクが理由の観察打ち切りであれば無情報とは言えず、「無情報打ち切りの仮定」は成り立ちません。

率と割合の1対1対応関係が崩れる

いま、競合するイベントが2種類あるとしましょう。イベント1の原因別ハザード関数を $h_1(t)$ とするとき、イベント1の累積発生割合は、

$\begin{aligned} F_1(t) &= \int_0^t S(s) h_1(s)ds\\ \end{aligned}$

となります。

このときのリスクセット（= 興味あるイベントが発生しうるとして観察する集団）を「どちらのイベントも発生していない人」とすると、生存関数は、

$\begin{aligned} S(t) &= e^{-H_1(t)-H_2(t)}\\ \end{aligned}$

となります。

ハザード関数はイベント1のことだけに注目しているのに、発生割合（あるいは生存割合）に変換しようとするとイベント2の発生状況も関係してくることがわかります。つまり、ハザード関数（= rate）と累積発生関数（=risk）が1対1対応ではありません。ハザードを対象としたモデルで関連がみられた共変量が、発生割合に対しても同じ方向の効果を持つとは限らない、ということになります。

CIFと1-KMの違い

競合イベントがない場合（例：全死亡）では、カプランマイヤー法で推定した生存時間を1から引いたもの（1-KM）と累積発生関数（CIF）は一致します。

競合イベントがある場合に、興味あるイベント以外のイベントを打ち切り扱いにしていると、1-KMは実際のCIFよりも大きくなります。既に他のイベントを発生してリスクセットから除くべき人を打ち切り扱いにするため、1-KMを用いてCIFを推定すると過大評価することになります。このため、複数のイベントの累積発生割合を1-KMで推定すると、

$\begin{aligned} \Sigma CIF_{KM} \geq 1 - KM_{all} \end{aligned}$

となり、その合計が1を超えてしまう状況が生じる可能性があります。

そこで、前述の原因別ハザード関数を使ってCIFを推定すれば、

$\begin{aligned} \Sigma CIF_{CS} = 1 - KM_{all} \end{aligned}$

となり、整合性が取れます。

競合リスクデータに対する解析方針

競合リスクを含んだデータを解析するとき、モデル化する対象は主に3つ考えられます。

イベント発生のハザード（rate）
イベント発生のリスク（risk）
イベント発生までの時間（time-to-event）

1. ハザードを対象としたモデル

先程説明した「原因別ハザード関数（cause-specific hazard function）」を使ったモデルで解析します。モデルから推定された共変量の係数を指数変換したもの $e^{\beta}$ を、ハザード比（hazard ratio, HR）として解釈する点はCox比例ハザードモデルと同じです。

このモデルでのリスクセットは「まだどのイベントも発生していない人」なので、共変量の係数を「ある時点までどのイベントも発生していなかった人において、瞬間発生率がどれくらい増加・減少するか」と解釈できます。（あとで説明する部分分布ハザードモデルと比べると）この状況を想定することは難しくなく、共変量とイベント1の因果関係として疫学的に解釈することがが容易にできます*2。

しかし、前述のとおりハザードと累積発生割合は1対1対応していません。ハザードを増加（あるいは減少）させる共変量が、発生割合に対しても同様の効果を持つかどうかはわかりません。このため「率（rate）を対象としたモデル（rate-regression model）」と表現してみました*3。

2. リスクを対象としたモデル

「部分分布ハザード関数（sub-distribution hazard function）」を使ったモデルで、提唱者の名前をとってFine-Gray regression modelとも呼ばれます。このモデルでは、次の式で表される部分分布ハザードを対象にモデル化します。

$\begin{aligned} h_1(t) &= \frac{dF_1(t)/df}{1-F_1(t)}\\ \end{aligned}$

先程の原因別ハザードでは、分母が $S(t) = e^{-H_1(t)-H_2(t)}$ でしたが、このモデルでは $1-F(t) = 1 - e^{-H_1(t)}$ を用います。つまり、リスクセットを「他のイベントの発生の有無は問わず、興味あるイベントを発生していない人」と考えます。原因別ハザードでは、興味あるイベント以外のイベントを発生した人をリスクセットから除いたのに対し、部分分布ハザードでは興味の対象のイベントを発生していない限り、これらの人もリスクセットに含める点が異なっています。

前述の部分分布ハザードの定義の両辺を積分して計算すると、

$\begin{aligned} \int \frac{dF_1(t)}{1-F_1(t)} &= \int h_1(s) ds \\ \ln(1 - F_1(t)) &= - \int h_1(s) ds \\ F_1(t) &= 1 - exp \left(-\int h_1(s) ds \right) \end{aligned}$

となり、イベント1についてハザードと累積発生割合が1対1で対応することになります。このモデルで推定された係数βは、累積発生割合との関連をモデル化していると考えられるため、リスク回帰モデル（risk-regression model）と呼ばれることがあります*4。「ハザード比」として記載されていても、実際はリスクとの関連を評価しているので、部分分布ハザード比（sub-distribution hazard ratio, SDHR）と区別する方がいいでしょう。

このモデルでは、実際には興味のイベントを発生しえない人（例：死亡が競合リスクの場合）もリスクセットに含むため、具体的な状況を想定することが難しく、因果推論などの疫学的な解釈には向いていません。代わりに、累積発生割合との関連を直接モデル化するので、発生割合をターゲットにした予測をしたいときに有用です。

3. イベント発生までの時間を対象としたモデル

「対象者が、それぞれのイベント発生までの時間を潜在的に持っているが、実際観察されるのは一部のみ」という考え方でイベント発生までの時間をモデル化する方法です。イベント間での相関関係をコピュラ（copula）という関数で表現します。

Joint frailty-copula modelについては、また別の機会に...（多分）

おわりに

半年ぶりくらいに書きました。競合リスクについて、あまり分かってなかったことが分かりました。
猫の恫喝にも負けず、ちゅ〜るは1週間に2回までにしています。

ランキング参加中

統計

参考資料

Andersen PK, Geskus RB, de Witte T, Putter H. Competing risks in epidemiology: possibilities and pitfalls. Int J Epidemiol. 2012 Jun;41(3):861-70. doi: 10.1093/ije/dyr213. Epub 2012 Jan 9. PMID: 22253319; PMCID: PMC3396320.
Austin PC, Lee DS, Fine JP. Introduction to the Analysis of Survival Data in the Presence of Competing Risks. Circulation. 2016 Feb 9;133(6):601-9. doi: 10.1161/CIRCULATIONAHA.115.017719. PMID: 26858290; PMCID: PMC4741409.
西川正子, 生存時間解析における競合リスクモデル, 計量生物学, 2008, 29 巻, 2 号, p. 141-170
部分分布ハザードモデルの報告方法に関しては、下の文献で解説されています。

pubmed.ncbi.nlm.nih.gov

共変量の係数βの大小比較して、どちらの共変量が発生に対して影響が大きいかは言えるが、大きさの違いなど定量的なことは言えない
異なったイベントに対して推定されたモデルで係数βを比較して、その共変量がどちらのイベントに対して影響が大きいかということは言えない
割合が小さい（<0.2）ときは、complementary log-log linkとlogit linkの差が小さいので、累積発生割合が低いなら $e^{\beta}$ をオッズ比と解釈できる

などがナルホドでした。

*1:定量的にも同じ効果と言えるわけではない（例：ハザード比=2だからといって、リスク比=2と言えるわけではない）

*2:因果関係と解釈する他の要件を満たしているかどうかは、とりあえず置いておいて...

*3:文献に記載があったわけではなく、risk-regression modelに対応させた表現として記載しただけです。

*4:CIFのcomplementary log-log modelの回帰係数と解釈できる