# Baseline OLS, descriptive output, and Figures 1--3 -------------------------

suppressPackageStartupMessages({library(ggplot2); library(dplyr); library(tidyr)})
df <- prepare_analysis_data()
write.csv(validate_baseline_data(df), file.path(OUTPUTS, "data_validation.csv"), row.names = FALSE)

# Table 2 uses the transformed but unwinsorized model frame.
summary_source <- df %>% mutate(size = size_unwinsorized, se = se_unwinsorized)
summary_variables <- c(size = "Size slope", se = "Standard error",
                       trust_index = "Generalized trust index", rule_of_law = "Rule of law",
                       placebo_index = "Family-trust placebo index")
table2 <- bind_rows(lapply(names(summary_variables), function(v) {
  x <- summary_source[[v]]
  data.frame(variable = summary_variables[[v]], mean = mean(x, na.rm = TRUE),
             sd = sd(x, na.rm = TRUE), min = min(x, na.rm = TRUE),
             p25 = quantile(x, .25, na.rm = TRUE), median = median(x, na.rm = TRUE),
             p75 = quantile(x, .75, na.rm = TRUE), max = max(x, na.rm = TRUE))
}))
write.csv(table2, file.path(TABLES, "Table_2.csv"), row.names = FALSE)

baseline_model <- lm(baseline_formula, data = df)
study <- cluster_table(baseline_model, df$idstudy)
country <- cluster_table(baseline_model, df$cn_id)
focal <- institutional_terms
table3 <- merge(study[study$term %in% focal, ], country[country$term %in% focal, ],
                by = "term", suffixes = c("_study", "_country"))
table3 <- table3[match(focal, table3$term), ]
write.csv(table3, file.path(TABLES, "Table_3.csv"), row.names = FALSE)
saveRDS(list(model = baseline_model, study = study, country = country, data = df),
        file.path(OUTPUTS, "baseline_ols.rds"))

# Figure 1: distributions by country, excluding only countries with <5 rows.
plot_df <- df %>% add_count(cn_id, name = "country_n") %>% filter(country_n >= 5) %>%
  group_by(cn_id) %>% mutate(country_median = median(size)) %>% ungroup() %>%
  mutate(cn_id = reorder(cn_id, country_median))
size_limits <- range(df$size, na.rm = TRUE)
p1 <- ggplot(plot_df, aes(x = size, y = cn_id)) +
  geom_vline(xintercept = 0, linewidth = 0.35, color = "grey35") +
  geom_boxplot(outlier.shape = NA, width = 0.62, fill = "#E7EEF4",
               color = "#27323A", linewidth = 0.35) +
  coord_cartesian(xlim = size_limits) +
  labs(x = "Reported size slope", y = NULL) +
  theme_minimal(base_size = 10) +
  theme(panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(),
        axis.text.y = element_text(size = 8), plot.margin = margin(5.5, 8, 5.5, 5.5))

# Covariance-aware marginal effects.  Both clustering choices are retained in
# the CSV; the plotted bands use the two manuscript covariance matrices.
marginal_grid <- function(base, interaction, values, covariance, cluster_label) {
  b <- coef(baseline_model); z <- qnorm(.95)
  variance <- covariance[base, base] + values^2 * covariance[interaction, interaction] +
    2 * values * covariance[base, interaction]
  data.frame(value = values, effect = b[[base]] + values * b[[interaction]],
             low = b[[base]] + values * b[[interaction]] - z * sqrt(variance),
             high = b[[base]] + values * b[[interaction]] + z * sqrt(variance),
             cluster = cluster_label)
}
vc_study <- sandwich::vcovCL(baseline_model, cluster = df$idstudy, type = "HC1")
vc_country <- sandwich::vcovCL(baseline_model, cluster = df$cn_id, type = "HC1")
rol_values <- seq(min(df$rule_of_law), max(df$rule_of_law), length.out = 200)
trust_values <- seq(min(df$trust_index), max(df$trust_index), length.out = 200)
me_trust <- bind_rows(marginal_grid("trust_index", "trust_rule_of_law", rol_values, vc_study, "Study-clustered"),
                      marginal_grid("trust_index", "trust_rule_of_law", rol_values, vc_country, "Country-clustered"))
me_rol <- bind_rows(marginal_grid("rule_of_law", "trust_rule_of_law", trust_values, vc_study, "Study-clustered"),
                    marginal_grid("rule_of_law", "trust_rule_of_law", trust_values, vc_country, "Country-clustered"))
write.csv(me_trust, file.path(OUTPUTS, "Figure_2a_data.csv"), row.names = FALSE)
write.csv(me_rol, file.path(OUTPUTS, "Figure_2b_data.csv"), row.names = FALSE)
p2a <- ggplot(me_trust, aes(value, effect)) +
  geom_hline(yintercept = 0, color = "grey35", linewidth = 0.4, linetype = "dashed") +
  geom_ribbon(aes(ymin = low, ymax = high), fill = "#3B7EA1", alpha = 0.18) +
  geom_line(color = "#1B5E7A", linewidth = 0.9) +
  geom_point(data = data.frame(value = df$rule_of_law, y = 0), aes(value, y),
             inherit.aes = FALSE, color = "black", alpha = 0.35, size = 1) +
  facet_wrap(~cluster, ncol = 1) +
  labs(x = "Rule of law", y = "Marginal effect of generalized trust on the reported size slope") +
  theme_minimal(base_size = 12)
p2b <- ggplot(me_rol, aes(value, effect)) +
  geom_hline(yintercept = 0, color = "grey35", linewidth = 0.4, linetype = "dashed") +
  geom_ribbon(aes(ymin = low, ymax = high), fill = "#A65E2E", alpha = 0.18) +
  geom_line(color = "#7A4018", linewidth = 0.9) +
  geom_point(data = data.frame(value = df$trust_index, y = 0), aes(value, y),
             inherit.aes = FALSE, color = "black", alpha = 0.35, size = 1) +
  facet_wrap(~cluster, ncol = 1) +
  labs(x = "Generalized trust", y = "Marginal effect of rule of law on the reported size slope") +
  theme_minimal(base_size = 12)

# Figure 3 isolates the contribution of the three institutional terms.
b <- coef(baseline_model)
surface <- expand.grid(rule_of_law = seq(min(df$rule_of_law), max(df$rule_of_law), length.out = 200),
                       trust_index = seq(min(df$trust_index), max(df$trust_index), length.out = 200)) %>%
  mutate(partial_effect = b[["trust_index"]] * trust_index + b[["rule_of_law"]] * rule_of_law +
           b[["trust_rule_of_law"]] * trust_index * rule_of_law)
write.csv(surface, file.path(OUTPUTS, "Figure_3_data.csv"), row.names = FALSE)
p3_limit <- max(abs(surface$partial_effect))
p3 <- ggplot(surface, aes(rule_of_law, trust_index)) +
  geom_tile(aes(fill = partial_effect)) +
  geom_point(data = df, aes(rule_of_law, trust_index), inherit.aes = FALSE,
             color = "black", size = 1.2, alpha = 0.55) +
  scale_fill_gradient2(low = "#A53A2A", mid = "#F5F1E8", high = "#2B6F8A",
                       midpoint = 0, limits = c(-p3_limit, p3_limit), name = "Partial effect") +
  labs(x = "Rule of law", y = "Generalized trust") + theme_minimal(base_size = 12)

save_plot <- function(plot, stem, width, height) {
  ggsave(file.path(FIG_PNG, paste0(stem, ".png")), plot, width = width, height = height, dpi = 300)
  ggsave(file.path(FIG_PDF, paste0(stem, ".pdf")), plot, width = width, height = height,
         device = grDevices::pdf, version = "1.5")
}
save_plot(p1, "Figure_1", 6.8, 5.6); save_plot(p2a, "Figure_2a", 4.0, 5.2)
save_plot(p2b, "Figure_2b", 4.0, 5.2); save_plot(p3, "Figure_3", 8.0, 7.0)

# Country-level, not observation-weighted, quantiles used in the manuscript.
country_institutions <- df %>% group_by(cn_id) %>%
  summarise(trust_index = mean(trust_index), rule_of_law = mean(rule_of_law), .groups = "drop")
table7 <- bind_rows(lapply(c("trust_index", "rule_of_law"), function(v) data.frame(
  variable = v, p25 = quantile(country_institutions[[v]], .25),
  median = median(country_institutions[[v]]), p75 = quantile(country_institutions[[v]], .75))))
write.csv(table7, file.path(TABLES, "Table_7.csv"), row.names = FALSE)

# Illustrative economic-effects calculation from the Discussion -------------
# This is a point-estimate scaling exercise in reported-slope space, not a
# structural counterfactual. Institutional quantiles are calculated with one
# observation per country so the heavily represented United States does not
# determine the comparison. The slope percentile calculations then return to
# the complete estimate-level sample, exactly as described in the manuscript.
trust_p25 <- table7$p25[table7$variable == "trust_index"]
trust_p75 <- table7$p75[table7$variable == "trust_index"]
rol_p25 <- table7$p25[table7$variable == "rule_of_law"]
rol_p75 <- table7$p75[table7$variable == "rule_of_law"]
trust_change <- trust_p75 - trust_p25

trust_marginal_effect <- function(rule_of_law_value) {
  coef(baseline_model)[["trust_index"]] +
    coef(baseline_model)[["trust_rule_of_law"]] * rule_of_law_value
}
low_rol_marginal_effect <- trust_marginal_effect(rol_p25)
high_rol_marginal_effect <- trust_marginal_effect(rol_p75)
low_rol_fitted_change <- trust_change * low_rol_marginal_effect
high_rol_fitted_change <- trust_change * high_rol_marginal_effect

slope_p25 <- as.numeric(quantile(df$size_unwinsorized, 0.25))
low_rol_implied_slope <- slope_p25 + low_rol_fitted_change
low_rol_implied_percentile <- mean(df$size_unwinsorized <= low_rol_implied_slope)
between_count <- sum(df$size_unwinsorized >= slope_p25 &
                       df$size_unwinsorized <= low_rol_implied_slope)
between_share <- between_count / nrow(df)
country_weighted_mean_slope <- df %>%
  group_by(cn_id) %>%
  summarise(country_mean_slope = mean(size_unwinsorized), .groups = "drop") %>%
  summarise(value = mean(country_mean_slope)) %>%
  pull(value)

economic_effects <- data.frame(
  quantity = c(
    "trust_country_p25", "trust_country_p75", "trust_p75_minus_p25",
    "rule_of_law_country_p25", "rule_of_law_country_p75",
    "baseline_trust_coefficient", "baseline_interaction_coefficient",
    "trust_marginal_effect_at_rule_of_law_p25", "fitted_slope_change_at_rule_of_law_p25",
    "unwinsorized_size_slope_p25", "implied_slope_from_p25_at_rule_of_law_p25",
    "implied_slope_percentile", "estimates_between_observed_p25_and_implied_slope",
    "share_between_observed_p25_and_implied_slope", "country_weighted_mean_unwinsorized_slope",
    "trust_marginal_effect_at_rule_of_law_p75", "fitted_slope_change_at_rule_of_law_p75"
  ),
  value = c(
    trust_p25, trust_p75, trust_change, rol_p25, rol_p75,
    coef(baseline_model)[["trust_index"]], coef(baseline_model)[["trust_rule_of_law"]],
    low_rol_marginal_effect, low_rol_fitted_change, slope_p25, low_rol_implied_slope,
    low_rol_implied_percentile, between_count, between_share, country_weighted_mean_slope,
    high_rol_marginal_effect, high_rol_fitted_change
  ),
  unit = c(
    rep("index units", 5), rep("reported-slope units per index unit", 3),
    "reported-slope units", "reported-slope units", "reported-slope units",
    "share of estimates", "number of estimates", "share of estimates",
    "reported-slope units", "reported-slope units per index unit", "reported-slope units"
  ),
  manuscript_reading = c(
    "Country-level trust 25th percentile", "Country-level trust 75th percentile",
    "Trust increase of about 1.45 index units", "Lower-rule-of-law environment",
    "Higher-rule-of-law environment", "Baseline OLS trust coefficient",
    "Baseline OLS trust-by-rule-of-law coefficient", "Trust slope at lower rule of law",
    "Reported size slope rises by about 0.18", "Starting reported-slope percentile value",
    "Implied reported size slope of about 0.08", "Implied slope is around the 93rd percentile",
    "Number of estimates between the observed and implied slopes",
    "Roughly two thirds of the estimation sample", "Country-weighted baseline mean slope",
    "Trust slope at higher rule of law", "Reported size slope changes by about -0.04"
  )
)
write.csv(economic_effects, file.path(OUTPUTS, "economic_effects_calculation.csv"), row.names = FALSE)
