5

生存分析を行う際にggplotまたはlatticeを利用する方法を知っている人はいますか？トレリスまたはファセットのようなサバイバルグラフを作成すると便利です。

結局、私は遊んで、カプラン・マイヤープロットの解決策を見つけました。リスト要素をデータフレームに取り込む際の厄介なコードについてお詫びしますが、別の方法を理解できませんでした。

注：2つのレベルの階層でのみ機能します。誰かが私がx<-length(stratum)これを行うために使用できる方法を知っているなら、私に知らせてください（Stataではマクロに追加できます-これがRでどのように機能するかわからない）。

ggkm<-function(time,event,stratum) {

    m2s<-Surv(time,as.numeric(event))

    fit <- survfit(m2s ~ stratum)

    f$time <- fit$time

    f$surv <- fit$surv

    f$strata <- c(rep(names(fit$strata[1]),fit$strata[1]),
            rep(names(fit$strata[2]),fit$strata[2])) 

    f$upper <- fit$upper

    f$lower <- fit$lower

    r <- ggplot (f, aes(x=time, y=surv, fill=strata, group=strata))
        +geom_line()+geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.3)

    return(r)
}
4

2 に答える 2

4

私はで次のコードを使用していますlattice。最初の関数は1つのグループのKM曲線を描画し、通常は関数として使用されpanel.groupますが、2番目の関数はパネル全体のログランク検定のp値を追加します。

 km.panel <- function(x,y,type,mark.time=T,...){
     na.part <- is.na(x)|is.na(y)
     x <- x[!na.part]
     y <- y[!na.part]
     if (length(x)==0) return()
     fit <- survfit(Surv(x,y)~1)
     if (mark.time){
       cens <- which(fit$time %in% x[y==0])
       panel.xyplot(fit$time[cens], fit$surv[cens], type="p",...)
      }
     panel.xyplot(c(0,fit$time), c(1,fit$surv),type="s",...)
}

logrank.panel <- function(x,y,subscripts,groups,...){
    lr <-  survdiff(Surv(x,y)~groups[subscripts])
    otmp <- lr$obs
    etmp <- lr$exp
    df <- (sum(1 * (etmp > 0))) - 1
    p <- 1 - pchisq(lr$chisq, df)
    p.text <- paste("p=", signif(p, 2))
    grid.text(p.text, 0.95, 0.05, just=c("right","bottom"))
    panel.superpose(x=x,y=y,subscripts=subscripts,groups=groups,...)
}

このコードが機能するには、打ち切りインジケーターが0-1である必要があります。使用法は次のようになります。

library(survival)
library(lattice)
library(grid)
data(colon)  #built-in example data set
xyplot(status~time, data=colon, groups=rx, panel.groups=km.panel, panel=logrank.panel)

'panel = panel.superpose'を使用するだけでは、p値は取得されません。

于 2010-06-11T13:49:26.017 に答える
1

私はあなたがあなたの更新された答えで使用するアプローチをほぼ正確に従うことから始めました。しかし、survfitについてイライラするのは、各ティックではなく、変更のみをマークすることです。たとえば、0-100％、1-100％、2-100ではなく0-100％、3-88％になります。 ％、3-88％。これをggplotにフィードすると、線はフラットのままで3で真っ直ぐ下がるのではなく、0から3に傾斜します。これは、アプリケーションと仮定によっては問題ない場合がありますが、従来のKMプロットではありません。これは私がさまざまな数の層を処理した方法です：

groupvec <- c()
for(i in seq_along(x$strata)){
    groupvec <- append(groupvec, rep(x = names(x$strata[i]), times = x$strata[i]))
}
f$strata <- groupvec

価値があるのは、これが私がやった方法です-しかし、これは実際にはKMプロットでもありません。これは、KM推定値自体を計算していないためです（ただし、打ち切りはありません。したがって、これは同等です。 .. 私は信じている）。

survcurv <- function(surv.time, group = NA) {
    #Must be able to coerce surv.time and group to vectors
    if(!is.vector(as.vector(surv.time)) | !is.vector(as.vector(group))) {stop("surv.time and group must be coercible to vectors.")}

    #Make sure that the surv.time is numeric
    if(!is.numeric(surv.time)) {stop("Survival times must be numeric.")}

    #Group can be just about anything, but must be the same length as surv.time
    if(length(surv.time) != length(group)) {stop("The vectors passed to the surv.time and group arguments must be of equal length.")}

    #What is the maximum number of ticks recorded?
    max.time <- max(surv.time)  

    #What is the number of groups in the data?
    n.groups <- length(unique(group))

    #Use the number of ticks (plus one for t = 0) times the number of groups to
    #create an empty skeleton of the results.
    curves <- data.frame(tick = rep(0:max.time, n.groups), group = NA, surv.prop = NA)

    #Add the group names - R will reuse the vector so that equal numbers of rows
    #are labeled with each group.
    curves$group <- unique(group)

    #For each row, calculate the number of survivors in group[i] at tick[i]
    for(i in seq_len(nrow(curves))){
      curves$surv.prop[i] <- sum(surv.time[group %in% curves$group[i]] > curves$tick[i]) /
          length(surv.time[group %in% curves$group[i]])
    }

  #Return the results, ordered by group and tick - easier for humans to read.
  return(curves[order(curves$group, curves$tick), ])   

}
于 2010-07-07T16:24:07.003 に答える