# Power Analysis Script
#
# Power analysis for the WIBM BCH auxiliary analyses.
# Maps onto one-way ANOVA with k = 4 groups; effect size = Cohen's f.
#
# TWO APPROACHES (following Quirk, 2020):
#   1. ANALYTICAL  — pwr::pwr.anova.test(): closed-form, exact under normal/
#      equal-n assumptions. Used for a priori power, sensitivity, power curve.
#   2. SIMULATION  — Monte Carlo via replicate(), 5,000 reps. Validates the
#      analytical result and extends it to actual unequal class n's.
#
# REFERENCES
# Cohen, J. (1988). Statistical power analysis for the behavioral sciences
#   (2nd ed.). Lawrence Erlbaum Associates.
# Quirk, C. (2020). A case for simulated power analyses.
#   https://colinquirk.com/simulated-power/
#
# 1. SETUP ---------------------------------------------------------------------
if (!requireNamespace("pwr",      quietly = TRUE)) install.packages("pwr")
if (!requireNamespace("ggplot2",  quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("dplyr",    quietly = TRUE)) install.packages("dplyr")
if (!requireNamespace("tidyr",    quietly = TRUE)) install.packages("tidyr")
if (!requireNamespace("scales",   quietly = TRUE)) install.packages("scales")
if (!requireNamespace("flextable",quietly = TRUE)) install.packages("flextable")
if (!requireNamespace("officer",  quietly = TRUE)) install.packages("officer")

library(pwr)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(flextable)
library(officer)

set.seed(123)

cat("pwr version  :", as.character(packageVersion("pwr")), "\n")
cat("Analysis date:", format(Sys.Date(), "%d %B %Y"), "\n\n")


# 2. OUTPUT PATHS --------------------------------------------------------------

base_project <- here::here("output", "power_analysis")

figures_power_folder <- file.path(base_project, "Figures")
tables_power_folder  <- file.path(base_project, "Tables")

dir.create(figures_power_folder, showWarnings = FALSE, recursive = TRUE)
dir.create(tables_power_folder,  showWarnings = FALSE, recursive = TRUE)

cat(sprintf("Figures -> %s\n", figures_power_folder))
cat(sprintf("Tables  -> %s\n\n", tables_power_folder))

fig_path   <- function(...) file.path(figures_power_folder, ...)
table_path <- function(...) file.path(tables_power_folder,  ...)

save_fig <- function(p, stem, width = 8, height = 5) {
  ggsave(fig_path(paste0(stem, ".png")),
         plot = p, width = width, height = height,
         dpi = 600, bg = "white")
  ggsave(fig_path(paste0(stem, ".tiff")),
         plot = p, width = width, height = height,
         dpi = 600, bg = "white", compression = "lzw")
  cat(sprintf("  [SAVED] %s.png / .tiff\n", stem))
  invisible(p)
}

save_tbl_docx <- function(df, stem, title = NULL, note = NULL) {
  ft <- flextable(df) |>
    theme_booktabs() |>
    autofit()
  if (!is.null(title)) ft <- set_caption(ft, caption = title)
  doc <- read_docx() |> body_add_flextable(ft)
  if (!is.null(note)) {
    doc <- doc |> body_add_par(paste0("Note. ", note), style = "Normal")
  }
  path <- table_path(paste0(stem, ".docx"))
  print(doc, target = path)
  cat(sprintf("  [SAVED] %s.docx\n", stem))
  invisible(path)
}


# 3. STUDY PARAMETERS ----------------------------------------------------------

k        <- 4
N_total  <- 1441
alpha    <- 0.05
target   <- 0.80
n_sim    <- 5000

n_IO <- 421; n_IB <- 478; n_AR <- 279; n_NR <- 263
class_ns <- c(n_IO, n_IB, n_AR, n_NR)

# Harmonic mean n/group: recommended approximation for unequal groups (Cohen, 1988, p.361)
n_harmonic <- k / sum(1 / class_ns)

# Population SD of group index vector — used to scale means to hit target f exactly.
# Cohen's f uses population SD (divide by N), not R's sample SD (divide by N-1).
pop_sd_index <- sd(0:(k - 1)) * sqrt((k - 1) / k)

cat("=== STUDY PARAMETERS ===\n")
cat(sprintf("  Groups (k)            : %d\n",   k))
cat(sprintf("  Total N               : %d\n",   N_total))
cat(sprintf("  Class sizes           : IO=%d  IB=%d  AR=%d  NR=%d\n",
            n_IO, n_IB, n_AR, n_NR))
cat(sprintf("  Harmonic mean n/group : %.1f\n", n_harmonic))
cat(sprintf("  Alpha                 : %.2f\n", alpha))
cat(sprintf("  Target power          : %.2f\n", target))
cat(sprintf("  Simulations (n_sim)   : %d\n\n", n_sim))


# 4. HELPER FUNCTIONS ----------------------------------------------------------

# Cohen's f from BCH-weighted group means and sigma_within
cohen_f_from_means <- function(means, ns, sigma_within) {
  stopifnot(length(means) == length(ns))
  n_tot    <- sum(ns)
  mu_grand <- sum(ns * means) / n_tot
  sqrt(sum(ns * (means - mu_grand)^2) / n_tot) / sigma_within |> round(4)
}

# Construct k evenly-spaced means whose population SD = f * sigma_within exactly
make_means_for_f <- function(f, sigma_within = 1.0, k = 4) {
  index <- 0:(k - 1)
  index * ((f * sigma_within) / pop_sd_index)
}

# Analytical power
analytical_power <- function(f, n = n_harmonic, grps = k, sig = alpha) {
  pwr.anova.test(k = grps, n = n, f = f, sig.level = sig)$power
}

# Simulation helpers (Quirk, 2020)
get_data <- function(means, sds, ns) {
  data.frame(
    y     = unlist(mapply(rnorm, n = ns, mean = means, sd = sds,
                          SIMPLIFY = FALSE)),
    group = factor(rep(seq_along(ns), times = ns))
  )
}

run_anova <- function(means, sds, ns, sig = alpha) {
  d   <- get_data(means, sds, ns)
  fit <- aov(y ~ group, data = d)
  summary(fit)[[1]][["Pr(>F)"]][1] < sig
}

sim_power <- function(means, sds, ns, nsim = n_sim, sig = alpha) {
  mean(replicate(nsim, run_anova(means, sds, ns, sig)))
}

# Run one outcome through both analytical and simulation methods
power_for_outcome <- function(var_label, means, ns, sigma_within, grps = 4,
                               nsim = n_sim) {
  f_obs   <- cohen_f_from_means(means, ns, sigma_within)
  n_hmean <- grps / sum(1 / ns)
  pwr_an  <- pwr.anova.test(k = grps, n = n_hmean,
                              f = f_obs, sig.level = alpha)$power
  pwr_sim <- sim_power(means, rep(sigma_within, grps), ns, nsim = nsim)
  tibble(
    variable       = var_label,
    k_groups       = grps,
    f              = f_obs,
    n_harmonic     = round(n_hmean, 1),
    power_analytic = round(pwr_an,  4),
    power_sim      = round(pwr_sim, 4)
  )
}

# 5. STEP 1 — VALIDATE: SIMULATION VS ANALYTICAL -------------------------------
# Validation on Cohen (1988) Ex. 8.1: k=4, n=20, f=0.28 (expected power ≈ 0.515).
# Confirms simulation matches analytical before applying to study data.

cat("=== STEP 1: VALIDATE SIMULATION vs ANALYTICAL ===\n")
cat("Cohen (1988) Ex. 8.1: k=4, n=20, f=0.28  (expected power ≈ 0.515)\n\n")

val_f     <- 0.28
val_n     <- 20
val_sigma <- 1.0
val_means <- make_means_for_f(val_f, val_sigma, k)
val_ns    <- rep(val_n, k)

val_f_check <- cohen_f_from_means(val_means, val_ns, val_sigma)
cat(sprintf("  Target f = %.4f  |  Constructed f = %.4f  (should match)\n\n",
            val_f, val_f_check))

analytic_val  <- pwr.anova.test(k = k, n = val_n,
                                 f = val_f, sig.level = alpha)$power
simulated_val <- sim_power(val_means, rep(val_sigma, k), val_ns)

cat(sprintf("  Analytical power    : %.4f\n", analytic_val))
cat(sprintf("  Simulated power     : %.4f\n", simulated_val))
cat(sprintf("  Absolute difference : %.4f\n", abs(analytic_val - simulated_val)))

tol <- 0.03
if (abs(analytic_val - simulated_val) <= tol) {
  cat(sprintf("  Result: within tolerance (<=%.0f pp) — VALIDATED \u2713\n\n", tol * 100))
} else {
  cat(sprintf("  Result: EXCEEDS tolerance (%.0f pp) — check n_sim or mean construction\n\n",
              tol * 100))
}


# 6. STEP 2 — A PRIORI POWER AT N = 1,441 --------------------------------------

cat("=== STEP 2: A PRIORI POWER AT N = 1,441 ===\n")
cat(sprintf("  Harmonic mean n/group = %.1f;  k = %d;  alpha = %.2f\n\n",
            n_harmonic, k, alpha))

f_benchmarks <- c(small = 0.10, medium = 0.25, large = 0.40)

apriori <- tibble(
  label          = names(f_benchmarks),
  f              = f_benchmarks,
  power_analytic = sapply(f_benchmarks, analytical_power),
  power_sim      = sapply(f_benchmarks, function(f) {
    m <- make_means_for_f(f, 1.0, k)
    sim_power(m, rep(1, k), rep(round(n_harmonic), k))
  })
)

print(apriori, n = Inf)
cat("\n")

save_tbl_docx(
  apriori |>
    mutate(across(where(is.numeric), ~ round(., 4))) |>
    rename(`Effect size` = label, `Cohen's f` = f,
           `Power (analytical)` = power_analytic,
           `Power (simulated)`  = power_sim),
  stem  = "PA_Table1_Apriori_Power",
  title = "Table 1\nA Priori Power at N = 1,441 by Effect Size Benchmark",
  note  = paste0(
    "k = 4 groups, harmonic mean n/group = ", round(n_harmonic, 1),
    ", alpha = ", alpha, ". Analytical = pwr::pwr.anova.test(). ",
    "Simulated = Monte Carlo, n_sim = ", n_sim, ", set.seed(123). ",
    "Both power columns rounded to 4 decimal places. ",
    "Benchmarks from Cohen (1988).")
)


# 7. STEP 3 — SENSITIVITY ANALYSIS ---------------------------------------------

cat("=== STEP 3: SENSITIVITY — MINIMUM DETECTABLE f AT 80% POWER ===\n")

sens     <- pwr.anova.test(k = k, n = n_harmonic,
                            sig.level = alpha, power = target)
f_min    <- sens$f
eta2_min <- f_min^2 / (1 + f_min^2)

cat(sprintf("  Minimum detectable Cohen's f : %.4f\n", f_min))
cat(sprintf("  Equivalent eta-squared (\u03b7\u00b2)  : %.4f\n", eta2_min))
cat("  (Benchmarks: small \u03b7\u00b2 = 0.01, medium = 0.06, large = 0.14)\n\n")

cat("=== Required total N for 80% power ===\n")
f_seq <- c(0.10, 0.12, 0.15, 0.18, 0.20, 0.25, 0.30, 0.40)
req_n <- tibble(
  `Cohen's f`   = f_seq,
  `n per group` = sapply(f_seq, function(f)
    ceiling(pwr.anova.test(k = k, f = f,
                            sig.level = alpha, power = target)$n)),
  `Total N`     = `n per group` * k
)
print(req_n, n = Inf)
cat("\n")

sens_tbl <- tibble(
  Parameter = c("Groups (k)", "Harmonic mean n/group", "Alpha",
                "Target power", "Min detectable f", "Equivalent \u03b7\u00b2"),
  Value     = c(as.character(k),
                as.character(round(n_harmonic, 1)),
                as.character(alpha),
                as.character(target),
                sprintf("%.3f", round(f_min, 3)),
                sprintf("%.3f", round(eta2_min, 3)))
)

save_tbl_docx(
  sens_tbl,
  stem  = "PA_Table2_Sensitivity",
  title = "Table 2\nSensitivity Analysis: Minimum Detectable Effect at 80% Power",
  note  = paste0(
    "Analytical solution via pwr::pwr.anova.test(). N = ", N_total,
    " total participants. Benchmarks: small \u03b7\u00b2 = 0.01, medium = 0.06, ",
    "large = 0.14 (Cohen, 1988).")
)

save_tbl_docx(
  req_n,
  stem  = "PA_Table3_Required_N",
  title = "Table 3\nRequired Total N for 80% Power Across Effect Sizes (Future Study Planning)",
  note  = paste0(
    "k = 4 groups, alpha = ", alpha,
    ". n per group rounded up to next integer. Total N = n per group x k.")
)


# 8. STEP 4 — POWER CURVE FIGURE -----------------------------------------------

n_range  <- seq(10, 500, by = 5)

curve_df <- expand.grid(
  n_per_grp = n_range,
  f         = c(0.10, 0.25, 0.40)
) |>
  as_tibble() |>
  mutate(
    power   = mapply(function(n, f)
      pwr.anova.test(k = k, n = n, f = f, sig.level = alpha)$power,
      n_per_grp, f),
    f_label = factor(paste0("f = ", f),
                     levels = c("f = 0.1", "f = 0.25", "f = 0.4"))
  )

fig_curve <- ggplot(curve_df,
                    aes(x = n_per_grp, y = power,
                        colour = f_label, group = f_label)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 0.80, linetype = "dashed",
             colour = "black", linewidth = 0.6) +
  geom_vline(xintercept = round(n_harmonic), linetype = "dotted",
             colour = "grey40", linewidth = 0.6) +
  annotate("text",
           x = round(n_harmonic) + 8, y = 0.04,
           label = paste0("Harmonic mean\nn/group = ", round(n_harmonic)),
           hjust = 0, size = 3, colour = "grey30") +
  annotate("text",
           x = 12, y = 0.82,
           label = "80% power threshold",
           hjust = 0, size = 3, colour = "black") +
  scale_colour_manual(
    values = c("f = 0.1"  = "#CC6677",
               "f = 0.25" = "#0D2B55",
               "f = 0.4"  = "#44AA99"),
    name = "Effect size") +
  scale_x_continuous(breaks = seq(0, 500, 50)) +
  scale_y_continuous(labels = percent_format(accuracy = 1),
                     limits = c(0, 1.01)) +
  labs(
    x       = "N per group",
    y       = "Power",
    caption = paste0(
      "One-way ANOVA, k = 4 groups, \u03b1 = 0.05. ",
      "Vertical dotted line = harmonic mean class size (n = ",
      round(n_harmonic), "). ",
      "Analytical solution via pwr::pwr.anova.test().")
  ) +
  theme_classic(base_size = 11) +
  theme(legend.position = "right",
        plot.caption    = element_text(size = 8, colour = "grey40", hjust = 0))

print(fig_curve)
save_fig(fig_curve, "fig_power_curve", width = 8, height = 5)


# 9. STEP 5 — OBSERVED POWER PER BCH OUTCOME -----------------------------------
# Cohen's f computed from BCH-weighted class means and actual class n's.
# sigma_within = scale_range / 4 = 1.0 for all 1-5 ordinal items (Cohen, 1988).
# Age uses observed total-sample SD (16.3 years).
# 3-group outcomes exclude No RMT class (structural NAs for RMT-specific items).

cat("=== STEP 5: OBSERVED POWER PER BCH OUTCOME ===\n")
cat(sprintf("  (Simulation: %d replications per outcome)\n\n", n_sim))

# 4-group outcomes (all four classes included)
outcomes_4grp <- list(

  # Demographics
  list(var   = "Age (years)",
       means = c(38.73, 34.65, 36.11, 41.58),
       ns    = c(380,   327,   470,   263),
       sigma = 16.3),

  list(var   = "Playing ability (0.5-5.0)",
       means = c(3.95, 4.22, 4.16, 3.81),
       ns    = c(380,  327,  469,  263),
       sigma = 1.0),

  list(var   = "Playing frequency (1-5)",
       means = c(4.11, 4.28, 4.17, 3.63),
       ns    = c(380,  327,  470,  263),
       sigma = 1.0),

  list(var   = "Years playing (1-5)",
       means = c(3.61, 3.63, 3.60, 3.92),
       ns    = c(380,  327,  470,  263),
       sigma = 1.0),

  # Breathlessness
  list(var   = "Breathing awareness (1-5)",
       means = c(3.58, 3.85, 3.67, 3.19),
       ns    = c(380,  326,  470,  262),
       sigma = 1.0),

  list(var   = "Worst symptom frequency / freqMAX (1-5)",
       means = c(2.55, 2.65, 2.78, 2.63),
       ns    = c(331,  298,  432,  223),
       sigma = 1.0),

  # Beliefs — opinions
  list(var   = "Playing enough (1-5)",
       means = c(3.92, 3.55, 3.59, 3.62),
       ns    = c(366,  322,  463,  245),
       sigma = 1.0),

  list(var   = "RMT improves performance (1-5)",
       means = c(4.13, 4.41, 4.31, 3.84),
       ns    = c(298,  294,  425,  174),
       sigma = 1.0),

  # Beliefs — body importance
  list(var   = "Oral/facial muscles (1-5)",
       means = c(4.34, 4.45, 4.35, 4.32),
       ns    = c(378,  326,  462,  260),
       sigma = 1.0),

  list(var   = "Airways & lungs (1-5)",
       means = c(4.54, 4.70, 4.62, 4.41),
       ns    = c(378,  325,  457,  260),
       sigma = 1.0),

  list(var   = "Breathing muscles general (1-5)",
       means = c(4.40, 4.64, 4.51, 4.18),
       ns    = c(374,  324,  459,  249),
       sigma = 1.0),

  list(var   = "Postural muscles (1-5)",
       means = c(4.00, 4.20, 4.14, 3.81),
       ns    = c(370,  317,  445,  242),
       sigma = 1.0),

  list(var   = "Diaphragm (1-5)",
       means = c(4.26, 4.25, 4.15, 4.40),
       ns    = c(368,  320,  452,  245),
       sigma = 1.0),

  list(var   = "Abdominal muscles (1-5)",
       means = c(3.95, 4.05, 3.95, 3.72),
       ns    = c(359,  317,  444,  227),
       sigma = 1.0),

  list(var   = "Rib cage muscles (1-5)",
       means = c(3.79, 3.97, 3.96, 3.57),
       ns    = c(352,  308,  442,  226),
       sigma = 1.0),

  list(var   = "Accessory muscles (1-5)",
       means = c(3.55, 3.61, 3.63, 3.24),
       ns    = c(237,  227,  334,  155),
       sigma = 1.0),

  # Receptivity
  list(var   = "Athlete identity pre-survey (1-5)",
       means = c(3.93, 4.08, 4.02, 3.63),
       ns    = c(375,  323,  465,  254),
       sigma = 1.0),

  list(var   = "Athlete identity post-survey (1-5)",
       means = c(3.81, 4.01, 3.95, 3.54),
       ns    = c(380,  327,  470,  263),
       sigma = 1.0),

  list(var   = "Interest with evidence / RCT arm (1-5)",
       means = c(3.14, 3.30, 3.40, 2.88),
       ns    = c(201,  155,  229,  136),
       sigma = 1.0),

  list(var   = "Interest without evidence / non-RCT arm (1-5)",
       means = c(2.98, 3.34, 3.34, 2.75),
       ns    = c(180,  172,  241,  127),
       sigma = 1.0)
)

# 3-group outcomes (No RMT excluded — structural NAs for RMT-specific items)
outcomes_3grp <- list(

  list(var   = "RMT frequency (1-5) [IO, IB, AR only]",
       means = c(3.70, 4.26, 4.06),
       ns    = c(379,  327,  469),
       sigma = 1.0),

  list(var   = "RMT Eff: Instrument (1-5) [IO, IB, AR only]",
       means = c(3.66, 3.96, 3.79),
       ns    = c(369,  324,  456),
       sigma = 1.0),

  list(var   = "RMT Eff: Body (1-5) [IO, IB, AR only]",
       means = c(3.11, 3.87, 3.71),
       ns    = c(253,  320,  431),
       sigma = 1.0),

  list(var   = "RMT Eff: Device (1-5) [IO, IB, AR only]",
       means = c(2.87, 3.02, 3.51),
       ns    = c(172,  119,  279),
       sigma = 1.0)
)

cat("Running 4-group outcomes...\n")
results_4grp <- bind_rows(lapply(outcomes_4grp, function(x)
  power_for_outcome(x$var, x$means, x$ns, x$sigma, grps = 4)))

cat("Running 3-group outcomes (No RMT excluded)...\n")
results_3grp <- bind_rows(lapply(outcomes_3grp, function(x)
  power_for_outcome(x$var, x$means, x$ns, x$sigma, grps = 3)))


# 9b. BINARY OUTCOMES (proportions) -------------------------------------------
# For binary variables, sigma_within = sqrt(p_pooled * (1 - p_pooled))
# where p_pooled is the BCH-weighted grand proportion across all included groups.
# This treats the binary outcome as continuous under the ANOVA framework (Cohen, 1988).
# Proportions taken directly from BCH descriptive tables.
# Note: Unsure response dummies are not included; the ANOVA power framework applies
# to primary outcomes. Binary sparse variables (ensemble, self-taught, sparse regions,
# Asia/Africa/LatAm, individual disorder subtypes with n<5%) are excluded.

sigma_binary <- function(props, ns) {
  p_pool <- sum(props * ns) / sum(ns)
  sqrt(p_pool * (1 - p_pool))
}

outcomes_binary_4grp <- list(

  # ── Demographics ─────────────────────────────────────────────────────────────
  list(var   = "Gender: Female (binary)",
       props = c(0.48, 0.51, 0.52, 0.47),
       ns    = c(255,  358,  305,  445)),

  list(var   = "Role: Leisure (binary)",
       props = c(0.63, 0.52, 0.44, 0.43),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Role: Performer (binary)",
       props = c(0.59, 0.62, 0.65, 0.64),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Role: Student (binary)",
       props = c(0.25, 0.36, 0.42, 0.38),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Role: Teacher (binary)",
       props = c(0.17, 0.28, 0.43, 0.42),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Education: Tertiary (binary)",
       props = c(0.34, 0.42, 0.55, 0.52),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Education: Non-tertiary (binary)",
       props = c(0.54, 0.48, 0.38, 0.40),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Instrument: Woodwind (binary)",
       props = c(0.32, 0.36, 0.39, 0.36),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Instrument: Brass (binary)",
       props = c(0.40, 0.41, 0.48, 0.49),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Instrument: Reed (binary)",
       props = c(0.50, 0.53, 0.44, 0.48),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Region of residence: Europe (binary)",
       props = c(0.46, 0.35, 0.22, 0.23),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Region of residence: N. America (binary)",
       props = c(0.32, 0.37, 0.56, 0.51),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Region of residence: Oceania (binary)",
       props = c(0.22, 0.25, 0.20, 0.23),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Region of education: Europe (binary)",
       props = c(0.47, 0.36, 0.23, 0.23),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Region of education: N. America (binary)",
       props = c(0.32, 0.38, 0.57, 0.52),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Region of education: Oceania (binary)",
       props = c(0.21, 0.24, 0.19, 0.23),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Any diagnosed disorder (binary)",
       props = c(0.45, 0.45, 0.52, 0.52),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Disorder: Mental health (binary; disorder subsample)",
       props = c(0.55, 0.60, 0.73, 0.67),
       ns    = c(127,  165,  164,  245)),

  list(var   = "Disorder: Respiratory (binary; disorder subsample)",
       props = c(0.27, 0.40, 0.30, 0.31),
       ns    = c(127,  165,  164,  245)),

  list(var   = "Disorder: Other (binary; disorder subsample)",
       props = c(0.40, 0.34, 0.34, 0.41),
       ns    = c(127,  165,  164,  245)),

  list(var   = "Performance income (binary; performer subsample)",
       props = c(0.17, 0.15, 0.20, 0.28),
       ns    = c(154,  230,  203,  283)),

  list(var   = "Teaching income (binary; teacher subsample)",
       props = c(0.50, 0.58, 0.59, 0.65),
       ns    = c(44,   105,  135,  190)),

  # ── Beliefs: influences ───────────────────────────────────────────────────────
  list(var   = "Influence: Music teacher (binary)",
       props = c(0.56, 0.71, 0.84, 0.69),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Influence: WI peers (binary)",
       props = c(0.20, 0.39, 0.53, 0.48),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Influence: Personal research (binary)",
       props = c(0.19, 0.38, 0.51, 0.50),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Influence: Music school (binary)",
       props = c(0.16, 0.29, 0.42, 0.36),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Influence: Non-musical education (binary)",
       props = c(0.12, 0.17, 0.27, 0.30),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Influence: Medical practitioner (binary)",
       props = c(0.08, 0.10, 0.18, 0.21),
       ns    = c(263,  380,  327,  470)),

  # ── Beliefs: perceived barriers (non-device users only) ──────────────────────
  list(var   = "Barrier: Accessibility issues (binary; non-device)",
       props = c(0.14, 0.25, 0.36, 0.22),
       ns    = c(263,  380,  327,   75)),

  list(var   = "Barrier: Education/knowledge (binary; non-device)",
       props = c(0.62, 0.56, 0.50, 0.32),
       ns    = c(263,  380,  327,   75)),

  list(var   = "Barrier: Better alternatives (binary; non-device)",
       props = c(0.36, 0.41, 0.44, 0.30),
       ns    = c(263,  380,  327,   75)),

  list(var   = "Barrier: Method issues (binary; non-device)",
       props = c(0.13, 0.09, 0.10, 0.07),
       ns    = c(263,  380,  327,   75)),

  # ── Dyspnoea prevalence ───────────────────────────────────────────────────────
  list(var   = "Dyspnoea: Breathlessness (binary)",
       props = c(0.40, 0.38, 0.45, 0.42),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Breathing discomfort (binary)",
       props = c(0.17, 0.20, 0.31, 0.30),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Breathing effort (binary)",
       props = c(0.21, 0.28, 0.37, 0.33),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Air hunger (binary)",
       props = c(0.33, 0.39, 0.48, 0.44),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Chest tightness (binary)",
       props = c(0.14, 0.17, 0.30, 0.32),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Mental breath effort (binary)",
       props = c(0.14, 0.20, 0.33, 0.32),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Unplanned breaths (binary)",
       props = c(0.39, 0.42, 0.48, 0.44),
       ns    = c(263,  380,  327,  470)),

  list(var   = "Dyspnoea: Unfinished phrases (binary)",
       props = c(0.59, 0.65, 0.70, 0.62),
       ns    = c(263,  380,  327,  470))
)

# 3-group binary outcomes (No RMT excluded — structural NAs or inapplicable)
outcomes_binary_3grp <- list(

  list(var   = "Teaches students: Instrument (binary; teacher subsample)",
       props = c(0.62, 0.98, 0.99, 0.98),
       ns    = c(45,   108,  140,  197)),

  list(var   = "Teaches students: Body (binary; teacher subsample)",
       props = c(0.38, 0.46, 0.91, 0.82),
       ns    = c(45,   108,  140,  197)),

  list(var   = "Teaches students: None (binary; teacher subsample)",
       props = c(0.36, 0.02, 0.01, 0.02),
       ns    = c(45,   108,  140,  197))
)

cat("Running 4-group binary outcomes...\n")
results_binary_4grp <- bind_rows(lapply(outcomes_binary_4grp, function(x) {
  sig <- sigma_binary(x$props, x$ns)
  power_for_outcome(x$var, x$props, x$ns, sig, grps = 4)
}))

cat("Running 3-group binary outcomes...\n")
results_binary_3grp <- bind_rows(lapply(outcomes_binary_3grp, function(x) {
  sig <- sigma_binary(x$props, x$ns)
  power_for_outcome(x$var, x$props, x$ns, sig, grps = 3)
}))

results_all <- bind_rows(results_4grp, results_3grp,
                         results_binary_4grp, results_binary_3grp) |>
  arrange(desc(f))

cat("\nResults (sorted by observed f, largest first):\n\n")
print(results_all, n = Inf)

cat("\n=== SUMMARY ===\n")
cat(sprintf("  Total outcomes assessed       : %d\n",    nrow(results_all)))
cat(sprintf("  Power \u2265 80%% (analytical)    : %d / %d\n",
            sum(results_all$power_analytic >= 0.80), nrow(results_all)))
cat(sprintf("  Power \u2265 80%% (simulated)     : %d / %d\n",
            sum(results_all$power_sim      >= 0.80), nrow(results_all)))
cat(sprintf("  Minimum detectable f (80%%)    : %.4f\n",  f_min))
cat(sprintf("  Smallest observed f            : %.4f (%s)\n",
            min(results_all$f),
            results_all$variable[which.min(results_all$f)]))
cat(sprintf("  Largest observed f             : %.4f (%s)\n\n",
            max(results_all$f),
            results_all$variable[which.max(results_all$f)]))

save_tbl_docx(
  results_all |>
    mutate(f = round(f, 4)) |>
    rename(
      `Outcome variable`   = variable,
      `Groups`             = k_groups,
      `Cohen's f`          = f,
      `n (harmonic mean)`  = n_harmonic,
      `Power (analytical)` = power_analytic,
      `Power (simulated)`  = power_sim
    ),
  stem  = "PA_Table4_Observed_Power_Per_Outcome",
  title = "Table 4\nObserved Cohen's f and Power Per BCH Outcome Variable",
  note  = paste0(
    "Cohen's f computed from BCH-weighted class means or proportions and n_eff values. ",
    "For ordinal outcomes: sigma_within = scale range / 4 = 1.0 for 1-5 items (Cohen, 1988); ",
    "age uses observed SD = 16.3 years. ",
    "For binary outcomes: sigma_within = sqrt(p_pooled x (1 - p_pooled)) where p_pooled is the ",
    "BCH-weighted grand proportion across included classes (Cohen, 1988 binary ANOVA framework). ",
    "Power (analytical) = pwr::pwr.anova.test() using harmonic mean n/group. ",
    "Power (simulated) = Monte Carlo with actual unequal class n's ",
    "(n_sim = ", n_sim, ", set.seed(123); Quirk, 2020). ",
    "Variables labelled [IO, IB, AR only] exclude No RMT class (structural NAs). ",
    "Unsure response dummies not included: the ANOVA framework applies to primary outcomes. ",
    "Excluded variables: ensemble (n = 31, 2.1% of sample; too sparse for stable BCH estimates) and ",
    "self-taught (captured within the ed_nonTertiary composite in the primary BCH analysis); ",
    "Asia, Africa, and Latin America regions excluded due to near-zero cell counts (< 1% each). ",
    "Outcomes sorted by Cohen's f, largest first.")
)


# 10. FIGURE: OBSERVED f AND SIMULATED POWER PER OUTCOME -----------------------

plot_df <- results_all |>
  mutate(
    variable  = factor(variable, levels = rev(variable)),
    powered   = power_sim >= 0.80,
    pwr_label = paste0(round(power_sim * 100, 0), "%")
  )

fig_obs <- ggplot(plot_df, aes(x = f, y = variable, colour = powered)) +
  geom_vline(xintercept = f_min, linetype = "dashed",
             colour = "grey55", linewidth = 0.6) +
  geom_point(size = 3.2) +
  geom_text(aes(label = pwr_label), hjust = -0.3, size = 2.5) +
  scale_colour_manual(
    values = c(`TRUE`  = "#0D2B55",
               `FALSE` = "#CC6677"),
    labels = c(`TRUE`  = "\u2265 80% power",
               `FALSE` = "< 80% power"),
    name = NULL) +
  scale_x_continuous(
    limits = c(0, max(plot_df$f) * 1.4),
    breaks = seq(0, 1, 0.05)) +
  annotate("text",
           x = f_min + 0.003, y = 0.7,
           label = paste0("Min detectable f\n= ", round(f_min, 3)),
           hjust = 0, vjust = 0, size = 2.9, colour = "grey40") +
  labs(
    x       = "Cohen's f  (computed from BCH-weighted class means or proportions)",
    y       = NULL,
    caption = paste0(
      "Power estimated via Monte Carlo simulation (", format(n_sim, big.mark = ","),
      " replications per outcome; Quirk, 2020).\n",
      "Labels = simulated power (%). ",
      "Dashed line = minimum detectable f at 80% power with N = ",
      format(N_total, big.mark = ","), ", k = ", k,
      " groups, \u03b1 = ", alpha, ".\n",
      "Ordinal outcomes: sigma_within = scale range / 4 = 1.0 (Cohen, 1988); age uses SD = 16.3 years.\n",
      "Binary outcomes: sigma_within = sqrt(p_pooled \u00d7 (1 - p_pooled)).")
  ) +
  theme_classic(base_size = 10) +
  theme(
    legend.position = "top",
    axis.text.y     = element_text(size = 7.5),
    plot.caption    = element_text(size = 7.5, colour = "grey40", hjust = 0)
  )

# Scale figure height with number of outcomes (0.22 inches per row, min 8)
fig_height <- max(8, nrow(plot_df) * 0.22 + 2)
print(fig_obs)
save_fig(fig_obs, "fig_power_per_outcome", width = 11, height = fig_height)


# 11. SESSION INFO -------------------------------------------------------------

cat("=== SESSION INFO ===\n")
cat(sprintf("R version : %s\n", R.version$version.string))
cat(sprintf("pwr       : %s\n", as.character(packageVersion("pwr"))))
cat(sprintf("ggplot2   : %s\n", as.character(packageVersion("ggplot2"))))
cat(sprintf("set.seed  : 123\n"))
cat(sprintf("n_sim     : %d\n", n_sim))
cat(sprintf("Figures   : %s\n", figures_power_folder))
cat(sprintf("Tables    : %s\n", tables_power_folder))
cat(sprintf("Date      : %s\n", format(Sys.Date(), "%d %B %Y")))