Joint Frailty-Copula Modelを使った生存時間解析 (1)：理論編

同一施設の患者は類似した傾向を持っている可能性があります。同一のクラスターに属する患者が共有する共通性を、フレイルティ（frailty）項*1としてモデルに含めることができます。これは一般化線形モデルでクラスター効果を考慮するときに導入される変量効果（random-effect）と同じようなものです。

比例ハザードモデルでは、個々の患者のハザード関数 $\lambda_{j}$ は下のように、ベースラインハザード $\lambda_0(t)$ と患者の特性によって決まる係数 $\exp(Z_j \beta)$ の積として表されます。「患者の特性によって、その患者のハザードがベースライン（= 全ての共変量が0のときのハザード）の何倍になるか」を決めているということになります。

$\begin{aligned} \lambda_{j} (t | Z_j) &= \lambda_0(t) \exp(Z_j \beta) \end{aligned}$

ここで、患者が属するクラスターによって決まるフレイルティ項 $b_i$ （ $b_i > 0$ ）を導入します。

$\begin{aligned} \lambda_{ij} (t | Z_{ij}) &= b_i \lambda_0(t) \exp(Z_{ij} \beta) \end{aligned}$

ハザードが高いクラスターでは $b_i$ > 1、低いクラスターでは 0 < $b_i$ < 1となります。同じクラスターに属する人は、同じフレイルティを共有するので、類似した傾向を持つことになるわけです。

また、混合効果モデルと同様、各クラスターの $b_i$ は直接推定されるのではなく、設定した確率分布に従う潜在変数として扱われます*2。混合効果モデルではクラスター効果は正規分布に従うと仮定されることが多いですが、フレイルティ $b_i$ ではガンマ分布が用いられることが多いとのこと。

コピュラ（copula）で2つの生存関数の相関をモデル化する

医学研究では1つの対象者について複数のイベントを測定することが珍しくありません。このようなデータを解析するとき、コピュラを用いて生存関数の依存関係（= 相関）をモデル化することで競合リスクを考慮することができます。

コピュラ（copula）, $C_{\theta}$ とは、2つの生存関数を結びつけ、1つの結合分布として扱うために用いられる関数です。添字にあるパラメータθは2つの変数間の依存度の強さを表すために用いられます。

よく用いられるコピュラはアルキメデス型コピュラと呼ばれるもので、[0, 1] の範囲を[0, ∞] に変換する生成関数（generator function） $\Phi_{\theta}$ とその逆関数 $\Phi_{\theta}^{-1}$ を用いて、2つの生存関数

$\begin{aligned} u &= Pr(X>x|Z)\\ v &= Pr(Y>y|Z) \end{aligned}$

を次のように結合させます。

$\begin{aligned} Pr(X>x, Y>y |Z) &= C_{\theta} (u, v) \\ &= \Phi_{\theta}^{-1} \left( \Phi_{\theta}(u) + \Phi_{\theta}(v)\right) \end{aligned}$

生存時間解析でよく使われるのがクレイトンコピュラ（Clayton copula）で、その生成関数および逆生成関数は、

$\begin{aligned} \Phi_{\theta} (t) &= \frac{1}{\theta} (t^{-\theta} - 1) \\ \Phi_{\theta}^{-1} (s) &= (1 + \theta s)^{-\frac{1}{\theta}} \end{aligned}$

となります（0<t<1, s>0, θ>0）。

ちなみに、それぞれの定義域・値域を確認すると以下のとおり：

t→0のとき $\Phi_{\theta} (t)$ → ∞
t→1のとき $\Phi_{\theta} (t)$ → 0
s→0のとき $\Phi_{\theta}^{-1} (s)$ → 1
s→∞のとき $\Phi_{\theta}^{-1} (s)$ → 0

前述の生成関数と逆生成関数から、クレイトンコピュラの式は下のようになることが分かります。

$\begin{aligned} C_{\theta} (u, v) &= \Phi_{\theta}^{-1} \left( \frac{1}{\theta} (u^{-\theta} - 1) + \frac{1}{\theta} (v^{-\theta} - 1)\right) \\ &= \Phi_{\theta}^{-1} \left( \frac{1}{\theta} (u^{-\theta} + v^{-\theta} - 2) \right) \\ &= \left( 1 + \theta \frac{1}{\theta} (u^{-\theta} + v^{-\theta} - 2) \right)^{-\frac{1}{\theta}} \\ &= \left(u^{-\theta} + v^{-\theta} - 1 \right)^{-\frac{1}{\theta}} \\ \end{aligned}$

依存度の指標1：ケンドールの順位相関係数

生存関数間の依存度を表す1つ目の指標は、ケンドールの順位相関係数（Kendall's $\tau$ ）です。

$\begin{aligned} \tau &= Pr \{(X_1 - X_2)(Y_1 - Y_2)>0\} - Pr \{(X_1 - X_2)(Y_1 - Y_2)<0\} \end{aligned}$

τは2つの生存関数からランダムに抽出した生存時間（X1, X2とY1, Y2）の順位が一致する確率から、一致しない確率を引いたもので、 順位が完全に一致、つまり完全に正の相関があればτ=1となり、順位が完全に逆であればτ=-1となります。

τはθによって決まり、クレイトンコピュラでは

$\begin{aligned} \tau &= \frac{\theta}{\theta + 2} \end{aligned}$

の関係が成り立ちます（θ=0に近いほど依存度が低く、θが大きくなるほど依存度が高くなる）。

依存度の指標2：クロスレシオ関数

2つ目の指標はクロスレシオ関数*3（cross-ratio function, CRF）です。こちらもパラメータθの関数になるので、θを添えて $R_{\theta}$ と書くみたいです。

CRFは下の式のように、コピュラ関数とその2階微分、偏微分で定義されます。

$\begin{aligned} R_\theta (u, v) &= \frac{\frac{\partial^2 C_\theta(u,v)}{\partial u \partial v}}{\frac{\partial C_\theta(u,v)}{\partial u} \frac{\partial C_\theta(u,v)}{\partial v}}C_\theta(u,v) \end{aligned}$

2階微分は同時分布の密度、偏微分はそれぞれの生存確率の周辺分布の密度を表しています（多分）。2つの変数が独立なら、同時分布は周辺分布の積になるので、「周辺分布に比べて、同時分布が大きいなら相関がある」という解釈になるのかと…。コピュラ関数をかけているのは標準化のためでしょうか。

一般的にCRFは生存確率の関数なので、「ある時点での局所依存度」を表しますが、 クレイトンコピュラでは $R_{\theta} = 1 + \theta$ であり、生存時間に寄らずCRFが一定になっています。

解釈は以下のとおり：

$R_{\theta}$ > 1のとき、局所における正の相関
0 < $R_{\theta}$ < 1のとき、局所における負の相関
$R_{\theta}$ = 1のとき、局所における独立

もう少し実践的には、「ある時点xでイベント1を発生した人は、その時点まで発生しなかった人に比べて、イベント2ののハザードが何倍か」と解釈できます。条件付きのハザード比・オッズ比として解釈することの詳細は教科書などを読んでみてください。

Joint Frailty-Copula Model

フライリティとコピュラを組み込んだ結合分布モデルを数式で書いておきます。

まず、2つのイベントのハザードは、

$\begin{aligned} \lambda_{xij} (t | b_i) &= b_i \lambda_{x0}(t) \exp(Z_1 \beta_1) \\ \lambda_{yij} (t | b_i) &= b_i^{\alpha} \lambda_{y0}(t) \exp(Z_1 \beta_1) \\ \end{aligned}$