outCppE2 <- function(n = 300, plot = F){
  
  A <- readRDS(paste0("Simulation Study 2/Results/simCpp_S2_n",n,".rds"))
  
  co <- complete.cases(A)
  
  M1 <- A[co,1:65]
  M2 <- A[co,66:130]
  M3 <- A[co,131:195]
  ti <- A[co,196]
  
  imu0 <- 1:8
  ipi0 <- 9:16
  imu1 <- 17:24
  ipi1 <- 25:32
  isg0 <- 33:40
  isg1 <- 41:48
  inu0 <- 49:56
  inu1 <- 57:64
  iR <- ncol(M1)
  
  MSEp <- mean((M1[,c(ipi0,ipi1)] - M2[,c(ipi0,ipi1)])^2)
  MSEm <- mean((M1[,c(imu0,imu1)] - M2[,c(imu0,imu1)])^2)
  MSEs <- mean((M1[,c(isg0,isg1)] - M2[,c(isg0,isg1)])^2)
  MSEn <- mean((M1[,c(inu0,inu1)] - M2[,c(inu0,inu1)])^2)
  MSER <- mean((M1[,iR] - M2[,iR])^2)
  
  ABBp <- mean(abs(colMeans(M1[,c(ipi0,ipi1)] - M2[,c(ipi0,ipi1)])))
  ABBm <- mean(abs(colMeans(M1[,c(imu0,imu1)] - M2[,c(imu0,imu1)])))
  ABBs <- mean(abs(colMeans(M1[,c(isg0,isg1)] - M2[,c(isg0,isg1)])))
  ABBn <- mean(abs(colMeans(M1[,c(inu0,inu1)] - M2[,c(inu0,inu1)])))
  ABBR <- mean(abs(colMeans(M1[,iR,drop = F] - M2[,iR,drop = F])))
  
  ACOp <- mean((M1[,c(ipi0,ipi1)] > M2[,c(ipi0,ipi1)] - qnorm(0.975)*M3[,c(ipi0,ipi1)]) & (M1[,c(ipi0,ipi1)] < M2[,c(ipi0,ipi1)] + qnorm(0.975)*M3[,c(ipi0,ipi1)]))
  ACOm <- mean((M1[,c(imu0,imu1)] > M2[,c(imu0,imu1)] - qnorm(0.975)*M3[,c(imu0,imu1)]) & (M1[,c(imu0,imu1)] < M2[,c(imu0,imu1)] + qnorm(0.975)*M3[,c(imu0,imu1)]))
  ACOs <- mean((M1[,c(isg0,isg1)] > M2[,c(isg0,isg1)] - qnorm(0.975)*M3[,c(isg0,isg1)]) & (M1[,c(isg0,isg1)] < M2[,c(isg0,isg1)] + qnorm(0.975)*M3[,c(isg0,isg1)]))
  ACOn <- mean((M1[,c(inu0,inu1)] > M2[,c(inu0,inu1)] - qnorm(0.975)*M3[,c(inu0,inu1)]) & (M1[,c(inu0,inu1)] < M2[,c(inu0,inu1)] + qnorm(0.975)*M3[,c(inu0,inu1)]))
  ACOR <- mean((M1[,iR] > (M2[,iR] - qnorm(0.975)*M3[,iR])) & (M1[,iR] < (M2[,iR] + qnorm(0.975)*M3[,iR])))
  
  COMPT <- mean(ti)
  
  if(plot){
    boxplot(M2, main = paste0("Simulation Study II (n: ", n,")"))
    points(colMeans(M1), col = 2, pch = 16)
    plot(apply(M2,2,sd), colMeans(M3), main = paste0("Simulation Study II (n: ", n,")"),
         xlab = "Empirical SEs", ylab = "Estimated SEs")
    abline(a = 0, b = 1, lty = 3, lwd = 2)
  }
  
  return(c(n,MSEp, MSEm, MSEs, MSEn, MSER,
           ABBp, ABBm, ABBs, ABBn, ABBR,
           ACOp, ACOm, ACOs, ACOn, ACOR, COMPT))
}
