Probabilidade e Estatística com R
Claudio Rodrigues Nunes
2025-05-09
Introdução
Este documento é o trabalhado final do estudo do Air Quality Data in India (Kaggle), mas especificamente as tabelas city_day.csv, station_day.csv e stations.csv.
Técnicamente eu preferi subir algumas bibliotecas e junções logo no inicio do script, algo do tipo cabeçalho, que neste é caso é o “setup”.
Recriei a página que o professor me informou com as orientações no formato de chek list, me facilitando dentro do processo em marcar os atendimentos.
Check-list do Trabalho de AQI na Índia
1. Leitura e limpeza dos dados
2. Análise exploratória (EDA)
3. Visualização (ggplot2)
Requisitos gerais
3.1 Scatter plots
3.2 Histogramas
3.3 Boxplots
3.4 Barras e colunas
3.5 Gráficos de linha
Hospedei também o código no meu repositório do git: https://github.com/ClaudioRodriguesNunes/probabilidade_estatistica/tree/main
No decorrer das seções eu passo para a fase de tratamento dos dados (tabelas), depois faço uma exploração bruta, tendo como intenção de atender boa parte do check list, e depois finalizo respondendo umas três perguntas.
Peguntas que este trabalho pretende responder:
❓ Pergunta 1: Quais cidades apresentam a maior média de AQI no período completo?
❓ Pergunta 2: Das cidades cuja média inicial de AQI estava acima da média geral, quais conseguiram reduzir para ficar abaixo da média nos dois últimos anos?
❓ Pergunta 3: O comportamento das medições de PM10 e PM2.5 é semelhante em todas as cidades?
Metodologia & Preparação do ambiente de desenvolvimento
Subi as seguintes bibliotecas:
📝 stations <- read_csv(“stations.csv”) %>% clean_names()
📝 station_day<- read_csv(“station_day.csv”)%>% clean_names()
Importamos stations.csv (contendo station_id, city, state) e station_day.csv (contendo medições diárias de AQI e poluentes por estação). Em seguida fizemos um left join para anexar city e state a cada registro diário de AQI.
| Arquivo | Descrição | Principais Variáveis | Período Coberto |
|---|---|---|---|
stations.csv |
Metadados das estações de coleta | station_id, city, state |
— |
station_day.csv |
Medições diárias de AQI e poluentes por estação | station_id, date, aqi,
pm2_5, pm10, no2, … |
2015-01-01 a 2020-12-31 |
city_day.csv |
Agregados diários de AQI por cidade | city, date, aqi,
pm2_5, pm10, … |
2015-01-01 a 2020-12-31 |
Resultados: Exploração de Dados
1) Gráfico de Dispersão (duas variáveis contínuas)
library(ggplot2)
# duas variáveis contínuas: pm2_5 e pm10
aqi_stations %>%
ggplot(aes(x = pm2_5, y = pm10)) +
geom_point(alpha = 0.3) +
labs(
title = "Scatter PM2.5 vs PM10 (todos os registros)",
x = "PM2.5 (µg/m³)",
y = "PM10 (µg/m³)"
)
Neste scatter plot, concluo que PM2.5 e PM10 têm uma relação muito forte e quase linear: à medida que PM2.5 aumenta, PM10 também cresce. A nuvem de pontos concentrada ao longo da reta de tendência reforça o alto grau de correlação entre esses poluentes.
2) Inpecionei a estrutura da tibble com glimpse() e gerei um sumário completo da tibble
glimpse(aqi_stations)## Rows: 108,035
## Columns: 20
## $ station_id <chr> "AP001", "AP001", "AP001", "AP001", "AP001", "AP001", "AP…
## $ date <date> 2017-11-24, 2017-11-25, 2017-11-26, 2017-11-27, 2017-11-…
## $ pm2_5 <dbl> 71.36, 81.40, 78.32, 88.76, 64.18, 72.47, 69.80, 73.96, 8…
## $ pm10 <dbl> 115.75, 124.50, 129.06, 135.32, 104.09, 114.84, 114.86, 1…
## $ no <dbl> 1.75, 1.44, 1.26, 6.60, 2.56, 5.23, 4.69, 4.58, 7.71, 0.9…
## $ no2 <dbl> 20.65, 20.50, 26.00, 30.85, 28.07, 23.20, 20.17, 19.29, 2…
## $ n_ox <dbl> 12.40, 12.08, 14.85, 21.77, 17.01, 16.59, 14.54, 13.97, 1…
## $ nh3 <dbl> 12.19, 10.72, 10.28, 12.91, 11.42, 12.25, 10.95, 10.95, 1…
## $ co <dbl> 0.10, 0.12, 0.14, 0.11, 0.09, 0.16, 0.12, 0.10, 0.10, 0.1…
## $ so2 <dbl> 10.76, 15.24, 26.96, 33.59, 19.00, 10.55, 14.07, 13.90, 1…
## $ o3 <dbl> 109.26, 127.09, 117.44, 111.81, 138.18, 109.74, 118.09, 1…
## $ benzene <dbl> 0.17, 0.20, 0.22, 0.29, 0.17, 0.21, 0.16, 0.17, 0.25, 0.2…
## $ toluene <dbl> 5.92, 6.50, 7.95, 7.63, 5.02, 4.71, 3.52, 2.85, 2.79, 3.8…
## $ xylene <dbl> 0.10, 0.06, 0.08, 0.12, 0.07, 0.08, 0.06, 0.04, 0.07, 0.0…
## $ aqi <dbl> NA, 184, 197, 198, 188, 173, 165, 191, 191, 227, 168, 198…
## $ aqi_bucket <chr> NA, "Moderate", "Moderate", "Moderate", "Moderate", "Mode…
## $ station_name <chr> "Secretariat, Amaravati - APPCB", "Secretariat, Amaravati…
## $ city <chr> "Amaravati", "Amaravati", "Amaravati", "Amaravati", "Amar…
## $ state <chr> "Andhra Pradesh", "Andhra Pradesh", "Andhra Pradesh", "An…
## $ status <chr> "Active", "Active", "Active", "Active", "Active", "Active…dfSummary(aqi_stations, varnumbers = FALSE, style = "grid")## Data Frame Summary
## aqi_stations
## Dimensions: 108035 x 20
## Duplicates: 0
##
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | Variable | Stats / Values | Freqs (% of Valid) | Graph | Valid | Missing |
## +==============+===============================+=======================+=====================+==========+=========+
## | station_id | 1. DL007 | 2009 ( 1.9%) | | 108035 | 0 |
## | [character] | 2. DL008 | 2009 ( 1.9%) | | (100.0%) | (0.0%) |
## | | 3. DL013 | 2009 ( 1.9%) | | | |
## | | 4. DL021 | 2009 ( 1.9%) | | | |
## | | 5. DL033 | 2009 ( 1.9%) | | | |
## | | 6. GJ001 | 2009 ( 1.9%) | | | |
## | | 7. KA003 | 2009 ( 1.9%) | | | |
## | | 8. KA009 | 2009 ( 1.9%) | | | |
## | | 9. MH005 | 2009 ( 1.9%) | | | |
## | | 10. TN001 | 2009 ( 1.9%) | | | |
## | | [ 100 others ] | 87945 (81.4%) | IIIIIIIIIIIIIIII | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | date | min : 2015-01-01 | 2009 distinct values | . : | 108035 | 0 |
## | [Date] | med : 2018-12-02 | | . . : : | (100.0%) | (0.0%) |
## | | max : 2020-07-01 | | . : : : : | | |
## | | range : 5y 6m 0d | | . . : : : : : | | |
## | | | | : : : : : : : : : : | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | pm2_5 | Mean (sd) : 80.3 (76.5) | 22395 distinct values | : | 86410 | 21625 |
## | [numeric] | min < med < max: | | : | (80.0%) | (20.0%) |
## | | 0 < 56 < 1000 | | : | | |
## | | IQR (CV) : 68 (1) | | : | | |
## | | | | : : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | pm10 | Mean (sd) : 158 (123.4) | 29575 distinct values | : | 65329 | 42706 |
## | [numeric] | min < med < max: | | : : | (60.5%) | (39.5%) |
## | | 0 < 122.1 < 1000 | | : : | | |
## | | IQR (CV) : 138.5 (0.8) | | : : : | | |
## | | | | : : : : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | no | Mean (sd) : 23.1 (34.5) | 11963 distinct values | : | 90929 | 17106 |
## | [numeric] | min < med < max: | | : | (84.2%) | (15.8%) |
## | | 0 < 10.3 < 470 | | : | | |
## | | IQR (CV) : 20.1 (1.5) | | : | | |
## | | | | : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | no2 | Mean (sd) : 35.2 (29.5) | 12050 distinct values | : | 91488 | 16547 |
## | [numeric] | min < med < max: | | : | (84.7%) | (15.3%) |
## | | 0 < 27.2 < 448 | | : | | |
## | | IQR (CV) : 31.8 (0.8) | | : . | | |
## | | | | : : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | n_ox | Mean (sd) : 41.2 (45.1) | 15608 distinct values | : | 92535 | 15500 |
## | [numeric] | min < med < max: | | : | (85.7%) | (14.3%) |
## | | 0 < 26.7 < 467.6 | | : | | |
## | | IQR (CV) : 36.5 (1.1) | | : | | |
## | | | | : : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | nh3 | Mean (sd) : 28.7 (24.9) | 9119 distinct values | : | 59930 | 48105 |
## | [numeric] | min < med < max: | | : | (55.5%) | (44.5%) |
## | | 0 < 23.6 < 418.9 | | : | | |
## | | IQR (CV) : 26.2 (0.9) | | : | | |
## | | | | : : | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | co | Mean (sd) : 1.6 (4.4) | 2352 distinct values | : | 95037 | 12998 |
## | [numeric] | min < med < max: | | : | (88.0%) | (12.0%) |
## | | 0 < 0.9 < 175.8 | | : | | |
## | | IQR (CV) : 0.9 (2.7) | | : | | |
## | | | | : | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | so2 | Mean (sd) : 12.3 (13) | 5801 distinct values | : | 82831 | 25204 |
## | [numeric] | min < med < max: | | : | (76.7%) | (23.3%) |
## | | 0 < 8.9 < 195.7 | | : | | |
## | | IQR (CV) : 9.9 (1.1) | | : | | |
## | | | | : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | o3 | Mean (sd) : 38.1 (39.1) | 11166 distinct values | : | 82467 | 25568 |
## | [numeric] | min < med < max: | | : | (76.3%) | (23.7%) |
## | | 0 < 30.8 < 963 | | : | | |
## | | IQR (CV) : 28.2 (1) | | : | | |
## | | | | : | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | benzene | Mean (sd) : 3.4 (11.2) | 3017 distinct values | : | 76580 | 31455 |
## | [numeric] | min < med < max: | | : | (70.9%) | (29.1%) |
## | | 0 < 1.2 < 455 | | : | | |
## | | IQR (CV) : 3.4 (3.3) | | : | | |
## | | | | : | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | toluene | Mean (sd) : 15.3 (29.3) | 8713 distinct values | : | 69333 | 38702 |
## | [numeric] | min < med < max: | | : | (64.2%) | (35.8%) |
## | | 0 < 4.3 < 454.9 | | : | | |
## | | IQR (CV) : 16.8 (1.9) | | : | | |
## | | | | : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | xylene | Mean (sd) : 2.4 (6.5) | 1892 distinct values | : | 22898 | 85137 |
## | [numeric] | min < med < max: | | : | (21.2%) | (78.8%) |
## | | 0 < 0.4 < 170.4 | | : | | |
## | | IQR (CV) : 2.1 (2.7) | | : | | |
## | | | | : | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | aqi | Mean (sd) : 179.7 (131.3) | 930 distinct values | : | 87025 | 21010 |
## | [numeric] | min < med < max: | | : | (80.6%) | (19.4%) |
## | | 8 < 132 < 2049 | | : | | |
## | | IQR (CV) : 168 (0.7) | | : : | | |
## | | | | : : . | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | aqi_bucket | 1. Good | 5510 ( 6.3%) | I | 87025 | 21010 |
## | [character] | 2. Moderate | 29417 (33.8%) | IIIIII | (80.6%) | (19.4%) |
## | | 3. Poor | 11493 (13.2%) | II | | |
## | | 4. Satisfactory | 23636 (27.2%) | IIIII | | |
## | | 5. Severe | 5207 ( 6.0%) | I | | |
## | | 6. Very Poor | 11762 (13.5%) | II | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | station_name | 1. Alandur Bus Depot, Chenna | 2009 ( 1.9%) | | 108035 | 0 |
## | [character] | 2. Bandra, Mumbai - MPCB | 2009 ( 1.9%) | | (100.0%) | (0.0%) |
## | | 3. BWSSB Kadabesanahalli, Be | 2009 ( 1.9%) | | | |
## | | 4. Central School, Lucknow - | 2009 ( 1.9%) | | | |
## | | 5. CRRI Mathura Road, Delhi | 2009 ( 1.9%) | | | |
## | | 6. DTU, Delhi - CPCB | 2009 ( 1.9%) | | | |
## | | 7. IHBAS, Dilshad Garden, De | 2009 ( 1.9%) | | | |
## | | 8. Lalbagh, Lucknow - CPCB | 2009 ( 1.9%) | | | |
## | | 9. Manali, Chennai - CPCB | 2009 ( 1.9%) | | | |
## | | 10. Maninagar, Ahmedabad - GP | 2009 ( 1.9%) | | | |
## | | [ 100 others ] | 87945 (81.4%) | IIIIIIIIIIIIIIII | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | city | 1. Delhi | 45360 (42.0%) | IIIIIIII | 108035 | 0 |
## | [character] | 2. Bengaluru | 11996 (11.1%) | II | (100.0%) | (0.0%) |
## | | 3. Hyderabad | 8752 ( 8.1%) | I | | |
## | | 4. Chennai | 6406 ( 5.9%) | I | | |
## | | 5. Lucknow | 6099 ( 5.6%) | I | | |
## | | 6. Mumbai | 5504 ( 5.1%) | I | | |
## | | 7. Kolkata | 3165 ( 2.9%) | | | |
## | | 8. Jaipur | 3089 ( 2.9%) | | | |
## | | 9. Gurugram | 2831 ( 2.6%) | | | |
## | | 10. Patna | 2678 ( 2.5%) | | | |
## | | [ 16 others ] | 12155 (11.3%) | II | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | state | 1. Delhi | 45360 (42.0%) | IIIIIIII | 108035 | 0 |
## | [character] | 2. Karnataka | 11996 (11.1%) | II | (100.0%) | (0.0%) |
## | | 3. Telangana | 8752 ( 8.1%) | I | | |
## | | 4. Tamil Nadu | 6792 ( 6.3%) | I | | |
## | | 5. Uttar Pradesh | 6099 ( 5.6%) | I | | |
## | | 6. Maharashtra | 5504 ( 5.1%) | I | | |
## | | 7. West Bengal | 3165 ( 2.9%) | | | |
## | | 8. Rajasthan | 3089 ( 2.9%) | | | |
## | | 9. Haryana | 2831 ( 2.6%) | | | |
## | | 10. Bihar | 2678 ( 2.5%) | | | |
## | | [ 11 others ] | 11769 (10.9%) | II | | |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+
## | status | 1. Active | 106102 (98.5%) | IIIIIIIIIIIIIIIIIII | 107711 | 324 |
## | [character] | 2. Inactive | 1609 ( 1.5%) | | (99.7%) | (0.3%) |
## +--------------+-------------------------------+-----------------------+---------------------+----------+---------+Com glimpse() tive uma visão rápida da estrutura e tipos de cada coluna, mas senti dificuldade em identificar frequências e estatísticas de forma imediata. Ao usar dfSummary(), obtenho um sumário completo em HTML que exibe claramente:
Tipos de variáveis,
Quantidade e percentual de NAs,
Medidas descritivas (mínimo, quartis, mediana, média, máximo),
Frequência das categorias em colunas categóricas. Isso me ajudou a planejar o tratamento de valores faltantes e compreender a distribuição geral dos dados.”
3) Gráfico de Barras + Facetas (três variáveis categóricas)
library(dplyr); library(lubridate); library(ggplot2)
# calcula AQI médio por state e season
aqi_sazonal <- aqi_stations %>%
mutate(
date = ymd(date),
month = month(date),
season = case_when(
month %in% 3:5 ~ "Primavera",
month %in% 6:8 ~ "Verão",
month %in% 9:11 ~ "Outono",
TRUE ~ "Inverno"
)
) %>%
group_by(state, season) %>%
summarise(aqi_medio = mean(aqi, na.rm = TRUE), .groups="drop")
# escolhe os 6 estados mais poluentes no geral
top6_states <- media_aqi %>%
group_by(state) %>%
summarise(aqi_estado = mean(media_aqi), .groups="drop") %>%
slice_max(order_by = aqi_estado, n = 6) %>%
pull(state)
# plota só esses
aqi_sazonal %>%
filter(state %in% top6_states) %>%
ggplot(aes(x = season, y = aqi_medio, fill = season)) +
geom_col(show.legend = FALSE) +
facet_wrap(~ state, scales = "free_y") +
labs(
title = "AQI médio sazonal nos 6 estados mais poluentes",
x = "Estação",
y = "AQI Médio"
) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Este gráfico facetado mostra o AQI médio por estação para os seis estados com maior poluição anual. Observa-se consistentemente que o Inverno costuma registrar os picos mais altos de AQI – possivelmente devido à menor dispersão atmosférica – enquanto no Verão os índices caem. Essa sazonalidade reforça a importância de políticas direcionadas aos meses frios
4) Gráfico de Linha (evolução ao longo do tempo)
Como tenho a variável date (contínua no tempo), posso
plotar o AQI médio mensal para todo o país – ou
restringir a algumas cidades/estados.
Série Temporal do AQI Médio Mensal (Requisito 3)
library(dplyr); library(lubridate); library(ggplot2)
aqi_stations %>%
mutate(date = ymd(date)) %>%
group_by(month = floor_date(date, "month")) %>%
summarise(aqi_medio = mean(aqi, na.rm = TRUE), .groups="drop") %>%
ggplot(aes(x = month, y = aqi_medio)) +
geom_line() +
labs(
title = "AQI médio mensal na Índia (todas as estações)",
x = "Mês",
y = "AQI Médio"
)
A série temporal mensal evidencia picos recorrentes de AQI logo no início de cada ano, sugerindo um padrão sazonal possivelmente associado a condições climáticas de dispersão reduzida ou eventos específicos (e.g., queima de resíduos). Há também uma leve tendência de queda após o primeiro trimestre.
5) Vou extrair o nome do mês a partir da data, mas usando stringr em vez de lubridate:
labels_pt <- c(
"Janeiro", "Fevereiro", "Março", "Abril", "Maio", "Junho",
"Julho", "Agosto", "Setembro", "Outubro", "Novembro", "Dezembro"
)
aqi_stations %>%
mutate(
mes_num = month(ymd(date)), # extrai número do mês
mes_texto = labels_pt[mes_num] # mapeia número → nome em PT
) %>%
select(date, mes_num, mes_texto) %>%
head(10)## # A tibble: 10 × 3
## date mes_num mes_texto
## <date> <dbl> <chr>
## 1 2017-11-24 11 Novembro
## 2 2017-11-25 11 Novembro
## 3 2017-11-26 11 Novembro
## 4 2017-11-27 11 Novembro
## 5 2017-11-28 11 Novembro
## 6 2017-11-29 11 Novembro
## 7 2017-11-30 11 Novembro
## 8 2017-12-01 12 Dezembro
## 9 2017-12-02 12 Dezembro
## 10 2017-12-03 12 Dezembroaqi_stations %>%
mutate(
mes_texto = str_to_title(str_extract(month(date), "[:digit:]+"))
) %>%
select(date, mes_texto) %>%
head(10)## # A tibble: 10 × 2
## date mes_texto
## <date> <chr>
## 1 2017-11-24 11
## 2 2017-11-25 11
## 3 2017-11-26 11
## 4 2017-11-27 11
## 5 2017-11-28 11
## 6 2017-11-29 11
## 7 2017-11-30 11
## 8 2017-12-01 12
## 9 2017-12-02 12
## 10 2017-12-03 12O número 11 neste caso representa o mês de novembro.
6) mutate() para converter unidades
aqi_stations %>%
mutate(
pm2_5_mg = pm2_5 / 1000,
pm10_mg = pm10 / 1000
) %>%
select(pm2_5, pm2_5_mg, pm10, pm10_mg) %>%
head(10)## # A tibble: 10 × 4
## pm2_5 pm2_5_mg pm10 pm10_mg
## <dbl> <dbl> <dbl> <dbl>
## 1 71.4 0.0714 116. 0.116
## 2 81.4 0.0814 124. 0.124
## 3 78.3 0.0783 129. 0.129
## 4 88.8 0.0888 135. 0.135
## 5 64.2 0.0642 104. 0.104
## 6 72.5 0.0725 115. 0.115
## 7 69.8 0.0698 115. 0.115
## 8 74.0 0.0740 114. 0.114
## 9 89.9 0.0899 140. 0.140
## 10 87.1 0.0871 131. 0.131Converti PM2.5 e PM10 de µg/m³ para mg/m³ dividindo por 1000. Isso facilita comparações caso eu queira trabalhar em escalas diferentes ou agregar valores em mg.
7) Contagem de valores diferentes
aqi_stations %>%
summarise(
estados = n_distinct(state),
cidades = n_distinct(city),
dias = n_distinct(date)
)## # A tibble: 1 × 3
## estados cidades dias
## <int> <int> <int>
## 1 21 26 2009Verifico que existem 21 estados, 26 cidades e 2009 datas únicas no conjunto de dados, o que confirma a diversidade espacial e temporal da amostra.
8) Contagem de NAs por coluna
aqi_stations %>%
summarise(across(
c(aqi, pm2_5, pm10, no2),
~ sum(is.na(.)),
.names = "n_na_{col}"
))## # A tibble: 1 × 4
## n_na_aqi n_na_pm2_5 n_na_pm10 n_na_no2
## <int> <int> <int> <int>
## 1 21010 21625 42706 16547Identifico que a coluna no2 possui o maior número de NAs, enquanto PM10 tem poucos faltantes. Isso me orienta a focar o tratamento de valores faltantes sobretudo em NO₂ ou a usar imputação quando necessário
9) Perguntas extras e resposta via EDA
Pergunta: “Qual o estado com maior variabilidade de PM2.5 ao longo do ano?”
# calcula desvio padrão de pm2_5 por estado
variabilidade <- aqi_stations %>%
group_by(state) %>%
summarise(sd_pm2_5 = sd(pm2_5, na.rm=TRUE)) %>%
arrange(desc(sd_pm2_5))
variabilidade %>% head(6)## # A tibble: 6 × 2
## state sd_pm2_5
## <chr> <dbl>
## 1 Haryana 90.4
## 2 Bihar 88.4
## 3 Delhi 88.1
## 4 Meghalaya 83.6
## 5 Uttar Pradesh 83.2
## 6 Assam 61.6Descubro que o estado de Haryana apresenta a maior variabilidade de PM2.5 (desvio-padrão), indicando maiores oscilações na poluição ao longo do ano. Esse insight pode guiar investigações sobre causas locais de alta variabilidade.
10) Scatter plots (e geom_smooth)
aqi_stations %>%
ggplot(aes(x = pm2_5, y = pm10)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "Relação PM2.5 × PM10 com reta de regressão",
subtitle = "Todos os registros",
x = "PM2.5 (µg/m³)",
y = "PM10 (µg/m³)"
)
Scatter PM2.5 vs PM10 com reta de regressão
A adição da reta de regressão (geom_smooth(method=‘lm’)) reforça visualmente a magnitude da correlação: o ajuste linear se encaixa bem aos pontos, confirmando que PM2.5 explica grande parte da variação de PM10.
11) Histogramas (facetados por poluente)
aqi_stations %>%
select(aqi, pm2_5, pm10) %>%
pivot_longer(everything(), names_to = "poluente", values_to = "valor") %>%
ggplot(aes(x = valor)) +
geom_histogram(bins = 30, alpha = 0.7) +
facet_wrap(~ poluente, scales = "free") +
labs(
title = "Distribuição de AQI, PM2.5 e PM10",
x = "Valor",
y = "Frequência"
) +
theme_minimal()
Histogramas de AQI, PM2.5 e PM10
Os histogramas mostram distribuição assimétrica de AQI, PM2.5 e PM10. Por exemplo, PM2.5 tem cauda longa à direita, sugerindo poucos dias com valores extremamente altos. Esses detalhes me ajudam a avaliar se transformações logarítmicas seriam adequadas.
12) Boxplots (comparando PM2.5 e PM10)
aqi_stations %>%
select(pm2_5, pm10) %>%
pivot_longer(everything(), names_to = "poluente", values_to = "valor") %>%
ggplot(aes(x = poluente, y = valor)) +
geom_boxplot(outlier.alpha = 0.3) +
labs(
title = "Boxplot de PM2.5 vs PM10",
x = "Poluente",
y = "Concentração (µg/m³)"
)
Boxplots de PM2.5 e PM10
Comparando boxplots, observo que PM10 apresenta outliers mais extremos que PM2.5, confirmando maior variabilidade e possíveis eventos episódicos de poluição grossa. A mediana de PM2.5 é consistentemente menor.
13) Barras e colunas (proporções e stacked)
a) Proporção de dias por estação
aqi_stations %>%
mutate(
month = month(ymd(date)),
season = case_when(
month %in% 3:5 ~ "Primavera",
month %in% 6:8 ~ "Verão",
month %in% 9:11 ~ "Outono",
TRUE ~ "Inverno"
)
) %>%
count(season) %>%
mutate(prop = n / sum(n)) %>%
ggplot(aes(x = season, y = prop, fill = season)) +
geom_col(show.legend = FALSE) +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Proporção de registros por estação do ano",
x = "Estação",
y = "Percentual de registros"
)
Proporção de registros por estação do ano
Vejo que a maior parte dos registros ocorre na primavera e inverno, refletindo o período de maior monitoração ou oscilações sazonais na coleta dos dados.
b) Stacked AQI médio por estação e top-3 estados
top3_states <- media_aqi %>%
group_by(state) %>%
summarise(aqi_estado = mean(media_aqi)) %>%
slice_max(aqi_estado, n = 3) %>%
pull(state)
aqi_stations %>%
filter(state %in% top3_states) %>%
mutate(mes = floor_date(ymd(date), "month")) %>%
group_by(mes, state) %>%
summarise(aqi_medio = mean(aqi, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(x = mes, y = aqi_medio, fill = state)) +
geom_col() +
labs(
title = "AQI médio mensal (stacked) – Top 3 estados",
x = "Mês",
y = "AQI Médio",
fill = "Estado"
)
AQI médio mensal empilhado dos 3 estados mais poluentes
No gráfico empilhado noto que Uttar Pradesh contribui constantemente com a maior fração do AQI mensal, ao passo que Gujarat mostra aumento acentuado entre 2018 e 2020. Isso aponta para diferentes trajetórias de poluição em cada estado.
14) Gráficos de linha (por cidade)
top3_cities <- media_aqi %>%
slice_max(media_aqi, n = 3) %>%
pull(city)
aqi_stations %>%
filter(city %in% top3_cities) %>%
mutate(mes = floor_date(ymd(date), "month")) %>%
group_by(mes, city) %>%
summarise(aqi_medio = mean(aqi, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(x = mes, y = aqi_medio, color = city)) +
geom_line() +
labs(
title = "AQI médio mensal nas 3 cidades mais poluentes",
x = "Mês",
y = "AQI Médio",
color = "Cidade"
)
Série temporal de AQI médio nas 3 cidades mais poluentes
Ao sobrepor as linhas das três cidades mais poluentes, percebo que cada uma tem seu padrão de picos:
“Patna apresenta maior amplitude no começo do ano;” “Delhi mostra picos adicionais no outono;” “Ludhiana tem flutuações mais irregulares.” “Essas diferenças sugerem impactos locais (clima urbano, fontes de emissão) distintos.”
📝❓ Pergunta 1: Cidades com maior média de AQI no período completo
a) Amostra dos dados brutos
head(aqi_stations)## # A tibble: 6 × 20
## station_id date pm2_5 pm10 no no2 n_ox nh3 co so2 o3
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AP001 2017-11-24 71.4 116. 1.75 20.6 12.4 12.2 0.1 10.8 109.
## 2 AP001 2017-11-25 81.4 124. 1.44 20.5 12.1 10.7 0.12 15.2 127.
## 3 AP001 2017-11-26 78.3 129. 1.26 26 14.8 10.3 0.14 27.0 117.
## 4 AP001 2017-11-27 88.8 135. 6.6 30.8 21.8 12.9 0.11 33.6 112.
## 5 AP001 2017-11-28 64.2 104. 2.56 28.1 17.0 11.4 0.09 19 138.
## 6 AP001 2017-11-29 72.5 115. 5.23 23.2 16.6 12.2 0.16 10.6 110.
## # ℹ 9 more variables: benzene <dbl>, toluene <dbl>, xylene <dbl>, aqi <dbl>,
## # aqi_bucket <chr>, station_name <chr>, city <chr>, state <chr>, status <chr>b) Linhas com problemas (AQI faltante)
aqi_stations %>%
filter(is.na(aqi)) %>%
head(10)## # A tibble: 10 × 20
## station_id date pm2_5 pm10 no no2 n_ox nh3 co so2 o3
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AP001 2017-11-24 71.4 116. 1.75 20.6 12.4 12.2 0.1 10.8 109.
## 2 AP001 2018-03-21 48.7 96.5 5.63 14.9 12.6 13.0 0.92 15.4 27.5
## 3 AP001 2018-04-06 29.9 71.8 5.08 9.62 9.25 9.17 0.87 10.3 20.6
## 4 AP001 2018-04-11 13 29 3.8 2.8 4.6 20.1 NA NA 19.6
## 5 AP001 2018-04-12 9.25 35.4 4.7 4.6 6.48 9.06 0.54 6.2 21.7
## 6 AP001 2018-04-18 17.8 54.0 4.51 6.55 7.31 8.37 0.73 10.5 28.5
## 7 AP001 2018-04-20 24.8 67.2 3.41 6.22 6.25 10.4 0.56 15.2 30.3
## 8 AP001 2018-04-23 25 61 3.1 8.8 7.2 7.5 NA NA 19.6
## 9 AP001 2018-04-24 23.4 91.8 4.6 7.1 7.6 7.8 0.62 11.2 35.0
## 10 AP001 2018-04-26 24.8 66.5 4.32 7.23 7.54 9.86 0.6 13.0 29.1
## # ℹ 9 more variables: benzene <dbl>, toluene <dbl>, xylene <dbl>, aqi <dbl>,
## # aqi_bucket <chr>, station_name <chr>, city <chr>, state <chr>, status <chr>c) Tabela tratada (remoção de NAs)
aqi_stations_clean <- aqi_stations %>%
filter(!is.na(aqi))
head(aqi_stations_clean)## # A tibble: 6 × 20
## station_id date pm2_5 pm10 no no2 n_ox nh3 co so2 o3
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AP001 2017-11-25 81.4 124. 1.44 20.5 12.1 10.7 0.12 15.2 127.
## 2 AP001 2017-11-26 78.3 129. 1.26 26 14.8 10.3 0.14 27.0 117.
## 3 AP001 2017-11-27 88.8 135. 6.6 30.8 21.8 12.9 0.11 33.6 112.
## 4 AP001 2017-11-28 64.2 104. 2.56 28.1 17.0 11.4 0.09 19 138.
## 5 AP001 2017-11-29 72.5 115. 5.23 23.2 16.6 12.2 0.16 10.6 110.
## 6 AP001 2017-11-30 69.8 115. 4.69 20.2 14.5 11.0 0.12 14.1 118.
## # ℹ 9 more variables: benzene <dbl>, toluene <dbl>, xylene <dbl>, aqi <dbl>,
## # aqi_bucket <chr>, station_name <chr>, city <chr>, state <chr>, status <chr>d) Cálculo da média de AQI por estado–cidade (dados tratados)
media_aqi <- aqi_stations_clean %>%
group_by(state, city) %>%
summarise(media_aqi = mean(aqi, na.rm = TRUE), .groups = "drop") %>%
arrange(desc(media_aqi))
head(media_aqi, 10)## # A tibble: 10 × 3
## state city media_aqi
## <chr> <chr> <dbl>
## 1 Gujarat Ahmedabad 452.
## 2 Delhi Delhi 239.
## 3 Uttar Pradesh Lucknow 216.
## 4 Bihar Patna 216.
## 5 Haryana Gurugram 211.
## 6 Odisha Talcher 173.
## 7 Jharkhand Jorapokhar 159.
## 8 Odisha Brajrajnagar 150.
## 9 Assam Guwahati 140.
## 10 Rajasthan Jaipur 135.e) Gráfico das Top 10 Cidades/Estados
library(ggplot2)
media_aqi %>%
slice_head(n = 10) %>%
mutate(label = paste(state, city, sep = " – ")) %>%
ggplot(aes(x = reorder(label, media_aqi), y = media_aqi)) +
geom_col(fill = "steelblue") +
coord_flip() +
labs(
title = "Top 10 Cidades/Estados por AQI Médio",
x = NULL,
y = "Média AQI"
)
f) Tabela final da Pergunta 1
media_aqi %>%
slice_head(n = 10) %>%
arrange(desc(media_aqi)) %>%
mutate(
Rank = row_number(),
`Média AQI` = round(media_aqi, 1)
) %>%
select(Rank, Estado = state, Cidade = city, `Média AQI`) %>%
knitr::kable(
col.names = c("Rank","Estado","Cidade","Média AQI"),
align = c("c","l","l","r")
)| Rank | Estado | Cidade | Média AQI |
|---|---|---|---|
| 1 | Gujarat | Ahmedabad | 452.1 |
| 2 | Delhi | Delhi | 238.7 |
| 3 | Uttar Pradesh | Lucknow | 216.5 |
| 4 | Bihar | Patna | 215.6 |
| 5 | Haryana | Gurugram | 211.0 |
| 6 | Odisha | Talcher | 172.9 |
| 7 | Jharkhand | Jorapokhar | 159.3 |
| 8 | Odisha | Brajrajnagar | 150.3 |
| 9 | Assam | Guwahati | 140.1 |
| 10 | Rajasthan | Jaipur | 134.8 |
📝❓ Pergunta 2: Cidades/estados que saíram de acima da média geral para abaixo nos últimos 2 anos
a) Amostra dos dados brutos
head(aqi %>% select(date, state, city, aqi, year))## # A tibble: 6 × 5
## date state city aqi year
## <date> <chr> <chr> <dbl> <dbl>
## 1 2017-11-24 Andhra Pradesh Amaravati NA 2017
## 2 2017-11-25 Andhra Pradesh Amaravati 184 2017
## 3 2017-11-26 Andhra Pradesh Amaravati 197 2017
## 4 2017-11-27 Andhra Pradesh Amaravati 198 2017
## 5 2017-11-28 Andhra Pradesh Amaravati 188 2017
## 6 2017-11-29 Andhra Pradesh Amaravati 173 2017b) Linhas “problemáticas” (sem dados de AQI ou sem comparativo)
aqi %>%
filter(is.na(aqi)) %>%
head(10)## # A tibble: 10 × 21
## station_id date pm2_5 pm10 no no2 n_ox nh3 co so2 o3
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AP001 2017-11-24 71.4 116. 1.75 20.6 12.4 12.2 0.1 10.8 109.
## 2 AP001 2018-03-21 48.7 96.5 5.63 14.9 12.6 13.0 0.92 15.4 27.5
## 3 AP001 2018-04-06 29.9 71.8 5.08 9.62 9.25 9.17 0.87 10.3 20.6
## 4 AP001 2018-04-11 13 29 3.8 2.8 4.6 20.1 NA NA 19.6
## 5 AP001 2018-04-12 9.25 35.4 4.7 4.6 6.48 9.06 0.54 6.2 21.7
## 6 AP001 2018-04-18 17.8 54.0 4.51 6.55 7.31 8.37 0.73 10.5 28.5
## 7 AP001 2018-04-20 24.8 67.2 3.41 6.22 6.25 10.4 0.56 15.2 30.3
## 8 AP001 2018-04-23 25 61 3.1 8.8 7.2 7.5 NA NA 19.6
## 9 AP001 2018-04-24 23.4 91.8 4.6 7.1 7.6 7.8 0.62 11.2 35.0
## 10 AP001 2018-04-26 24.8 66.5 4.32 7.23 7.54 9.86 0.6 13.0 29.1
## # ℹ 10 more variables: benzene <dbl>, toluene <dbl>, xylene <dbl>, aqi <dbl>,
## # aqi_bucket <chr>, station_name <chr>, city <chr>, state <chr>,
## # status <chr>, year <dbl>c) Tabela tratada
aqi_clean2 <- aqi %>%
filter(!is.na(aqi))
head(aqi_clean2)## # A tibble: 6 × 21
## station_id date pm2_5 pm10 no no2 n_ox nh3 co so2 o3
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AP001 2017-11-25 81.4 124. 1.44 20.5 12.1 10.7 0.12 15.2 127.
## 2 AP001 2017-11-26 78.3 129. 1.26 26 14.8 10.3 0.14 27.0 117.
## 3 AP001 2017-11-27 88.8 135. 6.6 30.8 21.8 12.9 0.11 33.6 112.
## 4 AP001 2017-11-28 64.2 104. 2.56 28.1 17.0 11.4 0.09 19 138.
## 5 AP001 2017-11-29 72.5 115. 5.23 23.2 16.6 12.2 0.16 10.6 110.
## 6 AP001 2017-11-30 69.8 115. 4.69 20.2 14.5 11.0 0.12 14.1 118.
## # ℹ 10 more variables: benzene <dbl>, toluene <dbl>, xylene <dbl>, aqi <dbl>,
## # aqi_bucket <chr>, station_name <chr>, city <chr>, state <chr>,
## # status <chr>, year <dbl>d) Cálculo das médias inicial e recente
média geral
media_geral <- mean(aqi_clean2$aqi, na.rm=TRUE)ano inicial e anos finais
ano_ini <- min(aqi_clean2$year)
ano_max <- max(aqi_clean2$year)
med_iniciais <- aqi_clean2 %>%
filter(year == ano_ini) %>%
group_by(state, city) %>%
summarise(media_inicial = mean(aqi, na.rm=TRUE), .groups="drop")
med_recentes <- aqi_clean2 %>%
filter(year >= ano_max-1) %>%
group_by(state, city) %>%
summarise(media_recente = mean(aqi, na.rm=TRUE), .groups="drop")e) Gráfico comparativo (inicial × recente)
library(tidyr)
library(ggplot2)
#### 1. Junta medias iniciais e recentes pelas chaves
comparativo <- med_iniciais %>%
inner_join(med_recentes, by = c("state","city")) %>%
#### 2. Filtra quem melhorou
filter(media_inicial > media_geral, media_recente < media_geral) %>%
#### 3. Prepara para pivot_longer
rename(Inicial = media_inicial, Recente = media_recente) %>%
select(state, city, Inicial, Recente)
#### 4. Empilha para fazer dodge bars
comparativo_long <- comparativo %>%
pivot_longer(
cols = c(Inicial, Recente),
names_to = "Período",
values_to = "AQI"
) %>%
mutate(label = paste(state, city, sep=" – "))
#### 5. Plota
ggplot(comparativo_long, aes(x = reorder(label, AQI), y = AQI, fill = Período)) +
geom_col(position="dodge") +
coord_flip() +
labs(
title = "Comparativo AQI Inicial vs Recente (quem melhorou)",
x = NULL,
y = "AQI",
fill = NULL
)
### f) Tabela final de resposta
med_iniciais %>%
inner_join(med_recentes, by=c("state","city")) %>%
filter(media_inicial > media_geral, media_recente < media_geral) %>%
arrange(desc(media_inicial)) %>%
slice_head(n = 10) %>%
mutate(
Rank = row_number(),
`Média Inicial` = round(media_inicial, 1),
`Média Recente` = round(media_recente, 1)
) %>%
select(Rank, Estado=state, Cidade=city, `Média Inicial`, `Média Recente`) %>%
knitr::kable(align=c("c","l","l","r","r"))| Rank | Estado | Cidade | Média Inicial | Média Recente |
|---|---|---|---|---|
| 1 | Bihar | Patna | 350.6 | 177.8 |
📝❓ Pergunta 3: Comportamento de PM10 vs PM2.5 em cada cidade
a) Amostra dos dados brutos
head(aqi_stations %>% select(pm2_5, pm10, state, city))## # A tibble: 6 × 4
## pm2_5 pm10 state city
## <dbl> <dbl> <chr> <chr>
## 1 71.4 116. Andhra Pradesh Amaravati
## 2 81.4 124. Andhra Pradesh Amaravati
## 3 78.3 129. Andhra Pradesh Amaravati
## 4 88.8 135. Andhra Pradesh Amaravati
## 5 64.2 104. Andhra Pradesh Amaravati
## 6 72.5 115. Andhra Pradesh Amaravatib) Linhas “problemáticas” (menos de 2 pares completos)
aqi_stations %>%
group_by(state, city) %>%
summarise(n_obs = sum(!is.na(pm2_5)&!is.na(pm10))) %>%
filter(n_obs < 2)## # A tibble: 1 × 3
## # Groups: state [1]
## state city n_obs
## <chr> <chr> <int>
## 1 Uttar Pradesh Lucknow 0c) Tabela tratada – cálculo de correlação
correlacoes2 <- aqi_stations %>%
group_by(state, city) %>%
summarise(
n_obs = sum(!is.na(pm2_5)&!is.na(pm10)),
cor_pm = ifelse(n_obs>1, cor(pm2_5, pm10, use="complete.obs"), NA_real_),
.groups="drop"
)
head(correlacoes2)## # A tibble: 6 × 4
## state city n_obs cor_pm
## <chr> <chr> <int> <dbl>
## 1 Andhra Pradesh Amaravati 892 0.923
## 2 Andhra Pradesh Visakhapatnam 1221 0.857
## 3 Assam Guwahati 501 0.784
## 4 Bihar Patna 723 0.655
## 5 Chandigarh Chandigarh 289 0.935
## 6 Delhi Delhi 31766 0.828d) Gráfico scatter (top 6 cidades)
top6 <- correlacoes2 %>%
filter(!is.na(cor_pm)) %>%
slice_max(cor_pm, n=6, with_ties=FALSE) %>%
select(state, city)
aqi_stations %>%
semi_join(top6, by=c("state","city")) %>%
ggplot(aes(x=pm2_5, y=pm10)) +
geom_point(alpha=0.3) +
facet_wrap(~ city, scales="free") +
labs(
title="Scatter PM2.5 vs PM10 nas 6 cidades de maior correlação",
x="PM2.5", y="PM10"
)
e) Tabela final de resposta
correlacoes2 %>%
filter(!is.na(cor_pm)) %>%
slice_max(cor_pm, n=10, with_ties=FALSE) %>%
mutate(
Rank = row_number(),
Correlação = round(cor_pm, 2)
) %>%
select(Rank, Estado=state, Cidade=city, Correlação) %>%
knitr::kable(align=c("c","l","l","r"))| Rank | Estado | Cidade | Correlação |
|---|---|---|---|
| 1 | Meghalaya | Shillong | 1.00 |
| 2 | Mizoram | Aizawl | 0.98 |
| 3 | Kerala | Ernakulam | 0.96 |
| 4 | West Bengal | Kolkata | 0.95 |
| 5 | Chandigarh | Chandigarh | 0.93 |
| 6 | Andhra Pradesh | Amaravati | 0.92 |
| 7 | Kerala | Thiruvananthapuram | 0.90 |
| 8 | Tamil Nadu | Coimbatore | 0.89 |
| 9 | Odisha | Brajrajnagar | 0.89 |
| 10 | Karnataka | Bengaluru | 0.88 |
Conclusão
Neste relatório de análise da qualidade do ar na Índia, respondemos três perguntas-chave:
- Cidades e estados mais poluentes
Top 10 Cidades/Estados Mais Poluentes (Pergunta 1)
| Rank | Estado | Cidade | Média AQI |
|---|---|---|---|
| 1 | Gujarat | Ahmedabad | 452.1 |
| 2 | Delhi | Delhi | 238.7 |
| 3 | Uttar Pradesh | Lucknow | 216.5 |
| 4 | Bihar | Patna | 215.6 |
| 5 | Haryana | Gurugram | 211.0 |
| 6 | Odisha | Talcher | 172.9 |
| 7 | Jharkhand | Jorapokhar | 159.3 |
| 8 | Odisha | Brajrajnagar | 150.3 |
| 9 | Assam | Guwahati | 140.1 |
| 10 | Rajasthan | Jaipur | 134.8 |
- Melhoria nos últimos dois anos
Das cidades/estados que começaram acima da média geral de AQI, as que conseguiram reduzir suas médias anuais para abaixo da média geral nos dois últimos anos foram:
| Estado | Cidade | Média Inicial | Média Recente |
|---|---|---|---|
| Bihar | Patna | 350.6 | 177.8 |
- Similaridade PM2.5 vs PM10
Calculamos o coeficiente de correlação entre PM2.5 e PM10 por cidade e verificamos que a maioria apresenta correlações muito altas (próximas de 1), indicando comportamento fortemente semelhante. Nas 6 cidades de maior correlação, o scatter plot mostra um alinhamento quase linear, confirmando essa relação.
| Estado | Cidade | Correlação |
|---|---|---|
| Meghalaya | Shillong | 1.00 |
| Mizoram | Aizawl | 0.98 |
| Kerala | Ernakulam | 0.96 |
| West Bengal | Kolkata | 0.95 |
| Chandigarh | Chandigarh | 0.93 |
| Andhra Pradesh | Amaravati | 0.92 |
| Kerala | Thiruvananthapuram | 0.90 |
| Tamil Nadu | Coimbatore | 0.89 |
| Odisha | Brajrajnagar | 0.89 |
| Karnataka | Bengaluru | 0.88 |
Limitações e próximos passos:
- Em algumas cidades faltaram dados completos de PM2.5/PM10, o que
excluiu esses casos da análise de correlação.
- Poderíamos estender o estudo agregando por região ou avaliando
sazonalidade mais fina (mensal/diária).
- Uma análise de séries temporais (ARIMA, decomposição sazonal) ajudaria
a entender tendências e picos de poluição.