OmnibusParadox

library(MASS)
library(ellipse)

#Helper function to make colors more transluscent
add.alpha <- function(col, alpha=1){#from https://www.r-bloggers.com/how-to-change-the-alpha-value-of-colours-in-r/
  if(missing(col))
    
    stop("Please provide a vector of colours.")
  
  apply(sapply(col, col2rgb)/255, 2, 
        
        function(x) 
          
          rgb(x[1], x[2], x[3], alpha=alpha))  
}

#Colors I like
salmon <- add.alpha("salmon2",.35)
maroon <- add.alpha("maroon4",.8)
Grey   <- add.alpha("grey20",.4)

#Helper function to compute mahalanobis dist. for vector v & matrix A
#Reduces to Euclidean dist. when A is proportional to an identity matrix
quad.form <- function(v,A){
  
  Q <- v%*%A%*%v
  
  return(Q)
}


#Now we can do our simulation
set.seed(1337)
r <- 0                   #correlation between tests, 0=indep.
k <- 3e4                 #replicates
n <- 2                   #number of tests
test.z <- c(1.75,1.85)   #example z values for r=0
#test.z <- c(2.05, 2.05)  #example z values for r=.5
#test.z <- c(2.25,2.25)    #example z values for r=.8


mu <- rep(0,n)                               #Mean vector
sigma <- diag(1-r,n)+ matrix(rep(r,n^2),n,n) #Covariance matrix
tau <- solve(sigma)                          #Invert cov. matrix
#tau <- 1/(1-r^2)*matrix(c(1,-r,-r,1),nrow=2)  #Precise tau for 2 Z tests 


Z <- mvrnorm(k,mu,sigma)         # generate our z-values
D2 <- apply(Z,1,quad.form,A=tau) # compute our omnibus test statistic 

crit <- qchisq(.95,df=n)         # critical value of chi-square

#Plot the sampling distribution and add ellipse/circle and box and arrow
plot(Z[,1],Z[,2],xlim=c(-3.2,3.2),ylim=c(-3.2,3.2),col=salmon,pch=20,bty="n",xlab="Z1",ylab="Z2")
lines(ellipse(sigma,level=.95),col=maroon,lwd=5,lty=1)
lines(c(-1.96,-1.96,-1.96,1.96,1.96,1.96,1.96,-1.96),
       c(-1.96,1.96,1.96,1.96,1.96,-1.96,-1.96,-1.96),lwd=3,col=Grey)
#abline(v=c(-1.96,1.96),col=Grey,lwd=3)         #Extend the lines of the box
#abline(h=c(-1.96,1.96),col=Grey,lwd=3)         #Extend the lines of the box
arrows(0,0,test.z[1],test.z[2],col="grey30",lwd=5)
points(0,0,pch=3,col="red",lwd=3)