library(pwr)
library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(car)
1. Importação e Limpeza dos Dados
petro_prod_total <- read_csv(
"C:/Users/Guilherme/Downloads/petro/oil-production-by-country.csv"
) |>
rename(
pais = Entity,
codigo = Code,
ano = Year,
producao_oil = Oil
) |>
drop_na() |>
filter(!is.na(codigo) & codigo != "")
gini_index <- read_csv(
"C:/Users/Guilherme/Downloads/GINI/API_SI.POV.GINI_DS2_en_csv_v2_288.csv",
skip = 4
) |>
select(-"Indicator Name", -"Indicator Code") |>
pivot_longer(
cols = -c("Country Name", "Country Code"),
names_to = "ano",
values_to = "gini"
) |>
rename(pais = "Country Name", codigo = "Country Code") |>
mutate(ano = as.integer(ano)) |>
filter(ano <= 2023) |>
drop_na(gini)
2. Filtro 2015 — Países em Comum
paises_comuns_gini <- inner_join(
gini_index |> filter(ano == 2015) |> select(codigo),
petro_prod_total |> filter(ano == 2015) |> select(codigo),
by = "codigo"
)
gini_2015 <- gini_index |> filter(ano == 2015, codigo %in% paises_comuns_gini$codigo)
petro_gini_2015 <- petro_prod_total |> filter(ano == 2015, codigo %in% paises_comuns_gini$codigo)
3. Teste de Hipótese — GINI vs Petróleo
dados_gini <- inner_join(gini_2015, petro_gini_2015, by = "codigo")
mediana_producao <- median(dados_gini$producao_oil)
dados_gini <- dados_gini |>
mutate(grupo = ifelse(producao_oil >= mediana_producao,
"Acima da mediana",
"Abaixo da mediana"))
levene_result <- leveneTest(gini ~ grupo, data = dados_gini)
levene_result
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3333 0.2515
## 84
var_igual <- levene_result$`Pr(>F)`[1] > 0.05
cat("Variâncias iguais?", ifelse(var_igual, "SIM → teste t padrão", "NÃO → Welch t-test"), "\n")
## Variâncias iguais? SIM → teste t padrão
table(dados_gini$grupo)
##
## Abaixo da mediana Acima da mediana
## 43 43
ggplot(dados_gini, aes(x = gini, fill = grupo)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = median(dados_gini$gini[dados_gini$grupo == "Acima da mediana"]),
color = "blue", linetype = "dashed") +
geom_vline(xintercept = median(dados_gini$gini[dados_gini$grupo == "Abaixo da mediana"]),
color = "red", linetype = "dashed") +
scale_x_continuous(limits = c(0, 80)) +
labs(
title = "Distribuição do GINI por Grupo de Produção de Petróleo (2015)",
x = "Índice GINI",
y = "Densidade",
fill = "Grupo"
) +
theme_minimal()
ggplot(dados_gini, aes(x = gini, fill = grupo)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = mean(dados_gini$gini[dados_gini$grupo == "Acima da mediana"]),
color = "blue", linetype = "dashed", show.legend = FALSE) +
geom_vline(xintercept = mean(dados_gini$gini[dados_gini$grupo == "Abaixo da mediana"]),
color = "red", linetype = "dashed", show.legend = FALSE) +
facet_wrap(~ grupo, ncol = 1) +
scale_x_continuous(limits = c(0, 80)) +
labs(
title = "Distribuição do GINI por Grupo (2015)",
x = "Índice GINI",
y = "Densidade",
fill = "Grupo"
) +
theme_minimal()
4. Estatísticas e Teste t
dados_gini |>
group_by(grupo) |>
summarise(
n = n(),
media_gini = mean(gini),
dp_gini = sd(gini)
)
## # A tibble: 2 × 4
## grupo n media_gini dp_gini
## <chr> <int> <dbl> <dbl>
## 1 Abaixo da mediana 43 37.7 8.52
## 2 Acima da mediana 43 35.9 6.81
resultado_teste_gini <- t.test(
gini ~ grupo,
data = dados_gini,
var.equal = var_igual
)
resultado_teste_gini
##
## Two Sample t-test
##
## data: gini by grupo
## t = 1.0749, df = 84, p-value = 0.2855
## alternative hypothesis: true difference in means between group Abaixo da mediana and group Acima da mediana is not equal to 0
## 95 percent confidence interval:
## -1.520062 5.096806
## sample estimates:
## mean in group Abaixo da mediana mean in group Acima da mediana
## 37.66977 35.88140
alfa <- 0.05
cat("Erro Tipo 1 (alfa):", alfa, "\n")
## Erro Tipo 1 (alfa): 0.05
media_acima <- mean(dados_gini$gini[dados_gini$grupo == "Acima da mediana"])
media_abaixo <- mean(dados_gini$gini[dados_gini$grupo == "Abaixo da mediana"])
dp_acima <- sd(dados_gini$gini[dados_gini$grupo == "Acima da mediana"])
dp_abaixo <- sd(dados_gini$gini[dados_gini$grupo == "Abaixo da mediana"])
dp_pooled <- sqrt((dp_acima^2 + dp_abaixo^2) / 2)
d_cohen <- (media_acima - media_abaixo) / dp_pooled
gl <- 2 * 49 - 2
t_critico <- qt(1 - 0.05 / 2, df = gl)
cat("Graus de liberdade:", gl, "\n")
## Graus de liberdade: 96
cat("Valor crítico (±t):", round(t_critico, 4), "\n")
## Valor crítico (±t): 1.985
cat("t observado:", round(resultado_teste_gini$statistic, 4), "\n")
## t observado: 1.0749
cat("p-valor:", round(resultado_teste_gini$p.value, 4), "\n")
## p-valor: 0.2855
cat("Rejeita H0 se t <", round(-t_critico, 4), "ou t >", round(t_critico, 4), "\n")
## Rejeita H0 se t < -1.985 ou t > 1.985
cat("t observado está na região de",
ifelse(abs(resultado_teste_gini$statistic) > t_critico, "REJEIÇÃO ❌", "NÃO REJEIÇÃO ✅"), "\n")
## t observado está na região de NÃO REJEIÇÃO ✅
x <- seq(-4, 4, length.out = 1000)
y <- dt(x, df = gl)
ggplot(data.frame(x, y), aes(x, y)) +
geom_line(size = 1) +
geom_area(data = subset(data.frame(x, y), x > t_critico), aes(x, y), fill = "red", alpha = 0.4) +
geom_area(data = subset(data.frame(x, y), x < -t_critico), aes(x, y), fill = "red", alpha = 0.4) +
geom_vline(xintercept = t_critico, color = "red", linetype = "dashed") +
geom_vline(xintercept = -t_critico, color = "red", linetype = "dashed") +
geom_vline(xintercept = resultado_teste_gini$statistic, color = "blue", size = 1) +
annotate("text", x = t_critico + 0.3, y = 0.35, label = paste("t crítico =", round( t_critico, 2)), color = "red") +
annotate("text", x = -t_critico - 0.3, y = 0.35, label = paste("t crítico =", round(-t_critico, 2)), color = "red") +
annotate("text", x = resultado_teste_gini$statistic, y = 0.38,
label = paste("t obs =", round(resultado_teste_gini$statistic, 2)), color = "blue") +
labs(
title = "Distribuição t — Teste GINI vs Produção de Petróleo (2015)",
x = "Valor t",
y = "Densidade"
) +
theme_minimal()
5. Teste Eliminando Não Produtores
Amostra_petro2 <- petro_gini_2015 |> filter(producao_oil > 0)
paises_comuns_gini_2 <- inner_join(
gini_2015 |> select(codigo),
Amostra_petro2 |> select(codigo),
by = "codigo"
)
dados_gini_2 <- inner_join(gini_2015, Amostra_petro2, by = "codigo")
mediana_producao_2 <- median(dados_gini_2$producao_oil)
dados_gini_2 <- dados_gini_2 |>
mutate(grupo = ifelse(producao_oil >= mediana_producao_2,
"Acima da mediana",
"Abaixo da mediana"))
dados_gini_2 |>
group_by(grupo) |>
summarise(
n = n(),
media_gini = mean(gini),
dp_gini = sd(gini)
)
## # A tibble: 2 × 4
## grupo n media_gini dp_gini
## <chr> <int> <dbl> <dbl>
## 1 Abaixo da mediana 25 34.0 5.82
## 2 Acima da mediana 26 36.9 7.17
resultado_teste_gini_2 <- t.test(
gini ~ grupo,
data = dados_gini_2,
var.equal = FALSE
)
resultado_teste_gini_2
##
## Welch Two Sample t-test
##
## data: gini by grupo
## t = -1.6148, df = 47.683, p-value = 0.1129
## alternative hypothesis: true difference in means between group Abaixo da mediana and group Acima da mediana is not equal to 0
## 95 percent confidence interval:
## -6.6184869 0.7231023
## sample estimates:
## mean in group Abaixo da mediana mean in group Acima da mediana
## 33.96000 36.90769
alfa_2 <- 0.05
cat("Erro Tipo 1 (alfa):", alfa_2, "\n")
## Erro Tipo 1 (alfa): 0.05
media_acima <- mean(dados_gini_2$gini[dados_gini_2$grupo == "Acima da mediana"])
media_abaixo <- mean(dados_gini_2$gini[dados_gini_2$grupo == "Abaixo da mediana"])
dp_acima <- sd(dados_gini_2$gini[dados_gini_2$grupo == "Acima da mediana"])
dp_abaixo <- sd(dados_gini_2$gini[dados_gini_2$grupo == "Abaixo da mediana"])
dp_pooled <- sqrt((dp_acima^2 + dp_abaixo^2) / 2)
d_cohen_2 <- (media_acima - media_abaixo) / dp_pooled
n_grupo <- min(table(dados_gini_2$grupo))
gl <- 2 * n_grupo - 2
t_critico <- qt(1 - 0.05 / 2, df = gl)
cat("Graus de liberdade:", gl, "\n")
## Graus de liberdade: 48
cat("Valor crítico (±t):", round(t_critico, 4), "\n")
## Valor crítico (±t): 2.0106
cat("t observado:", round(resultado_teste_gini_2$statistic, 4), "\n")
## t observado: -1.6148
cat("p-valor:", round(resultado_teste_gini_2$p.value, 4), "\n")
## p-valor: 0.1129
cat("Rejeita H0 se t <", round(-t_critico, 4), "ou t >", round(t_critico, 4), "\n")
## Rejeita H0 se t < -2.0106 ou t > 2.0106
cat("t observado está na região de",
ifelse(abs(resultado_teste_gini_2$statistic) > t_critico, "REJEIÇÃO", "NÃO REJEIÇÃO"), "\n")
## t observado está na região de NÃO REJEIÇÃO
ggplot(dados_gini_2, aes(x = gini, fill = grupo)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = median(dados_gini_2$gini[dados_gini_2$grupo == "Acima da mediana"]),
color = "blue", linetype = "dashed", show.legend = FALSE) +
geom_vline(xintercept = median(dados_gini_2$gini[dados_gini_2$grupo == "Abaixo da mediana"]),
color = "red", linetype = "dashed", show.legend = FALSE) +
scale_x_continuous(limits = c(0, 80)) +
labs(
title = "Distribuição do GINI — Apenas Produtores de Petróleo (2015)",
x = "Índice GINI",
y = "Densidade",
fill = "Grupo"
) +
theme_minimal()
ggplot(dados_gini_2, aes(x = gini, fill = grupo)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = mean(dados_gini_2$gini[dados_gini_2$grupo == "Acima da mediana"]),
color = "blue", linetype = "dashed", show.legend = FALSE) +
geom_vline(xintercept = mean(dados_gini_2$gini[dados_gini_2$grupo == "Abaixo da mediana"]),
color = "red", linetype = "dashed", show.legend = FALSE) +
facet_wrap(~ grupo, ncol = 1) +
scale_x_continuous(limits = c(0, 80)) +
labs(
title = "Distribuição do GINI — Apenas Produtores de Petróleo (2015)",
x = "Índice GINI",
y = "Densidade",
fill = "Grupo"
) +
theme_minimal()
cat("\nTamanho do efeito (d de Cohen):", round(d_cohen_2, 4), "\n")
##
## Tamanho do efeito (d de Cohen): 0.4514
cat("╔══════════════════════════════════════════╗\n")
## ╔══════════════════════════════════════════╗
cat("║ RESULTADO DO TESTE t ║\n")
## ║ RESULTADO DO TESTE t ║
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ t observado :", round(resultado_teste_gini_2$statistic, 4), "\n")
## ║ t observado : -1.6148
cat("║ t crítico : ±", round(t_critico, 4), "\n")
## ║ t crítico : ± 2.0106
cat("║ p-valor :", round(resultado_teste_gini_2$p.value, 4), "\n")
## ║ p-valor : 0.1129
cat("║ alfa : 0.05\n")
## ║ alfa : 0.05
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ d de Cohen :", round(d_cohen_2, 4), "\n")
## ║ d de Cohen : 0.4514
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ DECISÃO:",
ifelse(resultado_teste_gini_2$p.value < 0.05,
"REJEITA H0 ❌",
"NÃO REJEITA H0 ✅"), "\n")
## ║ DECISÃO: NÃO REJEITA H0 ✅
cat("╚══════════════════════════════════════════╝\n")
## ╚══════════════════════════════════════════╝