# Evangélicos, católicos e a queda de Lula entre 2022 e 2026.
# Primeiro post da seção Análises de dados-eleicoes.danielmariani.com.br.
# Pergunta: Lula caiu mais nos territórios católicos ou evangélicos?
# Cada "degrau" refaz a comparação respondendo uma objeção:
#   1. cidades católicas × evangélicas (bruto)   → "é diferença de UF"
#   2. cidades da mesma UF                       → "é renda"
#   3. áreas de mesma renda, na mesma UF         → "evangélicas tinham piso"
#   4. pares que votaram IGUAL em 2022           → "é urbano × rural"
#   5. pares do mesmo tipo de área               → "é idade, porte..."
#   6. pares iguais também em idade e porte de cidade
# Em cada degrau: a queda média de CADA grupo, quantas unidades de cada
# grupo caíram/subiram, e a diferença entre os grupos.
# Dados: dois arquivos públicos do repositório (licença CC BY 4.0):
#   municipios_eleicao.parquet (resultado OFICIAL por município, TSE)
#   apond_completo.parquet     (áreas de ponderação: religião, renda etc.)

library(arrow)    # ler parquet
library(dplyr)    # manipulação
library(ggplot2)  # gráficos
library(svglite)  # salvar svg

BASE <- "https://dados-eleicoes.danielmariani.com.br/dados"

# cache local: baixa cada arquivo uma vez e reaproveita nas execuções
baixa <- function(arquivo) {
  dir.create("cache", showWarnings = FALSE)
  destino <- file.path("cache", arquivo)
  if (!file.exists(destino))
    download.file(file.path(BASE, arquivo), destino, mode = "wb")
  destino
}

municipios_of <- arrow::read_parquet(baixa("municipios_eleicao.parquet"))
apond <- arrow::read_parquet(baixa("apond_completo.parquet"))

# ---------------------------------------------------------------- tema ----
CAT <- "#B08A3E"   # territórios mais católicos (ocre)
EV  <- "#27568B"   # mais evangélicos (azul)
CINZA <- "#777777"

tema <- theme_minimal(base_size = 13) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = "#ececec", linewidth = .4),
        axis.title = element_text(size = 10.5, colour = "#555555"),
        axis.text = element_text(colour = "#555555"),
        plot.margin = margin(6, 12, 6, 6),
        legend.position = "none")

salva <- function(g, nome, h) {
  ggsave(file.path("fig", paste0(nome, ".svg")), g,
         width = 7.2, height = h, device = svglite::svglite)
  ggsave(file.path("fig", paste0(nome, ".png")), g,
         width = 7.2, height = h, dpi = 150)   # prova local
}

virgula <- function(x, dec = 1) sub("\\.", ",", sprintf("%+.*f", dec, x))

# resumo-padrão de um degrau: queda média, quantas caíram e subiram
resumo <- function(dados, grupo_col, delta_col, rotulo) {
  r <- dados |>
    group_by(grupo = .data[[grupo_col]]) |>
    summarise(n = n(), caiu = sum(.data[[delta_col]] < 0),
              subiu = sum(.data[[delta_col]] > 0),
              delta = mean(.data[[delta_col]]), .groups = "drop")
  cat(rotulo, "\n")
  for (i in seq_len(nrow(r)))
    cat(sprintf("  %-17s Δ %s p.p. | caiu em %d de %d (%.0f%%), subiu em %d\n",
                r$grupo[i], virgula(r$delta[i]), r$caiu[i], r$n[i],
                100 * r$caiu[i] / r$n[i], r$subiu[i]))
  cat("  diferença (católicas − evangélicas):",
      virgula(r$delta[r$grupo == "mais católicas"] -
              r$delta[r$grupo == "mais evangélicas"]), "p.p.\n\n")
  r
}

# ------------------------------------------------- bases de comparação ----
# Nível cidade: o resultado OFICIAL do TSE por município (sem estimativa).
# Religião vem da amostra do Censo: a menor geografia é a ÁREA DE
# PONDERAÇÃO (AP); a religião municipal é a média das APs (peso = pop).
voto_mun <- municipios_of |>
  select(cd_mun, nm_mun, uf, pct_2022 = pct_lula_2022t1,
         pct_2026 = pct_lula_2026t1, delta = delta_lula_t1) |>
  filter(!is.na(delta))

religiao_mun <- apond |>
  filter(!is.na(pct_catolica)) |>
  group_by(cd_mun) |>
  summarise(pct_catolica = weighted.mean(pct_catolica, pop),
            pct_evangelica = weighted.mean(pct_evangelica, pop),
            pop = sum(pop), .groups = "drop")

municipios <- inner_join(voto_mun, religiao_mun, by = "cd_mun") |>
  mutate(grupo = case_when(
    ntile(pct_catolica, 5) == 5   ~ "mais católicas",
    ntile(pct_evangelica, 5) == 5 ~ "mais evangélicas",
    TRUE ~ NA))
grupos_mun <- filter(municipios, !is.na(grupo))
cat(nrow(municipios), "municípios;",
    sum(grupos_mun$grupo == "mais católicas"), "+",
    sum(grupos_mun$grupo == "mais evangélicas"), "nos grupos\n\n")

# ------------------------------------------------------------ degrau 1 ----
d1 <- resumo(grupos_mun, "grupo", "delta",
             "degrau 1: todas as cidades (bruto)")

g1 <- ggplot(grupos_mun, aes(delta, grupo, colour = grupo)) +
  geom_vline(xintercept = 0, colour = "#999999") +
  geom_jitter(height = .22, size = .55, alpha = .18) +
  # média destacada: anel branco + rótulo com fundo (legível na nuvem)
  stat_summary(aes(fill = grupo), fun = mean, geom = "point", size = 5.4,
               shape = 21, colour = "white", stroke = 1.6) +
  stat_summary(aes(label = paste("média", virgula(after_stat(x)))),
               fun = mean, geom = "label", vjust = -0.75, size = 3.8,
               fontface = "bold", fill = "white", label.size = 0,
               label.padding = unit(0.18, "lines")) +
  scale_fill_manual(values = c("mais católicas" = CAT,
                               "mais evangélicas" = EV)) +
  scale_colour_manual(values = c("mais católicas" = CAT,
                                 "mais evangélicas" = EV)) +
  scale_x_continuous(labels = \(x) sprintf("%+.0f", x)) +
  labs(x = "variação do voto em Lula 2022 → 2026 (1º turno), em pontos percentuais",
       y = NULL) +
  tema + theme(axis.text.y = element_text(size = 12, face = "bold",
                                          colour = "#121212"))
salva(g1, "d1_cidades", 3.1)

# ------------------------------------------------------------ degrau 2 ----
# "É diferença de estado." Compara só cidades da MESMA UF.
d2 <- grupos_mun |>
  group_by(uf, grupo) |>
  summarise(delta = mean(delta), n = n(), .groups = "drop") |>
  tidyr::pivot_wider(names_from = grupo, values_from = c(delta, n)) |>
  filter(!is.na(`delta_mais católicas`), !is.na(`delta_mais evangélicas`)) |>
  mutate(gap = `delta_mais católicas` - `delta_mais evangélicas`,
         n_cidades = `n_mais católicas` + `n_mais evangélicas`)
grupos_mun_d2 <- semi_join(grupos_mun, d2, by = "uf")
r2 <- resumo(grupos_mun_d2, "grupo", "delta",
             "degrau 2: cidades da mesma UF (gap ponderado abaixo)")
cat("  gap médio ponderado dentro das UFs:",
    virgula(weighted.mean(d2$gap, d2$n_cidades)),
    "p.p. | católicas caíram mais em", sum(d2$gap < 0), "de",
    nrow(d2), "UFs\n\n")

d2g <- arrange(d2, gap) |> mutate(uf = factor(uf, uf))
g2 <- ggplot(d2g) +
  geom_vline(xintercept = 0, colour = "#999999") +
  geom_segment(aes(x = `delta_mais evangélicas`,
                   xend = `delta_mais católicas`, y = uf, yend = uf),
               colour = "#cccccc", linewidth = 1.6) +
  geom_point(aes(`delta_mais evangélicas`, uf), colour = EV, size = 2.6) +
  geom_point(aes(`delta_mais católicas`, uf), colour = CAT, size = 2.6) +
  annotate("text", x = d2g$`delta_mais católicas`[1], y = nrow(d2g) + 1.6,
           label = "cidades mais\ncatólicas", colour = CAT, size = 3.4,
           fontface = "bold", lineheight = .95) +
  annotate("text", x = d2g$`delta_mais evangélicas`[1] + 3.5,
           y = nrow(d2g) + 1.6, label = "cidades mais\nevangélicas",
           colour = EV, size = 3.4, fontface = "bold", lineheight = .95) +
  coord_cartesian(clip = "off") +
  scale_y_discrete(expand = expansion(add = c(.6, 2.4))) +
  scale_x_continuous(labels = \(x) sprintf("%+.0f", x)) +
  labs(x = "variação do voto em Lula 2022 → 2026, em pontos percentuais (média das cidades da UF)",
       y = NULL) + tema
salva(g2, "d2_ufs", 5.9)

# ------------------------------------------------------------ degrau 3 ----
# "É renda." Desce para as ÁREAS DE PONDERAÇÃO e compara só áreas da
# mesma UF e do mesmo quinto de renda per capita (quintil_renda_br).
aps <- apond |>
  filter(!is.na(pct_catolica), !is.na(delta_lula_t1),
         !is.na(renda_pc_media_sm)) |>
  mutate(grupo = case_when(
    ntile(pct_catolica, 5) == 5   ~ "mais católicas",
    ntile(pct_evangelica, 5) == 5 ~ "mais evangélicas",
    TRUE ~ NA)) |>
  filter(!is.na(grupo))

d3 <- aps |>
  group_by(uf, quintil_renda_br, grupo) |>
  summarise(delta = mean(delta_lula_t1), n = n(), .groups = "drop") |>
  tidyr::pivot_wider(names_from = grupo, values_from = c(delta, n)) |>
  filter(!is.na(`delta_mais católicas`), !is.na(`delta_mais evangélicas`)) |>
  mutate(gap = `delta_mais católicas` - `delta_mais evangélicas`,
         n_aps = `n_mais católicas` + `n_mais evangélicas`)
aps_d3 <- semi_join(aps, d3, by = c("uf", "quintil_renda_br"))
r3 <- resumo(aps_d3, "grupo", "delta_lula_t1",
             "degrau 3: áreas da mesma UF e renda")
cat("  gap médio ponderado nas", nrow(d3), "células UF × renda:",
    virgula(weighted.mean(d3$gap, d3$n_aps)),
    "p.p. | católicas caíram mais em", sum(d3$gap < 0), "células\n\n")

# por quinto de renda, para o gráfico (médias das células, ponderadas)
d3q <- d3 |>
  group_by(quintil_renda_br) |>
  summarise(catolicas = weighted.mean(`delta_mais católicas`, n_aps),
            evangelicas = weighted.mean(`delta_mais evangélicas`, n_aps)) |>
  tidyr::pivot_longer(-quintil_renda_br, names_to = "grupo",
                      values_to = "delta")
rotulo_q <- c("1\n20% mais\npobres", "2", "3", "4", "5\n20% mais\nricas")
g3 <- ggplot(d3q, aes(factor(quintil_renda_br), delta, fill = grupo)) +
  geom_hline(yintercept = 0, colour = "#999999") +
  geom_col(position = position_dodge(.62), width = .56) +
  geom_text(aes(label = virgula(delta)),
            position = position_dodge(.62), vjust = 1.35, size = 3.3,
            fontface = "bold", colour = "#333333") +
  annotate("text", x = .62, y = -6.9, label = "áreas mais católicas",
           colour = CAT, hjust = 0, size = 3.6, fontface = "bold") +
  annotate("text", x = .62, y = -7.6, label = "áreas mais evangélicas",
           colour = EV, hjust = 0, size = 3.6, fontface = "bold") +
  scale_fill_manual(values = c(catolicas = CAT, evangelicas = EV)) +
  scale_x_discrete(labels = rotulo_q) +
  scale_y_continuous(labels = \(x) sprintf("%+.0f", x),
                     expand = expansion(mult = c(.12, .05))) +
  labs(x = "quinto de renda domiciliar per capita da área (Brasil)",
       y = "variação de Lula 2022 → 2026, p.p.") +
  tema + theme(panel.grid.major.x = element_blank(),
               axis.text.x = element_text(lineheight = .95))
salva(g3, "d3_renda", 4.2)

# ------------------------------------------------- pareamentos (4 a 6) ----
# Função de pareamento: pares católica–evangélica da mesma UF que passam
# nos filtros ("calipers"); cada área entra em no máximo um par,
# começando pelos pares mais parecidos.
pop_mun <- apond |> group_by(cd_mun) |> summarise(pop_mun = sum(pop))
aps <- left_join(aps, pop_mun, by = "cd_mun")

pareia <- function(aps, filtro) {
  candidatos <- inner_join(
      filter(aps, grupo == "mais católicas"),
      filter(aps, grupo == "mais evangélicas"),
      by = "uf", suffix = c("_cat", "_ev"),
      relationship = "many-to-many") |>
    filtro() |>
    mutate(dist = scale(renda_pc_media_sm_cat - renda_pc_media_sm_ev)^2 +
                  scale(pct_lula_2022t1_cat - pct_lula_2022t1_ev)^2 +
                  scale(pct_pop_urbana_cat - pct_pop_urbana_ev)^2) |>
    arrange(dist)
  usada <- new.env(); escolhidos <- logical(nrow(candidatos))
  for (i in seq_len(nrow(candidatos))) {
    a <- candidatos$apond_cat[i]; b <- candidatos$apond_ev[i]
    if (is.null(usada[[a]]) && is.null(usada[[b]])) {
      usada[[a]] <- TRUE; usada[[b]] <- TRUE; escolhidos[i] <- TRUE
    }
  }
  candidatos[escolhidos, ] |>
    mutate(dif = delta_lula_t1_cat - delta_lula_t1_ev)
}

resumo_pares <- function(pares, rotulo) {
  cat(rotulo, "\n")
  for (lado in c("cat", "ev")) {
    d <- pares[[paste0("delta_lula_t1_", lado)]]
    cat(sprintf("  %-17s Δ %s p.p. | caiu em %d de %d (%.0f%%), subiu em %d\n",
                ifelse(lado == "cat", "mais católicas", "mais evangélicas"),
                virgula(mean(d)), sum(d < 0), length(d),
                100 * mean(d < 0), sum(d > 0)))
  }
  cat(sprintf("  diferença média (cat − ev): %s p.p. | católica caiu mais em %d de %d pares (%.0f%%)\n\n",
              virgula(mean(pares$dif)), sum(pares$dif < 0), nrow(pares),
              100 * mean(pares$dif < 0)))
}

# degrau 4: "evangélicas já votavam pouco em Lula, não tinham de onde
# cair" (efeito-piso). Pares de mesma renda (±20%) que votaram IGUAL em
# 2022 (±5 p.p.): mesmo ponto de partida.
p4 <- pareia(aps, \(d) filter(d,
  abs(log(renda_pc_media_sm_cat / renda_pc_media_sm_ev)) <= log(1.2),
  abs(pct_lula_2022t1_cat - pct_lula_2022t1_ev) <= 5))
resumo_pares(p4, "degrau 4: pares de mesma UF e renda que votaram igual em 2022")

# degrau 5: "é urbano contra rural": + urbanização parecida (±15 p.p.)
# e mesmo lado da linha urbana
p5 <- pareia(aps, \(d) filter(d,
  abs(log(renda_pc_media_sm_cat / renda_pc_media_sm_ev)) <= log(1.2),
  abs(pct_lula_2022t1_cat - pct_lula_2022t1_ev) <= 5,
  abs(pct_pop_urbana_cat - pct_pop_urbana_ev) <= 15,
  (pct_pop_urbana_cat >= 50) == (pct_pop_urbana_ev >= 50)))
resumo_pares(p5, "degrau 5: ... e mesmo tipo de área (urbana/rural)")

# degrau 6: "é idade, porte de cidade": + estrutura etária parecida
# (60+ ±5 p.p.) e cidades DIFERENTES de porte parecido (população ≤2×)
p6 <- pareia(aps, \(d) filter(d,
  cd_mun_cat != cd_mun_ev,
  pmax(pop_mun_cat, pop_mun_ev) / pmin(pop_mun_cat, pop_mun_ev) <= 2,
  abs(log(renda_pc_media_sm_cat / renda_pc_media_sm_ev)) <= log(1.2),
  abs(pct_lula_2022t1_cat - pct_lula_2022t1_ev) <= 5,
  abs(pct_pop_urbana_cat - pct_pop_urbana_ev) <= 15,
  (pct_pop_urbana_cat >= 50) == (pct_pop_urbana_ev >= 50),
  abs(pct_60mais_cat - pct_60mais_ev) <= 5))
resumo_pares(p6, "degrau 6: ... e mesma idade, cidades diferentes de mesmo porte")
cat("  (verificação do pareamento final) voto 2022:",
    round(mean(p6$pct_lula_2022t1_cat), 1), "vs",
    round(mean(p6$pct_lula_2022t1_ev), 1), "% | renda:",
    round(mean(p6$renda_pc_media_sm_cat), 2), "vs",
    round(mean(p6$renda_pc_media_sm_ev), 2), "SM | urbana:",
    round(mean(p6$pct_pop_urbana_cat)), "vs",
    round(mean(p6$pct_pop_urbana_ev)), "% | 60+:",
    round(mean(p6$pct_60mais_cat), 1), "vs",
    round(mean(p6$pct_60mais_ev), 1), "%\n\n")

g6 <- ggplot(p6, aes(dif, 0)) +
  geom_vline(xintercept = 0, colour = "#999999") +
  geom_jitter(height = .3, size = 1.2, alpha = .45, colour = "#555555") +
  annotate("point", x = mean(p6$dif), y = 0, colour = CAT, size = 5) +
  annotate("text", x = mean(p6$dif), y = .5,
           label = paste("média:", virgula(mean(p6$dif)), "p.p."),
           colour = CAT, fontface = "bold", size = 3.7) +
  annotate("text", x = min(p6$dif), y = -.52, hjust = 0, size = 3.4,
           colour = CINZA, label = sprintf(
             "em %d%% dos pares, Lula caiu mais na área mais católica",
             round(100 * mean(p6$dif < 0))), lineheight = .95) +
  scale_x_continuous(labels = \(x) sprintf("%+.0f", x)) +
  scale_y_continuous(limits = c(-.68, .68)) +
  labs(x = "diferença da variação dentro do par (católica − evangélica), em pontos percentuais",
       y = NULL) +
  tema + theme(axis.text.y = element_blank(),
               panel.grid.major.y = element_blank())
salva(g6, "d6_pares", 2.9)

# ------------------------------------------------------------- escada -----
rotulos <- c("todas as cidades",
             "cidades da mesma UF",
             "áreas da mesma UF e renda",
             "pares que votaram igual em 2022",
             "... e mesmo tipo (urbana/rural)",
             "... e mesma idade e porte de cidade")
escada <- tibble::tibble(
  degrau = factor(rotulos, levels = rev(rotulos)),
  gap = c(d1$delta[d1$grupo == "mais católicas"] -
            d1$delta[d1$grupo == "mais evangélicas"],
          weighted.mean(d2$gap, d2$n_cidades),
          weighted.mean(d3$gap, d3$n_aps),
          mean(p4$dif), mean(p5$dif), mean(p6$dif)))
print(escada)

g7 <- ggplot(escada, aes(gap, degrau)) +
  geom_col(width = .58, fill = c(rep("#c9c3b8", 5), CAT)) +
  geom_text(aes(label = virgula(gap)), hjust = 1.15, size = 3.9,
            fontface = "bold", colour = "#333333") +
  geom_vline(xintercept = 0, colour = "#999999") +
  scale_x_continuous(expand = expansion(mult = c(.2, .02))) +
  labs(x = "pontos percentuais a mais de queda nos territórios católicos",
       y = NULL) +
  tema + theme(panel.grid.major.y = element_blank(),
               axis.text.y = element_text(size = 11, hjust = 1,
                                          colour = "#121212"))
salva(g7, "d7_escada", 3.6)

# números do post
jsonlite::write_json(list(
  municipios = nrow(municipios), d1 = d1, d2 = r2,
  d2_ufs = nrow(d2), d2_gap = weighted.mean(d2$gap, d2$n_cidades),
  d2_placar = sum(d2$gap < 0),
  d3 = r3, d3_cel = nrow(d3), d3_gap = weighted.mean(d3$gap, d3$n_aps),
  d3_placar = sum(d3$gap < 0), d3q = d3q,
  p4 = list(n = nrow(p4), gap = mean(p4$dif),
            placar = sum(p4$dif < 0)),
  p5 = list(n = nrow(p5), gap = mean(p5$dif),
            placar = sum(p5$dif < 0)),
  p6 = list(n = nrow(p6), gap = mean(p6$dif),
            placar = sum(p6$dif < 0)),
  escada = escada), "numeros.json", digits = 4, pretty = TRUE)
cat("figuras em fig/, números em numeros.json\n")