## =============================================================================
## Workshop CIPE | FGV EAESP | 19 de agosto de 2026
##
## Pedro Borrego
##
## Exercício R: dimensão da amostra (amostragem aleatória simples e por
## conglomerados) e análise de potência, para a estimação de uma média
##
## Mesmo caso de estudo dos Exercícios Excel 1 e 2 (fictício, inspirado em:
## Bruhn, M., Garber, G., Koyama, S., Zia, B. (2022). The Long-Term Impact of
## High School Financial Education: evidence from Brazil. WP 563.
## Banco Central do Brasil Working Papers):
## uma rede estadual com N = 45 000 alunos da 3ª série, em M = 300 escolas,
## quer estimar a pontuação média (escala 0-100) numa prova de literacia financeira.
## Nota: o estudo original é um estudo experimental (aleartorizado e controlado) 
## que compara tratamento vs. controlo, mas neste exercício apenas se
## estima um parâmetro único (a média), servindo o estudo só de inspiração
## para o contexto.
##
##
## Só usa packages base do R para os cálculos estatísticos, para aplicar as
## fórmulas dos slides. Em alternativa, podem ser explorados packages como
## 'pwr' ou 'powertools'.
##
##
## =============================================================================


## -----------------------------------------------------------------------
## Carregar pacotes necessários
## -----------------------------------------------------------------------

library(ggplot2)   # para gráficos
# se não tiver instalado, usar install.packages("ggplot2") na primeira vez que 
# correr este script, mas depois já náão vair ser necessário instalar novamente


## -----------------------------------------------------------------------
## 0. Parâmetros gerais do caso de estudo (comuns aos Exercícios 1 e 2)
## -----------------------------------------------------------------------
N     <- 45000   # populacao: n. de alunos do 3.o ano
M     <- 300     # populacao: n. de escolas
conf  <- 0.95     # nivel de confianca
alpha <- 1 - conf
z_conf <- qnorm(1 - alpha / 2)   # z(1-alfa/2) -- usa-se qnorm em vez de NORMSINV do Excel
d     <- 3.0      # margem de erro absoluta desejada (pontos, escala 0-100)

cat("z(1-alfa/2) para", conf * 100, "% de confianca:", round(z_conf, 3), "\n\n")


## =========================================================================
## PARTE 1 -- Dimensao da amostra: amostragem aleatoria simples sem reposição
## (replica o Exercicio Excel 1, para verificação cruzada)
## =========================================================================

# Estimativa de sigma -- Metodo A: regra do intervalo/6 (mais conservador)
amplitude_esperada <- 100
sigma_A <- amplitude_esperada / 6

# Estimativa de sigma -- Metodo B: amostra piloto (25 alunos)
piloto_individual <- c(49.6, 49.2, 50.2, 63.2, 50.0, 28.0, 57.4, 47.8, 48.6, 53.8,
                        55.8, 70.6, 62.6, 53.8, 40.2, 35.8, 56.0, 73.0, 52.6, 50.2,
                        60.6, 28.8, 47.0, 59.8, 66.0)
sigma_B <- sd(piloto_individual)   # sd() em R usa por omissao o divisor (n-1), tal como STDEV.S no Excel

cat("sigma (Metodo A, intervalo/6):", round(sigma_A, 2), "\n")
cat("sigma (Metodo B, amostra piloto, n=25):", round(sigma_B, 2), "\n\n")


# Cálculo de dimensão da amostra com sigma_A (desvio padrão obtido no Método A)
z <- qnorm(1 - (1 - conf) / 2)
n_aas<-(z * sigma_A / d)^2 / (1 + (z^2 * sigma_A^2 / (d^2*N)))
n_ass
# Arredonda para o número intereiro imediatamente acima
n_aas_final <- ceiling(n_aas)

cat("\n--> Dimensao de amostra (AAS, Metodo A):", n_aas_final, "alunos\n")
cat("    (no Exercicio Excel 1, para os mesmos pressupostos, o resultado foi 119 alunos)\n\n")


## =========================================================================
## PARTE 2 -- Dimensao da amostra: amostragem por conglomerados
## (Tiragem de clusters com probabilidades iguais e sem reposição)
## (replica o Exercicio Excel 2, para verificação cruzada)
## =========================================================================

# Amostra piloto ao nivel da escola: 12 escolas, todos os alunos avaliados
piloto_escolas <- data.frame(
  escola   = 1:12,
  N_g      = c(145, 150, 155, 142, 152, 148, 158, 144, 150, 156, 146, 151),
  media_g  = c(48.2, 55.6, 62.4, 39.0, 51.2, 44.0, 58.8, 41.6, 53.0, 66.2, 47.4, 56.0)
)
piloto_escolas$tau_g <- piloto_escolas$N_g * piloto_escolas$media_g

var_tau <- var(piloto_escolas$tau_g)   # var() no R usa por omissao o divisor (n-1) -- variância corrigida
cat("Amostra piloto (escolas):\n")
print(piloto_escolas)
cat("\nVariancia corrigida dos totais por escola (s^2_tau):", round(var_tau, 1), "\n\n")

# Funcao: dimensao da amostra por conglomerados
m_clusters <- function(var_tau, d, M, N, conf = 0.95) {
  z <- qnorm(1 - (1 - conf) / 2)
  m <- (z^2 * M^2 * var_tau) / (d^2 * N^2 + z^2 * M * var_tau)
  alunos_esperados <- ceiling(m) * (N / M)
  data.frame(m = m, m_arredondado = ceiling(m), alunos_esperados = alunos_esperados)
}

resultado_clusters <- m_clusters(var_tau, d, M, N, conf)
print(resultado_clusters)
cat("\n--> Dimensao de amostra (clusters):", resultado_clusters$m_arredondado,
    "escolas, ~", resultado_clusters$alunos_esperados, "alunos\n")
cat("    (no Exercicio Excel 2, para os mesmos pressupostos, o resultado foi 37 escolas / 5550 alunos)\n\n")

cat("COMPARAÇÃO (mesma precisao, d =", d, "pontos, confianca =", conf * 100, "%):\n")
cat("  Aleatoria simples: ", n_aas_final, " alunos\n", sep = "")
cat("  Conglomerados:     ", resultado_clusters$alunos_esperados,
    " alunos (", resultado_clusters$m_arredondado, " escolas)\n\n", sep = "")


## =========================================================================
## PARTE 3 -- Analise de potencia: teste-z bilateral para uma média
## 
## =========================================================================

# Funcao de potencia, teste bilateral, para H0: mu = mu0 vs H1: mu = mu0 + delta
# Potencia = 1 - beta = P(rejeitar H0 | H1 verdadeira)
potencia_z <- function(n, delta, sigma, alpha = 0.05) {
  z_alpha <- qnorm(1 - alpha / 2)
  termo <- sqrt(n) * abs(delta) / sigma
  pnorm(termo - z_alpha) + pnorm(-termo - z_alpha)
}

# Funcao inversa: n necessario para atingir uma potencia-alvo (1-beta), dado delta e sigma
n_para_potencia <- function(potencia_alvo, delta, sigma, alpha = 0.05) {
  z_alpha <- qnorm(1 - alpha / 2)
  z_beta  <- qnorm(potencia_alvo)
  n <- ((z_alpha + z_beta) * sigma / delta)^2
  ceiling(n)
}

# --- Verificação cruzada com o Exercicio Excel 1 (separador Potencia) ---
sigma_potencia <- sigma_A     # 16.67, tal como no Excel
delta_potencia <- 4           # diferenca minima de interesse (pontos)
alpha_potencia <- 0.05

potencia_119 <- potencia_z(n_aas_final, delta_potencia, sigma_potencia, alpha_potencia)
n_para_80pc  <- n_para_potencia(0.80, delta_potencia, sigma_potencia, alpha_potencia)

cat("Verificação cruzada com o Exercicio Excel 1 (potencia):\n")
cat("  Potencia com n =", n_aas_final, "e delta =", delta_potencia, "pontos:",
    round(potencia_119 * 100, 1), "%  (Excel: 74.5%)\n")
cat("  n necessario para 80% de potencia:", n_para_80pc, " (Excel: 137)\n\n")


## -----------------------------------------------------------------------
## 3.1 Explorar potencia vs. dimensao da amostra, para varios efeitos (Delta)
##
## -----------------------------------------------------------------------

n_grid <- seq(20, 400, by = 5)

# Tres cenarios de efeito minimo a detetar, em pontos (escala 0-100):
# - Delta = 2   : diferenca pequena
# - Delta = 4   : diferenca moderada (o mesmo cenario do exercício do Excel)
# - Delta = 0.25 * sigma_A : dimensao do efeito "pequeno a medio" de Cohen,
#   e tambem ordem de grandeza do efeito de 0.25 desvio-padrao encontrado por Bruhn et al. (2022)
#   no proprio programa de educacao financeira que inspira este exercicio
deltas <- c(2, 4, 0.25 * sigma_A)
labels_deltas <- c("Delta = 2 pontos", "Delta = 4 pontos",
                    paste0("Delta = 0,25 sigma"))

potencia_matriz <- sapply(deltas, function(d) potencia_z(n_grid, d, sigma_A, alpha))


# Gráfico com o package ggplot2
library(ggplot2)
grafico1_potencia_vs_n <- ggplot(data = data.frame(n = n_grid), aes(x = n)) +
  geom_line(aes(y = potencia_matriz[, 1], color = labels_deltas[1]), linewidth = 1) +
  geom_line(aes(y = potencia_matriz[, 2], color = labels_deltas[2]), linewidth = 1) +
  geom_line(aes(y = potencia_matriz[, 3], color = labels_deltas[3]), linewidth = 1) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = n_aas_final, linetype = "dotted", color = "#1F4E78") +
  annotate("text", x = max(n_grid) * 0.98, y = 0.82, label = "potencia = 80%", color = "grey40", hjust = 1, size = 3) +
  labs(title = "Potencia vs. n da amostra, para diferentes efeitos minimos",
       x = "Dimensao da amostra (n)", y = "Potencia (1 - beta)", color = "Cenarios de Delta") +
  scale_color_manual(values = c("#1F4E78", "#C00000", "#548235")) +
  theme_minimal()
# Visualiza o gráfico
grafico1_potencia_vs_n



# n necessario para 80% de potencia, em cada cenario de Delta
cat("\nn necessario para 80% de potencia, por cenario de Delta:\n")
for (i in seq_along(deltas)) {
  n_needed <- n_para_potencia(0.80, deltas[i], sigma_A, alpha)
  cat("  ", labels_deltas[i], ": n =", n_needed, "\n")
}


## -----------------------------------------------------------------------
## 3.2 Explorar o trade-off entre nivel de significancia (alfa) e potencia
##
## -----------------------------------------------------------------------

alphas <- c(0.10, 0.05, 0.01)
labels_alphas <- paste0("alfa = ", alphas)
delta_fixo <- 4   # pontos

potencia_matriz_alpha <- sapply(alphas, function(a) potencia_z(n_grid, delta_fixo, sigma_A, a))


# Gráfico com o package ggplot2
grafico2_potencia_vs_alfa <- ggplot(data = data.frame(n = n_grid), aes(x = n)) +
  geom_line(aes(y = potencia_matriz_alpha[, 1], color = labels_alphas[1]), linewidth = 1) +
  geom_line(aes(y = potencia_matriz_alpha[, 2], color = labels_alphas[2]), linewidth = 1) +
  geom_line(aes(y = potencia_matriz_alpha[, 3], color = labels_alphas[3]), linewidth = 1) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = n_aas_final, linetype = "dotted", color = "#1F4E78") +
  annotate("text", x = max(n_grid) * 0.98, y = 0.82, label = "potencia = 80%", color = "grey40", hjust = 1, size = 3) +
  labs(title = paste0("Trade-off entre alfa e potencia (Delta = ", delta_fixo, " pontos)"),
       x = "Dimensao da amostra (n)", y = "Potencia (1 - beta)", color = "Nivel de significancia") +
  scale_color_manual(values = c("#1F4E78", "#C00000", "#548235")) +
  theme_minimal()
# Visualiza o gráfico
grafico2_potencia_vs_alfa

## -----------------------------------------------------------------------
## Fim do script
## -----------------------------------------------------------------------