R/atomic.sim.R

In blebedenko/ME-AUC: Estimate AUC in ME Models - our functions

# This function is one atomic simulation
# of the normal setting
# which returns the following estimates for AUC:

# EBLUP
# REISER 2000 REML
# SIMEX "naive"
# SIMEX on averages
# SIMEX on EBLUPs

atomic.sim <- function(m,n,p,muX,muY,
                       sigmaX,sigmaY,sigmaEps,sigmaEta,
                       K,LAM1,LAM2)
  
{  
 
  
  X <- rnorm(m,muX,sigmaX)
  Y <- rnorm(n,muY,sigmaY)
  
  # additive normal ME's with var sigmaEps, sigmaEta
  eps <- matrix(rnorm(n = m*p,mean = 0, sd = sigmaEps),nrow = m,ncol = p)
  eta <- matrix(rnorm(n = n*p,mean = 0, sd = sigmaEta),nrow = n,ncol = p)
  
  x <- X + eps
  y <- Y + eta
  
  # eblup
  eblup.est <- eblup.auc.REML(x,y)
  
  e <- eblup.est[[1]]     # the MW estimate from EBLUPS
  sEps <- eblup.est[[3]][1] # REML estimate of sigma_eps
  sEta <- eblup.est[[3]][2] # REML estimate of sigma_eta
  ebX <- eblup.est$blups1
  ebY <- eblup.est$blups2
  sX <- eblup.est[[2]][1] # REML estimate of sigma_X
  sY <- eblup.est[[2]][2] # REML estimate of sigma_Y
  # simex
  simex <-  simex.auc.naive.RM(y,x,K=K,LAM1 = LAM1,LAM2 = LAM2,tau1.reml = sEta,tau2.reml = sEps)
  
  s1 <- simex[1] # linear 
  s2 <- simex[2] # quadratic
  
  # Reiser
  r <- reiser.2000.REML(x,y,sX,sY)
  
  # Simex on averages - ME sd is sigma.e/sqrt(p)
  
  simex.aver <- simex.auc.reml(rowMeans(x),rowMeans(y),K = K,
                               LAM1 = LAM1,LAM2 = LAM2,sigmaX = sEps/sqrt(p),sigmaY = sEta/sqrt(p))
  
  sa1 <- simex.aver[1]
  sa2 <- simex.aver[2]
  
  # simex on eblups
  
  simex.eblups <- simex.auc.reml(ebX,ebY,K = K,LAM1 = LAM1,LAM2 = LAM2,sigmaX = sEps,sigmaY = sEta)
  
  # auc based on averages, ignoring ME
  
  auc.ave <- auc.averages(x,y)
  
  vec <- c(e,s1,s2,sa1,sa2,simex.eblups,auc.ave)
  names(vec) <- c("EBLUP","SIMEX linear","SIMEX quad","SIMEX A linear","SIMEX A quad",
                  "SIMEX EBLUP linear","SIMEX EBLUP quad","Averages MW")
  return(vec)
  
}

# 
# m <- n <- 200
# p <- 3
# muX <- 1.8
# muY <- 0
# sigmaX <- sigmaY <- 1
# sigmaEps <- sigmaEta <- 1
# K <- 1000
# LAM1 <- LAM2 <- seq(0.1,1.9,by=0.2)
# 
# a <- Sys.time()
# ww <- atomic.sim(m=m,n=n,p=p,muX=muX,muY=muY,sigmaX = sigmaX,sigmaY = sigmaY,sigmaEps = sigmaEps,sigmaEta = sigmaEta,
#            K=10,LAM1 = LAM1,LAM2=LAM2)
# 
# Sys.time()-a