Relatório de Modelo: SARIMAX

SARIMA com Variáveis Exógenas - Análise de Desempenho

1. Introdução

1.1 Objetivo

Este relatório analisa o desempenho do modelo SARIMAX (Seasonal AutoRegressive Integrated Moving Average with eXogenous variables), que estende o SARIMA tradicional ao incorporar variáveis exógenas (regressores externos) como preditoras adicionais.

1.2 Metodologia

• Modelo: SARIMAX (Auto-ajustado)
• Estratégia: Rolling Replication (Janela deslizante) com recálculo de features por fold
• Horizontes de Previsão: 4, 6 e 8 semanas
• Comparativo: Comparação com SNaive, SARIMA e modelos anteriores
• Métricas: RMSE, MAE, MAPE, RMSLE, MASE
• Paralelismo: Execução paralela via furrr para otimização de tempo

1.3 Variáveis Exógenas: Conceitos Estatísticos

O modelo SARIMAX incorpora variáveis exógenas (\(\mathbf{X}_t\)) à estrutura tradicional do SARIMA:

\[Y_t = \boldsymbol{\beta}' \mathbf{X}_t + \eta_t\]

onde \(\eta_t\) segue um processo SARIMA:

\[\phi(B)\Phi(B^s)(1-B)^d(1-B^s)^D \eta_t = \theta(B)\Theta(B^s)\varepsilon_t\]

1.3.1 Tipos de Variáveis Exógenas

Categorias de Regressores

1. Climáticas Defasadas: Temperatura e precipitação com lags (4-8 semanas) para respeitar o ciclo biológico do vetor.
2. Tendência Climática: Médias móveis de temperatura para capturar persistência térmica.

Nota: Optamos por não incluir harmônicos de Fourier ou lags autoregressivos como regressores externos (xreg) nesta iteração, pois o componente SARIMA (Seasonal AR+MA) já modela a sazonalidade estocástica e a autocorrelação de forma eficiente. A inclusão duplicada poderia gerar colinearidade e instabilidade.

1.3.2 Importância do Lag Epidemiológico

O ciclo biológico do Aedes aegypti impõe um delay natural entre condições climáticas favoráveis e aumento de casos:

• Temperatura: Afeta desenvolvimento larval (7-14 dias) + período de incubação extrínseco (8-12 dias)
• Precipitação: Influência na disponibilidade de criadouros com atraso de 2-4 semanas

Por isso, utilizamos lags de 4-8 semanas para variáveis climáticas.

1.3.3 Prevenção de Data Leakage

Cuidado Crítico

Ao incluir variáveis derivadas do alvo (como lag_1, roll_mean_4), é essencial que estas sejam calculadas apenas com dados disponíveis até o momento da previsão. No backtesting, isso significa recalcular features a cada fold.

Para garantir isso, utilizamos a função run_backtest_parallel() que:

1. Recebe dados base (sem features calculadas)
2. Para cada origem de backtest, chama make_features_safe(data, config, cutoff_idx = origin)
3. Treina o modelo apenas com features calculadas até o ponto de corte

2. Dados e Preparação

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

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

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

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

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

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

# IMPORTANTE: Para SARIMAX, NÃO chamamos make_features() aqui!
# O backtest paralelo vai calcular features por fold para evitar leakage

# Dados BASE (sem features) para o backtesting
df_base <- df |>
    mutate(data_iniSE = as.Date(data_iniSE)) |>
    arrange(data_iniSE)

# Carregar config de features para passar ao backtester
config_features <- tryCatch(load_config("features"), error = function(e) list())
if (length(config_features) == 0) {
    config_features <- list(temporal = list(enabled = TRUE), seasonal = list(enabled = TRUE))
}

cat("Dados base carregados:", nrow(df_base), "linhas\n")

#> Dados base carregados: 625 linhas

cat("Features serão calculadas por fold no backtesting (sem leakage)\n")

#> Features serão calculadas por fold no backtesting (sem leakage)

# Visualizar estrutura
glimpse(df_base |> 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 SARIMAX estende o SARIMA ao incluir regressores externos:

\[ SARIMAX(p,d,q)(P,D,Q)_s + \boldsymbol{\beta}' \mathbf{X}_t \]

Os componentes são:

• \((p,d,q)\): Ordem AR, diferenciação, MA da parte não-sazonal
• \((P,D,Q)_s\): Ordem sazonal com período \(s = 52\) (semanal)
• \(\mathbf{X}_t\): Vetor de variáveis exógenas

# Carregar definição do modelo
model_name <- "sarimax"
model <- get_model(model_name)

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

#> Modelo: sarimax
#> Família: statistical
#> Descrição: SARIMAX - Seasonal ARIMA com regressores externos dinamicos

# Lista de regressores selecionados
cat("\n### Estratégia de Regressores:\n")

#> 
#> ### Estratégia de Regressores:

cat("Apenas variáveis climáticas foram selecionadas para complementar a estrutura SARIMA:\n")

#> Apenas variáveis climáticas foram selecionadas para complementar a estrutura SARIMA:

cat("- temperature_lag4: Temperatura média (lag 4 semanas)\n")

#> - temperature_lag4: Temperatura média (lag 4 semanas)

cat("- precipitation_lag4: Precipitação (lag 4 semanas)\n")

#> - precipitation_lag4: Precipitação (lag 4 semanas)

cat("- temp_roll_mean_4: Tendência de temperatura (média móvel 4 semanas)\n")

#> - temp_roll_mean_4: Tendência de temperatura (média móvel 4 semanas)

cat("\n*Foram excluídos termos de Fourier e lags da variável alvo para evitar redundância com os componentes Sazonal e AR do modelo SARIMA.*\n")

#> 
#> *Foram excluídos termos de Fourier e lags da variável alvo para evitar redundância com os componentes Sazonal e AR do modelo SARIMA.*

4. Backtesting

4.1 Prevenção de Data Leakage

Modo Seguro Ativado

O backtesting utiliza run_backtest_parallel() que recalcula todas as features a cada fold usando apenas dados disponíveis até aquele ponto. Isso garante que:

1. Lags são calculados apenas com observações passadas
2. Médias móveis não “vazam” informação do futuro
3. A avaliação reflete o desempenho real em produção

4.2 Execução com Paralelismo

Utilizamos a execução paralela para reduzir o tempo de processamento:

# Seleção Racional de Features (Baseada na Decomposição)
# A análise de decomposição mostrou forte sazonalidade e componente AR.
# O modelo SARIMA puro já captura bem a sazonalidade estocástica e a autocorrelação.
# Portanto, para o SARIMAX, devemos evitar redundância:
# - Remover 'lags' da variável alvo (deixar o componente AR do ARIMA cuidar disso)
# - Remover termos de Fourier (deixar o componente Sazonal do SARIMA cuidar disso)
# - Manter apenas variáveis CLIMÁTICAS que trazem informação nova

selected_xregs <- c(
    "temperature_lag4", # Efeito biológico da temperatura
    "precipitation_lag4", # Efeito de criadouros
    "temp_roll_mean_4" # Tendência térmica recente
)

# Configuração do Backtest
backtest_config <- list(
    backtest = list(
        horizons = c(4, 6, 8),
        initial_window = 52 * 5, # 5 anos de treino inicial
        step = 4 # Avançar a cada 4 semanas
    ),
    data = list(
        date_col = "data_iniSE",
        target = "est_inc100k"
    ),
    # Config de features para recálculo por fold
    features = config_features
)

# Configuração do modelo SARIMAX refinado
model_config <- list(
    target = "est_inc100k",
    season = 52,
    auto = TRUE,
    lambda = 0, # Box-Cox (Log)
    stepwise = TRUE,
    approximation = TRUE,
    xreg_cols = selected_xregs, # Forçar uso apenas das climáticas
    seasonal = list(D = 1) # Forçar diferenciação sazonal para evitar previsão "flat"
)

cat("### Estratégia de Modelagem Refinada\n")

#> ### Estratégia de Modelagem Refinada

cat("Selecionando apenas regressores climáticos para evitar colinearidade com estrutura SARIMA:\n")

#> Selecionando apenas regressores climáticos para evitar colinearidade com estrutura SARIMA:

print(selected_xregs)

#> [1] "temperature_lag4"   "precipitation_lag4" "temp_roll_mean_4"

# Verificar disponibilidade de furrr
has_parallel <- requireNamespace("furrr", quietly = TRUE) &&
    requireNamespace("future", quietly = TRUE)

if (has_parallel) {
    cat("\n✅ Pacotes furrr/future disponíveis - usando execução PARALELA\n")

    # Executar backtest paralelo
    bt_result <- run_backtest_parallel(
        model_name = model_name,
        data_base = df_base,
        config = backtest_config,
        model_config = model_config,
        verbose = TRUE
    )

    cat(sprintf("\n🚀 Backtest executado com %d workers\n", bt_result$n_workers))
} else {
    cat("\n⚠️ Pacotes furrr/future não disponíveis - usando execução SEQUENCIAL\n")

    # Fallback para versão sequencial segura
    bt_result <- run_backtest_safe(
        model_name = model_name,
        data_base = df_base,
        config = backtest_config,
        model_config = model_config,
        verbose = TRUE
    )
}

#> 
#> ✅ Pacotes furrr/future disponíveis - usando execução PARALELA
#> [2026-01-11 14:57:29] INFO: ========== INICIO: Backtest PARALLEL: sarimax ========== 
#> [2026-01-11 14:57:29] INFO: Origens: 90 | Horizontes: 4, 6, 8 | Workers: 15 
#> [2026-01-11 14:57:29] INFO: Modo: Features recalculadas por fold (sem leakage)

#> [2026-01-11 15:17:46] INFO: Previsoes geradas: 270 
#> [2026-01-11 15:17:46] INFO: Erros: 0 
#> [2026-01-11 15:17:46] INFO: Tempo total: 1216.5 segundos 
#> [2026-01-11 15:17:46] INFO: ========== FIM: Backtest PARALLEL: sarimax (1216.53s) ========== 
#> 
#> 🚀 Backtest executado com 15 workers

# Verificar modo seguro
if (isTRUE(bt_result$safe_mode)) {
    cat("✅ Backtest executado em modo SAFE (sem data leakage)\n")
}

#> ✅ Backtest executado em modo SAFE (sem data leakage)

cat(sprintf("⏱️ Tempo total: %.1f segundos\n", bt_result$duration))

#> ⏱️ Tempo total: 1216.5 segundos

# Resumo
print(bt_result)

#> 
#> === Resultado de Backtest ===
#> 
#> Modelo: sarimax 
#> Origens: 90 
#> Previsoes: 270 
#> Duracao: 1216.5 s
#> 
#> Metricas por Horizonte:
#>  h      mae     rmse     mape    smape     mase     rmsle  n
#>  4 34.78484 90.69557 43.82694 47.52683 1.027074 0.7821497 90
#>  6 38.68131 92.93367 54.77168 58.48330 1.142124 1.2348616 90
#>  8 41.54976 98.37409 65.16628 66.22393 1.226819 1.2545504 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(caption = "Métricas de Performance do SARIMAX por Horizonte") |>
    kable_styling(bootstrap_options = c("striped", "hover"))

Métricas de Performance do SARIMAX por Horizonte

Horizonte	RMSE	MAE	MAPE	RMSLE	MASE
4 semanas	90.6956	34.7848	43.8269	0.7821	1.0271
6 semanas	92.9337	38.6813	54.7717	1.2349	1.1421
8 semanas	98.3741	41.5498	65.1663	1.2546	1.2268

5.2 Comparação com Baselines

Comparação do SARIMAX com modelos anteriores (SNaive e SARIMA).

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

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

    # Pegar as últimas execuções de cada modelo
    snaive_metrics <- model_metrics |>
        filter(model_name == "seasonal_naive") |>
        arrange(desc(execution_date)) |>
        slice(1)

    sarima_metrics <- model_metrics |>
        filter(model_name == "sarima") |>
        arrange(desc(execution_date)) |>
        slice(1)

    if (nrow(snaive_metrics) > 0 && nrow(sarima_metrics) > 0) {
        # Criar tabela de comparação 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(
                sarima_metrics$horizon_4w_rmse,
                sarima_metrics$horizon_4w_mae,
                sarima_metrics$horizon_4w_mape,
                sarima_metrics$horizon_4w_rmsle
            ),
            SARIMAX = 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(
                SNaive = round(as.numeric(SNaive), 4),
                SARIMA = round(as.numeric(SARIMA), 4),
                SARIMAX = round(as.numeric(SARIMAX), 4),
                `Melhor` = case_when(
                    SARIMAX <= SARIMA & SARIMAX <= SNaive ~ "SARIMAX",
                    SARIMA <= SARIMAX & SARIMA <= SNaive ~ "SARIMA",
                    TRUE ~ "SNaive"
                )
            ) |>
            kable(caption = "Comparativo: SNaive vs SARIMA vs SARIMAX (Horizonte 4 semanas)") |>
            kable_styling(bootstrap_options = c("striped", "hover")) |>
            column_spec(5, bold = TRUE)
    } else {
        cat("Métricas completas de modelos anteriores não disponíveis.")
    }
} else {
    cat("Arquivo de métricas anteriores não encontrado. Execute os relatórios 01_seasonal_naive.qmd e 03_sarima.qmd primeiro.")
}

Comparativo: SNaive vs SARIMA vs SARIMAX (Horizonte 4 semanas)

Métrica	SNaive	SARIMA	SARIMAX	Melhor
RMSE (4 sem)	96.3101	36.3290	90.6956	SARIMA
MAE (4 sem)	42.7400	14.2967	34.7848	SARIMA
MAPE (4 sem)	66.1433	25.1256	43.8269	SARIMA
RMSLE (4 sem)	0.9683	0.2930	0.7821	SARIMA

5.3 Análise Visual Multi-Horizonte

# Extrair previsões de múltiplos 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 estático para ser convertido em interativo
p_backtest <- ggplot() +
    geom_line(data = preds_multi, aes(x = target_date, y = actual), color = "black", linewidth = 0.5) +
    geom_line(data = preds_multi, aes(x = target_date, y = predicted, color = horizon), linewidth = 0.5, alpha = 0.8) +
    scale_color_manual(values = c("4 semanas" = "#E74C3C", "6 semanas" = "#3498DB", "8 semanas" = "#2ECC71")) +
    labs(
        title = "Backtest SARIMAX: Previsões Multi-Horizonte",
        subtitle = "Features recalculadas por fold (sem data leakage)",
        x = "Data", y = "Incidência / 100k hab"
    ) +
    theme_minimal()

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

5.4 Erros por Regime Epidemiológico

Avaliação do desempenho em diferentes regimes de transmissão (baixo, médio, alto).

# Adicionar classificação de regime aos resultados
preds_all <- bt_result$predictions
preds_all$regime <- classify_regime(preds_all$actual)

# Calcular métricas 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"),
        caption = "Métricas por Regime Epidemiológico"
    )

# Visualizar erros por regime
preds_all |>
    mutate(error = predicted - actual) |>
    ggplot(aes(x = regime, y = error, fill = regime)) +
    geom_boxplot(alpha = 0.7) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
    scale_fill_manual(values = c("low" = "#2ECC71", "medium" = "#F39C12", "high" = "#E74C3C")) +
    labs(
        title = "Distribuição de Erros por Regime Epidemiológico",
        x = "Regime", y = "Erro (Previsto - Real)"
    ) +
    theme_minimal() +
    theme(legend.position = "none")

6. Previsão Futura

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

# Filtrar dados base para treino final
df_final_base <- df_base |>
    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

# Para DEPLOY: calcular features com dados filtrados
# (diferente do backtesting que recalcula por fold)
df_model <- make_features_for_prediction(df_final_base, config_features)

#> [2026-01-11 15:17:46] INFO: ========== INICIO: Feature Engineering ========== 
#> [2026-01-11 15:17:46] INFO: Adicionando features temporais... 
#> [2026-01-11 15:17:46] INFO: Adicionando features de sazonalidade... 
#> [2026-01-11 15:17:46] INFO: Adicionando features climaticas... 
#> [2026-01-11 15:17:46] INFO: Features criadas: 47 novas colunas 
#> [2026-01-11 15:17:46] INFO: ========== FIM: Feature Engineering ==========

# Treinar modelo final com todos os dados
model_fit <- model$fit(df_model, model_config)

# Exibir parâmetros 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[1:3], collapse = ",")))

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

cat(sprintf("- **Ordem Sazonal**: (%s)\n", paste(params$order[4:6], collapse = ",")))

#> - **Ordem Sazonal**: (1,1,0)

cat(sprintf("- **Período**: %s\n", params$order[7]))

#> - **Período**: 52

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

#> - **Box-Cox Lambda**: 0 (0 = Log)

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

#> - **AIC**: -592.50

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

#> - **BIC**: -566.49

# Regressores utilizados
if (!is.null(params$feature_names) && length(params$feature_names) > 0) {
    cat("\n**Regressores Exógenos Utilizados**:\n")
    for (feat in params$feature_names) {
        cat(sprintf("- %s\n", feat))
    }
}

#> 
#> **Regressores Exógenos Utilizados**:
#> - temperature_lag4
#> - precipitation_lag4
#> - temp_roll_mean_4

# Coeficientes
if (length(params$coefficients) > 0) {
    cat("\n**Coeficientes Estimados**:\n")
    coef_df <- data.frame(
        Parâmetro = names(params$coefficients),
        Estimativa = round(as.numeric(params$coefficients), 4)
    )
    print(kable(coef_df))
}

#> 
#> **Coeficientes Estimados**:
#> 
#> 
#> |Parâmetro          | Estimativa|
#> |:------------------|----------:|
#> |ar1                |     0.3280|
#> |sar1               |    -0.4490|
#> |temperature_lag4   |     0.0026|
#> |precipitation_lag4 |     0.0025|
#> |temp_roll_mean_4   |    -0.0386|

# Datas para previsão (seguindo padrão SARIMA: iniciando após o último dado disponível, T+1)
# O modelo foi treinado até T-2, mas projetamos 8 semanas à frente para cobrir 
# o período de nowcast (T-1, T) e futuro (T+1..T+6), rotulando como T+1..T+8
# Usando df_base para garantir data máxima correta
last_date <- max(df_base$data_iniSE)
future_dates <- seq(from = last_date + 7, by = "week", length.out = 8)

# Previsão
forecast_h8 <- model$predict(model_fit, h = 8, new_data = df_model)
intervals <- model$predict_interval(model_fit, h = 8, level = 0.95, new_data = df_model)

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 SARIMAX - Próximas 8 Semanas") |>
    kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Previsão SARIMAX - Próximas 8 Semanas

Data	Semana_Epidem	Previsao	Inferior_95	Superior_95
2025-12-28	53	20.43	15.50	26.91
2026-01-04	1	21.21	12.49	31.25
2026-01-11	2	22.11	19.36	65.05
2026-01-18	3	22.93	26.82	115.29
2026-01-25	4	23.08	30.83	164.11
2026-02-01	5	22.50	34.54	222.35
2026-02-08	6	22.22	39.24	300.28
2026-02-15	7	22.20	42.06	377.45

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

# Historico ignorado
history_ignored <- df_base |>
    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 (SARIMAX)"
)

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

p_future <- ggplot(viz_df, aes(x = Date, y = Value, color = Type)) +
    geom_line(linewidth = 0.8) +
    geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = Type), alpha = 0.2, color = NA) +
    scale_color_manual(values = c("Histórico" = "black", "Previsão (SARIMAX)" = "blue", "Provisório (Ignorado)" = "gray")) +
    scale_fill_manual(values = c("Histórico" = NA, "Previsão (SARIMAX)" = "blue", "Provisório (Ignorado)" = NA)) +
    labs(
        title = "Previsão SARIMAX - Próximas 8 Semanas",
        subtitle = "Com intervalos de confiança (95%)",
        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, "sarimax.rds"))
cat("Modelo salvo em data/models/sarimax.rds\n")

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

7.2 Atualizar Leaderboard

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

# Criar entrada base (compatível com estrutura existente)
new_entry <- data.frame(
    model_name = "sarimax",
    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)
    # Usar bind_rows que lida com colunas diferentes graciosamente
    model_metrics <- dplyr::bind_rows(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)") |>
    kable_styling(bootstrap_options = c("striped", "hover"))

Leaderboard de Modelos (Melhores Métricas Recentes)

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

8. Conclusão

O modelo SARIMAX combina a robustez da estrutura SARIMA com a flexibilidade de regressores externos, incluindo:

• Harmônicos sazonais para capturar padrões cíclicos
• Variáveis climáticas defasadas respeitando o ciclo biológico do vetor
• Lags autoregressivos atualizados dinamicamente nas previsões

Pontos de Destaque

1. Prevenção de Data Leakage: Features recalculadas por fold no backtesting
2. Paralelismo: Redução significativa do tempo de execução com furrr
3. Box-Cox: Transformação logarítmica garante previsões não-negativas
4. Previsão Recursiva: Lags são atualizados dinamicamente em previsões multi-passo

Pontos de Atenção

• A seleção automática de parâmetros (auto.arima) pode ser computacionalmente custosa
• Coeficientes de regressores devem ser monitorados para estabilidade
• O modelo assume relação linear entre regressores e a série diferenciada