In elong0527/simuCPP: Simulate Concordance Probability, AUC and Time-Dependent AUC Based on Copula

library(simuCCP)
library(mvtnorm)
# C-index functon defined in evacure::cindex
cindex <- function (risk, pi, time, delta, tau = NULL){
    risk <- as.numeric(risk)
    pi <- as.numeric(pi)
    n <- length(risk)
    t1 <- outer(time, time, "<")
    t2 <- rep(1, n)
    if (!is.null(tau)) 
        t2 <- time < tau
    r1 <- outer(risk, risk, function(x, y) (x > y) + 0.5 * (x == 
        y))
    p1 <- outer(pi, pi, "*")
    num <- delta * t1 * t2 * r1 * p1
    diag(num) <- 0
    dev <- delta * t1 * t2 * p1
    diag(dev) <- 0
    sum(num)/sum(dev)
}

Numerical Example:

C-index

We provide one example of simulating correlated markers with pre-defined C-index discussed in the paper. The code below consider two markers with Gaussian and Clayton copulas (family = c(1,3)) with C-index both equal to 0.75. The conditional copula between two markers is Gumbel (fam2 = 4) with parameter 1 (par2 = 1). To reduce running time, we only repeat simulation 10 times instead of 1,000 times as in the paper.

set.seed(123)
N_sim <- 10 # Number of Simulation
res <- list()


for(i in 1:N_sim){
  N = 100 # Sample size
  data <- simuCCP(N, value = 0.75, family = c(1,3), fam2 = 4, par2 = 1)


  #Assume a exponential marginal distribution of survival time T with mean 1.
  T <- exp( data$data[,1] )   

  # Assume a standard normal marginal distribution of marker X1 and marker X2
  X1 <- (- qnorm(data$data[,2]) ) 
  X2 <- (- qnorm(data$data[,3]))

  pi <- rep(1, N)
  delta <- rep(1,N)

  c_X1 <- cindex(risk = X1, pi = rep(1, length(T) ), time = T, delta = rep(1, length(T)) )
  c_X2 <- cindex(risk = X1, pi = rep(1, length(T) ), time = T, delta = rep(1, length(T)) )
  cor_X1_X2 <- cor(X1,X2, method = "kendall")
  res[[i]] <- c(c_X1, c_X2, cor_X1_X2)

}
res <- do.call(rbind, res)
res_mean <- apply(res, 2, mean)
res_sd   <- apply(res, 2, sd)

simu_cindex_res <- rbind(res_mean, res_sd)
colnames(simu_cindex_res) <- c("X1", "X2", "Correlation")
round(simu_cindex_res, 2)

Cumulative time dependent AUC

library(ROCR)
set.seed(123)
N_sim <- 10 # Number of Simulation
res <- list()


for(i in 1:N_sim){
  N = 100 # Sample size
  data <- simuCCP(N, value = 0.75, metric = "AUC", family = c(1,3), 
                  fam2 = 4, par2 = 1, c0 = 0.5)


  #Assume a exponential marginal distribution of survival time T with mean 1.
  T <- exp( data$data[,1] )   

  # Assume a standard normal marginal distribution of marker X1 and marker X2
  X1 <- (- qnorm(data$data[,2]) ) 
  X2 <- (- qnorm(data$data[,3]))

  pi <- rep(1, N)
  delta <- rep(1,N)
  auc_X1 <- ROCR::performance(prediction(X1, T < median(T)), "auc")@y.values
  auc_X2 <- ROCR::performance(prediction(X2, T < median(T)), "auc")@y.values
  cor_X1_X2 <- cor(X1, X2, method = "kendall")
  res[[i]] <- unlist(c(auc_X1, auc_X2, cor_X1_X2))

}
res <- do.call(rbind, res)
res_mean <- apply(res, 2, mean)
res_sd   <- apply(res, 2, sd)

simu_cindex_res <- rbind(res_mean, res_sd)
colnames(simu_cindex_res) <- c("X1", "X2", "Correlation")
round(simu_cindex_res, 2)

elong0527/simuCPP documentation built on March 29, 2021, 10:03 a.m.