Run the model search

best_models <- find_bestfit(
  data = epi_data,
  response = "y",
  group = "epi_id",
  time = "time",
  vars = c("tmean","rain","wetness"),
  max_lag = 85,
  df_var_grid = c(2,3,4),
  df_lag_grid = c(2,3,4),
  model_engine = "glmmTMB",
  family = "beta",
  keep_fits = TRUE)

plot_performance(
  object = best_models,
  model_id = c(5,12,19,26,28,33,35,40,42,47,49,54,56,61,63),
  metrics = c("CCC"),
  x = "observed",
  y = "predicted",
  scale_factor = 1,
  ncol = 5,
  x_lab = "Observed ",
  y_lab = "Predicted",
  size = 1.4
) +
  ggplot2::theme(
    text = ggplot2::element_text(
      size = 10,
      face = "bold"
    )
  )

Observed versus predicted performance for selected candidate models

Inspect results

Depending on the selected validation settings, the resulting object can include:

• a model identifier;
• the fitted predictor combination;
• exposure and lag spline settings;
• the modeling engine and response family;
• the Concordance Correlation Coefficient;
• the Root Mean Square Error;
• the Mean Absolute Error;
• bias-related metrics;
• the fitted model object when keep_fits = TRUE.

Top-performing models

The following example ranks models by decreasing CCC. This is one possible ranking criterion, not a universal rule. Candidate models should also be examined using error metrics, bias, convergence diagnostics, model complexity, and the scientific plausibility of the exposure-lag structure.

best_models %>%
  arrange(desc(CCC)) %>%
  select(
    model_id,
    vars,
    CCC,
    RMSE,
    MAE
  ) %>%
  head(10)
#> # A tibble: 10 × 5
#>    model_id vars                     CCC   RMSE    MAE
#>       <dbl> <chr>                  <dbl>  <dbl>  <dbl>
#>  1       63 tmean + rain + wetness 0.980 0.0529 0.0382
#>  2       42 tmean + rain + wetness 0.978 0.0550 0.0395
#>  3       56 tmean + rain + wetness 0.970 0.0642 0.0475
#>  4       35 tmean + rain + wetness 0.969 0.0656 0.0483
#>  5       49 tmean + rain + wetness 0.952 0.0803 0.0577
#>  6       28 tmean + rain + wetness 0.949 0.0825 0.0589
#>  7       40 tmean + wetness        0.859 0.132  0.0988
#>  8       61 tmean + wetness        0.858 0.132  0.0998
#>  9       19 tmean + wetness        0.851 0.134  0.102 
#> 10       47 tmean + wetness        0.848 0.136  0.104

Concordance Correlation Coefficient

CCC evaluates agreement between observed and predicted outcomes by combining precision and accuracy.

Values closer to 1 indicate stronger concordance. However, CCC should be interpreted together with error and bias metrics rather than used as the only selection criterion.

Root Mean Square Error

RMSE summarizes the magnitude of prediction errors while assigning greater influence to larger errors because residuals are squared before averaging.

Lower values indicate better predictive performance when models are evaluated on the same response scale and validation data.

Mean Absolute Error

MAE is the average absolute difference between observed and predicted outcomes.

Lower values indicate better predictive performance. Because MAE remains on the response scale, it is often easier to interpret directly than RMSE.

Inspect results

The ensemble summary can be extracted from the fitted ensemble object:

ens_unweighted$ensemble_summary

Inspect results

The ensemble summary can be extracted from the fitted ensemble object:

ens_weighted$ensemble_summary

Inspect results

The ensemble summary can be extracted from the fitted ensemble object:

ens_stacked$ensemble_summary

Why construct a lag ensemble?

Model-level ensembles combine predictions from multiple candidate models. However, models may also differ in how they distribute the contribution of an exposure across the fitted lag interval.

A lag ensemble combines the lag-specific epidemiological contribution indices obtained from the selected models. For lag \(l\), the ensemble contribution is defined as:

\[ ECI_{\mathrm{ensemble},l} = \sum_m \alpha_m ECI_{m,l}, \]

where \(ECI_{m,l}\) is the lag-specific weighted ECI from model \(m\), and \(\alpha_m\) is the model weight derived from out-of-fold model performance.

The model weights therefore come from cross-validation, whereas the lag-specific ECI values are calculated from the selected models refitted to the complete dataset.

The fitted models are not used to recalculate the cross-validation metrics. The model ranking and ensemble weights remain based on the out-of-fold predictions produced by find_bestfit(). The retained full-data fits are used only to obtain the final lag-specific decompositions for the selected models.

Constructing unweighted and weighted lag ensembles

The following example selects the three highest-ranked models and constructs two lag ensembles:

• "unweighted" gives equal weight to each selected model;
• "weighted" derives model weights from the selected performance metric.

To keep the example computationally efficient, the lag ensemble is demonstrated using one complete epidemic profile. Users may instead supply all available epidemic profiles, provided that each profile contains the complete max_lag + 1 observations required by the fitted models. The computational cost will increase with the number of epidemics, exposures, selected models, and lags. The results shown here are therefore specific to Epidemic 1.

epi_data1 = epi_data |> 
  filter(epi_id == 1)

lag_ensemble <- ensemble_bestfit(
  bestfit = best_models,
  data = epi_data1,
  group = "epi_id",
  var = c(
    "tmean",
    "rain",
    "wetness"
  ),
  ensemble_scope = "lag",
  method = c(
    "unweighted",
    "weighted"
  ),
  top_n = 3,
  weight_metric = NULL,
  weight_transform = "softmax",
  lag_group_cols = c(
    "epi_id",
    "var",
    "lag"
  ),
  lg_strategy = "requested_available",
  rl_weights = TRUE,
  compute_ecilag_args = list(
    uncertainty = FALSE
  ),
  verbose = TRUE
)

When weight_metric = NULL, ensemble_bestfit() uses the ranking metric stored by find_bestfit(). The direction of that metric is handled automatically: metrics such as CCC are maximized, whereas metrics such as RMSE are minimized.

The argument:

lg_strategy = "requested_available"

uses each requested exposure when that exposure is available in the selected model. Because candidate models may contain different combinations of exposures, not every model necessarily contributes to every exposure-specific lag ensemble.

With:

rl_weights = TRUE

the model weights are renormalized among the models available for each exposure-lag unit. This ensures that the available model weights sum to one within each lag-specific ensemble calculation.

Inspecting selected models and weights

The selected candidate models can be inspected directly:

lag_ensemble$selected_models

The model weights used by each ensemble method are stored separately:

lag_ensemble$model_weights

For the unweighted ensemble, the selected models receive equal lag weights. For the weighted ensemble, the weights are derived from the selected cross-validation performance metric.

These weights describe the relative influence of the selected models in the lag ensemble. They are not lag-specific DLNM coefficients.

Inspecting ensemble lag contributions

The combined lag-specific contributions are returned in:

lag_ensemble$ensemble_by_lag

Inspect the first rows:

The principal output columns are:

• var: exposure variable;
• lag: retrospective lag, where lag 0 is the most recent exposure;
• method: ensemble method;
• ECI_weighted_ens: signed ensemble ECI contribution at that lag;
• ECI_percent_ens: absolute percentage contribution of that lag within the corresponding exposure and ensemble method;
• n_models: number of models contributing to that lag unit;
• raw_model_weight_sum: total original model-weight mass available for that lag unit;
• weights_renormalized: whether weights were renormalized among the available models;
• models_used: selected models contributing to the lag unit.

The signed value ECI_weighted_ens preserves the direction of the ensemble contribution. In contrast, ECI_percent_ens is based on the absolute contribution magnitude and describes how the total lag contribution is distributed across the fitted lag interval.

Visualizing signed lag contributions

ensemble_by_lag = lag_ensemble$ensemble_by_lag

lag_ensemble_signed_plot <- ggplot(
  ensemble_by_lag,
  aes(
    x = lag,
    y = ECI_weighted_ens,
    color = method
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey50",
    linetype = "dashed",
    linewidth = 0.5
  ) +
  geom_line(
    linewidth = 1
  ) +
  facet_wrap(
    vars(var),
    ncol = 1,
    scales = "free_y",
    labeller = as_labeller(
      c(
        tmean = "Mean temperature",
        rain = "Rainfall",
        wetness = "Leaf wetness"
      )
    )
  ) +
  scale_x_reverse(
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      unweighted = "#4C72B0",
      weighted = "#C44E52"
    ),
    labels = c(
      unweighted = "Unweighted",
      weighted = "Weighted"
    )
  ) +
  theme_bw() +
  labs(
    x = "Retrospective lag",
    y = "Ensemble weighted ECI",
    color = "Ensemble method"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    strip.background = element_rect(
      fill = "white",
      color = "black"
    ),
    strip.text = element_text(
      face = "bold"
    ),
    legend.position = "top"
  )

lag_ensemble_signed_plot

The lag axis is displayed retrospectively. Lag 0 corresponds to the most recent exposure observation, whereas increasing lag values represent progressively older exposure conditions.

Positive and negative values describe the direction of the combined lag-specific contribution. The weighted and unweighted curves may differ when the selected models receive substantially different performance-based weights.

Visualizing relative lag contributions

The percentage contribution provides a complementary interpretation based on the absolute contribution magnitude:

lag_ensemble_percent_plot <- ggplot(
  ensemble_by_lag,
  aes(
    x = lag,
    y = ECI_percent_ens,
    color = method
  )
) +
  geom_line(
    linewidth = 1
  ) +
  facet_wrap(
    vars(var),
    ncol = 1,
    scales = "free_y",
    labeller = as_labeller(
      c(
        tmean = "Mean temperature",
        rain = "Rainfall",
        wetness = "Leaf wetness"
      )
    )
  ) +
  scale_x_reverse(
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      unweighted = "#4C72B0",
      weighted = "#C44E52"
    ),
    labels = c(
      unweighted = "Unweighted",
      weighted = "Weighted"
    )
  ) +
  theme_bw() +
  labs(
    x = "Retrospective lag",
    y = "Absolute lag contribution (%)",
    color = "Ensemble method"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    strip.background = element_rect(
      fill = "white",
      color = "black"
    ),
    strip.text = element_text(
      face = "bold"
    ),
    legend.position = "top"
  )

lag_ensemble_percent_plot

ECI_percent_ens is calculated from the absolute lag contributions within each exposure and ensemble method. Consequently, the percentages describe where the ensemble contribution is concentrated across lags, but they do not preserve the positive or negative direction of that contribution.

The signed and percentage plots should therefore be interpreted together:

• ECI_weighted_ens describes the direction and magnitude of the combined contribution;
• ECI_percent_ens describes the relative concentration of absolute contribution across lags.

Using strict exposure availability

The previous example used:

lg_strategy = "requested_available"

which allows different selected models to contribute to different exposures when their fitted variable sets differ.

For a stricter comparison, the user can require every selected model to contain all requested exposures:

lag_ensemble_strict <- ensemble_bestfit(
  bestfit = best_models,
  data = epi_data1,
  group = "epi_id",
  var = c(
    "tmean",
    "rain",
    "wetness"
  ),
  ensemble_scope = "lag",
  method = c(
    "unweighted",
    "weighted"
  ),
  top_n = 3,
  weight_transform = "softmax",
  lag_group_cols = c(
    "epi_id",
    "var",
    "lag"
  ),
  lg_strategy = "strict",
  rl_weights = TRUE,
  compute_ecilag_args = list(
    uncertainty = FALSE
  )
)

This strict mode stops when a selected model does not contain one of the requested exposures or when a lag unit is not available from every selected model. It is useful when the analyst wants the same model set to support every exposure-specific lag ensemble.

Lag ensembles should be interpreted as model-performance-weighted summaries of lag-specific ECI contributions. They do not replace external validation, do not represent raw DLNM coefficients, and should not be interpreted as causal attribution across lags.

Avoiding optimistic ensemble evaluation

Ensemble weights should be estimated from out-of-sample or cross-validated predictions whenever possible. Estimating weights and reporting performance on the same fitted predictions can produce optimistic results.

The final evaluation should use observations that were not reused to select candidate models or estimate ensemble weights. If a completely independent test set is unavailable, nested or repeated cross-validation may provide a more defensible assessment.