3

我正在尝试使用ggplot2和在以下位置找到的代码绘制一些 Kaplan-Meier 曲线:https ://github.com/kmiddleton/rexamples/blob/master/qplot_survival.R

我在不同的数据库中使用这个很棒的代码得到了很好的结果。但是,在这种情况下,它给了我以下错误......好像我的数据框中有空行:

Error en if (nrow(layer_data) == 0) return() : argument is of length zero.

以前有关此类错误的问题似乎对我没有用,因为在我的情况下,数据和函数的类型并不相同。

我对使用 R 进行统计分析比较陌生,而且我没有编程背景,所以我认为这一定是我的数据中的一个“愚蠢的错误”,但我找不到它在哪里……看来 ggplot2 绝对可以'找不到要绘制的行。拜托,你能以任何方式帮助我,提供线索、建议等吗?

这是我的数据和使用的代码,按顺序为控制台准备好-我在 knitr 脚本中尝试过-。最后,我发布了我的 sessionInfo:

library(splines)
library(survival)
library(abind)
library(ggplot2)
library(grid)

我创建了一个名为acbi30(真实数据)的数据框:

mort28day <- c(1,0,1,0,0,0,0,1,0,0,0,1,1,0,1,0,0,1,0,1,1,1,1,0,1,1,1,0,0,1)
daysurv <- c(4,29,24,29,29,29,29,19,29,29,29,3,9,29,15,29,29,11,29,5,13,20,22,29,16,21,9,29,29,15)
levo <- c(0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0)
acbi30 <- data.frame(mort28day, daysurv, levo)
save(acbi30, file="acbi30.rda")
acbi30

然后,我粘贴以下命令以使用 ggplot2 创建一个函数:

t.Surv <- Surv(acbi30$daysurv, acbi30$mort28day)
t.survfit <- survfit(t.Surv~1, data=acbi30)


#define custom function to create a survival data.frame#
createSurvivalFrame <- function(f.survfit){

#initialise frame variable#
f.frame <- NULL

#check if more then one strata#
if(length(names(f.survfit$strata)) == 0){

#create data.frame with data from survfit#
f.frame <- data.frame(time=f.survfit$time, n.risk=f.survfit$n.risk, n.event=f.survfit$n.event, n.censor = f.survfit
$n.censor, surv=f.survfit$surv, upper=f.survfit$upper, lower=f.survfit$lower)
#create first two rows (start at 1)#
f.start <- data.frame(time=c(0, f.frame$time[1]), n.risk=c(f.survfit$n, f.survfit$n), n.event=c(0,0),
n.censor=c(0,0), surv=c(1,1), upper=c(1,1), lower=c(1,1))
#add first row to dataset#
f.frame <- rbind(f.start, f.frame)
#remove temporary data#
rm(f.start)
}
else {
#create vector for strata identification#
f.strata <- NULL
for(f.i in 1:length(f.survfit$strata)){
#add vector for one strata according to number of rows of strata#
f.strata <- c(f.strata, rep(names(f.survfit$strata)[f.i], f.survfit$strata[f.i]))
}
#create data.frame with data from survfit (create column for strata)#
f.frame <- data.frame(time=f.survfit$time, n.risk=f.survfit$n.risk, n.event=f.survfit$n.event, n.censor = f.survfit
$n.censor, surv=f.survfit$surv, upper=f.survfit$upper, lower=f.survfit$lower, strata=factor(f.strata))
#remove temporary data#
rm(f.strata)
#create first two rows (start at 1) for each strata#
for(f.i in 1:length(f.survfit$strata)){
#take only subset for this strata from data#
f.subset <- subset(f.frame, strata==names(f.survfit$strata)[f.i])
#create first two rows (time: 0, time of first event)#
f.start <- data.frame(time=c(0, f.subset$time[1]), n.risk=rep(f.survfit[f.i]$n, 2), n.event=c(0,0),
n.censor=c(0,0), surv=c(1,1), upper=c(1,1), lower=c(1,1), strata=rep(names(f.survfit$strata)[f.i],
2))
#add first two rows to dataset#
f.frame <- rbind(f.start, f.frame)
#remove temporary data#
rm(f.start, f.subset)
}
#reorder data#
f.frame <- f.frame[order(f.frame$strata, f.frame$time), ]
#rename row.names#
rownames(f.frame) <- NULL
}
#return frame#
return(f.frame)
}


#define custom function to draw kaplan-meier curve with ggplot#
qplot_survival <- function(f.frame, f.CI="default", f.shape=3){
#use different plotting commands dependig whether or not strata's are given#
if("strata" %in% names(f.frame) == FALSE){
#confidence intervals are drawn if not specified otherwise#
if(f.CI=="default" | f.CI==TRUE ){
#create plot with 4 layers (first 3 layers only events, last layer only censored)#
#hint: censoring data for multiple censoring events at timepoint are overplotted#
#(unlike in plot.survfit in survival package)#
ggplot(data=f.frame) + geom_step(aes(x=time, y=surv), direction="hv") + geom_step(aes(x=time,
y=upper), directions="hv", linetype=2) + geom_step(aes(x=time,y=lower), direction="hv", linetype=2) +
geom_point(data=subset(f.frame, n.censor==1), aes(x=time, y=surv), shape=f.shape)
}
else {
#create plot without confidence intervals#
ggplot(data=f.frame) + geom_step(aes(x=time, y=surv), direction="hv") +
geom_point(data=subset(f.frame, n.censor==1), aes(x=time, y=surv), shape=f.shape)
}
}
else {
if(f.CI=="default" | f.CI==FALSE){
#without CI#
ggplot(data=f.frame, aes(group=strata, colour=strata)) + geom_step(aes(x=time, y=surv),
direction="hv") + geom_point(data=subset(f.frame, n.censor==1), aes(x=time, y=surv), shape=f.shape)
}
else {
#with CI (hint: use alpha for CI)#
ggplot(data=f.frame, aes(colour=strata, group=strata)) + geom_step(aes(x=time, y=surv),
direction="hv") + geom_step(aes(x=time, y=upper), directions="hv", linetype=2, alpha=0.5) +
geom_step(aes(x=time,y=lower), direction="hv", linetype=2, alpha=0.5) +
geom_point(data=subset(f.frame, n.censor==1), aes(x=time, y=surv), shape=f.shape)
}
}
}

绘制全局生存曲线(95% CI):

它没有给出任何错误:

# Kaplan-Meier plot, global survival (with CI)
t.survfit <- survfit(t.Surv~1, data=acbi30)
t.survframe <- createSurvivalFrame(t.survfit)
t.survfit
qplot_survival(t.survframe, TRUE, 20)

绘制分层生存曲线:

给出上述错误:

# Kaplan-Meier plot, stratified survival
t.survfit2 <- survfit(t.Surv~levo, data=acbi30)
t.survframe2 <- createSurvivalFrame(t.survfit2)
t.survfit2
qplot_survival(t.survframe2, TRUE, 20)

在没有 ggplot2 的情况下绘制结果:

t.survframe2 的结构对我来说似乎还可以,没有任何空行,所以它一定是 qplot_survival 读取我在t.survframe2中的数据的问题。制作一个简单的情节不会返回任何错误:

t.survframe2
plot(t.survfit2)

我的数据框的问题在哪里?创建的函数适用于其他数据集,但不适用于这个......

先感谢您,

马列维夫

会话信息:

sessionInfo()

R 版本 2.15.2 (2012-10-26) 平台:i386-w64-mingw32/i386(32 位)

locale:

[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
[3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
[5] LC_TIME=Spanish_Spain.1252    

attached base packages:
[1] grid      splines   stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
[1] ggplot2_0.9.3    abind_1.4-0      survival_2.36-14 knitr_0.8       

loaded via a namespace (and not attached):
 [1] colorspace_1.1-1   dichromat_1.2-4    digest_0.5.2      
 [4] evaluate_0.4.2     formatR_0.7        gtable_0.1.2      
 [7] labeling_0.1       MASS_7.3-22        munsell_0.4       
[10] plyr_1.8           proto_0.3-9.2      RColorBrewer_1.0-5
[13] reshape2_1.2.1     scales_0.2.3       stringr_0.6.1     
[16] tools_2.15.2    
4

2 回答 2

2

qplot_survival()我对你的功能做了一点整容手术。主要问题似乎是您在data =论点中的子集条件geom_point;在t.survframet.survframe2中,生成n.censor值 0、3 和 12 的表。通过将子集条件更改为n.censor > 0,我设法在所有情况下都得到了一个图。我也没有看到f.CI = "default",所以我将默认设置为 TRUE 并相应地修改了 if 条件。

qplot_survival <- function(f.frame, f.CI= TRUE, f.shape=3)
{
 # use different plotting commands depending whether 
 # or not strata are given#
if(!("strata" %in% names(f.frame)))
{
  #confidence intervals are drawn if not specified otherwise#
   if( isTRUE(f.CI) )
   {
      # create plot with 4 layers (first 3 layers only events, 
      # last layer only censored)#
      # hint: censoring data for multiple censoring events at 
      # timepoint are overplotted#
      # (unlike in plot.survfit in survival package)#
   ggplot(data=f.frame) + 
      geom_step(aes(x=time, y=surv), direction="hv") + 
      geom_step(aes(x=time, y=upper), direction ="hv", linetype=2) + 
      geom_step(aes(x=time,y=lower), direction="hv", linetype=2) +
      geom_point(data=subset(f.frame, n.censor > 0), 
                 aes(x=time, y=surv), shape=f.shape)
   } else {
  #create plot without confidence intervals#
   ggplot(data=f.frame) + 
      geom_step(aes(x=time, y=surv), direction="hv") +
      geom_point(data=subset(f.frame, n.censor > 0), 
                 aes(x=time, y=surv), shape=f.shape)
          }
} else {
  if( !(isTRUE(f.CI)) ){
#without CI#
  ggplot(data=f.frame, aes(group=strata, colour=strata)) + 
     geom_step(aes(x=time, y=surv), direction="hv") + 
     geom_point(data=subset(f.frame, n.censor > 0), 
                aes(x=time, y=surv), shape=f.shape)
} else {

#with CI (hint: use alpha for CI)#
  ggplot(data=f.frame, aes(x = time, colour=strata, group=strata)) + 
      geom_step(aes(y=surv), direction="hv") + 
      geom_step(aes(y=upper), direction="hv", 
                   linetype=2, alpha=0.5) +
      geom_step(aes(y=lower), direction="hv", 
                   linetype=2, alpha=0.5) +
      geom_point(data=subset(f.frame, n.censor > 0), 
                 aes(y=surv), shape=f.shape)
      }
   }
}

进行这些更改后,以下图表都对我有用:

qplot_survival(t.survframe2, TRUE, 20)
qplot_survival(t.survframe2, FALSE, 20)
qplot_survival(t.survframe, TRUE, 20)
qplot_survival(t.survframe, FALSE, 20)

一些评论:

  1. 函数内部的子集化可能很危险,因为有时,如本例所示,满足条件会返回零行数据帧。我会考虑该geom_point()层是否真的有必要。
  2. 在几个地方,你有directions = "hv"一个geom_step()电话。这个论点不是复数形式,上面已经改变了。
  3. 我认为这可以更有效地完成,但是从survfit对象中提取感兴趣的列的一种方法t.survfit是这样的:

(当存在地层时扩展组合)

comps <- c(2:6, 8, 10);
t.fit <- as.data.frame(do.call(cbind, lapply(comps, function(j) t.survfit[[j]])))
names(t.fit) <- names(t.survfit)[comps]
于 2013-01-05T06:13:03.587 回答
1

这是另一个版本,它也解释了数据中没有审查点的情况(@Dennis 的版​​本在这种情况下仍然失败)。这可以提高效率,可能是通过创建一个变量来预先存储整个数据帧中有多少审查点,然后重新使用它,而不是像我在每种情况下再次进行测试。

# define custom function to draw kaplan-meier curve with ggplot
qplot_survival <- function(f.frame, f.CI="default", f.shape=3){
  # use different plotting commands dependig whether or not strata's are given
  if("strata" %in% colnames(f.frame) == FALSE){
    # confidence intervals are drawn if not specified otherwise
    if(f.CI=="default" | f.CI==TRUE ){
      # create plot with 4 layers (first 3 layers only events, last layer only censored)
      # hint: censoring data for multiple censoring events at timepoint are overplotted



      # (unlike in plot.survfit in survival package)
      p <- ggplot(data=f.frame) + geom_step(aes(x=time, y=surv), direction="hv") + geom_step(aes(x=time, 
                                                                                            y=upper), directions="hv", linetype=2) + geom_step(aes(x=time,y=lower), direction="hv", linetype=2)
      if(nrow(subset(f.frame, n.censor > 0)) > 0){
        p+geom_point(data=subset(f.frame, n.censor > 0), aes(x=time, y=surv), shape=f.shape)
      }else{
        p
      }
    }
    else {
      # create plot without confidence intervalls
      p <- ggplot(data=f.frame) + geom_step(aes(x=time, y=surv), direction="hv")
      if(nrow(subset(f.frame, n.censor > 0)) > 0){
        p + geom_point(data=subset(f.frame, n.censor > 0), aes(x=time, y=surv), shape=f.shape)
      }else{
        p
      }
    }
  }
  else {
    if(f.CI=="default" | f.CI==FALSE){
      # without CI 
      p <- ggplot(data=f.frame, aes(group=strata, colour=strata)) + geom_step(aes(x=time, y=surv), 
                                                                         direction="hv") 
      if(nrow(subset(f.frame, n.censor > 0)) > 0){
        p +geom_point(data=subset(f.frame, n.censor > 0), aes(x=time, y=surv), shape=f.shape)
      }else{
        p
      }
    }
    else {
      # with CI (hint: use alpha for CI)
      p <- ggplot(data=f.frame, aes(colour=strata, group=strata)) + geom_step(aes(x=time, y=surv), 
                                                                         direction="hv") + geom_step(aes(x=time, y=upper), directions="hv", linetype=2, alpha=0.5) + 
        geom_step(aes(x=time,y=lower), direction="hv", linetype=2, alpha=0.5)
      if(nrow(subset(f.frame, n.censor > 0)) > 0){
        p + geom_point(data=subset(f.frame, n.censor > 0), aes(x=time, y=surv), shape=f.shape)
      }else{
        p
      }
    }
  }
}
于 2013-11-17T20:33:53.983 回答