# =============================================================================
# HYROX health-problem epidemiology
# Prospective 12-week OSTRC-H2 surveillance: full analysis
#
# All-injury PRIMARY analysis; HYROX-attributed injuries as a SENSITIVITY.
# Produces every number, table and figure in the manuscript and supplement.
#
# Requires: MASS, tidyverse, readxl, sandwich, broom, ggpubr
# Source the whole file; run order matters.
# =============================================================================

library(MASS)          # glm.nb - load BEFORE tidyverse so select() is masked back
library(tidyverse)
library(readxl)
library(sandwich)      # vcovHC, HC3 robust covariance
library(broom)
library(ggpubr)        # theme_pubr for the forest plot
select <- dplyr::select
set.seed(20260809)

filepath <- "data/HYROX_raw_data.xlsx"
out_dir  <- "figures"
dir.create(out_dir, showWarnings = FALSE)

REGIONS <- c("Head","Neck","Shoulder","Upper arm","Elbow","Forearm","Wrist",
             "Hand/fingers","Chest/ribs/upper back","Abdomen","Pelvis/lower back",
             "Hip/groin","Thigh","Knee","Lower leg/Achilles","Ankle","Foot","Other")
MODALITY <- c("competition_hrs","hyrox_conditioning_hours","strength_training_hours",
              "running_hours","stretching_hours")
TRAINING <- setdiff(MODALITY, "stretching_hours")   # denominator excluding recovery work
N_BOOT   <- 2000

theme_pub <- function(base_size = 13) {
  theme_bw(base_size = base_size) %+replace%
    theme(panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
          panel.grid.major = element_line(colour = "grey92", linewidth = 0.3),
          panel.grid.minor = element_blank(),
          axis.text = element_text(size = rel(0.9), colour = "black"),
          axis.title = element_text(size = rel(1.05), face = "bold"),
          axis.ticks = element_line(colour = "black", linewidth = 0.3),
          legend.background = element_blank(), legend.key = element_blank(),
          legend.position = "bottom", legend.text = element_text(size = rel(0.95)),
          plot.margin = margin(10, 10, 6, 6))
}

# Figure 1 bands: darkest innermost (time-loss), lightest outermost (all).
# NOTE: placeholder hex codes; substitute the published values before archiving.
pal_prev <- c("All problems" = "#F4A6C0", "Substantial problems" = "#C0392B",
              "Time-loss problems" = "#1A1A1A")

pois_rate <- function(k, hours) {
  ci <- poisson.test(k, hours / 1000)$conf.int
  sprintf("%.1f (%.1f to %.1f)", k / hours * 1000, ci[1], ci[2])
}

# =============================================================================
# 1. Import, de-duplicate, derive
#    Five participant-weeks were submitted twice (P045 wk5, P093 wks 1-3,
#    P105 wk1). Where duplicates disagree the record carrying a health problem
#    is retained, so no reported problem is discarded.
#    Five participants (P016, P061, P103, P107, P109) returned no baseline
#    questionnaire. They are retained for prevalence, incidence and exposure,
#    and are dropped only by the risk-marker model, which is the source of its
#    n of 84.
# =============================================================================
pss_raw    <- read_excel(filepath, sheet = "PSS")
weekly_raw <- read_excel(filepath, sheet = "Weekly")

dup_weeks <- weekly_raw %>% count(participant_id, week_number) %>% filter(n > 1)

pss <- pss_raw %>%
  mutate(sex = factor(sex, levels = c(1, 2), labels = c("Male", "Female")),
         age_cat = case_when(age_group %in% 1:2 ~ "16-29", age_group %in% 3:4 ~ "30-39",
                             age_group %in% 5:6 ~ "40-49", age_group >= 7 ~ "50+") %>%
           factor(levels = c("16-29","30-39","40-49","50+")),
         hyrox_exp_num = as.numeric(years_hyrox),
         hyrox_exp_cat = factor(years_hyrox, levels = 1:5,
                                labels = c("0-6 months","6 months - 1 year","1-2 years",
                                           "2-3 years","3+ years")),
         prev_injury_bin = factor(prev_injury_timeloss, levels = c(2, 1), labels = c("No","Yes")),
         no_baseline = is.na(age_group) & is.na(sex) & is.na(years_hyrox))

weekly <- weekly_raw %>%
  arrange(participant_id, week_number, desc(replace_na(any_problem, 0))) %>%
  distinct(participant_id, week_number, .keep_all = TRUE) %>%
  left_join(pss %>% select(participant_id, sex, age_cat, hyrox_exp_num, prev_injury_bin),
            by = "participant_id") %>%
  mutate(
    severity_wk = replace_na(OSTRC_severity_score, 0),
    hp_type     = factor(health_problem_type, levels = 1:3,
                         labels = c("Acute injury","Overuse injury","Illness")),
    any_wk      = replace_na(as.integer(any_problem == 1), 0L),
    injury_wk   = replace_na(as.integer(hp_type %in% c("Acute injury","Overuse injury") &
                                          any_problem == 1), 0L),
    illness_wk  = replace_na(as.integer(hp_type == "Illness" & any_problem == 1), 0L),
    acute_wk    = replace_na(as.integer(hp_type == "Acute injury" & any_problem == 1), 0L),
    overuse_wk  = replace_na(as.integer(hp_type == "Overuse injury" & any_problem == 1), 0L),
    hyrox_wk    = replace_na(as.integer(hyrox_injury == 1 & any_problem == 1), 0L),
    body_region = factor(injured_body_region, levels = 1:18, labels = REGIONS),
    mechanism   = factor(injury_mechanism, levels = 1:3,
                         labels = c("Strength","Running","HYROX station")),
    # Substantial: moderate/severe reduction in training (Q2) or performance (Q3),
    # i.e. item score >= 17, or any time-loss, so time-loss nests within substantial.
    tl_any     = replace_na(as.integer(timeloss_problem == 1), 0L),
    subst_flag = as.integer(replace_na(OSTRC_Q2_score, 0) >= 17 |
                              replace_na(OSTRC_Q3_score, 0) >= 17 | tl_any == 1),
    subst_any  = as.integer(any_wk == 1 & subst_flag == 1),
    hx_subst   = as.integer(hyrox_wk == 1 & subst_flag == 1),
    hx_tl      = as.integer(hyrox_wk == 1 & tl_any == 1),
    mod_hrs    = rowSums(across(all_of(MODALITY), ~replace_na(.x, 0))),
    train_hrs  = rowSums(across(all_of(TRAINING), ~replace_na(.x, 0)))
  ) %>%
  arrange(participant_id, week_number)

n_participants <- nrow(pss)
n_female       <- sum(pss$sex == "Female", na.rm = TRUE)
n_personweeks  <- nrow(weekly)
response_pct   <- round(n_personweeks / (n_participants * 12) * 100, 1)

# =============================================================================
# 2. Exposure denominators
#    Three totals coexist in the workbook. The participant-level summary is
#    retained as the primary denominator for continuity with the submitted
#    analysis; the weekly modality columns supply the modality-specific rates
#    and the recovery-excluded sensitivity.
# =============================================================================
exposure <- weekly %>%
  group_by(participant_id) %>%
  summarise(weeks = n(), across(all_of(MODALITY), ~sum(replace_na(.x, 0))),
            hrs_weekly_field = sum(replace_na(weekly_exposure_hrs, 0)),
            hrs_modality = sum(mod_hrs), hrs_training = sum(train_hrs), .groups = "drop") %>%
  right_join(pss %>% select(participant_id, hrs_pss = total_exposure_hrs, n_weeks_responded),
             by = "participant_id") %>%
  mutate(across(where(is.numeric), ~replace_na(.x, 0)))

denominators <- tibble(
  source = c("PSS total_exposure_hrs (primary)", "Weekly weekly_exposure_hrs",
             "Weekly modality columns", "Weekly modality excluding recovery work"),
  hours  = c(sum(exposure$hrs_pss), sum(exposure$hrs_weekly_field),
             sum(exposure$hrs_modality), sum(exposure$hrs_training)))

denominator_disagreements <- exposure %>%
  filter(abs(hrs_pss - hrs_modality) > 1e-6 | n_weeks_responded != weeks) %>%
  select(participant_id, weeks, n_weeks_responded, hrs_pss, hrs_weekly_field, hrs_modality)

modality_share <- weekly %>%
  summarise(across(all_of(MODALITY), ~sum(replace_na(.x, 0)))) %>%
  pivot_longer(everything(), names_to = "modality", values_to = "hours") %>%
  mutate(pct = round(hours / sum(hours) * 100, 1))

E_pub   <- sum(exposure$hrs_pss)
E_mod   <- sum(exposure$hrs_modality)
E_train <- sum(exposure$hrs_training)
ids_ir  <- exposure %>% filter(hrs_pss > 0) %>% pull(participant_id)

# =============================================================================
# 3. Episode reconstruction
#    (i)   An episode opens only at a week reporting a new injury onset.
#    (ii)  It continues through later injury-weeks that are not new onsets.
#    (iii) A returned week reporting no injury closes it.
#    (iv)  A change of body region closes it.
#    (v)   Missing weeks do not close it, unless max_gap is set.
#    Injury-weeks that cannot be attached to a recorded onset are counted as
#    orphan weeks and excluded from incidence. NOTE: new_injury is recorded as
#    a tick (1 or NA), so an unticked week cannot be distinguished from an
#    unanswered one; the orphan weeks may be unflagged recurrences. The
#    sensitivity analysis below counts them as such.
# =============================================================================
assign_episodes <- function(d, max_gap = Inf) {
  n <- nrow(d); id <- rep(NA_integer_, n); cur <- 0L; open <- FALSE
  for (i in seq_len(n)) {
    if (d$injury_wk[i] != 1) { open <- FALSE; next }
    if (isTRUE(d$new_injury[i] == 1)) { cur <- cur + 1L; open <- TRUE; id[i] <- cur; next }
    if (!open) next
    prev <- max(which(!is.na(id[seq_len(i - 1)])))
    same_region <- is.na(d$body_region[i]) || is.na(d$body_region[prev]) ||
      d$body_region[i] == d$body_region[prev]
    gap <- d$week_number[i] - d$week_number[prev] - 1
    if (same_region && gap <= max_gap) id[i] <- cur else open <- FALSE
  }
  id
}

tag_episodes <- function(w, max_gap = Inf) {
  w %>% group_by(participant_id) %>%
    group_modify(~mutate(.x, epi = assign_episodes(.x, max_gap))) %>% ungroup()
}

build_episodes <- function(w, max_gap = Inf) {
  d <- tag_episodes(w, max_gap)
  eps <- d %>% filter(injury_wk == 1, !is.na(epi)) %>%
    group_by(participant_id, epi) %>%
    summarise(body_region = first(body_region), hp_type = first(hp_type),
              mechanism = first(mechanism), sex = first(sex),
              hyrox = as.integer(first(hyrox_injury) == 1),
              cum_severity = sum(OSTRC_severity_score, na.rm = TRUE),
              total_days = sum(replace_na(as.numeric(timeloss_days), 0)),
              n_weeks = n(), first_week = min(week_number), last_week = max(week_number),
              .groups = "drop")
  attr(eps, "orphan_weeks") <- sum(d$injury_wk == 1 & is.na(d$epi))
  eps
}

injury_episodes <- build_episodes(weekly)
ep_ir <- injury_episodes %>% filter(participant_id %in% ids_ir) %>%
  mutate(sexg = ifelse(is.na(sex), "Unspecified", as.character(sex)))

# Illness episode: run of consecutive illness weeks, a gap of >1 week starting a new one
illness_episodes <- weekly %>% filter(illness_wk == 1) %>%
  group_by(participant_id) %>%
  mutate(ill_id = cumsum(as.integer(row_number() == 1 | (week_number - lag(week_number)) > 1))) %>%
  ungroup() %>% group_by(participant_id, ill_id) %>%
  summarise(total_days = sum(replace_na(as.numeric(timeloss_days), 0)),
            severity = sum(OSTRC_severity_score, na.rm = TRUE), n_weeks = n(), .groups = "drop")

# Sensitivity: alternative closure rules, plus the opposing assumption that
# unattached injury weeks are incident episodes.
episode_sensitivity <- map_dfr(
  list(`Missing weeks never close an episode` = Inf, `Gap > 2 weeks closes` = 2,
       `Gap > 1 week closes` = 1, `Any missing week closes` = 0),
  function(g) {
    e <- build_episodes(weekly, g)
    orph <- attr(e, "orphan_weeks")
    e <- filter(e, participant_id %in% ids_ir)
    tibble(episodes = nrow(e), time_loss = sum(e$total_days > 0), hyrox = sum(e$hyrox == 1),
           days_lost = sum(e$total_days), orphan_weeks = orph,
           incidence = pois_rate(nrow(e), E_pub),
           tl_incidence = pois_rate(sum(e$total_days > 0), E_pub))
  }, .id = "rule")

# Orphan runs that sit immediately before a flagged onset in the same region are
# late ticks and merge forward; the rest are treated as unflagged recurrences.
orphan_wk <- tag_episodes(weekly) %>% filter(injury_wk == 1, is.na(epi)) %>%
  select(participant_id, week_number, body_region)
onset_wk <- weekly %>% filter(injury_wk == 1, new_injury == 1) %>%
  select(participant_id, week_number, body_region)
late_tick <- orphan_wk %>%
  inner_join(onset_wk, by = "participant_id", suffix = c("", "_on")) %>%
  filter(week_number_on == week_number + 1,
         is.na(body_region) | is.na(body_region_on) | body_region == body_region_on) %>%
  distinct(participant_id, week_number)
orphan_runs <- orphan_wk %>% anti_join(late_tick, by = c("participant_id","week_number")) %>%
  group_by(participant_id) %>%
  summarise(runs = sum(c(TRUE, diff(week_number) != 1)), .groups = "drop")
episode_recurrence_rule <- tibble(
  episodes = nrow(ep_ir) + sum(orphan_runs$runs),
  late_ticks_merged = nrow(late_tick),
  incidence = pois_rate(nrow(ep_ir) + sum(orphan_runs$runs), E_pub))

# =============================================================================
# 4. Incidence
# =============================================================================
k_all <- nrow(ep_ir); k_tl <- sum(ep_ir$total_days > 0)
k_hx  <- sum(ep_ir$hyrox == 1); k_hx_tl <- sum(ep_ir$hyrox == 1 & ep_ir$total_days > 0)

incidence_summary <- tibble(
  measure = c("All injuries","All injuries (time-loss)",
              "HYROX-attributed","HYROX-attributed (time-loss)"),
  k = c(k_all, k_tl, k_hx, k_hx_tl)) %>%
  rowwise() %>% mutate(rate = pois_rate(k, E_pub)) %>% ungroup()

# Table 2: incidence by sex, Unspecified retained so rows sum to Combined
ir_by_sex <- exposure %>% filter(hrs_pss > 0) %>%
  left_join(pss %>% select(participant_id, sex), by = "participant_id") %>%
  mutate(sexg = ifelse(is.na(sex), "Unspecified", as.character(sex))) %>%
  group_by(sexg) %>% summarise(n = n(), hrs = sum(hrs_pss), .groups = "drop") %>%
  left_join(ep_ir %>% group_by(sexg) %>%
              summarise(new_all = n(), tl_all = sum(total_days > 0), .groups = "drop"),
            by = "sexg") %>%
  mutate(across(c(new_all, tl_all), ~replace_na(.x, 0))) %>%
  { bind_rows(., summarise(., sexg = "Combined", n = sum(n), hrs = sum(hrs),
                           new_all = sum(new_all), tl_all = sum(tl_all))) } %>%
  rowwise() %>%
  mutate(overall = pois_rate(new_all, hrs), timeloss = pois_rate(tl_all, hrs)) %>%
  ungroup() %>%
  mutate(sex = factor(sexg, levels = c("Male","Female","Unspecified","Combined"))) %>%
  arrange(sex) %>% select(sex, n, exposure_h = hrs, overall, timeloss)

# Supplementary Table S1: denominator and reconstruction sensitivity
incidence_by_denominator <- tibble(
  denominator = c("PSS total (primary)","Weekly modality total",
                  "Weekly modality excluding recovery work"),
  hours = c(E_pub, E_mod, E_train)) %>%
  rowwise() %>%
  mutate(all_injuries = pois_rate(k_all, hours), time_loss = pois_rate(k_tl, hours),
         hyrox_attributed = pois_rate(k_hx, hours)) %>% ungroup()

# Supplementary Table S2: exploratory modality-specific incidence.
# Numerator is the participant-reported mechanism; mechanism was missing for
# roughly half of episodes, and more often for injuries not attributed to HYROX.
mech_hours <- c(Running = sum(weekly$running_hours, na.rm = TRUE),
                Strength = sum(weekly$strength_training_hours, na.rm = TRUE),
                `HYROX station` = sum(weekly$hyrox_conditioning_hours, na.rm = TRUE))
incidence_by_modality <- tibble(mechanism = names(mech_hours), hours = as.numeric(mech_hours)) %>%
  left_join(ep_ir %>% count(mechanism, name = "episodes") %>%
              mutate(mechanism = as.character(mechanism)), by = "mechanism") %>%
  rowwise() %>% mutate(rate = pois_rate(episodes, hours)) %>% ungroup()

mechanism_missing <- tibble(
  episodes = nrow(ep_ir), missing = sum(is.na(ep_ir$mechanism)),
  pct_missing_hyrox = round(mean(is.na(ep_ir$mechanism[ep_ir$hyrox == 1])) * 100, 0),
  pct_missing_other = round(mean(is.na(ep_ir$mechanism[ep_ir$hyrox == 0])) * 100, 0))

# =============================================================================
# 5. Prevalence
#    Pooled across returned participant-weeks, with cluster bootstrap CIs
#    resampling participants with replacement.
# =============================================================================
prev_outcomes <- c("any_wk","subst_any","tl_any","hyrox_wk","hx_subst","hx_tl",
                   "overuse_wk","acute_wk","illness_wk")
stopifnot(all(prev_outcomes %in% names(weekly)))

by_pid <- split(weekly[, prev_outcomes, drop = FALSE], weekly$participant_id)
boot_prev <- map_dfr(seq_len(N_BOOT), function(i) {
  d <- bind_rows(by_pid[sample(names(by_pid), replace = TRUE)])
  summarise(d, across(everything(), ~mean(.x) * 100))
})
prevalence_pooled <- tibble(outcome = prev_outcomes,
                            estimate = map_dbl(prev_outcomes, ~mean(weekly[[.x]]) * 100),
                            ci_lo = map_dbl(prev_outcomes, ~quantile(boot_prev[[.x]], 0.025)),
                            ci_hi = map_dbl(prev_outcomes, ~quantile(boot_prev[[.x]], 0.975))) %>%
  mutate(across(where(is.numeric), ~round(.x, 1)))

n_with_problem <- weekly %>% group_by(participant_id) %>%
  summarise(any = max(any_wk), .groups = "drop") %>% summarise(sum(any)) %>% pull()

# =============================================================================
# 6. Non-response
# =============================================================================
resp <- pss %>% select(participant_id) %>%
  left_join(count(weekly, participant_id, name = "weeks"), by = "participant_id") %>%
  mutate(weeks = replace_na(weeks, 0L))

pattern <- weekly %>% group_by(participant_id) %>%
  summarise(k = n(), last = max(week_number), .groups = "drop") %>%
  mutate(monotone = (k == last))   # returned weeks form an unbroken run from week 1

response_pattern <- tibble(
  person_weeks = n_personweeks, duplicated_weeks = nrow(dup_weeks),
  response_rate_pct = response_pct, median_weeks = median(resp$weeks),
  iqr_low = quantile(resp$weeks, .25), iqr_high = quantile(resp$weeks, .75),
  completed_all_12 = sum(resp$weeks == 12), completed_9_or_more = sum(resp$weeks >= 9),
  completed_3_or_fewer = sum(resp$weeks <= 3),
  monotone_dropout = sum(pattern$monotone), of_n = nrow(pattern))

grid <- expand_grid(participant_id = pss$participant_id, week_number = 1:12) %>%
  left_join(weekly %>% select(participant_id, week_number, any_wk, subst_any, tl_any),
            by = c("participant_id","week_number")) %>%
  arrange(participant_id, week_number) %>% group_by(participant_id) %>%
  mutate(responded = as.integer(!is.na(any_wk)), resp_next = lead(responded)) %>% ungroup()

informative_missingness <- grid %>% filter(responded == 1, !is.na(resp_next)) %>%
  group_by(problem_in_week_t = any_wk == 1) %>%
  summarise(weeks = n(), responded_next_pct = round(mean(resp_next) * 100, 1), .groups = "drop")
nonresponse_test <- with(filter(grid, responded == 1, !is.na(resp_next)),
                         fisher.test(table(any_wk == 1, resp_next)))

smd <- function(x, g) {
  a <- x[g == 1]; b <- x[g == 0]; a <- a[!is.na(a)]; b <- b[!is.na(b)]
  s <- sqrt((var(a) + var(b)) / 2)
  if (is.na(s) || s == 0) NA_real_ else (mean(a) - mean(b)) / s
}
comp <- pss %>% left_join(resp, by = "participant_id") %>%
  mutate(high = as.integer(weeks >= median(resp$weeks)))

# Supplementary Table S3
compliance_table <- tibble(
  variable = c("Age group (1-9)","Female","HYROX experience (1-5)",
               "Previous time-loss injury","Baseline total training h/wk",
               "Baseline HYROX training h/wk","Competed in last 12 months"),
  value = list(comp$age_group, as.numeric(comp$sex == "Female"), comp$hyrox_exp_num,
               as.numeric(comp$prev_injury_bin == "Yes"), comp$total_training_hrs_wk,
               comp$hyrox_training_hrs_wk, as.numeric(comp$competed_last_12m == 1))) %>%
  rowwise() %>%
  mutate(high_compliance = round(mean(value[comp$high == 1], na.rm = TRUE), 2),
         low_compliance  = round(mean(value[comp$high == 0], na.rm = TRUE), 2),
         SMD = round(smd(value, comp$high), 2)) %>%
  ungroup() %>% select(-value)

# Supplementary Table S4. Interior-gap LOCF fills only gaps between two returned
# questionnaires; it does not extrapolate beyond a participant's final response.
last_wk <- weekly %>% group_by(participant_id) %>%
  summarise(last = max(week_number), .groups = "drop")
prevalence_sensitivity <- map_dfr(c("any_wk","subst_any","tl_any"), function(v) {
  cc <- mean(weekly[[v]]) * 100
  g <- grid %>% left_join(last_wk, by = "participant_id") %>%
    mutate(filled = .data[[v]]) %>% group_by(participant_id) %>%
    fill(filled, .direction = "down") %>% ungroup() %>%
    filter(week_number <= last, !is.na(filled))
  z <- mutate(grid, filled = replace_na(.data[[v]], 0))
  tibble(outcome = v, returned_weeks_only = round(cc, 1),
         interior_gap_locf = round(mean(g$filled) * 100, 1),
         nonresponse_assumed_healthy = round(mean(z$filled) * 100, 1))
})

# =============================================================================
# 7. Descriptives
# =============================================================================
med_iqr <- function(x) sprintf("%.1f (%.1f-%.1f)", median(x, na.rm = TRUE),
                               quantile(x, .25, na.rm = TRUE), quantile(x, .75, na.rm = TRUE))
n_pct   <- function(k) sprintf("%d (%.0f%%)", k, k / n_participants * 100)

# Table 1. Participants with no baseline questionnaire are reported as missing
# rather than absorbed into a category.
table1 <- bind_rows(
  tibble(characteristic = "Age group", value = ""),
  pss %>% count(age_cat, .drop = FALSE) %>%
    transmute(characteristic = ifelse(is.na(age_cat), "  Not reported", paste0("  ", age_cat)),
              value = n_pct(n)),
  tibble(characteristic = "Sex: female", value = n_pct(n_female)),
  tibble(characteristic = "  Sex not reported", value = n_pct(sum(is.na(pss$sex)))),
  tibble(characteristic = "HYROX experience", value = ""),
  pss %>% count(hyrox_exp_cat, .drop = FALSE) %>%
    transmute(characteristic = ifelse(is.na(hyrox_exp_cat), "  Not reported",
                                      paste0("  ", hyrox_exp_cat)), value = n_pct(n)),
  tibble(characteristic = c("Total training, h/wk","HYROX-specific training, h/wk",
                            "Total exposure, h","Weeks responded"),
         value = c(med_iqr(pss$total_training_hrs_wk), med_iqr(pss$hyrox_training_hrs_wk),
                   med_iqr(pss$total_exposure_hrs), med_iqr(resp$weeks))),
  tibble(characteristic = c("Previous injury (time-loss) in last 12 months",
                            "Competed in last 12 months"),
         value = c(n_pct(sum(pss$prev_injury_bin == "Yes", na.rm = TRUE)),
                   n_pct(sum(pss$competed_last_12m == 1, na.rm = TRUE)))))

exposure_variability <- tibble(
  measure = c("Weekly exposure, h","Total exposure, h",
              "Within-participant CV of weekly hours"),
  value = c(med_iqr(weekly$mod_hrs), med_iqr(exposure$hrs_modality),
            med_iqr(weekly %>% group_by(participant_id) %>% filter(n() >= 3) %>%
                      summarise(cv = sd(mod_hrs) / mean(mod_hrs), .groups = "drop") %>% pull(cv))))

# Table 3. Locations by type and IOC time-loss band
ioc_band <- function(days) cut(days, c(-.1, 0, 7, 28, Inf),
                               labels = c("Slight","Mild","Moderate","Severe"))
table3 <- ep_ir %>%
  mutate(region = fct_na_value_to_level(body_region, "Other/unspecified"),
         band = ioc_band(total_days)) %>%
  group_by(region, hp_type) %>%
  summarise(n = n(), slight = sum(band == "Slight"), mild = sum(band == "Mild"),
            moderate = sum(band == "Moderate"), severe = sum(band == "Severe"),
            days = sum(total_days), .groups = "drop") %>%
  arrange(hp_type, desc(n))
table3_illness <- illness_episodes %>%
  mutate(band = ioc_band(total_days)) %>%
  summarise(n = n(), slight = sum(band == "Slight"), mild = sum(band == "Mild"),
            moderate = sum(band == "Moderate"), severe = sum(band == "Severe"),
            days = sum(total_days))

region_freq <- ep_ir %>% count(body_region, sort = TRUE) %>%
  mutate(pct = round(n / sum(n) * 100, 1))
type_freq <- ep_ir %>% count(hp_type) %>% mutate(pct = round(n / sum(n) * 100, 1))
mech_freq <- ep_ir %>% filter(!is.na(mechanism)) %>% count(mechanism) %>%
  mutate(pct = round(n / sum(n) * 100, 1))

# =============================================================================
# 8. Burden
# =============================================================================
named_region <- function(d) filter(d, !is.na(body_region), body_region != "Other")

burden_by_region_type <- ep_ir %>% named_region() %>%
  group_by(body_region, hp_type) %>%
  summarise(n_new = n(), mean_sev = mean(cum_severity), .groups = "drop") %>%
  mutate(incidence = n_new / E_pub * 1000, burden = incidence * mean_sev,
         type = ifelse(hp_type == "Overuse injury", "Gradual-onset injury", "Acute injury"),
         label = as.character(body_region))
burden_all <- burden_by_region_type %>% filter(n_new >= 1) %>% slice_max(burden, n = 8)

burden_by_region <- ep_ir %>% named_region() %>%
  group_by(body_region) %>%
  summarise(n_new = n(), mean_sev = mean(cum_severity),
            days = sum(total_days), tl_eps = sum(total_days > 0),
            gradual_pct = round(mean(hp_type == "Overuse injury") * 100, 0),
            mech_recorded = sum(!is.na(mechanism)),
            mech_mix = paste(names(table(mechanism))[table(mechanism) > 0],
                             table(mechanism)[table(mechanism) > 0], collapse = ", "),
            .groups = "drop") %>%
  mutate(incidence = n_new / E_pub * 1000, burden = round(incidence * mean_sev, 0),
         pct_of_episodes = round(n_new / nrow(ep_ir) * 100, 0)) %>%
  arrange(desc(burden))

# Table 4 source rows: the six highest-burden regions
table4 <- burden_by_region %>% slice_head(n = 6) %>%
  select(body_region, n_new, pct_of_episodes, gradual_pct, days, mech_mix, mech_recorded, burden)

# How far the elbow ranking depends on individual episodes
elbow_episodes <- ep_ir %>% filter(body_region == "Elbow") %>%
  select(participant_id, hp_type, n_weeks, first_week, last_week, total_days,
         cum_severity, mechanism)

# =============================================================================
# 9. Risk-marker model
#    Negative binomial with log(exposure) offset and HC3 robust standard errors.
#    The five participants without a baseline questionnaire are dropped here.
# =============================================================================
cum_sev <- weekly %>% group_by(participant_id) %>%
  summarise(cum_severity_injury = sum(severity_wk * injury_wk),
            cum_severity_hyrox  = sum(severity_wk * hyrox_wk), .groups = "drop")

pss_rm <- pss %>% left_join(cum_sev, by = "participant_id") %>%
  left_join(exposure %>% select(participant_id, hrs_pss), by = "participant_id") %>%
  mutate(across(starts_with("cum_severity"), ~replace_na(.x, 0))) %>%
  filter(!is.na(sex), !is.na(age_cat), !is.na(hyrox_exp_num), !is.na(prev_injury_bin),
         hrs_pss > 0) %>%
  mutate(cum_sev_int = as.integer(round(cum_severity_injury)),
         cum_sev_int_hx = as.integer(round(cum_severity_hyrox)),
         ln_exposure_hrs = log(hrs_pss))

fit_nb <- function(outcome) {
  f <- as.formula(paste(outcome,
                        "~ sex + age_cat + prev_injury_bin + hyrox_exp_num + offset(ln_exposure_hrs)"))
  m <- tryCatch(glm.nb(f, data = pss_rm),
                error = function(e) glm(f, data = pss_rm, family = poisson))
  v <- sandwich::vcovHC(m, type = "HC3")
  tidy(m) %>% filter(term != "(Intercept)") %>%
    mutate(robust_se = sqrt(diag(v)[term]), RR = exp(estimate),
           ci_lo = exp(estimate - 1.96 * robust_se),
           ci_hi = exp(estimate + 1.96 * robust_se),
           p_val = 2 * pnorm(-abs(estimate / robust_se)))
}
nb_res    <- fit_nb("cum_sev_int")
nb_res_hx <- fit_nb("cum_sev_int_hx")

fmt_rr <- function(d, dp = 3) d %>%
  transmute(term, RR = round(RR, 2), ci_lo = round(ci_lo, dp),
            ci_hi = round(ci_hi, 2), p = round(p_val, dp))
regression_primary     <- fmt_rr(nb_res)
regression_sensitivity <- fmt_rr(nb_res_hx)

# =============================================================================
# 10. Figures
#     All three at base_size 13. The forest plot no longer encodes p < 0.05 in
#     the point fill, consistent with the exploratory framing of the model.
# =============================================================================

# ---- Figure 1: weekly prevalence, faceted bands ----
weekly_prev <- weekly %>% group_by(week_number) %>%
  summarise(n_resp = n(),
            prev_any = mean(any_wk) * 100,
            prev_any_sub = mean(subst_any) * 100,
            prev_any_tl = mean(tl_any) * 100,
            prev_overuse = mean(overuse_wk) * 100,
            prev_overuse_sub = mean(overuse_wk == 1 & subst_flag == 1) * 100,
            prev_overuse_tl = mean(overuse_wk == 1 & tl_any == 1) * 100,
            prev_acute = mean(acute_wk) * 100,
            prev_acute_sub = mean(acute_wk == 1 & subst_flag == 1) * 100,
            prev_acute_tl = mean(acute_wk == 1 & tl_any == 1) * 100,
            prev_illness = mean(illness_wk) * 100,
            prev_illness_sub = mean(illness_wk == 1 & subst_flag == 1) * 100,
            prev_illness_tl = mean(illness_wk == 1 & tl_any == 1) * 100, .groups = "drop")

fig1_long <- weekly_prev %>%
  pivot_longer(-c(week_number, n_resp), names_to = "var", values_to = "prevalence") %>%
  mutate(panel = case_when(str_detect(var, "^prev_any") ~ "All health problems",
                           str_detect(var, "^prev_overuse") ~ "Overuse injury",
                           str_detect(var, "^prev_acute") ~ "Acute injury",
                           str_detect(var, "^prev_illness") ~ "Illness"),
         level = case_when(str_detect(var, "_tl$") ~ "Time-loss problems",
                           str_detect(var, "_sub$") ~ "Substantial problems",
                           TRUE ~ "All problems")) %>%
  filter(!is.na(panel)) %>%
  mutate(panel = factor(panel, levels = c("All health problems","Overuse injury",
                                          "Acute injury","Illness")),
         level = factor(level, levels = c("All problems","Substantial problems",
                                          "Time-loss problems")))

# Spline-interpolate between weeks, then enforce nesting (interpolation overshoots)
fig1_wide <- fig1_long %>% group_by(panel, level) %>%
  reframe(week_smooth = seq(min(week_number), max(week_number), length.out = 300),
          prevalence = spline(week_number, prevalence, xout = week_smooth,
                              method = "natural")$y) %>%
  mutate(prevalence = pmax(prevalence, 0)) %>%
  pivot_wider(names_from = level, values_from = prevalence) %>%
  group_by(panel) %>%
  mutate(`Substantial problems` = pmin(`Substantial problems`, `All problems`),
         `Time-loss problems` = pmin(`Time-loss problems`, `Substantial problems`)) %>%
  ungroup() %>% arrange(panel, week_smooth)

# Explicit stacked bands, so there is no draw-order dependence
fig1_bands <- bind_rows(
  fig1_wide %>% transmute(panel, week_smooth, level = "Time-loss problems",
                          ymin = 0, ymax = `Time-loss problems`),
  fig1_wide %>% transmute(panel, week_smooth, level = "Substantial problems",
                          ymin = `Time-loss problems`, ymax = `Substantial problems`),
  fig1_wide %>% transmute(panel, week_smooth, level = "All problems",
                          ymin = `Substantial problems`, ymax = `All problems`)) %>%
  mutate(level = factor(level, levels = c("All problems","Substantial problems",
                                          "Time-loss problems")))

p_prev <- ggplot(fig1_bands, aes(x = week_smooth, ymin = ymin, ymax = ymax, fill = level)) +
  geom_ribbon() + facet_wrap(~panel, ncol = 1) +
  scale_x_continuous(breaks = 1:12, expand = expansion(mult = c(0.02, 0.02))) +
  scale_y_continuous(limits = c(0, 100), labels = function(x) paste0(x, "%"),
                     expand = expansion(mult = c(0, 0.02))) +
  scale_fill_manual(values = pal_prev, name = NULL,
                    breaks = c("All problems","Substantial problems","Time-loss problems")) +
  labs(x = "Week", y = "Prevalence") +
  theme_pub(13) +
  theme(strip.text = element_text(size = rel(1.05), face = "bold"),
        panel.spacing = unit(0.6, "lines"))
ggsave(file.path(out_dir, "fig1_prevalence.tiff"), p_prev, width = 6.9, height = 7.3, dpi = 300)

# ---- Figure 2: burden matrix ----
x_max <- max(burden_all$incidence) * 1.35
y_max <- max(burden_all$mean_sev) * 1.25
grid_bg <- expand.grid(x = seq(0, x_max, length.out = 300),
                       y = seq(0, y_max, length.out = 300)) %>% mutate(burden = x * y)

p_burden <- ggplot() +
  geom_raster(data = grid_bg, aes(x, y, fill = burden), interpolate = TRUE) +
  geom_contour(data = grid_bg, aes(x, y, z = burden), colour = "grey90", linewidth = 0.3,
               breaks = pretty(range(grid_bg$burden), n = 8)) +
  scale_fill_gradientn(colours = c("white","#FFF3C4","#FFD54F","#FFB300","#E65100"),
                       guide = "none") +
  geom_point(data = burden_all, aes(incidence, mean_sev, shape = type),
             size = 4, colour = "black", stroke = 0.4) +
  scale_shape_manual(values = c("Gradual-onset injury" = 16, "Acute injury" = 15), name = NULL) +
  geom_text(data = burden_all, aes(incidence, mean_sev, label = label),
            hjust = -0.12, vjust = -0.7, size = 4.1, colour = "grey15") +
  scale_x_continuous(expand = expansion(mult = c(0, 0))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0))) +
  coord_cartesian(xlim = c(0, x_max), ylim = c(0, y_max)) +
  labs(x = "Incidence (new injuries per 1000 h)",
       y = "Severity (mean cumulative severity score)") +
  theme_pub(13) + theme(panel.grid = element_blank())
ggsave(file.path(out_dir, "fig2_burden.tiff"), p_burden, width = 7, height = 5.5, dpi = 300)

# ---- Figure 3: forest plot, primary all-injury model ----
fig3_data <- bind_rows(
  tibble(RR = 1, ci_lo = NA_real_, ci_hi = NA_real_,
         label = c("Sex: Male (ref)","Age: 16-29 (ref)","Previous injury: No (ref)"),
         order = c(1, 3, 7)),
  nb_res %>% mutate(label = case_match(term,
                                       "sexFemale" ~ "Sex: Female",
                                       "prev_injury_binYes" ~ "Previous injury: Yes",
                                       "hyrox_exp_num" ~ "HYROX experience (per category)",
                                       "age_cat30-39" ~ "Age: 30-39", "age_cat40-49" ~ "Age: 40-49",
                                       "age_cat50+" ~ "Age: 50+", .default = term),
                    order = case_match(term, "sexFemale" ~ 2, "age_cat30-39" ~ 4,
                                       "age_cat40-49" ~ 5, "age_cat50+" ~ 6, "prev_injury_binYes" ~ 8,
                                       "hyrox_exp_num" ~ 9, .default = 10)) %>%
    select(RR, ci_lo, ci_hi, label, order)) %>%
  arrange(desc(order)) %>% mutate(label = factor(label, levels = label))

p_forest <- ggplot(fig3_data, aes(RR, label)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey60", linewidth = 0.4) +
  geom_errorbarh(aes(xmin = ci_lo, xmax = ci_hi), height = 0.25, colour = "grey40",
                 linewidth = 0.5, na.rm = TRUE) +
  geom_point(size = 3.4, shape = 21, fill = "white", colour = "black", stroke = 0.5,
             na.rm = TRUE) +
  scale_x_log10() + labs(x = "Rate ratio", y = NULL) +
  ggpubr::theme_pubr(base_size = 13, legend = "none") +
  theme(panel.grid.major.x = element_line(colour = "grey92", linewidth = 0.3))
ggsave(file.path(out_dir, "fig3_risk_markers.tiff"), p_forest, width = 7, height = 5.5, dpi = 300)

# =============================================================================
# 11. Output
# =============================================================================
cat(sprintf("n = %d (%d female, %d without baseline data); %d participant-weeks (%d duplicates removed); response %.1f%%\n",
            n_participants, n_female, sum(pss$no_baseline), n_personweeks,
            nrow(dup_weeks), response_pct))
cat(sprintf("%d of %d participants (%.0f%%) reported at least one health problem\n",
            n_with_problem, n_participants, n_with_problem / n_participants * 100))
cat(sprintf("Days lost: injuries %d, illness %d (illness %.0f%% of all days)\n\n",
            round(sum(ep_ir$total_days)), round(sum(illness_episodes$total_days)),
            sum(illness_episodes$total_days) /
              (sum(ep_ir$total_days) + sum(illness_episodes$total_days)) * 100))

print(table1)
print(denominators); print(modality_share); print(denominator_disagreements)
print(exposure_variability)
print(incidence_summary); print(ir_by_sex)
print(episode_sensitivity); print(episode_recurrence_rule); print(orphan_wk)
print(incidence_by_denominator); print(incidence_by_modality); print(mechanism_missing)
print(prevalence_pooled)
print(response_pattern); print(informative_missingness); print(nonresponse_test)
print(compliance_table); print(prevalence_sensitivity)
print(region_freq); print(type_freq); print(mech_freq)
print(table3); print(table3_illness)
print(burden_by_region); print(table4); print(elbow_episodes)
print(regression_primary); print(regression_sensitivity)

cat(sprintf("\nFigures written to %s/fig1_prevalence.tiff, %s/fig2_burden.tiff, %s/fig3_risk_markers.tiff\n",
            out_dir, out_dir, out_dir))