Joint Frailty-Copula Modelを使った生存時間解析 (2)：実践編（joint.Coxパッケージ） [R]

2つのイベントを同時に考えてLCVを最大化するのは計算負荷が大きいそうで、それぞれのイベントについてLCV1, LCV2を最大にするようにsmoothing parameterを決めます。このときsmoothing parameterを探索する範囲をkappa1, kappa2で指定できますが、正直どれくらいにしたらいいのか分からないのでデフォルトの設定を使うことにします。

では、joint.Coxパッケージに付属の卵巣腫瘍のデータ（dataOvarian）を使って試してみます。

library(tidyverse)
library(joint.Cox)
data(dataOvarian)

サンプルデータの中身：

> dataOvarian %>% glimpse
Rows: 1,003
Columns: 6
$ t.event <int> 1650, 30, 720, 450, 510, 1110, 210, 1380, 1800…
$ event   <dbl> 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0…
$ t.death <int> 1650, 30, 720, 780, 990, 1110, 600, 1410, 1800…
$ death   <dbl> 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0…
$ group   <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
$ CXCL12  <dbl> 1.30594160, 1.28621637, -1.36903145, 1.6132696…

これをjointCox.reg( )に渡します。計算の過程で "Error in integrate..." が出続けますが無視していいとのこと。

result <- jointCox.reg(
  t.event = dataOvarian$t.event,
  event = dataOvarian$event,
  t.death = dataOvarian$t.death,
  death = dataOvarian$death,
  Z1 = dataOvarian$CXCL12,
  Z2 = dataOvarian$CXCL12,
  group = dataOvarian$group,
  alpha = 0
)

出力結果

結果は以下の項目が一括で表示されます。少しずつ見ていきます。

まずは使用したデータのイベント数をまとめたもの。

> result
$count
   No.of samples No.of events No.of deaths No.of censors
4            110           76           46            64
8             58           48           36            22
11           278          185          113           165
14           557          266          290           267

次に、共変量（ここではCXCL12）に対する係数とその信頼区間。

$beta1
  estimate         SE      Lower      Upper 
0.19897418 0.03580646 0.12879353 0.26915484 

$beta2
  estimate         SE      Lower      Upper 
0.16540746 0.04234325 0.08241470 0.24840023

η（イータ）はフレイルティ項の分散の推定値。

$eta
   estimate          SE       Lower       Upper 
0.033275908 0.029043661 0.006014176 0.184112694

θはコピュラのパラメータです。どうやらクレイトンコピュラのみ使用可みたいです。 τはケンドールの順位相関係数で、クレイトンコピュラでは $\tau = \frac{\theta}{|theta + 2}$ で計算されます（この出力でもそのようになっていることを確認）。

$theta
 estimate        SE     Lower     Upper 
2.3632714 0.2354481 1.9440548 2.8728881 

$tau
  estimate     tau_se      Lower      Upper 
0.54162833 0.02473436 0.49290765 0.58956578

LCVを最大にするκとそのときのLCV（だと思います）。

$LCV1
           K1          LCV1 
 3.103448e+16 -4.591487e+03 

$LCV2
           K2          LCV2 
 3.448276e+16 -4.159568e+03

M-スプラインの係数です。gはeventのベースラインハザード関数、hはdeathのそれを表すのに使います（詳細は後ほど）。

$g
[1]  9.136547e-01  1.686345e+00  1.855953e-97  4.447143e-02 8.208009e-155

$h
[1] 2.121898e-01 1.083565e+00 9.906754e-01 1.763836e-01 3.368575e-94

$g_var
              [,1]          [,2]         [,3]          [,4]          [,5]
[1,]  0.0044487471 -0.0027831825  0.00000e+00 -0.0003813422  0.000000e+00
[2,] -0.0027831825  0.0248178870  0.00000e+00  0.0008556739  0.000000e+00
[3,]  0.0000000000  0.0000000000 3.44456e-190  0.0000000000  0.000000e+00
[4,] -0.0003813422  0.0008556739  0.00000e+00  0.0524971946  0.000000e+00
[5,]  0.0000000000  0.0000000000  0.00000e+00  0.0000000000 6.737141e-305

$h_var
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,]  9.506484e-04 -0.0013733137 -3.285041e-06 -2.162385e-04  0.000000e+00
[2,] -1.373314e-03  0.0099792406 -1.177544e-02 -1.839687e-04  0.000000e+00
[3,] -3.285041e-06 -0.0117754352  5.668603e-02 -9.863048e-05  0.000000e+00
[4,] -2.162385e-04 -0.0001839687 -9.863048e-05  4.356045e-02  0.000000e+00
[5,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00 1.134729e-183

最後に収束状況全般について。

$convergence
                 MPL                   DF                  LCV                 code     No.of.iterations No.of.randomizations 
           -8604.369               20.000            -8618.603                1.000               91.000               10.000 

$convergence.parameters
NULL

ベースラインハザード関数を描く

推定されたスプラインの係数（g, h）から、スプライン曲線を描くことができます。

ここで使われているM-スプラインは、負にならないスプラインです。M-スプラインでは、隣接する区域の基底関数を使って再帰的に目的の次元（ここでは3次）の基底関数が作られます。

本を読むと、

$\varepsilon_1$ ：最小のイベント発生時間
$\varepsilon_2$ ： $(\varepsilon_1 + \varepsilon_3)/2$
$\varepsilon_3$ ：最大のフォローアップ期間

の3つが内在するノットに設定されています。

内在ノット数3で3次のMスプラインを定義すると、基底関数は5つ必要になる*2そうで、この基底関数の係数がgやhとして与えられているようです。

joint.CoxパッケージにあるM.spline( )関数を使って、5つの基底関数を求めることができます。

M.spline(time, t_min, t_max)

以下の3つの引数を指定して使います。

time：対象区間の時点をベクトル形式で指定（例：time = c(1,2,3,4,5)など）。t_min, t_maxの範囲に収まっている必要あり。
t_min：時点の下限。non-terminal event発生までの時間の最小値を指定。
t_max：時点の上限。terminal event発生までの時間の最大値を指定*3。

例にならって下のように指定します。

t_min = min(dataOvarian$t.event)
t_max = max(dataOvarian$t.death)

ベースラインハザード関数を書くときは、M.spline( )に実際の時点を与えるのではなく、関数式として定義してcurve( )に渡します。

# 平均線
r1_func = function(t){as.numeric(M.spline(t, t_min, t_max)%*%(result$g))}

# 95%信頼区間の下端
r1_low_func = function(t){
  M_vec = M.spline(t, t_min, t_max)
  r1_V = M_vec%*%(result$g_var)%*%t(M_vec)
  as.numeric(M_vec%*%(result$g) - 1.96*sqrt(diag(r1_V)))
}

# 95%信頼区間の上端
r1_upp_func = function(t){
  M_vec = M.spline(t, t_min, t_max)
  r1_V = M_vec%*%(result$g_var)%*%t(M_vec)
  as.numeric(M_vec%*%(result$g) + 1.96*sqrt(diag(r1_V)))
}

平均の線はM.spline( )で作成した5つの基底関数と、モデルの当てはめで推定したスプラインの係数（g）との積を計算した式を定義しています。

信頼区間の場合は、スプラインの係数の共分散行列（g_var）を使って標準誤差を求めて、mean±1.96*SEとして計算した式を定義しています。

これを下のコードで重ね描きすると、

curve(r1_func, from=t_min, to=t_max)
curve(r1_low_func, from=t_min, to=t_max, lty="dotted", add=TRUE)
curve(r1_upp_func, from=t_min, to=t_max, lty="dotted", add=TRUE)

となります。軸ラベルや範囲は適宜指定してください。

jointCox.Weibull.reg( )でワイブル分布を用いる方法

ベースラインハザード関数にワイブル分布を使う場合は、jointCox.Weibull.reg( )を使います。指定する主な引数はjointCox.reg( )と同じです。

result <- jointCox.Weibull.reg(
  t.event = dataOvarian$t.event,
  event = dataOvarian$event,
  t.death = dataOvarian$t.death,
  death = dataOvarian$death,
  Z1 = dataOvarian$CXCL12,
  Z2 = dataOvarian$CXCL12,
  group = dataOvarian$group,
  alpha = 0
)

出力結果

出力結果を見ていきましょう（一部省略しています）。

Mスプラインを使ったときと、概ね同じようなβが推定されました。フライリティ項の分散は少し大きくなりました。

$beta1
  estimate         SE      Lower      Upper 
0.22058618 0.03912474 0.14390170 0.29727067 

$beta2
  estimate         SE      Lower      Upper 
0.17069401 0.04460975 0.08325891 0.25812912 

$eta
  Estimate         SE      Lower      Upper 
0.09406251 0.06685348 0.02335754 0.37879662

クレイトンコピュラのパラメータθと、そこから求められるケンドールのτも値が少し変わりました。

$theta
 Estimate        SE     Lower     Upper 
1.7017831 0.1951057 1.3592960 2.1305630 

$tau
  Estimate         SE      Lower      Upper 
0.45971983 0.02847593 0.40463716 0.51580451

イベント1, イベント2のそれぞれのベースラインハザードについて、ワイブル分布の尺度パラメータ（scale）と形状パラメータ（shape1）が推定されています。

$scale1
    Estimate           SE        Lower        Upper 
0.0006696365 0.0001630500 0.0004155048 0.0010792008 

$shape1
  Estimate         SE      Lower      Upper 
1.10866686 0.03383253 1.04429928 1.17700187 

$scale2
    Estimate           SE        Lower        Upper 
4.544557e-05 1.563874e-05 2.315104e-05 8.920982e-05 

$shape2
  Estimate         SE      Lower      Upper 
1.31537842 0.04662033 1.22710415 1.41000287

ベースラインハザード関数を描く

推定されたワイブル分布の尺度パラメータと形状パラメータを取り出します。

scale = result$scale1[1]
shape = result$shape1[1]

密度関数にワイブル分布を使ったときのハザード関数は、

$\begin{aligned} h(t) &= \left( \frac{\alpha}{\sigma} \right) \left( \frac{t}{\sigma} \right)^{\alpha - 1} \end{aligned}$

なので、下のコードのようにう関数を記述し、curve( )で描いてみます（dweibull( )とpweibull( )を使っても書けると思いますが、かえってコードが長くなりそうなので）。

r1_func_weibull <- function(t){(shape/scale)*(t/scale)^(shape-1)}
curve(r1_func_weibull, from=t_min, to=t_max)

出来たグラフは下のとおりなんですが...

M-スプラインで描いたときとあまりに違うので、どこか間違えているのかも...。

おわりに

「Rで実行する方法をまとめてみたいと思います！」と意気込んだ割に、コードの意味を表面的になぞるだけになってしまいました。
GJRMというパッケージを使うともう少し自由度が高い joint model も扱えるっぽいです。
NHK「おれ、ねこ」の投稿用写真を選んでます。療養食なので「いつものご飯」と「スペシャルご飯」が同じになってしまう...。

ランキング参加中

統計

参考資料

M-スプラインの話ではないですが、スプライン補間とは何かについて5回シリーズで詳しく解説してくださっています。私のような疫学分野の末端ユーザーでは、局面まで拡大する機会はないでしょうし、難しくて理解しきれなかったところもありますが、とても面白く読ませていただきました。

qiita.com

*1:書籍にはフレイリティのαの記載があったので、おそらくそのことだろう

*2:M-スプラインの基底関数の定義式をみると、当該区域とその隣区域の低次基底関数を使って、次の次数の基底関数を定義しているので、3次まで定義しようとすると5つよりも多くなる気がするんですが、私は理解できてません

*3:3番目の内在ノットは「最大のフォローアップ期間」と記載されているが、例では死亡発生までの時間の最大値を使用しているので、おそらくこれでよいのではないかと思う

ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

Joint Frailty-Copula Modelを使った生存時間解析 (2)：実践編（joint.Coxパッケージ） [R]

ベースラインハザード関数の表し方

jointCox.reg( )でM-スプラインを用いる方法

出力結果

ベースラインハザード関数を描く

jointCox.Weibull.reg( )でワイブル分布を用いる方法

出力結果

ベースラインハザード関数を描く

おわりに

参考資料