R/simex.auc.naive.RM.R

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

simex.auc.naive.RM <- function(Y2,Y1,
                               K=1000,LAM1=seq(0.1,1.9,0.2),LAM2=seq(0.1,1.9,0.2),
                               tau1.reml=NA,tau2.reml=NA){
  # Y1, Y2 are matrices of RM (Y2 is the diseased pop.)
  # K is the resampling size for the SIMEX procedure
  # LAM1,2 are the vectors of Lambda for simex resampling
  # can accept REML estimates of variance components
  # otherwise, uses ANOVA estimates
  n <- nrow(Y1)
  p <- ncol(Y1)
  # ANOVA estimates in case REML not supplied
  if(is.na(tau1.reml)){
    
 
  
  y1_bar <- rowMeans(Y1)
  y2_bar <- rowMeans(Y2)
  
  # sample estimates
  
  mu1 <- mean(Y1)
  mu2 <- mean(Y2)
  
  
  SSB1 <- p * sum( (y1_bar-mu1)^2)
  SSB2 <- p * sum( (y2_bar-mu2)^2)
  
  tau1_hat <-sqrt(( 1 / ((p-1)*n)) * sum((Y1-do.call(cbind,replicate(p,y1_bar,simplify = FALSE)))^2))
  tau2_hat <-sqrt(( 1 / ((p-1)*n)) * sum((Y2-do.call(cbind,replicate(p,y2_bar,simplify = FALSE)))^2))
  }
  
  # in case REML estimates are supplied
  else{
    tau1_hat <- tau1.reml
    tau2_hat <- tau2.reml
  }
  
  TC <- matrix(data = NA,nrow = 0,ncol = 3)
  for(lam1 in LAM1){
    for(lam2 in LAM2){
      
      w1 <- matrix(rnorm(K*p*n,0,sqrt(lam1)*tau1_hat), n, p*K )
      w2 <- matrix(rnorm(K*p*n,0,sqrt(lam2)*tau2_hat), n, p*K )
      
      Y1.rep <- do.call(cbind,replicate(K,Y1,simplify = FALSE))
      W1 <- Y1.rep + w1
      
      Y2.rep <- do.call(cbind,replicate(K,Y2,simplify = FALSE))
      W2 <- Y2.rep + w2
      
      TC <- rbind(TC,c(lam1,lam2, mw.outer.matrix(W1,W2)))
    }
  }
  
  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)
}