Relatório de Modelo: SARIMA Baseline

Análise de Desempenho e Comparação com SNaive

1. Introdução

1.1 Objetivo

Este relatório analisa o desempenho do modelo SARIMA (Seasonal AutoRegressive Integrated Moving Average), que serve como uma alternativa de baseline para o sistema de previsão. O SARIMA modela tanto a autocorrelação serial quanto a sazonalidade da série temporal.

1.2 Metodologia

• Modelo: SARIMA
• Estratégia: Rolling Replication (Janela deslizante)
• Horizontes de Previsão: 4, 6 e 8 semanas
• Comparativo: Comparação direta com o modelo Seasonal Naive
• Métricas: RMSE, MAE, MAPE, RMSLE, MASE

2. Dados e Preparação

# 1. Carregar dados brutos
df_raw <- load_raw_data()

#> [2026-01-11 14:46:57] INFO: Carregando dados de: /Users/caiosainvallio/ses-sp/forecast_dengue/data/raw/dengue.RData 
#> [2026-01-11 14:46:57] INFO: Dados carregados: 248865 linhas, 23 colunas

# 2. Agregar por estado
df_state <- aggregate_state(df_raw)

#> [2026-01-11 14:46:57] INFO: Agregando dados por estado... 
#> [2026-01-11 14:46:57] INFO: Dados agregados: 626 semanas

# 3. Preprocessamento (garante coluna 'semana')
df <- preprocess_data(df_state)

#> [2026-01-11 14:46:57] INFO: ========== INICIO: Preprocessamento ========== 
#> [2026-01-11 14:46:57] INFO: Removidas 3 linhas da semana 53 
#> [2026-01-11 14:46:57] WARN: Semanas faltantes detectadas: 2 
#> [2026-01-11 14:46:57] WARN: ATENCAO: Alvo nao sera imputado (sera mantido como NA) 
#> [2026-01-11 14:46:57] INFO: Imputados 8 NAs em mean_temp (fill) 
#> [2026-01-11 14:46:57] INFO: Imputados 8 NAs em mean_max_temp (fill) 
#> [2026-01-11 14:46:57] INFO: Imputados 8 NAs em mean_min_temp (fill) 
#> [2026-01-11 14:46:57] INFO: Imputados 8 NAs em mean_precip (fill) 
#> [2026-01-11 14:46:57] INFO: Dados preprocessados: 625 linhas 
#> [2026-01-11 14:46:57] INFO: ========== FIM: Preprocessamento ==========

# NOTA: SARIMA puro NÃO precisa de features externas
# Modela apenas a própria série temporal (AR + MA + Sazonalidade)
# Por isso, NÃO chamamos make_features() aqui

df_features <- df  # Usar dados preprocessados diretamente

# Verificar coluna semana
if (!"semana" %in% names(df_features)) {
    stop("Erro: Coluna 'semana' não encontrada no dataframe.")
}

# Visualizar estrutura
glimpse(df_features |> select(data_iniSE, est_inc100k, semana, ano))

#> Rows: 625
#> Columns: 4
#> $ data_iniSE  <date> 2014-01-05, 2014-01-12, 2014-01-19, 2014-01-26, 2014-02-0…
#> $ est_inc100k <dbl> 4.688448, 5.657748, 6.589646, 8.467320, 10.189317, 14.3069…
#> $ semana      <dbl> 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
#> $ ano         <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014…

3. Definição do Modelo

O modelo SARIMA estende o ARIMA inluindo termos para a parte sazonal:

\[ SARIMA(p,d,q)(P,D,Q)_s \]

Onde \(s\) é o período sazonal (52 semanas).

# Carregar definicao do modelo
model_name <- "sarima"
model <- get_model(model_name)

cat(sprintf(
    "Modelo: %s\nFamília: %s\nDescrição: %s",
    model$name, model$family, model$description
))

#> Modelo: sarima
#> Família: baseline
#> Descrição: SARIMA - Modelo autoregressivo integrado sazonal

4. Backtesting

Utilizando a mesma configuração de backtesting do relatório anterior para garantir comparabilidade justa.

# Configuracao do Backtest (Idêntica ao SNaive)
backtest_config <- list(
    backtest = list(
        horizons = c(4, 6, 8),
        initial_window = 52 * 5, # 5 anos de treino inicial
        step = 4 # Avancar a cada 4 semanas
    ),
    data = list(
        date_col = "data_iniSE",
        target = "est_inc100k"
    )
)

# Executar Backtest
bt_result <- run_backtest(
    model_name = model_name,
    data = df_features,
    config = backtest_config,
    verbose = FALSE
)

# Resumo
print(bt_result)

#> 
#> === Resultado de Backtest ===
#> 
#> Modelo: sarima 
#> Origens: 90 
#> Previsoes: 270 
#> Duracao: 476.6 s
#> 
#> Metricas por Horizonte:
#>  h      mae     rmse     mape    smape      mase     rmsle  n
#>  4 14.29671 36.32899 25.12562 23.19337 0.4221317 0.2929710 90
#>  6 24.03573 67.46882 37.59521 31.61564 0.7096908 0.4310768 90
#>  8 31.75510 84.90015 52.41051 39.06296 0.9376168 0.5568200 90

5. Avaliação de Desempenho

5.1 Métricas por Horizonte

bt_result$metrics |>
    mutate(
        across(where(is.numeric), \(x) round(x, 4)),
        h_desc = paste(h, "semanas")
    ) |>
    select(Horizonte = h_desc, RMSE = rmse, MAE = mae, MAPE = mape, RMSLE = rmsle, MASE = mase) |>
    kable() |>
    kable_styling(bootstrap_options = c("striped", "hover"))

Horizonte	RMSE	MAE	MAPE	RMSLE	MASE
4 semanas	36.3290	14.2967	25.1256	0.2930	0.4221
6 semanas	67.4688	24.0357	37.5952	0.4311	0.7097
8 semanas	84.9001	31.7551	52.4105	0.5568	0.9376

5.2 Comparação com Baseline (SNaive)

Abaixo comparamos o desempenho do SARIMA com o Seasonal Naive.

# Carregar metricas salvas
metrics_path <- file.path(PROJECT_ROOT, "data", "model_metrics.RData")

if (file.exists(metrics_path)) {
    load(metrics_path)

    # Pegar a ultima execucao do SNaive
    snaive_metrics <- model_metrics |>
        filter(model_name == "seasonal_naive") |>
        arrange(desc(execution_date)) |>
        slice(1)

    # Criar tabela de comparacao para horizonte 4 semanas
    comparison <- data.frame(
        Métrica = c("RMSE (4 sem)", "MAE (4 sem)", "MAPE (4 sem)", "RMSLE (4 sem)"),
        SNaive = c(
            snaive_metrics$horizon_4w_rmse,
            snaive_metrics$horizon_4w_mae,
            snaive_metrics$horizon_4w_mape,
            snaive_metrics$horizon_4w_rmsle
        ),
        SARIMA = c(
            bt_result$metrics[bt_result$metrics$h == 4, "rmse"],
            bt_result$metrics[bt_result$metrics$h == 4, "mae"],
            bt_result$metrics[bt_result$metrics$h == 4, "mape"],
            bt_result$metrics[bt_result$metrics$h == 4, "rmsle"]
        )
    )

    comparison |>
        mutate(
            Diferença = SARIMA - SNaive,
            Melhor = ifelse(SARIMA < SNaive, "SARIMA", "SNaive"),
            `% Melhoria` = ifelse(SNaive > 0, (SNaive - SARIMA) / SNaive * 100, 0)
        ) |>
        mutate(across(where(is.numeric), \(x) round(x, 4))) |>
        kable(caption = "Comparativo Direto: SNaive vs SARIMA (Horizonte 4 semanas)") |>
        kable_styling(bootstrap_options = "striped") |>
        column_spec(5, bold = TRUE, color = ifelse(comparison$SARIMA < comparison$SNaive, "green", "red"))
} else {
    cat("Arquivo de métricas anteriores não encontrado. Execute o relatório 01_seasonal_naive.qmd primeiro.")
}

Comparativo Direto: SNaive vs SARIMA (Horizonte 4 semanas)

Métrica	SNaive	SARIMA	Diferença	Melhor	% Melhoria
RMSE (4 sem)	96.3101	36.3290	-59.9812	SARIMA	62.2792
MAE (4 sem)	42.7400	14.2967	-28.4433	SARIMA	66.5496
MAPE (4 sem)	66.1433	25.1256	-41.0176	SARIMA	62.0133
RMSLE (4 sem)	0.9683	0.2930	-0.6754	SARIMA	69.7449

5.3 Análise Visual (4 semanas)

# Extrair previsoes de multiplos horizontes
preds_multi <- bind_rows(
  extract_predictions(bt_result, horizon = 4) |> mutate(horizon = "4 semanas"),
  extract_predictions(bt_result, horizon = 6) |> mutate(horizon = "6 semanas"),
  extract_predictions(bt_result, horizon = 8) |> mutate(horizon = "8 semanas")
)

# Plot estatico para ser convertido em interativo
p_backtest <- ggplot() +
  geom_line(data = preds_multi, aes(x = target_date, y = actual), color = "black", size = 0.5) +
  geom_line(data = preds_multi, aes(x = target_date, y = predicted, color = horizon), size = 0.5, alpha = 0.8) +
  scale_color_manual(values = c("4 semanas" = "#E74C3C", "6 semanas" = "#3498DB", "8 semanas" = "#2ECC71")) +
  labs(title = "Backtest: Previsões Multi-Horizonte",
       x = "Data", y = "Incidência") +
  theme_minimal()

ggplotly(p_backtest) |>
  layout(legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.2))

5.3 Erros por Regime Epidemiológico

É importante avaliar se o modelo erra mais em picos epidêmicos ou em períodos de baixa transmissão.

# Adicionar classificacao de regime aos resultados
preds_all <- bt_result$predictions
preds_all$regime <- classify_regime(preds_all$actual)

# Calcular metricas por regime
metrics_regime <- compute_metrics_by_regime(preds_all)

metrics_regime |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  datatable(options = list(dom = 't'))

6. Previsão Futura

# Definir atraso de notificação
delay_weeks <- 2
last_reliable_date <- max(df_features$data_iniSE) - (delay_weeks * 7)

# Filtrar dados para treino final
df_final <- df_features |>
    filter(data_iniSE <= last_reliable_date)

cat(sprintf("Ignorando últimas %d semanas (dados provisórios).\n", delay_weeks))

#> Ignorando últimas 2 semanas (dados provisórios).

cat(sprintf("Treinando modelo com dados até: %s\n", last_reliable_date))

#> Treinando modelo com dados até: 2025-12-07

# Treinar com dados confiáveis
model_fit <- model$fit(df_final, list(week_col = "semana"))

# Exibir parametros do modelo selecionado
params <- model$get_params(model_fit)

cat("### Especificação do Modelo Final\n")

#> ### Especificação do Modelo Final

cat(sprintf("- **Ordem ARIMA**: (%s)\n", paste(params$order, collapse = ",")))

#> - **Ordem ARIMA**: (2,1,0,0,1,0,52)

if (!is.null(params$lambda)) {
    cat(sprintf("- **Box-Cox Lambda**: %.4f (0 = Log)\n", params$lambda))
}
cat(sprintf("- **AIC**: %.2f\n", params$aic))

#> - **AIC**: 4522.55

cat(sprintf("- **BIC**: %.2f\n", params$bic))

#> - **BIC**: 4535.58

# Coeficientes
if (length(params$coefficients) > 0) {
    cat("\n**Coeficientes**:\n")
    print(params$coefficients)
}

#> 
#> **Coeficientes**:
#>        ar1        ar2 
#> 0.46294283 0.05842771

# Datas futuras
last_date <- max(df_features$data_iniSE)
future_dates <- seq(from = last_date + 7, by = "week", length.out = 8)

# Previsao
forecast_h8 <- model$predict(model_fit, h = 8, new_data = df_features)
intervals <- model$predict_interval(model_fit, h = 8, level = 0.95, new_data = df_features)

forecast_df <- data.frame(
    Data = future_dates,
    Semana_Epidem = lubridate::epiweek(future_dates),
    Previsao = round(forecast_h8, 2),
    Inferior_95 = round(intervals$lower, 2),
    Superior_95 = round(intervals$upper, 2)
)

forecast_df |>
    kable(caption = "Previsão SARIMA - Próximas 8 Semanas") |>
    kable_styling(full_width = FALSE)

Previsão SARIMA - Próximas 8 Semanas

Data	Semana_Epidem	Previsao	Inferior_95	Superior_95
2025-12-28	53	19.66	-5.58	44.91
2026-01-04	1	20.52	-24.22	65.26
2026-01-11	2	51.62	-11.00	114.24
2026-01-18	3	74.28	-4.43	152.99
2026-01-25	4	88.76	-4.42	181.95
2026-02-01	5	110.26	3.97	216.55
2026-02-08	6	121.75	3.50	240.00
2026-02-15	7	127.57	-1.70	256.84

# Plotar historico recente + previsao
# Historico usado no treino
history_train <- df_final |> 
  tail(104) |> 
  select(Date = data_iniSE, Value = est_inc100k) |>
  mutate(Type = "Histórico", Lower = NA, Upper = NA)

# Historico ignorado
history_ignored <- df_features |>
  filter(data_iniSE > last_reliable_date) |>
  select(Date = data_iniSE, Value = est_inc100k) |>
  mutate(Type = "Provisório (Ignorado)", Lower = NA, Upper = NA)

history <- bind_rows(history_train, history_ignored)

future_viz <- data.frame(
  Date = future_dates,
  Value = forecast_h8,
  Lower = intervals$lower,
  Upper = intervals$upper,
  Type = "Previsão (8 sem)"
)

# Combinar para plot
viz_df <- bind_rows(history, future_viz)

p_future <- ggplot(viz_df, aes(x = Date, y = Value, color = Type)) +
  geom_line() +
  geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = Type), alpha = 0.2, color = NA) +
  scale_color_manual(values = c("Histórico" = "black", "Previsão (8 sem)" = "blue", "Provisório (Ignorado)" = "gray")) +
  scale_fill_manual(values = c("Histórico" = NA, "Previsão (8 sem)" = "blue", "Provisório (Ignorado)" = NA)) +
  labs(title = "Previsão SARIMA - Próximas 8 Semanas",
       y = "Incidência / 100k hab") +
  theme_minimal()

ggplotly(p_future)

7. Salvamento e Registro

7.1 Salvar Modelo

models_dir <- file.path(PROJECT_ROOT, "data", "models")
if (!dir.exists(models_dir)) dir.create(models_dir, recursive = TRUE)

saveRDS(model_fit, file.path(models_dir, "sarima.rds"))
cat("Modelo salvo em data/models/sarima.rds\n")

#> Modelo salvo em data/models/sarima.rds

7.2 Atualizar Leaderboard

metrics_path <- file.path(PROJECT_ROOT, "data", "model_metrics.RData")

new_entry <- data.frame(
    model_name = "sarima",
    execution_date = Sys.time(),
    horizon_4w_rmse = bt_result$metrics[bt_result$metrics$h == 4, "rmse"],
    horizon_4w_mae = bt_result$metrics[bt_result$metrics$h == 4, "mae"],
    horizon_4w_mape = bt_result$metrics[bt_result$metrics$h == 4, "mape"],
    horizon_4w_rmsle = bt_result$metrics[bt_result$metrics$h == 4, "rmsle"],
    horizon_8w_rmse = bt_result$metrics[bt_result$metrics$h == 8, "rmse"],
    horizon_8w_mae = bt_result$metrics[bt_result$metrics$h == 8, "mae"],
    horizon_8w_mape = bt_result$metrics[bt_result$metrics$h == 8, "mape"],
    horizon_8w_rmsle = bt_result$metrics[bt_result$metrics$h == 8, "rmsle"],
    avg_rmse = mean(bt_result$metrics$rmse),
    avg_mape = mean(bt_result$metrics$mape)
)

if (file.exists(metrics_path)) {
    load(metrics_path)
    
    # Remove anterior do mesmo dia se houver (para evitar duplicação em desenvolvimento)
    # Mas aqui vou apenas adicionar
    model_metrics <- rbind(model_metrics, new_entry)
} else {
    model_metrics <- new_entry
}

save(model_metrics, file = metrics_path)

# Mostrar Leaderboard Atual
model_metrics |>
    group_by(model_name) |>
    slice(which.max(execution_date)) |>
    ungroup() |>
    select(Modelo = model_name, RMSE_4w = horizon_4w_rmse, MAPE_4w = horizon_4w_mape, RMSLE_4w = horizon_4w_rmsle) |>
    arrange(RMSE_4w) |>
    kable(caption = "Leaderboard de Modelos (Melhores Métricas Recentes)")

Leaderboard de Modelos (Melhores Métricas Recentes)

Modelo	RMSE_4w	MAPE_4w	RMSLE_4w
sarima	36.32899	25.12562	0.2929710
seasonal_mean	91.28455	65.00128	1.0641488
seasonal_naive	96.31014	66.14327	0.9683374

8. Conclusão

A comparação com o SNaive ajuda a entender se a complexidade do SARIMA traz benefícios reais.

Note

Modelos SARIMA são poderosos para capturar dinâmicas temporais, mas podem ser computacionalmente custosos para ajustar em backtest extensivo e sensíveis a mudanças estruturais na série que violem a estacionariedade assumida.