以前、競合するイベントがある場合の生存時間解析として、tidycmprskパッケージを使ってFine-Grayモデルを当てはめる方法を学んだ。
今回はtimeregパッケージを使って実行してみようと思う。
が、解説を読み始めて「そもそもflexible risk regressionって何？？」となったのでそこから勉強...
柔軟なリスク回帰モデル（flexible risk regression model）とは
参考文献*1から抜粋。は累積発生関数（cumulative incidence function, CIF）。
The regression parameters β= (β1, . . . , βK) in the following absolute risk regression model (ARR) have the desired interpretation:
Here F1,0 is an unspecified function of time which represents the cumulative incidence for subjects with X = 0. The logarithmic link function yields that predictor variable changes are translated into ratios of the corresponding cumulative incidences:
（中略）
A corresponding semi-parametric theory is available for the transformation model:
where g is a known differentiable function. The transformation model includes the model (2) as the special case where g(p) = log(p) and β0(t) = log F1,0(t). The logistic-link model corresponds to g(p) = log(p/(1 – p)) and the Fine-Gray model to the complementary log-log link: g(p) = log(–log(p)).
自分なりに噛み砕くと、
- （問題の背景として）ハザードをモデル化するのではなく、累積発生確率（=CIF; ここでは
）をモデル化したい
- 式(2)：「ある患者（あるいは集団）のCIFが、ベースラインCIF（
の何倍か」という形でモデル化してみよう。これをabsolute risk regression model (ARR)と呼ぶよ。
- 式(3)：こうすれば係数βは変数が1単位変化したときの累積発生確率比（=リスク比）になるよね
- 式(4)：もう少し一般化して書いたよ。CIFの変化を時間依存性の部分
と、共変量依存性の部分
に分けて説明するところはCox比例ハザードモデルとおんなじだね。
(4)のg( )は変換の関数で、
- absolute risk regression model：
とすれば(2)が(4)の特殊形だと分かる
- logistic-link model：
としたもの
- Fine-Gray model：
、つまり補対数対数変換したもの。
に対応していると説明されている。
打ち切りがある場合は、打ち切り確率を使った逆確率重み付けを行って推定するみたい。
色々なパッケージを使ってFine-Grayモデルを当てはめる
パッケージとデータセットの準備
データはtidycmprskパッケージに含まれるtrialというシミュレーションデータを使う。 使用する関数によって引数の指定に因子型を要求したり、数値型を要求したりとバラバラなので用意しておく。 使用する関数によっては、欠測が含まれているとエラーが出てしまうようなのであらかじめ除外した（N=200→173になった）。
prodlimパッケージはriskRegression( )で生存時間データを引数として指定するときに必要な関数Hist( )を使うために必要。
library(tidycmprsk) library(timereg) library(tidyverse) library(riskRegression) library(prodlim) trial_clean <- trial %>% mutate(fctr_death = as.factor(.$death), trt_B = if_else(trt=="Drug A",0,1)) %>% drop_na()
> trial_clean # A tibble: 173 × 11 trt age marker stage grade response death death_cr ttdeath fctr_death trt_B <chr> <dbl> <dbl> <fct> <fct> <int> <int> <fct> <dbl> <fct> <dbl> 1 Drug A 23 0.16 T1 II 0 0 censor 24 0 0 2 Drug B 9 1.11 T2 I 1 0 censor 24 0 1 3 Drug A 31 0.277 T1 II 0 0 censor 24 0 0 4 Drug A 51 2.77 T4 III 1 1 death other causes 16.4 1 0 5 Drug B 39 0.613 T4 I 0 1 death from cancer 15.6 1 1 6 Drug A 37 0.354 T1 II 0 0 censor 24 0 0 7 Drug A 32 1.74 T1 I 0 1 death other causes 18.4 1 0 8 Drug A 31 0.144 T1 II 0 0 censor 24 0 0 9 Drug B 34 0.205 T3 I 0 1 death from cancer 10.5 1 1 10 Drug B 42 0.513 T1 III 0 0 censor 24 0 1 # … with 163 more rows
trt：化学療法age：年齢marker：腫瘍マーカーstage：Tステージgrade：グレードresponse：腫瘍の反応性death：死亡death_cr：腫瘍関連死亡、腫瘍以外による死亡、打ち切りttdeath：死亡または打ち切りまでの時間trt_B：0=治療A, 1=治療Bfctr_death：deathを因子型にした。今回は使わないが、通常のCoxハザードモデル当てはめの確認のために作成しただけ。
以下の3つの関数を使って、Fine-Grayモデルの当てはめをやってみる。
- timeregパッケージ comp.risk( )
- riskRegressionパッケージ riskRegression( )
- tidycmprskパッケージ crr( )
方法1：timeregパッケージ comp.risk( )を使う
引数の指定方法の注意点は、以下のとおり。
Event( )：Surv(time, event)と同じだが、この関数ではEvent( )を使う必要あり。cens.code：Event( )の中で使用する。どれが打ち切りを示しているかを指定する。death_crではcensorを1として保持しているのでそのように指定。const( )：時間依存性効果（=時間によって係数が変化する）を許容した解析もできるので、ベースラインから固定された変数を指定するときにはconst( )で囲む。cause：この解析で興味があるイベント。イベントの状態を表すdeath_crの中で、2番目のレベルがdeath from cancerなのでcause=2と指定した。カテゴリー名ではダメみたい。model：="prop"を指定するとFine-Grayモデルになる。="fg"としても同じ。
fit_comprisk <- comp.risk(Event(ttdeath, death_cr, cens.code=1) ~ const(trt)+const(age)+const(marker)+const(stage)+const(grade), data = trial_clean, cause = 2, model = "prop")
結果は以下のとおり。ちょっと推定の問題があるようなコメントがあるが、とりあえずここではパラメトリック項の推定値のみ注目。
> summary(fit_comprisk) Competing risks Model No test for non-parametric terms Did not converge, allow more iterations Parametric terms : Coef. SE Robust SE z P-val lower2.5% upper97.5% const(trt)Drug B 0.9690 0.376 0.376 2.58 0.00990 0.23200 1.7100 const(age) 0.0266 0.012 0.012 2.22 0.02640 0.00308 0.0501 const(marker) -0.2550 0.214 0.214 -1.19 0.23400 -0.67400 0.1640 const(stage)T2 1.1400 0.666 0.666 1.71 0.08710 -0.16500 2.4500 const(stage)T3 0.9960 0.725 0.725 1.37 0.17000 -0.42500 2.4200 const(stage)T4 1.7000 0.638 0.638 2.66 0.00788 0.45000 2.9500 const(grade)II 1.0100 0.484 0.484 2.08 0.03720 0.06140 1.9600 const(grade)III 1.0100 0.455 0.455 2.22 0.02650 0.11800 1.9000
方法2：riskRegressionパッケージ riskRegression( )を使う
引数の指定方法の注意点は、以下のとおり。多分内部でcomp.risk( )を使っている？
Hist( )：Surv(time, event)と同じだが、この関数ではHist( )を使う必要あり。prodlimパッケージが必要。cens.code：Event( )と同様、打ち切りを示すHist( )内で指定する。となっているが、他のレベルを打ち切りに指定しても結果が変わらないので、要らないのかも。cause：この解析で興味があるイベント。今回は数字ではなくラベル名で指定する必要があるみたい。link：="prop"を指定するとFine-Grayモデルになる。
なぜかtrtは受け付けてくれなかった。ラベルに含まれるスペースがよくなかった？
fit_rskreg <- riskRegression(Hist(ttdeath, death_cr, cens.code="censor") ~ trt_B + age + marker + stage + grade, data = trial_clean, cause = "death from cancer", link = "prop")
結果はcomp.risk( )と同じ。
> summary(fit_rskreg) riskRegression: Competing risks regression model IPCW estimation. The weights are based on the Kaplan-Meier estimate for the censoring distribution. Link function: 'proportional' yielding sub-hazard ratios (Fine & Gray 1999), see help(riskRegression). Covariates with time-varying effects: Intercept (numeric) The effects of these variables depend on time.The column 'Intercept' is the baseline risk where all the covariates have value zero (Intercept) 5.3 "0.00015" 13 "0.00209" 17 "0.00409" 20 "0.00594" 24 "0.00828" Shown are selected time points, use plot.riskRegression to investigate the full shape. Covariates with time-constant effects: trt_B (numeric) age (numeric) marker (numeric) stageT2 (numeric) stageT3 (numeric) stageT4 (numeric) gradeII (numeric) gradeIII (numeric) Time constant regression coefficients: Factor Coef exp(Coef) StandardError z CI_95 Pvalue trt_B 0.969 2.634 0.376 2.579 [1.262; 5.499] 0.009902 age 0.0266 1.0269 0.0120 2.2198 [1.003; 1.051] 0.026432 marker -0.255 0.775 0.214 -1.191 [0.510; 1.179] 0.233616 stageT2 1.140 3.127 0.666 1.711 [0.847;11.538] 0.087055 stageT3 0.996 2.708 0.725 1.374 [0.654;11.218] 0.169508 stageT4 1.696 5.454 0.638 2.657 [1.561;19.059] 0.007882 gradeII 1.008 2.739 0.484 2.083 [1.061; 7.069] 0.037232 gradeIII 1.010 2.745 0.455 2.219 [1.125; 6.700] 0.026497 Note: The values exp(Coef) are sub-hazard ratios (Fine & Gray 1999)
方法3：tidycmprskパッケージ crr( )を使う
crr( )の使い方は以前説明したとおり。
failcode：この解析で興味があるイベント。ラベル名で指定する必要あり。- `cencode：打ち切りを表すラベル名を指定。
fit_crr <- crr(Surv(ttdeath, death_cr) ~ trt + age + marker + stage + grade, failcode = "death from cancer", cencode = "censor", data = trial_clean)
で、結果なのだが...
── crr() ────────────────────────────────────────────────────────────────────────────────────────────── • Call Surv(ttdeath, death_cr) ~ trt + age + marker + stage + grade • Failure type of interest "death from cancer" Variable Coef SE HR 95% CI p-value trtDrug B 0.390 0.287 1.48 0.84, 2.59 0.17 age 0.011 0.011 1.01 0.99, 1.03 0.34 marker -0.009 0.162 0.99 0.72, 1.36 0.96 stageT2 0.070 0.494 1.07 0.41, 2.82 0.89 stageT3 0.595 0.476 1.81 0.71, 4.61 0.21 stageT4 1.00 0.426 2.71 1.17, 6.25 0.019 gradeII 0.425 0.384 1.53 0.72, 3.25 0.27 gradeIII 0.484 0.343 1.62 0.83, 3.17 0.16
comp.risk( )およびriskRegressionと一致しない...。 興味あるイベントを指定し間違え？とも思ったり、他のデータセットでも試したけど、やっぱりダメ。競合リスクがない場合はcrr( )とcoxph( )は同じ結果になったんだけど。
そもそもcrr( )の使い方は表面的に見ただけなので、本当にこれでFine-Grayモデルが当てはめられているのか分からない。
参考資料にリンクをつけた文献*2に、次のようなことが書いてあった。
The cumulative incidence function based on the FG model is given by
,
where Λ1(t) is an unknown increasing function and β is a vector of regression coefficients. FG proposed using an inverse probability of censoring weighting technique to estimate β and Λ1(t). This approach is implemented in the crr() function in the cmprsk package (Gray 2010) for R (R Development Core Team 2010).
（中略）
We consider an extended version of the FG model
,
where some effects are proportional as in the FG model (γ) and some effects are allowed to change their effects on the cumulative incidence function over time (α(t)).
2番目のモデル式で、xは効果が時間依存性の変数、zは効果が固定の変数を表していて、x=1とすればFine-Grayモデルになる。そうすると式としては1番目と同じ形になるのだが...。
補対数対数変換について
Fine-Grayモデルは累積発生確率を補対数対数変換（complementary log-log transformation, cloglog）して線形予測子とつなぐモデルなので、補対数対数変換について整理しておく。
変換の目的は、ロジット変換、プロビット変換と同様、あらゆる実数を取りうる変数を、0〜1の範囲*3に収まる変数に変換すること。
- y：-∞〜+∞を取りうる変数
- p：0〜1の範囲しか取らない変数（Fine-Grayモデルでは累積発生確率）
補対数対数変換は、次の数式で表される。
逆関数で書けば、
の形。y=-∞のときはp=0, y=+∞のときはp=1となるような単調増加関数。
似たような変換で、
があって、これもy=-∞のときはp=0, y=+∞のときはp=1となるような単調増加関数。 なんていう名前なんだろうか。
おわりに
- リスク回帰モデルについてはriskRegression( )の使い方をもう少し勉強してみようと思います。
- 我が家の猫は、膝に乗るとおやつがもらえることを学びつつあります。
参考資料
- timeregについていろいろ。
- 統計ERさんの解析。riskRegression( )の存在はこちらで知りました。
timeregパッケージで実行していることを説明。
Scheike TH, Zhang MJ. Analyzing Competing Risk Data Using the R timereg Package. J Stat Softw. 2011 Jan;38(2):i02. PMID: 22707920; PMCID: PMC3375021.Flexible risk regressionについて。解析パッケージの使い方から離れて読めるのでよかった。
Gerds TA, Scheike TH, Andersen PK. Absolute risk regression for competing risks: interpretation, link functions, and prediction. Stat Med. 2012 Dec 20;31(29):3921-30. doi: 10.1002/sim.5459. Epub 2012 Aug 2. PMID: 22865706; PMCID: PMC4547456.補対数対数変換について。
CLOGLOG - Complementary Log-Log Transform
*1:Gerds TA, Scheike TH, Andersen PK. Absolute risk regression for competing risks: interpretation, link functions, and prediction. Stat Med. 2012 Dec 20;31(29):3921-30. doi: 10.1002/sim.5459. Epub 2012 Aug 2. PMID: 22865706; PMCID: PMC4547456.
*2:Scheike TH, Zhang MJ. Analyzing Competing Risk Data Using the R timereg Package. J Stat Softw. 2011 Jan;38(2):i02. PMID: 22707920; PMCID: PMC3375021.
*3:取りうる範囲の境界値を0, 1以外に一般化したものもあるがここでは扱わない
