カプランマイヤー曲線をggplot系の図で作りたい【R】【ggsurvplot】【 survminerパッケージ】

Rでカプランマイヤー曲線の図を作る機会があった。　

せっかくなので、デフォルトのsurvivalパッケージだけではなく， survminerパッケージのggsurvplotというggplot系コードを使って鮮やかにしてみたので、備忘録的に

今回使うデータとは異なるデータですが、こんな感じのが作れたりします

使うデータと分析

データはMASSパッケージのgehanを使います。

このページで示されてるのと同じ分析をggsurvplotでやってみます

https://www1.doshisha.ac.jp/~mjin/R/Chap_36/36.html

通常のsurvival パッケージでのカプランマイヤー

#パッケージの読み込み

library(survival)　#生存時間の関数を作るために呼びだし

library(survminer)　# ggsurvplotのパッケージ

library(MASS)　# gehanデータを使うために呼び出し

data(gehan) #データの呼び出し 

#変数リスト

#pair:  Pair index,  time: Remission time (weeks),  cens:  0=censored,  treat: Treatment: control or 6-MP

ge.sf <- survival::survfit(Surv(time,cens)~treat,data=gehan) # モデル
plot(ge.sf,lty=1:2)　#プロット

legend("topright" ,c("6-MP投与群","対照群"), lty=c(1,2))

ggsurvplotでカプランマイヤーやってみる

シンプルなカプランマイヤー曲線

ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,,xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray())

# sizeは曲線の太さ， ggthemeでは背景を様々に変えることができます。

ちょっとオプションを足すと・・・

#信頼区間やLog-rank testの結果を示す　pval や conf.int をTRUEに

ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,
pval =T ,conf.int = T, xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray())

図の設定をいじりたい

#ggplot系ですので、どんどん後ろにスクリプトを足していくことで、いろいろな設定を変えていけます。

#たとえば　さきほどのシンプルな図を　temp　に代入し，

temp <- ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,,xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray())

temp$plot + scale_x_continuous(breaks=seq(0,40,4))

#と付け足すと

x軸を4weekごとの目盛りを入れることができます。

ちょっとした注意点として，tempに代入後，temp$plot の続きに足していく要素を書かないとエラーが出てきてしまいます。

Risk Table?

あと社会学系では見たことないのですが、リスクテーブルなるものも表示できるようです。

ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,,xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray(),risk.table=T)

おまけ

最初の図はこんなコードでつくりました

fit <- survfit(Surv(month,event)~age311_cat,data=dat)
temp <- ggsurvplot(fit, data = dat, size=2,legend.title = "世帯主年齢", legend.labs=c("20未満","20代","30代","40代","50代","60代","70代","80歳以上"),legend=c(.2,.85),conf.int = F, ylab="すまい再建率", xlab="時間", censor=F,pval = F, ggtheme = theme_gray(), fun="event")
temp <- temp$plot + theme(axis.text=element_text(size=30))
temp <- temp + theme(axis.title.x = element_text(size=40),axis.title.y = element_text(size=40))
temp <- temp + scale_x_continuous(breaks=c(0,12,24,36,48,60,72), labels=c("2011/3", "2012/3","2013/3","2014/3","2015/3","2016/3","2017/3"))
temp <- temp+ylim(0,1)
temp + theme(legend.title = element_text(size=38),legend.text = element_text(size=38))

 # fun= "event" で　y軸を逆転できます。　あとはこのデータは，2011年3月観察開始で，仮設住宅を退居するまでの時間を分析する図なので，x軸は12ヶ月後ごとに設定，そのラベルとして年を示しています。