# 02_PostHocTests_MediumDataset.R
# PeerJ-CS, "Personality-based pair programming"
#
# This script:
#   1) (Re)fits 'model_interaction' & 'model_cluster' if needed.
#   2) Performs post-hoc tests on Role-based & Cluster-based models.
#   3) Computes Cohen's d for pairwise Role differences in the interaction model.
#   4) Outputs results to .txt files and optionally plots.
# Note: the warning “Results may be misleading due to involvement in interactions”
# indicates that we are requesting the marginal means (and effect sizes) for a 
# main effect of Role, even though Role interacts with the Big Five variables in
# our model. That’s not a fatal error, just a caution from emmeans.


# ------------------------------------------------
# 1) Install/load required libraries
# ------------------------------------------------
required_packages <- c("openxlsx", "nlme", "dplyr", "emmeans", "ggplot2")
new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
lapply(required_packages, require, character.only = TRUE)

# ------------------------------------------------
# 2) (Re)Fit Models (or load from previous script)
# ------------------------------------------------
model_interaction <- lme(
  fixed = IntrinsicMotivation ~ Role * (B5_O + B5_C + B5_E + B5_A + B5_N),
  random = ~1 | Student_ID,
  data = long_data,
  method = "REML"
)

model_cluster <- lme(
  fixed = IntrinsicMotivation ~ Role * PersonalityCluster,
  random = ~1 | Student_ID,
  data = long_data,
  method = "REML"
)

# ------------------------------------------------
# 3) Role-based pairwise comparisons (Interaction Model)
# ------------------------------------------------
emmeans_interaction_role <- emmeans(model_interaction, pairwise ~ Role, adjust = "tukey")
role_summary <- summary(emmeans_interaction_role, infer = TRUE)   # Means + pairwise

# Save the raw post-hoc results (with estimate, SE, t-ratio, p-value) to file
capture.output(role_summary, file = "PostHoc_Interaction_Role.txt", append = FALSE)

# ------------------------------------------------
# 4) Cohen's d for each pairwise difference
# ------------------------------------------------
# A) Emmeans of Role (main effect) -> pairwise contrasts
role_emm <- emmeans(model_interaction, ~ Role)
role_contr <- pairs(role_emm, adjust = "tukey")  # returns an emmGrid of differences

# B) Supply sigma & edf from our model
sigma_value <- model_interaction$sigma # residual SD
edf_value   <- 1014

# C) Compute effect sizes
role_effsize <- eff_size(role_emm, sigma = sigma_value, edf = edf_value, contrast = "pairs")

# D) Summarize effect sizes with inference & confidence intervals
role_effsize_summary <- summary(role_effsize, infer = TRUE)
role_effsize_confint <- confint(role_effsize)
role_effsize_summary
role_effsize_confint

# E) Export them to a text file
capture.output(
  list(
    "CohenD_Pairwise_Summary" = role_effsize_summary,
    "CohenD_Pairwise_CI"      = role_effsize_confint
  ),
  file = "PostHoc_Interaction_Effsize.txt",
  append = FALSE
)

# ------------------------------------------------
# 5) Cluster model post-hoc tests
# ------------------------------------------------
emmeans_cluster_role <- emmeans(model_cluster, pairwise ~ Role, adjust = "tukey")
emmeans_cluster_by_cluster <- emmeans(model_cluster, pairwise ~ Role | PersonalityCluster, adjust = "tukey")
#emmeans_cluster_role
#emmeans_cluster_by_cluster

capture.output(
  summary(emmeans_cluster_role),
  file = "PostHoc_Cluster_Role.txt"
)
capture.output(
  summary(emmeans_cluster_by_cluster),
  file = "PostHoc_Cluster_byCluster.txt"
)

# Optionally, plots
plot(emmeans_interaction_role$emmeans)
plot(emmeans_cluster_role$emmeans)

cat("\nPost-hoc tests complete. Results (including effect sizes) saved to text files,
    e.g., PostHoc_Interaction_Role.txt, PostHoc_Interaction_Effsize.txt, and others.\n")