Análise de Séries Temporais de PM2.5: Modelos Estatísticos e Padrões Sazonais
ASISA 2025 - Diego Ricardo Xavier
2025-07-06
Introdução
Este documento apresenta uma análise abrangente de séries temporais de concentrações de material particulado fino (PM2.5) utilizando diferentes abordagens estatísticas. O PM2.5 é um importante indicador de qualidade do ar, sendo monitorado devido aos seus impactos significativos na saúde pública e no meio ambiente.
A análise inclui decomposição de séries temporais, ajuste de diversos modelos estatísticos (regressão linear, modelos harmônicos, segmentados, ARIMA e GAM), além da comparação com limites estabelecidos pelo CONAMA e pela Organização Mundial da Saúde (OMS).
Configuração do Ambiente
Instalação e Carregamento de Bibliotecas
Primeiro, instalamos e carregamos todas as bibliotecas necessárias para a análise:
Análise Mensal e Padrões Sazonais
Importação e Preparação dos Dados
A análise inicia com a importação dos dados mensais de PM2.5 para o município com código IBGE 130260 (Manaus). Os dados são organizados em uma série temporal para facilitar a análise de padrões temporais.
# Definir caminho do arquivo
caminho <- "C:/Users/diego/Documents/ASISA/ASISA_2025/pm25_anomes_2003_2023/dados_pm25.csv"
# Importar o CSV e selecionar variáveis para Manaus (código 130260)
serie <- read_csv(caminho) %>%
filter(cod6 == "130260") %>%
dplyr::select(ano_mes, pm25)
# Ordenar os dados por data
serie_ord <- serie %>%
arrange(ano_mes)
# Extrair ano e mês da primeira data para configuração da série temporal
ano_inicio <- year(min(serie_ord$ano_mes))
mes_inicio <- month(min(serie_ord$ano_mes))
# Converter para objeto ts (série temporal com frequência mensal)
dados <- ts(serie_ord$pm25,
start = c(ano_inicio, mes_inicio),
frequency = 12)
# Verificar estrutura da série
print(dados)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2003 908 1903 698 1743 1952 1622 1215 648 2243 1107 1733 2488
## 2004 1235 1224 1239 137 618 523 814 2003 1851 1742 1675 1498
## 2005 1106 763 573 1657 832 1017 1033 1308 1787 1816 10273 1929
## 2006 886 727 1746 2755 1179 744 57 1014 2421 137 2994 1814
## 2007 1084 1988 1369 772 1039 1218 1119 1305 2081 2686 858 968
## 2008 961 2064 204 919 1449 1032 939 1486 1173 2802 1641 2858
## 2009 76 613 578 105 981 2089 616 1258 1612 2725 1124 1918
## 2010 731 1475 743 1204 672 701 988 13 1289 183 1271 1366
## 2011 1002 552 1286 2491 836 1213 1065 1294 1526 2213 1558 86
## 2012 96 1452 929 2058 1169 823 893 1046 1693 112 1767 1136
## 2013 772 1175 1408 1445 94 159 1167 688 1295 117 1228 1371
## 2014 1615 844 1446 2054 1682 707 1166 1238 1039 2297 3393 2325
## 2015 1564 2001 953 404 1098 998 1485 902 1615 2882 4219 2528
## 2016 3039 1674 2606 1026 1192 64 806 1439 1261 756 892 2465
## 2017 785 502 715 897 1073 1088 809 1928 1895 1698 1417 1545
## 2018 1162 1532 479 1612 684 989 92 937 8 2146 25 113
## 2019 1862 1805 4927 1112 1032 92 1843 1689 1258 1746 2682 2513
## 2020 1593 1171 783 943 1092 1189 95 1667 111 1164 1705 826
## 2021 815 1254 2017 685 642 814 897 1853 109 2611 1093 947
## 2022 1651 1936 132 1578 847 945 1345 1076 416 4645 1904 1223
## 2023 1122 311 729 475 977 964 812 755 149 1646 2344 2416Decomposição da Série Temporal
A decomposição clássica permite separar os componentes de tendência, sazonalidade e aleatoriedade da série temporal, fornecendo insights sobre os padrões subjacentes nos dados de PM2.5.
# Realizar decomposição clássica
decomp <- decompose(dados)
# Visualizar a decomposição completa
autoplot(decomp) +
ggtitle("Decomposição da Série Temporal (PM2.5)") +
theme_minimal() +
labs(y = "PM2.5 (μg/m³)")
# Extrair componentes individuais
tendencia <- decomp$trend
sazonalidade <- decomp$seasonal
aleatorio <- decomp$random
# Criar plots individuais dos componentes
par(mfrow = c(3, 1), mar = c(4, 4, 2, 1))
plot(tendencia, main = "Componente de Tendência",
ylab = "PM2.5 (μg/m³)", xlab = "Tempo")
plot(sazonalidade, main = "Componente Sazonal",
ylab = "PM2.5 (μg/m³)", xlab = "Tempo")
plot(aleatorio, main = "Componente Aleatório",
ylab = "PM2.5 (μg/m³)", xlab = "Tempo")
Modelos Estatísticos
Para esta seção, utilizaremos dados do município com código IBGE 355030 (São Paulo) para demonstrar diferentes abordagens de modelagem estatística aplicadas a séries temporais de PM2.5.
Preparação dos Dados para Modelagem
# Importar dados para São Paulo (código 355030)
dados <- read_csv(caminho) %>%
filter(cod6 == "355030") %>%
dplyr::select(ano_mes, pm25) %>%
mutate(ano_mes = as.Date(paste0(ano_mes, "-01")))
# Visualizar estrutura dos dados
head(dados)## # A tibble: 6 × 2
## ano_mes pm25
## <date> <dbl>
## 1 2003-01-01 9528
## 2 2003-02-01 5962
## 3 2003-03-01 5002
## 4 2003-04-01 6446
## 5 2003-05-01 15508
## 6 2003-06-01 15451## ano_mes pm25
## Min. :2003-01-01 Min. : 12
## 1st Qu.:2008-03-24 1st Qu.: 1849
## Median :2013-06-16 Median : 3671
## Mean :2013-06-16 Mean : 4545
## 3rd Qu.:2018-09-08 3rd Qu.: 6231
## Max. :2023-12-01 Max. :17304Modelo da Média
O modelo mais simples assume que todas as observações flutuam em torno de uma média constante, sem considerar tendências temporais ou padrões sazonais.
# Calcular média das concentrações de PM2.5
media_pm25 <- mean(dados$pm25, na.rm = TRUE)
# Visualizar dados observados e linha da média
plot(dados$ano_mes, dados$pm25, type = "p", pch = 20,
col = "blue",
xlab = "Tempo",
ylab = "PM2.5 (μg/m³)",
main = "Modelo da Média")
abline(h = media_pm25, col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c("Observados", "Média"),
col = c("blue", "red"),
pch = c(20, NA),
lty = c(NA, 2),
lwd = c(NA, 2),
bty = "n")
## Concentração média de PM2.5: 4545.01 μg/m³Modelo de Regressão Linear
A regressão linear simples modela a relação entre as concentrações de PM2.5 e o tempo, permitindo identificar tendências de longo prazo.
# Ajustar modelo de regressão linear
modelo_lm <- lm(pm25 ~ ano_mes, data = dados)
# Exibir resumo do modelo
summary(modelo_lm)##
## Call:
## lm(formula = pm25 ~ ano_mes, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4889.3 -2778.3 -861.9 1740.0 11947.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9024.61085 1558.70435 5.790 2.11e-08 ***
## ano_mes -0.28223 0.09726 -2.902 0.00404 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3419 on 250 degrees of freedom
## Multiple R-squared: 0.03258, Adjusted R-squared: 0.02871
## F-statistic: 8.42 on 1 and 250 DF, p-value: 0.004042# Calcular valores ajustados
dados$ajustados_lm <- fitted(modelo_lm)
# Visualizar ajuste do modelo
plot(dados$ano_mes, dados$pm25, type = "p", pch = 20,
col = "blue",
xlab = "Tempo",
ylab = "PM2.5 (μg/m³)",
main = "Modelo de Regressão Linear")
lines(dados$ano_mes, dados$ajustados_lm, col = "red", lwd = 2)
legend("topright",
legend = c("Observados", "Ajustados (LM)"),
col = c("blue", "red"),
pch = c(20, NA),
lty = c(NA, 1),
lwd = c(NA, 2),
bty = "n")
Modelo Harmônico
O modelo harmônico utiliza funções trigonométricas para capturar padrões sazonais e cíclicos nas concentrações de PM2.5, considerando diferentes frequências de oscilação.
# Preparar variáveis harmônicas para diferentes frequências
dados <- dados %>%
mutate(
ano_mes_num = seq_len(nrow(dados)),
media = mean(pm25, na.rm = TRUE),
# Componentes anuais (12 meses)
sin_12 = sin(2 * pi * ano_mes_num / 12),
cos_12 = cos(2 * pi * ano_mes_num / 12),
# Componentes semestrais (6 meses)
sin_6 = sin(2 * pi * ano_mes_num / 6),
cos_6 = cos(2 * pi * ano_mes_num / 6),
# Componentes trimestrais (3 meses)
sin_3 = sin(2 * pi * ano_mes_num / 3),
cos_3 = cos(2 * pi * ano_mes_num / 3)
)
# Ajustar modelo harmônico completo
modelo_harm <- lm(pm25 ~ sin_12 + cos_12 +
sin_6 + cos_6 +
sin_3 + cos_3 +
ano_mes + media,
data = dados)
# Exibir resumo do modelo
summary(modelo_harm)##
## Call:
## lm(formula = pm25 ~ sin_12 + cos_12 + sin_6 + cos_6 + sin_3 +
## cos_3 + ano_mes + media, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6388.4 -2158.6 -336.5 1624.4 10478.2
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.096e+03 1.420e+03 6.404 7.70e-10 ***
## sin_12 -1.191e+03 2.774e+02 -4.292 2.55e-05 ***
## cos_12 -1.482e+03 2.773e+02 -5.347 2.06e-07 ***
## sin_6 7.910e+02 2.773e+02 2.853 0.00471 **
## cos_6 2.534e+02 2.773e+02 0.914 0.36168
## sin_3 4.530e+01 2.773e+02 0.163 0.87036
## cos_3 -3.545e+02 2.773e+02 -1.279 0.20221
## ano_mes -2.867e-01 8.863e-02 -3.235 0.00138 **
## media NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3112 on 244 degrees of freedom
## Multiple R-squared: 0.2175, Adjusted R-squared: 0.1951
## F-statistic: 9.69 on 7 and 244 DF, p-value: 1.228e-10# Calcular valores ajustados
dados$ajustados_harm_media <- fitted(modelo_harm)
# Visualizar ajuste do modelo harmônico
plot(dados$ano_mes, dados$pm25, type = "p", pch = 20,
col = "blue",
xlab = "Tempo",
ylab = "PM2.5 (μg/m³)",
main = "Modelo Harmônico")
lines(dados$ano_mes, dados$ajustados_harm_media, col = "red", lwd = 2)
legend("topright",
legend = c("Observados", "Ajustados"),
col = c("blue", "red"),
pch = c(20, NA),
lwd = c(NA, 2),
bty = "n")
Modelo Segmentado
O modelo segmentado (ou regressão por partes) permite identificar pontos de mudança estrutural na série temporal, ajustando diferentes tendências lineares para diferentes períodos.
# Ajustar modelo segmentado com detecção automática de ponto de quebra
modelo_seg <- segmented(lm(pm25 ~ ano_mes, data = dados),
seg.Z = ~ ano_mes)
# Exibir resumo do modelo
summary(modelo_seg)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm(pm25 ~ ano_mes, data = dados), seg.Z = ~ano_mes)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.ano_mes 13422 639.741
##
## Coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29473.706 15828.866 1.862 0.0638 .
## ano_mes -1.849 1.242 -1.488 0.1379
## U1.ano_mes 1.714 1.249 1.372 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3405 on 248 degrees of freedom
## Multiple R-Squared: 0.04815, Adjusted R-squared: 0.03663
##
## Boot restarting based on 6 samples. Last fit:
## Convergence attained in 1 iterations (rel. change 1.5856e-06)# Visualizar ajuste do modelo segmentado
plot(dados$ano_mes, dados$pm25, type = "p", pch = 20,
col = "blue",
xlab = "Tempo",
ylab = "PM2.5 (μg/m³)",
main = "Modelo Segmentado")
lines(dados$ano_mes, fitted(modelo_seg), col = "red", lwd = 2)
# Adicionar linhas verticais nos pontos de inflexão
breakpoints <- modelo_seg$psi[, "Est."]
abline(v = breakpoints, col = "darkgreen", lty = 2, lwd = 2)
legend("topright",
legend = c("Observados", "Ajustados", "Ponto de Inflexão"),
col = c("blue", "red", "darkgreen"),
pch = c(20, NA, NA),
lty = c(NA, 1, 2),
bty = "n")
# Exibir pontos de quebra identificados
cat("Pontos de quebra identificados:", as.Date(breakpoints, origin = "1970-01-01"))## Pontos de quebra identificados: 13422Modelo ARIMA
O modelo ARIMA (AutoRegressive Integrated Moving Average) é especialmente adequado para séries temporais, capturando dependências temporais através de componentes autorregressivos e de médias móveis.
# Ajustar modelo ARIMA com seleção automática de parâmetros
modelo_arima <- auto.arima(dados$pm25)
# Exibir resumo do modelo
summary(modelo_arima)## Series: dados$pm25
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.0927 -0.9741
## s.e. 0.0649 0.0172
##
## sigma^2 = 11930005: log likelihood = -2401.51
## AIC=4809.03 AICc=4809.13 BIC=4819.6
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -390.0805 3433.363 2797.718 -425.1864 450.4232 0.8031944
## ACF1
## Training set -0.02454218# Calcular valores ajustados
valores_ajustados_arima <- fitted(modelo_arima)
# Visualizar ajuste do modelo ARIMA
plot(dados$ano_mes, dados$pm25, type = "p", pch = 20,
col = "blue",
xlab = "Tempo",
ylab = "PM2.5 (μg/m³)",
main = "Modelo ARIMA")
lines(dados$ano_mes, valores_ajustados_arima,
col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c("Observados", "Ajustados (ARIMA)"),
col = c("blue", "red"),
lty = c(NA, 2),
pch = c(20, NA),
bty = "n")
Modelo GAM (Generalized Additive Model)
O modelo GAM utiliza funções suaves (splines) para capturar relações não-lineares complexas entre as variáveis, oferecendo maior flexibilidade na modelagem de tendências.
# Ajustar modelo GAM com função suave do tempo
modelo_gam <- gam(pm25 ~ s(ano_mes_num), data = dados)
# Exibir resumo do modelo
summary(modelo_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## pm25 ~ s(ano_mes_num)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4545.0 213.7 21.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ano_mes_num) 3.622 4.488 2.798 0.0198 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0434 Deviance explained = 5.72%
## GCV = 1.1726e+07 Scale est. = 1.1511e+07 n = 252# Calcular predições
dados$pred_gam <- predict(modelo_gam, newdata = dados)
# Visualizar ajuste do modelo GAM
plot(dados$ano_mes, dados$pm25, type = "p", pch = 20,
col = "blue",
xlab = "Tempo",
ylab = "PM2.5 (μg/m³)",
main = "Modelo GAM")
lines(dados$ano_mes, dados$pred_gam,
col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c("Observados", "Ajustados (GAM)"),
col = c("blue", "red"),
lty = c(NA, 2),
pch = c(20, NA),
bty = "n")
# Diagnóstico de resíduos do modelo GAM
plot(modelo_gam)
# Adicionar linha de referência zero
abline(h = 0, col = "red", lty = 2, lwd = 2)
Comparação e Diagnóstico de Modelos
Critério de Informação de Akaike (AIC)
O AIC permite comparar diferentes modelos, favorecendo aqueles que oferecem melhor ajuste com menor complexidade.
# Comparar modelos usando AIC
aic_comparison <- AIC(modelo_lm, modelo_harm, modelo_seg, modelo_gam)
print(aic_comparison)## df AIC
## modelo_lm 3.000000 4820.182
## modelo_harm 9.000000 4778.717
## modelo_seg 5.000000 4820.094
## modelo_gam 5.622445 4818.940## df AIC
## modelo_harm 9.000000 4778.717
## modelo_gam 5.622445 4818.940
## modelo_seg 5.000000 4820.094
## modelo_lm 3.000000 4820.182Diagnóstico de Resíduos
A análise de resíduos é fundamental para verificar as suposições dos modelos e identificar possíveis problemas no ajuste.
## Diagnóstico do Modelo GAM:
##
## Breusch-Godfrey test for serial correlation of order up to 13
##
## data: Residuals from GCV
## LM test = 43.873, df = 13, p-value = 3.221e-05##
## Diagnóstico do Modelo ARIMA:
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)
## Q* = 31.481, df = 8, p-value = 0.0001153
##
## Model df: 2. Total lags used: 10##
## Diagnóstico do Modelo Harmônico:
##
## Breusch-Godfrey test for serial correlation of order up to 12
##
## data: Residuals
## LM test = 15.127, df = 12, p-value = 0.2346##
## Diagnóstico do Modelo Linear:
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 29.519, df = 10, p-value = 0.001026##
## Diagnóstico do Modelo Segmentado:
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 32.498, df = 10, p-value = 0.0003306Análise de Dados Diários e Conformidade com Padrões
Limites CONAMA e OMS
Esta seção analisa a conformidade das concentrações diárias de PM2.5 com os padrões estabelecidos pelo CONAMA (Conselho Nacional do Meio Ambiente) e pela OMS (Organização Mundial da Saúde).
# Definir caminho para dados diários
caminho_diario <- "C:/Users/diego/Documents/ASISA/ASISA_2025/pm25_anomes_2003_2023/dados_dia.csv"
# Importar dados diários e criar indicadores de conformidade
serie_diaria <- read_csv(caminho_diario) %>%
filter(cod6 == "130260") %>%
mutate(
med_anual = mean(med, na.rm = TRUE),
conama = max > 50, # Limite CONAMA: 50 μg/m³
oms = max > 15 # Limite OMS: 15 μg/m³
)
# Tabela de cruzamento entre limites CONAMA e OMS
cat("Tabela de Cruzamento - Violações dos Limites:\n")## Tabela de Cruzamento - Violações dos Limites:##
## FALSE TRUE
## FALSE 212 75
## TRUE 0 3# Calcular percentuais de violação
total_dias <- nrow(serie_diaria)
violacoes_conama <- sum(serie_diaria$conama, na.rm = TRUE)
violacoes_oms <- sum(serie_diaria$oms, na.rm = TRUE)
cat("\nEstatísticas de Conformidade:\n")##
## Estatísticas de Conformidade:## Total de dias analisados: 290cat("Violações CONAMA (>50 μg/m³):", violacoes_conama,
"(", round(100 * violacoes_conama / total_dias, 1), "%)\n")## Violações CONAMA (>50 μg/m³): 3 ( 1 %)cat("Violações OMS (>15 μg/m³):", violacoes_oms,
"(", round(100 * violacoes_oms / total_dias, 1), "%)\n")## Violações OMS (>15 μg/m³): 78 ( 26.9 %)Função de Visualização por Município
Criamos uma função reutilizável para visualizar dados de PM2.5 por município, incluindo valores máximos, mínimos e médios, além dos limites regulatórios.
# Função para plotar dados de PM2.5 por município
plot_municipio <- function(data, ibge_code) {
# Filtrar dados para o município específico
munic_data <- data %>%
filter(cod6 == ibge_code)
# Criar gráfico com ggplot2
ggplot() +
# Pontos para valores máximos, mínimos e médios
geom_point(data = munic_data, aes(x = ano_mes, y = max),
color = "red", alpha = 0.6, size = 1) +
geom_point(data = munic_data, aes(x = ano_mes, y = min),
color = "green", alpha = 0.6, size = 1) +
geom_point(data = munic_data, aes(x = ano_mes, y = med),
color = "orange", alpha = 0.6, size = 1) +
# Linhas de tendência suavizadas
geom_smooth(data = munic_data, aes(x = ano_mes, y = max),
color = "red", se = FALSE, size = 1.2) +
geom_smooth(data = munic_data, aes(x = ano_mes, y = min),
color = "green", se = FALSE, size = 1.2) +
geom_smooth(data = munic_data, aes(x = ano_mes, y = med),
color = "orange", se = FALSE, size = 1.2) +
# Linhas de referência para limites regulatórios
geom_texthline(yintercept = 15, label = "Limite OMS (15 μg/m³)",
hjust = 0.1, color = "blue", linetype = "dashed") +
geom_texthline(yintercept = 50, label = "Limite CONAMA (50 μg/m³)",
hjust = 0.1, color = "purple", linetype = "dashed") +
# Configurações do tema e labels
labs(
title = paste("Concentrações de PM2.5 - Município", ibge_code),
x = "Tempo",
y = "PM2.5 (μg/m³)",
caption = "Vermelho: Máximo | Verde: Mínimo | Laranja: Médio"
) +
theme_minimal() +
theme(
axis.line = element_line(size = 1),
axis.ticks = element_line(size = 0.8),
axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
plot.title = element_text(size = 14, hjust = 0.5),
legend.background = element_rect(colour = "black"),
legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
legend.position = "bottom"
)
}
Conclusões
Esta análise demonstrou diferentes abordagens para modelagem de séries temporais de PM2.5, cada uma com suas vantagens específicas:
- Modelo da Média: Fornece uma baseline simples, mas não captura variações temporais
- Regressão Linear: Identifica tendências de longo prazo de forma direta
- Modelo Harmônico: Captura padrões sazonais através de componentes trigonométricos
- Modelo Segmentado: Detecta mudanças estruturais e pontos de inflexão
- Modelo ARIMA: Considera dependências temporais e é robusto para previsões
- Modelo GAM: Oferece flexibilidade para capturar relações não-lineares complexas
A comparação através do AIC e análise de resíduos permite selecionar o modelo mais adequado para cada contexto específico. A análise de conformidade com padrões CONAMA e OMS fornece informações cruciais para gestão da qualidade do ar e proteção da saúde pública.