The lopensemble package provides an easy way to combine
predictions from individual time series or panel data models to an
ensemble. lopensemble stacks Models according to the
Continuous Ranked Probability Score (CRPS) over k-step ahead
predictions. It is therefore especially suited for timeseries and panel
data. A function for leave-one-out CRPS may be added in the future.
Predictions need to be predictive distributions represented by
predictive samples. Usually, these will be sets of posterior predictive
simulation draws generated by an MCMC algorithm.
Given some training data with true observed values as well as
predictive samples generated from different models,
crps_weights finds the optimal (in the sense of minimizing
expected cross-validation predictive error) weights to form an ensemble
from these models. Using these weights,
mixture_from_samples can then provide samples from the
optimal model mixture by drawing from the predictice samples of the
individual models in the correct proportion. This gives a mixture model
solely based on predictive samples and is in this regard superior to
other ensembling techniques like Bayesian Model Averaging.
Here is an example of how to use lopensemble with a toy
data set. The data contains true observed values as well as predictive
samples generated by different models.
library("data.table")
#>
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#>
#> %notin%
splitdate <- as.Date("2020-03-28")
print(example_data)
#> geography model sample_id date predicted observed
#> <char> <char> <num> <Date> <num> <num>
#> 1: Tatooine ARIMA 1 2020-03-14 1.719445 1.655068
#> 2: Tatooine ARIMA 2 2020-03-14 1.896555 1.655068
#> 3: Tatooine ARIMA 3 2020-03-14 1.766821 1.655068
#> 4: Tatooine ARIMA 4 2020-03-14 1.714007 1.655068
#> 5: Tatooine ARIMA 5 2020-03-14 1.726421 1.655068
#> ---
#> 103996: Coruscant Naive 496 2020-04-08 1.460421 1.543976
#> 103997: Coruscant Naive 497 2020-04-08 1.057846 1.543976
#> 103998: Coruscant Naive 498 2020-04-08 1.433936 1.543976
#> 103999: Coruscant Naive 499 2020-04-08 1.719357 1.543976
#> 104000: Coruscant Naive 500 2020-04-08 0.781818 1.543976Obtain weights based on trainin data, create mixture based on predictive samples in the testing data.
weights <- lopensemble::crps_weights(traindata)
test_mixture <- lopensemble::mixture_from_samples(testdata, weights = weights)
print(test_mixture)
#> Key: <geography, sample_id, date, observed>
#> geography sample_id date observed predicted model
#> <char> <num> <Date> <num> <num> <char>
#> 1: Coruscant 1 2020-03-29 1.5002089 1.2133070 CRPS_Mixture
#> 2: Coruscant 1 2020-03-30 1.4713551 2.3203675 CRPS_Mixture
#> 3: Coruscant 1 2020-03-31 1.4154616 1.2099242 CRPS_Mixture
#> 4: Coruscant 1 2020-04-01 1.3426560 2.3038299 CRPS_Mixture
#> 5: Coruscant 1 2020-04-02 1.3298251 1.7740069 CRPS_Mixture
#> ---
#> 10996: Tatooine 500 2020-04-04 0.6750499 0.3437742 CRPS_Mixture
#> 10997: Tatooine 500 2020-04-05 0.7225369 0.2328917 CRPS_Mixture
#> 10998: Tatooine 500 2020-04-06 0.7532836 0.3716259 CRPS_Mixture
#> 10999: Tatooine 500 2020-04-07 0.7598158 0.4109292 CRPS_Mixture
#> 11000: Tatooine 500 2020-04-08 0.7719879 0.5695866 CRPS_MixtureScore predictions using the scoringutils package
(install from CRAN using
install.packages("scoringutils")).
# combine data.frame with mixture with predictions from other models
score_df <- data.table::rbindlist(list(testdata, test_mixture), fill = TRUE)
# CRPS score for all predictions using scoringutils
score_df[, crps := scoringutils::crps(unique(observed), t(predicted)),
by = .(geography, model, date)
]
# summarise scores
score_df[, mean(crps), by = model][, data.table::setnames(.SD, "V1", "CRPS")]Display other metrics
Given a cumulative distribution function (CDF) with a finite first moment of a probabilistic forecast, the corresponding Continuous Ranked Probability Score can be written as \[crps(F,y)=\mathbb{E}_X|X-y|- \frac{1}{2}\mathbb{E}_{X,X^\prime}|X-X^\prime|.\tag{1}\]
The CRPS can be used to stack different models to obtain one ensemble mixture model. Let us assume we have data from \(T\) time points \(t = 1, \dots, T\) in \(R\) regions \(1, \dots, R\). Observations are denoted \(y_{tr}\). Predictive samples are generated from \(K\) different models \(k, \dots, K\). For every observation \(y_{tr}\) the \(S\) predictive samples are denoted \(x_{1ktr}, \dots, x_{Sktr}\).
Using these predictive draws, we can compute the CRPS of the \(k\)-th model for the observation \(y_{tr}\) at time \(t\) in region \(r\) as
\[ \widehat {crps}_{ktr}= \widehat {crps}_(x_{1ktr}, \dots, x_{Sktr},y_{tr})= \frac{1}{S} \sum_{s=1}^S |x_{sktr}-y_{tr}| - \frac{1}{2S^2} \sum_{s, j=1}^S |x_{sktr}- x_{jktr}|. \tag{2}\]
Now we want to aggregate predictions from these \(K\) models. When the prediction is a mixture of the \(K\) models with weights \(w_1, \dots, w_s\), the CRPS can be expressed as
\[ \widehat {crps}_{agg, tr} (w_1, \dots, w_K) = \frac{1}{S} \sum_{k=1}^K w_k \sum_{s=1}^S |x_{skt}-y_t| - \frac{1}{2S^2} (\sum_{k=1}^K \sum_{k, k'=1 }^K w_k w_{k'} \sum_{s, j=1}^S |x_{skt}- x_{jk't}| ). \tag{3}\]
This expression is quadratic in \(w\). We only need to compute \(\sum_{s=1}^S |x_{skt}-y_{tr}|\), \(\sum_{s, j=1}^S |x_{sktr}- x_{jktr}|\), and \(\sum_{s, j=1}^S |x_{sktr}- x_{jk'tr}|\) for all \(k, k'\) pairs once and store them for all weight values in the optimization.
The CRPS for the mixture of all models for all observations can simply obtained by summing up the individual results from equation 3 over all regions and time points. However, we can also weight predictions from different time points and regions differently. This makes sense for example if we have very little data in the beginning or if current older observations are less characteristic of the current and future dynamics. In this case we might want to downweight the early phase in the final model evaluation. Similarly, we might want to give more or less weight to certain regions. Mathematically we can introduce a time-varying weight \(\lambda_1, \dots, \lambda_T\), e.g. \(\lambda_t = 2-(1-t/T)^2\) to penalize earler estimates. Likewise we can introduce a region-specific weight \(\tau_r\).
To obtain the CRPS weights we finally solve a quadratic optimization:
\[\begin{align*} &\min_{w_1, \dots, w_K} \sum_{t=1}^T \sum_{r=1}^R\lambda_t\tau_r \widehat {crps}_{agg, tr} (w), \\ &s.t. ~{0\leq w_1, \dots, w_K \leq 1, \sum_{k=1}^K w_k=1}. \end{align*}\]