library(normtest)
library(nortest)
library(sn)
library(ggplot2)


#Comparing Jarque-Bera, Lilliefors, Anderson-Darling and Shapiro-Wilk normalitytests
ssize <- seq(from=100, to=500, by=100)
alpha <- 0.05
sims <- 200

pow.jb <- rep(NA, length(ssize))
pow.lillie <- rep(NA, length(ssize))
pow.ad <- rep(NA, length(ssize))
pow.shap <- rep(NA, length(ssize))

#Different sample sizes
for (i in 1:length(ssize)){
  prob.jb <- numeric(sims)
  prob.lillie <- numeric(sims)
  prob.ad <- numeric(sims)
  prob.shap <-numeric(sims)
  N <- ssize[i]

  #Run the simulations with the current sample size
  for (j in 1:sims){
    ##### Select one distribution ####
    #data <- rsn(N)
    data <- rchisq(N, df=5)
    #data <- rcauchy(N)
    #data <- rt(N,df=5)
    ##################################

    #Execute the normalitytests
    prob.jb[j] <- jb.norm.test(data)$p.value
    prob.lillie[j] <- lillie.test(data)$p.value
    prob.ad[j] <- ad.test(data)$p.value
    prob.shap[j] <- shapiro.test(data)$p.value
  }

  #Calculate the powers of the normalitytests
  pow.jb[i] <- mean(prob.jb < alpha)
  pow.lillie[i] <- mean(prob.lillie < alpha)
  pow.ad[i] <- mean(prob.ad < alpha)
  pow.shap[i] <- mean(prob.shap < alpha)
}

#Plot the results
df<- data.frame(samples=ssize,
                jb=pow.jb,
                lillie=pow.lillie,
                ad=pow.ad,
                shap=pow.shap)
ggplot(df, aes(samples)) +
  geom_line(aes(y = jb, colour = "jb")) +
  geom_line(aes(y = lillie, colour = "lillie")) +
  geom_line(aes(y = ad, colour = "ad")) +
  geom_line(aes(y = shap, colour = "shap"))