Análise: Poluição do Ar, Gasto com Saúde e Renda

---

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 == 2018) |>
  drop_na(densidade_pop) |>
  filter(!is.na(codigo) & codigo != "")

PIB per Capita

gdp_total <- read_csv(
  "C:/Users/Guilherme/Downloads/GDP/API_NY.GDP.PCAP.CD_DS2_en_csv_v2_245.csv",
  skip = 4
) |>
  select(-"Indicator Name", -"Indicator Code") |>
  pivot_longer(
    cols      = -c("Country Name", "Country Code"),
    names_to  = "ano",
    values_to = "gdp_per_capita"
  ) |>
  rename(pais = "Country Name", codigo = "Country Code") |>
  mutate(ano = as.integer(ano)) |>
  filter(ano == 2018) |>
  drop_na(gdp_per_capita)

---

2. Montagem da Base Principal

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")

cat("Países na base principal:", nrow(dados_trabalho), "\n")

## Países na base principal: 186

---

3. Teste de Hipótese — Poluição vs Gasto com Saúde

H₀: A média do gasto com saúde é igual entre países de alta e baixa poluição.
H₁: A média do gasto com saúde é diferente entre países de alta e baixa poluição.

# Divide pela mediana de poluição
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

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. Filtro — Segundo Quartil de Renda

Mantemos apenas países com renda per capita entre o percentil 25 e 50, isolando países de renda média-baixa.

q1 <- quantile(gdp_total$gdp_per_capita, 0.25)
q2 <- quantile(gdp_total$gdp_per_capita, 0.50)

gdp_q2 <- gdp_total |>
  filter(gdp_per_capita >= q1 & gdp_per_capita <= q2)

cat("Países no 2° quartil:", nrow(gdp_q2), "\n")

## Países no 2° quartil: 65

cat("Faixa de renda: US$", round(q1), "a US$", round(q2), "\n")

## Faixa de renda: US$ 2163 a US$ 6714

---

6. Montagem da Base de Regressão

paises_regressao <- inner_join(
  dados_trabalho |> select(codigo),
  gdp_q2         |> select(codigo),
  by = "codigo"
) |>
  inner_join(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(gdp_q2        |> select(codigo, gdp_per_capita), by = "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: 52

base_regressao |>
  count(continente) |>
  kable(caption = "Países por Continente na Base de Regressão")

Países por Continente na Base de Regressão

continente	n
Africa	12
America	11
Asia	14
Europe	6
Oceania	9

---

7. Regressão Linear — Impacto da Poluição no Gasto com Saúde

Modelo:

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

• Y = Gasto com saúde (% do PIB)
• β₁ = Efeito da poluição (PM2.5) — coeficiente de interesse
• Controles = Densidade, % idosos, PIB per capita e continente

Regressão Geral

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

summary(reg_inv_geral)

## 
## Call:
## lm(formula = gasto_saude_pct_gdp ~ pm25 + densidade_pop + pct_65_mais + 
##     gdp_per_capita + continente, data = base_regressao)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0038 -1.5835 -0.1004  1.0772  7.1853 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        1.2707063  2.0779563   0.612   0.5441   
## pm25              -0.0226328  0.0416185  -0.544   0.5894   
## densidade_pop      0.0120642  0.0039785   3.032   0.0041 **
## pct_65_mais        0.2049898  0.1796448   1.141   0.2602   
## gdp_per_capita     0.0007270  0.0003295   2.207   0.0327 * 
## continenteAmerica -0.5332345  1.2845522  -0.415   0.6801   
## continenteAsia    -1.6049577  1.1338313  -1.416   0.1641   
## continenteEurope  -1.8202986  2.3420856  -0.777   0.4413   
## continenteOceania  1.4089381  1.4123527   0.998   0.3241   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.589 on 43 degrees of freedom
## Multiple R-squared:  0.3929, Adjusted R-squared:  0.2799 
## F-statistic: 3.478 on 8 and 43 DF,  p-value: 0.003542

Regressão por Continente

continentes <- unique(base_regressao$continente)

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

  reg <- lm(
    gasto_saude_pct_gdp ~ pm25 + densidade_pop + pct_65_mais + gdp_per_capita,
    data = dados_cont
  )

  coef_pm25 <- coef(summary(reg))["pm25", ]

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

tabela_inv <- do.call(rbind, resultados_inv)
kable(tabela_inv, caption = "Beta da Poluição sobre Gasto com Saúde por Continente")

Beta da Poluição sobre Gasto com Saúde por Continente

	continente	n_paises	beta1	erro_padrao	t_valor	p_valor	significativo
Estimate	Africa	12	-0.0087	0.0469	-0.1862	0.8576	NAO
Estimate1	America	11	-0.0005	0.1094	-0.0048	0.9963	NAO
Estimate2	Asia	14	0.0966	0.0749	1.2898	0.2293	NAO
Estimate3	Europe	6	-0.0054	0.0531	-0.1022	0.9351	NAO
Estimate4	Oceania	9	-0.7128	0.6038	-1.1804	0.3032	NAO

---

8. Visualizações — Regressão

β₁ por Continente com Intervalo de Confiança

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

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

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

df_ic_inv <- do.call(rbind, ic_inv)

ic_geral_inv <- confint(reg_inv_geral)["pm25", ]
df_ic_inv <- rbind(df_ic_inv, data.frame(
  continente = "GERAL",
  beta       = coef(reg_inv_geral)["pm25"],
  ic_inf     = ic_geral_inv[1],
  ic_sup     = ic_geral_inv[2]
))

ggplot(df_ic_inv, 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 da Poluicao sobre Gasto com Saude por Continente — 2 Quartil de Renda (2018)",
    subtitle = "Controlando por densidade, % de idosos e renda per capita",
    x        = "Continente",
    y        = "Beta (PM2.5 -> Gasto com Saude % PIB)"
  ) +
  theme_minimal()

Scatter — Poluição vs Gasto com Saúde

ggplot(base_regressao, aes(x = pm25, y = gasto_saude_pct_gdp, 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 — 2 Quartil de Renda (2018)",
    x      = "PM2.5 (µg/m³)",
    y      = "Gasto com Saude (% do PIB)",
    color  = "Continente"
  ) +
  theme_minimal()

---

9. Conclusão

Com p-valor = 0 no teste t, rejeitamos H₀ — há evidência estatística de diferença no gasto com saúde entre países de alta e baixa poluição.

A regressão linear para países do segundo quartil de renda estima o β₁ — o efeito puro da poluição do ar sobre o gasto com saúde, removendo os efeitos da densidade populacional, percentual de idosos, renda per capita e continente.

---

Fontes:

• Poluição do Ar (PM2.5): World Bank
• Gasto com Saúde: World Bank
• População 65+: World Bank
• Densidade Populacional: World Bank
• PIB per Capita: World Bank