Zenn
🦔

Rによる生存時間解析

2022/02/25に公開
2件
R
統計
生存時間解析
survivalanalysis
生物統計
tech

はじめに

本記事では、医療分野における生存時間解析の基本的な解析のRでの実装を紹介する。生物統計の初学者よりは一定の知識がある者を対象とし数式やRの使用方法に関する記述を省いたためわかりづらい箇所があるかもしれないが、ご容赦いただければ幸いである。そのため生存時間解析をはじめて学習する方はチュートリアル論文や教科書で先に学ぶことを推奨する。

【参考文献】

【変更履歴】
2022/11/9 対比検定のc2の係数およびその結果を変更
2022/11/18 効果量の推定時の変数をログランク検定の章と整合性が取れるように修正

目次

  • Package
  • Byar & Greene prostate cancer data
    • Variables
    • データの読み込み
  • データの要約
  • 生存時間の記述
    • カプランマイヤー曲線
    • Median 生存時間
  • ログランク検定
    • ログランク検定（層別、重み付き）
    • 対比検定
  • 効果量の推定
    • Cox回帰
    • RMST
  • サンプルサイズ計算
    • 基本のサンプルサイズ計算
    • 登録期間を考慮したサンプルサイズ計算
    • 脱落を考慮したサンプルサイズ計算
    • （補足）HRの推定
  • まとめ
    • （メモ）今後、追加する可能性のある解析

Package

利用するPackageは以下の通り。

  • tidyverse
  • gtsummary
  • survival
  • survminer
  • MASS
  • survRM2
  • gsDesign
#package
library(tidyverse)
library(gtsummary)
library(survival)
library(survminer)
library(MASS) 
library(survRM2)
library(gsDesign)

Byar & Greene prostate cancer data

前立腺がん治療のランダム化比較試験（Byar and Green, 1980）のデータを利用した。前立腺がんのStage3、Stage4の患者を対象にエストロゲン治療群（0.2 mg, 1.0 mg, 5.0 mg）とPlacebo群の4群に割付が行われた。データはVanderbilt Biostatistics Datasetsの"Byar & Greene prostate cancer data"からダウンロードが可能である（2022-2-25時点）。

Variables

・データに含まれている変数のうち利用した変数

変数名 概要
age 年齢
stage Stage(3 or 4)
rx 割付群（placebo, 0.2 mg estrogen, 1.0 mg estrogen or 5.0 mg estrogen）
dtime 追跡期間（月）
status alive もしくは 死因


・新たに作成した変数

変数名 概要
rx_c rxを数値化
rx_c2 placebo群（0）とエストロゲン群（1）
status2 生存、前立腺がんによる死亡、心疾患による死亡、その他に分類
status3-5 status2のダミー変数（status5が前立腺がんによる死亡）
stage_c StageをFactorにした変数
dtime_y 追跡期間（年）

データの読み込み

#データの読み込み
prostate <- read.table("path/prostate.csv",
                   sep=",",
                   header=TRUE,
                   na.strings=c(''),
                   fill=TRUE)

#変数の作成
th1 <- prostate %>% 
  mutate(rx_c = case_when(
    rx == "placebo" ~ 0,
    rx == "0.2 mg estrogen" ~ 1,
    rx == "1.0 mg estrogen" ~ 2,
    rx == "5.0 mg estrogen" ~ 3
  )) %>% 
  mutate(rx_c2 = case_when(
    rx == "placebo" ~ 0,
    rx == "0.2 mg estrogen" ~ 1,
    rx == "1.0 mg estrogen" ~ 1,
    rx == "5.0 mg estrogen" ~ 1
  )) %>% 
  mutate(status2 = case_when(
    status == "alive" ~ 0,
    status == "dead - prostatic ca" ~ 1,
    status %in% c("dead - heart or vascular","dead - cerebrovascular") ~ 2,
    TRUE ~ 3
  )) %>% 
  mutate(status3 =
           ifelse(status2==0,0,
                  ifelse(status2==2,1,0))) %>% 
  mutate(status4 = 
           ifelse(status2==0,0,
                  ifelse(status2==3,1,0))) %>% 
  mutate(status5 = 
           ifelse(status2==1,1,0)) %>% 
  mutate(stage_c=as.factor(stage)) %>% 
  mutate(dtime_y=dtime/12)

データの要約

治療群ごとに、年齢、Stage、Status、追跡期間を集計した。0から順にplacebo, 0.2 mg estrogen, 1.0 mg estrogen, 5.0 mg estrogen群である。

th2 <- th1 %>% dplyr::select(rx_c,age,stage_c,status,dtime_y,status5) 
  
table1 <- tbl_summary(
  th2,
  by = rx_c,
  missing = "ifany") %>% 
  add_n() %>% 
  modify_header(label = "**Variable**")

生存時間の記述

カプランマイヤー曲線

生存時間解析でよく記述されるカプランマイヤー曲線をPlaceboとEstrogenの2群で作成した。また2年おきのAt Riskの人数も算出し、95％信頼区間もつけて図を作成した。補足として、生存割合の信頼区間はGreenwoodの公式から分散を計算して算出されるが、変数変換を行わない場合には生存割合として不適切な値をとり得るためdelta法によりLog変換やLog-Log変換を行い0から1に収まるようにした上での信頼区間がしばしば計算される。survfitではデフォルトではLog変換が行われているようである。[1]

#KM曲線
fit <- survfit( Surv(dtime_y, status5) ~ rx_c2, data = th1 )

ggsurvplot(
  fit,                     
  data = th1,             
  risk.table = TRUE,       
  pval = TRUE,             
  conf.int = TRUE,         
  xlab = "YEAR", 
  ggtheme = theme_light(), 
  risk.table.y.text.col = T, 
  risk.table.y.text = FALSE, 
  legend.labs =
   c("Placebo", "Estrogen")
)

カプランマイヤー曲線をPlacebo, 0.2 mg estrogen, 1.0 mg estrogen, 5.0 mg estrogenの4群ごとに作成した。

fit2 <- survfit( Surv(dtime_y, status5) ~ rx_c, data = th1 )

ggsurvplot(
  fit2,                     
  data = th1,             
  risk.table = TRUE,       
  pval = TRUE,             
  conf.int = FALSE,         
  xlab = "YEAR", 
  ggtheme = theme_light(), 
  risk.table.y.text.col = T, 
  risk.table.y.text = FALSE, 
  legend.labs =
    c("Placebo", "0.2 mg estrogen", "1.0 mg estrogen", "5.0 mg estrogen")
)

median 生存時間

PlaceboとEstrogen群の2群に分けた場合での、各群におけるmedian生存時間を算出した。ただしEstorogen群では観察期間内に生存割合が50%を下回っておらず算出できなかった。
Placebo群のmedian 生存時間は5.92年となった。

#median survival time
survfit(Surv(dtime_y, status5) ~ rx_c2, data = th1)

## Call: survfit(formula = Surv(dtime_y, status5) ~ rx_c2, data = th1)
## 
##           n events median 0.95LCL 0.95UCL
## rx_c2=0 127     37   5.92    5.75      NA
## rx_c2=1 375     93     NA      NA      NA

生存時間の中央値だけではなく、特定の生存割合の生存時間も出力が行える。以下では、生存割合が75％、50％、25%の順に生存時間を出力した。

quantile(fit, c(1-0.75,1-0.5,1-0.25))

## $quantile
##               25       50 75
## rx_c2=0 3.000000 5.916667 NA
## rx_c2=1 3.416667       NA NA
## 
## $lower
##               25   50 75
## rx_c2=0 2.166667 5.75 NA
## rx_c2=1 2.583333   NA NA
## 
## $upper
##               25 50 75
## rx_c2=0 4.083333 NA NA
## rx_c2=1 4.583333 NA NA

逆に、特定の生存時間のときの生存割合も出力できる。5年時点の生存割合を算出した。
5年生存割合は、Placebo群では62.6%、Estorogen群では66.4%であった。

#Survival probability at a certain time
summary(fit,times = 5)

## Call: survfit(formula = Surv(dtime_y, status5) ~ rx_c2, data = th1)
## 
##                 rx_c2=0 
##         time       n.risk      n.event     survival      std.err lower 95% CI 
##        5.000       28.000       35.000        0.626        0.052        0.532 
## upper 95% CI 
##        0.737 
## 
##                 rx_c2=1 
##         time       n.risk      n.event     survival      std.err lower 95% CI 
##       5.0000      81.0000      89.0000       0.6638       0.0311       0.6055 
## upper 95% CI 
##       0.7278

ログランク検定

ログランク検定(層別、重み付き)

ログランク検定は生存時間に仮定を置く必要がないノンパラメトリックな手法であり、よく生存時間解析では用いられる。単純なログランク検定はsurvdiffで容易に行える。

※KM曲線をみると交差している箇所もあり、ログランク検定が適切であるかは議論の余地があるが、本記事はあくまでRでの実装を主旨としているため、ここでは妥当性については言及しない。

re <- survdiff(Surv(dtime_y, status5) ~ rx_c2, data = th1)
1 - pchisq(re$chisq, length(re$n) - 1) #p-value

## [1] 0.3328387

層別ログランク検定や、重み付きログランク検定は以下のコードで実装できる。層別ログランク検定は交絡調整の方法の1つであり、各層ごとでの結果を統合して最終的な結果を得る手法である。交絡調整としては、後述するCox回帰もしばしば利用される。また重み付きログランク検定はログランク検定が全時点に等しい重みづけを行っているのに対して、例えば試験開始時点の方が人数が多いため開始時に重みを付けたり、薬剤の性質を考慮して適切な重みを付けられる手法である。比例ハザード性の下ではログランク検定の検出力が最大である。

ただし、重み付きの層別解析を実装する方法はわからず、記載できていない。もしご存知の方がいれば、ご連絡頂きたい。

#startified log-rank test
survdiff(Surv(dtime_y, status5) ~ rx_c2 + strata(stage) , data = th1)

#Weighted log-rank test
fit <- survfit( Surv(dtime_y, status5) ~ rx_c2, data = th1 )

##log-rank
ggsurvplot(fit, data = th1, pval = TRUE, pval.method = TRUE)

##Gehan-Breslow
ggsurvplot(fit, data = th1, pval = TRUE, pval.method = TRUE,log.rank.weights = "n")

##Tharone-Ware
ggsurvplot(fit, data = th1, pval = TRUE, pval.method = TRUE,log.rank.weights = "sqrtN")

##Peto-Peto
ggsurvplot(fit, data = th1, pval = TRUE, pval.method = TRUE,log.rank.weights = "S1")
survdiff(Surv(dtime, status5) ~ rx_c2, data = th1, rho = 1)

##Peto-Peto(modified)
ggsurvplot(fit, data = th1, pval = TRUE, pval.method = TRUE,log.rank.weights = "S2")

##Fleming-Harrington
ggsurvplot(fit, data = th1, pval = TRUE, pval.method = TRUE,log.rank.weights = "FH_p=1_q=1")

対比検定

ログランク検定の対比検定を行う。対比検定は3群以上かつ用量反応性のような群間の関連が想定される場合に利用される。多群の検定を単にログランク検定で行うと用量反応性は考慮されず、「全ての群で等しい治療効果」という帰無仮説に対する検定を行うことになる。packageなどで実装する方法が明らかでなかったため、行列計算で対比検定を実装した。

uはログランクスコア、Vがその分散共分散行列になる。c1はPlaceboに-1、Estorogenの各群に1/3をあてており、これはrx_c2のようにPlacebo群とEstrogen群の2群に分類するのと同義の対比係数の割当である。c2は線形な用量反応性あてはめており、Estorogenの投与量が多いほど効果が大きいことが期待される状況での解析となる。

pvalはsurvdiffによって出力されたログランク検定の結果、pval_Lはsurvdiffによって出力されたログランクスコアと分散共分散行列からカイ二乗分布でのp値を出力した結果、pval_c1は対比係数にc1を用いた結果、val_c2は対比係数にc2を用いた結果になる。

結果をみると、pvalとpval_Lはきちんと等しく、pval_c2はそれよりも小さくなった。つまり、用量反応性があることが示唆されるような結果である。また、pval_c1は先ほどPlacebo群とEstorogen群の2群でログランク検定を行ったときのp値と等しくなっていることが確認できた。

#対比検定
re <- survdiff(Surv(dtime_y, status5) ~ rx_c, data = th1)

u <- re[["obs"]]- re[["exp"]]
V <- re[["var"]]
c1 <- c(-1,1/3,1/3,1/3)
c2 <- c(-3/2,-1/2,1/2,3/2)

kai2_L<- t(u) %*% ginv(V) %*% u
kai2_c1 <- (t(c1) %*% u)^2 / (t(c1) %*% V %*% c1)
kai2_c2 <- (t(c2) %*% u)^2 / (t(c2) %*% V %*% c2)

pval <- 1 - pchisq(re$chisq, length(re$n) - 1) 
pval_L = 1-pchisq(kai2_L,3)
pval_c1 = 1-pchisq(kai2_c1,1)
pval_c2 = 1-pchisq(kai2_c2,1)

th3 <- c(pval,pval_L,pval_c1,pval_c2)
names(th3) <- c("pval","pval_L","pval_c1","pval_c2")
print(th3)

##       pval     pval_L    pval_c1    pval_c2 
## 0.03308962 0.03308962 0.33283869 0.04999783

効果量の推定

ログランク検定はノンパラメトリックな仮定の下で行われる検定であり、差の有無は判断できるが効果の大きさは明らかでない。ログランクスコアと分散共分散行列からログランク検定に基づくPeto法で効果量（HR）を推定することもできるが、ここでは比例ハザード性の仮定の下で行われるCox回帰と、比例ハザード性が成り立たない場合の代替指標の1つとしてRMSTの推定を行った。

Cox回帰

Stageで調整を行ったCox回帰を実施した。coefは対数ハザード比になるため、ハザード比はexp(coef)を確認すれば良い。

cox_reg <-coxph( Surv(dtime_y, status5) ~ rx_c2 + stage_c + rx_c2:stage_c, data=th1)
summary(cox_reg)

## Call:
## coxph(formula = Surv(dtime_y, status5) ~ rx_c2 + stage_c + rx_c2:stage_c, 
##     data = th1)
## 
##   n= 502, number of events= 130 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)    
## rx_c2          -0.5567    0.5731   0.3451 -1.613  0.10674    
## stage_c4        1.3010    3.6730   0.3399  3.827  0.00013 ***
## rx_c2:stage_c4  0.3187    1.3753   0.4206  0.758  0.44865    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## rx_c2             0.5731     1.7449    0.2914     1.127
## stage_c4          3.6730     0.2723    1.8866     7.151
## rx_c2:stage_c4    1.3753     0.7271    0.6031     3.137
## 
## Concordance= 0.707  (se = 0.021 )
## Likelihood ratio test= 68.77  on 3 df,   p=8e-15
## Wald test            = 57.52  on 3 df,   p=2e-12
## Score (logrank) test = 70.26  on 3 df,   p=4e-15

RMST

3年で打ち切る条件の下でのRMSTを推定した。Estorogen群での制限付き生存時間は2.51年、Placebo群では2.63年であった。RMSTについて製薬協の資料が参考になる[2]。

rmst <- rmst2(th1$dtime_y, th1$status5, th1$rx_c2, tau=3)
print(rmst)

## 
## The truncation time: tau = 3  was specified. 
## 
## Restricted Mean Survival Time (RMST) by arm 
##               Est.    se lower .95 upper .95
## RMST (arm=1) 2.656 0.041     2.575     2.736
## RMST (arm=0) 2.626 0.073     2.482     2.770
## 
## 
## Restricted Mean Time Lost (RMTL) by arm 
##               Est.    se lower .95 upper .95
## RMTL (arm=1) 0.344 0.041     0.264     0.425
## RMTL (arm=0) 0.374 0.073     0.230     0.518
## 
## 
## Between-group contrast 
##                       Est. lower .95 upper .95     p
## RMST (arm=1)-(arm=0) 0.030    -0.135     0.194 0.724
## RMST (arm=1)/(arm=0) 1.011     0.950     1.076 0.725
## RMTL (arm=1)/(arm=0) 0.921     0.588     1.443 0.719

plot(rmst, xlab="TIME", ylab="Probability", density=60)

サンプルサイズ計算

2群での生存時間解析のサンプルサイズ計算を行った。基本的には、コントロール群のハザード及び、ハザード比からサンプルサイズは計算される。生存時間に指数分布の仮定を置いているため、ハザードは一定となる。

ハザードの推定方法については、

  • n年生存割合からの推定
  • median生存時間からの推定
  • 人年法による推定（ポアソン回帰）

を紹介するが、いづれも指数分布の仮定の下で行われている。指数分布の仮定の下で S(t),f(t),λ(t) の関連から説明できる。

nSurv でサンプルサイズを計算するが、lambdaCがコントロール群のハザード、hrはハザード比、Tは全期間、minfupは最後に登録された試験対象者の追跡期間になる。そのため、T - minfup が登録期間に該当する。

基本のサンプルサイズ計算

#5年生存割合からの推定
pic <- 0.65 
pie <- 0.80
followup <- 5
lambda_c <- -log(pic)/followup
lambda_e <- -log(pie)/followup
hr <- lambda_e / lambda_c

nSurv(lambdaC = lambda_c, hr = hr, T = 5 + .00001 ,minfup = 5,alpha=.05, beta= .2, sided = 2) 

## Fixed design, two-arm trial with time-to-event
## outcome (Lachin and Foulkes, 1986).
## Solving for:  Accrual rate 
## Hazard ratio                  H1/H0=0.518/1
## Study duration:                   T=5
## Accrual duration:                   0
## Min. end-of-study follow-up: minfup=5
## Expected events (total, H1):        73.5703
## Expected sample size (total):       267.5281
## Accrual rates:
##     Stratum 1
## 0-0  26752812
## Control event rates (H1):
##       Stratum 1
## 0-Inf    0.0862
## Censoring rates:
##       Stratum 1
## 0-Inf         0
## Power:                 100*(1-beta)=80%
## Type I error (1-sided):   100*alpha=5%
## Equal randomization:          ratio=1

登録期間を考慮したサンプルサイズ計算

#登録期間2年
nSurv(lambdaC = lambda_c, hr = hr, T = 7,minfup = 5,alpha=.05, beta= .2, sided = 2) 

## Fixed design, two-arm trial with time-to-event
## outcome (Lachin and Foulkes, 1986).
## Solving for:  Accrual rate 
## Hazard ratio                  H1/H0=0.518/1
## Study duration:                   T=7
## Accrual duration:                   2
## Min. end-of-study follow-up: minfup=5
## Expected events (total, H1):        73.3359
## Expected sample size (total):       230.0435
## Accrual rates:
##     Stratum 1
## 0-2  115.0218
## Control event rates (H1):
##       Stratum 1
## 0-Inf    0.0862
## Censoring rates:
##       Stratum 1
## 0-Inf         0
## Power:                 100*(1-beta)=80%
## Type I error (1-sided):   100*alpha=5%
## Equal randomization:          ratio=1

脱落を考慮したサンプルサイズ計算

試験期間中の脱落を考慮したサンプルサイズ設計では、etaでコントロール群の脱落率、etaEで治療群の脱落率を指定することができる。

#脱落を考慮
nSurv(lambdaC = lambda_c, hr = hr, T = 7,minfup = 5,eta = 0.05, etaE = 0.05, 
      alpha=.05, beta= .2, sided = 2) 

## Fixed design, two-arm trial with time-to-event
## outcome (Lachin and Foulkes, 1986).
## Solving for:  Accrual rate 
## Hazard ratio                  H1/H0=0.518/1
## Study duration:                   T=7
## Accrual duration:                   2
## Min. end-of-study follow-up: minfup=5
## Expected events (total, H1):        73.4166
## Expected sample size (total):       263.958
## Accrual rates:
##     Stratum 1
## 0-2   131.979
## Control event rates (H1):
##       Stratum 1
## 0-Inf    0.0862
## Censoring rates:
##       Stratum 1
## 0-Inf      0.05
## Power:                 100*(1-beta)=80%
## Type I error (1-sided):   100*alpha=5%
## Equal randomization:          ratio=1

（補足）HRの推定

#HRの推定 いづれも指数分布の仮定を置く
##median生存時間
median_time <-5
lambda <- log(2)/median_time
print(lambda)

## [1] 0.1386294

##人年法（ポアソン回帰）
th1 %>% 
  filter(dtime != 0) %>% 
  mutate(log_dtime_y = log(dtime_y)) -> th4

re <- glm(status5 ~ rx_c , data=th4, family="poisson", offset = log_dtime_y)
v <- re[["coefficients"]]
lambda_c <- exp(v[1])
lambda_e <- exp(v[1]+v[2])
hr <- lambda_e / lambda_c
hr_2 <- exp(v[2])

th5 <- c(lambda_c,lambda_e,hr,hr_2)
names(th5) <- c("lambda_c","lambda_e","hr","hr_2")            
print(th5)

##   lambda_c   lambda_e         hr       hr_2 
## 0.10664357 0.08997709 0.84371786 0.84371786

まとめ

カプランマイヤー曲線やMedian生存時間などの生存時間の記述から、ログランク検定、Cox回帰、RMST、サンプルサイズ計算までの生存時間解析での基本的な解析を紹介した。本記事ではPlacebo群とEstorogen群の2群に分類したが、例えばエストロゲン0.2 mgでは治療効果がないと考えPlacebo群と0.2 mgをあわせた群を対照に、エストロゲン治療群（1.0 mg, 5.0 mg）の効果を推定するような解析も考えられる。ちなみにこの分類であればKM曲線はきれいに差がつき交わりは見られなくなる。またRに不慣れな点も多く誤りがあるかもしれないが、もし誤りを見つけられた方は、ご連絡頂ければ幸いである。

（メモ）今後、追加する可能性のある解析

  • パラメトリックな解析（生存時間に仮定を置くモデル）
  • RMSTの回帰分析
  • 中間解析を考慮したサンプルサイズ設計
  • 競合リスクを考慮した解析（累積発生関数：Gray検定）

サポートして頂けるとモチベーションに繋がりますのでぜひ宜しくお願いします。データ解析や臨床研究でのご相談があれば、お気軽にTwiiterもしくはメールにてご連絡下さい。

作成者：Masahiro Kondo
作成日：2022/2/25
連絡先：m.kondo1042(at)gmail.com

脚注

  1. https://www.rdocumentation.org/packages/survival/versions/2.11-4/topics/survfit ↩︎

  2. 生存時間型応答の評価指標（第2版） ↩︎

Discussion

Michel TakahashiMichel Takahashi

非常に勉強になりました。
・比例ハザード性の検定(二重対数プロット、累積シェーンフェルド残差の図示)
・累積発症関数の図示
・Fine & Gray検定のやり方
・交互作用の検討のやり方
などご教授いただけましたら幸いです。

KONDO MasahiroKONDO Masahiro

ご質問ありがとうございます。いずれも重要な点なのですが、私の時間が確保できず記事に記載できておりません。またどれもプログラム例等の提示には時間がかかり、すぐにはお答えすることが難しいです。大変申し訳ありません。