Rでの生存時間解析


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


#生存時間解析とは

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

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

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



#打ち切りの種類

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

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

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

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


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



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

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


```{r}
library(tidyverse)
library(survival)
survival::leukemia
```


```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5e7193e0-d69a-0de1-c638-b8d318a91d82.png)


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


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



```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/2ad6b132-ce4f-70a6-c8bc-5e45636a0faa.png)


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

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


```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/8ed5f19d-bdd0-07e6-9894-806a0a898844.png)


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


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


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



```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e73dabcf-93d6-65ff-6913-46c6baad94f1.png)



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

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

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


```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/9f3b42e6-9f6a-7a60-e6d1-8a3a23a76830.png)



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


```{math}
P(T<t+Δ | T>t)
```


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



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

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

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



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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/aaa8a50b-89b3-742e-6211-039cabf4dcde.png)


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



```{r}
df3 %>%
  ggplot()+
  aes(x=time,y=suv_prod)+
  geom_step()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/6c090541-bfd3-966b-18ea-879f94a51a35.png)


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

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


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

```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/68416d05-2a2b-4a2a-c32d-cac43442b431.png)



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

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


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



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


```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/062b7d47-d95c-2617-003a-5d34ffe7fe06.png)




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



```{r}
fit
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/15d849bf-baf9-87ea-fd42-ace83493c8f6.png)



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


```{r}
summary(fit)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/7d0b2a26-9bb2-c402-31b5-52e4c108287c.png)


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

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


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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/7d5ddb99-c5d8-30ea-99b4-f52e0170af26.png)





#信頼区間・分散の推定

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

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

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



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

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




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

```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/48bfc530-596f-f440-ff76-1d45b575a436.png)




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


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


```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/858a83fc-a63e-932c-d5db-6a7255c14e55.png)



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


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


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



```{r}
survdiff(Surv(leukemia$time,leukemia$status)~leukemia$x,rho = 0)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5ebf9d2e-2134-0e74-24b8-5cd76682b4a1.png)



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

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

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


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4989057b-895d-795a-cafa-5637c1ba4476.png)


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


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



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d40bc105-52e1-535e-f1f7-ee880bd04a8b.png)



#plotの補助「surveminer」



```{r}
library(survminer)

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

ggsurvplot(fit, data = lung)
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/807560ec-91ef-592e-4b31-09b5357d01c3.png)





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





![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5059b6c9-7a1a-6147-5b70-764e409b0b2a.png)



累積イベント率


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



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/99024f30-cd0e-816b-2ffd-ee190cfcae89.png)



#plotの補助「ggrandomforest」





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




```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4bb8536e-b551-89d4-ce07-97efee97e646.png)





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


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



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

plot(gg_dta)

```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/a130aa6f-652a-e1d6-17d9-503b803a7c15.png)








#coxphの変数選択

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


```{r}
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からハザード関数を出力


```{r}
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b568c4b7-a29e-019e-691e-49228fd0889e.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](https://qiita.com/MorinibuTakeshi/items/f637e22d5cf216d9974b)
[生存時間分析の色々なアルゴリズムをまとめてみました](https://qiita.com/saltcooky/items/409329485be499a5b270)
[Survival Analysis | Concepts and Implementation in R](https://www.youtube.com/playlist?list=PLqzoL9-eJTNDdnKvep_YHIwk2AMqHhuJ0)
[カプラン・マイヤー法: 生存時間解析の基本手法 (統計学One Point)](https://www.amazon.co.jp/%E3%82%AB%E3%83%97%E3%83%A9%E3%83%B3%E3%83%BB%E3%83%9E%E3%82%A4%E3%83%A4%E3%83%BC%E6%B3%95-%E7%94%9F%E5%AD%98%E6%99%82%E9%96%93%E8%A7%A3%E6%9E%90%E3%81%AE%E5%9F%BA%E6%9C%AC%E6%89%8B%E6%B3%95-%E7%B5%B1%E8%A8%88%E5%AD%A6One-%E8%A5%BF%E5%B7%9D-%E6%AD%A3%E5%AD%90/dp/4320112628)
[R doc survival](https://www.rdocumentation.org/packages/survival/versions/3.2-7)
[R doc ggRandomForests](https://www.rdocumentation.org/packages/ggRandomForests/versions/2.0.1)
[R doc survminer](https://www.rdocumentation.org/packages/survminer/versions/0.4.8)

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


その他紹介

[R Lecture lesson5 生存時間分析](https://rpubs.com/Ryohei/567609)
[R: Kaplan-Meier曲線を作成する](http://www.yamamotor.jp/?p=2697)
[カプランマイヤー法（Kaplan-Meier method)によるイベント発生率の点推定とNumber at riskの算出](http://rcommanderdeigakutoukeikaiseki.com/Kaplan_Meier.html)
[Kaplan-Meierの生存曲線をRで理解する](https://mmbiostats.com/kaplan_meier_r)
[第11章　生命表解析](http://www.snap-tck.com/room04/c01/stat/stat11/stat1101.html)
[臨床泌尿器，40(12)，985-990，1986](http://halbau.world.coocan.jp/life.html)
[競合リスクデータの解析](https://www.icrweb.jp/course/view.php?id=341)
[Cox比例ハザードモデルの変数選択](https://www.nag-j.co.jp/naglib/tr/SurvivalAnalysis/VariableSelection.htm)
[PUBG - Survival Analysis (Kaplan-Meier)](https://www.kaggle.com/datark1/pubg-survival-analysis-kaplan-meier)
[ログ・ランク検定](http://jspt.japanpt.or.jp/ebpt_glossary/log-rank-test.html)
[Rと生存時間分析（１）](https://www1.doshisha.ac.jp/~mjin/R/Chap_36/36.html)