R Functions for Competing Risk Analysis and Stratified Analysis in Reduce

Stratified logrank Test

library(survival)
f <- survdiff(S ~ treatment + strata(cluster))
chisq <- f$chisq
Pval <- 1 - pchisq(chisq, 1)

Cumulative Incidence Function Plots for Competing Risks

library(cmprsk)
set.seed(1)
n <- 200
treat <- sample(c('a','b'), n, TRUE)
event <- sample(c('death','heart attack','censored'), n, TRUE)
table(event)
time.to.event <- runif(n, .5, 10)
f <- cuminc(time.to.event, event, treat, cencode='censored', rho=0) 
# also has strata argument which doesn't affect plots but affects the test

# Modification of plot.cuminc to allow vector lwd argument
plot.cuminc <- function (x, main = " ", curvlab, ylim = c(0, 1), xlim, wh = 2, 
    xlab = "Years", ylab = "Probability", lty = 1:length(x), 
    color = 1, lwd = rep(1, length(x)), ...) 
{
    if (!is.null(x$Tests)) 
        x <- x[names(x) != "Tests"]
    nc <- length(x)
    if (length(lty) < nc) 
        lty <- rep(lty[1], nc)
    else lty <- lty[1:nc]
    if (length(color) < nc) 
        color <- rep(color[1], nc)
    else color <- color[1:nc]
    lwd <- rep(lwd, length.out=nc)
    if (missing(curvlab)) {
        if (mode(names(x)) == "NULL") {
            curvlab <- as.character(1:nc)
        }
        else curvlab <- names(x)[1:nc]
    }
    if (missing(xlim)) {
        xmax <- 0
        for (i in 1:nc) {
            xmax <- max(c(xmax, x[[i]][[1]]))
        }
        xlim <- c(0, xmax)
    }
    plot(x[[1]][[1]], x[[1]][[2]], type = "n", ylim = ylim, xlim = xlim, 
        main = main, xlab = xlab, ylab = ylab, bty = "l", ...)
    if (length(wh) != 2) {
        wh <- c(xlim[1], ylim[2])
    }
    legend(wh[1], wh[2], legend = curvlab, col = color, lty = lty, lwd=lwd,
        bty = "n", bg = -999999, ...)
    for (i in 1:nc) {
        lines(x[[i]][[1]], x[[i]][[2]], lty = lty[i], col = color[i], 
              lwd=lwd[i],...)
    }
}
plot(f) # use all defaults
plot(f, xlab='Years', ylab='Cumulative Incidence', color=gray(c(0, .7, 0, .7)), lwd=c(.5, 3, .5, 3),
     lty=c(1,1,3,3))
Pvalues <- f$Tests[,'pv']
Pvalues
Pvalues['heart attack']
f$Tests[,'stat']
g <- timepoints(f, c(2,4))  # get estimates at 2 and 4 years

# Get approximate 0.95 confidence intervals for 2- and 4-y cumulative incidence
g$est - 1.96*sqrt(g$var)
g$est + 1.96*sqrt(g$var)
Topic attachments
I Attachment Action Size Date Who Comment
cuminc.pngpng cuminc.png manage 8.0 K 19 Sep 2006 - 17:15 FrankHarrell Example cumulative incidence plot, 2 events, 2 treatments
Topic revision: r1 - 19 Sep 2006, FrankHarrell
 

This site is powered by FoswikiCopyright © 2013-2020 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding Vanderbilt Biostatistics Wiki? Send feedback