################################################################################
#                                                                              #
#  HARMONIC SUSCEPTIBILITY INDEX (HSI) - REPLICATION SCRIPT                     #
#                                                                              #
#  This single, self-contained script reproduces the core numerical results    #
#  reported in the paper "The Harmonic Susceptibility Index" from the deposited #
#  analytical datasets. It is a faithful reduction of the full research         #
#  pipeline to the steps that generate the published figures.                   #
#                                                                              #
#  INPUTS (provided alongside this script):                                     #
#    - data_individual_items.csv : individual-level data, 21 normalised EVS-WVS #
#         items + country code (cntry_AN, ISO-2). N = 156,658 respondents.      #
#    - data_country_HOI.csv      : country-level Harmonic Outcome Index (HOI),  #
#         the independent well-being validation measure (92 countries).         #
#    - (optional) external benchmark CSVs for the discriminant-validity,       #
#         external-outcome and spatial sections - see PARTS 6, 10, 11.          #
#    - (optional) data_individual_welzel.csv : 12 recoded emancipative-values   #
#         items (Welzel 2013) + country code, for the PART 12 comparison.       #
#                                                                              #
#  OUTPUT:                                                                       #
#    Console report of all reproduced quantities, each printed next to the      #
#    value reported in the paper, so that agreement can be checked directly.    #
#                                                                              #
#  REQUIREMENTS: base R (>= 4.0). No external packages required.                #
#                                                                              #
#  Run:  Rscript replicate_HSI.R                                                #
#                                                                              #
################################################################################

options(stringsAsFactors = FALSE)

# ------------------------------------------------------------------------------
# 0. Locate the data folder
# ------------------------------------------------------------------------------
# The script reads all data files from a single folder. It first tries the
# current working directory; if the files are not there, it asks the user to
# pick ANY file from the data folder (via a dialog) and uses that folder for all
# files. This works both in automated runs (correct working directory) and in
# interactive sessions (no need to set the working directory by hand).
DATA_DIR <- getwd()

.files_present <- function(dir) {
  all(file.exists(file.path(dir, c("data_individual_items.csv",
                                   "data_country_HOI.csv"))))
}

if (!.files_present(DATA_DIR)) {
  if (interactive()) {
    cat("Data files not found in the working directory.\n")
    cat("Please select ANY file inside the replication folder...\n")
    picked  <- file.choose()              # opens a dialog; pick e.g. the README or a CSV
    DATA_DIR <- dirname(picked)
  }
}

if (!.files_present(DATA_DIR)) {
  stop("Could not find 'data_individual_items.csv' and 'data_country_HOI.csv'. ",
       "Place all files in one folder and either run the script from that ",
       "folder or select a file inside it when prompted.")
}

PATH_ITEMS   <- file.path(DATA_DIR, "data_individual_items.csv")
PATH_COUNTRY <- file.path(DATA_DIR, "data_country_HOI.csv")

cat("\n==============================================================\n")
cat(" HSI REPLICATION - core results\n")
cat("==============================================================\n")
cat(sprintf(" Data folder: %s\n", DATA_DIR))

# ------------------------------------------------------------------------------
# 1. Item-to-dimension map (S.P.A. taxonomy, 7 dimensions x 3 sub-components)
#    Each dimension is the robust mean of its three items (Spirit/Thought/Action).
#    Inverse-polarity items already carry the "_inv" tag in the deposited data:
#    every item is oriented so that higher values = higher predisposition.
# ------------------------------------------------------------------------------
triples <- list(
  D1_Objectives   = c("A041_norm", "Y002_norm", "A080_01_norm"),
  D2_Origin       = c("A006_inv_norm", "Etica_inv_norm", "E025_norm"),
  D3_Quality      = c("B008_norm", "D3_Qualita_Pensiero_norm", "E026_norm"),
  D4_Context      = c("A165_norm", "A002_inv_norm", "A066_dic"),
  D5_Method       = c("A035_norm", "E039_norm", "A080_02_norm"),
  D6_Organisation = c("A042_inv_norm", "E235_norm", "E027_norm"),
  D7_Approach     = c("E037_inv_norm", "A032_norm", "A029_norm")
)
dim_names <- names(triples)
# Pole partition (D3 Quality is transversal and excluded from the poles)
systemic_pole  <- c("D1_Objectives", "D2_Origin", "D5_Method")
relational_pole <- c("D4_Context", "D6_Organisation", "D7_Approach")

# Certified penalty weights (lambda), paper Section 3.4
LAMBDA <- c(diss7 = 0.60, weak7 = 0.35, blockA = 0.20, blockB = 0.20, d3 = 0.25)

# ------------------------------------------------------------------------------
# 2. Load individual-level items and build the 7 dimensions per respondent
#    Rule "at least 2 of 3": a dimension is NA if fewer than 2 items are present.
# ------------------------------------------------------------------------------
di <- read.csv(PATH_ITEMS, check.names = FALSE, fileEncoding = "UTF-8")
cat(sprintf("\nLoaded individual data: %d respondents, %d items\n",
            nrow(di), ncol(di) - 1))

build_dimension <- function(df, items) {
  M  <- as.matrix(df[, items, drop = FALSE])
  vn <- rowSums(!is.na(M))
  mv <- rowMeans(M, na.rm = TRUE)
  mv[vn < 2] <- NA
  mv
}
dim_ind <- as.data.frame(lapply(triples, build_dimension, df = di))
dim_ind$cntry <- di$cntry_AN

# ------------------------------------------------------------------------------
# 3. Aggregate to country level (mean over respondents), keep complete profiles
# ------------------------------------------------------------------------------
agg <- aggregate(. ~ cntry, data = dim_ind, FUN = function(x) mean(x, na.rm = TRUE))
agg <- agg[stats::complete.cases(agg[, dim_names]), ]
cat(sprintf("Countries with complete 7-dimension profile: %d\n", nrow(agg)))

# ------------------------------------------------------------------------------
# 4. HSI formula (certified, paper Section 3.4)
#    HSI_raw = mean(7D)
#             - lambda_diss7 * sd(7D)
#             - lambda_weak7 * (1 - min(7D))
#             - lambda_blockA * sd(systemic pole)
#             - lambda_blockB * sd(relational pole)
#             - lambda_d3    * (1 - D3)
#    then min-max normalised to [0,1] over the sample.
# ------------------------------------------------------------------------------
compute_hsi_raw <- function(D) {
  X  <- as.matrix(D[, dim_names])
  mean7 <- rowMeans(X)
  diss7 <- apply(X, 1, stats::sd)
  weak7 <- apply(X, 1, min)
  dissA <- apply(D[, systemic_pole],  1, stats::sd)
  dissB <- apply(D[, relational_pole], 1, stats::sd)
  d3    <- D[["D3_Quality"]]
  mean7 - LAMBDA["diss7"] * diss7 -
          LAMBDA["weak7"] * (1 - weak7) -
          LAMBDA["blockA"] * dissA -
          LAMBDA["blockB"] * dissB -
          LAMBDA["d3"] * (1 - d3)
}
minmax <- function(x) (x - min(x, na.rm = TRUE)) / diff(range(x, na.rm = TRUE))

agg$Integral_Score <- rowMeans(as.matrix(agg[, dim_names]))   # S (non-penalised mean)
agg$HSI_raw <- as.numeric(compute_hsi_raw(agg))
agg$HSI     <- minmax(agg$HSI_raw)
agg$rank    <- rank(-agg$HSI, ties.method = "min")

# ------------------------------------------------------------------------------
# 5. Harmonic classes (FIXED cut-offs on the min-max-normalised HSI, paper 4.1:
#    High >= 0.75, Good [0.50, 0.75), Medium [0.25, 0.50), Low < 0.25.
#    Applied to the observed distribution this yields 12 / 29 / 42 / 9 countries)
#    and quadrants from sample-standardised poles.
# ------------------------------------------------------------------------------
classify_hsi <- function(x) {
  as.character(cut(x, breaks = c(-Inf, 0.25, 0.50, 0.75, Inf),
                   labels = c("Low", "Medium", "Good", "High"), right = FALSE))
}
agg$class <- classify_hsi(agg$HSI)

z <- function(x) (x - mean(x)) / stats::sd(x)
PS <- rowMeans(scale(agg[, systemic_pole]))
PR <- rowMeans(scale(agg[, relational_pole]))
agg$quadrant <- ifelse(PS > 0 & PR > 0, "Q1",
                ifelse(PS > 0 & PR < 0, "Q2",
                ifelse(PS < 0 & PR < 0, "Q3", "Q4")))

# ------------------------------------------------------------------------------
# 6. External validation against the HOI (independent well-being measure)
# ------------------------------------------------------------------------------
cl <- read.csv(PATH_COUNTRY, check.names = FALSE, fileEncoding = "UTF-8")
m  <- merge(agg[, c("cntry", "HSI", "Integral_Score")],
            cl[, c("iso2", "HOI")], by.x = "cntry", by.y = "iso2")

cor_hsi_hoi  <- cor(m$HSI, m$HOI)
fit_hoi_hsi  <- lm(HOI ~ HSI, data = m)
r2_hsi       <- summary(fit_hoi_hsi)$r.squared
fit_hoi_s    <- lm(HOI ~ Integral_Score, data = m)
r2_s         <- summary(fit_hoi_s)$r.squared

# ------------------------------------------------------------------------------
# 6b. Discriminant validity (5 external non-survey benchmarks) and
#     convergent validation against the World Happiness Report life ladder.
#     External files are optional: if absent, this part is skipped.
# ------------------------------------------------------------------------------
# ISO-2 (HSI data) -> ISO-3 (benchmark files) lookup, restricted to the
# countries present in the analysis (same crosswalk logic as the full pipeline).
iso2_to_iso3 <- c(
  AD="AND",AL="ALB",AM="ARM",AR="ARG",AT="AUT",AU="AUS",AZ="AZE",BA="BIH",BD="BGD",
  BG="BGR",BO="BOL",BR="BRA",BY="BLR",CA="CAN",CH="CHE",CL="CHL",CN="CHN",CO="COL",
  CY="CYP",CZ="CZE",DE="DEU",DK="DNK",EC="ECU",EE="EST",EG="EGY",ES="ESP",ET="ETH",
  FI="FIN",FR="FRA",GB="GBR",GE="GEO",GR="GRC",GT="GTM",HK="HKG",HR="HRV",HU="HUN",
  ID="IDN",IN="IND",IQ="IRQ",IR="IRN",IS="ISL",IT="ITA",JO="JOR",JP="JPN",KE="KEN",
  KG="KGZ",KR="KOR",KZ="KAZ",LB="LBN",LT="LTU",LV="LVA",LY="LBY",MA="MAR",ME="MNE",
  MK="MKD",MM="MMR",MN="MNG",MO="MAC",MV="MDV",MX="MEX",MY="MYS",NG="NGA",NI="NIC",
  NIR="GBR",NL="NLD",NO="NOR",NZ="NZL",PE="PER",PH="PHL",PK="PAK",PL="POL",PR="PRI",
  PT="PRT",RO="ROU",RS="SRB",RU="RUS",SE="SWE",SG="SGP",SI="SVN",SK="SVK",TH="THA",
  TJ="TJK",TN="TUN",TR="TUR",TW="TWN",UA="UKR",US="USA",UY="URY",UZ="UZB",VE="VEN",
  VN="VNM",ZW="ZWE"
)
agg$iso3 <- iso2_to_iso3[agg$cntry]

# Some HSI units share an ISO-3 code (e.g. Great Britain 'GB' and Northern
# Ireland 'NIR' both map to 'GBR'). External benchmark files carry a single
# row per ISO-3, so before merging we collapse HSI to one row per ISO-3,
# averaging the HSI value. This reproduces the country coverage of the full
# pipeline (e.g. 86 countries matched for the WHR life ladder) and prevents
# duplicate rows from inflating the sample.
agg_iso3 <- aggregate(HSI ~ iso3, data = agg, FUN = mean)

# Robust CSV reader: detect separator (',' or ';'), handle comma decimals.
read_benchmark <- function(path, score_col) {
  first <- readLines(path, n = 1, warn = FALSE, encoding = "latin1")
  sep <- if (lengths(regmatches(first, gregexpr(";", first))) >
             lengths(regmatches(first, gregexpr(",", first)))) ";" else ","
  # latin1 + quote='' avoids truncation on stray bytes / unbalanced quotes in
  # free-text note columns, which can otherwise abort reading mid-file.
  df <- read.csv(path, sep = sep, check.names = FALSE,
                 encoding = "latin1", quote = "\"", fill = TRUE)
  names(df) <- sub("^\ufeff", "", names(df))          # strip BOM from 1st name
  names(df) <- trimws(names(df))
  sc <- df[[score_col]]
  if (is.character(sc)) sc <- as.numeric(gsub(",", ".", sc))  # comma decimals
  data.frame(iso3 = toupper(trimws(as.character(df[["iso3"]]))), score = sc)
}

benchmarks <- list(
  HDI      = list(file = "HDI_2023.csv",          col = "HDI_score_2023"),
  LPI      = list(file = "LPI_2023.csv",          col = "lpi_score_2023"),
  TEHA_GAI = list(file = "TEHA_GAI_2025.csv",      col = "gai_score_2025"),
  WGI      = list(file = "WGI_2024.csv",           col = "wgi_composite_score_2024"),
  GII      = list(file = "GII_2025.csv",           col = "gii_score_2025")
)
PAPER_DR2 <- c(HDI = 0.273, LPI = 0.242, TEHA_GAI = 0.282, WGI = 0.265, GII = 0.309)

disc_results <- list()
have_benchmarks <- all(sapply(benchmarks, function(b) file.exists(file.path(DATA_DIR, b$file))))
if (have_benchmarks) {
  # HSI (one row per ISO-3) + HOI, joined via the ISO-2 country code
  hsi_hoi <- merge(agg[, c("cntry", "iso3", "HSI")],
                   cl[, c("iso2", "HOI")], by.x = "cntry", by.y = "iso2")
  # collapse to one row per ISO-3 (HSI and HOI averaged where codes coincide)
  hsi_hoi <- aggregate(cbind(HSI, HOI) ~ iso3, data = hsi_hoi, FUN = mean)
  for (bid in names(benchmarks)) {
    b <- benchmarks[[bid]]
    bench <- read_benchmark(file.path(DATA_DIR, b$file), b$col)
    d <- merge(hsi_hoi, bench, by = "iso3")
    d <- d[is.finite(d$HSI) & is.finite(d$HOI) & is.finite(d$score), ]
    r2_bench <- summary(lm(HOI ~ score, data = d))$r.squared
    r2_full  <- summary(lm(HOI ~ score + HSI, data = d))$r.squared
    disc_results[[bid]] <- c(n = nrow(d), dR2 = r2_full - r2_bench)
  }
}

# World Happiness Report life ladder (convergent, independent source)
whr_r <- NA_real_; whr_n <- NA_integer_
if (file.exists(file.path(DATA_DIR, "WHR_LifeLadder_2020.csv"))) {
  whr <- read_benchmark(file.path(DATA_DIR, "WHR_LifeLadder_2020.csv"), "life_ladder_2020")
  dw <- merge(agg_iso3, whr, by = "iso3")
  dw <- dw[is.finite(dw$HSI) & is.finite(dw$score), ]
  whr_r <- cor(dw$HSI, dw$score); whr_n <- nrow(dw)
}

# ------------------------------------------------------------------------------
# 7. Robustness (paper Section 4.5)
#    (a) variance-based dimensional weights vs equal weighting
#    (b) geometric vs arithmetic aggregation
#    (c) joint Monte Carlo perturbation of the five lambda (+/-50%)
# ------------------------------------------------------------------------------
hsi_with_weights <- function(D, w = rep(1, 7)) {
  X <- as.matrix(D[, dim_names]); w <- w / sum(w)
  mean7 <- as.numeric(X %*% w)
  diss7 <- apply(X, 1, stats::sd); weak7 <- apply(X, 1, min)
  dissA <- apply(D[, systemic_pole], 1, stats::sd)
  dissB <- apply(D[, relational_pole], 1, stats::sd); d3 <- D[["D3_Quality"]]
  raw <- mean7 - LAMBDA["diss7"]*diss7 - LAMBDA["weak7"]*(1-weak7) -
         LAMBDA["blockA"]*dissA - LAMBDA["blockB"]*dissB - LAMBDA["d3"]*(1-d3)
  minmax(raw)
}
hsi_geometric <- function(D) {
  X <- pmax(as.matrix(D[, dim_names]), 1e-9)
  geo <- exp(rowMeans(log(X)))
  diss7 <- apply(X, 1, stats::sd); weak7 <- apply(X, 1, min)
  dissA <- apply(D[, systemic_pole], 1, stats::sd)
  dissB <- apply(D[, relational_pole], 1, stats::sd); d3 <- D[["D3_Quality"]]
  raw <- geo - LAMBDA["diss7"]*diss7 - LAMBDA["weak7"]*(1-weak7) -
         LAMBDA["blockA"]*dissA - LAMBDA["blockB"]*dissB - LAMBDA["d3"]*(1-d3)
  minmax(raw)
}
var_w   <- apply(agg[, dim_names], 2, var)
rho_varw <- cor(hsi_with_weights(agg, var_w), agg$HSI, method = "spearman")
rho_geom <- cor(hsi_geometric(agg),           agg$HSI, method = "spearman")

set.seed(20260525)
NB <- 5000
rho_mc <- numeric(NB)
top12_base <- order(-agg$HSI)[1:12]
stab12 <- numeric(NB)
for (b in seq_len(NB)) {
  lam_b <- LAMBDA * runif(5, 0.5, 1.5); names(lam_b) <- names(LAMBDA)
  X <- as.matrix(agg[, dim_names])
  raw <- rowMeans(X) - lam_b["diss7"]*apply(X,1,stats::sd) -
         lam_b["weak7"]*(1-apply(X,1,min)) -
         lam_b["blockA"]*apply(agg[,systemic_pole],1,stats::sd) -
         lam_b["blockB"]*apply(agg[,relational_pole],1,stats::sd) -
         lam_b["d3"]*(1-agg[["D3_Quality"]])
  hsi_b <- minmax(raw)
  rho_mc[b] <- cor(hsi_b, agg$HSI, method = "spearman")
  stab12[b] <- length(intersect(order(-hsi_b)[1:12], top12_base)) / 12
}
rho_mc_med <- median(rho_mc)
rho_mc_ci  <- quantile(rho_mc, c(0.05, 0.95))
stab12_mean <- mean(stab12)

# ------------------------------------------------------------------------------
# 8. Zero-penalty comparison (paper Section 4.5, Fig. 5)
#    HSI vs the index obtained with all five lambda set to zero (the min-max-
#    normalised Integral Score). Same fixed class cut-offs applied to both.
# ------------------------------------------------------------------------------
zp <- minmax(agg$Integral_Score)
rho_zp <- cor(agg$HSI, zp, method = "spearman")
cls_lv <- c(Low = 1, Medium = 2, Good = 3, High = 4)
zp_class  <- classify_hsi(zp)
n_change  <- sum(agg$class != zp_class)
# Direction of movement, stated on both metrics for full transparency:
#  - CLASS: under the penalised HSI vs the zero-penalty index (higher / lower)
#  - RANK (among the class-changers): rank improves / worsens once penalised
n_cls_up   <- sum(cls_lv[agg$class] > cls_lv[zp_class])
n_cls_down <- sum(cls_lv[agg$class] < cls_lv[zp_class])
rk_h <- rank(-agg$HSI, ties.method = "min"); rk_z <- rank(-zp, ties.method = "min")
chg  <- agg$class != zp_class
n_rk_impr <- sum(chg & rk_h < rk_z); n_rk_wors <- sum(chg & rk_h > rk_z)
high_demoted <- sort(agg$cntry[agg$class == "High" & zp_class == "Good"])

# ------------------------------------------------------------------------------
# 9. Further robustness reported in the paper (Sections 3.5, 4.2, 4.5)
# ------------------------------------------------------------------------------
# 9a. Individual perturbation of each lambda: baseline x {0.50, 0.75, 1.25, 1.50}
#     (20 specifications). Spearman/Kendall vs baseline; top-10 set; High class.
specs_rho_s <- c(); specs_rho_k <- c(); top10_same <- c(); high_exch25 <- c(); high_exch50 <- c()
top10_base <- agg$cntry[order(-agg$HSI)][1:10]
high_base  <- agg$cntry[order(-agg$HSI)][1:12]   # High class = top 12 by rank
for (li in names(LAMBDA)) for (f in c(0.50, 0.75, 1.25, 1.50)) {
  lam <- LAMBDA; lam[li] <- lam[li] * f
  X <- as.matrix(agg[, dim_names])
  raw <- rowMeans(X) - lam["diss7"]*apply(X,1,stats::sd) -
         lam["weak7"]*(1-apply(X,1,min)) -
         lam["blockA"]*apply(agg[,systemic_pole],1,stats::sd) -
         lam["blockB"]*apply(agg[,relational_pole],1,stats::sd) -
         lam["d3"]*(1-agg[["D3_Quality"]])
  h <- minmax(raw)
  specs_rho_s <- c(specs_rho_s, cor(h, agg$HSI, method = "spearman"))
  specs_rho_k <- c(specs_rho_k, cor(h, agg$HSI, method = "kendall"))
  top10_same  <- c(top10_same, length(setdiff(top10_base, agg$cntry[order(-h)][1:10])))
  # High-class membership by rank (top 12), as in the certified robustness protocol
  n_ex <- length(setdiff(high_base, agg$cntry[order(-h)][1:12]))
  if (f %in% c(0.75, 1.25)) high_exch25 <- c(high_exch25, n_ex) else high_exch50 <- c(high_exch50, n_ex)
}

# 9b. Removal of the most value-laden items (religiosity A006, obedience A042),
#     singly and jointly, plus systematic leave-one-item-out (21 items).
#     Rule as in the pipeline: a dimension needs at least ceiling(2/3 of its
#     items); with a reduced 2-item triple both items are required.
east_asia <- c("JP","KR","CN","TW","HK","MO","SG","MN")
hsi_quad_from_triples <- function(tr) {
  dims <- as.data.frame(lapply(tr, function(items) {
    M <- as.matrix(di[, items, drop = FALSE])
    vn <- rowSums(!is.na(M)); mv <- rowMeans(M, na.rm = TRUE)
    mv[vn < max(1, ceiling(length(items) * 2/3))] <- NA
    mv
  }))
  dims$cntry <- di$cntry_AN
  dims <- dims[stats::complete.cases(dims[, names(tr)]), ]
  a <- aggregate(. ~ cntry, data = dims, FUN = function(x) mean(x, na.rm = TRUE))
  X <- as.matrix(a[, names(tr)])
  raw <- rowMeans(X) - LAMBDA["diss7"]*apply(X,1,stats::sd) -
         LAMBDA["weak7"]*(1-apply(X,1,min)) -
         LAMBDA["blockA"]*apply(a[, systemic_pole],1,stats::sd) -
         LAMBDA["blockB"]*apply(a[, relational_pole],1,stats::sd) -
         LAMBDA["d3"]*(1-a[["D3_Quality"]])
  PSx <- rowMeans(scale(a[, systemic_pole])); PRx <- rowMeans(scale(a[, relational_pole]))
  data.frame(cntry = a$cntry, HSI = minmax(raw),
             Q = ifelse(PSx > 0 & PRx > 0, "Q1", ifelse(PSx > 0 & PRx < 0, "Q2",
                 ifelse(PSx < 0 & PRx < 0, "Q3", "Q4"))), stringsAsFactors = FALSE)
}
rho_vs_base <- function(alt) {
  mm <- merge(agg[, c("cntry","HSI")], alt[, c("cntry","HSI")], by = "cntry")
  cor(mm$HSI.x, mm$HSI.y, method = "spearman")
}
q4_asia <- function(a) sum(a$Q[a$cntry %in% east_asia] == "Q4")
tr_noRel <- triples; tr_noRel$D2_Origin       <- setdiff(triples$D2_Origin, "A006_inv_norm")
tr_noObb <- triples; tr_noObb$D6_Organisation <- setdiff(triples$D6_Organisation, "A042_inv_norm")
tr_no2   <- tr_noRel; tr_no2$D6_Organisation  <- tr_noObb$D6_Organisation
r_noRel <- hsi_quad_from_triples(tr_noRel); r_noObb <- hsi_quad_from_triples(tr_noObb)
r_no2   <- hsi_quad_from_triples(tr_no2)
rho_val_min <- min(rho_vs_base(r_noRel), rho_vs_base(r_noObb), rho_vs_base(r_no2))
loo <- sapply(unlist(triples, use.names = FALSE), function(it) {
  rho_vs_base(hsi_quad_from_triples(lapply(triples, function(v) setdiff(v, it))))
})
ref_full <- hsi_quad_from_triples(triples)

# 9c. Inter-dimensional correlation structure (paper 4.5) and pole support (3.5)
Md <- cor(agg[, dim_names])
pairs_r <- Md[upper.tri(Md)]
pole_intra <- c(Md[systemic_pole, systemic_pole][upper.tri(diag(3))],
                Md[relational_pole, relational_pole][upper.tri(diag(3))])
pole_inter <- as.vector(Md[systemic_pole, relational_pole])

# 9d. Bootstrap CI for the HOI-on-HSI slope (paper 4.2), 10,000 country resamples
set.seed(20260525)
slopes <- replicate(10000, {
  i <- sample.int(nrow(m), replace = TRUE)
  coef(lm(HOI ~ HSI, data = m[i, ]))[2]
})
slope_ci <- quantile(slopes, c(0.025, 0.975))

# ------------------------------------------------------------------------------
# 10. Spatial / regional non-independence (paper Section 4.5)
#     Moran's I on the residuals of HOI ~ HSI with a macro-region contiguity
#     matrix (row-normalised), and region-clustered (CR1) standard errors.
#     Sample: countries with HSI + HOI matched to at least one benchmark, as in
#     the paper. Macro-regions from a fixed ISO-2 -> continent lookup (NIR, the
#     only non-ISO unit, assigned to Europe as in the pipeline).
# ------------------------------------------------------------------------------
REGION <- c(
  EG="Africa",ET="Africa",KE="Africa",LY="Africa",MA="Africa",NG="Africa",TN="Africa",ZW="Africa",
  AR="Americas",BO="Americas",BR="Americas",CA="Americas",CL="Americas",CO="Americas",EC="Americas",
  GT="Americas",MX="Americas",NI="Americas",PE="Americas",PR="Americas",US="Americas",UY="Americas",VE="Americas",
  AM="Asia",AZ="Asia",BD="Asia",CN="Asia",CY="Asia",GE="Asia",HK="Asia",ID="Asia",IN="Asia",IQ="Asia",
  IR="Asia",JO="Asia",JP="Asia",KG="Asia",KR="Asia",KZ="Asia",LB="Asia",MM="Asia",MN="Asia",MO="Asia",
  MV="Asia",MY="Asia",PH="Asia",PK="Asia",SG="Asia",TH="Asia",TJ="Asia",TR="Asia",TW="Asia",UZ="Asia",VN="Asia",
  AD="Europe",AL="Europe",AT="Europe",BA="Europe",BG="Europe",BY="Europe",CH="Europe",CZ="Europe",DE="Europe",
  DK="Europe",EE="Europe",ES="Europe",FI="Europe",FR="Europe",GB="Europe",GR="Europe",HR="Europe",HU="Europe",
  IS="Europe",IT="Europe",LT="Europe",LV="Europe",ME="Europe",MK="Europe",NIR="Europe",NL="Europe",NO="Europe",
  PL="Europe",PT="Europe",RO="Europe",RS="Europe",RU="Europe",SE="Europe",SI="Europe",SK="Europe",UA="Europe",
  AU="Oceania",NZ="Oceania"
)
moran_res <- NULL
if (have_benchmarks) {
  bench_iso3_all <- unique(unlist(lapply(names(benchmarks), function(bid) {
    b <- benchmarks[[bid]]
    read_benchmark(file.path(DATA_DIR, b$file), b$col)$iso3
  })))
  sp <- merge(agg[, c("cntry","iso3","HSI")], cl[, c("iso2","HOI")],
              by.x = "cntry", by.y = "iso2")
  sp <- sp[sp$iso3 %in% bench_iso3_all & is.finite(sp$HSI) & is.finite(sp$HOI), ]
  sp$region <- REGION[sp$cntry]
  sp <- sp[!is.na(sp$region), ]
  fit_sp <- lm(HOI ~ HSI, data = sp)
  beta_sp <- coef(fit_sp)["HSI"]
  se_ols  <- summary(fit_sp)$coefficients["HSI","Std. Error"]
  # region-clustered SE (CR1)
  X <- model.matrix(fit_sp); u <- residuals(fit_sp)
  XtXi <- solve(t(X) %*% X); G <- length(unique(sp$region)); n <- nrow(sp)
  meat <- matrix(0, ncol(X), ncol(X))
  for (g in unique(sp$region)) {
    sg <- t(X[sp$region == g, , drop = FALSE]) %*% u[sp$region == g]
    meat <- meat + sg %*% t(sg)
  }
  Vcl <- XtXi %*% meat %*% XtXi * (G/(G-1)) * ((n-1)/(n-ncol(X)))
  se_cl <- sqrt(diag(Vcl))["HSI"]
  # Moran's I with region-block contiguity weights (row-normalised)
  W <- outer(sp$region, sp$region, "=="); diag(W) <- FALSE
  W <- W / pmax(rowSums(W), 1)
  zres <- u - mean(u); S0 <- sum(W)
  I_obs <- (n/S0) * sum(zres * (W %*% zres)) / sum(zres^2)
  set.seed(20260525)
  I_perm <- replicate(999, { zp2 <- sample(zres)
    (n/S0) * sum(zp2 * (W %*% zp2)) / sum(zp2^2) })
  moran_res <- list(n = n, beta = beta_sp, t_ols = beta_sp/se_ols,
                    t_cl = beta_sp/se_cl, I = I_obs, EI = -1/(n-1),
                    p = (sum(I_perm >= I_obs) + 1) / 1000)
}

# ------------------------------------------------------------------------------
# 11. Discriminant validity against an EXTERNAL outcome (paper Section 4.6)
#     The WHR life ladder (independent source) replaces HOI as the dependent
#     variable, removing the shared-survey-source confound: incremental R2 of
#     the HSI over each benchmark, and the reverse contribution of each
#     benchmark over the HSI (complementarity reading).
# ------------------------------------------------------------------------------
PAPER_CANTRIL <- list(HDI = c(0.113, 0.380), LPI = c(0.042, 0.409),
                      TEHA_GAI = c(0.139, 0.160), WGI = c(0.048, 0.325),
                      GII = c(0.076, 0.228))
cantril_res <- list()
if (have_benchmarks && file.exists(file.path(DATA_DIR, "WHR_LifeLadder_2020.csv"))) {
  whr2 <- read_benchmark(file.path(DATA_DIR, "WHR_LifeLadder_2020.csv"), "life_ladder_2020")
  names(whr2)[2] <- "cantril"
  zs <- function(x) (x - mean(x)) / stats::sd(x)
  for (bid in names(benchmarks)) {
    b <- benchmarks[[bid]]
    bench <- read_benchmark(file.path(DATA_DIR, b$file), b$col)
    d <- merge(merge(agg_iso3, whr2, by = "iso3"), bench, by = "iso3")
    d <- d[stats::complete.cases(d[, c("HSI","cantril","score")]), ]
    d$Cz <- zs(d$cantril); d$Hz <- zs(d$HSI); d$Bz <- zs(d$score)
    dr2_fwd <- summary(lm(Cz ~ Bz + Hz, d))$r.squared - summary(lm(Cz ~ Bz, d))$r.squared
    dr2_rev <- summary(lm(Cz ~ Hz + Bz, d))$r.squared - summary(lm(Cz ~ Hz, d))$r.squared
    cantril_res[[bid]] <- c(n = nrow(d), fwd = dr2_fwd, rev = dr2_rev)
  }
}

# ------------------------------------------------------------------------------
# 12. Discriminant validity vs Welzel's emancipative values (paper Section 4.6)
#     Reads data_individual_welzel.csv (12 recoded items, emancipative direction,
#     [0,1]; VOICE already composed from the post-materialist priority items).
#     Only 2 of the 12 items (independence, obedience) also enter the HSI.
#     Decomposition: the Integral Score is the LEVEL component of the HSI; the
#     Compensation Delta (= sum of the penalty terms = Integral Score - raw HSI)
#     is its CONFIGURATIONAL component.
# ------------------------------------------------------------------------------
welzel_res <- NULL
PATH_WELZEL <- file.path(DATA_DIR, "data_individual_welzel.csv")
if (file.exists(PATH_WELZEL)) {
  wz <- read.csv(PATH_WELZEL, check.names = FALSE, fileEncoding = "UTF-8")
  sub_mean <- function(cols) rowMeans(wz[, cols, drop = FALSE], na.rm = TRUE)
  autonomy <- sub_mean(c("w_independence","w_imagination","w_obedience_rev"))
  equality <- sub_mean(c("w_jobs","w_pol_leaders","w_university"))
  choice   <- sub_mean(c("w_homosexuality","w_abortion","w_divorce"))
  voice    <- wz$w_voice
  wz$emanc <- rowMeans(cbind(autonomy, equality, choice, voice), na.rm = TRUE)
  wc <- aggregate(emanc ~ cntry_AN, data = wz, FUN = function(x) mean(x, na.rm = TRUE))
  dw <- merge(agg, wc, by.x = "cntry", by.y = "cntry_AN")
  dw$Delta <- dw$Integral_Score - dw$HSI_raw   # sum of penalty terms
  welzel_res <- list(
    n     = nrow(dw),
    r_p   = cor(dw$HSI, dw$emanc),
    r_s   = cor(dw$HSI, dw$emanc, method = "spearman"),
    r_lvl = cor(dw$emanc, dw$Integral_Score),
    r_cfg = cor(dw$emanc, dw$Delta)
  )
}

# ==============================================================================
# REPORT: reproduced value vs value reported in the paper
# ==============================================================================
chk <- function(label, got, paper) {
  cat(sprintf("  %-46s  computed = %-9s | paper = %s\n", label, got, paper))
}
cat("\n--------------------------------------------------------------\n")
cat(" REPRODUCED RESULTS  (computed vs reported in paper)\n")
cat("--------------------------------------------------------------\n")

top3 <- agg[order(-agg$HSI), ][1:3, ]
bot1 <- agg[order(agg$HSI), ][1, ]
chk("Top country HSI", sprintf("%s %.3f", top3$cntry[1], top3$HSI[1]), "SE 1.000")
chk("2nd / 3rd HSI",
    sprintf("%.3f/%.3f", top3$HSI[2], top3$HSI[3]), "0.946 / 0.943")
chk("Bottom country HSI", sprintf("%s %.3f", bot1$cntry, bot1$HSI), "EG 0.000")
it <- agg[agg$cntry == "IT", ]
chk("Italy HSI / rank", sprintf("%.3f / #%d", it$HSI, it$rank), "0.498 / #42")

ct <- table(factor(agg$class, levels = c("High","Good","Medium","Low")))
chk("Harmonic classes (H/G/M/L)",
    paste(ct, collapse = "/"), "12/29/42/9")
qt <- table(factor(agg$quadrant, levels = c("Q1","Q2","Q3","Q4")))
chk("Quadrants (Q1/Q2/Q3/Q4)",
    paste(qt, collapse = "/"), "22/24/29/17")

chk("HSI-HOI Pearson r", sprintf("%.3f", cor_hsi_hoi), "0.521")
chk("HSI-HOI R-squared",  sprintf("%.3f", r2_hsi), "0.271")
chk("Integral Score R-squared (S only)", sprintf("%.3f", r2_s), "0.188")

chk("Robustness: variance weights (rho)", sprintf("%.3f", rho_varw), ">= 0.99")
chk("Robustness: geometric mean (rho)",   sprintf("%.3f", rho_geom), ">= 0.99")
chk("Robustness: Monte Carlo median rho", sprintf("%.3f", rho_mc_med), "~0.997")
chk("  Monte Carlo 90% CI",
    sprintf("[%.3f,%.3f]", rho_mc_ci[1], rho_mc_ci[2]), "[0.992,0.999]")
chk("  High class stability (MC)", sprintf("%.1f%%", 100*stab12_mean), "~98.7%")

cat("\n  Discriminant validity - incremental R-squared of HSI over each benchmark:\n")
if (length(disc_results) > 0) {
  for (bid in names(disc_results)) {
    chk(sprintf("  dR2 over %s (n=%d)", bid, disc_results[[bid]]["n"]),
        sprintf("%.3f", disc_results[[bid]]["dR2"]),
        sprintf("%.3f", PAPER_DR2[bid]))
  }
} else {
  cat("  (benchmark files not found - section skipped)\n")
}
chk("HSI-WHR life ladder Pearson r",
    if (is.na(whr_r)) "NA" else sprintf("%.3f (n=%d)", whr_r, whr_n), "0.585")

cat("\n  Zero-penalty comparison (Section 4.5, Fig. 5):\n")
chk("  Spearman rho (HSI vs zero-penalty)", sprintf("%.3f", rho_zp), "0.870")
chk("  Countries changing class", sprintf("%d", n_change), "33")
chk("  Class direction (higher/lower under HSI)",
    sprintf("%d / %d", n_cls_up, n_cls_down), "28 / 5")
chk("  Rank among changers (improves/worsens)",
    sprintf("%d / %d", n_rk_impr, n_rk_wors), "25 / 8")
chk("  High members demoted without penalties",
    paste(high_demoted, collapse = ","), "CA,GB,NIR,US,UY (5)")

cat("\n  Penalty-weight robustness, individual perturbation (Section 4.5):\n")
chk("  min Spearman / min Kendall (20 specs)",
    sprintf("%.3f / %.3f", min(specs_rho_s), min(specs_rho_k)), ">=0.99 / >=0.94")
chk("  Top-10 max set difference (20 specs)",
    sprintf("%d", max(top10_same)), "0 (see note 3)")
chk("  High-class exchanges (25% | 50%)",
    sprintf("%d | %d", max(high_exch25), max(high_exch50)), "0 | <= 1")

cat("\n  Value-laden items and leave-one-item-out (Section 4.5):\n")
chk("  min rho (religion/obedience/joint removal)", sprintf("%.4f", rho_val_min), ">= 0.989")
chk("  LOIO min / median rho", sprintf("%.3f / %.3f", min(loo), median(loo)), ">= 0.989")
chk("  East Asia in Q4 (base / joint removal)",
    sprintf("%d/8 / %d/8", q4_asia(ref_full), q4_asia(r_no2)), "8/8 / core preserved (5/8)")

cat("\n  Dimensional structure (Sections 3.5, 4.5):\n")
chk("  Inter-dimensional r: range / median",
    sprintf("[%.2f, %.2f] / %.2f", min(pairs_r), max(pairs_r), median(pairs_r)),
    "[-0.29, 0.59] / 0.36")
chk("  Pairs |r|<0.30 / negative pairs",
    sprintf("%d/21 / %d", sum(abs(pairs_r) < 0.30), sum(pairs_r < 0)), "9/21 / 1")
chk("  Pole support: within / between",
    sprintf("%.2f / %.2f", mean(pole_intra), mean(pole_inter)), "0.46 / 0.15")
chk("HOI-on-HSI slope bootstrap 95% CI",
    sprintf("[%.3f, %.3f]", slope_ci[1], slope_ci[2]), "[0.090, 0.184]")

if (!is.null(moran_res)) {
  cat("\n  Spatial / regional non-independence (Section 4.5):\n")
  chk("  beta HSI / t OLS / t cluster-region",
      sprintf("%.3f / %.2f / %.2f", moran_res$beta, moran_res$t_ols, moran_res$t_cl),
      "0.141 / 5.85 / 7.08")
  chk("  Moran's I (E[I]) / permutation p",
      sprintf("%.3f (%.3f) / %.3f (n=%d)", moran_res$I, moran_res$EI, moran_res$p, moran_res$n),
      "0.078 (-0.011) / 0.023 (n=89)")
}

if (length(cantril_res) > 0) {
  cat("\n  Discriminant validity, EXTERNAL outcome = WHR life ladder (Section 4.6):\n")
  cat("  (forward: dR2 of HSI over benchmark; reverse: dR2 of benchmark over HSI)\n")
  for (bid in names(cantril_res)) {
    v <- cantril_res[[bid]]
    chk(sprintf("  %s fwd / rev (n=%d)", bid, v["n"]),
        sprintf("%.3f / %.3f", v["fwd"], v["rev"]),
        sprintf("%.3f / %.3f", PAPER_CANTRIL[[bid]][1], PAPER_CANTRIL[[bid]][2]))
  }
}

if (!is.null(welzel_res)) {
  cat("\n  Discriminant validity vs Welzel emancipative values (Section 4.6):\n")
  chk("  HSI ~ Welzel: Pearson / Spearman",
      sprintf("%.3f / %.3f (n=%d)", welzel_res$r_p, welzel_res$r_s, welzel_res$n),
      "0.654 / 0.530")
  chk("  Welzel ~ Integral Score (level)", sprintf("%.3f", welzel_res$r_lvl), "0.777")
  chk("  Welzel ~ Compensation Delta (config.)", sprintf("%.3f", welzel_res$r_cfg), "-0.434")
} else {
  cat("\n  (data_individual_welzel.csv not found - Welzel section skipped)\n")
}

cat("\n  NOTES ON MARGINAL DIFFERENCES\n")
cat("  1. Discriminant-validity, WHR, spatial and external-outcome values may\n")
cat("     differ marginally from the paper (typically by ~0.003-0.015): this\n")
cat("     reflects small differences in country coverage between this script -\n")
cat("     which uses a transparent standard ISO crosswalk - and the full\n")
cat("     pipeline, whose three-tier matching resolves a few additional\n")
cat("     countries. Region-clustered t statistics are the most sensitive to\n")
cat("     this (few clusters, G = 5). Substantive conclusions are unaffected.\n")
cat("  2. Bootstrap confidence intervals depend on the resampling protocol and\n")
cat("     seed; bounds may differ from the paper by ~0.01 without affecting\n")
cat("     the conclusion (interval clearly bounded away from zero).\n")
cat("  3. In the two most extreme penalty-weight specifications (x1.50) a\n")
cat("     single swap at the rank-10/11 boundary may appear in the top-10 set,\n")
cat("     reflecting third-decimal differences between this reduced rebuild\n")
cat("     and the certified pipeline values. The top-12 High class is\n")
cat("     preserved in all twenty specifications.\n")

cat("\n==============================================================\n")
cat(" Done.\n")
cat("==============================================================\n")
