R/AtomicSim.US.R

In blebedenko/thescript2: Boris's thesis scripts

#' Atomic simulation of AUC coverage experiment with undersmoothing
#' 
#' @param distribution one of the following strings : "normal", "gamma.skewed",  "gamma.symmetrical" ,
#' "log.normal", "bimodal1", "bimodal2", "t3", "laplace", "cauchy"
#' @param n number of observations in each group
#' @param AUC the probability P(X>Y) of correctly ranking a pair of random observations (diseased,healthy).
#' @param theta_squared contamination parameter (noise/signal)
#' @param var_eq logical - if TRUE VX=VY=1, if FALSE VX=2VY=2 ratio of the biomarker variances
#' @param ME_var_ratio (defaults to 1) the ration between ME variances, V
#' @param R bootstrap replications (defaults to 1000)
#' @param US_ratio vector of multipliers for the undersmoothing of bootstrap samples (defaults to 0.5,0.25)
#' @details This is the `doOne` function required by `simsalapar` for parallel simulations.
#' @examples AtomicSim("gamma.skewed",250,0.8,0.5,TRUE,1,R=1000)
#' @return  a vector of the point estimate and percentile CI and also the original observations,all concatenated in a named vector.
#' @export



 AtomicSim.US <- function(distribution,
                          n,
                          AUC,
                          theta_squared,
                          var_eq,
                          ME_var_ratio,
                          R=1000,US_ratio = c(0.5,0.25)){
   tictoc::tic()
   variances <- FindVariances(var_eq=var_eq,
                              ME_var_ratio=ME_var_ratio,
                              theta_squared=theta_squared)

   VX <- variances["VX"]
   VY <- variances["VY"]
   Veps <- variances["Veps"]
   Veta <- variances["Veta"]


   obs <-  GenerateData(distribution=distribution,
                      n = n ,
                      AUC = AUC,
                      theta_squared = theta_squared,
                      var_eq = var_eq,
                      ME_var_ratio = ME_var_ratio)

   x <- obs$x
   y <- obs$y
   
   Y_decon_cdf <- QPdecon::QPdecon(y,f_Z = sqrt(Veta),Pdf = FALSE) # F_Y
   X_decon_pdf <- QPdecon::QPdecon(x,f_Z = sqrt(Veps)) # f_X
   A.hat_original <- EstimateAUC.QP(x,y,MESD_X = sqrt(Veps),MESD_Y = sqrt(Veta))
   lambda_original_XY <- c(X_decon_pdf$lambda,Y_decon_cdf$lambda)
   
   Xboot <- replicate(R,sample(x,replace = TRUE))
   Yboot <- replicate(R,sample(y,replace = TRUE))
   res <- matrix(NA ,nrow = length(US_ratio), ncol = R)
   tictoc::tic()
   for(u in 1:length(US_ratio)){
      for(i in 1:R){
         
         Y_decon_cdf <- QPdecon::QPdecon(Yboot[ ,i],
                                f_Z = sqrt(Veta),
                                Pdf = FALSE,
                                lambda = lambda_original_XY[2]* US_ratio[u]) # F_Y
         X_decon_pdf <- QPdecon::QPdecon(Xboot[, i],
                                f_Z = sqrt(Veps),
                                lambda = lambda_original_XY[1]* US_ratio[u]) # f_X
         
         x <- X_decon_pdf$x
         fx <- X_decon_pdf$f_X
         delta <- mean(diff(x)) #the values vary very slightly so we just take the mean
         SUM_FYx <- Decon2CDF(qpobj = Y_decon_cdf,quant = x)
         A.hat <- sum(fx * SUM_FYx*delta)
         res[u,i] <- A.hat
      }
      
   }
   percentile <- apply(res,1,quantile,c(0.025,0.975))
   colnames(percentile) <- US_ratio
   rownames(percentile) <- c("Lp","Up")
   bca <- apply(res,1,coxed::bca)
   colnames(bca) <- colnames(percentile)   #return a vector of the point estimate and percentile CI
   #and also the original observations
   rownames(bca) <- c("Ub","Up")
   return(rbind(percentile,bca))

 }
# doCheck(AtomicSim.US,varList,nChks = 1)
# tictoc::tic()
# ll <- AtomicSim.US("gamma.skewed",250,0.8,0.5,TRUE,1,R=1000)
# tictoc::toc()
# QA passed