# SIMULATION FOR ERP N2-P2 COMPLEX AMPLITUDE

# Loading libraries
library("Superpower")
library("ggplot2")
library("esc")

# Initialize the main variables
# Based on Colloca et al., 2009 (in Tables 1, Fig.3B)

# ----
# DESIGN PARAMETERS
# Taken from Colloca, 2009 [absolute value of the difference between Baseline-PostTreatment]:
MeanUnCond = 22.54-18.10 # Mean of unconditioned expectations group 
MeanCond = 25.58-18.15 #  Mean of conditioned expectations group 
nG1 = 16 # N of unconditioned expectations group 
nG2 = 16 # N of conditioned expectations group 
AbsoluteES = (MeanCond-MeanUnCond) * 0.5 # half of the absolute effect size
MeanCond = AbsoluteES + MeanUnCond  # Recalculate the MeanCond value 

# Taken from Ribordy Lambert et al., 2020
SDUnCond = 2.6 # SD of unconditioned expectations group 
SDCond = 2.6 # SD of conditioned expectations group  


# SIMULATIONS PARAMETERS
Increment = 2 # Factor to increment parameter N (sample size)
N1 = 50 # Minimum N
N2 = 80 # Maximum N
Nseq = seq(N1,N2,Increment) # Sequence of N to test
Alpha = 0.05 # alpha level 
Power = NULL # Empty vector to store power values for each simulated N
counter = 0 # Initialize loop counter
Seed = 1234 # Seed for reproducibility
Nsim = 10000 # Number of simulations


# Effect size (hedge's gs + correction of bias for small sample sizes)
ESC = esc_mean_sd(grp1m = MeanUnCond, grp1sd = SDUnCond, grp1n = nG1,
            grp2m = MeanCond, grp2sd = SDCond, grp2n = nG2, es.type = "g")
HedgesG = abs(ESC$es)

# ----
# COMPUTING POWER ANALYSES [For each N]

for (k in 1:length(Nseq)){
  
  # Build design
  Design = ANOVA_design(design = "2b", n = Nseq[k], mu = c(MeanUnCond, MeanCond), sd = c(SDUnCond, SDCond), 
                        labelnames=c("Groups","Unconditioned expectations","Conditioned expectations"))
  
  # Build simulated datasets
  Datasets = ANOVA_power(Design, alpha_level = Alpha, correction = "none", nsims = Nsim, seed = Seed, verbose = F) 
  
  # Store results for each iteration
  # pc_results = Power analysis results for pairwise comparisons (in this case it is the same as in main_results, expect for effect size)
  Power[k] = Datasets$pc_results$power
  
  # Counter 
  counter = counter + 1
  Msg = sprintf("Computation status: %g%%", round((counter/length(Nseq)*100), 2))
  print(Msg)
}

# ----
# FIGURE
# Subtitle for the figure
Subtitle = sprintf("MeanUnCond = %.1f, MeanCond = %.1f, SDUnCond = %.1f, SDCond = %.1f, Alpha = %.2f, g = %.2f",MeanUnCond,MeanCond,SDUnCond,SDCond,Alpha,HedgesG)
Table = cbind.data.frame(Power,Nseq) 
g = ggplot(Table, aes(x = Nseq, y = Power)) + 
  geom_point() + geom_line() + 
  scale_x_continuous(limits = c(N1,N2), breaks = Nseq) +
  labs(title = "Sample size estimation",subtitle = Subtitle, x = "Total sample size", y = "Power") + 
  theme_bw() +
  geom_hline(yintercept=90, linetype="dashed", color = "red")

print(g)

# Save plot
# ggsave("EffectSize.png", device = "png", dpi = 300)
ggsave("Fig1EffectSize.tiff", device = "tiff", dpi = 300)