--- title: "Comparing Multiple Forecasters with SMCS" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Comparing Multiple Forecasters with SMCS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4 ) ``` ```{r setup} library(seqcomp) ``` ## Overview While `compare_forecasts()` evaluates exactly two models, real-world applications often require comparing an arbitrary number of candidate models ($m \ge 2$) simultaneously. Evaluating multiple models sequentially introduces a multiple-testing problem: if you run enough pairwise comparisons over enough time steps, you will eventually find false "statistically significant" differences by pure chance. To solve this, `seqcomp` implements **Sequential Model Confidence Sets (SMCS)** (Arnold et al., 2026). The SMCS is the set of models that have not yet been confidently beaten by any other model in the candidate pool. The package guarantees that the true best model(s) will remain in the SMCS with probability at least $1 - \alpha$ over the entire evaluation period. *(If you are only comparing two forecasters, we recommend starting with the [introductory vignette](seqcomp.html).)* ## Three Notions of Superiority When comparing multiple models, `seqcomp` tracks three different definitions of what makes a model "superior", following Arnold, Gavrilopoulos, Schulz and Ziegel (2026): 1. **The Strong Null (`smcs_strong(method = "betting")`)**: A model is strongly superior when, at every time step, its conditional expected score is at least as large as that of every competitor. The strong-null SMCS uses product-form betting e-processes and permanent exclusion. 2. **The Uniformly Weak Null (`smcs_strong(method = "mixture")`)**: A model is uniformly weakly superior if its running average performance against every competitor, taken over the *entire* evaluation horizon so far, never falls behind. This uses the same closed-testing machinery as the strong null above (just built from exponential-mixture e-processes instead of the product-form martingale) so it also maintains a running intersection: once a model is excluded, it stays excluded. 3. **The (Time-Varying) Weak Null (`smcs_weak()`)**: A model is weakly superior *at time t* if its average performance against every competitor, taken over the horizon up to t, is not behind. Unlike the two notions above, this is re-evaluated fresh at each t from Bonferroni-adjusted pairwise confidence sequences, rather than from closed testing over e-processes. This is what lets a model exit and later re-enter the set as its average performance recovers. The strong and uniformly weak nulls never allow re-entry. ## A simple multi-model example Let's generate a sequence of binary outcomes and three candidate forecasters: * `M1`: A highly accurate forecaster. * `M2`: A mediocre forecaster that guesses randomly. * `M3`: A terrible forecaster that is consistently wrong. ```{r} set.seed(2026) T_sim <- 300 y <- rbinom(T_sim, size = 1, prob = 0.7) # Create 3 forecasters forecasts <- matrix(0, nrow = T_sim, ncol = 3) colnames(forecasts) <- c("M1", "M2", "M3") forecasts[, "M1"] <- ifelse(y == 1, 0.8, 0.2) # Highly accurate forecasts[, "M2"] <- runif(T_sim, 0.4, 0.6) # Mediocre/random forecasts[, "M3"] <- ifelse(y == 1, 0.2, 0.8) # Consistently wrong head(forecasts) ``` This toy example is deliberately simple and outcome-conditioned. In a real application, forecasts would come from two forecasting models, analysts, or institutions, before observing `y`. ## Compare the forecasts The easiest workflow for evaluating three or more models is to use `smcs_compare()`. ```{r} multi_cmp <- smcs_compare( forecasts = forecasts, outcomes = y, scoring_rule = "brier", alpha = 0.05 ) # Printing the object provides a clean summary of which models # were excluded and when they first dropped out of the set. multi_cmp ``` While the printed summary gives a quick overview, the underlying object is a list containing the full time-indexed matrices for further analysis (alongside metadata such as `alpha` and `scoring_rule`): * `multi_cmp$scores`: The evaluated pointwise scores for each model. * `multi_cmp$smcs_strong`: A logical matrix tracking whether each model belongs to the strong-null SMCS at time $t$. * `multi_cmp$smcs_uniform_weak`: A logical matrix tracking whether each model belongs to the uniformly weak SMCS at time $t$. * `multi_cmp$smcs_weak`: A logical matrix tracking whether each model belongs to the weak-null SMCS at time $t$. ## Interpreting the SMCS output `smcs_strong()` maintains a **running intersection** over time: once a model is shown to violate the null against some competitor, it is excluded permanently. This matches the strong/uniformly-weak SMCS construction of Arnold, Gavrilopoulos, Schulz and Ziegel (2026, Section 3.2). The underlying hypothesis of `smcs_weak()` (the *time-varying* weak null, Section 3.3) is defined *at* each time $t$, not cumulatively, so a model's average performance can legitimately recover after a bad stretch. `smcs_weak()` therefore returns `smcs`, which can re-admit models that were previously excluded. Let's look at the strong SMCS at the end of the evaluation: ```{r} tail(multi_cmp$smcs_strong) ``` We can visualize how the confidence set shrinks over time by plotting whether a model is `TRUE` (Included) or `FALSE` (Excluded): ```{r} par(mfrow = c(3, 1), mar = c(2, 4, 2, 1)) colors <- c("blue", "gray", "red") for (i in 1:3) { plot( 1:T_sim, multi_cmp$smcs_strong[, i], type = "s", col = colors[i], lwd = 2, ylim = c(-0.1, 1.1), yaxt = "n", ylab = "In SMCS?", main = paste("Model:", colnames(forecasts)[i]) ) axis(2, at = c(0, 1), labels = c("Excluded", "Included"), las = 2) } par(mfrow = c(1, 1)) ``` Notice how quickly the set shrinks: 1. `M3` (the terrible model) is confidently beaten by `M1` almost immediately and drops out of the set. 2. `M2` (the mediocre model) survives a bit longer, but as evidence accumulates, it too is permanently excluded. 3. `M1` (the true best model) remains in the SMCS for the entire duration. ## Using lower-level functions directly Just like the two-model case, `seqcomp` exposes the underlying multi-model primitives if you need custom bounds, different scoring rules, or adaptive betting fractions. You can compute your own score matrices and pass them directly to `smcs_strong()` or `smcs_weak()`. ```{r} # Manually compute Brier scores scores_mat <- matrix(0, nrow = T_sim, ncol = 3) for(i in 1:3) scores_mat[, i] <- brier_score(forecasts[, i], y) # Construct the Weak SMCS directly # Brier score differences are bounded in [-1, 1], so we use c_param = 2 res_weak <- smcs_weak( scores = scores_mat, alpha = 0.05, cs_method = "bernstein", c_param = 2 ) tail(res_weak$smcs) ``` ## Adaptive betting fractions for the strong null The strong-null e-process used by `smcs_strong(method = "betting")` is a product-form betting martingale, $E_t = \prod_{r \le t}(1 + \lambda_r d_r)$. Arnold et al. (2026) only require that $\lambda_r$ be predictable and lie in $[0, 1/c_r]$ at every round; `smcs_strong()`'s default, if `lambda_param` is left `NULL`, is the conservative fixed fraction $\lambda_r = 1/(2c_r)$. `seqcomp` also provides two adaptive betting-fraction rules, `lambda_betting_agrapa()` and `lambda_betting_ons()`, which adapt the aGRAPA and ONS-m algorithms of Waudby-Smith and Ramdas (2024) to Arnold's bounded strong-null setting. These adaptations are **not** given in Arnold et al. (2026) itself; they are original extensions implemented in this package (see the function documentation for the derivation. Additionally, see [the adaptive betting vignette](adaptive_betting.html) for a comparison of performances under different temporal structures). `build_agrapa_betting_array()` and `build_ons_betting_array()` build the full $T \times m \times m$ `lambda_param` array these rules require: ```{r} lam_agrapa <- build_agrapa_betting_array(scores_mat, c_mat = 2) res_strong_agrapa <- smcs_strong( scores_mat, alpha = 0.05, method = "betting", c_param = 2, lambda_param = lam_agrapa ) tail(res_strong_agrapa$smcs) ``` `build_ons_betting_array()` is a drop-in alternative using ONS-m instead of aGRAPA. Both accept a `period` argument: if competitors differ systematically by a known periodic covariate (for example, day-of-week seasonality), passing `period` conditions each rule on its own interleaved sub-stream rather than pooling across the full history. Pooling otherwise dilutes the adaptive rule's ability to react to the periodic pattern, since the pooled running mean/variance stays governed by the majority (non-anomalous) observations: ```{r, eval = FALSE} lam_agrapa_p7 <- build_agrapa_betting_array(scores_mat, c_mat = 2, period = 7) ``` This `period` currently only works with the aGRAPA and ONS-m builders. ## Conditionally bounded scores (Tick loss) Some scoring rules, like tick loss for quantile forecasting, are unbounded globally but bounded *conditionally* based on the distance between the forecasts. For these rules, the strong SMCS requires a dynamically updating 3D array of bounds and betting fractions. The wrapper `smcs_compare()` handles this automatically when `scoring_rule = "tick"`, or you can construct the arrays manually using `build_quantile_betting_arrays()`. ```{r} set.seed(7) tau <- 0.5 truth <- rnorm(T_sim) q_forecasts <- cbind( M1 = truth + rnorm(T_sim, sd = 0.1), # accurate M2 = rep(0, T_sim), # uninformative M3 = truth + rnorm(T_sim, sd = 0.1) + 1 # biased ) tick_cmp <- smcs_compare( forecasts = q_forecasts, outcomes = truth, scoring_rule = "tick", tau = tau, alpha = 0.05 ) tail(tick_cmp$smcs_strong) ``` `smcs_uniform_weak` and `smcs_weak` are `NULL` here, since tick-loss score differences are only conditionally bounded, not uniformly bounded, while the current mixture and confidence-sequence implementations require a fixed uniform bound.