# Load necessary packages
if (!require("lme4")) install.packages("lme4")
if (!require("simr")) install.packages("simr")
if (!require("readxl")) install.packages("readxl")
if (!require("performance")) install.packages("performance")

library(lme4)
library(simr)
library(readxl)
library(performance)

# Set your file path
file_path <- "data_in_brief_submitted.xlsx"

# Import the data
data <- read_excel(file_path, sheet = "SIMR_Spanish")

# Convert 'ID' to a factor if it's not already
data$ID <- as.factor(data$ID)

# Convert 'Test' to a factor if it's not already
data$Test <- as.factor(data$Test)

# Relevel 'Test' with 'a' as the reference level
data$Test <- relevel(data$Test, ref = "a")

#Create MODEL0 (including pitch, formant, and risetime acuity) 
MODEL0 <- lmer(DV ~ FA_pitch*Test + FA_formant*Test + FA_risetime*Test + (1|ID), data = data)

# Print the summary of the model
summary(MODEL0)

# Print R squared
performance::r2(MODEL0)

#POWER ANALYSIS
power_sim_pitch <- powerSim(MODEL0, test = fixed("FA_pitch", "z"), nsim = 100)

print(power_sim_pitch)

# Extend the number of levels for 'ID' to n = 200 (*** new n has to be more than simulated dataset)
MODEL0_ID <- extend(MODEL0, along="ID", n=200)  

power_sim_pitch_ID <- powerSim(MODEL0_ID, test = fixed("FA_pitch", "z"), nsim = 10)
print(power_sim_pitch_ID)

# Define power curve for 'ID' for each predictor *if* N = 200 

pc_ID_pitch <- powerCurve(MODEL0_ID, along = "ID", test = fixed("FA_pitch", "z"), nsim = 100)

#Print power curve for 'ID' for each predictor
print(pc_ID_pitch)
plot(pc_ID_pitch)