######################################################################
# study2_analysis.R
# Study 2 confirmatory + secondary analysis (preregistered specification)
# Manuscript: "Regaining Trust in AI: How Trust Erodes and Recovers
#              After Temporary and Permanent Forecasting Errors"
#
# Built from analysis_master.R (Study 1 OLS/ME/figures/Excel) and
# AI_Advise_Revision.R (participant-cluster bootstrap, BIC/AIC + LRT,
# breakpoint/dummy robustness). Specification follows the OSF
# preregistration and the 2026-09-04 transparent-changes note:
#   - Confirmatory sample: record_source = "connect", status = completed,
#     30/30 rounds, consented. Manipulation-check failures stay in the
#     primary sample and get a labeled sensitivity analysis.
#   - Phases: C1 = task rounds 19-25 (erosion), C2 = rounds >= 26 (recovery);
#     A-series uses a single post-shock dummy D1 = round >= 19.
#     Robustness: erosion from round 18.
#   - Direct replication subset: theta-matched AO vs LS-D4 (2 x 2).
#     A1-A3, B1-B3 with (a) OLS, (b) LME random intercept, (c) participant-
#     cluster bootstrap 95% CI (B = 1000). Plus the confirmatory
#     Advisor x Shock x Phase LME.
#   - LS extension: LME Advisor x Profile(D4/D8/D12) x Phase for WOA_adj,
#     trust_pre, trust_post, trust_change. Overall 8-cell model too.
#   - Two-sided tests, alpha = .05.
#
# USAGE (Rscript, run from C:/Temp/AIAdvisor):
#   Rscript study2_analysis.R extract connect     # pull DB -> locked CSVs
#   Rscript study2_analysis.R analyze connect     # analyze locked CSVs
#   Rscript study2_analysis.R extract simulation  # pipeline check on agents
#   Rscript study2_analysis.R analyze simulation
#   Optional 3rd arg: bootstrap replications (default 1000).
# Outputs -> study2_results/<source>/ : CSV snapshot, log, tables (xlsx),
#   figures (png).
######################################################################

.libPaths(c("C:/Users/meals/AppData/Local/R/win-library/4.6", .libPaths()))
options(warn = -1, width = 140)
suppressMessages({
  library(dplyr); library(tidyr); library(ggplot2)
  library(lme4);  library(lmerTest); library(openxlsx); library(jsonlite)
})

args   <- commandArgs(trailingOnly = TRUE)
MODE   <- if (length(args) >= 1) args[[1]] else "analyze"
SOURCE <- if (length(args) >= 2) args[[2]] else "connect"
B_BOOT <- if (length(args) >= 3) as.integer(args[[3]]) else 1000L
stopifnot(MODE %in% c("extract", "analyze"), SOURCE %in% c("connect", "simulation", "mturk"))

OUT_DIR <- file.path("study2_results", SOURCE)
dir.create(OUT_DIR, recursive = TRUE, showWarnings = FALSE)
F_PART  <- file.path(OUT_DIR, sprintf("study2_%s_participants.csv", SOURCE))
F_ROUND <- file.path(OUT_DIR, sprintf("study2_%s_rounds.csv", SOURCE))
F_FLOW  <- file.path(OUT_DIR, sprintf("study2_%s_participant_flow.csv", SOURCE))

######################################################################
# 0. EXTRACT: DB -> locked CSV snapshot (participant keys replaced by pid)
######################################################################
if (MODE == "extract") {
  library(DBI); library(RMariaDB)
  sec <- fromJSON("C:/Users/meals/.secret/ENVIRONMENT_VARIABLES.json")$aiforalab_server
  con <- dbConnect(RMariaDB::MariaDB(), host = sec$AIFORALAB_IP, port = as.integer(sec$AIFORALAB_SQL_PORT),
                   user = sec$AIFORALAB_DB_ID, password = sec$AIFORALAB_DB_PW, dbname = "aiadvice")
  on.exit(dbDisconnect(con), add = TRUE)

  # Participant-flow accounting (all statuses, the requested source)
  flow <- dbGetQuery(con, sprintf("
    SELECT advisor, shock, persistence, status, COUNT(*) AS n
    FROM study2_participant WHERE record_source = '%s'
    GROUP BY advisor, shock, persistence, status ORDER BY 1,2,3,4", SOURCE))
  write.csv(flow, F_FLOW, row.names = FALSE)

  part <- dbGetQuery(con, sprintf("
    SELECT p.participant_key, p.advisor, p.shock, p.persistence, p.status, p.consented,
           p.gender, p.age, p.statistics_familiarity, p.initial_trust, p.manipulation_check,
           p.agent_type, p.task_started_at, p.completed_at,
           (SELECT COUNT(*) FROM study2_round r WHERE r.participant_key = p.participant_key) AS n_rounds
    FROM study2_participant p
    WHERE p.record_source = '%s' AND p.status = 'completed'
    ORDER BY p.completed_at, p.participant_key", SOURCE))
  part$pid <- seq_len(nrow(part))
  rounds <- dbGetQuery(con, sprintf("
    SELECT r.participant_key, r.task_round, r.global_period, r.actual, r.advice,
           r.initial_forecast, r.final_forecast, r.stated_trust, r.trust_pre_error,
           r.trust_post_error, r.forecast_error, r.error_rate, r.raw_woa, r.woa_capped_0_2, r.woa_adj
    FROM study2_round r JOIN study2_participant p ON p.participant_key = r.participant_key
    WHERE p.record_source = '%s' AND p.status = 'completed'
    ORDER BY r.participant_key, r.task_round", SOURCE))
  rounds <- rounds %>% inner_join(part[, c("participant_key", "pid")], by = "participant_key")
  part$participant_key <- NULL; rounds$participant_key <- NULL
  write.csv(part,   F_PART,  row.names = FALSE)
  write.csv(rounds, F_ROUND, row.names = FALSE)
  cat(sprintf("Extracted %s: %d completed participants, %d round rows -> %s\n",
              SOURCE, nrow(part), nrow(rounds), OUT_DIR))
  cat("Snapshot time (KST):", format(Sys.time(), tz = "Asia/Seoul"), "\n")
  quit(save = "no")
}

######################################################################
# 1. LOAD LOCKED SNAPSHOT AND BUILD ANALYSIS VARIABLES
######################################################################
log_file <- file.path(OUT_DIR, sprintf("study2_%s_analysis_log.txt", SOURCE))
sink(log_file, split = TRUE)
cat("Study 2 analysis | source =", SOURCE, "| run", format(Sys.time(), tz = "Asia/Seoul"), "KST\n")
cat("Bootstrap replications:", B_BOOT, "\n\n")

part   <- read.csv(F_PART,  stringsAsFactors = FALSE)
rounds <- read.csv(F_ROUND, stringsAsFactors = FALSE)

# Inclusion: consented, completed, exactly 30 rounds
part <- part %>% filter(consented == 1, status == "completed", n_rounds == 30)
part <- part %>% mutate(
  cell    = ifelse(shock == "AO", paste0(advisor, "_AO"), paste0(advisor, "_LS_", persistence)),
  Profile = factor(ifelse(shock == "AO", "AO", persistence), levels = c("AO", "D4", "D8", "D12")),
  Advisor_Human = as.integer(advisor == "Human"),
  Shock_LS      = as.integer(shock == "LS"),
  mc_pass = manipulation_check == advisor)

df <- rounds %>% inner_join(part %>% select(pid, cell, advisor, shock, persistence, Profile,
                                           Advisor_Human, Shock_LS, mc_pass), by = "pid") %>%
  rename(Round = task_round) %>%
  arrange(pid, Round)

stopifnot(all(table(df$pid) == 30))
cat("Participants:", n_distinct(df$pid), " rows:", nrow(df), "\n")

# Outcome definitions (prereg): WOA_adj = max(0, 1 - |WOA - 1|); recheck DB value
df$WOA_adj_recalc <- pmax(0, 1 - abs(df$raw_woa - 1))
mism <- sum(abs(df$WOA_adj_recalc - df$woa_adj) > 1e-6, na.rm = TRUE)
cat("WOA_adj recomputation mismatches vs DB:", mism, "\n")
df$WOA_adj      <- df$woa_adj
df$WOA_raw      <- df$raw_woa
df$trust_pre    <- df$trust_pre_error
df$trust_post   <- df$trust_post_error
df$trust_change <- df$trust_post - df$trust_pre

# Phase coding (prereg): D1 = round >= 19 (A-series); C1 = 19..25, C2 = >= 26 (B-series)
df$D1 <- as.integer(df$Round >= 19)
df$C1 <- as.integer(df$Round >= 19 & df$Round <= 25)
df$C2 <- as.integer(df$Round >= 26)
df$D1_r18 <- as.integer(df$Round >= 18)
df$C1_r18 <- as.integer(df$Round >= 18 & df$Round <= 25)
df$phase <- factor(ifelse(df$Round < 19, "baseline", ifelse(df$Round <= 25, "erosion", "recovery")),
                   levels = c("baseline", "erosion", "recovery"))

######################################################################
# 2. HELPERS
######################################################################
stars_fn <- function(p) ifelse(p < .001, "***", ifelse(p < .01, "**", ifelse(p < .05, "*", "")))
fmt_p    <- function(p) ifelse(p < .001, "<.001", sprintf("%.3f", p))

coef_table <- function(m, label) {
  s <- summary(m)$coefficients
  data.frame(model = label, term = rownames(s), est = s[, 1], se = s[, 2],
             p = s[, ncol(s)], stars = stars_fn(s[, ncol(s)]), row.names = NULL)
}
fit_lme <- function(f, d) lmer(update(f, . ~ . + (1 | pid)), data = d, REML = FALSE)
icc_of  <- function(m) { vc <- as.data.frame(VarCorr(m)); vc$vcov[1] / sum(vc$vcov) }

# Participant-cluster bootstrap for a list of OLS formulas on data d
cluster_boot <- function(formulas, d, B = B_BOOT, seed = 20260904) {
  set.seed(seed)
  pids <- unique(d$pid); pid_data <- split(d, d$pid)
  res <- lapply(formulas, function(f) list())
  for (b in seq_len(B)) {
    samp <- sample(pids, length(pids), replace = TRUE)
    bd <- do.call(rbind, pid_data[as.character(samp)])
    bd$pid <- rep(seq_along(samp), each = 30)
    for (nm in names(formulas)) res[[nm]][[b]] <- coef(lm(formulas[[nm]], data = bd))
    if (b %% 250 == 0) cat(sprintf("  bootstrap %d/%d\n", b, B))
  }
  lapply(res, function(l) { m <- do.call(rbind, l)
    data.frame(term = colnames(m), boot_mean = colMeans(m), boot_se = apply(m, 2, sd),
               ci_lo = apply(m, 2, quantile, .025), ci_hi = apply(m, 2, quantile, .975), row.names = NULL) })
}

# Full A/B pipeline (OLS + LME + bootstrap + fit indices) for a 2 x 2 subset
run_ab_pipeline <- function(d, label, erosion_from = 19, do_boot = TRUE) {
  if (erosion_from == 18) { d$D1 <- d$D1_r18; d$C1 <- d$C1_r18 }
  forms <- list(
    A1 = WOA_adj ~ D1,
    A2 = WOA_adj ~ Advisor_Human + D1 + Advisor_Human:D1,
    A3 = WOA_adj ~ Shock_LS + D1 + Shock_LS:D1,
    B1 = WOA_adj ~ C1 + C2,
    B2 = WOA_adj ~ Advisor_Human + C1 + Advisor_Human:C1 + C2 + Advisor_Human:C2,
    B3 = WOA_adj ~ Shock_LS + C1 + Shock_LS:C1 + C2 + Shock_LS:C2)
  ols <- lapply(forms, lm, data = d)
  lme <- lapply(forms, fit_lme, d = d)
  cat(sprintf("\n===== %s (erosion from round %d) =====\n", label, erosion_from))
  tab <- bind_rows(
    lapply(names(forms), function(n) coef_table(ols[[n]], n) %>% mutate(estimator = "OLS")),
    lapply(names(forms), function(n) coef_table(lme[[n]], n) %>% mutate(estimator = "LME")))
  if (do_boot) {
    cat("Participant-cluster bootstrap...\n")
    bt <- cluster_boot(forms, d)
    bt_tab <- bind_rows(lapply(names(bt), function(n) bt[[n]] %>% mutate(model = n)))
    tab <- tab %>% left_join(bt_tab %>% select(model, term, ci_lo, ci_hi, boot_se), by = c("model", "term"))
  }
  fit <- data.frame(model = names(forms),
    k = sapply(ols, function(m) length(coef(m))),
    F = sapply(ols, function(m) summary(m)$fstatistic[1]),
    R2 = sapply(ols, function(m) summary(m)$r.squared),
    OLS_AIC = sapply(ols, AIC), OLS_BIC = sapply(ols, BIC),
    ME_logLik = sapply(lme, function(m) as.numeric(logLik(m))),
    ME_AIC = sapply(lme, AIC), ME_BIC = sapply(lme, BIC),
    ICC = sapply(lme, icc_of), row.names = NULL)
  lrt <- rbind(
    data.frame(test = "B1 -> B2 (advisor interactions)", as.data.frame(anova(lme$B1, lme$B2))[2, c("Chisq", "Df", "Pr(>Chisq)")]),
    data.frame(test = "B1 -> B3 (shock interactions)",   as.data.frame(anova(lme$B1, lme$B3))[2, c("Chisq", "Df", "Pr(>Chisq)")]))
  print(tab %>% filter(estimator == "OLS") %>% select(-estimator) %>% mutate(across(where(is.numeric), ~ round(., 4))))
  cat("\nLME (random intercept):\n")
  print(tab %>% filter(estimator == "LME") %>% select(model, term, est, se, p, stars) %>% mutate(across(where(is.numeric), ~ round(., 4))))
  cat("\nFit indices:\n"); print(fit %>% mutate(across(where(is.numeric), ~ round(., 3))))
  cat("\nNested LRT:\n"); print(lrt)
  list(coefs = tab, fit = fit, lrt = lrt, ols = ols, lme = lme)
}

# Phase-interaction LME: outcome ~ Advisor * Factor * (C1 + C2) + (1|pid)
run_phase_lme <- function(d, outcome, factor_var, label) {
  f <- as.formula(sprintf("%s ~ Advisor_Human * %s * (C1 + C2) + (1 | pid)", outcome, factor_var))
  m <- lmer(f, data = d, REML = FALSE)
  s <- summary(m)$coefficients
  tab <- data.frame(analysis = label, outcome = outcome, term = rownames(s), est = s[, 1], se = s[, 2],
                    df = s[, 3], t = s[, 4], p = s[, 5], stars = stars_fn(s[, 5]), row.names = NULL)
  cat(sprintf("\n--- %s | outcome = %s | ICC = %.3f ---\n", label, outcome, icc_of(m)))
  print(tab %>% select(term, est, se, p, stars) %>% mutate(across(where(is.numeric), ~ round(., 4))))
  list(model = m, table = tab)
}

######################################################################
# 3. PARTICIPANT FLOW, DEMOGRAPHICS, DESCRIPTIVES
######################################################################
cat("\n===== PARTICIPANT FLOW (all statuses, from DB extract) =====\n")
if (file.exists(F_FLOW)) print(read.csv(F_FLOW))

cat("\n===== CELL SIZES (analysis sample) =====\n")
cells <- part %>% count(cell) %>% arrange(cell); print(cells)
cat("Manipulation-check failures:", sum(!part$mc_pass), "of", nrow(part), "\n")

demo <- part %>% group_by(cell) %>% summarise(
  N = n(), Female_pct = round(100 * mean(gender == "Female", na.rm = TRUE), 1),
  M_Age = round(mean(age, na.rm = TRUE), 1), SD_Age = round(sd(age, na.rm = TRUE), 2),
  Stats_pct = round(100 * mean(statistics_familiarity == "Yes", na.rm = TRUE), 1),
  MC_fail = sum(!mc_pass), .groups = "drop")
overall_demo <- part %>% summarise(N = n(), Female_pct = round(100 * mean(gender == "Female", na.rm = TRUE), 1),
  M_Age = round(mean(age, na.rm = TRUE), 1), SD_Age = round(sd(age, na.rm = TRUE), 2),
  Stats_pct = round(100 * mean(statistics_familiarity == "Yes", na.rm = TRUE), 1), MC_fail = sum(!mc_pass))
cat("\nDemographics by cell:\n"); print(demo); cat("\nOverall:\n"); print(overall_demo)

phase_means <- df %>% group_by(cell, phase) %>% summarise(
  WOA_adj = mean(WOA_adj), WOA_raw = mean(WOA_raw), trust_pre = mean(trust_pre),
  trust_post = mean(trust_post), trust_change = mean(trust_change), .groups = "drop")
cat("\nPhase means by cell:\n"); print(phase_means %>% mutate(across(where(is.numeric), ~ round(., 3))), n = 50)

round_means <- df %>% group_by(cell, advisor, Profile, Round) %>%
  summarise(WOA_adj = mean(WOA_adj), trust_pre = mean(trust_pre), trust_post = mean(trust_post),
            trust_change = mean(trust_change), .groups = "drop")

# Advice-use classification (descriptive): reject (raw WOA <= .05), adopt (>= .95), compromise
df$advice_use <- ifelse(df$WOA_raw <= 0.05, "reject", ifelse(df$WOA_raw >= 0.95, "adopt", "compromise"))
use_tab <- df %>% count(cell, phase, advice_use) %>% group_by(cell, phase) %>%
  mutate(pct = round(100 * n / sum(n), 1)) %>% ungroup()
cat("\nAdvice-use classification (% by cell x phase):\n")
print(use_tab %>% select(-n) %>% pivot_wider(names_from = advice_use, values_from = pct), n = 50)

######################################################################
# 4. DIRECT REPLICATION: theta-matched AO vs LS-D4 (2 x 2)
######################################################################
rep_df <- df %>% filter(Profile %in% c("AO", "D4"))
cat(sprintf("\nReplication subset: %d participants (%s)\n", n_distinct(rep_df$pid),
            paste(sprintf("%s=%d", names(table(rep_df$cell[!duplicated(rep_df$pid)])),
                          table(rep_df$cell[!duplicated(rep_df$pid)])), collapse = ", ")))
REP <- run_ab_pipeline(rep_df, "DIRECT REPLICATION AO vs LS-D4", erosion_from = 19, do_boot = TRUE)

cat("\n===== CONFIRMATORY: Advisor x Shock x Phase LME (AO/LS-D4 subset) =====\n")
X3 <- run_phase_lme(rep_df, "WOA_adj", "Shock_LS", "Advisor x Shock x Phase")

# Linear contrasts used in Study 1 reporting: total recovery for Human (B2) and LS (B3)
contrast_total <- function(m, base, inter) {
  V <- vcov(m); b <- coef(m)
  est <- b[base] + b[inter]; se <- sqrt(V[base, base] + V[inter, inter] + 2 * V[base, inter])
  p <- 2 * pt(-abs(est / se), df.residual(m)); c(est = unname(est), se = unname(se), p = unname(p)) }
cat("\nB2 Human total recovery (gamma0 + gamma1):"); print(round(contrast_total(REP$ols$B2, "C2", "Advisor_Human:C2"), 4))
cat("B3 LS total recovery (gamma0 + gamma1):");     print(round(contrast_total(REP$ols$B3, "C2", "Shock_LS:C2"), 4))

######################################################################
# 5. ROBUSTNESS: erosion from round 18; raw WOA outcome; MC sensitivity
######################################################################
REP18 <- run_ab_pipeline(rep_df, "ROBUSTNESS erosion from round 18", erosion_from = 18, do_boot = FALSE)

cat("\n===== ROBUSTNESS: raw WOA (capped 0-2) as outcome, A3/B3 =====\n")
raw_df <- rep_df %>% mutate(WOA_adj = woa_capped_0_2)
RAW <- list(A3 = lm(WOA_adj ~ Shock_LS * D1, raw_df), B3 = lm(WOA_adj ~ Shock_LS * (C1 + C2), raw_df))
raw_tab <- bind_rows(coef_table(RAW$A3, "A3_rawWOA"), coef_table(RAW$B3, "B3_rawWOA")); print(raw_tab %>% mutate(across(where(is.numeric), ~ round(., 4))))

cat("\n===== SENSITIVITY: manipulation-check passers only =====\n")
mc_df <- rep_df %>% filter(mc_pass)
cat("Participants:", n_distinct(mc_df$pid), "\n")
MC <- run_ab_pipeline(mc_df, "MC-PASS ONLY AO vs LS-D4", erosion_from = 19, do_boot = FALSE)
MC3 <- run_phase_lme(mc_df, "WOA_adj", "Shock_LS", "Advisor x Shock x Phase (MC pass)")

######################################################################
# 6. LS EXTENSION: Advisor x Profile(D4/D8/D12) x Phase, four outcomes
######################################################################
ls_df <- df %>% filter(Profile %in% c("D4", "D8", "D12")) %>% mutate(Profile = droplevels(Profile))
cat(sprintf("\nLS extension subset: %d participants\n", n_distinct(ls_df$pid)))
EXT <- lapply(c("WOA_adj", "trust_pre", "trust_post", "trust_change"), function(o)
  run_phase_lme(ls_df, o, "Profile", "LS extension Advisor x Profile x Phase"))
names(EXT) <- c("WOA_adj", "trust_pre", "trust_post", "trust_change")

# Secondary trust outcomes in the replication subset (prereg: Advisor x Profile x Phase)
TRUST_REP <- lapply(c("trust_pre", "trust_post", "trust_change"), function(o)
  run_phase_lme(rep_df, o, "Shock_LS", "Secondary trust, AO/LS-D4 Advisor x Shock x Phase"))
names(TRUST_REP) <- c("trust_pre", "trust_post", "trust_change")

# Overall 8-cell model (does not create AO-D8/AO-D12 cells: Profile has 4 levels, AO ref)
ALL8 <- run_phase_lme(df, "WOA_adj", "Profile", "Overall 8-cell Advisor x Profile x Phase")

######################################################################
# 7. FIGURES
######################################################################
theme_set(theme_minimal())
vl <- list(geom_vline(xintercept = 17, linetype = "dashed", color = "red"),
           geom_vline(xintercept = 19, linetype = "dotted", color = "grey40"),
           geom_vline(xintercept = 26, linetype = "dotted", color = "blue"))
p1 <- ggplot(round_means, aes(Round, WOA_adj, color = Profile, linetype = advisor)) + geom_line(linewidth = .8) + vl +
  labs(title = "Study 2: mean WOA_adj by round, all 8 cells", y = "WOA_adj", linetype = "Advisor") + theme(legend.position = "bottom")
ggsave(file.path(OUT_DIR, "fig_woa_8cells.png"), p1, width = 9, height = 5, dpi = 150)
p2 <- ggplot(round_means %>% filter(Profile %in% c("AO", "D4")), aes(Round, WOA_adj, color = Profile, linetype = advisor)) +
  geom_line(linewidth = .9) + geom_point(size = 1.4) + vl +
  labs(title = "Theta-matched AO vs LS-D4 by advisor", y = "WOA_adj", linetype = "Advisor") + theme(legend.position = "bottom")
ggsave(file.path(OUT_DIR, "fig_woa_replication.png"), p2, width = 9, height = 5, dpi = 150)
p3 <- ggplot(round_means %>% filter(Profile != "AO") %>% group_by(Profile, Round) %>% summarise(WOA_adj = mean(WOA_adj), .groups = "drop"),
             aes(Round, WOA_adj, color = Profile)) + geom_line(linewidth = .9) + geom_point(size = 1.4) + vl +
  labs(title = "LS profiles D4/D8/D12 (collapsed across advisor)", y = "WOA_adj")
ggsave(file.path(OUT_DIR, "fig_woa_ls_profiles.png"), p3, width = 8, height = 4.5, dpi = 150)
tr_long <- round_means %>% select(Profile, advisor, Round, trust_pre, trust_post, trust_change) %>%
  pivot_longer(c(trust_pre, trust_post, trust_change), names_to = "measure")
p4 <- ggplot(tr_long %>% group_by(Profile, Round, measure) %>% summarise(value = mean(value), .groups = "drop"),
             aes(Round, value, color = Profile)) + geom_line(linewidth = .8) + vl +
  facet_wrap(~ measure, scales = "free_y", ncol = 1) + labs(title = "Stated trust by round and profile", y = NULL)
ggsave(file.path(OUT_DIR, "fig_trust_profiles.png"), p4, width = 8, height = 8, dpi = 150)
cat("\nFigures saved to", OUT_DIR, "\n")

######################################################################
# 8. EXCEL OUTPUT
######################################################################
wb <- createWorkbook()
add <- function(name, x) { addWorksheet(wb, name); writeData(wb, name, as.data.frame(x)) }
add("Cells", cells); add("Demographics", bind_rows(demo, overall_demo %>% mutate(cell = "ALL")))
if (file.exists(F_FLOW)) add("Participant_flow", read.csv(F_FLOW))
add("Phase_means", phase_means); add("Round_means", round_means); add("Advice_use", use_tab)
add("Rep_AB_coefs", REP$coefs); add("Rep_fit", REP$fit); add("Rep_LRT", REP$lrt)
add("Rep_3way_LME", X3$table)
add("Robust_r18_coefs", REP18$coefs); add("Robust_r18_fit", REP18$fit); add("Robust_rawWOA", raw_tab)
add("MC_pass_AB", MC$coefs); add("MC_pass_3way", MC3$table)
add("LS_ext_LME", bind_rows(lapply(EXT, function(e) e$table)))
add("Trust_rep_LME", bind_rows(lapply(TRUST_REP, function(e) e$table)))
add("All8_LME", ALL8$table)
saveWorkbook(wb, file.path(OUT_DIR, sprintf("study2_%s_results.xlsx", SOURCE)), overwrite = TRUE)
cat("Saved:", file.path(OUT_DIR, sprintf("study2_%s_results.xlsx", SOURCE)), "\n")
cat("\nDone", format(Sys.time(), tz = "Asia/Seoul"), "KST\n")
sink()