#!/usr/bin/env Rscript
# ---------------------------------------------------------------------------
# z93_ftdna_register_analysis.R (2026-08-18)
# Independent replication pass for the /r1a-z93-distribution post: re-run the
# temporal and geographic tests on the FTDNA Discover ancient register
# (364 sequence-level re-calls, disjoint pipeline from the AADR v66 audit)
# and test the modern inversion formally on the FTDNA haplotree census.
# Era bins are identical to z93_distribution_analysis.R (BP breaks).
# ---------------------------------------------------------------------------
D <- "/home/cdr/projects/sohi-qpadm/runs/discover-ftdna-2026-08-18"
reg <- read.csv(file.path(D, "ftdna-register-arm-year.csv"), stringsAsFactors = FALSE)
reg$bp <- 1950 - as.numeric(reg$year)
era <- function(bp) cut(bp, breaks = c(-Inf, 1000, 1800, 2800, 3500, 4300, Inf),
                        labels = c("6 Medieval+", "5 Migration-era 150-950CE",
                                   "4 Scythian-Xiongnu", "3 LBA-Iron",
                                   "2 Sintashta-Andronovo", "1 pre-Sintashta"))
reg$era <- factor(era(reg$bp), levels = rev(levels(era(0))))

cat("# FTDNA Discover ancient register — replication analysis (2026-08-18)\n\n")
cat("Register: ", nrow(reg), " ancient Z93 males with a usable date, classified into\n")
cat("arms by subtree membership in the FTDNA haplotree (Y3 = below F2597;\n")
cat("Z2124 = below Z2124; else basal/other). Pipeline disjoint from AADR v66.\n\n")

cat("## 1. Era x arm counts (register)\n\n")
t1 <- table(reg$era, reg$arm)
print(t1)
share <- prop.table(t1, 1)[, "Z2124 arm"]
cat("\nZ2124-arm share by era:\n"); print(round(share, 3))
cat("\nAADR v66 audit shares for comparison (published): 0 / .51 / .54 / .67 / .70 / —\n\n")

cat("## 2. Temporal trend, register only\n\n")
reg$z2124 <- as.integer(reg$arm == "Z2124 arm")
reg$kyr <- reg$bp / 1000
m <- glm(z2124 ~ kyr, data = reg, family = binomial)
co <- summary(m)$coefficients
cat("logistic Z2124 ~ kyr-BP: coef =", round(co[2, 1], 3),
    " p =", signif(co[2, 4], 3), "\n")
cat("(negative coef = Z2124 share RISES toward the present; same caveats as the\n")
cat("audit: cemetery patrilines violate independence, era bins are choices.)\n\n")

cat("## 3. Y3-arm ancients in the register\n\n")
y3 <- reg[reg$arm == "Y3 arm", c("code", "year", "country", "study")]
y3 <- y3[order(as.numeric(y3$year)), ]
print(y3, row.names = FALSE)
cat("\nPre-CE Y3-arm men: ", sum(as.numeric(y3$year) < 0), " of ", nrow(y3), "\n\n")

cat("## 4. Modern inversion test (FTDNA haplotree census)\n\n")
mod <- read.csv(file.path("/home/cdr/projects/sohi-qpadm/runs",
                          "ftdna-modern-country-x-arm-2026-08-18.csv"),
                check.names = FALSE)
pick <- function(cs) colSums(mod[mod$country %in% cs,
                                 c("Y3 arm (toward L657)", "Z2124 arm")])
sa <- pick(c("India", "Pakistan", "Sri Lanka", "Bangladesh"))
ee <- pick(c("Ukraine", "Poland", "Lithuania", "Belarus", "Russia"))
tab <- rbind(SouthAsia = sa, EasternEurope = ee)
print(tab)
ft <- fisher.test(tab)
cat("\nFisher exact, South Asia vs Eastern Europe, Y3 vs Z2124:\n")
cat("odds ratio =", signif(ft$estimate, 4), " p =", signif(ft$p.value, 3), "\n")
cat("(self-reported tester countries; ~63% of kits report no country)\n")