--- title: "`modelimportance`: Evaluating model importance within a multi-model ensemble in R" author: "Minsu Kim, Li Shandross, Evan L. Ray, and Nicholas G. Reich" bibliography: references.bib output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{`modelimportance`: Evaluating model importance within a multi-model ensemble in R} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4, warning = FALSE, message = FALSE ) oldpar <- par(no.readonly = TRUE) oldopts <- options() library(dplyr) library(tidyr) library(ggplot2) library(kableExtra) library(modelimportance) ``` ## Abstract Ensemble forecasts are commonly used to support decision-making and policy planning across various fields because they often offer improved accuracy and stability compared to individual models. As each model has its own unique characteristics, understanding and measuring the value of each constituent model can support the construction of effective ensembles. The [R]{.proglang} package `modelimportance` provides tools to quantify how each component model contributes to the accuracy of ensemble performance for both point and probabilistic forecasts. The package supports multiple ensemble methods and multiple model importance metrics. Additionally, the software offers customizable options for handling missing values. These features enable the package to serve as a versatile tool for researchers and practitioners. It helps not only in constructing an effective ensemble model across a wide range of forecasting tasks, but also in understanding the role of each model within the ensemble and gaining insights into individual models themselves. This package follows the ['hubverse'](https://hubverse.io/) framework, which is a collection of open-source software, tools and data standards developed to promote collaborative modeling hub efforts and simplify their setup and operation. Doing so enables seamless integration and flexibility with other forecasting tools and systems, allowing many analyses to be performed on existing hubs. ## 1. Introduction {#sec:intro} Ensemble forecasting is a method to produce a single, consolidated prediction by combining forecasts generated from different models. While individual models' strengths and weaknesses are pronounced within ensembles, the combination of many models tend to offset each other and create an ensemble forecast that is more robust and accurate than any single component model [@gneiting2005weather; @hastie01statisticallearning]. Specifically, ensembles effectively mitigate the bias and variance arising from the predictions of individual models by averaging them out, and aggregating in this way can reduce prediction errors and improve overall performance. Enhanced prediction accuracy and robustness enable the achievement of more reliable predictions, thereby improving decision-making. For this reason, ensemble forecasting is widely used across various domains such as weather forecasting [@Guerra_2020; @gneiting2005weather], financial modeling [@SUN2020101160; @math11041054], and infectious disease outbreak forecasting [@ray_prediction_2018; @reich_accuracy_2019; @lutz_applying_2019; @viboud_rapidd_2018] to name a few. For example, throughout the COVID-19 pandemic, the US COVID-19 Forecast Hub collected individual models developed by over 90 different research groups and built a probabilistic ensemble forecasting model for COVID-19 cases, hospitalizations, and deaths in the US, based on those models’ predictions. The ensemble model served as the official short-term forecasts for the US Centers for Disease Control and Prevention (CDC) [@cramer2022united]. The quality of forecasts is assessed by evaluating their errors, biases, sharpness, and/or calibration using different scoring metrics. The selection of scoring metrics depends on the type of forecast: point forecasts (e.g., mean, median) and probabilistic forecasts (e.g., quantiles, probability mass function). Commonly used assessment tools for point forecasts are the mean absolute error (MAE) and the mean squared error (MSE), which calculate the average magnitude of forecast errors. Scoring metrics for probabilistic forecasts consider the uncertainty and variability in predictions and provide concise evaluations through numerical scores [@gneiting_strictly_2007]. Some examples include the weighted interval score (WIS), the continuous ranked probability score (CRPS), and the log score [@bracher_evaluating_2021]. (Note that CRPS is a general scoring rule that can be computed either analytically in closed form or numerically from samples, and WIS is a quantile-based approximation of CRPS.) Several [R]{.proglang} packages have been developed to evaluate forecast quality across both point and probabilistic settings. To name a few, the `fable` package [@Rpackage-fable] is widely used for univariate time series forecasting and includes functions for accuracy measurement. The `Metrics` [@Rpackage-Metrics] and `MLmetrics` [@Rpackage-MLmetrics] provide a wide range of performance metrics specifically designed for evaluating machine learning models. The `scoringRules` [@Rpackage-scoringRules] package offers a comprehensive set of proper scoring rules for evaluating probabilistic forecasts and supports both univariate and multivariate settings. The `scoringutils` [@bosse2022evaluating] package offers additional features to the functionality provided by `scoringRules`, which makes it more useful for certain tasks, such as summarizing, comparing, and visualizing forecast performance. These packages have been valuable to evaluate individual models as independent entities, using performance metrics selected for each specific situation or problem type. However, they do not measure the individual models' contributions to the enhanced predictive accuracy when used as part of an ensemble. Building on our prior methodological study [@kim2026beyond], we emphasize that strong standalone model performance does not automatically imply a large positive contribution once the model is used within an ensemble. The `modelimportance` package operationalizes this idea in software, providing tools to evaluate the role of each model as an ensemble member within an ensemble model, rather than focusing on the individual predictive performance per se. In ensemble forecasting, certain models contribute more significantly to the overall predictions than others. Assessing the impact of each component model on ensemble predictions is methodologically similar to determining variable importance in traditional regression and machine learning models, where variable importance measures evaluate how much individual variables improve the accuracy of the model's predictive performance or reduce the average loss. [R]{.proglang} packages such as `randomForest` [@Rpackage-randomForest], `caret` [@Rpackage-caret], `xgboost` [@Rpackage-xgboost], and `gbm` [@Rpackage-gbm] implement these functions for different types of models: random forest models, general machine learning models, extreme gradient boosting models, and generalized boosted regression models, respectively. These packages focus on feature-level importance within a single model and do not measure the contribution of individual models within an ensemble. In contrast, the tools in `modelimportance` quantify how each component model helps enhance the ensemble model's predictive performance. They assign numerical scores to each model using a forecast accuracy metric selected based on the forecast type. Our methods are based on the concept of Shapley values in cooperative game theory, which measure a player’s average contribution to the game's overall outcome (@Shapley1953). There are several approaches that utilize Shapley values in black-box machine learning models to understand how each feature affects the model’s predictive power (@lundberg2017unified; @lundberg2020local; @covert2020understanding). Related work by @lipiecki2024postprocessing used Shapley-based attribution to study contributions of component models in a combined ensemble after postprocessing point forecasts into probabilistic forecasts. That methodology was later implemented in `PostForecasts.jl` [@Jlpackage-lipiecki2025postforecast]. `PostForecasts.jl` is conceptually aligned with our work in that both approaches leverage Shapley values to evaluate component forecasters. However, while `PostForecasts.jl` was developed for [Julia]{.proglang} users originally focusing on energy economics, our [R]{.proglang} package focuses on epidemiological forecasts and offers a broader range of capabilities, such as handling missing forecasts, compatibility with the hubverse forecasting ecosystem, and rich tools for summarizing and visualizing importance scores across tasks. These capabilities are particularly useful for hub organizers, who oversee collaborative forecasting systems that combine submissions from many teams into a single ensemble forecast [@shandross2026multi]. Examples include efforts coordinated by the US CDC and the European Centre for Disease Prevention and Control. By quantifying each component model's contribution, `modelimportance` can support evidence-based decisions about ensemble design and maintenance. The package follows conventions defined by the 'hubverse', a community-maintained ecosystem of open software and data standards for collaborative forecasting hubs [@hubverse_docs; @kerr2025coordinating]. Using these shared data structures allows our package to be seamlessly integrated with existing hub pipelines and deployed across multiple active hubs. We note that there are 31 hubs as of March 2026, including model development and training hubs. We highlight some development practices we employed, such as unit testing of individual functions, object-oriented programming using S3 classes, continuous integration testing on different operating systems, and independent code review by peer developers. This emphasis on quality control is a key strength of this work. The paper proceeds as follows. [Section 2](#sec:data) describes (a) how the `modelimportance` package relates to the hubverse framework, including its dependencies, (b) the model output formats defined within hubverse, and (c) the structure of data presentation for both forecasts and actual observations. [Section 3](#sec:algorithms) presents two algorithms implemented in `modelimportance` for calculating the model importance metric: leave-one-model-out and leave-all-subsets-of-models-out. [Section 4](#sec:main-function) demonstrates the various functionalities `modelimportance` supports. [Section 5](#sec:s3-infrastructure) describes the S3 class structure and related methods implemented in the package, followed by examples of their usage in [Section 6](#sec:examples). [Section 7](#sec:computational-complexity) discusses computational complexity and strategies for efficient computation. [Section 8](#sec:implementation-and-availability) highlights our quality assurance measures and the code's availability, and then we close this paper with a summary and a discussion of possible extensions. ## 2. Data {#sec:data} ### 2.1 Dependencies and related software {#subsec:dependences} The `modelimportance` package is designed to work with the hubverse framework and, accordingly, depends on several packages in the hubverse ecosystem, such as `hubUtils` (@Rpackage-hubUtils), `hubEnsembles` (@Rpackage-hubEnsembles), and `hubEvals` (@Rpackage-hubEvals). `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](#subsec:model_output_format)). 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](#sec:examples)) from the `hubExamples` (@Rpackage-hubExamples) 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` [@Rpackage-ggplot2] functions can be applied. Parallel computing is supported by the `furrr` [@Rpackage-furrr] and `future` [@Rpackage-future] 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 {#subsec: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 [@shandross2026multi]. 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. ```{r, results='asis'} #| echo: false #| label: tbl-example-model_output #| tbl-cap: "Example of the model output for incident influenza #| hospitalizations (top 10 rows) extracted from `forecast_data_example` data in the #| `modelimportance` package, which is originally from `hubExamples` package." forecast_data_example |> dplyr::filter( .data$reference_date == "2022-12-17", .data$horizon != 0 ) |> head(10) |> knitr::kable( format = "html", caption = "Table 1: Example of the model output for incident influenza hospitalizations (top 10 rows) extracted from `forecast_data_example` bundled in the `modelimportance` package, which is originally from `hubExamples` package.") |> kableExtra::kable_styling( font_size = 12, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = TRUE ) ``` 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. ```{r, fig.pos="t"} #| echo: false #| label: fig-example-model_output #| fig.cap: > #| 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. forecast_data_example |> filter( reference_date == "2022-12-17", horizon != 0, output_type_id %in% c(0.05, 0.5, 0.95) ) |> pivot_wider(names_from = output_type_id, values_from = value) |> rename(lower = "0.05", upper = "0.95", value = "0.5") |> ggplot(aes(x = target_end_date)) + facet_grid(~model_id) + geom_point(aes(y = value, color = "medians"), size = 2) + geom_line(aes(y = value, color = "medians"), linewidth = 1) + geom_ribbon( aes( ymin = lower, ymax = upper, fill = "#3388FF" ), alpha = 0.5 ) + geom_point( data = target_data_example |> filter(target_end_date <= "2022-12-17"), aes(y = observation, group = 1, color = "obs") ) + geom_line( data = target_data_example |> filter(target_end_date <= "2022-12-17"), aes(y = observation, group = 1, color = "obs") ) + geom_point( data = target_data_example |> filter(target_end_date > "2022-12-17"), aes(y = observation, group = 1, color = "truth"), shape = 1, alpha = 1 ) + geom_line( data = target_data_example |> filter(target_end_date > "2022-12-17"), aes(y = observation, group = 1, color = "truth"), alpha = 0.75 ) + # coord_cartesian(ylim = c(0, 500)) + scale_x_date(breaks = target_data_example$target_end_date, date_labels = "%Y-%m-%d") + labs( y = "Weekly Hospitalization", x = "Date" ) + scale_color_manual( name = "", values = c( "medians" = "DodgerBlue", "obs" = "Black", "truth" = "Black" ), labels = c( "Forecast (Predictive median)", "Observed data before forecasting", "Eventually observed value" ) ) + scale_fill_manual("", values = "#3388FF", labels = "90% Prediction interval" ) + theme( axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 8), axis.text.x = element_text(size = 6, angle = 90), axis.text.y = element_text(size = 8), strip.text = element_text(size = 8), legend.title = element_blank(), legend.text = element_text(size = 8), legend.position = "bottom", legend.direction = "vertical", legend.box = "horizontal" ) ``` ### 2.3 Forecast data representation {#subsec:model_output} 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]{.pkg} 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]{.pkg} 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. ```{r} #| echo: false #| label: tbl-pair-output-scoringrule #| tbl-cap: "Pairs of output types and their associated scoring rules for evaluating prediction performance." data.frame( "Output Type" = c("mean", "median", "quantile", "pmf"), "Scoring Rule" = c("SE", "AE", "WIS", "Log Score"), Description = c( "Squared error (SE): the squared difference between the predicted value and the observed value", "Absolute error (AE): the absolute difference between the predicted value and the observed value", paste0( "Weighted interval score (WIS): a quantile-based approximation of the continuous ranked", " probability score (CRPS) for evaluating quantile forecasts" ), "Logarithm of the probability assigned to the true outcome (LogScore)" ), check.names = FALSE ) |> knitr::kable( format = "html", caption = paste0( "Table 2: Pairs of output types and their associated scoring rules ", "for evaluating prediction performance." ) ) |> kableExtra::kable_styling( font_size = 13, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = TRUE ) ``` ### 2.4 Oracle output data {#subsec: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 [@hubverse_docs; @kerr2025coordinating]. 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. Method description and algorithms {#sec:algorithms} This section provides a brief description of the leave one model out (LOMO) and leave all subsets of models out (LASOMO) algorithms, which are used to compute the model importance score. The basic idea of measuring the importance of each component model is to evaluate the change in ensemble performance when that model is included or excluded in the ensemble construction. More specifically, we compare the performance of an ensemble with and without a specific model for a specific task, and consider the difference in performance as the importance of that model for that task (Figure 2). We apply this idea across many tasks and then average the task-level values to obtain model-level summaries. The full theoretical development and additional experiments are provided in @kim2026beyond. ```{r} #| echo: false #| label: fig-concept #| fig-cap: > #| Figure 2: Conceptual illustration of measuring model importance in a three-model setting. #| Each circle represents a forecasting model with component model(s) shown inside. #| Stars indicate the performance score measured by a positively oriented scoring rule (e.g., $-$WIS), #| with more stars indicating better accuracy. LOMO computes a single performance difference #| by removing the target model 1 from the full ensemble. LASOMO leverages multiple performance #| differences across all subsets containing the model 1 and aggregates them. #| out.width: "80%" #| fig-align: "center" knitr::include_graphics("figure-concept.png") ``` LOMO involves creating an ensemble by excluding one component model from the entire set of models. Let ${A}$ be a set of $n$ models and $F^i$ be a forecast produced by model $i$, where $i = 1,2, \dots, n.$ Each ensemble excludes exactly one model while including all the others. $F^{A^{-i}}$ denotes the ensemble forecast constructed using all forecasts $F^A$ except $F^i$. Model $i$'s importance score using LOMO is calculated as the difference in accuracy, as measured by a specific scoring rule, between $F^{A^{-i}}$ and $F^A$ (Algorithm 1). For example, when evaluating model 1 within an ensemble of three models ($n=3$), LOMO creates an ensemble forecast $F^{\{2,3\}}$ using only $F^2$ and $F^3$. The performance of this reduced ensemble is then compared to the full ensemble forecast $F^{\{1,2,3\}}$, which incorporates all three models. We note that a model can make an ensemble better or worse, and thus the importance score for model 1 can be positive or negative accordingly. ```{r} #| echo: false #| label: algorithm-lomo #| out.width: '90%' #| fig.cap: "Algorithm 1: Leave one model out (LOMO)" knitr::include_graphics("algorithm-lomo.jpg") ``` On the other hand, LASOMO involves ensemble constructions from all possible subsets of models. For each subset $S$ that does not contain the model $i$, $S \cup \{i\}$ plays the role of ${A}$ in the LOMO; the score associated with the subset $S$ is the difference of measures between $F^S$ and $F^{S \cup \{i\}}$. Then, all scores are aggregated across all possible subsets that the model $i$ does not belong to (Algorithm 2). For example, using the earlier setup of three forecast models, LASOMO considers three subsets, which we denote by $S_1=\{2\}$, $S_2=\{3\}$, and $S_3=\{2, 3\}$, to calculate the importance score of model 1 (excluding all subsets that include model 1). The ensemble forecasts $F^{\{2\}}, F^{\{3\}}$, and $F^{\{2,3\}}$ are then compared to $F^{\{1,2\}}, F^{\{1,3\}}$, and $F^{\{1,2,3\}}$, respectively. The performance differences attributable to model 1's inclusion are aggregated, which results in the importance score for model 1. We note that the subsets (e.g., $S_1, S_2,$ and $S_3$) may have different weights during the aggregation process. ```{r} #| echo: false #| label: algorithm-lasomo #| out.width: '90%' #| fig.cap: "Algorithm 2: Leave all subsets of models out (LASOMO)" knitr::include_graphics("algorithm-lasomo.jpg") ``` The `modelimportance` package offers two weighting options for subsets: one assigns equal (uniform) weights to all subsets, and the other assigns weights based on their size, similar to the definition of Shapley values, where a player’s average contribution is aggregated over all possible coalitions (or, equivalently, over all permutations of players) (@Shapley1953). Users can choose one to evaluate the contribution of each model in a manner suited to their preferred framework. Uniform weighting may be preferred when each subset size is of equal analytical interest, regardless of how many such subsets there are, while size-based weighting may be preferred when users want to maintain the original Shapley value interpretation, preventing numerous medium-sized subsets from dominating the importance scores. A detailed discussion of the differences between the two weighting schemes follows in the next section. ### 3.1 Comparison of weighting schemes in LASOMO {#sec:comparison-weighting-schemes} 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](#sec: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). ```{r, fig.pos="t"} #| echo: false #| label: fig-lasomo-weights #| fig.cap: > #| 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. data.frame(n = 2:10) |> mutate( w_eq = 1 / (2^(n - 1) - 1), w_perm_min = 1 / ((n - 1) * choose(n - 1, floor((n - 1) / 2))), w_perm_max = 1 / (n - 1) ) |> ggplot(aes(x = n)) + geom_line(aes(y = w_eq, color = "eq"), size = 0.75, linetype = "longdash") + geom_point(aes(y = w_eq, color = "eq"), size = 2, shape = 15) + geom_line(aes(y = w_perm_min, color = "perm_min"), size = 0.75) + geom_point(aes(y = w_perm_min, color = "perm_min"), size = 2, shape = 16) + geom_line(aes(y = w_perm_max, color = "perm_max"), size = 0.75) + geom_point(aes(y = w_perm_max, color = "perm_max"), size = 2, shape = 17) + scale_color_manual( name = "Weights", values = c( "eq" = "#F8766D", "perm_min" = "#619CFF", "perm_max" = "#00BA38" ), labels = c( eq = expression(w^{ eq }), perm_min = expression(w^{ perm - min }), perm_max = expression(w^{ perm - max }) ) ) + scale_x_continuous(breaks = 2:10) + labs( x = "Number of models (n)", y = "Weight assigned to a subset", color = "Weighting Scheme" ) + theme( axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 8), axis.text.y = element_text(size = 8), legend.title = element_text(size = 9), legend.position = "bottom", legend.text = element_text(size = 9), legend.spacing.x = unit(0.5, "mm"), legend.key.width = unit(1, "cm") ) ``` ## 4. Evaluating models with the model_importance() function{#sec:main-function} In this section, we describe the usage of the function `model_importance()`, where multiple options are available to customize the evaluation framework (Table 3). 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. The function requires a minimum of two models per task, as with only one model, the importance score is not defined. The output of the function is 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( ) 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. ```{r code-model_importance} #| eval: false 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](#subsec:oracle_output_data). The `ensemble_fun` argument specifies the ensemble method to use when evaluating model importance, which relies on implementations in the `hubEnsembles` package (@Rpackage-hubEnsembles). 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 [@brooks2018nonmechanistic; @reich_accuracy_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. ```{r} #| echo: false #| label: tbl-arguments1 #| tbl-cap: "Table 3: Description of the arguments for the #| `model_importance()` function, including their purpose, possible values, #| and default settings." tbl_argument1 <- data.frame( "Argument" = c( "`forecast_data`", "`oracle_output_data`", "`ensemble_fun`", "`importance_algorithm`", "`subset_wt`", "`min_log_score`", "`...`" ), "Description" = c( "Forecasts", "Ground truth data", "Ensemble method", "Algorithm to calculate importance", "Method for assigning weight to subsets when using LASOMO algorithm", "Minimum value to replace for log score", "Optional arguments for `'simple_ensemble'`" ), "Possible Values" = c( "Table of model output", "Table of oracle output", "`'simple_ensemble'`, `'linear_pool'`", "`'lomo', 'lasomo'`", "`'equal', 'perm_based'`", "Non-positive numeric", "Varies" ), Default = c( "N/A", "N/A", "`'simple_ensemble'`", "`'lomo'`", "`'equal'`", -10, "`agg_fun='mean'`" ), check.names = FALSE ) tbl_argument1 |> knitr::kable(format = "html", escape = FALSE, caption = "Table 3: Description of the arguments for the `model_importance()` function, including their purpose, possible values, and default settings.") |> kableExtra::kable_styling( font_size = 13, bootstrap_options = c("striped", "condensed", "responsive"), full_width = FALSE ) ``` ## 5. S3 class and methods {#sec:s3-infrastructure} We have defined a custom S3 class, `model_imp_tbl` that represents the output of the `model_importance()` function, which extends the base `data.frame` class. Objects of this class enable dispatch for class-specific methods for printing, summarizing, and aggregating them across tasks. ### 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). ```{r} #| echo: false #| label: tbl-arguments2 #| tbl-cap: > #| Description of the arguments for the `aggregate()` function for `model_imp_tbl` objects, #| including their purpose, possible values, and default settings." tbl_argument2 <- data.frame( "Argument" = c( "`importance_scores`", "`by`", "`na_action`", "`fun`", "`...`" ), "Description" = c( "Model importance scores produced by `model_importance()`", "Grouping variable(s) for summarization", "Method to handle `NA` values", "Function to summarize importance scores", "Optional arguments for `\"fun\"`" ), "Possible Values" = c( "data frame", "grouping variable(s)", "`\"drop\", \"worst\", \"average\"`", "summary function", "depends on `fun`" ), Default = c( "N/A", "`\"model_id\"`", "`\"drop\"`", "`mean`", "N/A" ), check.names = FALSE ) tbl_argument2 |> knitr::kable( format = "html", escape = FALSE, caption = paste0( "Table 4: Description of the arguments for the `aggregate()` function for ", "`model_imp_tbl` objects, including their purpose, possible values, and default settings." ) ) |> kableExtra::kable_styling( font_size = 13, bootstrap_options = c("striped", "condensed", "responsive"), full_width = FALSE ) ``` 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_`, where `` 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_`. ## 6. Examples {#sec:examples} In this section, we illustrate how to use the `model_importance()` function to evaluate the importance of component models within an ensemble. The examples show various combinations of the arguments described in [Section 4](#sec:main-function). Example forecast and target data are originally sourced from `hubExamples` package, which provides sample datasets for multiple modeling hubs in the hubverse format. ### 6.1 Example data {#sec: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 `NA`s 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. ```{r code-forecast_data} #| eval: false forecast_data ``` ```{r forecast_data, tidy = FALSE} #| echo: false forecast_data <- readRDS("vignette-example-forecast_data.rds") forecast_data |> knitr::kable( format = "html" ) |> kableExtra::kable_styling( font_size = 12, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = TRUE ) ``` The corresponding target data contains the observed hospitalization counts for these dates and locations. ```{r code-target_data} #| eval: false target_data ``` ```{r target_data} #| echo: false target_data <- readRDS("vignette-example-oracle_data.rds") target_data |> knitr::kable( format = "html" ) |> kableExtra::kable_styling( font_size = 12, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = TRUE ) ``` 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). ```{r, fig.pos="t"} #| echo: false #| label: fig-example-median-lomo #| fig.cap: > #| 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. target_data |> mutate(model_id = "Observed") |> rename(value = oracle_value) |> rbind( forecast_data |> select("target_end_date", "target", "location", "value", "model_id") ) |> ggplot(aes(x = target_end_date)) + geom_point(aes(y = value, color = model_id, shape = model_id), size = 2) + geom_point( data = target_data, aes(y = oracle_value, color = "Observed", shape = "Observed"), size = 3 ) + facet_wrap(~location, scales = "free_y", labeller = labeller(location = function(x) paste0("Location: ", x)) ) + scale_x_date( breaks = target_data$target_end_date, date_labels = "%Y-%m-%d", expand = expansion(add = c(5, 5)) ) + scale_color_manual( name = "model_id", values = c( "Flusight-baseline" = "#F8766D", "MOBS-GLEAM_FLUH" = "#00BA38", "PSI-DICE" = "#619CFF", "Observed" = "black" ), limits = c( "Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE", "Observed" ) ) + scale_shape_manual( name = "model_id", values = c( "Flusight-baseline" = 16, "MOBS-GLEAM_FLUH" = 17, "PSI-DICE" = 15, "Observed" = 18 ), limits = c( "Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE", "Observed" ) ) + labs( y = "Weekly Hospitalization", x = "Date" ) + theme( axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9), strip.text = element_text(size = 9), legend.title = element_blank(), legend.text = element_text(size = 8), legend.position = "bottom", legend.box = "horizontal", legend.spacing.x = unit(0.25, "mm"), legend.key.width = unit(0.4, "cm") ) ``` ### Evaluation using LOMO algorithm {#sec:example-lomo} 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. ```{r code-model_importance-lomo} #| eval: false 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. ```r 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.) ```{r} #| label: scores_lomo #| echo: false scores_lomo <- model_importance( forecast_data = forecast_data, oracle_output_data = target_data, ensemble_fun = "simple_ensemble", importance_algorithm = "lomo" ) ``` ```{r} #| eval: false 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 ) ) ``` ```{r} #| label: print_scores_lomo #| echo: false 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 )) ``` 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. ```{r} #| eval: false summary(scores_lomo) ``` ```{r} #| label: summary_scores_lomo #| echo: false summary(scores_lomo) ``` As indicated in the output, more details about the summary are available through the summary object's elements as follows. ```{r} #| eval: false s <- summary(scores_lomo) s$all_tasks ``` ```{r} #| label: summary_scores_lomo_all_tasks #| echo: false s <- summary(scores_lomo) s$all_tasks ``` Each row represents a unique combination of task IDs, from which we verify that four different tasks were evaluated. ```{r} #| eval: false s$model_summary ``` ```{r} #| label: summary_scores_lomo_model_summary #| echo: false s$model_summary ``` 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. ```{r} #| eval: false s$task_winners ``` ```{r} #| label: summary_scores_lomo_task_winners #| echo: false s$task_winners ``` 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. ```{r} #| eval: false 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" ) ``` ```{r} #| label: plot_scores_lomo #| echo: false #| fig-align: "center" #| fig-width: 6 #| fig-height: 4 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"`. ```{r} #| eval: false aggregate(scores_lomo, by = "model_id", na_action = "drop", fun = mean) ``` ```{r} #| label: aggregate_scores_lomo_drop #| echo: false aggregate(scores_lomo, by = "model_id", na_action = "drop", fun = mean) ``` 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. ```{r} #| eval: false aggregate(scores_lomo, by = "model_id", na_action = "worst", fun = mean) ``` ```{r} #| label: aggregate_scores_lomo_worst #| echo: false aggregate(scores_lomo, by = "model_id", na_action = "worst", fun = mean) ``` 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. ```{r} #| eval: false aggregate(scores_lomo, by = "model_id", na_action = "average", fun = mean) ``` ```{r} #| label: aggregate_scores_lomo_avg #| echo: false aggregate(scores_lomo, by = "model_id", na_action = "average", fun = mean) ``` ### Evaluation using LASOMO algorithm {#sec:example-lasomo} 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](#sec:example-lomo), 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. ```{r code-model_importance-summary-lasomo-equal} #| eval: false 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) ``` ```{r} #| label: lasomo-NAdrop-equal #| echo: false scores_lasomo_eq <- suppressMessages(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) ``` ```{r code-model_importance-summary-lasomo-perm} #| eval: false 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) ``` ```{r} #| label: lasomo-NAdrop-perm #| echo: false scores_lasomo_perm <- suppressMessages(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) ``` \ \ 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](#sec:comparison-weighting-schemes). 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 [@kim2026beyond]. 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 [@kim2026beyond] shows that many models, including the 'MOBS-GLEAM_FLUH' model, are assessed as more important than the baseline model. ## 7. Computational complexity {#sec:computational-complexity} Content coming soon. This section will describe computational complexity of LOMO and LASOMO algorithms implemented in the package, depending the numbers of models and tasks. This section describes the computational complexity of LOMO and LASOMO algorithms implemented in the package, depending on the numbers of ensemble component models and prediction tasks. We conducted a computational experiment using simulated data with a point forecast ('median'). The execution time was measured for both algorithms across varying numbers of models and tasks. We also compared the execution times between sequential and multisession computing environments to demonstrate the efficiency of parallelization in handling computationally intensive tasks. We performed all experiments on a machine running macOS 15 with an Apple M4 chip and 16 GB RAM under R version 4.4.3. Parallel computations were implemented using the `future` (`multisession` backend) with four worker processes. The execution times were recorded in seconds. Figure 5 illustrates the computational runtime for four different numbers of prediction tasks (10, 20, 50, and 100) and varying numbers of models (ranging from 2 to 10) for both the LOMO and LASOMO algorithms under sequential and parallel execution modes. The scale of the $y$-axis approximately ranges from 0 to 12 seconds for LOMO and from 0 to 1110 seconds for LASOMO. Overall, the time needed to compute model importance scores increases with the number of models and tasks for both algorithms. Specifically, when the number of tasks is 10 or 20, the difference in computation time between sequential and parallel executions is not significant. However, as the number of prediction tasks rises to 50, the difference becomes more noticeable. The speedup from parallel execution relative to sequential execution is approximately 4-fold for LOMO in these settings of 10, 20, 50, and 100 prediction tasks, which is consistent with the use of four worker processes. Similarly, for LASOMO, a speedup of approximately 3-fold is observed across these settings, which is less than the 4-fold speedup observed with LOMO due to the heavier computational intensity of LASOMO. Nevertheless, the efficiency of parallel execution is more pronounced for LASOMO than for LOMO, particularly when eight or more models are involved in the ensemble construction and the number of prediction tasks is large. ```{r} #| echo: false #| label: fig-feasibility-comparison #| fig.cap: > #| Figure 5: Runtime to compute model importance scores for point predictions by number of ensemble #| component models (ranging from 2 to 10) and number of prediction tasks (10, 20, 50, and 100) #| under sequential and parallel execution modes (multisession backend, with four workers) #| for LOMO and LASOMO algorithms. Each point represents the elapsed time (in seconds) #| for a given number of models, while lines (solid/dashed) show the overall trend #| fitted to the observations. The scale of the $y$-axis is different between the two algorithms. #| out.width: "100%" #| fig-align: "center" knitr::include_graphics("runtime_comparison-fitted.png") ``` Theoretically, the computational complexity of LOMO and LASOMO algorithms is $O(t \cdot n)$ and $O(t\cdot 2^n)$, respectively, where $n$ is the number of models and $t$ is the number of prediction tasks. However, the observations from our computational experiment in the sequential execution mode exhibited slight discrepancies from these theoretical complexities, under a limited setting involving between 2 and 10 models and 10, 20, 50, and 100 prediction tasks. While the observed runtimes showed linear growth with respect to the number of tasks for both algorithms, consistent with the theoretical complexities, the growth with respect to the number of models is sublinear for LOMO and slower than the expected exponential for LASOMO. This deviation from the theoretical patterns is likely due to the practical efficiency of the implementations, including shared computation. Specifically, once the ensemble predictions are computed, the importance scores are derived from these predictions without needing to recompute the ensemble for each model. It should be noted that the scaling increment and ratio between adjacent values of $n$ in LOMO and LASOMO, respectively, is gradually increasing, which implies the empirical complexity is approaching the theoretical complexity with the number of models. However, the execution time scales much more rapidly with the number of models than the theoretical complexities, notably in LASOMO, as larger sub-ensembles require more expensive evaluation and scoring operations across all prediction tasks as well as larger intermediate objects. That is, the theoretical complexity does not factor in subset size, which does, in fact, impact the computational speed in practice. This phenomenon is also observed even in the parallel execution setting (Figure 6b). We further explored computational feasibility, focusing on the parallel execution mode (Figure 6). We observe that, in LOMO, the elapsed time increases linearly with the number of models and the rate of increase depends on the number of tasks, with a steeper increase as the number of tasks grows. In contrast, in LASOMO, the elapsed time increases much more rapidly as the ensemble size grows. In particular, when the number of models exceeds 12, it takes more than 20 minutes to compute the importance scores, even with fixed 5 prediction tasks. It is because the number of subsets of models that need to be evaluated in LASOMO grows exponentially with the number of models, leading to a rapidly rising computational intensity. This observation highlights the trade-off between the comprehensiveness of the LASOMO algorithm and its computational challenges in practice, particularly with a large number of models. ```{r} #| echo: false #| label: fig-feasibility-parallel-all_in_one #| fig-cap: > #| Figure 6: Runtime to compute model importance scores for point predictions by number of ensemble #| component models and prediction tasks during parallel execution via multisession backend. #| Each dot represents the elapsed time (in seconds) for a given number of models. #| For LOMO, the number of models evaluated varies from 2 to 100 for each 5, 20, 50, and #| 100 prediction tasks. The fitted lines illustrate the trend of the observations. #| For LASOMO, the number of models evaluated varies from 2 to 14 while fixing 5 prediction tasks #| due to exponentially growing computational intensity with the number of models. #| The exact elapsed time is shown for LASOMO, which captures the steep increase in runtime #| as the number of models grows, even before the high jump occurs with models over 12. #| out.width: "100%" #| fig-align: "center" knitr::include_graphics("runtime_parallel-fitted.png") ``` ## 8. Implementation and availability {#sec:implementation-and-availability} The `modelimportance` package is implemented in [R]{.proglang} and distributed via CRAN and GitHub, under the MIT license. We conducted unit tests using the `testthat` package [@Rpackage-testthat] to verify that all inputs and outputs are properly formatted and ensure that all functions work correctly as expected, including those used internally. We also performed continuous integration testing using GitHub Actions to maintain functionality across platforms, including Windows, macOS, and Linux. Integrated GitHub Action, we employed the `lintr` package [@Rpackage-lintr] to maintain code quality and detect potential issues, and the Air formatter [@air_formatter] to ensure consistent code style across the codebase. All code changes were systematically reviewed by fellow team members, and this enhanced reliability. ## Summary and discussion {#sec:discussion} Multi-model ensemble forecasts often provide better accuracy and robustness than single models, and are widely used in decision-making and policy planning across various domains. The contribution of each component model to the accuracy of the ensemble depends on its own unique characteristics. The `modelimportance` package enables the quantification of the value that each component model adds to the ensemble performance in different evaluation contexts. In the example analysis, we showed the package workflow using a simple example dataset specifically designed solely to demonstrate our package's capabilities; thus the results should not be taken as a definitive assessment of the models. A more complete model importance analysis is provided in our companion methodological paper [@kim2026beyond]. The primary `model_importance()` function returns an S3 object of class `model_imp_tbl` of tabular data containing component models and their importance metrics. Users are given a choice of ensemble methods, model importance algorithms (either LOMO or LASOMO), and options to handle missing values. These features enable the package to serve as a versatile tool to aid collaborative efforts to construct an effective ensemble model across a wide range of forecasting tasks. We note that unit testing with continuous integration ensures the reliability of all functions and the overall quality of code across multiple platforms. The [modelimportance]{.pkg} package still has several areas in which it may be enhanced. Namely, although the package currently supports four different output types ('mean', 'median', 'quantile, and 'pmf'), other output types are widely used in practice. For example, the 'sample' output type is commonly used in the US Flu Scenario Modeling Hub [@FluSMH]. This format includes multiple simulated values (samples) from the forecast distribution. Support for the 'sample' output type is under consideration for future releases, and, in general, extensions to support more output types would aim to broaden the scope of applications in real-world forecasting tasks. ## Acknowledgements {.unnumbered} We acknowledge Zhian N. Kamvar for debugging and resolving coding issues while developing the package. We are also grateful to Matthew Cornell for his advice on unit testing, which greatly helped us improve the structure and testing our code with a solid understanding of unit testing practices. Moreover, we would like to thank the hubverse development team for their data standards, on which our package is based. ## Appendix {.unnumbered} ### Weights for subsets in LASOMO {#sec:appendix} In the LASOMO algorithm, two weighting schemes are available for subsets of models in the calculation of model importance scores: equal weights and permutation-based weights. Let $n$ be the total number of models and $k$ be the size of a subset that does not include the model being evaluated. The formulas for the weights under each scheme are as follows: \begin{align*} w^{\text{eq}} &= \frac{1}{2^{n-1}-1}, \\ w^{\text{perm}} &= \frac{1}{(n-1)\binom{n-1}{k}}, \end{align*} where the superscripts "eq" and "perm" denote the equal and permutation-based weighting schemes, respectively. The maximum weight under the permutation-based scheme occurs when $k=n-1$, which is ${1}/{(n-1)}$. The minimum weight occurs when the subset size is around ${(n-1)}/{2}$ (i.e., $k=\lfloor (n-1)/2 \rfloor$), which is approximately $\displaystyle\frac{\sqrt{\pi(n-1)/2}}{(n-1)2^{n-1}}$ by Stirling's approximation. Given a fixed mid-sized subset, as $n$ increases, the weight assigned to this subset under the equal weighting scheme decreases at a rate of $O({1}/{2^n})$, while under the permutation-based scheme, it decreases at a much faster rate of $O({1}/({\sqrt{n}\,2^n}))$. This indicates that as the number of models grows, that mid-sized subset becomes significantly less influential in determining model importance scores when using the permutation-based weighting scheme compared to the equal weighting scheme. On the other hand, for subsets of extreme sizes (e.g., $k=1$ or $k=n-1$), the weights under permutation-based weighting scheme decrease only at $O({1}/{n})$, much slower under the equal weighting scheme. This implies that in scenarios with a large number of models, the contributions of these extreme-sized subsets play a relatively larger role in the calculation of model importance scores when using permutation-based weights compared to the equal weighting approach. ```{r, include = FALSE} options(oldopts) par(oldpar) ``` ## References