R/CoxMI.R

In SurvMI: Multiple Imputation Method in Survival Analysis

#Imports:
#  survival
#'
#' This function calculate the Cox model with MI when event uncertainty presents
#'
#' @param data the data list which has been transformed by the transformation function
#' @param nMI number of imputation
#' @param covariates list of covariates on the RHS of Cox model
#' @return A class object with MI Cox estimation for given data
#' @export
#' @examples

##main function
##Derive the time and status variables from the long formatted data generated by the data transform function
CoxMI<-function(data_list,nMI=1000,covariates=NULL,id=NULL,...){
  form1="Surv(time,cens)~"
  form2=paste(covariates, collapse="+")


if (is.null(id)){
  form<-as.formula(paste("survival::",paste(form1,form2)))
}
if (!is.null(id)){
  form1="survival::Surv(start,end,cens)~"
form<-as.formula(paste(form1,form2))
}

data1<-data_list$data_uc

formc=as.formula(paste("~",form2))
x_beta<-model.matrix(formc,data1)[,-1]
if (!is.null(dimnames(x_beta)[[2]])){
  p1=length(dimnames(x_beta)[[2]])
}

if (is.null(dimnames(x_beta)[[2]])){
  p1=length(covariates)
  x_beta<-matrix(x_beta)
}

n<-nrow(data1)

  betamat=matrix(0,nrow=nMI,ncol=p1)
  vararray=array(0,dim=c(nMI,p1,p1))
  betamean=rep(0,p1)
  censor=rep(0,n)
  varmean=matrix(0,p1,p1)
  expected=0
  for (isim in 1:nMI) {
      if (is.null(id)){
        ##two derived variables used in Cox model
        data1$time=rep(0,n)
        data1$cens=rep(0,n)
    for (i in 1:n) {
      r<-runif(1)
      #generate censor indicator, if cum prob > runif[1] then event
      j=which(cumsum(data_list$weights[[i]])>r)[1]
      censor[i]=j
      data1$time[i]=data_list$time[[i]][j]

      if (censor[i]<data_list$e[i] ) {data1$cens[i]=1}
      else if (censor[i]>=data_list$e[i]) {data1$cens[i]=0}
    }

        coxph1=survival::coxph(form,data=data1, ...)
        expected=sum(data1$cens)+expected
      }
      ##for multiple event case
      if (!is.null(id)){
        start=end=censor=list()
        for (i in 1:n) {
          r<-runif(data_list$e[[i]])
          #generate censor indicator, if cum prob > runif[1] then event
          j=which(data_list$prob[[i]]>r)
          end_temp=c(data_list$time[[i]][j],data_list$time[[i]][data_list$e[[i]]])
          start_temp=c(0,end_temp[-length(end_temp)]+1e-06)
          censor_temp=c(rep(1,length(j)),0)
          data1$ntimes[i]=as.numeric(length(j)+1)
        start=c(start,list(start_temp))
        end=c(end,list(end_temp))
        censor=c(censor,list(censor_temp))
        }

        data2 <- as.data.frame(lapply(data1, rep, data1$ntimes))

        length<-length(unlist(start))
          #trt_long<-c(rep(0,length(unlist(x))))
          data2$start<-NA
          data2$end<-NA
          data2$cens<-NA
          cum<-0
          for (i in (1:n)){

            index<-cum+length(unlist(start[[i]]))
            data2$start[(cum+1):index]=matrix(unlist(start[[i]]))
            data2$end[(cum+1):index]=matrix(unlist(end[[i]]))
            data2$cens[(cum+1):index]=matrix(unlist(censor[[i]]))
            cum=index
          }
          id_var=data2[,c(id)]
          #coxph1=survival::coxph(form, cluster = id_var, data=data2)
          coxph1=survival::coxph(form, cluster = id_var, data=data2, ...)
          expected=sum(data2$cens)+expected
      }


    betamat[isim,]=coxph1$coef
    vararray[isim,,]=coxph1$var
    betamean=betamean+coxph1$coef
    varmean=varmean+coxph1$var
  }

  varmean=varmean/nMI
  x=betamean/nMI
  B=matrix(0,p1,p1)
  for (i in 1:nMI) {
    Bi=matrix(betamat[i,]-x,nrow=p1,ncol=1)
    B=B+Bi%*%t(Bi)
  }
  B=B/(nMI-1)
  var_MI<-varmean+((nMI+1)/nMI)*B
  p=as.numeric(2*(1-pnorm(abs(x/sqrt(diag(var_MI))))))

  res<-list(est=x,pvalue=p,var=var_MI,est_matrix=betamat,Var_mat=vararray,
            between_var=B, within_var=varmean,nMI=nMI,
            column_name=dimnames(x_beta)[[2]],n=n,en=expected/(nMI))
  attr(res,"est")=x
  attr(res,"var")=var_MI
  attr(res,"betamat")=betamat
  attr(res,"Var_matrix")=vararray
  attr(res,"Between Var")=B
  attr(res,"Within Var")=varmean
  attr(res,"nMI")=nMI
  attr(res,"pvalue")=p
  class(res) <- c("coxMI", class(res))
  return(res)
  }