# MIXED-EFFECTS REGRESSION

# Load necessary libraries
library(ggplot2)
library(dplyr)
library(lme4)
library(emmeans)
library(car)

# Prompt the user to select a file (LongHFScores.csv <- available in the Supplementary Materials)
selected_file <- file.choose()

# Load the dataframe file (LongHFScores.csv <- available in the Supplementary Materials)
data <- read.csv(selected_file)


# Ensure the data is correctly loaded
data$CONDITION <- factor(data$CONDITION, levels = c(1, 2, 3), labels = c("ADHD", "ASD", "ADHD&ASD"))
data$DIAGNOSIS <- factor(data$DIAGNOSIS, levels = c(0, 1), labels = c("No Diagnosis", "Formal Diagnosis"))
data$THERAPY <- factor(data$THERAPY, levels = c(0, 1), labels = c("No Therapy", "Therapy"))
data$SOCIAL_ISOLATION <- factor(data$SOCIAL_ISOLATION, levels = c(0, 1), labels = c("Before", "During"))

# Fit a mixed-effects model with main effects and interactions
model_main_effects <- lmer(TOTAL_SCORE ~ SOCIAL_ISOLATION * THERAPY + SOCIAL_ISOLATION * DIAGNOSIS + CONDITION +
                           (1 | ID), data = data, REML = TRUE)

# Display model summary
summary(model_main_effects)

# Perform Type III fixed effects tests
Anova(model_main_effects, type = "III")

# Calculate and display estimated marginal means (EMMs) for main effects and interactions
emm_main_effects <- emmeans(model_main_effects, ~ SOCIAL_ISOLATION * THERAPY * DIAGNOSIS + CONDITION)
pairs(emm_main_effects, adjust = "bonferroni")

# Check multicollinearity for main effects only, excluding interaction terms
vif_model_main_no_interaction <- lm(TOTAL_SCORE ~ SOCIAL_ISOLATION + THERAPY + DIAGNOSIS + CONDITION, data = data)
vif(vif_model_main_no_interaction)



# Ensure the factors are correctly labelled for plotting
data$CONDITION <- factor(data$CONDITION, levels = c("ADHD", "ASD", "ADHD&ASD"))
data$SOCIAL_ISOLATION <- factor(data$SOCIAL_ISOLATION, levels = c("Before", "During"))
data$THERAPY <- factor(data$THERAPY, levels = c("No Therapy", "Therapy"))
data$DIAGNOSIS <- factor(data$DIAGNOSIS, levels = c("No Diagnosis", "Formal Diagnosis"))

# Fit the mixed-effects model
model <- lmer(TOTAL_SCORE ~ SOCIAL_ISOLATION * CONDITION + SOCIAL_ISOLATION * DIAGNOSIS + SOCIAL_ISOLATION * THERAPY + 
              (1 | ID), data = data, REML = TRUE)

# Obtain Estimated Marginal Means (EMMs) for each interaction

# 1. EMMs for SOCIAL_ISOLATION by CONDITION
emm_condition <- emmeans(model, ~ SOCIAL_ISOLATION | CONDITION)

# 2. EMMs for SOCIAL_ISOLATION by DIAGNOSIS
emm_diagnosis <- emmeans(model, ~ SOCIAL_ISOLATION | DIAGNOSIS)

# 3. EMMs for SOCIAL_ISOLATION by THERAPY
emm_therapy <- emmeans(model, ~ SOCIAL_ISOLATION | THERAPY)

# Convert EMM results to data frames for plotting
df_condition <- as.data.frame(emm_condition)
df_diagnosis <- as.data.frame(emm_diagnosis)
df_therapy <- as.data.frame(emm_therapy)

# Define position dodge for slight horizontal shift
dodge_position <- position_dodge(width = 0.3)

# Plot 1: Estimated Mean HF Total Score by CONDITION and SOCIAL_ISOLATION
plot_condition <- ggplot(df_condition, aes(x = SOCIAL_ISOLATION, y = emmean, color = CONDITION, group = CONDITION)) +
  geom_point(size = 3, position = dodge_position) +
  geom_line(position = dodge_position) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = dodge_position) +
  labs(title = "Estimated Mean HF Total Score by Condition and Social Isolation",
       x = "Social Isolation", y = "Estimated Mean HF Total Score") +
  theme_minimal() +
  ylim(64, 75)  # Set y-axis limits

# Display the first plot
print(plot_condition)

# Plot 2: Estimated Mean HF Total Score by DIAGNOSIS and SOCIAL_ISOLATION
plot_diagnosis <- ggplot(df_diagnosis, aes(x = SOCIAL_ISOLATION, y = emmean, color = DIAGNOSIS, group = DIAGNOSIS)) +
  geom_point(size = 3, position = dodge_position) +
  geom_line(position = dodge_position) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = dodge_position) +
  labs(title = "Estimated Mean HF Total Score by Diagnosis and Social Isolation",
       x = "Social Isolation", y = "Estimated Mean HF Total Score") +
  theme_minimal() +
  ylim(64, 75)  # Set y-axis limits

# Display the second plot
print(plot_diagnosis)

# Plot 3: Estimated Mean HF Total Score by THERAPY and SOCIAL_ISOLATION
plot_therapy <- ggplot(df_therapy, aes(x = SOCIAL_ISOLATION, y = emmean, color = THERAPY, group = THERAPY)) +
  geom_point(size = 3, position = dodge_position) +
  geom_line(position = dodge_position) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = dodge_position) +
  labs(title = "Estimated Mean HF Total Score by Therapy and Social Isolation",
       x = "Social Isolation", y = "Estimated Mean HF Total Score") +
  theme_minimal() +
  ylim(64, 75)  # Set y-axis limits

# Display the third plot
print(plot_therapy)