Plotting survival curves after multiple imputation

3 messages
Open this post in threaded view
|

Plotting survival curves after multiple imputation

 I am working with some survival data with missing values. I am using the mice package to do multiple imputation. I have found code in this thread which handles pooling of the MI results: https://stat.ethz.ch/pipermail/r-help/2007-May/132180.htmlNow I would like to plot a survival curve using the pooled results. Here is a reproducible example: require(survival) require(mice) set.seed(2) dt <- colon fit <- coxph(Surv(time,etype)~rx + sex + age, data=colon) dummy <- data.frame(sex=c(1,1,1),rx=c("Obs","Lev","Lev+5FU"),age=c(40,40,40)) plot(survfit(fit, newdata=dummy) ) # now create some missing values in the data dt <- colon dt\$rx[sample(1:nrow(dt),50)] <- NA dt\$sex [sample(1:nrow(dt),50)] <- NA dt\$age[sample(1:nrow(dt),50)] <- NA imp <-mice(dt) fit.imp <- coxph.mids(Surv(time,etype)~rx + sex + age,imp) # Note, this function is defined below... imputed=summary.impute(pool.impute(fit.imp)) print(imputed) # now, how to plot a survival curve with the pooled results ? ########## begin code from linked thread above coxph.mids <- function (formula, data, ...) {      call <- match.call()      if (!is.mids(data)) stop("The data must have class mids")      analyses <- as.list(1:data\$m)      for (i in 1:data\$m) {        data.i        <- complete(data, i)        analyses[[i]] <- coxph(formula, data = data.i, ...)      }      object <- list(call = call, call1 = data\$call,                     nmis = data\$nmis, analyses = analyses)      return(object) } pool.impute <- function (object, method = "smallsample") {      if ((m <- length(object\$analyses)) < 2)        stop("At least two imputations are needed for pooling.\n")      analyses <- object\$analyses      k     <- length(coef(analyses[[1]]))      names <- names(coef(analyses[[1]]))      qhat  <- matrix(NA, nrow = m, ncol = k, dimnames = list(1:m,names))      u     <- array(NA, dim = c(m, k, k),                     dimnames = list(1:m, names, names))      for (i in 1:m) {        fit       <- analyses[[i]]        qhat[i, ] <- coef(fit)        u[i, , ]  <- vcov(fit)      }      qbar <- apply(qhat, 2, mean)      ubar <- apply(u, c(2, 3), mean)      e <- qhat - matrix(qbar, nrow = m, ncol = k, byrow = TRUE)      b <- (t(e) %*% e)/(m - 1)      t <- ubar + (1 + 1/m) * b      r <- (1 + 1/m) * diag(b/ubar)      f <- (1 + 1/m) * diag(b/t)      df <- (m - 1) * (1 + 1/r)^2      if (method == "smallsample") {        if( any( class(fit) == "coxph" ) ){          ### this loop is the hack for survival analysis ###          status   <- fit\$y[ , 2]          n.events <- sum(status == max(status))          p        <- length( coefficients( fit )  )          dfc      <- n.events - p        } else {          dfc <- fit\$df.residual        }        df <- dfc/((1 - (f/(m + 1)))/(1 - f) + dfc/df)      }      names(r) <- names(df) <- names(f) <- names      fit <- list(call = call, call1 = object\$call, call2 = object\$call1,                  nmis = object\$nmis, m = m, qhat = qhat, u = u,                  qbar = qbar, ubar = ubar, b = b, t = t, r = r, df = df,                  f = f)      return(fit) } summary.impute <- function(object){       if (!is.null(object\$call1)){         cat("Call: ")         dput(object\$call1)       }       est  <- object\$qbar       se   <- sqrt(diag(object\$t))       tval <- est/se       df   <- object\$df       pval <- 2 * pt(abs(tval), df, lower.tail = FALSE)       coefmat <- cbind(est, se, tval, pval)       colnames(coefmat) <- c("Estimate", "Std. Error",                                            "t value", "Pr(>|t|)")       cat("\nCoefficients:\n")       printCoefmat( coefmat, P.values=T, has.Pvalue=T, signif.legend=T )       cat("\nFraction of information about the coefficients                                       missing due to nonresponse:","\n")       print(object\$f)       ans <- list( coefficients=coefmat, df=df,                    call=object\$call1, fracinfo.miss=object\$f )       invisible( ans )   } ### end code from linked thread above ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.