library(SuppDists)
library(dplyr)
library(tidyr)

setwd("/home/ladmin/Documents/DND/EOHI/eohi3")

between_vars <- c("perspective", "temporalDO")
within_vars_MEAN <- c(
  "past_pref_MEAN", "past_pers_MEAN", "past_val_MEAN",
  "fut_pref_MEAN",  "fut_pers_MEAN",  "fut_val_MEAN"
)
within_vars_DGEN <- c(
  "past_pref_DGEN", "past_pers_DGEN", "past_val_DGEN",
  "fut_pref_DGEN",  "fut_pers_DGEN",  "fut_val_DGEN"
)

df <- read.csv("eohi3.csv", stringsAsFactors = FALSE, check.names = FALSE, na.strings = "NA")

anova_data_MEAN <- df %>%
  select(pID, all_of(between_vars), all_of(within_vars_MEAN)) %>%
  filter(!is.na(perspective), perspective != "",
         !is.na(temporalDO), temporalDO != "")

long_data_MEAN <- anova_data_MEAN %>%
  pivot_longer(
    cols = all_of(within_vars_MEAN),
    names_to = "variable",
    values_to = "MEAN_SCORE"
  ) %>%
  mutate(
    time   = ifelse(grepl("^past_", variable), "past", "fut"),
    domain = case_when(
      grepl("_pref_MEAN$", variable) ~ "pref",
      grepl("_pers_MEAN$", variable) ~ "pers",
      grepl("_val_MEAN$",  variable) ~ "val",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    TIME       = factor(time,   levels = c("past", "fut")),
    DOMAIN     = factor(domain, levels = c("pref", "pers", "val")),
    perspective = factor(perspective),
    temporalDO  = factor(temporalDO)
  ) %>%
  select(pID, perspective, temporalDO, TIME, DOMAIN, MEAN_SCORE) %>%
  filter(!is.na(MEAN_SCORE))

cell_vars_MEAN <- long_data_MEAN %>%
  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
  summarise(
    n        = n(),
    variance = var(MEAN_SCORE, na.rm = TRUE),
    .groups  = "drop"
  )

fmax_by_cell_MEAN <- cell_vars_MEAN %>%
  group_by(TIME, DOMAIN) %>%
  summarise(
    Fmax_observed = max(variance, na.rm = TRUE) / min(variance, na.rm = TRUE),
    df_min        = min(n) - 1L,
    .groups       = "drop"
  )

k <- 4

fmax_table_MEAN <- fmax_by_cell_MEAN %>%
  rowwise() %>%
  mutate(
    alpha_0.05 = SuppDists::qmaxFratio(0.95, df = df_min, k = k),
    alpha_0.01 = SuppDists::qmaxFratio(0.99, df = df_min, k = k)
  ) %>%
  ungroup() %>%
  mutate(
    Fmax_observed = round(Fmax_observed, 4),
    alpha_0.05    = round(alpha_0.05, 4),
    alpha_0.01    = round(alpha_0.01, 4)
  ) %>%
  select(TIME, DOMAIN, Fmax_observed, alpha_0.05, alpha_0.01)

# ---- MEAN: Print observed Hartley ratios ----
cat("\n--- Hartley ratios (MEAN) ---\n")
fmax_table_MEAN %>%
  mutate(across(where(is.numeric), ~ format(round(., 4), nsmall = 4))) %>%
  print()

# ---- DGEN: Observed Hartley ratios ----
anova_data_DGEN <- df %>%
  select(pID, all_of(between_vars), all_of(within_vars_DGEN)) %>%
  filter(!is.na(perspective), perspective != "",
         !is.na(temporalDO), temporalDO != "")

long_data_DGEN <- anova_data_DGEN %>%
  pivot_longer(
    cols = all_of(within_vars_DGEN),
    names_to = "variable",
    values_to = "DGEN_SCORE"
  ) %>%
  mutate(
    time   = ifelse(grepl("^past_", variable), "past", "fut"),
    domain = case_when(
      grepl("_pref_DGEN$", variable) ~ "pref",
      grepl("_pers_DGEN$", variable) ~ "pers",
      grepl("_val_DGEN$",  variable) ~ "val",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    TIME       = factor(time,   levels = c("past", "fut")),
    DOMAIN     = factor(domain, levels = c("pref", "pers", "val")),
    perspective = factor(perspective),
    temporalDO  = factor(temporalDO)
  ) %>%
  select(pID, perspective, temporalDO, TIME, DOMAIN, DGEN_SCORE) %>%
  filter(!is.na(DGEN_SCORE))

cell_vars_DGEN <- long_data_DGEN %>%
  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
  summarise(
    n        = n(),
    variance = var(DGEN_SCORE, na.rm = TRUE),
    .groups  = "drop"
  )

fmax_by_cell_DGEN <- cell_vars_DGEN %>%
  group_by(TIME, DOMAIN) %>%
  summarise(
    Fmax_observed = max(variance, na.rm = TRUE) / min(variance, na.rm = TRUE),
    df_min        = min(n) - 1L,
    .groups       = "drop"
  )

fmax_table_DGEN <- fmax_by_cell_DGEN %>%
  rowwise() %>%
  mutate(
    alpha_0.05 = SuppDists::qmaxFratio(0.95, df = df_min, k = k),
    alpha_0.01 = SuppDists::qmaxFratio(0.99, df = df_min, k = k)
  ) %>%
  ungroup() %>%
  mutate(
    Fmax_observed = round(Fmax_observed, 4),
    alpha_0.05    = round(alpha_0.05, 4),
    alpha_0.01    = round(alpha_0.01, 4)
  ) %>%
  select(TIME, DOMAIN, Fmax_observed, alpha_0.05, alpha_0.01)

cat("\n--- Hartley ratios (DGEN) ---\n")
fmax_table_DGEN %>%
  mutate(across(where(is.numeric), ~ format(round(., 4), nsmall = 4))) %>%
  print()