r - error message when ploting subjects at risk with survplot -
i following error message when trying plot subjects @ risk along x-axis in survplot:
error in text.default(tt[-1], yy, nri[-1], cex = cex.n.risk, adj = 1) : zero-length 'labels' specified
any help? rather new survival analysis, , failed find explanation of error. code in general seems fine, except when add n.risk=true option plot, error comes up. hint appreciated. many thanks.
below data code used.
here data
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")
the creation of survival object
library(survival) library(rms) duration$survobj <- with(duration, surv(peacedur, peacefail==1))
the fit + converstion npsurv
km.duration.totalps <- survfit(survobj ~ totalps, data = duration, conf.type = "log-log") class(km.duration.totalps) <- c(class(km.duration.totalps), "npsurv")
the plot:
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)
totalps=0
has 1 event. n.risk
doesn't that.
running survplot(km.duration.totalps[-1], ...)
works setting 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)
digging deeper root of problem, line in rms:::survplot.npsurv
text(tt[-1], yy, nri[-1], cex = cex.n.risk, adj = 1)
your nri totalps=0 vector of length 1, r attempting plot
text(1, 1, integer(0))
try , same error. remedy that, either use solutions above (plotting totalps=0 not interesting anyway since straight line), or can edit source code below, inserting if statement near end of code trick if (length(nri) > 1)
so can use new function full table/plot without errors (i wouldn't, though, cause can see, label covers @ risk table)
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)
code:
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 not have number @ risk\nis parametric model\nn.risk set 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) (i in 1:ns) { st <- stemp == 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, @ = if (logt) pretty(xlim) else seq(xlim[1], max(pretty(xlim)), time.inc), labels = true) mgp.axis(2, @ = pretty(ylim)) if (dots || length(grid)) { xlm <- pretty(xlim) xlm <- c(xlm[1], xlm[length(xlm)]) xp <- seq(xlm[1], xlm[2], = 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, = -yi) if (dots) (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 ", 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")) (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) }
Comments
Post a Comment