「統計学関連なんでもあり」の過去ログ--- 046

No.20748　時間依存性ROC曲線のAUCの95%信頼区間と差の検定の仕方について　　【BOZ】　2014/01/20(Mon) 18:10

Cox回帰分析による2つの予後予測モデル(モデル1とモデル2）のうち，どちらが正確に予後を推測出来るかを調べたいと思っています。
http://d.hatena.ne.jp/kaz_yos/20130113/1358094223のサイトを参考にして，RのrisksetROCを使い，2つのモデルの5年後生存を予測する時間依存性ROC曲線の曲線下面積（AUC）を出したのですが，AUCの95%信頼区間と差の検定の仕方がわからず困っています。
調べた範囲ではブートストラップ法で求められるようなのですが，Rのどのパッケージでどのようにコマンドを入力すればよいのかわからないでおります。
勉強不足で申し訳ございませんが，御教授いただけたら幸いです。

以下，risksetROCで入力したコマンドです。

> library(risksetROC)

> #モデル1のAUCを算出。
> ROC.model1=risksetROC(Stime=data$time, status=data$status,
+ marker=data$cox.model1, predict.time=365*5, method="Cox",
+ main="ROC Curve", lty=2, col="red")
> ROC.mode1$AUC
[1] 0.7947956

> #モデル2のAUCを算出。
> ROC.model2=risksetROC(Stime=data$time, status=data$status,
+ marker=data$cox.model2, predict.time=365*5, method="Cox",
+ main="ROC Curve", lty=2, col="red")
> ROC.model2$AUC
[1] 0.5778788

No.20800　Re: 時間依存性ROC曲線のAUCの95%信頼区間と差の検定の仕方について　　【taipapa】　2014/01/27(Mon) 21:09

遅くなりましたが，どなたも反応されないようなので．．．
pROCでできますよ．
http://rcommanderdeigakutoukeikaiseki.com/ROC_AUC_NRI_and_IDI.html
このページが分かりやすいです．
但し，このページはbootstrapに関して触れていませんが，method="bootstrap"のオプションをつければ出来ます．

No.20806　Re: 時間依存性ROC曲線のAUCの95%信頼区間と差の検定の仕方について　　【BOZ】　2014/01/29(Wed) 11:51

taipapa様

コメント大変ありがとうございます。
具体的なコマンドとしては，以下のようにするということでしょうか？

例えばモデル1のAUCの信頼区間を求めたい場合だと，
library(pROC)
ROC1 <- roc(response=data$status, predictor=data$cox.model1)
ci(ROC1$auc, method="bootstrap")

ただこれだと，risksetROCで入力したpredict.time=365*5（5年後の生存を予測するために入力するコマンドです）という情報入ってないので，正しくないのではと思ったのですが，いかがでしょうか？検討違いな事を言っていたら申し訳ございません。

No.20807　Re: 時間依存性ROC曲線のAUCの95%信頼区間と差の検定の仕方について　　【BOZ】　2014/01/29(Wed) 14:12

あれから懸命に調べ，risksetROCパッケージとbootパッケージの二つを用い下記のようにコマンドを入力したところ，モデル1のAUCの95%信頼区間は0.71～0.86となりました。明らかな間違い等ありましたら教えていただきたいところですが・・。
周りに聞ける人がいないので正直困っております(^_^;)

library(risksetROC)

#以下，ブートストラップの対象となる数式(for.boot)を作成
for.boot<-function(x,indices)
{
　Data<-x[indices,]　#無作為復元抽出するためのコマンド
　
#Cox回帰モデルを作成。Data$timeは生存時間，Data$statusは生死の情報。　
　fit0 <- coxph(Surv(Data$time, Data$status)~ Data$factor1 + Data$factor2)
　eta <- fit0$linear.predictor　
　
#以下，http://cran.r-project.org/web/packages/risksetROC/risksetROC.pdfの　　risksetROCにデータの当てはめ。生存予測を5年後(365*5)に設定するため，predict.timeを365*5に設定。
　ROC.model1=risksetROC(Stime=Data$time, status=Data$status,
　marker=eta, predict.time=365*5, method="Cox")　
　
　return(ROC.model1$AUC)　#5年後生存率予測能を示すAUCが算出
}

上記の数式(for.boot)を使って，5000回のブートストラップを施行。
library(boot)
result<-boot(data=Data, statistic=for.boot, R=5000)
quantile(result$t, c(.025,.975))
2.5% 97.5%
0.7172453 0.8643236

No.20808　Re: 時間依存性ROC曲線のAUCの95%信頼区間と差の検定の仕方について　　【taipapa】　2014/01/30(Thu) 00:29

Cox回帰モデルのROCだったんですね．よく読んでませんでした…(^^;;;
Cox回帰モデルなら，rms packageのcphとHmisc packageのrcorr.censが使えそうです（bootstrapではないですが）．
risksetROCって，知りませんでしたが，これのマニュアルの例を使うと，

> fit0.cph <- cph( Surv(survival.time,survival.status) ~ score + cell.type + tx + age, na.action=na.omit, x=TRUE, y=TRUE, surv=TRUE )
> rcorr.cens(x=fit0.cph$linear.predictors, S=Surv(survival.time,survival.status))
       C Index           Dxy           S.D.              n        missing …  
  2.621050e-01  -4.757900e-01   4.175206e-02   1.370000e+02   0.000000e+00   
 
上記のDxyは，Dxy = 2(c-0.5)
このcがAUCにあたります．ですから，
> -0.4758/2 + 0.5
[1] 0.2621

で，上記のC Indexと合います．
S.D.はDxyの標準誤差ですので，
CstatisticCI <- function(x) # x is object of rcorr.cens.
{
se <- x["S.D."]/2
Low95 <- x["C Index"] - 1.96*se
Upper95 <- x["C Index"] + 1.96*se

cbind(x["C Index"], Low95, Upper95)
}
とでも定義しておけば，
> CstatisticCI(rcorr.cens(x=fit0.cph$linear.predictors, S=Surv(survival.time,survival.status)))
 
                     Low95    Upper95
C Index 0.262105   0.221188   0.303022
となります．0.5未満なので，ひっくり返して，AUC = 0.74 といったほうが良いかもしれませんね．

とりあえず，これで95%信頼区間は出せます．ただ，BOZさんの場合は，時間依存性が肝で，5年後の値がほしいんですよね．これだと全体の評価ですから違いますね．risksetROCのソースを見ると，Coxの場合は，predict.timeは，最初のcoxphの部分では使われず，関数 CxWeightsに渡されて，その値以下のcoxphのlinear.predictorsのみを選択するために使われてますね．
違うような気がするけど，全部5年後で打ち切りにして上記の関数を使ってみるとか．．．

bootstrapに関しては，違うもので確かめたいですよね．rms packageのvalidateが参考になるかも．．．
> validate(fit0.cph, dxy=TRUE, B=100)
          index.orig     training    test      optimism   index.corrected   n
Dxy       0.4758        0.4931       0.4574     0.0357      0.4401         100
R2        0.3574        0.3929       0.3318     0.0611      0.2963         100
Slope     1.0000        1.0000       0.8682     0.1318      0.8682         100
D          0.0593        0.0679      0.0540     0.0139      0.0453         100
U         -0.0020       -0.0020      0.0047   -0.0067       0.0047         100
Q         0.0612         0.0699      0.0493     0.0207     0.0406          100
g          0.9616       1.0664       0.9065    0.1599      0.8017          100

No.20809　Re: 時間依存性ROC曲線のAUCの95%信頼区間と差の検定の仕方について　　【BOZ】　2014/01/30(Thu) 17:56

丁寧にコマンドまで記載していただきまして，大変ありがとうございます！rms packageやHmisc packageのことは知らなかったので，やってみたいと思います。
まずは取り急ぎお礼申し上げます。