#========================================================================= # Supplementary Code: Translated CES and Epistemic Deduction # for Non-Compensatory Composite Indicators (Online Resource 1) #========================================================================= # Article Title: Normative Translation and Pre-Screening to Overcome Domain Collapse... # Journal Name: Quality & Quantity # R Version: >= 4.1.0 (Utilizes native pipe |> and shorthand lambda \()) # Dependencies: 'TruncatedNormal', 'Matrix', 'matrixStats', 'withr', 'sensitivity' #========================================================================= # --- 0. Dependency Verification & Constants (Fail-Fast Best Practice) --- EPSILON <- 1e-6 PD_TOLERANCE <- 1e-8 VAR_FLOOR <- 1e-12 required_pkgs <- c("TruncatedNormal", "Matrix", "matrixStats", "withr", "sensitivity") pkg_status <- required_pkgs |> vapply(requireNamespace, logical(1), quietly = TRUE) missing_pkgs <- required_pkgs[!pkg_status] if (length(missing_pkgs) > 0) { stop(sprintf("Missing required packages: %s. Please install them before running.", paste(missing_pkgs, collapse = ", ")), call. = FALSE) } #------------------------------------------------------------------- # 1. Core function: Estimate baseline expectation and SD of lower-tail aggregation #------------------------------------------------------------------- #' @title Estimate Baseline Expectation and SD of Lower-Tail Aggregation #' @description Core function to estimate baseline expectation and SD of lower-tail aggregation. null_expectation_sd <- function(sigma = NULL, k, n_sim = 1e4, censor_low = NULL, censor_high = NULL, Sigma = NULL, boundary_type = "truncate", force_pd = FALSE, seed = NULL) { run_sim <- function() { if (!is.null(Sigma) && !is.null(sigma)) { warning("Both sigma and Sigma provided. Prioritizing Sigma matrix.", call. = FALSE) } if (is.null(Sigma)) { stopifnot(!is.null(sigma)) K <- length(sigma) stopifnot(k >= 1, k <= K) sim_vals <- matrix(rnorm(n_sim * K), nrow = n_sim, ncol = K) |> sweep(2, sigma, "*") } else { stopifnot(is.matrix(Sigma)) K <- nrow(Sigma) stopifnot(ncol(Sigma) == K, k >= 1, k <= K) is_pd <- all(eigen(Sigma, symmetric = TRUE, only.values = TRUE)$values > PD_TOLERANCE) if (!is_pd) { if (force_pd) { orig_vars <- diag(Sigma) Sigma_pd <- as.matrix(Matrix::nearPD(Sigma, corr = FALSE)$mat) pd_vars <- diag(Sigma_pd) D_mat <- diag(sqrt(pmax(orig_vars, VAR_FLOOR) / pmax(pd_vars, VAR_FLOOR))) Sigma <- D_mat %*% Sigma_pd %*% D_mat if (any(eigen(Sigma, symmetric = TRUE, only.values = TRUE)$values < PD_TOLERANCE)) { Sigma <- as.matrix(Matrix::nearPD(Sigma, corr = FALSE)$mat) } } else { stop("Sigma not positive-definite. Set force_pd=TRUE to auto-correct.", call. = FALSE) } } U <- chol(Sigma) Z <- matrix(rnorm(n_sim * K), nrow = n_sim, ncol = K) sim_vals <- Z %*% U } boundary_type <- match.arg(boundary_type, choices = c("censor", "truncate")) if (!is.null(censor_low) || !is.null(censor_high)) { if (boundary_type == "censor") { if (!is.null(censor_low)) sim_vals[sim_vals < censor_low] <- censor_low if (!is.null(censor_high)) sim_vals[sim_vals > censor_high] <- censor_high } else if (boundary_type == "truncate") { if (!is.null(Sigma)) { lower_bounds <- rep(ifelse(is.null(censor_low), -Inf, censor_low), K) upper_bounds <- rep(ifelse(is.null(censor_high), Inf, censor_high), K) sim_vals <- tryCatch({ TruncatedNormal::mvrandn(l = lower_bounds, u = upper_bounds, Sig = Sigma, n = n_sim, mu = rep(0, K)) }, error = function(e) { warning("mvrandn failed. Falling back to independent marginal truncation.", call. = FALSE) NULL }) if (is.null(sim_vals)) { sds <- sqrt(diag(Sigma)) low_p <- ifelse(is.null(censor_low), 0, pnorm(censor_low, sd = sds)) high_p <- ifelse(is.null(censor_high), 1, pnorm(censor_high, sd = sds)) sim_vals <- matrix(runif(n_sim * K), nrow = n_sim, ncol = K) |> sweep(2, high_p - low_p, "*") |> sweep(2, low_p, "+") |> as.vector() |> qnorm(sd = rep(sds, each = n_sim)) |> matrix(nrow = n_sim, ncol = K) warning("Silent Correlation Destruction: Fell back to independent marginal truncation.", call. = FALSE) } else { if (nrow(sim_vals) != n_sim) sim_vals <- t(sim_vals) } } else { low_p <- ifelse(is.null(censor_low), 0, pnorm(censor_low, sd = sigma)) high_p <- ifelse(is.null(censor_high), 1, pnorm(censor_high, sd = sigma)) sim_vals <- matrix(runif(n_sim * K), nrow = n_sim, ncol = K) |> sweep(2, high_p - low_p, "*") |> sweep(2, low_p, "+") |> as.vector() |> qnorm(sd = rep(sigma, each = n_sim)) |> matrix(nrow = n_sim, ncol = K) } } } calc_trimmed_tail <- function(row_vals) { idx <- order(row_vals)[seq_len(k)] tail_vals <- row_vals[idx] if (k >= 3) { trim_count <- max(1, floor(k * 0.20)) trimmed_vals <- tail_vals[seq_len(k - trim_count)] mean(trimmed_vals) } else { mean(tail_vals) } } raw_scores <- apply(sim_vals, 1, calc_trimmed_tail) list(null_exp = mean(raw_scores), null_sd = sd(raw_scores), sim_raw_scores = raw_scores) } if (!is.null(seed)) { withr::with_seed(seed, run_sim()) } else { run_sim() } } #------------------------------------------------------------------- # Helper function: Calculate Normative Translation Constants (C_global) #------------------------------------------------------------------- calculate_normative_C_global <- function(normative_floors_vector, epsilon = EPSILON) { max(abs(normative_floors_vector)) + epsilon } #------------------------------------------------------------------- # Helper function: Pessimistic Single Imputation (Max-Min MEU compliant) #------------------------------------------------------------------- pessimistic_impute <- function(mu_pred, sigma_pred, percentile = 0.10) { z_score <- qnorm(percentile) mu_pred + (z_score * sigma_pred) } #------------------------------------------------------------------- # Helper function: Multivariate Pessimistic Imputation (Gaussian Copula) #------------------------------------------------------------------- #' @title Multivariate Pessimistic Imputation via Gaussian Copula #' @description Approximates the deterministic conditional expectation of a #' multivariate truncated normal distribution to handle systemic MNAR missingness. pessimistic_impute_multivariate <- function(mu_pred_vec, sigma_pred_vec, R_miss, percentile = 0.10, n_draws = 1000) { if (!requireNamespace("TruncatedNormal", quietly = TRUE)) { stop("Requires 'TruncatedNormal' package for multivariate imputation.", call. = FALSE) } K <- length(mu_pred_vec) if (length(sigma_pred_vec) != K || nrow(R_miss) != K || ncol(R_miss) != K) { stop("Dimension mismatch between mu, sigma, and correlation matrix R_miss.", call. = FALSE) } Sigma <- diag(sigma_pred_vec) %*% R_miss %*% diag(sigma_pred_vec) is_pd <- all(eigen(Sigma, symmetric = TRUE, only.values = TRUE)$values > PD_TOLERANCE) if (!is_pd) { Sigma <- as.matrix(Matrix::nearPD(Sigma, corr = FALSE)$mat) } z_bound <- qnorm(percentile) upper_bounds <- mu_pred_vec + z_bound * sigma_pred_vec draws <- tryCatch({ TruncatedNormal::mvrandn(l = rep(-Inf, K), u = upper_bounds, Sig = Sigma, n = n_draws, mu = mu_pred_vec) }, error = function(e) { stop("mvrandn failed in multivariate imputation. The MEU joint worst-case assumption cannot be satisfied with independent fallbacks. Check your correlation matrix R_miss and bounds.", call. = FALSE) }) return(rowMeans(draws)) } #------------------------------------------------------------------- # 2. Lower-Tail Trimmed Mean (Null-centering bypassed to preserve absolute score) #------------------------------------------------------------------- lower_tail_owa <- function(observed, sigma = NULL, k, n_sim = 1e4, censor_low = NULL, censor_high = NULL, Sigma = NULL, boundary_type = "truncate", force_pd = FALSE, seed = NULL, calculate_shrinkage_diagnostic = FALSE, external_lambda = NULL, precomputed_null_exp = NULL) { if (anyNA(observed)) { stop("Domain Error: NA detected in indicators. Rigorous imputation required.", call. = FALSE) } K_valid <- length(observed) if (K_valid < k) stop("Insufficient indicators to evaluate lower-tail window k.", call. = FALSE) sorted_indices <- order(observed) tail_indices <- sorted_indices[seq_len(k)] tail_obs <- observed[tail_indices] if (k >= 3) { trim_count <- max(1, floor(k * 0.20)) sorted_tail <- sort(tail_obs) trimmed_tail <- sorted_tail[seq_len(k - trim_count)] raw_score <- mean(trimmed_tail) } else { raw_score <- mean(tail_obs) } if (is.null(precomputed_null_exp)) { warning("Computing null expectation on the fly. Pre-compute outside region loops.", call. = FALSE) ne <- null_expectation_sd(sigma = sigma, k = k, n_sim = n_sim, censor_low = censor_low, censor_high = censor_high, Sigma = Sigma, boundary_type = boundary_type, force_pd = force_pd, seed = seed) null_exp <- ne$null_exp null_sd <- ne$null_sd percentile_val <- mean(ne$sim_raw_scores <= raw_score) } else { null_exp <- precomputed_null_exp if (calculate_shrinkage_diagnostic) { ne <- null_expectation_sd(sigma = sigma, k = k, n_sim = n_sim, censor_low = censor_low, censor_high = censor_high, Sigma = Sigma, boundary_type = boundary_type, force_pd = force_pd, seed = seed) null_sd <- ne$null_sd } else { null_sd <- NA_real_ } percentile_val <- NA_real_ } lambda <- 0 corrected_score <- raw_score if (calculate_shrinkage_diagnostic) { if (k < 3) { warning("k < 3: Empirical null variance is statistically unstable. Setting lambda = 0.", call. = FALSE) var_tail <- NA } else { var_tail <- null_sd^2 sigma_tail_sq <- mean(sigma[tail_indices]^2) lambda <- max(0, 1 - (sigma_tail_sq / var_tail)) } } list(raw = raw_score, null_exp = null_exp, lambda = lambda, null_sd = null_sd, corrected = corrected_score, percentile = percentile_val) } #------------------------------------------------------------------- # 3. Min-Ratio Stabilized Global Normative-Anchor Translated CES #------------------------------------------------------------------- ces_translated_stable <- function(scores, weights, rho, C_vec, sigma_eps_vec = NULL, Z_score = 0, L_min = NULL) { if (anyNA(scores)) stop("Domain Error: NA detected in scores.", call. = FALSE) stopifnot(rho < 0, length(scores) == length(weights), length(C_vec) == length(scores)) if (!is.null(L_min)) { scores <- pmax(scores, L_min) } if (abs(sum(weights) - 1) > 1e-6) { stop("Domain Error: Weights do not sum to 1. Explicitly normalize weights prior to aggregation.", call. = FALSE) } C_global <- max(C_vec) C_vec <- rep(C_global, length(C_vec)) X <- scores + C_vec if (any(X <= 0)) stop("Domain Collapse: X_d <= 0 detected. The Normative C_global buffer is insufficient.", call. = FALSE) X_min <- min(X) inner_sum_X <- sum(weights * ((X / X_min) ^ rho)) Y_inner_unadjusted <- X_min * (inner_sum_X ^ (1 / rho)) Y_baseline_unadjusted <- Y_inner_unadjusted - C_global if (!is.null(sigma_eps_vec) && Z_score > 0) { floor_bound <- if (is.null(L_min)) -C_global + EPSILON else L_min + EPSILON S_star_unbounded <- scores - Z_score * sigma_eps_vec C_unbounded <- C_global + Z_score * max(sigma_eps_vec) + EPSILON X_star_unbounded <- S_star_unbounded + C_unbounded X_star_unbounded_min <- min(X_star_unbounded) inner_sum_X_star_unbounded <- sum(weights * ((X_star_unbounded / X_star_unbounded_min) ^ rho)) Y_ces_unbounded <- X_star_unbounded_min * (inner_sum_X_star_unbounded ^ (1 / rho)) - C_unbounded S_star <- pmax(S_star_unbounded, floor_bound) X_star <- S_star + C_vec X_star_min <- min(X_star) inner_sum_X_star <- sum(weights * ((X_star / X_star_min) ^ rho)) Y_ces_final <- X_star_min * (inner_sum_X_star ^ (1 / rho)) - C_global endogenous_weights <- (weights * ((X_star / X_star_min) ^ rho)) / inner_sum_X_star bottleneck_unadj <- which.min(X) bottleneck_star <- which.min(X_star) epistemic_shift_flag <- FALSE # Calculate the unbounded pre-screened score for the new bottleneck S_star_unbounded_bottleneck <- scores[bottleneck_star] - Z_score * sigma_eps_vec[bottleneck_star] # Strictly enforce the manuscript's definition: must cross the zero-point into deprivation if (bottleneck_unadj != bottleneck_star && scores[bottleneck_star] >= 0 && S_star_unbounded_bottleneck < 0) { epistemic_shift_flag <- TRUE } list(Y_ces = Y_ces_final, endogenous_weights = endogenous_weights, Y_unadjusted = Y_baseline_unadjusted, Y_unbounded = Y_ces_unbounded, epistemic_bottleneck_flag = epistemic_shift_flag) } else { endogenous_weights <- (weights * ((X / X_min) ^ rho)) / inner_sum_X list(Y_ces = Y_baseline_unadjusted, endogenous_weights = endogenous_weights, Y_unadjusted = Y_baseline_unadjusted, Y_unbounded = Y_baseline_unadjusted, epistemic_bottleneck_flag = FALSE) } } #------------------------------------------------------------------- # 4. Certainty-Equivalent Epistemic Deduction #------------------------------------------------------------------- apply_epistemic_penalty <- function(Y_ces_adjusted, Y_ces_unadjusted, Y_ces_unbounded, L_min = NULL) { if (is.null(L_min)) stop("Normative Floor (L_min) must be explicitly defined.", call. = FALSE) penalty <- Y_ces_unadjusted - Y_ces_adjusted Y_final <- max(L_min, Y_ces_adjusted) list(Y_final = Y_final, Y_unbounded = Y_ces_unbounded, penalty_component = penalty) } #------------------------------------------------------------------- # 5. Monte Carlo Rank-Stability Validation (Fully Vectorized R 4.x) #------------------------------------------------------------------- monte_carlo_rank_stability <- function(n_regions = 10000, D = 3, rho = -50, seed = 42) { run_mc <- function() { raw_gamma <- rgamma(n_regions * D, shape = 2, scale = 20) true_scores <- matrix(-pmin(raw_gamma, 100), nrow = n_regions, ncol = D) noise_sd <- matrix(runif(n_regions * D, 5, 25), nrow = n_regions, ncol = D) noise <- matrix(rnorm(n_regions * D, 0, noise_sd), nrow = n_regions, ncol = D) observed_scores <- pmax(true_scores + noise, -100) C_global <- 100 + EPSILON L_min <- -100 weights <- rep(1/D, D) Z_score <- 1.96 latent_capacity <- runif(n_regions, 10, 30) metadata_proxy <- latent_capacity + rnorm(n_regions * D, 0, 5) sigma_eps_mat <- pmax(metadata_proxy, 1) metadata_proxy_adversarial <- 0.3 * latent_capacity + 0.7 * runif(n_regions * D, 5, 25) sigma_eps_mat_adv <- pmax(metadata_proxy_adversarial, 1) vec_ces <- function(scores_mat, sigma_eps_mat) { S_star <- pmax(scores_mat - Z_score * sigma_eps_mat, L_min + EPSILON) X_star <- S_star + C_global X_star_min <- matrixStats::rowMins(X_star) ratio_X <- X_star / X_star_min inner_sum_star <- rowSums(sweep(ratio_X^rho, 2, weights, "*")) Y_ces <- X_star_min * (inner_sum_star^(1/rho)) - C_global pmax(Y_ces, L_min) } Y_norm_epistemic <- vec_ces(observed_scores, sigma_eps_mat) Y_norm_epistemic_adv <- vec_ces(observed_scores, sigma_eps_mat_adv) C_noise_mat <- 100 + (noise_sd * 2) X_noise <- observed_scores + C_noise_mat X_noise_min <- matrixStats::rowMins(X_noise) ratio_X_noise <- X_noise / X_noise_min inner_noise <- rowSums(sweep(ratio_X_noise^rho, 2, weights, "*")) min_C_noise <- matrixStats::rowMins(C_noise_mat) Y_noise_buffered <- X_noise_min * (inner_noise^(1/rho)) - min_C_noise rank_ontological <- rank(matrixStats::rowMins(true_scores)) cor_noise <- cor(rank_ontological, rank(Y_noise_buffered), method = "spearman") cor_norm <- cor(rank_ontological, rank(Y_norm_epistemic), method = "spearman") cor_norm_adv <- cor(rank_ontological, rank(Y_norm_epistemic_adv), method = "spearman") list( spearman_noise_buffered = cor_noise, spearman_norm_epistemic = cor_norm, spearman_norm_epistemic_adversarial = cor_norm_adv, message = sprintf("Noise-Buffered (r_s vs True): %.3f | Norm+Epistemic (Correlated r_s vs True): %.3f | Norm+Epistemic (Adversarial r_s vs True): %.3f", cor_noise, cor_norm, cor_norm_adv) ) } if (!is.null(seed)) { withr::with_seed(seed, run_mc()) } else { run_mc() } } #------------------------------------------------------------------- # 6. Morris Sensitivity Analysis Example #------------------------------------------------------------------- morris_sensitivity_example <- function(r = 100, base_C_global = 100, base_rho = -50, base_Z_score = 1.96, L_min = -100, sigma_eps_vec = c(15, 25, 10), weights = c(1/3, 1/3, 1/3), mock_scores = c(-20, -40, -10)) { if (!requireNamespace("sensitivity", quietly = TRUE)) { stop("Requires 'sensitivity' package.", call. = FALSE) } ces_model <- function(X) { vapply(seq_len(nrow(X)), function(i) { C_glob <- X[i, 1] rho_val <- X[i, 2] Z_val <- X[i, 3] C_vec <- rep(C_glob, length(mock_scores)) ces_out <- ces_translated_stable(scores = mock_scores, weights = weights, rho = rho_val, C_vec = C_vec, sigma_eps_vec = sigma_eps_vec, Z_score = Z_val, L_min = L_min) penalty_res <- apply_epistemic_penalty(ces_out$Y_ces, ces_out$Y_unadjusted, ces_out$Y_unbounded, L_min = L_min) penalty_res$Y_final }, numeric(1)) } set.seed(42) x <- sensitivity::morris(model = NULL, factors = 3, r = r, design = list(type = "oat", levels = 10, grid.jump = 5)) X_scaled <- x$X X_scaled[, 1] <- X_scaled[, 1] * 45 + 105 X_scaled[, 2] <- X_scaled[, 2] * 50 - 100 X_scaled[, 3] <- X_scaled[, 3] * (2.58 - 1.64) + 1.64 y <- ces_model(X_scaled) x <- sensitivity::tell(x, y) mu_star_scaled <- apply(x$ee, 2, function(col) mean(abs(col), na.rm = TRUE)) sigma_scaled <- apply(x$ee, 2, sd, na.rm = TRUE) range_C <- 150 - 105 range_rho <- -50 - (-100) range_Z <- 2.58 - 1.64 mu_star_orig <- mu_star_scaled / c(range_C, range_rho, range_Z) sigma_orig <- sigma_scaled / c(range_C, range_rho, range_Z) cat("\n--- Original-Scale Morris Mean Absolute Effects (Copy these to Table 2) ---\n") print(data.frame( Parameter = c("C_global", "rho", "Z"), mu_star_orig = round(mu_star_orig, 3), sigma_orig = round(sigma_orig, 3) )) return(x) } #------------------------------------------------------------------- # 7. Reproduce Stylized Numerical Trace #------------------------------------------------------------------- reproduce_stylized_trace <- function() { rho <- -50 C_global <- 100 + EPSILON L_min <- -100 Z_score <- 1.96 weights <- c(1/3, 1/3, 1/3) scores_X <- c(-40, -50, -45); sigma_eps_X <- c(5, 5, 5) scores_Y <- c(-30, -60, -20); sigma_eps_Y <- c(5, 25, 5) scores_Z <- c(60, 70, 65); sigma_eps_Z <- c(5, 5, 5) calc_noise <- function(scores, sigma_eps) { C_d <- 100 + (sigma_eps * 2) X_noise <- scores + C_d X_min <- min(X_noise) inner <- sum(weights * ((X_noise / X_min) ^ rho)) Y_ces <- X_min * (inner ^ (1 / rho)) - min(C_d) return(Y_ces) } Y_noise_X <- calc_noise(scores_X, sigma_eps_X) Y_noise_Y <- calc_noise(scores_Y, sigma_eps_Y) Y_noise_Z <- calc_noise(scores_Z, sigma_eps_Z) res_X <- ces_translated_stable(scores_X, weights, rho, rep(C_global, 3), sigma_eps_X, Z_score, L_min) res_Y <- ces_translated_stable(scores_Y, weights, rho, rep(C_global, 3), sigma_eps_Y, Z_score, L_min) res_Z <- ces_translated_stable(scores_Z, weights, rho, rep(C_global, 3), sigma_eps_Z, Z_score, L_min) message("--- Stylized Numerical Trace Reproduction ---") message(sprintf("Region X (Noise-Buffered): %.1f | Normative+Epistemic: %.1f", Y_noise_X, res_X$Y_ces)) message(sprintf("Region Y (Noise-Buffered): %.1f | Normative+Epistemic: %.1f", Y_noise_Y, res_Y$Y_ces)) message(sprintf("Region Z (Noise-Buffered): %.1f | Normative+Epistemic: %.1f", Y_noise_Z, res_Z$Y_ces)) message("---------------------------------------------") list( Region_X = list(Noise_Buffered = Y_noise_X, Normative_Epistemic = res_X$Y_ces), Region_Y = list(Noise_Buffered = Y_noise_Y, Normative_Epistemic = res_Y$Y_ces), Region_Z = list(Noise_Buffered = Y_noise_Z, Normative_Epistemic = res_Z$Y_ces) ) } #------------------------------------------------------------------- # 8. Derivation of Equation 2 Multiplier (Simulation Sweep) #------------------------------------------------------------------- derive_rho_bound_multiplier <- function(n_regions = 10000, D = 3, seed = 42) { if (!requireNamespace("matrixStats", quietly = TRUE)) stop("Requires 'matrixStats'.") set.seed(seed) message("Deriving Equation 2 Safe Harbor Multiplier...") raw <- rgamma(n_regions * D, shape = 2, scale = 20) scores <- matrix(-pmin(raw, 100), nrow = n_regions, ncol = D) C_global <- 100 + EPSILON sigma_S <- 200 / sqrt(12) sigma_S_sweep <- c(sigma_S * 0.8, sigma_S, sigma_S * 1.2) true_min <- matrixStats::rowMins(scores) weights <- rep(1/D, D) mults <- c(5, 8, 10, 12, 15) for (s in sigma_S_sweep) { message(sprintf("\n--- Sensitivity Sweep: sigma_S = %.2f ---", s)) for (m in mults) { rho_test <- -1 * m * (C_global / s) X <- scores + C_global X_min <- matrixStats::rowMins(X) ratio_X <- X / X_min inner <- rowSums(sweep(ratio_X^rho_test, 2, weights, "*")) Y_ces <- X_min * (inner^(1/rho_test)) - C_global drift_abs <- mean(abs(Y_ces - true_min)) drift_rel <- mean(abs(Y_ces - true_min) / C_global) message(sprintf("Multiplier: %2d | |rho|: %7.2f | Abs Drift: %.2f pts | Rel Drift (vs C_global): %.4f%%", m, abs(rho_test[1]), drift_abs, drift_rel * 100)) } } message("\nDerivation sweep complete. Verify that Rel Drift is < 3% at your chosen multiplier.") }