[R] Plotting survival curves after multiple imputation

[Frank Harrell]

I haven't seen anyone solve this.  I think it would be reasonable to do a
time point by time point averaging (over multiple imputations) of the
underlying survival curve, although there is some question about whether to
"freeze" the centering constant (sum of beta times covariate mean).  What
will be harder to get is confidence intervals for predicted probabilities
after multiple imputation.  I hope someone will take up the challenge.
Frank

Robert Long wrote
> 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.html
> 
> Now 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
> 
> ______________________________________________

> R-help@

>  mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.





-----
Frank Harrell
Department of Biostatistics, Vanderbilt University
--
View this message in context: http://r.789695.n4.nabble.com/Plotting-survival-curves-after-multiple-imputation-tp4658552p4659505.html
Sent from the R help mailing list archive at Nabble.com.