# Robustness summary (Table 5) and few-cluster inference (Table 8) ----------

suppressPackageStartupMessages({library(dplyr); library(tidyr)})
df <- prepare_analysis_data()
focal <- institutional_terms

# Return the three focal estimates and country-clustered p-values in one row.
focal_row <- function(data, formula = baseline_formula, label, weights = NULL,
                      term_names = institutional_terms, reading = "") {
  data$.analysis_weight <- if (is.null(weights)) 1 else weights
  model <- lm(formula, data = data, weights = .analysis_weight)
  tab <- cluster_table(model, data$cn_id)
  get <- function(term, column) tab[match(term, tab$term), column]
  data.frame(check = label,
    focal_coefficient = get(term_names[1], "estimate"), focal_p_value = get(term_names[1], "p_value"),
    formal_coefficient = get(term_names[2], "estimate"), formal_p_value = get(term_names[2], "p_value"),
    interaction_coefficient = get(term_names[3], "estimate"), interaction_p_value = get(term_names[3], "p_value"),
    main_reading = reading)
}

rows <- list()
flat <- prepare_analysis_data(file.path(ROOT, "data", "analysis_ready_persistent_trust.rds"))
rows[[1]] <- focal_row(flat, label = "Flat-data OLS", reading = "Larger interaction when trust is treated as persistent")

# Residualized trust removes its linear projection on law and development.
df_resid <- df %>% mutate(trust_residual = resid(lm(trust_index ~ rule_of_law + gdp_pc + gdp_growth, data = .)),
                          residual_interaction = trust_residual * rule_of_law)
resid_formula <- update(baseline_formula, . ~ . - trust_index - trust_rule_of_law + trust_residual + residual_interaction)
rows[[2]] <- focal_row(df_resid, resid_formula, "Residualized trust", term_names = c("trust_residual", "rule_of_law", "residual_interaction"),
                       reading = "Larger interaction after removing linear overlap")

reduced <- c("se", "any_trim", "indstock", "jan", "value", "momentum", "price_group", "mkt_char", "gdp_pc")
reduced_formula <- as.formula(paste("size ~", paste(c(reduced, institutional_terms), collapse = " + ")))
reduced_df <- prepare_analysis_data(sample_controls = reduced)
rows[[3]] <- focal_row(reduced_df, reduced_formula, "Reduced-control OLS", reading = "Same signs using high-PIP controls only")

unwin <- df %>% mutate(size = size_unwinsorized, se = se_unwinsorized)
rows[[4]] <- focal_row(unwin, label = "No winsorization OLS", reading = "Pattern is not created by tail treatment")
win95 <- df %>% mutate(size = winsorize(size_unwinsorized, c(.05, .95)), se = winsorize(se_unwinsorized, c(.05, .95)))
rows[[5]] <- focal_row(win95, label = "5/95 winsorized OLS", reading = "Aggressive compression attenuates the result")
rows[[6]] <- focal_row(df, label = "Precision-weighted limiting diagnostic", weights = 1 / df$se^2,
                       reading = "Coefficients collapse under concentrated weights")

placebo_df <- df %>% drop_na(placebo_index) %>% mutate(placebo_interaction = placebo_index * rule_of_law)
placebo_formula <- update(baseline_formula, . ~ . - trust_index - trust_rule_of_law + placebo_index + placebo_interaction)
rows[[7]] <- focal_row(placebo_df, placebo_formula, "Family-trust placebo",
                       term_names = c("placebo_index", "rule_of_law", "placebo_interaction"),
                       reading = "Particularized trust does not reproduce generalized trust")

for (formal in c("regulatory_quality", "government_effectiveness", "control_of_corruption")) {
  alt <- df %>% drop_na(all_of(formal))
  int <- paste0("trust_", formal); alt[[int]] <- alt$trust_index * alt[[formal]]
  f <- as.formula(paste("size ~", paste(c(controls, "trust_index", formal, int), collapse = " + ")))
  label <- c(regulatory_quality = "Regulatory-quality substitute", government_effectiveness = "Government-effectiveness substitute",
             control_of_corruption = "Control-of-corruption substitute")[[formal]]
  rows[[length(rows) + 1]] <- focal_row(alt, f, label, term_names = c("trust_index", formal, int),
                                        reading = "Rule-of-law-specific pattern")
}

median_law <- median(df$rule_of_law)
split_model <- function(data, label) {
  data <- data %>% mutate(trust_centered = trust_index - mean(trust_index), low_law = as.numeric(rule_of_law <= median_law),
                          trust_low_law = trust_centered * low_law)
  f <- as.formula(paste("size ~", paste(c(controls, "trust_centered", "low_law", "trust_low_law"), collapse = " + ")))
  focal_row(data, f, label, term_names = c("trust_centered", "low_law", "trust_low_law"), reading = "Threshold sensitivity diagnostic")
}
rows[[length(rows) + 1]] <- split_model(df, "Low-rule-of-law median split, 1/99 winsorized")
rows[[length(rows) + 1]] <- split_model(win95, "Low-rule-of-law median split, 5/95 winsorized")
rows[[length(rows) + 1]] <- focal_row(filter(df, cn_id != "US"), label = "Non-U.S. OLS",
                                      reading = "Continuous interaction collapses")

# The manuscript's collapsed row contains two interaction estimates; retain
# both in explicit columns rather than placing two numbers in one CSV cell.
collapsed_country <- df %>% group_by(cn_id) %>% summarise(across(c(size, trust_index, rule_of_law), mean), .groups = "drop") %>%
  mutate(trust_rule_of_law = trust_index * rule_of_law)
collapsed_sc <- df %>% group_by(idstudy, cn_id) %>% summarise(across(c(size, trust_index, rule_of_law), mean), .groups = "drop") %>%
  mutate(trust_rule_of_law = trust_index * rule_of_law)
cc <- summary(lm(size ~ trust_index + rule_of_law + trust_rule_of_law, collapsed_country))$coef
cs <- summary(lm(size ~ trust_index + rule_of_law + trust_rule_of_law, collapsed_sc))$coef
collapsed <- data.frame(check = "Collapsed country and study-country evidence", focal_coefficient = NA, focal_p_value = NA,
  formal_coefficient = NA, formal_p_value = NA, interaction_coefficient = cc["trust_rule_of_law", "Estimate"],
  interaction_p_value = cc["trust_rule_of_law", "Pr(>|t|)"], main_reading = "No strong collapsed cross-country relationship",
  study_country_interaction = cs["trust_rule_of_law", "Estimate"],
  study_country_p_value = cs["trust_rule_of_law", "Pr(>|t|)"])
table5 <- bind_rows(c(rows, list(collapsed)))
write.csv(table5, file.path(TABLES, "Table_5.csv"), row.names = FALSE)

# Table 8: conventional CR1, CR2/Satterthwaite, and null-imposed wild cluster
# bootstrap p-values.  The explicit implementation below is retained because
# the RNG behavior of fwildclusterboot changed after the manuscript results
# were frozen.  It reproduces the null-imposed Rademacher procedure exactly.
model <- lm(baseline_formula, data = df)
cr1 <- cluster_table(model, df$cn_id)
cr2 <- clubSandwich::coef_test(model, vcov = "CR2", cluster = df$cn_id, test = "Satterthwaite")
cluster_vcov_matrix <- function(x, residuals, cluster, rank, xtx_inv = solve(crossprod(x))) {
  m <- length(unique(cluster)); n <- length(cluster)
  correction <- (m / (m - 1)) * ((n - 1) / (n - rank))
  scores <- rowsum(x * as.numeric(residuals), group = cluster, reorder = FALSE)
  out <- correction * xtx_inv %*% crossprod(scores) %*% xtx_inv
  dimnames(out) <- list(colnames(x), colnames(x)); out
}
wild_cluster_p <- function(term, reps = 4999L) {
  restricted <- lm(update(baseline_formula, paste(". ~ . -", term)), data = df)
  x <- model.matrix(model); xtx_inv <- solve(crossprod(x)); index <- match(term, colnames(x))
  observed_v <- cluster_vcov_matrix(x, resid(model), df$cn_id, model$rank, xtx_inv)
  observed_t <- coef(model)[[term]] / sqrt(observed_v[term, term])
  clusters <- sort(unique(df$cn_id)); cluster_index <- match(df$cn_id, clusters)
  boot_t <- numeric(reps)
  for (i in seq_len(reps)) {
    weights <- sample(c(-1, 1), length(clusters), replace = TRUE)
    fit <- lm.fit(x, fitted(restricted) + resid(restricted) * weights[cluster_index])
    v <- cluster_vcov_matrix(x, fit$residuals, df$cn_id, model$rank, xtx_inv)
    boot_t[[i]] <- fit$coefficients[[index]] / sqrt(v[term, term])
  }
  mean(abs(boot_t) >= abs(observed_t), na.rm = TRUE)
}
set.seed(20260607)
# The original diagnostics script evaluated the three unwinsorized terms
# immediately before the manuscript's 1/99 specification.  Advance the RNG by
# the identical Rademacher draws so the frozen bootstrap results are bitwise
# reproducible without needlessly refitting the unused preliminary models.
for (i in seq_len(3L * 4999L)) sample(c(-1, 1), length(unique(df$cn_id)), replace = TRUE)
wild_p <- vapply(focal, wild_cluster_p, numeric(1))
table8 <- data.frame(term = focal,
  cr1_p_value = cr1$p_value[match(focal, cr1$term)],
  cr2_satterthwaite_p_value = cr2$p_Satt[match(focal, rownames(cr2))],
  wild_bootstrap_p_value = wild_p)
write.csv(table8, file.path(TABLES, "Table_8.csv"), row.names = FALSE)
