R package "JM"を使用してJoint modelを実装する。 - 統計学と疫学と時々、助教生活

前回は生存解析においてベースラインだけでなく、縦断的に測定している（繰り返し測定、反復測定値ともいう）値の影響も組み込んだ「時間依存性共変量」についてまとめました。今回は、さらに進化した（個人的見解です）Joint modelについて、簡単な紹介とRでの実装を試みたいと思います。

1. 生存時間解析や縦断解析とは
2. Joint modelとは何を「Joint」しているのか
3. Joint modelの数学的な発想
4. R package "JM"を使用したJoint modelの実装
5. Joint model の実装パッケージ紹介
6. まとめ
7. 参考ページ・図書

1. 生存時間解析や縦断解析とは

まず、生存時間解析は前回も引用したのですが、特定の曝露要因の有無と何かのイベント（死亡や発症）との関連について時間の概念も加えた解析になります。縦断解析ですが、今回は、同じ個人で繰り返し同じ項目について測定しているデータについて扱っているものを指しています（場合によっては、それとは異なるものをイメージすることもありますので、言葉のチョイスにはご注意）。

2. Joint modelとは何を「Joint」しているのか

Joint modelの本質でもありますが、このモデルで対象としているのが、大きく分けて2つのポイントであると論文には記載してあります。

繰り返し測定したデータの混合効果について
時間を含んだアウトカム（疾患の発症や死亡など）への曝露要因の効果について

Joint modelではこの両者の推定を同時に行うことが最大の売りとなっています。
対象者個人について、縦断的に測定しているデータ（血液中のバイオマーカーなど）とそのアウトカムデータを保有していれば、この解析が実行できます。すなわち、time-dependent covariate analysisで扱ったデータと同じですので、Joint modelと比較するときに、time-dependent covariateを使用したCox比例ハザードモデルがよく出てきます。

3. Joint modelの数学的な発想

難しいので、飛ばしたい人は4.の実装に行ってください。

さて、2で説明した内容を数式を使用して説明します。まず1つ目の式、 $Y_{ij}$ というのは、 $i$ 番目の個人の（たくさんらう曝露要因のうちの） $j$ 番目の曝露要因です。そして、この値というのは誤差を含んでいます。 $Y^{*}_{i} (t_{ij})$ というのは、 $t$ 時点の $Y_{ij}$ の本来の値で、 $Z_{ij}$ は測定誤差を示しています。さらに、 $Y^{*}_{i} (t_{ij})$ を2つの部分に分解すると、固定効果の $X_{ij}{\beta}$ と変量効果の $W_{ij}{\beta}$ になります。この先は、混合効果モデルの説明に譲りたいと思います。とにかく1つ目の式は、曝露要因の変量効果、測定誤差を考慮したモデルになっていることが理解できればOKです。

$Y_{ij} = Y^{*}_{i} (t_{ij}) + Z_{ij} = X_{ij}{\beta} + W_{ij}{\beta} + Z_{ij}$

さて、2つ目の式は、時間依存性共変量の ${\lambda}_i ( t | K_i ) = {\lambda}_0(t) exp ( K_i{\alpha}_1)$ が分かっている人には非常に簡単かもしれませんが、1式で求めた $Y^{*}_{i}(t)$ を組み込んだだけです。つまり、 ${\gamma}_2$ は、測定誤差を除いた $Y^{*}_{i}(t)$ とイベントとの関連を示したものということになります。

${\lambda}_i \Bigl( t | K_i, Y^{*}_{i}(t) \Bigr) = {\lambda}_0(t) exp \Bigl( K_i{\gamma}_1 + Y^{*}_{i}(t){\gamma}_2 \Bigr)$

再度まとめますと、この2つの式で与えられたバラメータを同時に推定するJoint modelは、まず全ての時間 $t$ における $Y^{*}_{i}(t)$ を推定することにより、 ${\gamma}_2$ を推定することを可能にしています。分かりにくいという方は、Asar öらのInt J Epi, 2015;44:334-344をご参照下さい。

4. R package "JM"を使用したJoint modelの実装

1. データの準備と研究の目的

今回の目的は、薬剤2剤（drugのddCおよびddI）の死亡リスクへの影響です。

データは、下記の2つのデータセットを使用します。

縦断解析用：「aids」という1405行（467人の個人について）×12変数のデータ（繰り返しデータ）
生存解析用：「aids.id」という467行×12変数のデータ（ベースラインの情報のみ）

library(JM) #パッケージの読み込み 
head(aids) 
  patient  Time death       CD4 obstime drug gender prevOI         AZT start  stop event
1       1 16.97     0 10.677078       0  ddC   male   AIDS intolerance     0  6.00     0
2       1 16.97     0  8.426150       6  ddC   male   AIDS intolerance     6 12.00     0
3       1 16.97     0  9.433981      12  ddC   male   AIDS intolerance    12 16.97     0
4       2 19.00     0  6.324555       0  ddI   male noAIDS intolerance     0  6.00     0
5       2 19.00     0  8.124038       6  ddI   male noAIDS intolerance     6 12.00     0
6       2 19.00     0  4.582576      12  ddI   male noAIDS intolerance    12 18.00     0

head(aids.id)
  patient  Time death       CD4 obstime drug gender prevOI         AZT start stop event
1       1 16.97     0 10.677078       0  ddC   male   AIDS intolerance     0  6.0     0
2       2 19.00     0  6.324555       0  ddI   male noAIDS intolerance     0  6.0     0
3       3 18.53     1  3.464102       0  ddI female   AIDS intolerance     0  2.0     0
4       4 12.70     0  3.872983       0  ddC   male   AIDS     failure     0  2.0     0
5       5 15.13     0  7.280110       0  ddI   male   AIDS     failure     0  2.0     0
6       6  1.90     1  4.582576       0  ddC female   AIDS     failure     0  1.9     1

2. 比較のために、time-dependent covariateを組み込んだCox proportional hazard modelを実行する。

td.cox <- coxph(Surv(start, stop, event) ~ drug + sqrt(CD4), data = aids)
summary(td.cox)
Call:
coxph(formula = Surv(start, stop, event) ~ drug + sqrt(CD4), 
    data = aids)

  n= 1405, number of events= 188 

              coef exp(coef) se(coef)      z Pr(>|z|)    
drugddI    0.32678   1.38650  0.14708  2.222   0.0263 *  
sqrt(CD4) -0.72302   0.48528  0.07997 -9.042   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
drugddI      1.3865     0.7212    1.0393    1.8498
sqrt(CD4)    0.4853     2.0606    0.4149    0.5676

Concordance= 0.696  (se = 0.018 )
Likelihood ratio test= 86.14  on 2 df,   p=<2e-16
Wald test            = 83.51  on 2 df,   p=<2e-16
Score (logrank) test = 83.25  on 2 df,   p=<2e-16

3. 今回の主題であるJoint modelを実装する（下記の順に従う）。

lme関数を使用して、観察時間がCD4陽性細胞に与える影響を混合効果モデルで推定する（fitLMEに格納）。
Surv関数を使用して、薬剤と死亡までの時間をモデリングする（fitSURVに格納）。
JMパッケージのjointModel関数を使用して、1と2の結果からパラメータ推定を行う（fit.JMに格納）。

fitLME <- lme(sqrt(CD4) ~ obstime + obstime:drug, random = ~ obstime | patient, data = aids)
fitSURV <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
fit.JM <- jointModel(fitLME, fitSURV, timeVar = "obstime", method = "piecewise-PH-GH")
summary(fit.JM)

Call:
jointModel(lmeObject = fitLME, survObject = fitSURV, timeVar = "obstime", 
    method = "piecewise-PH-GH")

Data Descriptives:
Longitudinal Process		Event Process
Number of Observations: 1405	Number of Events: 188 (40.3%)
Number of Groups: 467

Joint Model Summary:
Longitudinal Process: Linear mixed-effects model
Event Process: Relative risk model with piecewise-constant
		baseline risk function
Parameterization: Time-dependent 

   log.Lik      AIC      BIC
 -2107.647 4247.295 4313.636

Variance Components:
             StdDev    Corr
(Intercept)  0.8660  (Intr)
obstime      0.0388  0.0680
Residual     0.3754        

Coefficients:
Longitudinal Process
                  Value Std.Err z-value p-value
(Intercept)      2.5558  0.0372 68.7961 <0.0001
obstime         -0.0423  0.0046 -9.1931 <0.0001
obstime:drugddI  0.0051  0.0065  0.7821  0.4342

Event Process
            Value Std.Err z-value p-value
drugddI    0.3511  0.1537  2.2839  0.0224
Assoct    -1.1016  0.1180 -9.3388 <0.0001
log(xi.1) -1.6489  0.2498 -6.6000        
log(xi.2) -1.3393  0.2394 -5.5940        
log(xi.3) -1.0231  0.2861 -3.5758        
log(xi.4) -1.5802  0.3736 -4.2299        
log(xi.5) -1.4722  0.3500 -4.2069        
log(xi.6) -1.4383  0.4283 -3.3584        
log(xi.7) -1.4780  0.5455 -2.7094        

Integration:
method: Gauss-Hermite
quadrature points: 15 

Optimization:
Convergence: 0

4. 結果の解釈

薬剤の影響についてのcoefficientですが、0.33から0.35へと少し強い関連になっています。注目すべきは、CD4陽性細胞の影響についてですが、-0.72から-1.10へとかなり大きく変化しています。CD4陽性細胞は、測定誤差を多く含んでいるため、それらを組み込んでいる今回のモデルの方が時間依存性共変量を組み込んだモデルより絶対値が大きくなっていると考えられます。

DIsadvantageですが、一点。
非常に計算能力を要するので、コンピュータの性能に依存しますが、サンプルデータですら少々時間がかかります。さらに、データの行数（個人と測定回数）が増えると、より時間がかかると予想されますので、要注意です。

5. Joint model の実装パッケージ紹介

今回は、Joint modelについて本も執筆されているDr.RIzopoulos DのJMを用いましたが、様々な統計ソフトおよびパッケージで提供されています。長所や短所については、使用していないので把握しておりません。

R - JMBayes (https://cran.r-project.org/web/packages/JMbayes/JMbayes.pdf)
R - joineR (https://cran.r-project.org/web/packages/joineR/joineR.pdf)
SAS - JMfit (https://www.jstatsoft.org/article/view/v071i03)
SAS - JM (https://www.jm-macro.com/)
Stata - STJM (https://www.stata-journal.com/article.html?article=st0289)

6. まとめ

これからのデータサンプリング方法（よりリアルタイムになる）を考えると、このようなモデリングというのは非常に魅力的だなと思いました。とりわけ、私が扱っている血液データというのは測定誤差が生じやすいので、今後研究での実装も早々に検討していきたいところです。縦断的なデータかつイベントまでの時間データを持っている場合には、前回紹介した時間依存性共変量と今回のJoint modelを実装してみてはいかがでしょうか。何かおかしなポイントや間違っているところがあれば、お気軽にコメントください。

7. 参考ページ・図書

なかなかこの内容だけでは、理解できないと思いますので、下記の文献を繰り返し読む事をお勧めします。

1. JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data.
https://www.jstatsoft.org/article/view/v035i09

2. Joint modelling of repeated measurement and time-to-event data: an introductory tutorial.
https://academic.oup.com/ije/article/44/1/334/657852

3. Red Blood Cell Distribution Width Is Longitudinally Associated With Mortality and Anemia in Heart Failure Patients.
https://www.jstage.jst.go.jp/article/circj/78/2/78_CJ-13-0630/_article/-char/ja

4. Survival Analysis with R
https://rviews.rstudio.com/2017/09/25/survival-analysis-with-r/

5. 混合効果モデルのイメージについて（An Introduction to Hierarchical Modeling）
http://mfviz.com/hierarchical-models/

20200428
RF