Optimised pseudosurv2() function

pseudosurv2.R

pseudosurv2 <- function (time, event, tmax) 
{
  if (any(is.na(time))) 
    stop("missing values in 'time' vector")
  if (any(time < 0)) 
    stop("'time' must be nonnegative")
  if (any(is.na(event))) 
    stop("missing values in 'event' vector")
  if (any(event != 0 & event != 1)) 
    stop("'event' must be a 0/1 variable (alive/dead)")
  if (missing(tmax)) 
    tmax <- unique(time[event == 1])
  if (any(is.na(tmax))) 
    stop("missing values in 'tmax' vector")
  if (sum(tmax > max(time)) > 0) 
    stop("tmax greater than largest observation time")
  
  tmax <- sort(tmax)
  ltmax <- length(tmax)
  howmany <- length(time)
  pseudo <- data.frame(id = 1:howmany, time = time, event = event)
  pseudo <- pseudo[order(pseudo$time, -pseudo$event), ]
  
  tu <- unique(pseudo$time[pseudo$event == 1])
  
  inx <- vapply(tmax,function(x) sum(x >= tu),numeric(1))
  
  KM.omit <- surv.omit2(pseudo,inx)
  KM.tot <- surv.tot2(pseudo)[inx]
  KM.tot <- matrix(KM.tot, byrow = TRUE, nrow = howmany, ncol = length(tmax))
  pseu <- howmany * KM.tot - (howmany - 1) * KM.omit
  pseu <- as.data.frame(pseu)
  row.names(pseu) <- pseudo$id
  names(pseu) <- paste("time", 1:length(tmax), sep = "_")
  pseu <- pseu[order(pseudo$id),,drop=F]
  return(pseu)
}

surv.omit2 <- function (pseudo,inx) 
{
  td <- pseudo$time[pseudo$event == 1]
  lt.temp <- c(td[-1], td[length(td)] + 1)
  lt <- which(td != lt.temp)
  
  km_extract <- which(pseudo$event)[lt][inx]
  keep <- 1:nrow(pseudo) %in% km_extract
  
  km <- C_survomitkm(pseudo$event,keep)
  km
}


surv.tot2 <- function (pseudo, tmax) 
{
  howmany <- nrow(pseudo)
  td <- pseudo$time[pseudo$event == 1]
  lt.temp <- c(td[-1], td[length(td)] + 1)
  lt <- which(td != lt.temp)
  Y <- howmany:1
  N <- pseudo$event
  kmji <- (Y - N)/Y
  km <- cumprod(kmji)
  if (!missing(tmax)) {
    tt <- pseudo$time[pseudo$event == 1]
    tt <- tt[lt]
    tt <- c(0, tt, tmax)
    tt <- diff(tt)
  }
  km <- km[pseudo$event == 1]
  km <- km[lt]
  if (!missing(tmax)) {
    km <- sum(c(1, km) * tt)
  }
  km
}





Rcpp::cppFunction(
  'NumericMatrix C_survomitkm(LogicalVector event,
                               LogicalVector keep){
  int L = event.length();
  int K = 0;
  int M = 0;
  
  for(int i = 0; i < L; i++){
    if(keep[i] == TRUE){
      K++;
      M = i;
    }
  }
  
  
  NumericVector current_vec(L);
  NumericMatrix res(L,K);
  int iK = 0;
  
  for(int i = 0; i < M+1; i++){
      
    if((i == 0) & (event[0] == TRUE)) {
      current_vec[0] = 1;
      for(int j = 1; j < L;j++){
        current_vec[j] = 1 - 1/(1.0*L-1);
      }
    } else if((i == 0) & (event[0] == FALSE)) {
      for(int j = 0; j < L;j++){
        current_vec[j] = 1;
      }
    } else if((i == L-1) & (event[i] == TRUE)){
      for(int j = 0; j < L-1; j ++){
        current_vec[j] = 0;
      }
    } else if(event[i] == TRUE) {
      for(int j = 0;j < i; j++){
        current_vec[j] = current_vec[j]*(1-1/(1.0*L-1.0*i));
      }
      for(int j = i+1; j < L; j++){
        current_vec[j] = current_vec[j]*(1-1/(1.0*L-1.0*i-1));
      }
    }
    if(keep[i] == TRUE){
      res(_,iK) = current_vec;
      iK++;
    }
  }
  return res;
}')