傾向スコアと生存分析
UPDATE: 2022-10-23 12:26:06
はじめに
このノートではRを使って因果推論に関する基礎的な理論から分析を実行するための方法をまとめている。ここでは傾向スコアと生存分析についてまとめておく。
傾向スコアの数理
傾向スコアの数理については下記が詳しいので、ここでは扱わない。
傾向スコア
傾向スコアとバランシング
傾向スコアとIPW
傾向スコアの層別解析
傾向スコアとマッチング
準備
必要なライブラリを読み込んでおく。
library(tidyverse)
library(survival)
library(survminer)
library(MatchIt)
library(tableone)
library(cobalt)ここで使用するサンプルデータは、みんなの医療統計 多変量解析編で利用されているサンプルデータをお借りする。このデータは右心カテーテル検査(RHC)の有用性に関する論文で使用されたものをサンプリングしたもの。
傾向スコアで有名な研究論文「The effectiveness of right heart catheterization in the initial care of critically ill patients. SUPPORT Investigators.」のデータはこちらよりダウンロードできる。
このデータをもとに、書籍では紹介されていないRでの実行方法をまとめておく。
# ptid: 患者ID
# swang1: 1 is RHC, 0 is NoRHC
# dth30.num: 1 is death in 30days, 0 is alive
# t3d30: died x days entring ICU
# aps1: 集中治療開始時の重篤度スコア
# alb1: 集中治療開始時のアルブミン値
# resp1: 集中治療開始時の呼吸数
# hrt1: 集中治療開始時の心拍数
df <- read_csv('~/Desktop/rhc651.csv', show_col_types = FALSE) %>%
select(ptid, swang1, dth30.num, t3d30, age, aps1, alb1, resp1, hrt1)
df## # A tibble: 651 × 9
## ptid swang1 dth30.num t3d30 age aps1 alb1 resp1 hrt1
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 1 23 60.3 56 2 28 124
## 2 2 1 1 11 66.1 84 2 28 164
## 3 3 1 1 19 72.7 41 2.40 15 135
## 4 4 0 0 30 84.4 43 3.5 10 62
## 5 5 0 0 30 72.3 61 2 47 144
## 6 6 0 0 30 65.1 43 3.5 24 66
## 7 7 0 0 30 77.1 50 3.10 40 60
## 8 8 1 1 3 36.7 74 3.5 14 145
## 9 9 0 0 30 73.0 57 3.5 30 138
## 10 10 1 0 30 41.3 79 3.20 36 120
## # … with 641 more rows
## # ℹ Use `print(n = ...)` to see more rows傾向スコアマッチング
傾向スコアで有名な研究があるので、それをもとに傾向スコアの概念を説明する。ちなみに私は医学系ではないので、細かい内容は知らない。
RESULTS: By case-matching analysis, patients with RHC had an increased 30-day mortality (odds ratio, 1.24; 95% confidence interval, 1.03-1.49). The mean cost (25th, 50th, 75th percentiles) per hospital stay was $49 300 ($17 000, $30 500, $56 600) with RHC and $35 700 ($11 300, $20 600, $39 200) without RHC. Mean length of stay in the ICU was 14.8 (5, 9, 17) days with RHC and 13.0 (4, 7, 14) days without RHC. These findings were all confirmed by multivariable modeling techniques. Subgroup analysis did not reveal any patient group or site for which RHC was associated with improved outcomes. Patients with higher baseline probability of surviving 2 months had the highest relative risk of death following RHC. Sensitivity analysis suggested that a missing covariate would have to increase the risk of death 6-fold and the risk of RHC 6-fold for a true beneficial effect of RHC to be misrepresented as harmful.
CONCLUSION: In this observational study of critically ill patients, after adjustment for treatment selection bias, RHC was associated with increased mortality and increased utilization of resources. The cause of this apparent lack of benefit is unclear. The results of this analysis should be confirmed in other observational studies. These findings justify reconsideration of a randomized controlled trial of RHC and may guide patient selection for such a study.
これはICU入室から24時間以内に右心カテーテルによるモニタリングが行われたグループと行われなかったグループ間で傾向スコアマッチングを行い、30日後死亡率などが比較された研究。結果、右心カテーテルが予後に悪影響(死亡率が上がる)があることがわかり、この結果から、右心カテーテルに対するランダム化比較試験実施する必要性が検討された。
実施にはランダム化比較試験が実施され、同程度の結果となり、右心カテーテルの利用は低下した。という日常的な診療に大きな影響を与えた傾向スコアを使った研究。
ちなみに右心カテーテルとは、近畿中央呼吸器センター診療部によると、
心臓、肺の病気や症状があるとき、特に、心不全や肺高血圧症などの可能性が疑われる場合に受けていただき、診断の確定、今後の治療方針の決定に必要な情報を得る上で有用な検査です。
とのこと。では、傾向スコアの内容についてみていく。下記は、上記の論文The effectiveness of right heart catheterization in the initial care of critically ill patients. SUPPORT Investigators.の患者情報の部分の一部である。他の患者情報はリンク先で確認。
様々な患者の背景情報がまとめられている。年齢、性別、人種、病状などがあるが、RHC(right heart catheterization)を実施したグループと実施していないグループでは、背景情報がバラバラな状態である。これを揃える必要がある。そのために傾向スコアを計算する。
いずれかの治療を選択する場合、個人が一方の治療を受ける傾向をスコア化したのが傾向スコアである。RHCなのかNoRHCを比較する場合、RHCを受ける確率は、個人の背景情報をもとに推測が可能である。例えばRHC=1でNoRHC=0として、交絡の原因となる背景情報を説明変数にして、ロジスティック回帰を実行すれば、RHCを受ける確率が計算できる。つまり、複数の要因が治療に与える影響度合いを1つの傾向スコアという値に集約したことになる。基本的には、傾向スコアが同じ個人を比較すると、背景情報が似ていることになる。
そして、傾向スコアをもとに傾向スコアマッチングを行えば、背景情報が等しい個人をマッチさせることになるので、グループ間の背景情報を均質化できる。この傾向スコアマッチングを行った上で、術後のアウトカムを比較すれば、RHCが予後に影響を与えるかどうかを推測できる。
RHCを受けた患者ではday30の死亡率が高くなっていた(OR1.24[1.03-1.49])。RHCを受けた患者の入院1日当りの費用は49,300ドルであり、受けなかった患者の入院1日当りの費用は35,700ドルであった。RHCを受けた患者さんの ICUの平均滞在期間は 14.8日(5, 9, 17)、受けなかった患者のICUの平均滞在期間は13日であった。
というようにRHCを見直すべきような結果が傾向スコアマッチングの分析から発見された。
傾向スコアと生存分析
まずはバランシングせずにカプランマイヤー曲線を計算する。サンプリングしたデータとはいえ、RHCが行われたグループは術後予後が良くなく、RHCが行われていないグループのほうが1ヶ月以内の生存率が高い。
fit_by_swang <- survfit(Surv(t3d30, dth30.num) ~ swang1, data = df)
survminer::ggsurvplot(fit_by_swang,
fun = 'pct',
risk.table = TRUE,
tables.height = 0.3,
linetype = 'strata',
conf.int = TRUE,
pval = TRUE,
palette = c('#E86670', '#749FC6'),
legend = 'bottom',
break.x.by = 5)
Cox回帰を使って、ハザード比を計算しておく。ハザード比は1.49なので、RHCが行われていないグループを1とすると、RHCが行われたグループは、死亡リスクが1.49倍になることになる。
fit_coxph <- survival::coxph(Surv(t3d30, dth30.num) ~ swang1, data = df, method = 'breslow')
summary(fit_coxph)## Call:
## survival::coxph(formula = Surv(t3d30, dth30.num) ~ swang1, data = df,
## method = "breslow")
##
## n= 651, number of events= 204
##
## coef exp(coef) se(coef) z Pr(>|z|)
## swang1 0.4014 1.4939 0.1405 2.856 0.00429 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## swang1 1.494 0.6694 1.134 1.968
##
## Concordance= 0.552 (se = 0.018 )
## Likelihood ratio test= 8.02 on 1 df, p=0.005
## Wald test = 8.16 on 1 df, p=0.004
## Score (logrank) test = 8.27 on 1 df, p=0.004ここで、各群のベースラインを見てみると、治療開始時の重篤度スコアであるavg_aps1、治療開始時の心拍avg_hrt1に大きな差があるため、交絡が起こっている。そのため、RHCが予後に影響したのかどうかがわからない。つまり、検査を受けた人は、そもそも重篤なことが多く、カテーテル検査が行われたのであれば、単純にカテーテル検査が予後の生存と関わっているとは言えない。
# CreateTableOne(
# data = df,
# vars = c("dth30.num", "t3d30", "age", "aps1", "alb1", "resp1", "hrt1"),
# factorVars = c("dth30.num"),
# strata = "swang1",
# test = FALSE
# )
df %>%
group_by(swang1) %>%
summarise(
n = n(),
avg_age = mean(age),
avg_aps1 = mean(aps1),
avg_alb1 = mean(alb1),
avg_resp1 = mean(resp1),
avg_hrt1 = mean(hrt1)
) %>%
pivot_longer(cols = n:avg_hrt1,
names_to = 'name',
values_to = 'value') %>%
pivot_wider(names_from = swang1, names_prefix = 'flg_') %>%
mutate(diff = abs(flg_1 - flg_0))## # A tibble: 6 × 4
## name flg_0 flg_1 diff
## <chr> <dbl> <dbl> <dbl>
## 1 n 401 250 151
## 2 avg_age 62.1 60.9 1.22
## 3 avg_aps1 50.7 61.4 10.7
## 4 avg_alb1 3.12 2.98 0.141
## 5 avg_resp1 28.9 27.0 1.90
## 6 avg_hrt1 113. 121. 7.57交絡を解消するため傾向スコアマッチングを行う。傾向スコアマッチングはRCTができていない状態のデータにおいて、ベースラインを揃えることで、同じような背景の人のみを解析に入れ込むことで、ランダム化したようなデータが作れる。その状態で介入である治療の有無ごとにアウトカムを比較することで、介入による違いを明らかにできる。
傾向スコアは「各個人の背景情報をもとに、右心カテーテル検査(介入)を受ける傾向(確率)」をスコア化したもの。つまり介入に暴露する確率のこと。傾向スコアを並び替えてざっくりと確認すると、傾向スコアが低いとaps1も低く、傾向スコアが高いとaps1も高くなる。
全分析がそうではないことには注意が必要だが、今回の場合であれば、傾向スコアが重篤度合いを代替表現していると言えるので、傾向スコアで並び替えると、傾向スコアが高い患者ほど重篤で、傾向スコアが低いほど、軽症であるといえる。この状態のデータをバランシングせずに分析すれば、RHCを受けたグループのほうが、重篤な患者が多く偏っていたために、予後が悪くなってしまう。
バランシングすることで、RHCが行われているグループの重篤な患者、RHCが行われていないグループの軽症な患者がマッチせず、RHCが行われているグループの中等症の患者、RHCが行われていないグループの中等症の患者がマッチさせることができる。つまり似たものマッチングすることで、交絡を防げることになる。傾向スコアを使用したからと言って、必ずしも交絡しないわけではなく、必要な共変量が含まれていればの話である。
matching <- matchit(swang1 ~ age + aps1 + alb1 + resp1 + hrt1,
data = df,
method = "nearest")
summary(matching)##
## Call:
## matchit(formula = swang1 ~ age + aps1 + alb1 + resp1 + hrt1,
## data = df, method = "nearest")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.4423 0.3477 0.6066 1.3917 0.1786
## age 60.9120 62.1349 -0.0817 0.6982 0.0349
## aps1 61.4120 50.7406 0.5169 1.1562 0.1065
## alb1 2.9795 3.1206 -0.1936 1.2757 0.0382
## resp1 27.0040 28.9077 -0.1282 1.0579 0.0340
## hrt1 120.5480 112.9751 0.1885 0.9725 0.0447
## eCDF Max
## distance 0.2841
## age 0.1100
## aps1 0.2246
## alb1 0.0964
## resp1 0.1033
## hrt1 0.1031
##
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.4423 0.4054 0.2363 1.4100 0.0514
## age 60.9120 61.3043 -0.0262 0.6585 0.0359
## aps1 61.4120 58.1280 0.1591 1.2801 0.0353
## alb1 2.9795 3.0131 -0.0461 1.1730 0.0200
## resp1 27.0040 28.7600 -0.1182 1.0318 0.0381
## hrt1 120.5480 120.6520 -0.0026 0.9118 0.0163
## eCDF Max Std. Pair Dist.
## distance 0.164 0.2371
## age 0.112 1.2444
## aps1 0.108 0.6268
## alb1 0.064 0.9364
## resp1 0.128 1.0463
## hrt1 0.080 1.0899
##
## Sample Sizes:
## Control Treated
## All 401 250
## Matched 250 250
## Unmatched 151 0
## Discarded 0 0cobaltパッケージのlove.plot()で共変量の調整前後の状態が確認できる。aps1、resp1はうまくバランシングできておらず、\(|0.1|\)を超えている。
love.plot(matching, thresholds = 0.1)
match.data()でマッチングしたデータのみを取り出して、各群の共変量の差を計算すると、さきほどよりも小さくなっている。
matched_df <- match.data(matching)
matched_df %>%
group_by(swang1) %>%
summarise(
n = n(),
avg_age = mean(age),
avg_aps1 = mean(aps1),
avg_alb1 = mean(alb1),
avg_resp1 = mean(resp1),
avg_hrt1 = mean(hrt1)
) %>%
pivot_longer(cols = n:avg_hrt1,
names_to = 'name',
values_to = 'value') %>%
pivot_wider(names_from = swang1, names_prefix = 'flg_') %>%
mutate(diff = abs(flg_1 - flg_0))## # A tibble: 6 × 4
## name flg_0 flg_1 diff
## <chr> <dbl> <dbl> <dbl>
## 1 n 250 250 0
## 2 avg_age 61.3 60.9 0.392
## 3 avg_aps1 58.1 61.4 3.28
## 4 avg_alb1 3.01 2.98 0.0336
## 5 avg_resp1 28.8 27.0 1.76
## 6 avg_hrt1 121. 121. 0.104再度カプランマイヤー曲線を計算する。大きく結果が変わるわけではないが、曲線の差が縮まっていることがわかる。つまり、交絡が取り除かれていることによる影響。
fit_by_matched_swang <- survfit(Surv(t3d30, dth30.num) ~ swang1, data = matched_df)
survminer::ggsurvplot(fit_by_matched_swang,
fun = 'pct',
risk.table = TRUE,
tables.height = 0.3,
linetype = 'strata',
pval = TRUE,
conf.int = TRUE,
palette = c('#E86670', '#749FC6'),
legend = 'bottom',
break.x.by = 5)
Cox回帰のハザード比も1.32となっており、最初の結果よりも死亡リスクが小さい。
fit_coxph_by_matched_swang <- survival::coxph(Surv(t3d30, dth30.num) ~ swang1,
data = matched_df,
method = 'breslow')
summary(fit_coxph_by_matched_swang)## Call:
## survival::coxph(formula = Surv(t3d30, dth30.num) ~ swang1, data = matched_df,
## method = "breslow")
##
## n= 500, number of events= 170
##
## coef exp(coef) se(coef) z Pr(>|z|)
## swang1 0.2840 1.3285 0.1543 1.841 0.0657 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## swang1 1.328 0.7528 0.9818 1.798
##
## Concordance= 0.539 (se = 0.019 )
## Likelihood ratio test= 3.41 on 1 df, p=0.06
## Wald test = 3.39 on 1 df, p=0.07
## Score (logrank) test = 3.41 on 1 df, p=0.06他の変数をモデルに組み込んでみると、ハザード比は1.26となりさらに小さく調整される。いずれにせよ、RHCを行うと死亡リスクが高くなる。
fit_coxph_by_matched_cov <- survival::coxph(Surv(t3d30, dth30.num) ~ swang1 + age + aps1 + alb1 + resp1 + hrt1,
data = matched_df,
method = 'breslow')
summary(fit_coxph_by_matched_cov)## Call:
## survival::coxph(formula = Surv(t3d30, dth30.num) ~ swang1 + age +
## aps1 + alb1 + resp1 + hrt1, data = matched_df, method = "breslow")
##
## n= 500, number of events= 170
##
## coef exp(coef) se(coef) z Pr(>|z|)
## swang1 0.238902 1.269854 0.156390 1.528 0.1266
## age 0.009573 1.009619 0.004824 1.985 0.0472 *
## aps1 0.016430 1.016566 0.004117 3.991 6.58e-05 ***
## alb1 0.009394 1.009438 0.110797 0.085 0.9324
## resp1 0.001705 1.001707 0.005817 0.293 0.7694
## hrt1 -0.002527 0.997476 0.001869 -1.352 0.1764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## swang1 1.2699 0.7875 0.9346 1.725
## age 1.0096 0.9905 1.0001 1.019
## aps1 1.0166 0.9837 1.0084 1.025
## alb1 1.0094 0.9906 0.8124 1.254
## resp1 1.0017 0.9983 0.9904 1.013
## hrt1 0.9975 1.0025 0.9938 1.001
##
## Concordance= 0.621 (se = 0.021 )
## Likelihood ratio test= 26.61 on 6 df, p=2e-04
## Wald test = 26.89 on 6 df, p=2e-04
## Score (logrank) test = 26.85 on 6 df, p=2e-04