Análise de Séries Temporais - Dados PM2.5

Aplicação Interativa em Quarto

Introdução

Esta aplicação apresenta uma análise completa de séries temporais utilizando dados de PM2.5 (material particulado) e dados diários de qualidade do ar. A análise é baseada nos conceitos e métodos apresentados nos scripts originais fornecidos.

Objetivos

• Analisar padrões temporais nos dados de PM2.5
• Identificar tendências e sazonalidade
• Aplicar diferentes modelos de séries temporais
• Comparar valores observados com limites regulamentares (OMS e CONAMA)
• Fornecer visualizações interativas dos resultados

Dados Utilizados

• dados_pm25.csv: Dados mensais de PM2.5 por município (2003-2023)
• dados_dia.csv: Dados diários com valores mínimo, médio e máximo

Exploração dos Dados

Municípios Disponíveis

Primeiro, vamos explorar quais municípios estão disponíveis nos nossos conjuntos de dados:

# Listando municípios disponíveis
municipios <- listar_municipios()

=== CÓDIGOS DE MUNICÍPIOS DISPONÍVEIS ===
Códigos em dados_pm25.csv: 5570 municípios
Primeiros 10: 110001 110002 110003 110004 110005 110155 110006 110007 110008 110009 

Códigos em dados_dia.csv: 2275 municípios
Primeiros 10: 110001 110002 110003 110004 110005 110155 110006 110007 110008 110009 

Códigos em ambos os arquivos: 2275 municípios
Primeiros 10: 110001 110002 110003 110004 110005 110155 110006 110007 110008 110009 

cat("Total de municípios em dados PM2.5:", length(municipios$pm25), "\n")

Total de municípios em dados PM2.5: 5570 

cat("Total de municípios em dados diários:", length(municipios$diarios), "\n")

Total de municípios em dados diários: 2275 

cat("Municípios em ambos os conjuntos:", length(municipios$comuns), "\n\n")

Municípios em ambos os conjuntos: 2275 

cat("Primeiros 20 códigos de municípios disponíveis:\n")

Primeiros 20 códigos de municípios disponíveis:

print(head(municipios$comuns, 20))

 [1] 110001 110002 110003 110004 110005 110155 110006 110007 110008 110009
[11] 110010 110011 110012 110013 110014 110160 171870 171875 171880 171884

Seleção do Município para Análise

Para esta demonstração, vamos analisar o município com código 110001. Você pode modificar este código para analisar outros municípios disponíveis.

# Definindo município para análise
codigo_municipio <- "110001"
cat("Município selecionado:", codigo_municipio, "\n")

Município selecionado: 110001 

Carregamento e Preparação dos Dados

# Carregando dados do município selecionado
dados_pm25 <- carregar_dados_pm25(codigo_municipio = codigo_municipio)

Carregando dados PM2.5...

dados_diarios <- carregar_dados_diarios(codigo_municipio = codigo_municipio)

Carregando dados diários...

# Verificando se os dados foram carregados com sucesso
if (!is.null(dados_pm25)) {
  cat("✓ Dados PM2.5 carregados com sucesso\n")
  cat("  Período:", as.character(min(dados_pm25$dados$data)), "a", 
      as.character(max(dados_pm25$dados$data)), "\n")
  cat("  Total de observações:", nrow(dados_pm25$dados), "\n\n")
} else {
  cat("✗ Erro ao carregar dados PM2.5\n")
}

✓ Dados PM2.5 carregados com sucesso
  Período: 2003-01-01 a 2023-12-01 
  Total de observações: 252 

if (!is.null(dados_diarios)) {
  cat("✓ Dados diários carregados com sucesso\n")
  cat("  Período:", as.character(min(dados_diarios$data)), "a", 
      as.character(max(dados_diarios$data)), "\n")
  cat("  Total de observações:", nrow(dados_diarios), "\n")
} else {
  cat("✗ Erro ao carregar dados diários\n")
}

✓ Dados diários carregados com sucesso
  Período: 2023-01-01 a 2023-10-17 
  Total de observações: 290 

Estatísticas Descritivas

Dados PM2.5 (Mensais)

if (!is.null(dados_pm25)) {
  cat("Estatísticas dos dados PM2.5:\n")
  print(summary(dados_pm25$dados$pm25))
  
  cat("\nInformações adicionais:\n")
  cat("Desvio padrão:", round(sd(dados_pm25$dados$pm25, na.rm = TRUE), 2), "\n")
  cat("Coeficiente de variação:", 
      round(sd(dados_pm25$dados$pm25, na.rm = TRUE) / 
            mean(dados_pm25$dados$pm25, na.rm = TRUE) * 100, 2), "%\n")
}

Estatísticas dos dados PM2.5:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    3.0   629.5   908.5  2029.1  1552.2 32904.0 

Informações adicionais:
Desvio padrão: 3879.65 
Coeficiente de variação: 191.2 %

Dados Diários

if (!is.null(dados_diarios)) {
  cat("Estatísticas dos dados diários:\n\n")
  
  cat("Valores Mínimos:\n")
  print(summary(dados_diarios$min))
  
  cat("\nValores Médios:\n")
  print(summary(dados_diarios$med))
  
  cat("\nValores Máximos:\n")
  print(summary(dados_diarios$max))
}

Estatísticas dos dados diários:

Valores Mínimos:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.040   4.715   6.360   8.108   9.283  34.230 

Valores Médios:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.510   5.678   7.560   9.876  11.275  46.050 

Valores Máximos:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.010   6.875   8.950  12.450  12.898  79.640 

Visualizações e Análise Temporal

Série Temporal PM2.5

Vamos visualizar a série temporal dos dados PM2.5 ao longo do tempo:

if (!is.null(dados_pm25)) {
  # Plotando série temporal PM2.5
  plot(dados_pm25$dados$data, dados_pm25$dados$pm25, 
       type = "l", col = "blue", lwd = 2,
       main = paste("Série Temporal PM2.5 - Município", codigo_municipio),
       xlab = "Data", ylab = "PM2.5 (μg/m³)",
       cex.main = 1.2)
  
  # Adicionando linha de tendência
  dados_temp <- dados_pm25$dados
  dados_temp$tempo_num <- as.numeric(dados_temp$data - min(dados_temp$data))
  modelo_trend <- lm(pm25 ~ tempo_num, data = dados_temp)
  lines(dados_temp$data, fitted(modelo_trend), col = "red", lwd = 2, lty = 2)
  
  legend("topright", 
         legend = c("Dados Observados", "Tendência Linear"),
         col = c("blue", "red"), lty = c(1, 2), lwd = 2)
  
  # Estatísticas da tendência
  cat("\nAnálise de Tendência:\n")
  cat("Coeficiente de tendência:", round(coef(modelo_trend)[2], 6), "μg/m³ por dia\n")
  cat("R² do modelo linear:", round(summary(modelo_trend)$r.squared, 4), "\n")
  
  if (coef(modelo_trend)[2] > 0) {
    cat("Interpretação: Tendência de AUMENTO ao longo do tempo\n")
  } else {
    cat("Interpretação: Tendência de DIMINUIÇÃO ao longo do tempo\n")
  }
}

Análise de Tendência:
Coeficiente de tendência: -0.175504 μg/m³ por dia
R² do modelo linear: 0.0101 
Interpretação: Tendência de DIMINUIÇÃO ao longo do tempo

Decomposição da Série Temporal

A decomposição permite separar os componentes de tendência, sazonalidade e ruído:

if (!is.null(dados_pm25) && length(dados_pm25$serie_ts) >= 24) {
  # Decomposição da série
  decomp <- decompose(dados_pm25$serie_ts, type = "additive")
  
  # Plotando decomposição
  par(mfrow = c(4, 1), mar = c(4, 4, 2, 1))
  
  plot(decomp$x, main = "Série Original", ylab = "PM2.5", col = "blue")
  plot(decomp$trend, main = "Tendência", ylab = "Tendência", col = "red")
  plot(decomp$seasonal, main = "Sazonalidade", ylab = "Sazonal", col = "green")
  plot(decomp$random, main = "Resíduos (Ruído)", ylab = "Resíduos", col = "orange")
  
  par(mfrow = c(1, 1))
  
  # Análise dos componentes
  cat("\nAnálise dos Componentes:\n")
  cat("Variância da série original:", round(var(decomp$x, na.rm = TRUE), 2), "\n")
  cat("Variância da tendência:", round(var(decomp$trend, na.rm = TRUE), 2), "\n")
  cat("Variância sazonal:", round(var(decomp$seasonal, na.rm = TRUE), 2), "\n")
  cat("Variância dos resíduos:", round(var(decomp$random, na.rm = TRUE), 2), "\n")
  
} else {
  cat("Série muito curta para decomposição sazonal (mínimo 24 observações)\n")
}

Análise dos Componentes:
Variância da série original: 15051714 
Variância da tendência: 1115560 
Variância sazonal: 4236544 
Variância dos resíduos: 9918741 

Análise de Autocorrelação

A autocorrelação nos ajuda a entender a dependência temporal dos dados:

if (!is.null(dados_pm25)) {
  par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
  
  # Função de Autocorrelação (ACF)
  acf_result <- acf(dados_pm25$serie_ts, main = "Função de Autocorrelação (ACF)", 
                    lag.max = 24, plot = TRUE)
  
  # Função de Autocorrelação Parcial (PACF)
  pacf_result <- pacf(dados_pm25$serie_ts, main = "Função de Autocorrelação Parcial (PACF)", 
                      lag.max = 24, plot = TRUE)
  
  par(mfrow = c(1, 1))
  
  # Interpretação
  cat("\nInterpretação da Autocorrelação:\n")
  
  # Verificando autocorrelação significativa
  acf_significativo <- sum(abs(acf_result$acf[-1]) > 0.2)
  pacf_significativo <- sum(abs(pacf_result$acf) > 0.2)
  
  cat("Lags com ACF > 0.2:", acf_significativo, "\n")
  cat("Lags com PACF > 0.2:", pacf_significativo, "\n")
  
  if (acf_significativo > 0) {
    cat("A série apresenta dependência temporal significativa\n")
  } else {
    cat("A série apresenta baixa dependência temporal\n")
  }
}

Interpretação da Autocorrelação:
Lags com ACF > 0.2: 4 
Lags com PACF > 0.2: 2 
A série apresenta dependência temporal significativa

Modelagem de Séries Temporais

Modelos Básicos

Vamos comparar diferentes abordagens de modelagem:

Modelo 1: Média Simples

if (!is.null(dados_pm25)) {
  dados <- dados_pm25$dados
  
  # Modelo de média simples
  media_pm25 <- mean(dados$pm25, na.rm = TRUE)
  
  # Plotando
  plot(dados$data, dados$pm25, type = "p", col = "blue", pch = 20,
       main = "Modelo 1: Média Simples", xlab = "Data", ylab = "PM2.5 (μg/m³)")
  abline(h = media_pm25, col = "red", lwd = 3, lty = 2)
  
  legend("topright", 
         legend = c("Dados Observados", paste("Média =", round(media_pm25, 2))),
         col = c("blue", "red"), lty = c(NA, 2), pch = c(20, NA), lwd = c(NA, 3))
  
  # Calculando erro
  erro_media <- sqrt(mean((dados$pm25 - media_pm25)^2, na.rm = TRUE))
  cat("Erro Quadrático Médio (RMSE):", round(erro_media, 2), "\n")
}

Erro Quadrático Médio (RMSE): 3871.95 

Modelo 2: Regressão Linear

if (!is.null(dados_pm25)) {
  dados <- dados_pm25$dados
  dados$tempo_num <- as.numeric(dados$data - min(dados$data))
  
  # Modelo linear
  modelo_linear <- lm(pm25 ~ tempo_num, data = dados)
  
  # Plotando
  plot(dados$data, dados$pm25, type = "p", col = "blue", pch = 20,
       main = "Modelo 2: Regressão Linear", xlab = "Data", ylab = "PM2.5 (μg/m³)")
  lines(dados$data, fitted(modelo_linear), col = "red", lwd = 3)
  
  legend("topright", 
         legend = c("Dados Observados", "Ajuste Linear"),
         col = c("blue", "red"), lty = c(NA, 1), pch = c(20, NA), lwd = c(NA, 3))
  
  # Estatísticas do modelo
  cat("Resultados do Modelo Linear:\n")
  cat("R²:", round(summary(modelo_linear)$r.squared, 4), "\n")
  cat("RMSE:", round(sqrt(mean(residuals(modelo_linear)^2)), 2), "\n")
  cat("Coeficiente de tendência:", round(coef(modelo_linear)[2], 6), "por dia\n")
  
  # Teste de significância da tendência
  p_value <- summary(modelo_linear)$coefficients[2, 4]
  cat("P-valor da tendência:", round(p_value, 4), "\n")
  
  if (p_value < 0.05) {
    cat("A tendência é estatisticamente significativa (p < 0.05)\n")
  } else {
    cat("A tendência NÃO é estatisticamente significativa (p ≥ 0.05)\n")
  }
}

Resultados do Modelo Linear:
R²: 0.0101 
RMSE: 3852.4 
Coeficiente de tendência: -0.175504 por dia
P-valor da tendência: 0.112 
A tendência NÃO é estatisticamente significativa (p ≥ 0.05)

Modelo 3: Harmônico (Sazonal)

if (!is.null(dados_pm25)) {
  dados <- dados_pm25$dados
  n <- nrow(dados)
  
  # Criando variáveis harmônicas
  dados$tempo_num <- seq_len(n)
  dados$sin_12 <- sin(2 * pi * dados$tempo_num / 12)
  dados$cos_12 <- cos(2 * pi * dados$tempo_num / 12)
  dados$sin_6 <- sin(2 * pi * dados$tempo_num / 6)
  dados$cos_6 <- cos(2 * pi * dados$tempo_num / 6)
  
  # Modelo harmônico
  modelo_harm <- lm(pm25 ~ sin_12 + cos_12 + sin_6 + cos_6 + tempo_num, data = dados)
  
  # Plotando
  plot(dados$data, dados$pm25, type = "p", col = "blue", pch = 20,
       main = "Modelo 3: Harmônico (Sazonal)", xlab = "Data", ylab = "PM2.5 (μg/m³)")
  lines(dados$data, fitted(modelo_harm), col = "red", lwd = 2)
  
  legend("topright", 
         legend = c("Dados Observados", "Ajuste Harmônico"),
         col = c("blue", "red"), lty = c(NA, 1), pch = c(20, NA), lwd = c(NA, 2))
  
  # Estatísticas do modelo
  cat("Resultados do Modelo Harmônico:\n")
  cat("R²:", round(summary(modelo_harm)$r.squared, 4), "\n")
  cat("RMSE:", round(sqrt(mean(residuals(modelo_harm)^2)), 2), "\n")
  cat("R² ajustado:", round(summary(modelo_harm)$adj.r.squared, 4), "\n")
  
  # Significância dos componentes sazonais
  coefs <- summary(modelo_harm)$coefficients
  cat("\nSignificância dos componentes sazonais:\n")
  cat("sin_12 p-valor:", round(coefs["sin_12", "Pr(>|t|)"], 4), "\n")
  cat("cos_12 p-valor:", round(coefs["cos_12", "Pr(>|t|)"], 4), "\n")
  cat("sin_6 p-valor:", round(coefs["sin_6", "Pr(>|t|)"], 4), "\n")
  cat("cos_6 p-valor:", round(coefs["cos_6", "Pr(>|t|)"], 4), "\n")
}

Resultados do Modelo Harmônico:
R²: 0.2613 
RMSE: 3327.75 
R² ajustado: 0.2463 

Significância dos componentes sazonais:
sin_12 p-valor: 0 
cos_12 p-valor: 0.0022 
sin_6 p-valor: 3e-04 
cos_6 p-valor: 3e-04 

Comparação de Modelos

if (!is.null(dados_pm25)) {
  # Calculando métricas para comparação
  dados <- dados_pm25$dados
  
  # Modelo média
  media_pm25 <- mean(dados$pm25, na.rm = TRUE)
  rmse_media <- sqrt(mean((dados$pm25 - media_pm25)^2, na.rm = TRUE))
  
  # Modelo linear (já calculado)
  rmse_linear <- sqrt(mean(residuals(modelo_linear)^2))
  r2_linear <- summary(modelo_linear)$r.squared
  
  # Modelo harmônico (já calculado)
  rmse_harmonico <- sqrt(mean(residuals(modelo_harm)^2))
  r2_harmonico <- summary(modelo_harm)$r.squared
  
  # Tabela de comparação
  cat("COMPARAÇÃO DE MODELOS:\n")
  cat("========================\n")
  cat(sprintf("%-15s %8s %8s\n", "Modelo", "RMSE", "R²"))
  cat("------------------------\n")
  cat(sprintf("%-15s %8.2f %8s\n", "Média", rmse_media, "-"))
  cat(sprintf("%-15s %8.2f %8.4f\n", "Linear", rmse_linear, r2_linear))
  cat(sprintf("%-15s %8.2f %8.4f\n", "Harmônico", rmse_harmonico, r2_harmonico))
  
  # Determinando melhor modelo
  rmses <- c(rmse_media, rmse_linear, rmse_harmonico)
  modelos <- c("Média", "Linear", "Harmônico")
  melhor_modelo <- modelos[which.min(rmses)]
  
  cat("\nMelhor modelo (menor RMSE):", melhor_modelo, "\n")
}

COMPARAÇÃO DE MODELOS:
========================
Modelo              RMSE      R²
------------------------
Média           3871.95        -
Linear           3852.40   0.0101
Harmônico       3327.75   0.2613

Melhor modelo (menor RMSE): Harmônico 

Análise dos Dados Diários

Visualização dos Dados Diários

if (!is.null(dados_diarios)) {
  # Plotando dados diários
  plot(dados_diarios$data, dados_diarios$max, type = "p", col = "red", pch = 20,
       main = paste("Dados Diários - Município", codigo_municipio),
       xlab = "Data", ylab = "PM2.5 (μg/m³)",
       ylim = c(min(dados_diarios$min, na.rm = TRUE), 
                max(dados_diarios$max, na.rm = TRUE) * 1.1))
  
  points(dados_diarios$data, dados_diarios$med, col = "orange", pch = 20)
  points(dados_diarios$data, dados_diarios$min, col = "green", pch = 20)
  
  legend("topright", 
         legend = c("Máximo", "Médio", "Mínimo"),
         col = c("red", "orange", "green"),
         pch = 20)
}

Comparação com Limites Regulamentares

Vamos comparar os valores máximos diários com os limites recomendados pela OMS (Organização Mundial da Saúde) e pelo CONAMA (Conselho Nacional do Meio Ambiente).

• Limite OMS (24h): 15 μg/m³
• Limite CONAMA (24h): 50 μg/m³

if (!is.null(dados_diarios)) {
  # Definindo limites
  limite_oms <- 15
  limite_conama <- 50
  
  # Plotando com limites
  plot(dados_diarios$data, dados_diarios$max, type = "p", col = "red", pch = 20,
       main = "Valores Máximos Diários vs Limites Regulamentares",
       xlab = "Data", ylab = "PM2.5 (μg/m³)",
       ylim = c(min(dados_diarios$min, na.rm = TRUE), 
                max(c(dados_diarios$max, limite_conama * 1.1), na.rm = TRUE)))
  
  abline(h = limite_oms, col = "blue", lwd = 2, lty = 2)
  abline(h = limite_conama, col = "purple", lwd = 2, lty = 2)
  
  legend("topright", 
         legend = c("Máximo Diário", "Limite OMS (15)", "Limite CONAMA (50)"),
         col = c("red", "blue", "purple"),
         lty = c(NA, 2, 2), pch = c(20, NA, NA), lwd = c(NA, 2, 2))
  
  # Calculando violações
  violacoes_oms <- sum(dados_diarios$max > limite_oms, na.rm = TRUE)
  violacoes_conama <- sum(dados_diarios$max > limite_conama, na.rm = TRUE)
  total_dias <- nrow(dados_diarios)
  
  cat("\nAnálise de Violações:\n")
  cat("======================\n")
  cat("Total de dias analisados:", total_dias, "\n")
  cat("Dias acima do limite OMS:", violacoes_oms, "(", 
      round(violacoes_oms / total_dias * 100, 2), "%)\n")
  cat("Dias acima do limite CONAMA:", violacoes_conama, "(", 
      round(violacoes_conama / total_dias * 100, 2), "%)\n")
}

Análise de Violações:
======================
Total de dias analisados: 290 
Dias acima do limite OMS: 65 ( 22.41 %)
Dias acima do limite CONAMA: 4 ( 1.38 %)

Conclusão

Esta análise interativa demonstrou como explorar, visualizar e modelar dados de séries temporais de qualidade do ar. Os resultados podem ser utilizados para identificar tendências, padrões sazonais e avaliar a conformidade com os limites regulamentares.

Para análises futuras, recomenda-se:

• Utilizar modelos mais avançados (ARIMA, GAM) se as bibliotecas estiverem disponíveis.
• Incorporar outras variáveis explicativas (meteorologia, dados socioeconômicos).
• Realizar análises espaciais para comparar diferentes municípios.