かなり久しぶりの記事になりました。世間的にはCOVID-19が流行っておりますが、いかがお過ごしでしょうか。今日テーマにするのは、医学分野でいうと「生存時間解析」という解析です。おそらく、経済学や工学などの解析でも用いられることがある汎用性の高いモデルだと思います。
1. 生存解析（survival analysis）とは
医学的な意義で言えば、旧治療群と新治療群における死亡を比較してその治療法の効果を検証する。というとイメージがつきやすいでしょうか。それでも分からない人は、他に詳しく、かつ丁寧にご説明されているこちらのページをご参照ください1,2。
ご理解できたでしょうか？他にもたくさんの資料があるので今回は割愛。
2. 時間依存性共変量（time-dependent covariates）
生存解析においては、時間がキーワードであることは理解できたと思います。追跡時間を使うということは、ベースライン（初回のデータ収集ポイント）からイベントもしくは打ち切りがあるまでに時間があることになります。つまり、その間に興味のある変数が変わっていますことが大いにあり得ます。例えば、「血圧の心血管疾患罹患への影響を見たい」と思う場合ですが、心血管疾患を発症するまでに血圧のデータはベースラインだけでなく、何度も測定していることが臨床的には多いと思います。つまり、もともとベースラインで高血圧だった人も徐々に改善していることもありますし、逆にベースラインでは正常血圧だった人も徐々に血圧が上がっていたなんていうことも追跡期間には起こります（下図）。
このように、時間と共に変化しうる変数を「時間依存性共変量（time-dependent covariates）」というように呼んでいます。追跡している間に複数回データを収集している場合に有効な解析法となります。
主な事例としては、病院の電子カルテデータ、経済の顧客マーケティングデータ、動物実験の継時的なバイオモニタリングデータ、企業の毎期行われる人事評価データなどでしょうか。
3. R "survival"による実践
じゃあ、この変数をどのように評価するかということですが、今回はRのsurvivalという定番パッケージを使用して実行していきたいと思います。数学的な定義や方法論の制約は、参考図書をご参考にいただければと思います。
ちなみにこちらの英語資料を日本語に直しつつ、アレンジも加えております。
Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model
まずは、使用するデータですが、survivalパッケージ内のPBCデータセットを使用したいと思います。これは、1974年〜1984年まで米国・メイヨークリニックで行われた原発性胆汁性胆管炎に対する治療薬のRCT（ランダム化比較試験）のデータです。418名のデータのうち、初めの312名が研究に同意しているということでそのデータを今回も使います。
| 使用データ | PBC、pbcseq |
|---|---|
| 解析に使う人数 | 312名 |
| 目的変数 | 死亡 |
| 説明変数 | 血清ビリルビン値（連続値）、プロトロンビン時間（連続値） |
| 調整する変数 | 性別、年齢、血清アルブミン値（連続値） |
データを概観するとこんな感じです。
head(PBC) id time status trt age sex ascites hepato spiders edema bili 1 1 400 2 1 58.76523 f 1 1 1 1.0 14.5 2 2 4500 0 1 56.44627 f 0 1 1 0.0 1.1 3 3 1012 2 1 70.07255 m 0 0 0 0.5 1.4 4 4 1925 2 1 54.74059 f 0 1 1 0.5 1.8 5 5 1504 1 2 38.10541 f 0 1 1 0.0 3.4 6 6 2503 2 2 66.25873 f 0 1 0 0.0 0.8 chol albumin copper alk.phos ast trig platelet protime stage 1 261 2.60 156 1718.0 137.95 172 190 12.2 4 2 302 4.14 54 7394.8 113.52 88 221 10.6 3 3 176 3.48 210 516.0 96.10 55 151 12.0 4 4 244 2.54 64 6121.8 60.63 92 183 10.3 4 5 279 3.53 143 671.0 113.15 72 136 10.9 3 6 248 3.98 50 944.0 93.00 63 NA 11.0 3
ここから解析に必要な変数を選定していきます。
# ベースラインのデータになります。IDは312よりも小さい人に限定しています。 temp <- subset(pbc, id <= 312, select=c(id:sex, stage)) head(temp) id time status trt age sex stage 1 1 400 2 1 58.76523 f 4 2 2 4500 0 1 56.44627 f 3 3 3 1012 2 1 70.07255 m 4 4 4 1925 2 1 54.74059 f 4 5 5 1504 1 2 38.10541 f 3 6 6 2503 2 2 66.25873 f 3
このデータを元にして、時間依存性共変量を解析できるようなフォーマットに変換していきます。この時に有用なコマンドがtmergeですので、覚えておきましょう。pbcseqにこの対象者の生化学検査のデータ（血液検査の一部です）（繰り返しあり）があるのでこちらを結合していきます。
pbc2 <- tmerge(temp, temp, id=id, death = event(time, status)) pbc2 <- tmerge(pbc2, pbcseq, id=id, ascites = tdc(day, ascites), bili = tdc(day, bili), albumin = tdc(day, albumin), protime = tdc(day, protime), alk.phos = tdc(day, alk.phos)) head(pbc2) id time status trt age sex stage tstart tstop death ascites bili 1 1 400 2 1 58.76523 f 4 0 192 0 1 14.5 2 1 400 2 1 58.76523 f 4 192 400 2 1 21.3 3 2 4500 0 1 56.44627 f 3 0 182 0 0 1.1 4 2 4500 0 1 56.44627 f 3 182 365 0 0 0.8 5 2 4500 0 1 56.44627 f 3 365 768 0 0 1.0 6 2 4500 0 1 56.44627 f 3 768 1790 0 0 1.9 albumin protime alk.phos 1 2.60 12.2 1718 2 2.94 11.2 1612 3 4.14 10.6 7395 4 3.60 11.0 2107 5 3.55 11.6 1711 6 3.92 10.6 1365
PBCとpbc2は、何が違うかわかりますか？注目するのは、行数です。PBCは一行に一人が与えられていますが、時間依存性共変量の解析に使うpbc2は一個人が数行にまたがっています。これがポイントなのです。つまり、一個人が複数の追跡期間のデータを持っていることになるのです。8、9行目にあるstartとstopという列を見てもらうと分かりますね。
id = 1の人は、トータルで400日間の追跡がありましたが、192日のデータ（0〜192日までとして扱う）、400日のデータ（192〜400日までとして扱う）という2つの行が形成されています。このようにデータを整形することで、一個人でも複数の測定データを使用して解析ができるようになるのです。
A. ベースラインの情報のみモデル
さて、まずはベースラインの情報のみを使用した解析の結果を見てみようと思います。
fit1 <- coxph(Surv(time, status==2) ~ log(bili) + log(protime) + albumin + age + sex, pbc) summary(fit1) ## Call: ## coxph(formula = Surv(time, status == 2) ~ log(bili) + log(protime) + ## albumin + age + sex, data = pbc) ## ## n= 416, number of events= 160 ## (2 observations deleted due to missingness) ## ## coef exp(coef) se(coef) z Pr(>|z|) ## log(bili) 0.893512 2.443696 0.082101 10.883 < 2e-16 *** ## log(protime) 2.764261 15.867304 0.716816 3.856 0.000115 *** ## albumin -0.932484 0.393575 0.198010 -4.709 2.49e-06 *** ## age 0.039990 1.040801 0.007775 5.143 2.70e-07 *** ## sexf -0.120672 0.886325 0.232742 -0.518 0.604124 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## exp(coef) exp(-coef) lower .95 upper .95 ## log(bili) 2.4437 0.40922 2.0805 2.8703 ## log(protime) 15.8673 0.06302 3.8936 64.6634 ## albumin 0.3936 2.54081 0.2670 0.5802 ## age 1.0408 0.96080 1.0251 1.0568 ## sexf 0.8863 1.12825 0.5617 1.3986
B. 時間依存性共変量としたモデル
次に、ビリルビン値とプロトロンビン時間を時間依存性共変量として組み込んだモデルでは、どうなるか見てみます。
fit2 <- coxph(Surv(tstart, tstop, death==2) ~ log(bili) + log(protime) + albumin + age + sex, pbc2) summary(fit2) ## Call: ## coxph(formula = Surv(tstart, tstop, death == 2) ~ log(bili) + ## log(protime) + albumin + age + sex, data = pbc2) ## ## n= 1807, number of events= 125 ## ## coef exp(coef) se(coef) z Pr(>|z|) ## log(bili) 1.21547 3.37187 0.12493 9.729 < 2e-16 *** ## log(protime) 2.94747 19.05760 0.62497 4.716 2.40e-06 *** ## albumin -1.58789 0.20436 0.21013 -7.557 4.13e-14 *** ## age 0.04298 1.04392 0.01022 4.205 2.61e-05 *** ## sexf -0.05912 0.94259 0.25681 -0.230 0.818 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## exp(coef) exp(-coef) lower .95 upper .95 ## log(bili) 3.3719 0.29657 2.6395 4.3074 ## log(protime) 19.0576 0.05247 5.5988 64.8702 ## albumin 0.2044 4.89344 0.1354 0.3085 ## age 1.0439 0.95793 1.0232 1.0650 ## sexf 0.9426 1.06090 0.5698 1.5593
4. 結果
この結果から、log(bili)については、exp(coef)がそれぞれ2.44と3.37でした。また、log(protime)については、それぞれ15.87と19.06でした。どちらも時間共変量を考慮したモデルの方が推定値が大きかったということになります。
この参考にした文書によると、この疾患自体が段階的に肝臓の機能低下をもたらすことから、肝臓の機能指標であるこれらの値が高くなることは、PBCによる死亡リスクと正の関連を示しており、時間依存性共変量を考慮した方がより大きい推定値を示したことが考えられます。
5. 最後に
今回使用したケースでは、病態と関連し徐々に一定方向へ変化する指標を時間依存性共変量として組み込んだことでこのような結果をもたらしたと考えられます。そのため、どの場合でもこの解析が有効とは限りませんが、継時的なデータを持っている場合には、単純にベースライン値のみを使用した解析だけでなく、このような解析にも試みるとまた違う結果・考察をえるかもしれません。
6. 参考ページ・図書
1. Therneau T et al. Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model
https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf
2. 和歌山県立医大・下川先生による医学統計セミナー
https://waidai-csc.jp/updata/2019/05/20191114_%E7%AC%AC7%E5%9B%9E_%E8%B3%87%E6%96%99.pdf
3. Kooさんによる「生存解析のすヽめ：カプランマイヤー法とコックス比例ハザードモデル」
https://www.medi-08-data-06.work/entry/surv_analysis
4. Survival Analysis with R
https://rviews.rstudio.com/2017/09/25/survival-analysis-with-r/
5. Data Blueさんによる「Rでの生存時間解析（時間変化共変量の処理）」
Rでの生存時間解析（時間変化共変量の処理） - Data Blue
20200319
RF
