R/simex.auc.reml.thin.R

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

# This is the basic simex auc function
# gets input two _vectors_ and "known" ME variances
# returns the SIMEX estimate of P(X>Y)
# This thin version

simex.auc.reml <- function(X,Y,sigmaX,sigmaY,K=1000,LAM1=seq(0.1,1.9,0.2),LAM2=seq(0.1,1.9,0.2))
{                           
  
  m <- length(Y)
  n <- length(X)
  
  
  TC <- matrix(data = NA,nrow = 0,ncol = 3)
  for(lam1 in LAM1){
    for(lam2 in LAM2){
      
      # additive ME
      eX <- matrix(rnorm(K*n,0,sqrt(lam1)*sigmaX), n, K )
      eY <- matrix(rnorm(K*n,0,sqrt(lam2)*sigmaY), n, K )
      
      X.resamp <- do.call(cbind,replicate(K,X,simplify = FALSE))
      X.resamp <- X.resamp + eX
      
      Y.resamp <- do.call(cbind,replicate(K,Y,simplify = FALSE))
      Y.resamp <- Y.resamp + eY
      
      TC <- rbind(TC,c(lam1,lam2, mw.outer.matrix(X.resamp,Y.resamp)))
    }
  }
  
  TC <- as.data.frame(TC)
  names(TC) <- c("L1","L2","E_theta_N")
  
  # Now we obtain the intercept from the linear regression model
  TC$L1 <- 1+TC$L1
  TC$L2 <- 1+TC$L2
  
  TC$L1sq <- TC$L1^2
  TC$L2sq <- TC$L2^2
  
  fit1 <- lm(formula = E_theta_N ~ L1+L2,data = TC)
  fit2 <- lm(formula = E_theta_N ~ L1+L2+L1sq+L2sq,data = TC)
  
  s1 <- as.numeric(fit1$coefficients[1])
  s2 <- as.numeric(fit2$coefficients[1])
  
  # trimming of values >1
  s1 <- min(c(s1,1))
  s2 <- min(c(s2,1))
  
  v <- c(s1,s2)
  names(v) <- c("Linear","Quadratic")
  return(v)
  
  
  
}


# 
# X <- rnorm(m,muX,sigmaX)+rnorm(m,0,sigmaEps)
# Y <- rnorm(n,muY,sigmaY)+rnorm(n,0,sigmaEta)
# 
# 
# K <- 100
# a <- Sys.time()
# simex.auc.reml(X,Y,1,1)
# Sys.time()-a
# pnorm((muX-muY)/(sqrt(sigmaX^2+sigmaY^2)))
# 
# microbenchmark::microbenchmark(simex.auc.reml(X,Y,1,1),times = 3)