modelimportance: Evaluating model importance within a multi-model ensemble in R

2.1 Dependencies and related software

The modelimportance package is designed to work with the hubverse framework and, accordingly, depends on several packages in the hubverse ecosystem, such as hubUtils (Krystalli and Shandross (2025)), hubEnsembles (Shandross et al. (2025)), and hubEvals (Reich et al. (2025)). modelimportance uses a model_out_tbl S3 class as the model output format defined in hubUtils, which consists of utility functions to standardize prediction files and data formats (details in Section 2.2). Ensembling predictions from multiple models relies on hubEnsembles, which offers a broadly applicable framework to construct multi-model ensembles using various ensemble methods. Calculation of forecast accuracy using various metrics is based on hubEvals, which internally leverages scoringutils.

We derived example datasets used for testing and demonstration purposes (see Section 6) from the hubExamples (Ray et al. (2025)) package, which provides example datasets in the hubverse format. The dataset were locally stored for reproducibility within the package distribution, instead of requiring the hubExamples package as a dependency.

For visualization of model importance scores, standard ggplot2 (Wickham 2016) functions can be applied. Parallel computing is supported by the furrr (Vaughan and Dancho 2022) and future (Bengtsson 2021) packages, which allow users to specify the number of cores to use for parallel processing when calculating model importance scores across many tasks.

2.2 Model output format

Model outputs are structured in a tabular format designed specifically for predictions, which is a formal S3 object called model_out_tbl. In the hubverse standard, each row represents an individual prediction or a component of a prediction for a single task. More details about that prediction or prediction component are described in multiple columns through which one can identify the unique label assigned to each forecasting model, task characteristics, prediction representation type, and predicted values (Shandross et al. 2026). To elaborate on the task characteristics, each prediction task corresponds to a specific forecasting problem and it can be described by a set of task ID variables. Examples of such variables include a date on which forecasts are generated, the target to predict (e.g., flu-related incident deaths, cases, or hospitalizations), and the prediction horizon, which is the length of time into the future from the point when a model generate its forecast, for a specific location on a certain target date. Table 1 illustrates short-term forecasts of weekly incident influenza hospitalizations in the US for Massachusetts, generated by the model ‘Flusight-baseline’ on December 17, 2022, in the model_out_tbl format. The model_id column lists a uniquely identified model name. The reference_date, target, horizon, location, and target_end_date columns are collectively referred to as the task ID variables, which together defines the task characteristics. Note that the forecast generation date and the target date for which the prediction is made are mapped to the reference_date and target_end_date columns, respectively, and the location is represented based on the FIPS code (e.g., ‘25’ for Massachusetts). The time length to the target_end_date, which is the number of weeks ahead from the reference_date, is indicated in the horizon column. The prediction representation is specified as ‘quantile’ in the output_type column, and details are represented in the output_type_id column with seven quantiles of 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, and 0.95 for each target end date. The predicted value corresponding to each quantile is recorded in the value column.

Figure 1 visualizes the information on the prediction task provided by Table 1 for three models. For each model, the quantile-based forecasts are shown for the target end dates of December 24, 2022 (horizon 1), December 31, 2022 (horizon 2), and January 07, 2023 (horizon 3), which were made on December 17, 2022 based on the historical data available as of that date. The prediction intervals defined by the lowest and highest quantiles (0.05 and 0.95) represent the uncertainty of the predictions. To give a brief interpretation, the Flusight-baseline model under-predicted the outcomes for the first two target dates (horizon 1 and 2), but it over-predicted the outcome for the last target date (horizon 3). Its prediction intervals are narrow compared to the other two models, which indicates that it is more confident about its predictions. However, two of three prediction intervals (horizons 1 and 2) failed to cover the eventually observed values, implying that the model was overconfident and/or biased.

Figure 1: Example plot of three distributional forecasts corresponding to the model output for incident influenza hospitalizations shown in Table 1. Solid black dots indicate historically available data as of the forecast generation date, and open black circles indicate the eventually observed values. The blue dots represent predictive medians and the blue shaded area represents the corresponding 90% prediction interval defined by the 0.05 and 0.95 quantiles.

2.3 Forecast data representation

Generally, quantitative forecasts can be categorized as being either point forecasts or probabilistic forecasts. For a specific prediction task, point forecasts, represented by a single predicted value, provide a clear and concise prediction, making them easy to interpret and communicate. Probabilistic forecasts, on the other hand, provide a probability distribution over possible future values, which inherently involves uncertainty. They are represented in various ways, such as probability mass functions (pmf), cumulative distribution functions (cdf), samples, or probability quantiles (or intervals). Note that there might be other ways, such as named distributions (e.g. Normal(1,2)) that are not currently supported by hubverse.

The output_type and output_type_id columns in the hubverse model output format specify the forecast structure. In the modelimportance package, model outputs may only contain one output_type of ‘mean’, ‘median’, ‘quantile’, or ‘pmf’: ‘mean’ or ‘median’ for point forecasts and ‘quantile’ or ‘pmf’ for probabilistic forecasts.
As aforementioned, output_type_id column identifies additional detailed information, such as specific quantile levels (e.g., “0.1”, “0.25”, “0.5”, “0.75”, and “0.9”) for the ‘quantile’ output type and categorical values (e.g., “low”, “moderate”, “high”, and “very high”) for the ‘pmf’ output type. The predicted values for pmf are constrained to be between 0 and 1, indicating the probability at each categorical level, while they are unbounded numeric otherwise. Different output types correspond to different scoring rules for evaluating a model’s prediction performance. Table 2 presents the output types and their associated scoring rules supported by the modelimportance package. We note that while the scoring rules are listed in the table using their conventional names, the package uses their negative values for evaluation (e.g., \(-\text{RSE}, -\text{WIS}\)) so that higher scores indicate better performance. This positively oriented scoring rule facilitates the interpretation of importance scores: positive values indicate that the model’s inclusion in the ensemble improves ensemble performance, whereas negative values indicate that it worsens it. Further, using positive values to indicate improved ensemble performance aligns with the convention of Shapely values.

2.4 Oracle output data

The oracle_output_data is a data frame that contains the ground truth values for the variables used to define modeling targets (Consortium of Infectious Disease Modeling Hubs 2024; Kerr et al. 2025). Its name originates from the concept of an oracle making a perfect prediction formatted similarly to model output. This type of data must follow the oracle output format defined in the hubverse standard, which includes independent task ID columns (e.g., location, target_date), the output_type column specifying the output type of the predictions and an oracle_value column for the observed values. As in the forecast data, if the output_type is either "quantile" or "pmf", the output_type_id column is often required to provide further identifying information.

The model_out_tbl and oracle_output_data must have the same task ID columns and output_type, including output_type_id if necessary, as the oracle_output_data is joined with the model_out_tbl based on these fields in order to score the forecast performance.

3.1 Comparison of weighting schemes in LASOMO

The differences in how the two weighting schemes influence the importance scores become more pronounced as the number of models increases. As described in Algorithm 2, the formulas for their subset weights are \[\frac{1}{2^{n-1}-1} \quad\text{and}\quad \frac{1}{(n-1)\binom{n-1}{k}},\] where \(k\) is the size of each subset and \(n\) is the total number of models. The equal scheme (left formula) treats all subsets equally, so medium-sized subsets have considerable influence in the final result, as there are many such subsets. In contrast, the permutation-based scheme (right formula) adjusts the weights according to the subset size, giving the greatest weight to both the smallest and largest subsets while assigning small weights to the mid-sized subsets. Moreover, the weights assigned to the mid-sized subsets under the permutation-based approach decrease much faster with \(n\) than those under the equal weighting scheme (see Appendix for details). Consequently, when \(n\) is large,middle-sized subsets play a dominant role in determining the importance scores under the equal weighting scheme, whereas extreme-sized subsets primarily drive the scores under the permutation-based weighting approach.

Overall, the difference between the two weighting schemes is likely to arise mainly from the the extreme-sized subsets when \(n\) is large. This is because the weights given to the mid-sized subsets become increasingly similar, which are very small values on the order of \(10^{-3}\) even when \(n=8\), while the weights assigned to the smallest and largest subsets remain substantially different (Figure 3).

Figure 3: Comparison of weights assigned to a subset. The plot shows the weights assigned to a subset as the number of models \(n\) increases from 2 to 10. The red line represents the weights under the equal weighting scheme , while the blue and green lines represent the minimum and maximum weights, respectively, under the permutation-based weighting scheme. The minimum weight occurs when the subset size is around \((n-1)/2\), and the maximum weight occurs when the subset size is \(n-1\). As the number of models increases, the weights assigned by the two schemes become increasingly similar for mid-sized subsets whereas substantial differences remain for extreme-sized subsets.

model_importance( )

The model_importance() function calculates the importance scores of ensemble component models based on their contributions to improving ensemble prediction accuracy for each prediction task. It returns a single data frame of importance scores combined across all tasks. If a model missed predictions for a specific task, an NA value will be assigned for that task.

model_importance(
  forecast_data, oracle_output_data, ensemble_fun,
  importance_algorithm, subset_wt, min_log_score,
  ...
)

The forecast_data parameter is a data frame of predictions that is automatically coerced to a model_out_tbl format (if it is not one already), which is the standard S3 class model output format defined by the hubverse convention. If coercion fails, users may need to manually convert their data using the as_model_out_tbl() function from hubUtils.

The oracle_output_data is a data frame containing the actual observed values of the variables used to specify modeling targets. Further details on formatting and purpose are provided in Section 2.4.

The ensemble_fun argument specifies the ensemble method to use when evaluating model importance, which relies on implementations in the hubEnsembles package (Shandross et al. (2025)). The currently supported methods are "simple_ensemble" and "linear_pool". The "simple_ensemble" method returns the average of the predicted values from all component models per prediction task defined by task IDs, output_type, and output_type_id columns. The default aggregation function for this method is "mean", but it can be customized by specifying additional arguments through ..., such as agg_fun="median". When "linear_pool" is selected, the ensemble is a linear opinion pool of its component models. This method supports only an output_type of "mean", "quantile", or "pmf".

The importance_algorithm argument specifies the algorithm for model importance calculation, which can be either "lomo" (leave one model out) or "lasomo" (leave all subsets of models out). The subset_wt argument is employed only for the "lasomo" algorithm. This argument has two options: "equal" assigns equal weight to all subsets and "perm_based" assigns weight averaged over all possible permutations, as in the formula of Shapley values (Algorithm 2). The default values of importance_algorithm and subset_wt are "lomo" and "equal", respectively.

The min_log_score argument is relevant only for the output_type of "pmf", which uses Log Score as a scoring rule. It sets a minimum threshold for log scores to avoid issues with extremely low probabilities assigned to the true outcome, which can lead to undefined or negative infinite log scores. Any probability lower than this threshold will be adjusted to this minimum value before calculating the importance metric based on the log score. The default value is set to -10, following the CDC FluSight thresholding convention (Brooks et al. 2018; Reich et al. 2019). Users may choose a different value based on their practical needs.

The model_importance() function returns an object of S3 class model_imp_tbl of model important scores, with columns model_id, reference_date, output_type, and importance, along with any task ID columns (e.g., location, horizon, and target_end_date) present in the input forecast_data. Regardless of the original column name for the forecast generation date in the forecast_data, it is standardized to reference_date in the output. This standardization enables consistent handling of the forecast generation date across datasets, simplifying downstream processing. The importance column contains the calculated importance scores for each model and specific task, which are derived from the specified algorithm and ensemble method.

5.1 Print method

print() provides all computation results in a clear and organized manner. It displays model_id and importance by grouped task IDs (e.g., reference_date, location, and target_end_date) to facilitate easy interpretation of the importance scores for each model across different tasks.

5.2 Summary method

summary() provides a concise summary of the importance scores, including key statistics such as the number of models and tasks evaluated. It also displays top-scoring models for a subset of tasks. Additional summary details are available by specifying individual elements of the summary object as follows:

• .$all_tasks lists all the unique tasks evaluated, which are defined by the combinations of task ID columns present in the input forecast_data. This allows users to identify the scope of the evaluation.
• .$model_summary provides each model’s performance summary across tasks (e.g., the number of tasks each model submitted its forecast for, and the range of scores it achieved). This helps users understand the consistency and variability of each model’s importance across different tasks.
• .$task_winners identifies which model is the best for each task based on the highest importance score across the full set of tasks. From this information, users can quickly identify the most important model in each task.

5.3 Aggregate method

aggregate() allows users to obtain an overall importance score for each model by aggregating its importance scores across all evaluated tasks. This model_imp_tbl class-specific method provides arguments that users can specify to customize the aggregation process, including the handling of NA values and the choice of summary function (Table 4).

The by argument specifies the grouping variable(s) for summarization. Its default is "model_id", in which the aggregate() method summarizes importance scores for each model. Valid values for by include any combination of columns present in the importance_scores data frame.

The na_action argument allows for specifying how to handle NA values generated during importance score calculation for each task; these values occur when a model did not contribute to the ensemble prediction for a given task by missing its forecast submission. Three options are available: "worst", "average", and "drop". In each specific prediction task, if a model has any missing predictions, the "worst" option replaces the NA values with the smallest value among other models’ importance metrics, while the "average" option replaces them with the average of the other models’ importance metrics in that task. The "drop" option removes the NA values, which results in the exclusion of the model from the evaluation for that task.

The fun argument specifies a function used to summarize importance scores. An arithmetic mean (fun = mean) is the default, but other summary functions (e.g., fun = median) or user-defined functions may also be used. Additional arguments for the summary function fun can also be passed via ... if needed (e.g., fun = quantile, probs = 0.25 for a quartile summary).

The output returned by aggregate() is a data frame with columns model_id and importance_score_<fun>, where <fun> is the name of the summary function used (e.g., importance_score_mean when fun = mean). The output is sorted in descending order of the summarized importance scores, importance_score_<fun>.

6.1 Example data

Our example forecast data contains short-term predictions of weekly incident influenza hospitalizations in the US for Massachusetts (FIPS code 25) and Texas (FIPS code 48), generated on November 19, 2022. These forecasts are made for two target end dates, November 26, 2022 (horizon 1), and December 10, 2022 (horizon 3), and were produced by three models: ‘Flusight-baseline’, ‘MOBS-GLEAM_FLUH’, and ‘PSI-DICE’. The output type is median and the output_type_id column has NAs as no further specification is required for this output type. We have modified the example data slightly by removing some forecasts to demonstrate the handling of missing values. Therefore, MOBS-GLEAM_FLUH’s forecast for Massachusetts on November 26, 2022, and PSI-DICE’s forecast for Texas on December 10, 2022, are missing. We emphasize that this modification is made randomly and artificially to illustrate the impact of different handling approaches for missing forecasts and is not intended to provide a formal evaluation of model quality or performance.

forecast_data

The corresponding target data contains the observed hospitalization counts for these dates and locations.

target_data

When comparing the ground truth data and model predictions, we can see that forecasts for December 10, 2022 show larger deviations from the observed values compared to those for November 26, 2022. Thus, as we expect, prediction errors increase at longer horizons due to greater uncertainty. Additionally, the forecasts for Massachusetts are relatively more accurate compared to those for Texas, which tend to have higher errors (Figure 4).

Figure 4: Plot of three point forecasts (median) and the eventually observed values from the forecast_data and target_data for weekly incident influenza hospitalizations in Massachusetts (FIPS code 25) and Texas (FIPS code 48). Colored dots indicate the forecasts by three models, generated on November 19, 2022. Open black circles indicate the eventually observed values. MOBS-GLEAM_FLUH’s forecast for Massachusetts on November 26, 2022, and PSI-DICE’s forecast for Texas on December 10, 2022, are not shown. These forecasts were manually excluded from the original example dataset for demonstration purposes.

Evaluation using LOMO algorithm

We quantify the contribution of each model within the ensemble using the model_importance() function. The following code evaluates the importance of each ensemble member in the simple mean ensemble using the LOMO algorithm.

model_importance(
  forecast_data = forecast_data,
  oracle_output_data = target_data,
  ensemble_fun = "simple_ensemble",
  importance_algorithm = "lomo"
)

This call also generates informative messages that summarize the input data, including the number of dates on which forecasts were produced, the number of models and their ids, and whether the prediction tasks meet the minimum model requirement, alongside a note on how to use parallel processing to speed up the computation when there are many models and tasks. A print out of these messages is shown below.

Evaluating forecasts from 2022-11-19 to 2022-11-19  (a total of 1 forecast date(s)).

The available model IDs are:
     Flusight-baseline
     MOBS-GLEAM_FLUH
     PSI-DICE 
(a total of 3 models)

Note: This function uses 'furrr' and 'future' for parallelization.
To enable parallel execution, please set future::plan(multisession).

All tasks meet the minimum model requirement of 2 models.

The function output is a data frame containing model ids and their corresponding importance scores for each prediction task, along with task id columns. (Note that we rounded the importance scores and renamed the columns for better readability in the output below, but the original column names are retained in the actual output data frame.)

print(
  scores_lomo |>
    mutate(importance = round(importance, 2)) |>
    rename(
      ref_date = reference_date, h = horizon,
      loc = location, t_end_date = target_end_date,
      o_type = output_type, imp = importance
    )
)

#> Model importance result by task
#> ---------------------------------
#>             model_id   ref_date          target h loc t_end_date o_type     imp
#> 1  Flusight-baseline 2022-11-19 wk inc flu hosp 1  25 2022-11-26 median  -19.50
#> 2    MOBS-GLEAM_FLUH 2022-11-19 wk inc flu hosp 1  25 2022-11-26 median      NA
#> 3           PSI-DICE 2022-11-19 wk inc flu hosp 1  25 2022-11-26 median   19.50
#> 4  Flusight-baseline 2022-11-19 wk inc flu hosp 1  48 2022-11-26 median  -32.33
#> 5    MOBS-GLEAM_FLUH 2022-11-19 wk inc flu hosp 1  48 2022-11-26 median  -22.33
#> 6           PSI-DICE 2022-11-19 wk inc flu hosp 1  48 2022-11-26 median   54.67
#> 7  Flusight-baseline 2022-11-19 wk inc flu hosp 3  25 2022-12-10 median  -16.67
#> 8    MOBS-GLEAM_FLUH 2022-11-19 wk inc flu hosp 3  25 2022-12-10 median  -20.67
#> 9           PSI-DICE 2022-11-19 wk inc flu hosp 3  25 2022-12-10 median   37.33
#> 10 Flusight-baseline 2022-11-19 wk inc flu hosp 3  48 2022-12-10 median  182.00
#> 11   MOBS-GLEAM_FLUH 2022-11-19 wk inc flu hosp 3  48 2022-12-10 median -182.00
#> 12          PSI-DICE 2022-11-19 wk inc flu hosp 3  48 2022-12-10 median      NA

For models that missed forecasts for certain tasks, NA values were assigned in the importance column for those tasks.

Calling summary() shows that three models were used and four tasks were evaluated, along with a preview of the top-performing model for each task.

summary(scores_lomo)

#> === Summary of importance scores by task ===
#> Number of models: 3 
#> Number of tasks: 4 
#> 
#> === Top scoring model by task for a subset of tasks ======================================== 
#>           target horizon location target_end_date top_model importance
#>  wk inc flu hosp       1       25      2022-11-26  PSI-DICE      19.50
#>  wk inc flu hosp       1       48      2022-11-26  PSI-DICE      54.67
#>  wk inc flu hosp       3       25      2022-12-10  PSI-DICE      37.33
#> --------------------------------------------
#> * More details are available in the summary object (e.g., $all_tasks, $model_summary, $task_winners).

As indicated in the output, more details about the summary are available through the summary object’s elements as follows.

s <- summary(scores_lomo)
s$all_tasks

#>            target horizon location target_end_date
#> 1 wk inc flu hosp       1       25      2022-11-26
#> 2 wk inc flu hosp       1       48      2022-11-26
#> 3 wk inc flu hosp       3       25      2022-12-10
#> 4 wk inc flu hosp       3       48      2022-12-10

Each row represents a unique combination of task IDs, from which we verify that four different tasks were evaluated.

s$model_summary

#>            model_id n_tasks min_importance max_importance n_NA
#> 1 Flusight-baseline       4         -32.33         182.00    0
#> 2   MOBS-GLEAM_FLUH       4        -182.00         -20.67    1
#> 3          PSI-DICE       4          19.50          54.67    1

We observe that ‘Flusight-baseline’ submitted forecasts for all four tasks (n_NA = 0), while ‘MOBS-GLEAM_FLUH’ and ‘PSI-DICE’ submitted forecasts for only three tasks due to one missing forecast (n_NA = 1). Each model’s importance scores vary across tasks. ‘Flusight-baseline’ shows the largest range of scores that includes a negative minimum value and a positive maximum value, while ‘MOBS-GLEAM_FLUH’ and ‘PSI-DICE’ have scores that are all positive or all negative across the three tasks they submitted forecasts for.

s$task_winners

#>            target horizon location target_end_date         top_model max_score
#> 1 wk inc flu hosp       1       25      2022-11-26          PSI-DICE     19.50
#> 2 wk inc flu hosp       1       48      2022-11-26          PSI-DICE     54.67
#> 3 wk inc flu hosp       3       25      2022-12-10          PSI-DICE     37.33
#> 4 wk inc flu hosp       3       48      2022-12-10 Flusight-baseline    182.00

Models with the highest importance scores for each task are identified in the top_model column with their importance score in the max_score column. ‘PSI-DICE’ is the best model for three out of the four tasks, while ‘Flusight-baseline’ is the best for the remaining task.

Visualization can be performed using functions from ggplot2. The following example shows a bar plot of importance scores across models and tasks, with panels faceted by combinations of task ID values.

ggplot(scores_lomo, aes(x = model_id, y = importance, fill = model_id)) +
  geom_col() +
  coord_flip() +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.25) +
  facet_grid(
    cols = vars(target, horizon, location, target_end_date),
    scales = "free_x"
  ) +
  labs(
    x = "Model ID", y = "Importance Score",
    title = "Model Importance by Task"
  ) +
  scale_x_discrete(labels = function(x) gsub("[-_]", "-\n", x)) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    panel.spacing.x = unit(0.5, "lines"), legend.position = "none"
  )

We aggregate the importance scores for each model by averaging across all tasks. NA values are removed during the averaging process by setting the na_action argument to "drop".

aggregate(scores_lomo, by = "model_id", na_action = "drop", fun = mean)

#> Overall model importance across tasks
#> ---------------------------------------- 
#> # A tibble: 3 × 2
#>   model_id          importance_score_mean
#>   <chr>                             <dbl>
#> 1 PSI-DICE                           37.2
#> 2 Flusight-baseline                  28.4
#> 3 MOBS-GLEAM_FLUH                   -75

The results show that, overall, the model ‘PSI-DICE’ has the highest importance score, followed by ‘Flusight-baseline’ and ‘MOBS-GLEAM_FLUH’. That is, ‘PSI-DICE’ contributes the most to improving the ensemble’s predictive performance, whereas ‘MOBS-GLEAM_FLUH’, which has a negative score, detracts from the ensemble’s performance. The low importance score of ‘MOBS-GLEAM_FLUH’ is mainly due to a substantially larger prediction error for Texas on the target end date of December 10, 2022, compared to other models, while its missing forecast for Massachusetts for November 26, 2022, was not factored into the evaluation. This single large error significantly affected its contribution score.

Another approach to handling missing values is to use the "worst" option for na_action, which replaces missing values with the worst (i.e., minimum) score among the other models for the same task.

aggregate(scores_lomo, by = "model_id", na_action = "worst", fun = mean)

#> Overall model importance across tasks
#> ---------------------------------------- 
#> # A tibble: 3 × 2
#>   model_id          importance_score_mean
#>   <chr>                             <dbl>
#> 1 Flusight-baseline                  28.4
#> 2 PSI-DICE                          -17.6
#> 3 MOBS-GLEAM_FLUH                   -61.1

The results show that the importance score of ‘Flusight-baseline’ is unchanged because it has no missing forecasts. We observe that the importance score of ‘PSI-DICE’, which was previously positive, has now decreased to a negative value when compared to the evaluation using the "drop" option for na_action. Moreover, ‘MOBS-GLEAM_FLUH’ still ranks the lowest, but the importance score has increased. This change is related to the varying forecast accuracy across different tasks. For the target end date of November 26, 2022, in Massachusetts, most forecasts are relatively accurate. Thus, even if the ‘MOBS-GLEAM_FLUH’ is assigned the worst value of importance score for its missing forecast, including this value in the averaging is not detrimental to the overall importance metric; rather, it is more beneficial than excluding it. In contrast, for the target end date of December 10, 2022, in Texas, the forecasts have much larger errors across the board, and assigning the worst value of importance score to the missing forecast of ‘PSI-DICE’ in this task has a detrimental effect on averaging importance scores. This is because the scale of the importance scores is influenced by the magnitude of the prediction errors: for tasks with small errors, the scores remain moderate, while tasks with large errors can yield importance scores of much greater magnitude.

It is also possible to impute the missing scores with intermediate values by assigning the average importance scores of other models in the same task. This strategy may offer a more balanced trade-off by mitigating the influence of the missing data without overly penalizing or overlooking them.

aggregate(scores_lomo, by = "model_id", na_action = "average", fun = mean)

#> Overall model importance across tasks
#> ---------------------------------------- 
#> # A tibble: 3 × 2
#>   model_id          importance_score_mean
#>   <chr>                             <dbl>
#> 1 Flusight-baseline                  28.4
#> 2 PSI-DICE                           27.9
#> 3 MOBS-GLEAM_FLUH                   -56.2

Evaluation using LASOMO algorithm

We now demonstrate the use of the LASOMO algorithm for evaluating model importance. Since we explored the difference of na_action options in the previous LOMO example Section 6.2, we focus on options for subset_wt, which specifies how weights are assigned to subsets of models when calculating importance scores, with na_action fixed to "drop".

The following code and corresponding outputs illustrate the evaluation using each weighting scheme.

scores_lasomo_eq <- model_importance(
  forecast_data = forecast_data,
  oracle_output_data = target_data,
  ensemble_fun = "simple_ensemble",
  importance_algorithm = "lasomo",
  subset_wt = "equal"
)
aggregate(scores_lasomo_eq, by = "model_id", na_action = "drop", fun = mean)

#> Overall model importance across tasks
#> ---------------------------------------- 
#> # A tibble: 3 × 2
#>   model_id          importance_score_mean
#>   <chr>                             <dbl>
#> 1 PSI-DICE                           47.4
#> 2 Flusight-baseline                  24.3
#> 3 MOBS-GLEAM_FLUH                   -79.8

scores_lasomo_perm <- model_importance(
  forecast_data = forecast_data,
  oracle_output_data = target_data,
  ensemble_fun = "simple_ensemble",
  importance_algorithm = "lasomo",
  subset_wt = "perm_based"
)
aggregate(scores_lasomo_perm, by = "model_id", na_action = "drop", fun = mean)

#> Overall model importance across tasks
#> ---------------------------------------- 
#> # A tibble: 3 × 2
#>   model_id          importance_score_mean
#>   <chr>                             <dbl>
#> 1 PSI-DICE                           44.8
#> 2 Flusight-baseline                  25.3
#> 3 MOBS-GLEAM_FLUH                   -78.6

In this example, there are only three models (\(n = 3\)), and the weights do not differ significantly between the two weighting schemes. Therefore, the resulting outputs show little difference. However, in general, with a larger number of models, the two weighting schemes may yield quite different importance scores for each model, as discussed in Section 3.1.

In this section, we explored each component model’s contribution to the ensemble accuracy with only three models. An extensive application in more complex scenarios with larger ensembles is presented in our companion methodological paper (Kim et al. 2026).

It should be noted that the example presented here is designed for illustration purposes to demonstrate the use of the proposed software, and the results should not be interpreted as an authoritative evaluation of model performance. A comprehensive analysis in our prior methodological study (Kim et al. 2026) shows that many models, including the ‘MOBS-GLEAM_FLUH’ model, are assessed as more important than the baseline model.