Go to Qiita Advent Calendar 2024 Top
LoginSignup
4
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@Ringa_hyj(Ringa hyj)

Rでの生存時間解析

Last updated at Posted at 2022-02-23

だいぶ前にRで生存時間分析するとどんな感じなのか学びながら触ってみたメモが残っていたので備忘録として記事化しました。
内容はあまり期待しないでください。ご指摘大歓迎です。

#生存時間解析とは

イベントが起こるまでの時間や起こる確率を分析する手法。
薬効の確認、術後の生存日数、部品の故障までの時間、などをデータから分析する。
イベントが起こったかどうかはバイナリで持つと考える。

一般的に死亡をイベントと考えるが、起きて欲しいイベントが観測されるまでの期待的な時間という考え方で使う場合もある。
図にすると、観察開始からガントチャートを書いて死亡した場合に印をつける。
スタートをそろえてロリポップplotを行う場合もある。

生存時間解析の特徴的な考えに「打ち切り(censoring)」がある。
観察が途中で中止になった個体、イベントが発生せず試験期間が終わってしまった場合、が打ち切りとして扱われる。
途中で中止した場合などがある(右側打ち切り)。
主に扱う必要があるのはこの右側打ち切りになるので後述。
左側打ち切りは、スタート時期が分からないが観察が始まったデータのこと。
左側打ち切りは薬効を確かめたくて被験者を集めたが、既にその薬を服用していた対象者がいる場合などが当てはまる。

#打ち切りの種類

基本的に打ち切りの扱いは無情報打ち切りが想定される

  • 無情報打ち切り
    • 試験途中に追跡不可能(消息不明,来院しない)になった場合
    • 試験終了まで問題なく進んでいた場合(タイプ1打ち切り：時間打ち切り)
    • 観測したい個数のイベントが確認できたため試験を終了した(タイプ2打ち切り：個数打ち切り)
    • 観察対象によって試験期間がランダムに決まり、短くても試験終了となる(タイプ3打ち切り：ランダム打ち切り)

ただし、上記のように無情報打ち切りに見えても理由を調べると情報を持っている場合もある

  • 無情報打ち切りでない場合
    • 病状が悪化したためにリタイア
    • 悪化して病院を移った
    • 逆に病状が改善したために試験への参加をやめた
    • 上記のような「生存確率」に影響を与える情報を持ったリタイア

打ち切りは観察を中止することから観察打ち切りとも呼ばれる
本来観察したいイベントに対して、それ以外の別イベントが発生する場合もある
病気の生存を確認するはずが事故で死亡したり、自殺や他殺が発生するなど
これらを「競合リスク(competing risk)」と呼び、これを考慮しないで同じイベントとして分析すると誤った結果になる
競合リスクを考えた手法までは紹介できなかったが、
木構造モデルを使った手法での生存時間解析(ルール・アンサンブル法)が手法として知られているため、
「ggrandomforest」パッケージの紹介をしておく。

#生存者数の減少を見るグラフ

生存者が徐々に減っていく過程を図にしたいと考える。

library(tidyverse)
library(survival)
survival::leukemia

data <- leukemia %>% mutate(id = row_number())

ggplot(data, aes(x=id, y=time)) +
  geom_segment( aes(x=id, xend=id, y=0, yend=time), color="skyblue") +
  geom_point(aes(x=id,y=time,shape=as.factor(status)), color="blue", size=4, alpha=0.6) +
  scale_shape("event")+
  labs(title = "leukemia data")+
  theme_light()+
  coord_flip()

image.png

・Tはイベントが発生する瞬間
・生存関数とはS(t)と表す
・tがTよりも小さい範囲においての生存確率を出力する関数を生存関数S(t)と考える
・S(t) = P(T > t)

・打ち切りとか厄介な問題を一旦置いておく
・スタートした時点の総数からイベントにより死亡していくことが想像できる
・スタートが100人(100%)として考えると、以下のような図がかける

leukemia %>% 
  group_by(time) %>% 
  summarise(duplicate=n()) %>% 
  arrange(desc(time)) %>% 
  mutate(id = row_number(),prob=1/n(),cum=cumsum(prob)) %>%
  arrange(time)%>%
  ggplot() +
  geom_step(aes(x=time, y=cum))

image.png

これが生存関数の根本的な考え方。

50%の人がまだ生きている時刻を調べたければ,y=0.5の交点を調べれば時間が分かる

leukemia %>% 
  group_by(time) %>% 
  summarise(duplicate=n()) %>% 
  arrange(desc(time)) %>% 
  mutate(id = row_number(),prob=1/n(),cum=cumsum(prob)) %>%
  arrange(time)%>%
  ggplot() +
  geom_step(aes(x=time, y=cum))+
  geom_hline(yintercept = 0.5)

image.png

この(ニセ)生存関数について、ある時刻時点で重複している場合を考えたり、重複しているが片方は無情報ウチキリで片方は死亡の場合どうするか、
などのパターンを考慮した生存関数がカプラン・マイヤー法。

#カプラン・マイヤー法の計算手順と描画

まず、「イベントor打ち切り」と「時間」に対して重複をまとめる

df <- leukemia %>%
  arrange(time) %>% 
  mutate(num = 1,ranking=row_number()) %>% 
  group_by(time,status) %>% 
  summarise(events = sum(num),non_sum_rank=min(ranking)) %>% 
  ungroup() %>%
  pivot_wider(names_from = status, values_from = events, names_glue = "dupli_{status}",values_fill=0) %>% 
  group_by(time) %>% 
  summarise(sum_dupli_1=sum(dupli_1),sum_dupli_0=sum(dupli_0),non_sum_rank_pivot=min(non_sum_rank)) %>%
  ungroup() %>%
  mutate(after_sum_rank=row_number())

df

image.png

重複(順位タイ)を考慮した順位と、考慮せず集計後にランキングをつけた列を作った.

ここで、一番早いtime=5の時について考える。
5の時点でイベント、つまり死亡が2件起きている。
5になる時間的に微小時間だけ前には全数23人が生存していた、と考えると、
「死亡するかもしれないリスクのある集団」は23人いた。この人数を「リスク集合」と呼ぶ。

リスク集合は、この23人からイベントもしくはウチキリの数を順次引いていくことで求まる。
(リタイアしていくと考えると、残った集団がリスクを持っている)

num<-23

for(i in 1:(nrow(df)-1) ){
  num<-c(num,num[i]-(df$sum_dupli_1[i]+df$sum_dupli_0[i]))
}

df2 <- df %>% mutate(risk = num)
df2

image.png

ハザード関数とは、HAZ

ある任意の不定の時点tまで生きていたとき、次の瞬間死ぬ(イベントが起こる)確率を出力するもの。

P(T<t+Δ | T>t)

イベントはtの微小区間だけ先に発生すると考える

上記で作成したデータから各時点の死亡数を、各時点のリスク集団の大きさで割る。これがHAZ関数の値。
(その時点でイベントが起こる確率を集団数から割合として表現していると解釈)

逆に生存する確率も各時点で考える事ができる。
生存する確率、は
「イベントが起きない(打ち切り) = 1 or 1-イベントが起きる確率」
で数値化できる。

この生存率を累積積していく。
確率は始まりの1から順次減っていくので、リスク集合が減っていくという条件のもとでの生存率として累積積を計算している
生存率の列はより正確には「ある時点での条件付き生存率」と考えることができる。

df3<- df2 %>% mutate(haz = sum_dupli_1/risk, suv=1-haz, suv_prod=cumprod(suv))
df3

image.png

このようにして累積生存確率(生存関数)をつくる方法がカプラン・マイヤー法と呼ばれている。
ハザードが起こった時に降下する。
一般的に打ち切り地点に縦棒で印をつけ、リスク集合が減少していることを表現する。

df3 %>%
  ggplot()+
  aes(x=time,y=suv_prod)+
  geom_step()

image.png

もちろんこんな単純な話だけでカプランマイヤー法が成り立っているわけではないので、
これが全てと思わないでもらいたい。

しっかりパッケージになっているのでそちらを使う。

  • カプランマイヤーはデータから生存関数を求める手法であり。ノンパラメトリックな方法。
  • より滑らかに表現するために(負の)指数関数を仮定して生存関数を近似するのがexp生存関数モデル、パラメトリック。exp(-1 * haz * t)
  • この中間くらいに位置しているのがcoxのモデル。セミパラメトリック。
    • Cox proportional hazard model:(コックス比例ハザードモデル)
    • 指数関数ではなめらかすぎて急にイベントが大量発生してから、おちつくようなシグモイドカーブ様を再現できない
    • coxph()ならば表現することもできる

#コックス比例ハザードモデルcoxph

fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian) 
plot(survfit(fit, newdata=data.frame(age=60)),
     xscale=365.25, xlab = "Years", ylab="Survival")

image.png

王道として「survival」が生存時間分析のパッケージであり、
これの分析や結果表示・作図をサポートするものに「survminer」がある。

survivalでもplot()を使うことで作図はできるが、ggsurvplot()の方が機能的に便利。
実は「ggRandomForests」にもgg_survival()がこっそり潜んでいるのでこれも最後に紹介しておく。

もしかしたらsurvivalパッケージと分けて信頼区間を表示したい人がいるかもしれない。
その場合には「km.ci」というパッケージがsurvivalの解析結果を支援しているので紹介。

まず、単一の生存関数をカプランマイヤーで描いてもらいたいので、「Surv object」をfomulaのyとして渡し
影響を与えるグループは考えないとしてチルダ右側に1を入力

fit <- survfit(Surv(leukemia$time,leukemia$status)~1,type="kaplan-meier")
plot(fit, xlab = "time", ylab = "survival probability",main="S(t) - kaplan-meier model",conf.int= FALSE,mark.time=TRUE)

image.png

50%生存日数は27日と分かる
その信頼区間も確認できる

fit

image.png

先ほど計算したテーブルが返ってくる

summary(fit)

image.png

ここで、生存関数が求まったので、その時点で生存している確率が求まった。
しかし、今は確率と時点の点でしか求められない。
時点での考えられる生存確率の範囲として、時点ごとの分散を求めたい。

以下のように95%程度に信頼できる範囲を示したいとする。

fit <- survfit(Surv(leukemia$time,leukemia$status)~1)
plot(fit, xlab = "time", ylab = "survival probability",conf.int=0.95)

image.png

#信頼区間・分散の推定

時刻ごとの生存確率がどれだけ不確実なものかを確認することが出来た。

しかし、手元のデータは少なく、各時点ごとに約一点ずつしか得られないのになぜばらつきが求まるのか？
これは理論的に近似する分散から考える場合と、ブートストラップ法によって求める場合の二つが広く使われている。

理論的な近似方法として「グリーンウッド(Greenwood)の近似式」が有名

Var(\hat{S}(t)) = \hat{S}(t)^2 \sum_{i=t} \frac{d_i}{N_i(N_i - d_i)}

この平方根を標準誤差として扱い、正規分布の範囲で信頼区間を構築する。
ただし、確率の上限である1を超えず、0よりも下回らないように描画する。

外側のパッケージとして「km.ci」の紹介
実行は各自で

library(survival)
library(km.ci)
library(tidymodels)

#直腸がん:打ち切りや生存者はほとんど含まれていない
data(rectum.dat)

fit <- survfit(Surv(time, status) ~ 1, data=rectum.dat)
plot(fit)
fit%>%tidy

fit2 <- km.ci(fit)
fit2%>%tidy
plot(fit2)

km.ci(fit, conf.level=0.95, tl=NA, tu=NA, method="rothman")

image.png

#二種類の生存関数をつくる

leukemiaデータには化学療法を行ったかどうかを示すx列が存在する。
化学療法を使った場合と使わなかった場合で生存関数はどう変わるのか？
survfitのformula右式にグループ分けの条件を書き足す。

fit <- survfit(Surv(leukemia$time,leukemia$status)~leukemia$x)
p<- plot(fit, xlab = "time", ylab = "survival probability",conf.int=FALSE,col=c("skyblue","pink"))

text(p$x[1],p$y[1],labels="Maintained")
text(p$x[2],p$y[2],labels="Nonmaintained")

image.png

#二つの生存関数は統計的に違うものなのか検定

処置群と非処置群で生存時間が異なると言えるのであれば、延命を望む人には処置をしたほうがいい。
この根拠を検定して確かめる方法を紹介する。
二群で生存関数に差があるのかを確かめる場合には「Log-rank検定」が最も単純。

R言語では「survdiff()」によって検定を行うことができる。(Harrington and Fleming (1982)の実装版)
rho = 0の場合：log-rank test または Mantel-Haenszel 検定であり，
rho = 1の場合：Gehan-Wilcoxon 検定(Peto & Peto 補正)
分けたい群の中にNAが入っている場合、例えばコントロール群であった場合などには事前に埋めておく。

survdiff(Surv(leukemia$time,leukemia$status)~leukemia$x,rho = 0)

image.png

p値は0.07で95%信頼区間には届かないため異なる生存関数とは言えなさそう。

複数の変数を考慮してモデルを作る場合、変数を+で繋ぐ

mod <- survfit(Surv(time, status) ~ pat.karno + inst, data=lung)
plot(mod)

image.png

明示的にstrata()を指定することで、「層別log-rank検定」も実行できる

survdiff(Surv(time, status) ~ pat.karno + strata(inst), data=lung)

image.png

#plotの補助「surveminer」

library(survminer)

fit<- survfit(Surv(time, status) ~ sex, data = lung)

ggsurvplot(fit, data = lung)

image.png

ggsurvplot(fit, 
           pval = TRUE, 
           surv.median.line="hv",
           ggtheme = theme_light(),
           risk.table=T
           )

image.png

累積イベント率

ggsurvplot(fit, pval = TRUE, 
           xlab="time", 
           surv.median.line="hv",
           ggtheme = theme_light(),
           risk.table=T,
           fun="event")

image.png

#plotの補助「ggrandomforest」

生存時間解析では木構造法で生存関数を当てはめる場合もある。
ggrandomforestのplotを使う場合を紹介しておく。

library(survival)
library(cowplot)
library(randomForestSRC)
library(ggRandomForests)

data(pbc, package="randomForestSRC")
pbc$time <- pbc$days/364.25

gg_dta <- kaplan(interval="time", censor="status", 
                      data=pbc)
#plot(gg_dta, error="none")

gg_dta_s <- gg_survival(interval="time", censor="status", 
                      data=pbc,type="kaplan")
#plot(gg_dta_s)

plot_grid(plot(gg_dta), plot(gg_dta_s))

typeでの指定もデフォルトはカプランマイヤー

image.png

Nelson-Aalenの推定量を使いたい場合はtypeで
type = c("kaplan", "nelson")
調節するか
nelson(interval, censor, data, by = NULL, weight = NULL, ...)
を使用する。

二群に分けたければbyで指定

gg_dta <- gg_survival(interval="time", censor="status", 
                      data=pbc, by="treatment",, conf.int=0.5)

plot(gg_dta)

image.png

#coxphの変数選択

モデル内の「Likelihood Ration Test」を元に、変数を追加した場合の挙動を比較して、
変数を増やしたほうがいいかどうかを考える

cox_1 <- coxph(Surv(time,status)~treatment,data=pbc[pbc$status<2,])
cox_2 <- coxph(Surv(time,status)~treatment+hepatom,data=pbc[pbc$status<2,])
anova(cox_1,cox_2,test="LRT)

変数を増やしてもp値が大きいなら追加することの有意性を感じない。

#coxphからハザード関数を出力

coxph_model <- coxph(Surv(time,status)~trt,data=pbc[pbc$status<2,])
data<-basehaz(coxph_model, centered=TRUE)
plot(data$time,data$hazard,type="l")

image.png

#生存時間分析に使えるデータセット

data detail pkg
data(heart) stanford大学,心臓移植 survival
stanford2 stanford大学,心臓移植 survival
bladder 膀胱癌,bladder1と2もある survival
data("gbsg") 乳がん survival
cgd 慢性顆粒球症 survival
cgd0 慢性肉芽腫疾患 survival
data("retinopathy") 糖尿病網膜症 survival
data("diabetic") 糖尿病網膜症 survival
leukemia,aml 急性骨髄性白血病 survival
colon ステージB/Cの大腸がん survival
mgus,mgus1 モノクローナル・ガンモパチー(MGUS),mgus2は重大性不明のもの survival
lung,cancer NCCTG肺がんデータ survival
data("transplant") 肝移植 survival
veteran 肺がん survival
data("udca") 原発性胆汁性肝硬変(PBC),ウルソデオキシコール酸(UDCA) survival
nafld ノンアルコール脂肪肝疾患 survival
myeloid 急性骨髄性白血病 survival
pbc 胆汁性肝硬変 survival
pbcseq 胆汁性肝硬変 survival
kidney 腎臓,カテーテル挿入時の感染症 survival
data(flchain) 血清遊離軽鎖(FLC) survival
rats,rats2 ラット,発がん性物質 survival
ovarian 卵巣癌 survival
data("rotterdam") 乳がん survival
nwtco 腫瘍研究 survival
data(uspop2) 米国の年齢 2000年～2020年 survival

#参考

生存時間解析　ハザード関数からCox比例ハザードモデルまで　by using Python
生存時間分析の色々なアルゴリズムをまとめてみました
Survival Analysis | Concepts and Implementation in R
カプラン・マイヤー法: 生存時間解析の基本手法 (統計学One Point)
R doc survival
R doc ggRandomForests
R doc survminer

生存時間研究におけるルールアンサンブル法の開発
競合リスクを伴う生存時間データに対するルール・アンサンブル法の開発
生存時間解析入門 汪　金芳
医学統計勉強会 東北大学病院循環器内科
R による生存時間解析
Cox の比例ハザードモデルについて 土居

その他紹介

R Lecture lesson5 生存時間分析
R: Kaplan-Meier曲線を作成する
カプランマイヤー法（Kaplan-Meier method)によるイベント発生率の点推定とNumber at riskの算出
Kaplan-Meierの生存曲線をRで理解する
第11章　生命表解析
臨床泌尿器，40(12)，985-990，1986
競合リスクデータの解析
Cox比例ハザードモデルの変数選択
PUBG - Survival Analysis (Kaplan-Meier)
ログ・ランク検定
Rと生存時間分析（１）

4
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?