R/reiser.2000.R

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

reiser.2000 <- function(Y1,Y2){
  
  m <- nrow(Y1)
  n <- nrow(Y2)
  p <- ncol(Y1)
  
  y1_bar <- rowMeans(Y1)
  y2_bar <- rowMeans(Y2)
  
  mu1.hat <- mean(Y1)
  mu2.hat <- mean(Y2)
  
  
  SSB1 <- p * sum( (y1_bar-mu1.hat)^2)
  SSB2 <- p * sum( (y2_bar-mu2.hat)^2)
  
  tau1_hat <-sqrt(( 1 / ((p-1)*m)) * 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))
  
  # Reiser 2000, just above 3.2
  
  sig1 <- sqrt(((1/(n-1)) * SSB1 - tau1_hat^2)/p)
  sig2 <- sqrt (((1/(n-1)) * SSB2 - tau2_hat^2)/p)
  
  delta.hat <- pnorm((mu2.hat-mu1.hat)/(sqrt(sig1^2+sig2^2)))
  delta.hat
  
}