#### data preparation ##
data <- read.csv2 ("data_music.csv",na="NA", dec=",")
head(data)
describe(data)

# check skewness
library(moments)  # For skewness and kurtosis
skewness(data$global_psqi, na.rm = TRUE)
skewness(data$sleep_habit_numeric, na.rm = TRUE)
skewness(data$stress, na.rm = TRUE)
skewness(data$mood, na.rm = TRUE)
skewness(data$anxiety, na.rm = TRUE)

# scale variables
data$global_psqi_z <- data$global_psqi %>% scale()
data$sleep_habit_numeric_z <- data$sleep_habit_numeric %>% scale()
data$anxiety_z <- data$anxiety %>% scale()
data$mood_z <- data$mood %>% scale()
data$stress_z <- data$stress %>% scale()
data$spae_01_z <- data$spae_01 %>% scale()

# correlation matrix 
library(Hmisc)
data_corr <- subset(data, select = c(global_psqi, sleep_habit_numeric, anxiety, mood, stress, spae_01))

cor_results <- rcorr(as.matrix(data_corr))  # Compute correlations and p-values
cor_matrix <- round(cor_results$r, 2)  # Extract and round correlation coefficients
p_values <- round(cor_results$P, 3)  # Extract and round p-values
View(cor_matrix)  # Correlation coefficients
View(p_values)  # Corresponding p-values
significant_cor <- cor_matrix
significant_cor[p_values >= 0.05] <- NA  # Set non-significant correlations to NA
View(significant_cor)


#### research question 1 ####
# Create interaction terms manually
data$sleep_habit_anxiety_interaction <- data$sleep_habit_numeric_z * data$anxiety_z
data$sleep_habit_mood_interaction <- data$sleep_habit_numeric_z * data$mood_z
data$sleep_habit_stress_interaction <- data$sleep_habit_numeric_z * data$stress_z
data$sleep_habit_fitness_interaction <- data$sleep_habit_numeric_z * data$spae_01_z

# define the model
moderation_model_RQ1 <- '
  global_psqi_z ~ sleep_habit_numeric_z + anxiety_z + mood_z + stress_z + spae_01_z +
                sleep_habit_anxiety_interaction + sleep_habit_mood_interaction + sleep_habit_stress_interaction + 
                sleep_habit_fitness_interaction

  # Correlations
  anxiety_z ~~ mood_z
  anxiety_z ~~ stress_z
  mood_z ~~ stress_z
'

moderation_model_RQ1_fit <- sem(data = data, model = moderation_model_RQ1, estimator = "ML", missing = "FIML", fixed.x = FALSE)
varTable(moderation_model_RQ1_fit)
parameterEstimates(moderation_model_RQ1_fit)
summary(moderation_model_RQ1_fit, fit.measures = TRUE, rsq = TRUE, standardized = TRUE)

lavaanPlot(model = moderation_model_RQ1_fit, node_options = list(shape = "box", fontname = "Helvetica"), 
           edge_options = list(color = "grey"), coefs = TRUE, sig = 0.05)


### simple slopes analysis ###
library(interactions)  
library(ggplot2)
library(emmeans)

# Build a linear model with the interaction
model <- lm(global_psqi_z ~ sleep_habit_numeric_z * anxiety_z + mood_z + stress_z + spae_01_z, data = data)

# Get simple slopes
emtrends(model, ~ sleep_habit_numeric_z, var = "anxiety_z", at = list(sleep_habit_numeric_z = c(-1, 0, 1)))


### create figure ### 
# Create predicted values
library(interactions)
plot <- interact_plot(model, pred = anxiety_z, modx = sleep_habit_numeric_z,
              modx.values = c(-1, 0, 1),
              legend.main = "Listening to Pre-Sleep Music",
              modx.labels = c("-1 SD", "+/-0 SD", "+1 SD"),
              plot.points = FALSE,
              interval = TRUE) +
  labs(x = "Anxiety", y = "Poorer Sleep Quality", color = "Music Listening") +
  theme_classic()

# Save as a PNG
ggsave("interaction_plot.png", plot, width = 8, height = 6, dpi = 300)



### research question 2 ### 
# create dummy variable for gender
data <- data %>%
  mutate(
    gender_female = ifelse(gender == "1", 1, 0),
    gender_diverse = ifelse(gender == "3", 1, 0)
    # No need to create gender_male; it's the reference category
  )

data$age_z <- data$age %>% scale()
data$spae_01_z <- data$spae_01 %>% scale()

# Create interaction terms manually
data$sleep_habit_age_interaction <- data$sleep_habit_numeric_z * data$age_z
data$sleep_habit_gender_interaction <- data$sleep_habit_numeric_z * data$gender_female
data$sleep_habit_spae_interaction <- data$sleep_habit_numeric_z * data$spae_01_z

# define the model
moderation_model_RQ2 <- '
  global_psqi_z ~ sleep_habit_numeric_z + age_z + gender_female + spae_01_z + 
                sleep_habit_age_interaction + sleep_habit_gender_interaction + sleep_habit_spae_interaction
'
moderation_model__RQ2_fit <- sem(data = data, model = moderation_model_RQ2, estimator = "ML", missing = "FIML", fixed.x = FALSE)
parameterEstimates(moderation_model__RQ2_fit)
summary(moderation_model__RQ2_fit, fit.measures = TRUE, rsq = TRUE, standardized = TRUE)

lavaanPlot(model = moderation_model_RQ2, node_options = list(shape = "box", fontname = "Helvetica"), 
           edge_options = list(color = "grey"), coefs = TRUE, sig = 0.05)