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

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

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

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

# 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 
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)

#Create MODEL1 
MODEL1 <- lmer(DV ~ FA_pitch + FA_formant + FA_risetime + FA_attention + FA_integration + Test + (1|ID), data = data)

# Print the summary of the model
summary(MODEL1)

# Print R squared
performance::r2(MODEL1)

#Create MODEL1_interaction to check any interaction effects 

MODEL1_interaction <- lmer(DV ~ FA_pitch*Test + FA_formant*Test + FA_risetime*Test + FA_attention*Test + FA_integration*Test + (1|ID), data = data)

# Print the summary of the model
summary(MODEL1_interaction)

# Print R squared
performance::r2(MODEL1_interaction)


#Create MODEL2
# First, develop a linear model with control variables as predictors
control_model_experience <- lm(DV ~ AOA + LOR + L2useAVE + LengthofEFL + MusicTrainingYN, data = data)

# Create a new dependent variable by subtracting the predicted values from the original dependent variable
data$DV_residuals_experience <- resid(control_model_experience)

# Create a new mixed effects model using the new dependent variable
MODEL2 <- lmer(DV_residuals_experience ~ FA_pitch + FA_formant + FA_risetime + FA_attention + FA_integration + Test + (1|ID), data = data)


# Print the summary of the model
summary(MODEL2)
# Print R squared
performance::r2(MODEL2)

#Create MODEL3
# First, develop a linear model with control variables as predictors
control_model_wm <- lm(DV ~ digit_span_forwards + digit_span_backwards, data = data)

# Create a new dependent variable by subtracting the predicted values from the original dependent variable
data$DV_residuals_wm <- resid(control_model_wm)

# Create a new mixed effects model using the new dependent variable
MODEL3 <- lmer(DV_residuals_wm ~ FA_pitch + FA_formant + FA_risetime + FA_attention + FA_integration + Test + (1|ID), data = data)

# Print the summary of the model
summary(MODEL3)
# Print R squared
performance::r2(MODEL3)

# Compare MODEL0 vs. MODEL1
anova(MODEL0, MODEL1)

# Compare MODEL2 vs. MODEL0

# For the model comparison, MODEL0 was revised using residuals, with experience factors controlled for.

MODEL0_experience <- lmer(DV_residuals_experience ~ FA_pitch*Test + FA_formant*Test + FA_risetime*Test + (1|ID), data = data)

anova(MODEL2, MODEL0_experience)

# Compare MODEL3 with MODEL0

# For the model comparison, MODEL0 was revised using residuals, with working memory factors controlled for.

MODEL0_wm <- lmer(DV_residuals_wm ~ FA_pitch*Test + FA_formant*Test + FA_risetime*Test + (1|ID), data = data)

anova(MODEL3, MODEL0_wm)