r - 使用 survplot 绘制有风险的主题时出现错误消息

Question

尝试在 survplot 中沿 x 轴绘制有风险的受试者时，我收到以下错误消息：

Error in text.default(tt[-1], yy, nri[-1], cex = cex.n.risk, adj = 1) : zero-length 'labels' specified

有什么帮助吗？我对生存分析相当陌生，未能找到对此错误的任何解释。总体而言，代码看起来不错，除非我为绘图添加 n.risk=TRUE 选项，否则会出现错误。任何提示将不胜感激。非常感谢。

下面是数据以及使用的代码。

这里的数据

Duration <-  structure(list(conflict = c("Angola 75-89", "Angola 89-91", "Angola 92-94", 
    "Azerb (N-K) 89-94", "Bosnia 92-95", "Cambodia 70-91", "Chad 79-79", 
    "Chad 89-96", "Chechnya 94-96", "Colombia 48-57", "Croatia 91-91 (?)", 
    "Croatia 95-95", "DomRep 65-65", "El Salv 79-91", "GeorgA 89-92", 
    "GeorgB 92-94", "Guatem 68-96", "India 46-48", "Iraq 61-70", 
    "Laos 59-73", "Lebanon 58-58", "Lebanon 75-89", "Liberia 89-93", 
    "Malaysia 48-56", "Moldova 92-92", "Mozamb 81-92", "Nicara 81-89", 
    "Phil 72-96", "Rwanda 90-93", "SieLeo 91-96", "Stafrica 83-91", 
    "Sudan 63-72", "Tajik 92-97", "Yemen 62-70", "Zimbab 72-79", 
    "Guinea-Bissau June - November 1998", "Liberia 94-96", "Papua New Guinea 1990 - 2001", 
    "Afghanistan 1978 - 2001", "Ethiopia 1961-1993", "Indonesia (Aceh) 1976 - 2005", 
    "Kenya 2007- 2008", "Nepal 1996 - 2006", "Somalia 1991 - 2008", 
    "Bangladesh 1997", "Burundi 1993-2005", "Cote d'Ivoire 2002-2007", 
    "Democratic Republic of Congo 98-03", "Northern Ireland (68-98)", 
    "Darfur, Sudan 2003-2010", "Sudan 83-05", "Liberia 1999-2003"
    ), peacedur = c(2, 17, 58, 175, 157, 206, 7, 117, 34, 322, 43, 
    157, 520, 204, 192, 171, 144, 0.100000001490116, 48, 25, 199, 
    230, 12, 626, 89, 195, 232, 148, 8, 6, 204, 141, 138, 357, 348, 
    122, 40, 23, 0.100000001490116, 60, 40, 8, 24, 0.100000001490116, 
    133, 28, 22, 69, 128.5, 3, 71, 83), peacefail = c(1, 1, 1, 0, 
    0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 
    0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 
    0, 1, 0, 1, 1, 0), totalps = c(0L, 2L, 3L, 2L, 3L, 2L, 1L, 2L, 
    4L, 2L, 1L, 1L, 1L, 3L, 2L, 2L, 2L, 2L, 3L, 1L, 2L, 2L, 1L, 2L, 
    3L, 3L, 4L, 4L, 3L, 3L, 4L, 3L, 2L, 4L, 3L, 2L, 1L, 1L, 2L, 2L, 
    4L, 1L, 3L, 1L, 3L, 3L, 3L, 2L, 2L, 4L, 3L, 3L), year_end = c(1989L, 
    1991L, 1994L, 1994L, 1995L, 1991L, 1979L, 1996L, 1996L, 1957L, 
    1991L, 1995L, 1965L, 1991L, 1992L, 1994L, 1996L, 1948L, 1970L, 
    1973L, 1958L, 1989L, 1993L, 1956L, 1992L, 1992L, 1989L, 1996L, 
    1993L, 1996L, 1991L, 1972L, 1997L, 1970L, 1979L, 1998L, 1996L, 
    2001L, 2001L, 1993L, 2005L, 2008L, 2006L, 2008L, 1997L, 2005L, 
    2007L, 2003L, 1998L, 2010L, 2005L, 2003L), peacedur.year = c(1, 
    2, 5, 15, 14, 18, 1, 10, 3, 27, 4, 14, 44, 17, 16, 15, 12, 1, 
    4, 3, 17, 20, 1, 53, 8, 17, 20, 13, 1, 1, 17, 12, 12, 30, 29, 
    11, 4, 2, 1, 5, 4, 1, 2, 1, 12, 3, 2, 6, 11, 1, 6, 7), SurvObj = structure(c(2, 
    17, 58, 175, 157, 206, 7, 117, 34, 322, 43, 157, 520, 204, 192, 
    171, 144, 0.100000001490116, 48, 25, 199, 230, 12, 626, 89, 195, 
    232, 148, 8, 6, 204, 141, 138, 357, 348, 122, 40, 23, 0.100000001490116, 
    60, 40, 8, 24, 0.100000001490116, 133, 28, 22, 69, 128.5, 3, 
    71, 83, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 
    1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 
    1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0), .Dim = c(52L, 2L), .Dimnames = list(
        NULL, c("time", "status")), type = "right", class = "Surv")), .Names = c("conflict", 
    "peacedur", "peacefail", "totalps", "year_end", "peacedur.year", 
    "SurvObj"), row.names = c("1", "2", "3", "4", "5", "6", "7", 
    "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", 
    "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", 
    "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", 
    "41", "42", "43", "44", "45", "46", "47", "48", "49", "51", "52", 
        "53"), class = "data.frame")

生存对象的创建

library(survival)
library(rms)
Duration$SurvObj <- with(Duration, Surv(peacedur, peacefail==1))

适合 + 转换为 npsurv

KM.Duration.totalps <- survfit(SurvObj ~ totalps, data = Duration, conf.type = "log-log")
class(KM.Duration.totalps) <- c(class(KM.Duration.totalps), "npsurv")

剧情：

survplot(KM.Duration.totalps, 
         xlab="duration in months", ylab="survival prob",
         conf="none",
         label.curves = TRUE,                    
         time.inc=12,
         levels.only  = FALSE,                   
         n.risk=TRUE)

Answer 1

totalps=0只有 1 个事件。不n.risk喜欢那样。

运行survplot(KM.Duration.totalps[-1], ...)与设置一样n.risk = FALSE

library(survival)
library(rms)
Duration$SurvObj <- with(Duration, Surv(peacedur, peacefail==1))

KM.Duration.totalps <- survfit(SurvObj ~ totalps, data = Duration, conf.type = "log-log")
class(KM.Duration.totalps) <- c(class(KM.Duration.totalps), "npsurv")

summary(KM.Duration.totalps)

# Call: survfit(formula = SurvObj ~ totalps, data = Duration, conf.type = "log-log")
# 
# totalps=0 
# time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
# 2            1            1            0          NaN           NA           NA 
# 
# totalps=1 
# time n.risk n.event survival std.err lower 95% CI upper 95% CI
# 0.1     10       1    0.900  0.0949       0.4730        0.985
# 7.0      9       1    0.800  0.1265       0.4087        0.946
# ...



par(mfrow = c(2,1))
survplot(KM.Duration.totalps,
         xlab="duration in months", ylab="survival prob",
         conf="none",
         label.curves = TRUE,                    
         time.inc=12,
         levels.only  = FALSE,                   
         n.risk=FALSE)

survplot(KM.Duration.totalps[-1],
         xlab="duration in months", ylab="survival prob",
         conf="none",
         label.curves = TRUE,                    
         time.inc=12,
         levels.only  = FALSE,                   
         n.risk=TRUE)

深入挖掘问题的根源，这条线rms:::survplot.npsurv

text(tt[-1], yy, nri[-1], cex = cex.n.risk, adj = 1)

您的 totalps=0 的 nri 是长度为 1 的向量，因此 R 正在尝试绘制

text(1, 1, integer(0))

试试看，你会得到同样的错误。所以要解决这个问题，要么使用上面的解决方案（绘制 totalps=0 反正不是很有趣，因为它只是一条直线），或者你可以像下面这样编辑源代码，只在结尾附近插入一个 if 语句代码可以解决问题if (length(nri) > 1)

所以现在你可以使用你的新函数来获得完整的表格/情节没有错误（不过，我不会，因为你可以看到，标签掩盖了你的风险表）

survplot2(KM.Duration.totalps,
          xlab="duration in months", ylab="survival prob",
          conf="none",
          label.curves = TRUE,                    
          time.inc=12,
          levels.only  = FALSE,                   
          n.risk=TRUE)

代码：

survplot2 <- function (fit, xlim, ylim, xlab, ylab, time.inc, conf = c("bands", "bars", "diffbands", "none"),
          add = FALSE, label.curves = TRUE, 
          abbrev.label = FALSE, levels.only = FALSE, lty, lwd = par("lwd"), 
          col = 1, col.fill = gray(seq(0.95, 0.75, length = 5)), loglog = FALSE, 
          fun, n.risk = FALSE, logt = FALSE, dots = FALSE, dotsize = 0.003, 
          grid = NULL, srt.n.risk = 0, sep.n.risk = 0.056, adj.n.risk = 1, 
          y.n.risk, cex.n.risk = 0.6, pr = FALSE, ...) {
  conf <- match.arg(conf)
  polyg <- ordGridFun(grid = FALSE)$polygon
  conf.int <- fit$conf.int
  if (!length(conf.int) | conf == "none") 
    conf.int <- 0
  opar <- par(c("mar", "xpd"))
  on.exit(par(opar))
  fit.orig <- fit
  units <- fit$units
  if (!length(units)) 
    units <- "Day"
  maxtime <- fit$maxtime
  if (!length(maxtime)) 
    maxtime <- max(fit$time)
  mintime <- min(fit$time, 0)
  pret <- pretty(c(mintime, maxtime))
  maxtime <- max(pret)
  mintime <- min(pret)
  if (missing(time.inc)) {
    time.inc <- switch(units, Day = 30, Month = 1, Year = 1, 
                       (maxtime - mintime)/10)
    if (time.inc > maxtime) 
      time.inc <- (maxtime - mintime)/10
  }
  if (n.risk && !length(fit$n.risk)) {
    n.risk <- FALSE
    warning("fit does not have number at risk\nIs probably from a parametric model\nn.risk set to F")
  }
  trans <- loglog | !missing(fun)
  if (missing(ylab)) {
    if (loglog) 
      ylab <- "log(-log Survival Probability)"
    else if (trans) 
      ylab <- ""
    else ylab <- "Survival Probability"
  }
  if (loglog) 
    fun <- function(w) logb(-logb(ifelse(w == 0 | w == 1, 
                                         NA, w)))
  else if (!trans) 
    fun <- function(w) w
  if (missing(xlab)) {
    if (logt) 
      xlab <- paste("log Survival Time in ", units, "s", 
                    sep = "")
    else xlab <- if (units == " ") 
      ""
    else paste(units, "s", sep = "")
  }
  if (missing(xlim)) 
    xlim <- if (logt) 
      logb(c(maxtime/100, maxtime))
  else c(mintime, maxtime)
  if (trans) {
    fit$surv <- fun(fit$surv)
    fit$surv[is.infinite(fit$surv)] <- NA
    if (conf.int > 0) {
      fit$lower <- fun(fit$lower)
      fit$upper <- fun(fit$upper)
      fit$lower[is.infinite(fit$lower)] <- NA
      fit$upper[is.infinite(fit$upper)] <- NA
      if (missing(ylim)) 
        ylim <- range(c(fit$lower, fit$upper), na.rm = TRUE)
    }
    else if (missing(ylim)) 
      ylim <- range(fit$surv, na.rm = TRUE)
  }
  else if (missing(ylim)) 
    ylim <- c(0, 1)
  if (length(grid)) {
    dots <- FALSE
    if (is.logical(grid)) 
      grid <- if (grid) 
        gray(0.8)
    else NULL
  }
  if (logt | trans) {
    dots <- FALSE
    grid <- NULL
  }
  olev <- slev <- names(fit$strata)
  if (levels.only) 
    slev <- gsub(".*=", "", slev)
  sleva <- if (abbrev.label) 
    abbreviate(slev)
  else slev
  ns <- length(slev)
  slevp <- ns > 0
  labelc <- is.list(label.curves) || label.curves
  if (!slevp) 
    labelc <- FALSE
  ns <- max(ns, 1)
  y <- 1:ns
  stemp <- if (ns == 1) 
    rep(1, length(fit$time))
  else rep(1:ns, fit$strata)
  if (n.risk | (conf.int > 0 & conf == "bars")) {
    stime <- seq(mintime, maxtime, time.inc)
    v <- summary(fit, times = stime, print.it = FALSE)
    vs <- if (ns > 1) 
      as.character(v$strata)
  }
  xd <- xlim[2] - xlim[1]
  yd <- ylim[2] - ylim[1]
  if (n.risk && !add) {
    mar <- opar$mar
    if (mar[4] < 4) {
      mar[4] <- mar[4] + 2
      par(mar = mar)
    }
  }
  lty <- if (missing(lty)) 
    seq(ns + 1)[-2]
  else rep(lty, length = ns)
  lwd <- rep(lwd, length = ns)
  col <- rep(col, length = ns)
  if (labelc || conf == "bands") 
    curves <- vector("list", ns)
  Tim <- Srv <- list()
  par(xpd = NA)
  for (i in 1:ns) {
    st <- stemp == i
    time <- fit$time[st]
    surv <- fit$surv[st]
    if (logt) 
      time <- logb(time)
    s <- !is.na(time) & (time >= xlim[1])
    if (i == 1 & !add) {
      plot(time, surv, xlab = xlab, xlim = xlim, ylab = ylab, 
           ylim = ylim, type = "n", axes = FALSE)
      mgp.axis(1, at = if (logt) 
        pretty(xlim)
        else seq(xlim[1], max(pretty(xlim)), time.inc), labels = TRUE)
      mgp.axis(2, at = pretty(ylim))
      if (dots || length(grid)) {
        xlm <- pretty(xlim)
        xlm <- c(xlm[1], xlm[length(xlm)])
        xp <- seq(xlm[1], xlm[2], by = time.inc)
        yd <- ylim[2] - ylim[1]
        if (yd <= 0.1) 
          yi <- 0.01
        else if (yd <= 0.2) 
          yi <- 0.025
        else if (yd <= 0.4) 
          yi <- 0.05
        else yi <- 0.1
        yp <- seq(ylim[2], ylim[1] + if (n.risk && missing(y.n.risk)) 
          yi
          else 0, by = -yi)
        if (dots) 
          for (tt in xp) symbols(rep(tt, length(yp)), 
                                 yp, circles = rep(dotsize, length(yp)), inches = dotsize, 
                                 add = TRUE)
        else abline(h = yp, v = xp, col = grid, xpd = FALSE)
      }
    }
    tim <- time[s]
    srv <- surv[s]
    if (conf.int > 0 && conf == "bands") {
      blower <- fit$lower[st][s]
      bupper <- fit$upper[st][s]
    }
    if (max(tim) > xlim[2]) {
      srvl <- srv[tim <= xlim[2] + 1e-06]
      s.last <- srvl[length(srvl)]
      k <- tim < xlim[2]
      tim <- c(tim[k], xlim[2])
      srv <- c(srv[k], s.last)
      if (conf.int > 0 && conf == "bands") {
        low.last <- blower[time <= xlim[2] + 1e-06]
        low.last <- low.last[length(low.last)]
        up.last <- bupper[time <= xlim[2] + 1e-06]
        up.last <- up.last[length(up.last)]
        blower <- c(blower[k], low.last)
        bupper <- c(bupper[k], up.last)
      }
    }
    if (logt) {
      if (conf %nin% c("bands", "diffbands")) 
        lines(tim, srv, type = "s", lty = lty[i], col = col[i], 
              lwd = lwd[i])
      if (labelc || conf %in% c("bands", "diffbands")) 
        curves[[i]] <- list(tim, srv)
    }
    else {
      xxx <- c(mintime, tim)
      yyy <- c(fun(1), srv)
      if (conf %nin% c("bands", "diffbands")) 
        lines(xxx, yyy, type = "s", lty = lty[i], col = col[i], 
              lwd = lwd[i])
      if (labelc || conf %in% c("bands", "diffbands")) 
        curves[[i]] <- list(xxx, yyy)
    }
    if (pr) {
      zest <- rbind(time[s], surv[s])
      dimnames(zest) <- list(c("Time", "Survival"), rep("", 
                                                        sum(s)))
      if (slevp) 
        cat("\nEstimates for ", slev[i], "\n\n")
      print(zest, digits = 3)
    }
    if (conf.int > 0) {
      if (conf == "bands") {
        if (logt) 
          polyg(x = c(tim, max(tim), rev(tim)), y = c(blower, 
                                                      rev(bupper), max(bupper)), col = col.fill[i], 
                type = "s")
        else polyg(x = c(max(tim), tim, rev(c(tim, max(tim)))), 
                   y = c(fun(1), blower, rev(c(fun(1), bupper))), 
                   col = col.fill[i], type = "s")
      }
      else if (conf == "diffbands") 
        survdiffplot(fit.orig, conf = conf, fun = fun)
      else {
        j <- if (ns == 1) 
          TRUE
        else vs == olev[i]
        tt <- v$time[j]
        ss <- v$surv[j]
        lower <- v$lower[j]
        upper <- v$upper[j]
        if (logt) 
          tt <- logb(ifelse(tt == 0, NA, tt))
        tt <- tt + xd * (i - 1) * 0.01
        errbar(tt, ss, upper, lower, add = TRUE, lty = lty[i], 
               col = col[i])
      }
    }
    if (n.risk) {
      j <- if (ns == 1) 
        TRUE
      else vs == olev[i]
      tt <- v$time[j]
      nrisk <- v$n.risk[j]
      tt[1] <- xlim[1]
      if (missing(y.n.risk)) 
        y.n.risk <- ylim[1]
      yy <- y.n.risk + yd * (ns - i) * sep.n.risk
      nri <- nrisk
      nri[tt > xlim[2]] <- NA
      text(tt[1], yy, nri[1], cex = cex.n.risk, adj = adj.n.risk, 
           srt = srt.n.risk)


      ## add condition here
      if (length(nri) > 1)
        text(tt[-1], yy, nri[-1], cex = cex.n.risk, adj = 1)


      if (slevp) 
        text(xlim[2] + xd * 0.025, yy, adj = 0, sleva[i], 
             cex = cex.n.risk)
    }
  }
  if (conf %in% c("bands", "diffbands")) 
    for (i in 1:ns) lines(curves[[i]][[1]], curves[[i]][[2]], 
                          lty = lty[i], lwd = lwd[i], col = col[i], type = "s")
  if (labelc) 
    labcurve(curves, sleva, type = "s", lty = lty, lwd = lwd, 
             opts = label.curves, col. = col)
  invisible(slev)
}