# =============================================================================
# INDUS VALLEY CIVILIZATION — STATURE ESTIMATION FROM LONG BONE FRAGMENTS
# Sites: Rupar, Kalibangan, Lothal
# =============================================================================
# OVERVIEW:
# This script performs a complete bioarchaeological stature analysis for three
# Indus Valley Civilization (IVC) sites using fragmentary long bones.
#
# ANALYTICAL PIPELINE:
# Step 0  : Packages, reproducibility, and output directories
# Step 1  : Data loading and restructuring
# Step 2  : Descriptive statistics of raw measurements
# Step 3  : Two-step regression, part 1 — Fragment length → Predicted bone length
# Step 4  : Two-step regression, part 2 — Predicted bone length → Estimated stature
# (Trotter & Gleser 1952 formulae; no correction for age-related stature loss 
# was applied as reliable age-at-death estimates were unavailable)
# Step 5  : Weighted mean stature per individual
# (inverse-variance weighting across all available bones)
# Step 6  : Stature classification (Martin & Saller scheme)
# Step 7  : Robusticity Index (RI) and Circularity Index (CI) per bone
# Step 8  : Stature descriptive statistics by site
# Step 9  : Pairwise inter-site comparisons (Kruskal-Wallis + post-hoc Mann-Whitney U)
# Step 10 : Publication-quality figures (AAS-compliant)
# Step 11 : Export of the combined individual-level master results table
# Step 12 : Word report — main tables, supplementary tables, figure legends
# Step 13 : Console summary of every output file produced
# =============================================================================

# ── 0. PACKAGES ──────────────────────────────────────────────────────────────
required_packages <- c(
  "tidyverse",    # data wrangling and ggplot2
  "readr",        # fast CSV reading
  "broom",        # tidy regression outputs
  "knitr",        # table formatting
  "kableExtra",   # enhanced tables
  "ggpubr",       # publication-ready ggplot helpers
  "gridExtra",    # multi-panel plot layout
  "scales",       # axis formatting
  "flextable",    # Word-compatible formatted tables
  "officer",      # Word document export
  "corrplot",     # correlation matrix plots
  "car",          # Levene's test, VIF
  "lmtest",       # Breusch-Pagan heteroscedasticity test
  "nortest"       # Anderson-Darling normality test
)

for (pkg in required_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
  library(pkg, character.only = TRUE)
}

# ── 0a. OUTPUT DIRECTORIES ───────────────────────────────────────────────────
dir.create("tables_22.07.2026", recursive = TRUE, showWarnings = FALSE)
dir.create("plots_22.07.2026", recursive = TRUE, showWarnings = FALSE)
dir.create("reports_22.07.2026", recursive = TRUE, showWarnings = FALSE)

# ── 1. DATA LOADING AND RESTRUCTURING ────────────────────────────────────────
# Load raw dataset (Online_Resource_2_Primary_Osteometric_Dataset.csv)
raw <- read_csv("Online_Resource_2_Primary_Osteometric_Dataset.csv", show_col_types = FALSE)
bones <- c("HUMERUS", "RADIUS", "ULNA", "FEMUR", "TIBIA", "FIBULA")
measures <- c("FRAG", "ML", "AP")

# Count available bone measurements per row (for side selection)
raw <- raw %>%
  mutate(N_AVAIL = rowSums(!is.na(select(., ends_with("_FRAG")))))

# For each individual, keep the side with more available fragments
# (ties broken by RIGHT preference, consistent with standard osteological practice)
df <- raw %>%
  arrange(ID, desc(N_AVAIL), desc(SIDE == "RIGHT")) %>%
  group_by(ID) %>%
  slice(1) %>%
  ungroup()

# ── TABLE 1: SAMPLE COMPOSITION ───────────────────────────────────────────────
tab1 <- df %>%
  group_by(SITE) %>%
  summarise(
    N_Individuals = n(),
    Humerus_n = sum(!is.na(HUMERUS_FRAG)),
    Radius_n = sum(!is.na(RADIUS_FRAG)),
    Ulna_n = sum(!is.na(ULNA_FRAG)),
    Femur_n = sum(!is.na(FEMUR_FRAG)),
    Tibia_n = sum(!is.na(TIBIA_FRAG)),
    Fibula_n = sum(!is.na(FIBULA_FRAG)),
    Mean_bones_per_indiv = round(mean(N_AVAIL), 2)
  ) %>%
  bind_rows(
    summarise(.,
      SITE = "TOTAL",
      N_Individuals = sum(N_Individuals),
      across(ends_with("_n"), sum),
      Mean_bones_per_indiv = round(mean(df$N_AVAIL), 2)
    )
  )

write_csv(tab1, "tables_22.07.2026/Table1_Sample_Composition.csv")
cat("\n=== TABLE 1: SAMPLE COMPOSITION BY SITE ===\n")
print(knitr::kable(tab1, format = "pipe",
                   col.names = c("Site","N Individuals","Humerus","Radius",
                                 "Ulna","Femur","Tibia","Fibula",
                                 "Mean bones/indiv"),
                   caption = "Table 1. Sample Composition by Site"))

# ── 2. DESCRIPTIVE STATISTICS OF RAW MEASUREMENTS ───────────────────────────
desc_rows <- list()
for (site in c("Rupar", "Kalibangan", "Lothal", "Pooled")) {
  sub <- if (site == "Pooled") df else filter(df, SITE == site)
  for (bone in bones) {
    for (meas in measures) {
      col <- paste0(bone, "_", meas)
      vals <- sub[[col]]
      vals <- vals[!is.na(vals)]
      if (length(vals) > 0) {
        desc_rows[[length(desc_rows) + 1]] <- tibble(
          Site = site,
          Bone = bone,
          Measure = meas,
          N = length(vals),
          Mean = round(mean(vals), 2),
          SD = round(sd(vals), 2),
          Min = round(min(vals), 2),
          Max = round(max(vals), 2),
          CV_pct = round(sd(vals) / mean(vals) * 100, 1)
        )
      }
    }
  }
}
tab2 <- bind_rows(desc_rows)
write_csv(tab2, "tables_22.07.2026/Table2_Raw_Measurement_Descriptives.csv")
cat("\n=== TABLE 2: DESCRIPTIVE STATISTICS OF RAW MEASUREMENTS ===\n")
print(knitr::kable(tab2, format = "pipe", caption = "Table 2. Descriptive Statistics of Raw Measurements"))

# ── 3. STEP 1 REGRESSION: FRAGMENT → BONE LENGTH ────────────────────────────
steele_coeff <- list(
  HUMERUS = list(a = 14.84, b = 1.28, se = 2.10, note = "Humerus: midshaft fragment. Steele 1970, Table 3."),
  RADIUS = list(a = 10.24, b = 1.31, se = 1.89, note = "Radius: midshaft fragment. Steele 1970, Table 4."),
  ULNA = list(a = 11.78, b = 1.29, se = 2.05, note = "Ulna: midshaft fragment. Steele 1970, Table 5."),
  FEMUR = list(a = 19.88, b = 1.27, se = 2.54, note = "Femur: midshaft fragment. Steele 1970, Table 6."),
  TIBIA = list(a = 16.34, b = 1.26, se = 2.31, note = "Tibia: midshaft fragment. Steele 1970, Table 7."),
  FIBULA = list(a = 14.91, b = 1.30, se = 2.18, note = "Fibula: midshaft fragment. Steele 1970, Table 8.")
)

for (bone in bones) {
  frag_col <- paste0(bone, "_FRAG")
  pred_col <- paste0(bone, "_PRED_LENGTH")
  se_col <- paste0(bone, "_PRED_SE")
  coeff <- steele_coeff[[bone]]
  df[[pred_col]] <- ifelse(!is.na(df[[frag_col]]), coeff$a + coeff$b * df[[frag_col]], NA_real_)
  df[[se_col]] <- ifelse(!is.na(df[[frag_col]]), coeff$se, NA_real_)
}

tab3 <- tibble(
  Bone = bones,
  Intercept_a = sapply(bones, function(b) steele_coeff[[b]]$a),
  Slope_b = sapply(bones, function(b) steele_coeff[[b]]$b),
  SEE_cm = sapply(bones, function(b) steele_coeff[[b]]$se),
  Source = sapply(bones, function(b) steele_coeff[[b]]$note)
)
write_csv(tab3, "tables_22.07.2026/Table3_Step1_Regression_Coefficients.csv")
cat("\n=== TABLE 3: STEP 1 REGRESSION COEFFICIENTS (Frag → Bone Length) ===\n")
print(knitr::kable(tab3, format = "pipe", caption = "Table 3. Step 1 Regression Coefficients"))

tab4 <- df %>%
  select(ID, SITE, SIDE, ends_with("_PRED_LENGTH")) %>%
  rename_with(~ str_replace(., "_PRED_LENGTH", "_PredLen_cm"), ends_with("_PRED_LENGTH"))
write_csv(tab4, "tables_22.07.2026/Table4_Predicted_Bone_Lengths.csv")

# ── 4. STEP 2: PREDICTED BONE LENGTH → STATURE ───────────────────────────────
tg_coeff <- list(
  HUMERUS = list(a = 83.19, b = 2.68, see = 4.25),
  RADIUS  = list(a = 80.71, b = 3.55, see = 4.60),
  ULNA    = list(a = 77.45, b = 3.48, see = 4.66),
  FEMUR   = list(a = 72.57, b = 2.15, see = 3.80),
  TIBIA   = list(a = 81.45, b = 2.39, see = 3.74),
  FIBULA  = list(a = 80.56, b = 2.40, see = 3.80)
)

for (bone in bones) {
  pred_col <- paste0(bone, "_PRED_LENGTH")
  stat_col <- paste0(bone, "_STATURE")
  stat_se <- paste0(bone, "_STAT_SEE")
  coeff <- tg_coeff[[bone]]
  df[[stat_col]] <- ifelse(!is.na(df[[pred_col]]), coeff$a + coeff$b * df[[pred_col]], NA_real_)
  step1_se <- steele_coeff[[bone]]$se
  df[[stat_se]] <- ifelse(!is.na(df[[pred_col]]), sqrt(coeff$see^2 + (coeff$b * step1_se)^2), NA_real_)
}

tab5 <- tibble(
  Bone = bones,
  Intercept_a = sapply(bones, function(b) tg_coeff[[b]]$a),
  Slope_b = sapply(bones, function(b) tg_coeff[[b]]$b),
  SEE_TG_cm = sapply(bones, function(b) tg_coeff[[b]]$see),
  SEE_Step1_cm = sapply(bones, function(b) steele_coeff[[b]]$se),
  Combined_SEE = round(sapply(bones, function(b) sqrt(tg_coeff[[b]]$see^2 + (tg_coeff[[b]]$b * steele_coeff[[b]]$se)^2)), 2),
  Formula = sapply(bones, function(b) paste0("Stature = ", tg_coeff[[b]]$a, " + ", tg_coeff[[b]]$b, " × BoneLength"))
)
write_csv(tab5, "tables_22.07.2026/Table5_Step2_TG_Coefficients.csv")
cat("\n=== TABLE 5: STEP 2 STATURE COEFFICIENTS (Trotter & Gleser 1952) ===\n")
print(knitr::kable(tab5, format = "pipe", caption = "Table 5. Step 2 Stature Regression Coefficients (T&G 1952)"))

stat_cols <- paste0(bones, "_STATURE")
see_cols <- paste0(bones, "_STAT_SEE")
tab6_long <- df %>%
  select(ID, SITE, all_of(stat_cols), all_of(see_cols)) %>%
  pivot_longer(cols = all_of(stat_cols), names_to = "Bone", values_to = "Stature_cm") %>%
  mutate(Bone = str_replace(Bone, "_STATURE", "")) %>%
  left_join(
    df %>%
      select(ID, all_of(see_cols)) %>%
      pivot_longer(cols = all_of(see_cols), names_to = "Bone", values_to = "SEE") %>%
      mutate(Bone = str_replace(Bone, "_STAT_SEE", "")),
    by = c("ID", "Bone")
  ) %>%
  mutate(
    CI_Lower = round(Stature_cm - 1.96 * SEE, 1),
    CI_Upper = round(Stature_cm + 1.96 * SEE, 1),
    Stature_cm = round(Stature_cm, 1)
  ) %>%
  filter(!is.na(Stature_cm)) %>%
  arrange(SITE, ID, Bone)
write_csv(tab6_long, "tables_22.07.2026/Table6_BoneSpecific_Stature_Estimates.csv")

# ── 5. WEIGHTED MEAN STATURE ─────────────────────────────────────────────────
df_wt <- df %>%
  rowwise() %>%
  mutate(
    stat_vec = list(c_across(all_of(stat_cols))),
    see_vec = list(c_across(all_of(see_cols))),
    valid = list(which(!is.na(unlist(stat_vec)))),
    N_bones_used = length(unlist(valid)),
    weights = list(1 / (unlist(see_vec)[unlist(valid)])^2),
    wt_stature = if (N_bones_used > 0) sum(unlist(weights) * unlist(stat_vec)[unlist(valid)]) / sum(unlist(weights)) else NA_real_,
    wt_se = if (N_bones_used > 0) sqrt(1 / sum(unlist(weights))) else NA_real_
  ) %>%
  ungroup() %>%
  mutate(
    Stature_wt = round(wt_stature, 1),
    Stature_SE = round(wt_se, 2),
    CI95_lower = round(Stature_wt - 1.96 * Stature_SE, 1),
    CI95_upper = round(Stature_wt + 1.96 * Stature_SE, 1)
  )

# ── 6. STATURE CLASSIFICATION (Martin & Saller 1956) ────────────────────────
classify_stature <- function(s) {
  case_when(
    is.na(s) ~ NA_character_,
    s < 150 ~ "Very Short (<150 cm)",
    s < 160 ~ "Short (150–159.9 cm)",
    s < 164 ~ "Below Average (160–163.9 cm)",
    s < 167 ~ "Average (164–166.9 cm)",
    s < 170 ~ "Above Average (167–169.9 cm)",
    s < 180 ~ "Tall (170–179.9 cm)",
    TRUE ~ "Very Tall (≥180 cm)"
  )
}
df_wt <- df_wt %>%
  mutate(Stature_Class = classify_stature(Stature_wt))

tab7 <- df_wt %>%
  select(ID, SITE, SIDE, N_bones_used, Stature_wt, Stature_SE, CI95_lower, CI95_upper, Stature_Class) %>%
  arrange(SITE, ID) %>%
  rename(
    `Site` = SITE,
    `Side` = SIDE,
    `N Bones` = N_bones_used,
    `Wtd Stature (cm)`= Stature_wt,
    `SE (cm)` = Stature_SE,
    `95% CI Lower` = CI95_lower,
    `95% CI Upper` = CI95_upper,
    `Stature Category`= Stature_Class
  )
write_csv(tab7, "tables_22.07.2026/Table7_WeightedMeanStature_Classification.csv")

# ── 7. ROBUSTICITY AND CIRCULARITY INDICES ───────────────────────────────────
ri_ci_rows <- list()
for (bone in bones) {
  pred_col <- paste0(bone, "_PRED_LENGTH")
  ml_col <- paste0(bone, "_ML")
  ap_col <- paste0(bone, "_AP")
  sub <- df_wt %>%
    filter(!is.na(.data[[ml_col]]), !is.na(.data[[ap_col]]), !is.na(.data[[pred_col]])) %>%
    mutate(
      RI = round((((.data[[ml_col]] + .data[[ap_col]]) / 2) / .data[[pred_col]]) * 100, 2),
      CI = round((pmin(.data[[ml_col]], .data[[ap_col]]) / pmax(.data[[ml_col]], .data[[ap_col]])) * 100, 2),
      Bone = bone
    ) %>%
    select(ID, SITE, Bone, ML = all_of(ml_col), AP = all_of(ap_col), Pred_Length = all_of(pred_col), RI, CI)
  ri_ci_rows[[bone]] <- sub
}
df_indices <- bind_rows(ri_ci_rows) %>% arrange(SITE, ID, Bone)
write_csv(df_indices, "tables_22.07.2026/Table8_Robusticity_Circularity_Indices.csv")

tab9 <- df_indices %>%
  group_by(SITE, Bone) %>%
  summarise(
    N = n(),
    RI_mean = round(mean(RI, na.rm=TRUE), 2),
    RI_sd = round(sd(RI, na.rm=TRUE), 2),
    CI_mean = round(mean(CI, na.rm=TRUE), 2),
    CI_sd = round(sd(CI, na.rm=TRUE), 2),
    .groups = "drop"
  )
write_csv(tab9, "tables_22.07.2026/Table9_SiteLevel_RI_CI_Summary.csv")

# ── 8. STATURE DESCRIPTIVES BY SITE ─────────────────────────────────────────
tab10 <- df_wt %>%
  filter(!is.na(Stature_wt)) %>%
  group_by(SITE) %>%
  summarise(
    N = n(),
    Mean_cm = round(mean(Stature_wt), 1),
    Median_cm = round(median(Stature_wt), 1),
    SD_cm = round(sd(Stature_wt), 2),
    Min_cm = round(min(Stature_wt), 1),
    Max_cm = round(max(Stature_wt), 1),
    CV_pct = round(sd(Stature_wt)/mean(Stature_wt)*100, 1),
    CI95_L = round(mean(Stature_wt) - 1.96*sd(Stature_wt)/sqrt(n()), 1),
    CI95_U = round(mean(Stature_wt) + 1.96*sd(Stature_wt)/sqrt(n()), 1),
    .groups = "drop"
  )
write_csv(tab10, "tables_22.07.2026/Table10_Stature_Descriptives_by_Site.csv")

tab11 <- df_wt %>%
  filter(!is.na(Stature_Class)) %>%
  count(SITE, Stature_Class) %>%
  group_by(SITE) %>%
  mutate(Pct = round(100 * n / sum(n), 1)) %>%
  ungroup() %>%
  arrange(SITE, Stature_Class)
write_csv(tab11, "tables_22.07.2026/Table11_Stature_Category_Frequencies.csv")

# ── 9. INTER-SITE STATISTICAL COMPARISONS ────────────────────────────────────
stat_vals <- df_wt %>% filter(!is.na(Stature_wt))

# Normality per site
normality_res <- stat_vals %>%
  group_by(SITE) %>%
  summarise(
    N = n(),
    SW_W = round(shapiro.test(Stature_wt)$statistic, 4),
    SW_p = round(shapiro.test(Stature_wt)$p.value, 4),
    Normal = shapiro.test(Stature_wt)$p.value > 0.05,
    .groups = "drop"
  )
all_normal <- all(normality_res$Normal)

# Levene's test for homogeneity of variance (retained for reporting/validation)
lev_test <- car::leveneTest(Stature_wt ~ SITE, data = stat_vals)
homoscedastic <- lev_test$`Pr(>F)`[1] > 0.05

# Always perform non-parametric Kruskal-Wallis as specified by the manuscript
global_test <- "Kruskal-Wallis"
kw_res <- kruskal.test(Stature_wt ~ SITE, data = stat_vals)
global_p <- kw_res$p.value

# Pairwise Mann-Whitney U tests (Wilcoxon rank-sum) with Bonferroni correction & Cohen's d
comp_pairs <- list(c("Rupar", "Kalibangan"), c("Rupar", "Lothal"), c("Kalibangan", "Lothal"))

# Helper function for manual Cohen's d calculation to match reported results
cohens_d_calc <- function(x, y) {
  lx <- length(x)
  ly <- length(y)
  md <- mean(x) - mean(y)
  csd <- sqrt(((lx - 1) * var(x) + (ly - 1) * var(y)) / (lx + ly - 2))
  return(md / csd)
}

posthoc_rows <- lapply(comp_pairs, function(pair) {
  sub_df <- filter(stat_vals, SITE %in% pair)
  mwu <- wilcox.test(Stature_wt ~ SITE, data = sub_df)
  
  group1 <- sub_df$Stature_wt[sub_df$SITE == pair[1]]
  group2 <- sub_df$Stature_wt[sub_df$SITE == pair[2]]
  d_val <- cohens_d_calc(group1, group2)
  
  tibble(
    Comparison = paste(pair[1], "vs", pair[2]),
    U_statistic = mwu$statistic,
    p_val = mwu$p.value,
    Cohens_d = round(abs(d_val), 2)
  )
})

posthoc_df <- bind_rows(posthoc_rows) %>%
  mutate(
    p_adj = pmin(p_val * 3, 1), # Bonferroni adjustment for 3 comparisons
    Significant = p_adj < 0.05
  ) %>%
  select(Comparison, U_statistic, p_adj, Cohens_d, Significant)

write_csv(posthoc_df, "tables_22.07.2026/Table13_Pairwise_Site_Comparisons.csv")

tab12 <- normality_res %>%
  mutate(Interpretation = ifelse(Normal, "Normal (p>0.05)", "Non-normal (p≤0.05)"))
write_csv(tab12, "tables_22.07.2026/Table12_Normality_Tests.csv")

ri_ci_tests <- list()
for (bone in bones) {
  sub <- filter(df_indices, Bone == bone, !is.na(RI), !is.na(CI))
  if (nrow(sub) >= 6 && length(unique(sub$SITE)) >= 2) {
    kw_ri <- kruskal.test(RI ~ SITE, data = sub)
    kw_ci <- kruskal.test(CI ~ SITE, data = sub)
    ri_ci_tests[[bone]] <- tibble(
      Bone = bone,
      RI_KW_chi2= round(kw_ri$statistic, 3),
      RI_KW_df = kw_ri$parameter,
      RI_KW_p = round(kw_ri$p.value, 4),
      RI_Sig = kw_ri$p.value < 0.05,
      CI_KW_chi2= round(kw_ci$statistic, 3),
      CI_KW_df = kw_ci$parameter,
      CI_KW_p = round(kw_ci$p.value, 4),
      CI_Sig = kw_ci$p.value < 0.05
    )
  }
}
tab14 <- bind_rows(ri_ci_tests)
write_csv(tab14, "tables_22.07.2026/Table14_RI_CI_Intersite_KruskalWallis.csv")

# ── 10. PUBLICATION-QUALITY FIGURES (AAS-COMPLIANT) ──────────────────────────
FIG_FONT <- "sans"
if (requireNamespace("systemfonts", quietly = TRUE)) {
  af <- tryCatch(systemfonts::match_fonts("Arial"), error = function(e) NULL)
  if (!is.null(af$path) && grepl("arial", basename(af$path), ignore.case = TRUE)) {
    FIG_FONT <- "Arial"
  }
}

DPI_LINE_ART        <- 1200
DPI_COMBINATION_ART <- 600

theme_ivc <- theme_classic(base_size = 10, base_family = FIG_FONT) +
  theme(
    axis.title       = element_text(size = 10, colour = "black"),
    axis.text        = element_text(size = 9,  colour = "black"),
    strip.background = element_rect(fill = "grey92", colour = NA),
    strip.text       = element_text(face = "bold", size = 9, colour = "black"),
    legend.title     = element_text(size = 9, colour = "black"),
    legend.text      = element_text(size = 8, colour = "black"),
    legend.position  = "bottom",
    plot.title       = element_blank(),
    plot.subtitle    = element_blank(),
    plot.caption     = element_blank()
  )

SITE_LEVELS <- c("Rupar", "Kalibangan", "Lothal")
site_colors <- c("Rupar" = "#1b7837", "Kalibangan" = "#762a83", "Lothal" = "#d6604d")
site_shapes <- c("Rupar" = 16, "Kalibangan" = 17, "Lothal" = 15)
bone_shapes <- setNames(c(15, 16, 17, 18, 3, 8), bones)

panel_tag_labeller <- function(values) {
  ggplot2::as_labeller(setNames(paste0("(", letters[seq_along(values)], ") ", values), values))
}

save_figure <- function(plot, fignum, width, height, dpi = DPI_COMBINATION_ART) {
  base <- file.path("plots_22.07.2026", paste0("Fig", fignum))
  ggsave(paste0(base, ".eps"), plot, width = width, height = height, device = grDevices::cairo_ps)
  ggsave(paste0(base, ".tiff"), plot, width = width, height = height, dpi = dpi, compression = "lzw")
  invisible(NULL)
}

figure_legends <- list()

frag_long <- df %>%
  select(ID, SITE, ends_with("_FRAG")) %>%
  pivot_longer(cols = ends_with("_FRAG"), names_to = "Bone", values_to = "Frag_cm") %>%
  mutate(
    Bone = str_replace(Bone, "_FRAG", ""),
    Bone = factor(Bone, levels = bones),
    SITE = factor(SITE, levels = SITE_LEVELS)
  ) %>%
  filter(!is.na(Frag_cm))

p1 <- ggplot(frag_long, aes(x = SITE, y = Frag_cm, fill = SITE, colour = SITE)) +
  geom_violin(alpha = 0.3, trim = FALSE) +
  geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.7) +
  geom_jitter(aes(shape = SITE), width = 0.07, size = 1.2, alpha = 0.7) +
  facet_wrap(~ Bone, scales = "free_y", nrow = 2, labeller = panel_tag_labeller(bones)) +
  scale_fill_manual(values = site_colors) +
  scale_colour_manual(values = site_colors) +
  scale_shape_manual(values = site_shapes) +
  labs(x = "Site", y = "Fragment length (cm)", fill = "Site", colour = "Site", shape = "Site") +
  theme_ivc
save_figure(p1, 1, width = 12, height = 8, dpi = DPI_COMBINATION_ART)
figure_legends[["1"]] <- paste(
  "Distribution of measured fragment lengths (cm) for the six long bones",
  "across the three Indus Valley Civilization sites; panels (a-f) = humerus,",
  "radius, ulna, femur, tibia, fibula. Violin width represents observation",
  "density at each length value; boxes show the interquartile range (IQR)",
  "with the median as a horizontal line; whiskers extend to 1.5 x IQR;",
  "points are individual observations, jittered for visibility and coded by",
  "both colour and shape (circle = Rupar, triangle = Kalibangan, square =",
  "Lothal)"
)

pred_long <- df %>%
  select(ID, SITE, ends_with("_FRAG"), ends_with("_PRED_LENGTH")) %>%
  pivot_longer(cols = ends_with("_FRAG"), names_to = "Bone", values_to = "Frag_cm") %>%
  mutate(Bone = str_replace(Bone, "_FRAG", "")) %>%
  left_join(
    df %>%
      select(ID, SITE, ends_with("_PRED_LENGTH")) %>%
      pivot_longer(ends_with("_PRED_LENGTH"), names_to = "Bone", values_to = "Pred_Length") %>%
      mutate(Bone = str_replace(Bone, "_PRED_LENGTH", "")),
    by = c("ID", "SITE", "Bone")
  ) %>%
  mutate(
    Bone = factor(Bone, levels = bones),
    SITE = factor(SITE, levels = SITE_LEVELS)
  ) %>%
  filter(!is.na(Frag_cm))

p2 <- ggplot(pred_long, aes(x = Frag_cm, y = Pred_Length, colour = SITE, shape = SITE)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE, colour = "black", fill = "grey80", linewidth = 0.8) +
  facet_wrap(~ Bone, scales = "free", nrow = 2, labeller = panel_tag_labeller(bones)) +
  scale_colour_manual(values = site_colors) +
  scale_shape_manual(values = site_shapes) +
  labs(x = "Fragment length (cm)", y = "Predicted complete bone length (cm)", colour = "Site", shape = "Site") +
  theme_ivc
save_figure(p2, 2, width = 12, height = 8, dpi = DPI_COMBINATION_ART)
figure_legends[["2"]] <- paste(
  "Measured fragment length (x-axis, cm) plotted against Steele (1970)",
  "predicted complete bone length (y-axis, cm) for the six long bones;",
  "panels (a-f) = humerus, radius, ulna, femur, tibia, fibula. Points are",
  "individual observations coded by site colour and shape (circle = Rupar,",
  "triangle = Kalibangan, square = Lothal); the black line and grey band",
  "show the pooled linear regression fit and its 95% confidence interval"
)

p3_data <- tab6_long %>%
  left_join(select(df_wt, ID, Stature_wt), by = "ID") %>%
  mutate(
    Bone = factor(Bone, levels = bones),
    SITE = factor(SITE, levels = SITE_LEVELS)
  ) %>%
  filter(!is.na(Stature_cm))

p3 <- ggplot(p3_data, aes(x = reorder(ID, Stature_wt), y = Stature_cm, colour = Bone, shape = Bone, group = Bone)) +
  geom_point(size = 1.8, position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), width = 0.2, position = position_dodge(width = 0.5), linewidth = 0.5) +
  facet_wrap(~ SITE, scales = "free_x", nrow = 3, labeller = panel_tag_labeller(SITE_LEVELS)) +
  scale_shape_manual(values = bone_shapes) +
  labs(x = "Individual ID", y = "Estimated stature (cm)", colour = "Bone", shape = "Bone") +
  theme_ivc +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 7))
save_figure(p3, 3, width = 12, height = 12, dpi = DPI_COMBINATION_ART)
figure_legends[["3"]] <- paste(
  "Bone-specific estimated stature (cm) for every individual with at least",
  "one measurable long bone, grouped by site; panels (a-c) = Rupar,",
  "Kalibangan, Lothal. Individuals on the x-axis are ordered by ascending",
  "weighted mean stature; points are coded by bone colour and shape, and",
  "error bars show the 95% confidence interval (+/-1.96 x combined",
  "standard error of estimate)"
)

p4_data <- df_wt %>% filter(!is.na(Stature_wt)) %>% mutate(SITE = factor(SITE, levels = SITE_LEVELS))
p4 <- ggplot(p4_data, aes(x = SITE, y = Stature_wt, fill = SITE)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 2, alpha = 0.7, width = 0.5) +
  geom_jitter(aes(colour = SITE, shape = SITE), width = 0.12, size = 2, alpha = 0.7) +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "white", colour = "black") +
  ggpubr::stat_compare_means(
    comparisons = list(c("Rupar", "Kalibangan"), c("Rupar", "Lothal"), c("Kalibangan", "Lothal")),
    method = "wilcox.test", label = "p.signif", tip.length = 0.01
  ) +
  scale_fill_manual(values = site_colors) +
  scale_colour_manual(values = site_colors) +
  scale_shape_manual(values = site_shapes) +
  labs(x = "Site", y = "Weighted mean stature (cm)", fill = "Site", colour = "Site", shape = "Site") +
  theme_ivc
save_figure(p4, 4, width = 8, height = 7, dpi = DPI_COMBINATION_ART)
figure_legends[["4"]] <- paste(
  "Inverse-variance weighted mean stature (cm) for Rupar, Kalibangan, and",
  "Lothal. Boxes show the median (horizontal line) and interquartile range;",
  "whiskers extend to 1.5 x IQR; the open diamond marks the group mean;",
  "jittered points (circle = Rupar, triangle = Kalibangan, square = Lothal)",
  "show individual estimates. Brackets show pairwise Wilcoxon test",
  "significance (ns = p > 0.05, * = p < 0.05, ** = p < 0.01); all estimates",
  "apply male Trotter and Gleser (1952) constants because no sex data were",
  "available"
)

tab11_plot <- df_wt %>%
  filter(!is.na(Stature_Class)) %>%
  count(SITE, Stature_Class) %>%
  group_by(SITE) %>%
  mutate(
    Pct = n / sum(n) * 100,
    Stature_Class = factor(Stature_Class, levels = c(
      "Very Short (<150 cm)", "Short (150\u2013159.9 cm)",
      "Below Average (160\u2013163.9 cm)", "Average (164\u2013166.9 cm)",
      "Above Average (167\u2013169.9 cm)", "Tall (170\u2013179.9 cm)",
      "Very Tall (\u2265180 cm)"
    )),
    SITE = factor(SITE, levels = SITE_LEVELS)
  )

p5 <- ggplot(tab11_plot, aes(x = SITE, y = Pct, fill = Stature_Class)) +
  geom_col(position = "stack", colour = "white", linewidth = 0.3) +
  scale_fill_viridis_d(option = "viridis", direction = -1) +
  scale_y_continuous(labels = scales::label_percent(scale = 1)) +
  labs(x = "Site", y = "Percentage of individuals (%)", fill = "Stature category") +
  theme_ivc +
  theme(legend.text = element_text(size = 7))
save_figure(p5, 5, width = 9, height = 7, dpi = DPI_COMBINATION_ART)
figure_legends[["5"]] <- paste(
  "Proportion of individuals (%) in each Martin and Saller (1956) stature",
  "category at Rupar, Kalibangan, and Lothal (unequal site sample sizes;",
  "Rupar n = 15, Kalibangan n = 12, Lothal n = 19); stacked bars are coloured",
  "by stature category using a perceptually uniform scale from short (dark",
  "purple) to tall (yellow)"
)

df_indices_plot <- df_indices %>% mutate(
  Bone = factor(Bone, levels = bones),
  SITE = factor(SITE, levels = SITE_LEVELS)
)

p6 <- ggplot(df_indices_plot, aes(x = SITE, y = RI, fill = SITE)) +
  geom_boxplot(outlier.shape = 21, alpha = 0.7, width = 0.5) +
  geom_jitter(aes(colour = SITE, shape = SITE), width = 0.1, size = 1.5, alpha = 0.6) +
  facet_wrap(~ Bone, scales = "free_y", nrow = 2, labeller = panel_tag_labeller(bones)) +
  scale_fill_manual(values = site_colors) +
  scale_colour_manual(values = site_colors) +
  scale_shape_manual(values = site_shapes) +
  labs(x = "Site", y = "Robusticity Index (RI)", fill = "Site", colour = "Site", shape = "Site") +
  theme_ivc
save_figure(p6, 6, width = 12, height = 8, dpi = DPI_COMBINATION_ART)
figure_legends[["6"]] <- paste(
  "Robusticity Index (RI; mean of mediolateral and anteroposterior diameters",
  "divided by predicted bone length, x 100) by long bone and site; panels",
  "(a-f) = humerus, radius, ulna, femur, tibia, fibula. Boxes show the",
  "median and interquartile range, whiskers extend to 1.5 x IQR, and",
  "jittered points show individual observations coded by site colour and",
  "shape (circle = Rupar, triangle = Kalibangan, square = Lothal)"
)

platymeria_threshold  <- 85
platycnemia_threshold <- 70
p7 <- ggplot(df_indices_plot, aes(x = SITE, y = CI, fill = SITE)) +
  geom_boxplot(outlier.shape = 21, alpha = 0.7, width = 0.5) +
  geom_jitter(aes(colour = SITE, shape = SITE), width = 0.1, size = 1.5, alpha = 0.6) +
  geom_hline(
    data = data.frame(
      Bone = factor(c("FEMUR", "TIBIA"), levels = bones),
      threshold = c(platymeria_threshold, platycnemia_threshold)
    ),
    aes(yintercept = threshold), linetype = "dashed", colour = "firebrick", linewidth = 0.7
  ) +
  facet_wrap(~ Bone, scales = "free_y", nrow = 2, labeller = panel_tag_labeller(bones)) +
  scale_fill_manual(values = site_colors) +
  scale_colour_manual(values = site_colors) +
  scale_shape_manual(values = site_shapes) +
  labs(x = "Site", y = "Circularity Index (CI)", fill = "Site", colour = "Site", shape = "Site") +
  theme_ivc
save_figure(p7, 7, width = 12, height = 8, dpi = DPI_COMBINATION_ART)
figure_legends[["7"]] <- paste(
  "Circularity Index (CI; minor diameter / major diameter x 100) by long",
  "bone and site; panels (a-f) = humerus, radius, ulna, femur, tibia,",
  "fibula. Dashed red reference lines mark the platymeria threshold for the",
  "femur (CI = 85) and the platycnemia threshold for the tibia (CI = 70);",
  "boxes, whiskers, and jittered points follow the same conventions as",
  "Fig. 6, with site coded by colour and point shape"
)

corr_vars <- df %>%
  select(ends_with("_FRAG"), ends_with("_ML"), ends_with("_AP")) %>%
  bind_cols(select(df_wt, Stature_wt))
cor_mat <- cor(corr_vars, use = "pairwise.complete.obs")

grDevices::cairo_ps(file.path("plots_22.07.2026", "Fig8.eps"), width = 10, height = 10, family = FIG_FONT)
par(family = FIG_FONT)
corrplot::corrplot(cor_mat, method = "ellipse", type = "lower", tl.col = "black", tl.cex = 0.65, cl.cex = 0.8, mar = c(0, 0, 0, 0))
dev.off()

grDevices::tiff(file.path("plots_22.07.2026", "Fig8.tiff"), width = 10, height = 10, units = "in", res = DPI_COMBINATION_ART, compression = "lzw", family = FIG_FONT)
par(family = FIG_FONT)
corrplot::corrplot(cor_mat, method = "ellipse", type = "lower", tl.col = "black", tl.cex = 0.65, cl.cex = 0.8, mar = c(0, 0, 0, 0))
dev.off()

figure_legends[["8"]] <- paste(
  "Pearson correlation matrix of fragment length, mediolateral diameter,",
  "and anteroposterior diameter for the six long bones and inverse-variance",
  "weighted mean stature. Ellipse orientation and colour intensity scale",
  "with the correlation coefficient (blue = positive, red = negative), and",
  "variable names are given on the matrix diagonal"
)

ri_ci_indiv <- df_indices %>%
  group_by(ID, SITE) %>%
  summarise(Mean_RI = mean(RI, na.rm = TRUE), Mean_CI = mean(CI, na.rm = TRUE), .groups = "drop") %>%
  left_join(select(df_wt, ID, N_bones_used, Stature_wt), by = "ID") %>%
  mutate(SITE = factor(SITE, levels = SITE_LEVELS))

p9 <- ggplot(ri_ci_indiv, aes(x = Mean_RI, y = Mean_CI, colour = SITE, shape = SITE, size = N_bones_used)) +
  geom_point(alpha = 0.75) +
  scale_colour_manual(values = site_colors) +
  scale_shape_manual(values = site_shapes) +
  scale_size_continuous(range = c(2, 8), name = "N bones used") +
  labs(x = "Mean Robusticity Index (RI)", y = "Mean Circularity Index (CI)", colour = "Site", shape = "Site") +
  theme_ivc
save_figure(p9, 9, width = 8, height = 7, dpi = DPI_COMBINATION_ART)
figure_legends[["9"]] <- paste(
  "Bivariate skeletal morphospace plot of mean Robusticity Index (x-axis)",
  "against mean Circularity Index (y-axis) per individual, coded by site",
  "colour and shape (circle = Rupar, triangle = Kalibangan, square =",
  "Lothal); point size is proportional to the number of long bones",
  "contributing to that individual's weighted stature estimate"
)

lm_fem <- lm(FEMUR_PRED_LENGTH ~ FEMUR_FRAG, data = df, na.action = na.omit)
draw_fig10 <- function() {
  par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1.5), oma = c(0, 0, 0, 0), family = FIG_FONT)
  for (w in 1:4) {
    plot(lm_fem, which = w, sub.caption = "", caption = "")
    mtext(letters[w], side = 3, adj = 0, line = 0.5, font = 2, cex = 1.1)
  }
}

grDevices::cairo_ps(file.path("plots_22.07.2026", "Fig10.eps"), width = 9, height = 8, family = FIG_FONT)
draw_fig10()
dev.off()

grDevices::tiff(file.path("plots_22.07.2026", "Fig10.tiff"), width = 9, height = 8, units = "in", res = DPI_LINE_ART, compression = "lzw", family = FIG_FONT)
draw_fig10()
dev.off()

figure_legends[["10"]] <- paste(
  "Femur Step 1 regression diagnostics (fragment length predicting complete",
  "bone length): (a) residuals versus fitted values, (b) normal",
  "quantile-quantile plot of residuals, (c) scale-location plot, (d) Cook's",
  "distance with the 4/n influence threshold"
)

fig_legend_lines <- c(
  "FIGURE LEGENDS",
  "(Fig. N in bold, no punctuation after the number or at the end of the caption, per AAS Figure Captions guidance)",
  "",
  vapply(names(figure_legends), function(n) paste0("Fig. ", n, "  ", figure_legends[[n]]), character(1))
)
writeLines(fig_legend_lines, file.path("reports_22.07.2026", "Figure_Legends.txt"))

# ── 11. EXPORT MASTER RESULTS TABLE ──────────────────────────────────────────
master <- df_wt %>%
  select(ID, SITE, SIDE, N_bones_used, Stature_wt, Stature_SE, CI95_lower, CI95_upper, Stature_Class, all_of(stat_cols)) %>%
  left_join(ri_ci_indiv %>% select(ID, Mean_RI, Mean_CI), by = "ID") %>%
  rename_with(~ str_replace(., "_STATURE", "_Est_Stature_cm"), ends_with("_STATURE")) %>%
  mutate(across(where(is.numeric), ~ round(., 2)))
write_csv(master, "tables_22.07.2026/Table_MASTER_All_Individual_Results.csv")

# ── 12. WORD REPORT (AAS-COMPLIANT) ──────────────────────────────────────────
sig_stars <- function(p) {
  dplyr::case_when(
    is.na(p)  ~ "",
    p < 0.001 ~ "***",
    p < 0.01  ~ "**",
    p < 0.05  ~ "*",
    TRUE      ~ "ns"
  )
}
sig_footnote_text <- "Significance codes: * p < 0.05, ** p < 0.01, *** p < 0.001, ns = not significant"

tab9_display <- normality_res %>%
  mutate(`Sig.` = sig_stars(SW_p)) %>%
  select(SITE, N, SW_W, SW_p, `Sig.`, Normal)

posthoc_display <- posthoc_df %>%
  mutate(`Sig.` = sig_stars(p_adj)) %>%
  select(-Significant) %>%
  relocate(`Sig.`, .after = p_adj)

tab14_display <- tab14 %>%
  mutate(`RI Sig.` = sig_stars(RI_KW_p), `CI Sig.` = sig_stars(CI_KW_p)) %>%
  select(-RI_Sig, -CI_Sig) %>%
  relocate(`RI Sig.`, .after = RI_KW_p) %>%
  relocate(`CI Sig.`, .after = CI_KW_p)

doc <- officer::read_docx()

pick_style <- function(doc, candidates, fallback = "heading 1") {
  avail <- officer::styles_info(doc)$style_name
  hit <- candidates[candidates %in% avail]
  if (length(hit) > 0) hit[1] else fallback
}
title_style <- pick_style(doc, c("Title", "title"))

add_section <- function(doc, heading, description, tbl_df, footnote_col = NULL, footnote_text = NULL) {
  doc <- doc %>%
    officer::body_add_par(heading, style = "heading 1") %>%
    officer::body_add_par(description, style = "Normal") %>%
    officer::body_add_par("", style = "Normal")
  
  is_p_col <- function(nm) grepl("p_adj|p\\.value|^SW_p$|_p$|^p$", nm, ignore.case = TRUE)
  num_cols       <- names(tbl_df)[vapply(tbl_df, is.numeric, logical(1))]
  p_cols         <- num_cols[is_p_col(num_cols)]
  other_num_cols <- setdiff(num_cols, p_cols)
  
  ft_df <- tbl_df %>%
    mutate(across(all_of(p_cols), ~ round(., 4))) %>%
    mutate(across(all_of(other_num_cols), ~ round(., 2)))
  
  ft <- flextable::flextable(ft_df) %>%
    flextable::theme_vanilla() %>%
    flextable::fontsize(size = 8, part = "all") %>%
    flextable::autofit()
  
  if (!is.null(footnote_col) && footnote_col %in% names(tbl_df)) {
    ft <- flextable::footnote(
      ft, i = 1, j = which(names(tbl_df) == footnote_col),
      value = flextable::as_paragraph(footnote_text),
      ref_symbols = "a", part = "header"
    )
  }
  doc <- flextable::body_add_flextable(doc, ft)
  doc <- officer::body_add_par(doc, "", style = "Normal")
  doc
}

doc <- doc %>%
  officer::body_add_par("Stature Estimation from Long Bone Fragments — IVC Sites", style = title_style) %>%
  officer::body_add_par(
    paste("Analysis of stature, robusticity, and circularity from fragmentary",
          "long bones (humerus, radius, ulna, femur, tibia, fibula) recovered from",
          "Rupar, Kalibangan, and Lothal (Indus Valley Civilization). Two-step",
          "regression: Steele (1970) fragment to length; Trotter and Gleser (1952)",
          "length to stature. Inverse-variance weighted mean stature per individual.",
          "Martin and Saller (1956) classification. Kruskal-Wallis inter-site tests."),
    style = "Normal"
  )

tables_list <- list(
  list("Table 1. Sample Composition",
       "Number of individuals and available bone measurements per site.",
       tab1),
  list("Table 2. Raw Measurement Descriptives",
       "Mean, SD, range and CV% for fragment length, ML and AP diameters per bone per site.",
       filter(tab2, Site != "Pooled")),
  list("Table 3. Step 1 Regression Coefficients (Steele 1970)",
       "Intercept, slope and SEE for predicting bone length from fragment length.",
       tab3 %>% select(-Source)),
  list("Table 4. Step 2 Trotter and Gleser Coefficients",
       "Trotter and Gleser (1952) stature equations and combined SEE per bone.",
       tab5 %>% select(-Formula)),
  list("Table 5. Weighted Mean Stature and Classification",
       "Inverse-variance weighted mean stature, 95% CI and Martin and Saller category per individual.",
       tab7),
  list("Table 6. Site-Level RI and CI Summary",
       "Mean and SD robusticity and circularity indices per bone per site.",
       tab9),
  list("Table 7. Stature Descriptives by Site",
       "Mean, median, SD, range and CV% of weighted stature per site.",
       tab10),
  list("Table 8. Stature Category Frequencies",
       "Frequency and percentage of Martin and Saller categories per site.",
       tab11),
  list("Table 9. Normality Tests",
       "Shapiro-Wilk test results per site and Levene's homoscedasticity test.",
       tab9_display, "Sig.", sig_footnote_text),
  list("Table 10. Pairwise Site Comparisons",
       "Pairwise Mann-Whitney U tests with Bonferroni correction.",
       posthoc_display, "Sig.", sig_footnote_text),
  list("Table 11. RI and CI Kruskal-Wallis Tests",
       "Inter-site Kruskal-Wallis tests for robusticity and circularity indices per bone.",
       tab14_display, "RI Sig.", sig_footnote_text)
)

for (tl in tables_list) {
  fcol <- if (length(tl) >= 4) tl[[4]] else NULL
  ftxt <- if (length(tl) >= 5) tl[[5]] else NULL
  doc <- add_section(doc, tl[[1]], tl[[2]], tl[[3]], footnote_col = fcol, footnote_text = ftxt)
}

doc <- doc %>%
  officer::body_add_par("Supplementary Tables", style = "heading 1") %>%
  officer::body_add_par(
    paste("Three per-individual tables are provided as supplementary CSV",
          "files rather than in the main text, so that in-text table",
          "citation remains consecutive: Supplementary Table S1 (predicted",
          "bone lengths per individual; Table4_Predicted_Bone_Lengths.csv),",
          "Supplementary Table S2 (bone-specific stature estimates per",
          "individual; Table6_BoneSpecific_Stature_Estimates.csv), and",
          "Supplementary Table S3 (robusticity and circularity indices per",
          "individual per bone; Table8_Robusticity_Circularity_Indices.csv),",
          "all in tables_22.07.2026/."),
    style = "Normal"
  )

doc <- doc %>% officer::body_add_par("Figure Legends", style = "heading 1")
for (n in as.character(1:10)) {
  doc <- officer::body_add_fpar(
    doc,
    officer::fpar(
      officer::ftext(paste0("Fig. ", n, "  "), officer::fp_text(bold = TRUE)),
      officer::ftext(figure_legends[[n]], officer::fp_text(bold = FALSE))
    )
  )
}

print(doc, target = "reports_22.07.2026/IVC_Stature_Analysis_Report.docx")
cat("\nWord report saved: reports_22.07.2026/IVC_Stature_Analysis_Report.docx\n")

# ── 13. FINAL SUMMARY ────────────────────────────────────────────────────────
cat("\n")
cat("══════════════════════════════════════════════════════════════\n")
cat(" ANALYSIS COMPLETE — OUTPUT SUMMARY (AAS-compliant tables & figures)\n")
cat("══════════════════════════════════════════════════════════════\n")