R/uno.R

In Alanocallaghan/surverify: Calculation of survival concordance measures.

# X = min(time, censor)
# Z = covariate
# D / delta = censoring (1 = event)

# tau <- max(test_surv[, 1])

# D_i km(x_i)^(-2) i(x_i < x_j, x_i < tau) * i(pred_i > pred_j) / 
#   D_i km(x_i)^(-2) i(x_i < x_j, x_i < tau)



uno_c <- function(
    training,
    test,
    predictors,
    upper_limit = max(test[, 1])
  ) {

  stopifnot(inherits(training, "Surv"))
  stopifnot(inherits(test, "Surv"))

  time <- training[, 1]
  event <- 1 - training[, 2]

  time_new <- test[, 1]
  event_new <- test[, 2]
  
  n <- length(time)
  n_new <- length(time_new)
  n_lp <- length(predictors)
  n_tau <- length(upper_limit)

  stopifnot(length(test) == length(predictors))

  surv <- km(time, event)
  g <- stepfun(time_new, surv, time)

  num <- denom <- 0
  # todo: vectorise for multiple taus
  for (i in 1:n_new) {
    for (j in 1:n_new) {
      if (time_new[i] < time_new[j] && g[i] == 1) {
        denom <- denom + (1 / g[i]^2) * 
          event_new[i] *
          (time_new[i] < upper_limit)

        num <- num + (1 / g[i]^2) *
          event_new[i] *
          (time_new[i] < upper_limit) *
          (predictors[i] > predictors[j])
      }
    }
  }
  unname(num / denom)
}


stepfun <- function(time_new, surv, time) {
  surv_new <- numeric(length(time_new))
  for (j in 1:length(time_new)) {
    optim <- 1
    for (i in length(surv):1) {
      if (optim && time_new[j] >= time[i]) {
        surv_new[j] <- surv[i];
        optim <- 0;
      }
      if (optim) {
        surv_new[j] <- 1;
      }
    }
  }
  surv_new
}

km <- function(time, status) {
  n <- length(time)
  surv <- numeric(n)
  current <- 1
  order <- order(time)
  time <- time[order]
  status <- status[order]
  for (i in 1:n) {
    dead <- 0
    at_risk <- 0
    for (j in 1:n) {
      at_risk <- at_risk + (time[i] <= time[j])
      dead <- dead + (time[i] == time[j]) && (status[i])
    }
    current = current * (1.0 - dead / at_risk)
    surv[i] = current
  }
  surv
}



# Uno specificity
# param spec double, vector, which length is length(thres)*(length(t)+1)
# param thres double, vector  
# param t vector of the time points
# param marker - vector, linear predictor
# param new_data - vector of the new data
# param n.th - length of the vector thres
# param n.t - length of the vector t
# param n_new_data - length of the vector new_data

spec_uno <- function(new_data, marker, t) {

  thresh <- sort(unique(marker))
  n_th <- length(thresh)
  n_t <- length(t)
  d1 <- new_data[, 1]
  n_new_data <- length(new_data)
  spec <- numeric(n_t * (n_th + 1))
  for (k in 2:(n_th + 1)) {
    for (j in 1:n_t) {
      Ivec_zsp=0.0
      Ivec_nsp=0.0
      for (i in 1:n_new_data) {
        Ivec_zsp <- Ivec_zsp +
          (marker[i] <= thresh[k - 1]) *
          (t[j] < d1[i])
        Ivec_nsp <- Ivec_nsp +
          (t[j] < d1[i])
      }
      if (Ivec_nsp != 0.0) {
        spec[((k - 1) * (n_t) + j)] <- Ivec_zsp / Ivec_nsp
      } else {
        spec[((k - 1) * (n_t) + j)] <- 0.0
      }
    }
  }
  matrix(spec, nrow = n_t, ncol = n_th + 1)
}

# Uno sensitivity
# param sens double, vector, which length is length(thres)*(length(t)+1)
# param surv double, surival vector
# param surv_time double, surival vector
# param thres double, vector  
# param t vector of the time points
# param marker - vector, linear predictor
# param new_data - vector of the new data
# param n_th - length of the vector thres
# param n_t - length of the vector t
# param n_new_data - length of the vector new_data
sens_uno <- function(
    surv,
    surv_new,
    marker,
    times) {

  thresh <- sort(unique(marker))
  n_th <- length(thresh)
  n_t <- length(times)
  n_new_data <- length(surv_new)
  n_surv <- length(surv)
  surv_time <- surv[, 1]
  status <- surv[, 2]
  new_time <- surv_new[, 1]
  new_status <- surv_new[, 2]
  
  tmp_Ivec_nse=0.0

  time <- surv[, 1]
  event <- 1 - surv[, 2]

  surv <- km(time, event)
  g <- stepfun(new_time, surv, time)
  sens <- rep(1, n_t * (n_th + 1))

  for (k in 2:(n_th + 1)) {
    for (j in 1:n_t) {
      Ivec_zse = 0.0
      Ivec_nse = 0.0
      for (i in 1:n_new_data) {
        if (times[j] >= new_time[i]) {
          if (marker[i] > thresh[k - 1]){
            Ivec_zse <- Ivec_zse + new_status[i] / g[i];
          }
          Ivec_nse <- Ivec_nse + new_status[i] / g[i];
        }
      }
      if (Ivec_nse > .Machine$double.eps) {
        sens[(k - 1) * (n_t) + j] = Ivec_zse / Ivec_nse;
      } else {
        sens[(k - 1) * (n_t) + j] = 0.0;
      }
    }
  }
  matrix(sens, n_t, n_th+1)
}