1. Importação e Limpeza dos Dados

Poluição do Ar (PM2.5)

air_pollution <- read_csv2(
  "C:/Users/Guilherme/Downloads/air polution/base_com_continentes.csv"
) |>
  select(-"Indicator Name", -"Indicator Code", -"Unnamed: 70") |>
  pivot_longer(
    cols      = -c("Continente", "Country Name", "Country Code"),
    names_to  = "ano",
    values_to = "pm25"
  ) |>
  rename(continente = Continente, pais = `Country Name`, codigo = `Country Code`) |>
  mutate(ano = as.integer(ano)) |>
  filter(ano == 2018) |>
  drop_na(pm25) |>
  filter(continente != "Outros")

air_pollution |> count(continente) |> kable(caption = "Países por Continente")
Países por Continente
continente n
Africa 54
America 32
Asia 46
Europe 43
Oceania 14

Gasto com Saúde

gasto_saude <- read_csv(
  "C:/Users/Guilherme/Downloads/petro/Gasto com saúde.csv",
  skip = 4
) |>
  select(-"Indicator Name", -"Indicator Code") |>
  pivot_longer(cols = -c("Country Name", "Country Code"), names_to = "ano", values_to = "gasto_saude_pct_gdp") |>
  rename(pais = "Country Name", codigo = "Country Code") |>
  mutate(ano = as.integer(ano)) |>
  filter(ano == 2018) |>
  drop_na(gasto_saude_pct_gdp)

Percentual de Idosos (65+)

pop_idosa <- read_csv(
  "C:/Users/Guilherme/Downloads/% 65+/Population (age group as % of total population).csv"
) |>
  rename(
    pais        = Economy,
    ano         = Year,
    codigo      = `Economy Code`,
    pct_65_mais = `Population ages 65 and above (% of total population)`
  ) |>
  drop_na() |>
  filter(!is.na(codigo) & codigo != "")

Densidade Populacional

pop_densidade <- read_csv(
  "C:/Users/Guilherme/Downloads/densidade populacional/API_EN.POP.DNST_DS2_en_csv_v2_1453.csv",
  skip = 4
) |>
  select(-"Indicator Name", -"Indicator Code") |>
  pivot_longer(cols = -c("Country Name", "Country Code"), names_to = "ano", values_to = "densidade_pop") |>
  rename(pais = "Country Name", codigo = "Country Code") |>
  mutate(ano = as.integer(ano)) |>
  filter(ano == 2023) |>
  drop_na(densidade_pop) |>
  filter(!is.na(codigo) & codigo != "")

2. Montagem da Base

paises_comuns_air <- inner_join(
  air_pollution |> select(codigo),
  gasto_saude   |> select(codigo),
  by = "codigo"
)

Air_2018       <- air_pollution |> filter(codigo %in% paises_comuns_air$codigo)
Health_2018    <- gasto_saude   |> filter(codigo %in% paises_comuns_air$codigo)
dados_trabalho <- inner_join(Air_2018, Health_2018, by = "codigo")

mediana_pm25 <- median(dados_trabalho$pm25)

dados_trabalho <- dados_trabalho |>
  mutate(grupo = ifelse(pm25 >= mediana_pm25, "Alta poluição", "Baixa poluição"))

table(dados_trabalho$grupo)
## 
##  Alta poluição Baixa poluição 
##             93             93

3. Teste de Hipótese

H0: A média do gasto com saúde é igual nos dois grupos de poluição. H1: A média do gasto com saúde é diferente nos dois grupos de poluição.

Teste de Levene (Homogeneidade de Variâncias)

levene_result <- leveneTest(gasto_saude_pct_gdp ~ grupo, data = dados_trabalho)
levene_result
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   1  2.9271 0.08879 .
##       184                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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

Teste t

resultado_teste <- t.test(
  gasto_saude_pct_gdp ~ grupo,
  data      = dados_trabalho,
  var.equal = var_igual
)
resultado_teste
## 
##  Two Sample t-test
## 
## data:  gasto_saude_pct_gdp by grupo
## t = -5.427, df = 184, p-value = 1.791e-07
## alternative hypothesis: true difference in means between group Alta poluição and group Baixa poluição is not equal to 0
## 95 percent confidence interval:
##  -2.933378 -1.369211
## sample estimates:
##  mean in group Alta poluição mean in group Baixa poluição 
##                     5.328214                     7.479508

Estatísticas Descritivas

dados_trabalho |>
  group_by(grupo) |>
  summarise(
    n           = n(),
    media_saude = round(mean(gasto_saude_pct_gdp), 2),
    dp_saude    = round(sd(gasto_saude_pct_gdp), 2),
    media_pm25  = round(mean(pm25), 2)
  ) |>
  kable(caption = "Estatísticas Descritivas por Grupo")
Estatísticas Descritivas por Grupo
grupo n media_saude dp_saude media_pm25
Alta poluição 93 5.33 2.44 38.60
Baixa poluição 93 7.48 2.95 13.93

Poder do Teste

media_alta  <- mean(dados_trabalho$gasto_saude_pct_gdp[dados_trabalho$grupo == "Alta poluição"])
media_baixa <- mean(dados_trabalho$gasto_saude_pct_gdp[dados_trabalho$grupo == "Baixa poluição"])
dp_alta     <- sd(dados_trabalho$gasto_saude_pct_gdp[dados_trabalho$grupo == "Alta poluição"])
dp_baixa    <- sd(dados_trabalho$gasto_saude_pct_gdp[dados_trabalho$grupo == "Baixa poluição"])
dp_pooled   <- sqrt((dp_alta^2 + dp_baixa^2) / 2)
d_cohen     <- (media_alta - media_baixa) / dp_pooled

n_grupo   <- min(table(dados_trabalho$grupo))
gl        <- 2 * n_grupo - 2
t_critico <- qt(1 - 0.05 / 2, df = gl)

poder <- pwr.t.test(n = n_grupo, d = d_cohen, sig.level = 0.05, type = "two.sample")

cat("╔══════════════════════════════════════════╗\n")
## ╔══════════════════════════════════════════╗
cat("║     RESULTADO — POLUIÇÃO vs SAÚDE        ║\n")
## ║     RESULTADO — POLUIÇÃO vs SAÚDE        ║
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ Levene p-valor  :", round(levene_result$`Pr(>F)`[1], 4), "\n")
## ║ Levene p-valor  : 0.0888
cat("║ Variâncias      :", ifelse(var_igual, "Iguais", "Diferentes"), "\n")
## ║ Variâncias      : Iguais
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ t observado     :", round(resultado_teste$statistic, 4), "\n")
## ║ t observado     : -5.427
cat("║ t crítico       : ±", round(t_critico, 4), "\n")
## ║ t crítico       : ± 1.9729
cat("║ p-valor         :", round(resultado_teste$p.value, 4), "\n")
## ║ p-valor         : 0
cat("║ alfa            : 0.05\n")
## ║ alfa            : 0.05
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ d de Cohen      :", round(d_cohen, 4), "\n")
## ║ d de Cohen      : -0.7959
cat("║ Poder do teste  :", round(poder$power, 4), "\n")
## ║ Poder do teste  : 0.9997
cat("║ Erro Tipo 2     :", round(1 - poder$power, 4), "\n")
## ║ Erro Tipo 2     : 3e-04
cat("╠══════════════════════════════════════════╣\n")
## ╠══════════════════════════════════════════╣
cat("║ DECISÃO:", ifelse(resultado_teste$p.value < 0.05, "REJEITA H0", "NÃO REJEITA H0"), "\n")
## ║ DECISÃO: REJEITA H0
cat("╚══════════════════════════════════════════╝\n")
## ╚══════════════════════════════════════════╝

4. Visualizações — Teste t

Distribuição do Gasto com Saúde por Grupo

ggplot(dados_trabalho, aes(x = gasto_saude_pct_gdp, fill = grupo)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = mean(dados_trabalho$gasto_saude_pct_gdp[dados_trabalho$grupo == "Alta poluição"]),
             color = "red", linetype = "dashed", show.legend = FALSE) +
  geom_vline(xintercept = mean(dados_trabalho$gasto_saude_pct_gdp[dados_trabalho$grupo == "Baixa poluição"]),
             color = "blue", linetype = "dashed", show.legend = FALSE) +
  scale_x_continuous(limits = c(0, 30)) +
  labs(
    title = "Distribuição do Gasto com Saúde por Grupo de Poluição (2018)",
    x     = "Gasto com Saúde (% do PIB)",
    y     = "Densidade",
    fill  = "Grupo de Poluição"
  ) +
  theme_minimal()

Distribuições Separadas Centradas na Média

dados_trabalho <- dados_trabalho |>
  group_by(grupo) |>
  mutate(saude_centrado = gasto_saude_pct_gdp - mean(gasto_saude_pct_gdp)) |>
  ungroup()

ggplot(dados_trabalho, aes(x = saude_centrado, fill = grupo)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = 0, color = "black", linetype = "dashed", show.legend = FALSE) +
  facet_wrap(~ grupo, ncol = 1) +
  labs(
    title = "Distribuição do Gasto com Saúde Centrada na Média por Grupo (2018)",
    x     = "Desvio em relação à média do grupo (% do PIB)",
    y     = "Densidade",
    fill  = "Grupo de Poluição"
  ) +
  theme_minimal()

Distribuição t com Regiões Críticas

x <- seq(-4, 4, length.out = 1000)
y <- dt(x, df = gl)

ggplot(data.frame(x, y), aes(x, y)) +
  geom_line(linewidth = 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$statistic, color = "blue", linetype = "solid", linewidth = 1) +
  annotate("text", x =  t_critico + 0.4, y = 0.35, label = paste("t crítico =",  round( t_critico, 2)), color = "red") +
  annotate("text", x = -t_critico - 0.4, y = 0.35, label = paste("t crítico =", round(-t_critico, 2)), color = "red") +
  annotate("text", x = resultado_teste$statistic, y = 0.38,
           label = paste("t obs =", round(resultado_teste$statistic, 2)), color = "blue") +
  labs(title = "Distribuição t com Regiões Críticas — Poluição vs Saúde (2018)", x = "Valor t", y = "Densidade") +
  theme_minimal()

5. Regressão Linear — Isolando o Beta

Montagem da Base de Regressão

paises_regressao <- inner_join(
  dados_trabalho |> select(codigo),
  pop_densidade  |> select(codigo),
  by = "codigo"
) |>
  inner_join(pop_idosa |> select(codigo), by = "codigo")

base_regressao <- dados_trabalho |>
  filter(codigo %in% paises_regressao$codigo) |>
  inner_join(pop_densidade |> select(codigo, densidade_pop), by = "codigo") |>
  inner_join(pop_idosa     |> select(codigo, pct_65_mais),   by = "codigo")

cat("Países na regressão:", nrow(base_regressao), "\n")
## Países na regressão: 186

Regressão Geral

Modelo: PM2.5 = β₀ + β₁ × Gasto_Saúde + β₂ × Densidade + β₃ × Idosos + β₄ × Continente

reg_geral <- lm(
  pm25 ~ gasto_saude_pct_gdp + densidade_pop + pct_65_mais + continente,
  data = base_regressao
)

summary(reg_geral)
## 
## Call:
## lm(formula = pm25 ~ gasto_saude_pct_gdp + densidade_pop + pct_65_mais + 
##     continente, data = base_regressao)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.675  -8.387  -1.418   6.057  51.031 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.340e+01  2.787e+00  15.573  < 2e-16 ***
## gasto_saude_pct_gdp -5.770e-01  4.147e-01  -1.391 0.165837    
## densidade_pop        4.197e-04  7.415e-04   0.566 0.572059    
## pct_65_mais         -1.138e+00  3.168e-01  -3.592 0.000425 ***
## continenteAmerica   -1.036e+01  3.485e+00  -2.971 0.003374 ** 
## continenteAsia       1.785e+00  2.876e+00   0.621 0.535555    
## continenteEurope    -3.518e+00  5.192e+00  -0.678 0.498952    
## continenteOceania   -2.053e+01  4.334e+00  -4.737 4.42e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.62 on 178 degrees of freedom
## Multiple R-squared:  0.4162, Adjusted R-squared:  0.3932 
## F-statistic: 18.13 on 7 and 178 DF,  p-value: < 2.2e-16

Regressão por Continente

continentes <- unique(base_regressao$continente)

resultados_continente <- lapply(continentes, function(cont) {
  dados_cont <- base_regressao |> filter(continente == cont)
  if (nrow(dados_cont) < 5) return(NULL)

  reg <- lm(pm25 ~ gasto_saude_pct_gdp + densidade_pop + pct_65_mais, data = dados_cont)
  coef_saude <- coef(summary(reg))["gasto_saude_pct_gdp", ]

  data.frame(
    continente    = cont,
    n_paises      = nrow(dados_cont),
    beta1         = round(coef_saude["Estimate"],   4),
    erro_padrao   = round(coef_saude["Std. Error"], 4),
    t_valor       = round(coef_saude["t value"],    4),
    p_valor       = round(coef_saude["Pr(>|t|)"],   4),
    significativo = ifelse(coef_saude["Pr(>|t|)"] < 0.05, "SIM", "NAO")
  )
})

tabela_betas <- do.call(rbind, resultados_continente)
kable(tabela_betas, caption = "Beta do Gasto com Saude por Continente")
Beta do Gasto com Saude por Continente
continente n_paises beta1 erro_padrao t_valor p_valor significativo
Estimate Africa 54 -1.6477 0.9879 -1.6679 0.1016 NAO
Estimate1 America 30 -0.9008 0.4899 -1.8388 0.0774 NAO
Estimate2 Asia 45 0.8379 1.1859 0.7066 0.4838 NAO
Estimate3 Europe 43 -1.3894 0.5840 -2.3789 0.0224 SIM
Estimate4 Oceania 14 -0.2474 0.0883 -2.8032 0.0187 SIM

6. Visualizações — Regressão

Beta por Continente com Intervalo de Confiança

ic_continente <- lapply(continentes, function(cont) {
  dados_cont <- base_regressao |> filter(continente == cont)
  if (nrow(dados_cont) < 5) return(NULL)

  reg  <- lm(pm25 ~ gasto_saude_pct_gdp + densidade_pop + pct_65_mais, data = dados_cont)
  ic   <- confint(reg)["gasto_saude_pct_gdp", ]
  beta <- coef(reg)["gasto_saude_pct_gdp"]

  data.frame(continente = cont, beta = beta, ic_inf = ic[1], ic_sup = ic[2])
})

df_ic <- do.call(rbind, ic_continente)

ic_geral <- confint(reg_geral)["gasto_saude_pct_gdp", ]
df_ic <- rbind(df_ic, data.frame(
  continente = "GERAL",
  beta       = coef(reg_geral)["gasto_saude_pct_gdp"],
  ic_inf     = ic_geral[1],
  ic_sup     = ic_geral[2]
))

ggplot(df_ic, aes(x = reorder(continente, beta), y = beta, color = continente == "GERAL")) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = ic_inf, ymax = ic_sup), width = 0.2, linewidth = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  coord_flip() +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "steelblue"), guide = "none") +
  labs(
    title    = "Beta do Gasto com Saude sobre Poluicao do Ar por Continente (2018)",
    subtitle = "Controlando por densidade populacional e % de idosos",
    x        = "Continente",
    y        = "Beta (Gasto com Saude -> PM2.5)"
  ) +
  theme_minimal()

Scatter Geral com Linha de Regressão

ggplot(base_regressao, aes(x = gasto_saude_pct_gdp, y = pm25, color = continente)) +
  geom_point(alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 1) +
  labs(
    title  = "Poluicao do Ar vs Gasto com Saude (2018)",
    x      = "Gasto com Saude (% do PIB)",
    y      = "PM2.5 (µg/m³)",
    color  = "Continente"
  ) +
  theme_minimal()

7. Conclusão

A regressão geral controlando por densidade populacional, percentual de idosos e continente isola o efeito do gasto com saúde sobre a poluição do ar.

O β₁ representa quanto a poluição (PM2.5) muda a cada 1% a mais de gasto com saúde como percentual do PIB, removendo os efeitos confundidores das outras variáveis.

Fontes: