Relatório de Modelo: TBATS

Trigonometric Box-Cox ARMA Trend Seasonal

Authors

UIVS - Unidade de Inteligência em Vigilância em Saúde

Caio Sain Vallio

Published

11/01/2026

1. Introdução

1.1 Objetivo

Este relatório avalia o modelo TBATS (Trigonometric Box-Cox ARMA Trend Seasonal) para previsão de incidência de dengue. O TBATS é especialmente adequado para séries com sazonalidade complexa que muda de forma ao longo do tempo.

1.2 Metodologia

Modelo: TBATS (Auto-ajustado)
Estratégia: Rolling Replication (Janela deslizante)
Horizontes de Previsão: 4, 6 e 8 semanas
Comparativo: SNaive, SARIMA, SARIMAX
Métricas: RMSE, MAE, MAPE, RMSLE, MASE

1.3 Sobre o TBATS

O TBATS é uma extensão do modelo de suavização exponencial que inclui:

Componentes do TBATS

Trigonometric: Representação de sazonalidade via séries de Fourier
Box-Cox: Transformação para estabilizar variância
ARMA: Erros modelados como processo ARMA
Trend: Componente de tendência (possivelmente amortecida)
Seasonal: Múltiplos padrões sazonais simultâneos

1.3.1 Vantagens e Ajustes para o Novo Regime

Sazonalidade Flexível: Captura padrões que mudam de intensidade.
Adaptação a Surtos: Ao permitir que o modelo escolha entre tendência amortecida ou não, possibilitamos a captura de crescimentos explosivos (mudança de regime 2024).
Transformação Box-Cox: Essencial para lidar com a magnitude sem precedentes (“explosão de variância”).
Sem Necessidade de Features: Foca na dinâmica interna da série.

2. Dados e Preparação

Code

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

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

Code

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

#> [2026-01-11 15:18:58] INFO: Agregando dados por estado... 
#> [2026-01-11 15:18:58] INFO: Dados agregados: 626 semanas

Code

# 3. Preprocessamento
df <- preprocess_data(df_state)

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

Code

# TBATS não precisa de features - modela sazonalidade internamente
df_base <- df |>
    mutate(data_iniSE = as.Date(data_iniSE)) |>
    arrange(data_iniSE)

cat("Dados carregados:", nrow(df_base), "observações\n")

#> Dados carregados: 625 observações

Code

cat("Período:", min(df_base$data_iniSE), "a", max(df_base$data_iniSE), "\n")

#> Período: 16075 a 20443

Code

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 TBATS representa a série como:

\[y_t^{(\lambda)} = l_{t-1} + \phi b_{t-1} + \sum_{i=1}^{T} s_{t-1}^{(i)} + d_t\]

Onde:

\(y_t^{(\lambda)}\) = série transformada (Box-Cox)
\(l_t\) = nível
\(b_t\) = tendência (com amortecimento \(\phi\))
\(s_t^{(i)}\) = componente sazonal \(i\) (representado por Fourier)
\(d_t\) = resíduo ARMA

Code

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

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

#> Modelo: tbats
#> Família: statistical
#> Descrição: TBATS - Trigonometric Box-Cox ARMA Trend Seasonal

Code

cat("\n### Configuração Refinada:\n")

#> 
#> ### Configuração Refinada:

Code

cat("- Período sazonal: 52 semanas\n")

#> - Período sazonal: 52 semanas

Code

cat("- Transformação Box-Cox: Habilitada (Estabilização de variância)\n")

#> - Transformação Box-Cox: Habilitada (Estabilização de variância)

Code

cat("- Tendência e Amortecimento: Automático (AIC decide)\n")

#> - Tendência e Amortecimento: Automático (AIC decide)

Code

cat("  *Nota: Restrições relaxadas para capturar mudança de regime*\n")

#>   *Nota: Restrições relaxadas para capturar mudança de regime*

Code

cat("- Erros ARMA: Habilitados\n")

#> - Erros ARMA: Habilitados

4. Backtesting

Utilizando a mesma configuração de backtesting dos relatórios anteriores.

Code

# 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"
    )
)

# Configuração do modelo TBATS
# Configuração do modelo TBATS refinada
# Relaxamos as restrições de tendência (use.trend=NULL, use.damped.trend=NULL)
# para permitir que o modelo capture o "bom" de crescimento explosivo recente
# sem forçar um amortecimento artificial.
model_config <- list(
    target = "est_inc100k",
    seasonal.periods = c(52),
    use.box.cox = TRUE, # Manter Box-Cox (lidar com variância crescente)
    use.trend = NULL, # Deixar AIC decidir (Pode escolher tendência linear explosiva)
    use.damped.trend = NULL # Deixar AIC decidir (Não forçar amortecimento)
)

# TBATS não precisa de features, usar run_backtest normal
# Mas podemos usar run_backtest_parallel para velocidade
has_parallel <- requireNamespace("furrr", quietly = TRUE) &&
    requireNamespace("future", quietly = TRUE)

if (has_parallel) {
    cat("✅ Usando execução PARALELA\n")

    # Para TBATS, não precisamos recalcular features por fold
    # Mas usamos a versão parallel para consistência
    bt_result <- run_backtest_parallel(
        model_name = model_name,
        data_base = df_base,
        config = backtest_config,
        model_config = model_config,
        verbose = TRUE
    )
} else {
    cat("⚠️ Usando execução SEQUENCIAL\n")

    bt_result <- run_backtest(
        model_name = model_name,
        data = df_base,
        config = backtest_config,
        model_config = model_config,
        verbose = TRUE
    )
}

#> ✅ Usando execução PARALELA

#> [2026-01-11 15:18:58] INFO: ========== INICIO: Backtest PARALLEL: tbats ========== 
#> [2026-01-11 15:18:58] INFO: Origens: 90 | Horizontes: 4, 6, 8 | Workers: 15 
#> [2026-01-11 15:18:58] INFO: Modo: Features recalculadas por fold (sem leakage)

#> [2026-01-11 15:19:26] INFO: Previsoes geradas: 270 
#> [2026-01-11 15:19:26] INFO: Erros: 0 
#> [2026-01-11 15:19:26] INFO: Tempo total: 27.2 segundos 
#> [2026-01-11 15:19:26] INFO: ========== FIM: Backtest PARALLEL: tbats (27.24s) ==========

Code

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

#> ⏱️ Tempo total: 27.2 segundos

Code

print(bt_result)

#> 
#> === Resultado de Backtest ===
#> 
#> Modelo: tbats 
#> Origens: 90 
#> Previsoes: 270 
#> Duracao: 27.2 s
#> 
#> Metricas por Horizonte:
#>  h      mae     rmse     mape    smape      mase     rmsle  n
#>  4 10.39904 20.46878 23.29374 23.06774 0.3070472 0.3153439 90
#>  6 15.21271 28.22691 34.09684 33.03823 0.4491781 0.4554588 90
#>  8 20.02781 40.10261 47.91038 41.05286 0.5913510 0.5812118 90

5. Avaliação de Desempenho

5.1 Métricas por Horizonte

Code

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 TBATS por Horizonte") |>
    kable_styling(bootstrap_options = c("striped", "hover"))

Métricas de Performance do TBATS por Horizonte
Horizonte	RMSE	MAE	MAPE	RMSLE	MASE
4 semanas	20.4688	10.3990	23.2937	0.3153	0.3070
6 semanas	28.2269	15.2127	34.0968	0.4555	0.4492
8 semanas	40.1026	20.0278	47.9104	0.5812	0.5914

5.2 Comparação com Modelos Anteriores

Code

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

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

    # Últimas execuções de cada modelo
    comparison_data <- model_metrics |>
        group_by(model_name) |>
        slice(which.max(execution_date)) |>
        ungroup() |>
        select(model_name, horizon_4w_rmse, horizon_4w_mae, horizon_4w_mape, horizon_4w_rmsle)

    # Adicionar TBATS atual
    tbats_row <- data.frame(
        model_name = "tbats",
        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"]
    )

    all_models <- comparison_data |>
        filter(model_name != "tbats") |>
        dplyr::bind_rows(tbats_row) |>
        arrange(horizon_4w_rmse)

    all_models |>
        mutate(across(where(is.numeric), \(x) round(x, 4))) |>
        rename(
            Modelo = model_name,
            RMSE = horizon_4w_rmse,
            MAE = horizon_4w_mae,
            MAPE = horizon_4w_mape,
            RMSLE = horizon_4w_rmsle
        ) |>
        kable(caption = "Comparativo de Modelos - Horizonte 4 semanas") |>
        kable_styling(bootstrap_options = c("striped", "hover")) |>
        row_spec(which(all_models$model_name == "tbats"), bold = TRUE, background = "#E8F4F8")
} else {
    cat("Arquivo de métricas anteriores não encontrado.")
}

Comparativo de Modelos - Horizonte 4 semanas
Modelo	RMSE	MAE	MAPE	RMSLE
tbats	20.4688	10.3990	23.2937	0.3153
sarima	36.3290	14.2967	25.1256	0.2930
sarimax	90.6956	34.7848	43.8269	0.7821
seasonal_mean	91.2846	39.0308	65.0013	1.0641
seasonal_naive	96.3101	42.7400	66.1433	0.9683

5.3 Análise Visual Multi-Horizonte

Code

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")
)

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 TBATS: Previsões Multi-Horizonte",
        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

Code

preds_all <- bt_result$predictions
preds_all$regime <- classify_regime(preds_all$actual)

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")

6. Previsão Futura

Code

# Definir atraso de notificacao (semanas a ignorar)
delay_weeks <- 2

# Filtrar dados para treino final (removendo ultimas semanas incompletas)
last_reliable_date <- max(df_base$data_iniSE) - (delay_weeks * 7)
df_final <- df_base |>
    filter(data_iniSE <= last_reliable_date)

cat(sprintf("Ignorando últimas %d semanas devido ao atraso de notificação.\n", delay_weeks))

#> Ignorando últimas 2 semanas devido ao atraso de notificação.

Code

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

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

Code

# Treinar modelo final com dados confiáveis
model_fit <- model$fit(df_final, model_config)

# Parâmetros do modelo
params <- model$get_params(model_fit)

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

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

Code

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

#> - **Lambda (Box-Cox)**: 0.0090

Code

cat(sprintf("- **Alpha (nível)**: %.4f\n", params$alpha))

#> - **Alpha (nível)**: 0.5689

Code

if (!is.null(params$beta)) {
    cat(sprintf("- **Beta (tendência)**: %.4f\n", params$beta))
}
if (!is.null(params$damping)) {
    cat(sprintf("- **Damping**: %.4f\n", params$damping))
}
cat(sprintf("- **AIC**: %.2f\n", params$aic))

#> - **AIC**: 4759.61

Code

# Decomposição do TBATS
decomp <- model$decompose(model_fit)

# Visualizar componentes
if (!is.null(decomp$components_df)) {
    decomp_long <- decomp$components_df |>
        mutate(t = row_number()) |>
        pivot_longer(cols = -t, names_to = "component", values_to = "value")

    ggplot(decomp_long, aes(x = t, y = value)) +
        geom_line(color = "steelblue") +
        facet_wrap(~component, scales = "free_y", ncol = 1) +
        labs(title = "Decomposição TBATS", x = "Tempo", y = "Valor") +
        theme_minimal()
}

Code

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

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

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

Previsão TBATS - Próximas 8 Semanas
Data	Semana_Epidem	Previsao	Inferior_95	Superior_95
2025-12-28	53	24.54	19.53	30.84
2026-01-04	1	26.85	19.03	37.86
2026-01-11	2	30.88	19.75	48.22
2026-01-18	3	36.71	21.21	63.39
2026-01-25	4	43.95	22.85	84.20
2026-02-01	5	51.80	24.15	110.54
2026-02-08	6	59.53	24.84	141.71
2026-02-15	7	67.03	25.09	177.57

Code

# Dados historicos usados no treino
history_train <- df_final |>
    tail(104) |>
    select(Date = data_iniSE, Value = est_inc100k) |>
    mutate(Type = "Histórico", Lower = NA, Upper = NA)

# Dados ignorados (Provisorio)
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 (TBATS)"
)

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 (TBATS)" = "purple", "Provisório (Ignorado)" = "gray")) +
    scale_fill_manual(values = c("Histórico" = NA, "Previsão (TBATS)" = "purple", "Provisório (Ignorado)" = NA)) +
    labs(
        title = "Previsão TBATS - 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

Code

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, "tbats.rds"))
cat("Modelo salvo em data/models/tbats.rds\n")

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

7.2 Atualizar Leaderboard

Code

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

new_entry <- data.frame(
    model_name = "tbats",
    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)
    model_metrics <- dplyr::bind_rows(model_metrics, new_entry)
} else {
    model_metrics <- new_entry
}

save(model_metrics, file = metrics_path)

# 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") |>
    kable_styling(bootstrap_options = c("striped", "hover"))

Leaderboard de Modelos
Modelo	RMSE_4w	MAPE_4w	RMSLE_4w
tbats	20.46878	23.29374	0.3153439
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 TBATS oferece uma abordagem robusta para capturar sazonalidades complexas sem necessidade de features externas.

Pontos de Destaque

Modelagem Automática: TBATS seleciona automaticamente os melhores parâmetros
Sazonalidade Flexível: Séries de Fourier capturam padrões que mudam de forma
Box-Cox: Transformação garante previsões não-negativas
Tendência Amortecida: Evita extrapolação irreal de tendências

Limitações

Não incorpora variáveis exógenas (clima, eventos)
Pode ser lento para séries muito longas
Assume que padrões passados continuam no futuro