survival
データ
ここでは、データを作成します。
censor = 1 とは、アウトカムのイベント発生を意味します。
dfSurvival <- data.frame(
patientid = c(1,2,3,4,5,6,7,8,9,10,11,12,13),
year = c(10,9,9,8,7,6,5,4,4,4,3,2,1),
censor = c(1,1,1,1,0,1,0,1,1,0,1,0,1),
group = c("exposure","control","exposure","control","exposure","control","exposure","control","exposure","control","exposure","control","exposure")
)
dfSurvival## patientid year censor group
## 1 1 10 1 exposure
## 2 2 9 1 control
## 3 3 9 1 exposure
## 4 4 8 1 control
## 5 5 7 0 exposure
## 6 6 6 1 control
## 7 7 5 0 exposure
## 8 8 4 1 control
## 9 9 4 1 exposure
## 10 10 4 0 control
## 11 11 3 1 exposure
## 12 12 2 0 control
## 13 13 1 1 exposure変数
- dfSurvival: データフレームであることを明示するために、dfを先頭につけた変数名としています。
関数
- data.frame: base 関数。データフレームを作成する。他のプログラミング言語だと . は意味を持つことが多いが、R ではあまり意味はない。
実際には、File > Import Dataset などを使って、エクセルファイルまたはCSVファイルから読み込みます。 この時に、変数名を指定します。 変数名は、dataset や Dataset がデフォルトになっていることが多いようです。
生データからピボット（標準）
研究の生データは、通常は上の表のようなデータではありません。
実際には、毎年生存の確認をする場合、以下のような表を作成するでしょう。
## Warning in get_current_theme(): This functionality requires shiny v1.6 or higherこのようなデータ形式を、 long フォーマットといいます。 このデータを、表１のようなフォーマット（wide フォーマット）に成型します。
reshape2 を用いる場合
略
tidyr を用いる場合
library(tidyr)
dfTemp <- tidyr::pivot_wider(data = dfSurvivalEntry, names_from = "year", values_from = "censor")
#datatable(dfTemp)
dfTemp## # A tibble: 13 x 12
## patientid `2000` `2001` `2002` `2003` `2004` `2005` `2006` `2007` `2008`
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 1 0 0 0 0 0 0 0 0 0
## 2 2 0 0 0 0 0 0 0 0 0
## 3 3 0 0 0 0 0 0 0 0 0
## 4 4 NA NA 0 0 0 0 0 0 0
## 5 5 0 0 0 0 0 0 0 0 NA
## 6 6 NA 0 0 0 0 0 0 1 NA
## 7 7 NA 0 0 0 0 0 0 NA NA
## 8 8 NA 0 0 0 0 1 NA NA NA
## 9 9 NA 0 0 0 0 1 NA NA NA
## 10 10 NA 0 0 0 0 1 NA NA NA
## 11 11 0 0 0 1 NA NA NA NA NA
## 12 12 NA 0 0 0 NA NA NA NA NA
## 13 13 NA NA 0 1 NA NA NA NA NA
## # … with 2 more variables: 2009 <int>, 2010 <int>変数
- dfSurvivalWider: データフレーム（tibble）形式
関数
- tidyr::pivot_wider: ライブラリ tidyr が提供する関数
これは、確かに見やすいですが、表１の形式とは異なります。
そこで、
dfSurvivalWider <- tidyr::pivot_wider(data = dfSurvivalEntry, id_cols = "patientid", names_from = "censor", values_from = "year", values_fn = min)
names(dfSurvivalWider)[c(2,3)] <- c("alive","dead") # 列名を変えます。
dfSurvivalWider## # A tibble: 13 x 3
## patientid alive dead
## <int> <int> <int>
## 1 1 2000 2010
## 2 2 2000 2009
## 3 3 2000 2009
## 4 4 2002 2010
## 5 5 2000 NA
## 6 6 2001 2007
## 7 7 2001 NA
## 8 8 2001 2005
## 9 9 2001 2005
## 10 10 2001 2005
## 11 11 2000 2003
## 12 12 2001 NA
## 13 13 2002 2003関数
c: base 関数。ベクトルを返す。combination の c。
names: base 関数。列名を表示するのが主な機能。ここでは列名変換に使用した。
このままですと、イベント発生しなかった人が NA のままになってしまいます。
dfSurvivalWider2 <- tidyr::pivot_wider(data = dfSurvivalEntry, id_cols = "patientid", names_from = "censor", values_from = "year", values_fn = max)
names(dfSurvivalWider2)[2] <- "alivelast" # 列名を変えます。
dfSurvivalWider2## # A tibble: 13 x 3
## patientid alivelast `1`
## <int> <int> <int>
## 1 1 2009 2010
## 2 2 2008 2009
## 3 3 2008 2009
## 4 4 2009 2010
## 5 5 2007 NA
## 6 6 2006 2007
## 7 7 2006 NA
## 8 8 2004 2005
## 9 9 2004 2005
## 10 10 2004 2005
## 11 11 2002 2003
## 12 12 2003 NA
## 13 13 2002 2003library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondfSurvivalWider <- inner_join(dfSurvivalWider, dfSurvivalWider2, by = "patientid")
dfSurvivalWider## # A tibble: 13 x 5
## patientid alive dead alivelast `1`
## <int> <int> <int> <int> <int>
## 1 1 2000 2010 2009 2010
## 2 2 2000 2009 2008 2009
## 3 3 2000 2009 2008 2009
## 4 4 2002 2010 2009 2010
## 5 5 2000 NA 2007 NA
## 6 6 2001 2007 2006 2007
## 7 7 2001 NA 2006 NA
## 8 8 2001 2005 2004 2005
## 9 9 2001 2005 2004 2005
## 10 10 2001 2005 2004 2005
## 11 11 2000 2003 2002 2003
## 12 12 2001 NA 2003 NA
## 13 13 2002 2003 2002 2003次に、alive 列 と dead 列 から、duration 列 を計算します。
dfSurvivalWider$Duration <- dfSurvivalWider$alivelast - dfSurvivalWider$alive
dfSurvivalWider$Duration <- dfSurvivalWider$dead - dfSurvivalWider$alive
dfSurvivalWider[dfSurvivalWider$patientid == 1, "Group"] <- "exposure"
dfSurvivalWider[dfSurvivalWider$patientid == 2, "Group"] <- "control"
dfSurvivalWider$Group <- as.factor(dfSurvivalWider$Group)
dfSurvivalWider## # A tibble: 13 x 7
## patientid alive dead alivelast `1` Duration Group
## <int> <int> <int> <int> <int> <int> <fct>
## 1 1 2000 2010 2009 2010 10 exposure
## 2 2 2000 2009 2008 2009 9 control
## 3 3 2000 2009 2008 2009 9 <NA>
## 4 4 2002 2010 2009 2010 8 <NA>
## 5 5 2000 NA 2007 NA NA <NA>
## 6 6 2001 2007 2006 2007 6 <NA>
## 7 7 2001 NA 2006 NA NA <NA>
## 8 8 2001 2005 2004 2005 4 <NA>
## 9 9 2001 2005 2004 2005 4 <NA>
## 10 10 2001 2005 2004 2005 4 <NA>
## 11 11 2000 2003 2002 2003 3 <NA>
## 12 12 2001 NA 2003 NA NA <NA>
## 13 13 2002 2003 2002 2003 1 <NA>dplyr を駆使する方法もありますが、これより簡単には思えません。
最後の Group 列に exposure または control を入力する方法は、美しい方法ではありません。 しかし、エクセルなどの票に格納して inner_join などをするよりは、わかりやすく、また情報漏洩やミスのリスクも減るような気がします。
検定
library(survival)
objSurv <- survival::survfit( survival::Surv(year,censor) ~ group, data = dfSurvival)
summary(objSurv)## Call: survfit(formula = survival::Surv(year, censor) ~ group, data = dfSurvival)
##
## group=control
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 5 1 0.800 0.179 0.5161 1
## 6 3 1 0.533 0.248 0.2142 1
## 8 2 1 0.267 0.226 0.0507 1
## 9 1 1 0.000 NaN NA NA
##
## group=exposure
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 7 1 0.857 0.132 0.633 1
## 3 6 1 0.714 0.171 0.447 1
## 4 5 1 0.571 0.187 0.301 1
## 9 2 1 0.286 0.223 0.062 1
## 10 1 1 0.000 NaN NA NA変数
objSurv: survival パッケージのオブジェクト。
dfSurvival: 上で作ったデータフレーム。
関数
survfit: survival パッケージの関数。生存時間データのカプランマイヤー推定で用いる関数
Surv: survival パッケージの関数。生存時間データオブジェクトを作成する。
summary: base 関数。
２群の比較 log-rank 検定
ログ・ランク (log-rank)検定は、群ごとのイベントありとなしに集計したクロス表のカイ２乗検定。R の survival パッケージが自動的に計算する。
survdiff(Surv(year,censor)~group, data = dfSurvival) ## Call:
## survdiff(formula = Surv(year, censor) ~ group, data = dfSurvival)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## group=control 6 4 3.58 0.0486 0.105
## group=exposure 7 5 5.42 0.0321 0.105
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7プロット
下準備
plot
plot(objSurv)
いろいろなオプションを試してみましょう。
plot(objSurv, conf.int = T, col = c("red","blue"), xlab= "year", ylab = "event")
ggplot2
survminer
library(survminer, warn.conflicts=FALSE, quietly=TRUE)
ggsurvplot(objSurv, data = dfSurvival, risk.table = TRUE, pval = TRUE, conf.int = TRUE, risk.table.y.text.col = TRUE)
ggfortify
library(ggplot2)
library(ggfortify)
autoplot(objSurv, data = dfSurvival )
参考
- 東北大学 医学統計勉強会 : 東北大学病院 - 循環器内科 第７回 「生存時間解析 生存曲線 Cox比例ハザードモデル」 https://www.cardio.med.tohoku.ac.jp/2005/newmember/medical_statistics.html