# Load Packages
library(openxlsx)
library(readxl)
library(performance)
library(lmtest)
library(sandwich)
library(car)

# Read DF
df2 <- read_excel("/Users/carolineseton/Desktop/df2.xlsx")

# Moderation Analyses Unadjusted Models
# Define variables
dependent_vars <- c("GAD", "GDS", "WEMWBS")
berq_subscales <- c("BERQ_ActivelyApproach", "BERQ_Distraction", "BERQ_Withdrawal",
                    "BERQ_SocialSupport", "BERQ_Ignoring")
cerq_subscales <- c("CERQ_SelfBlame", "CERQ_Acceptance", "CERQ_Rumination",
                    "CERQ_Refocus", "CERQ_Planning", "CERQ_Reappraisal",
                    "CERQ_Perspective", "CERQ_Catastrophising", "CERQ_OtherBlame")
predictors <- c("CFRS", "PCA_RC1", "PCA_RC2")

# Combine all ER subscales
emotion_regulation_strategies <- c(berq_subscales, cerq_subscales)

# Loop through all combinations
for (dependent in dependent_vars) {
  for (emotion_regulation in emotion_regulation_strategies) {
    for (predictor in predictors) {
      
      # Construct the formula (no covariates)
      formula <- as.formula(paste(dependent, "~", emotion_regulation, "*", predictor))
      
      # Fit the model
      outcome_model <- lm(formula, data = df2)
      
      # Build interaction term label
      interaction_term <- paste(emotion_regulation, ":", predictor, sep = "")
      
      # Check if interaction term exists in the model
      if (interaction_term %in% rownames(summary(outcome_model)$coefficients)) {
        interaction_pval <- summary(outcome_model)$coefficients[interaction_term, "Pr(>|t|)"]
        
        # If interaction is significant, print results
        if (interaction_pval < 0.05) {
          cat("\nSignificant Interaction Model:", emotion_regulation, "*", predictor, "→", dependent, "\n")
          print(summary(outcome_model))
          
          # Assumption testing using the performance package
          cat("\nAssumption Testing:\n")
          homoscedasticity_test <- check_heteroscedasticity(outcome_model)
          normality_test <- check_normality(outcome_model)
          outlier_test <- check_outliers(outcome_model)
          collinearity_test <- check_collinearity(outcome_model)
          linearity_test <- raintest(outcome_model)
          
          # Print assumption diagnostics
          print(homoscedasticity_test)
          print(normality_test)
          print(outlier_test)
          print(collinearity_test)
          print(linearity_test)
          
          # Shapiro-Wilk Normality Test (residual check)
          resids <- residuals(outcome_model)
          if (all(is.finite(resids))) {
            shapiro_test <- shapiro.test(resids)
            cat("\nShapiro-Wilk Test for Normality:\n")
            print(shapiro_test)
          }
          
          # Ramsey RESET for Linearity
          reset_test <- resettest(outcome_model)
          cat("\nRamsey RESET Test for Linearity:\n")
          print(reset_test)
          
          # Robust standard errors
          robust_se_outcome <- vcovHC(outcome_model, type = "HC1")
          cat("\nRobust Standard Errors (HC1):\n")
          print(coeftest(outcome_model, vcov = robust_se_outcome))
        }
      }
    }
  }
}