Worked example: combining three flu hospitalisation forecasts

This page runs every combination operation and weight estimator in the package. The forecasts are a real hubverse slice bundled with ForecastEnsembles.jl: three models from the example-complex-forecast-hub, each predicting weekly flu hospitalisations on 2022-12-17 at horizon 1 across five US locations (national plus CA, FL, NY, TX), at the standard 23 quantile levels.

Every block below is executed when the documentation is built, so what you read is what the current code does.

julia

using ForecastEnsembles, CSV, DataFrames, Distributions, Random

flu = CSV.read(joinpath(pkgdir(ForecastEnsembles), "data", "flu_forecasts.csv"),
    DataFrame; types = Dict(:output_type_id => Float64, :location => String))

ft = ForecastTable(flu; task_id_cols = [:reference_date, :target_end_date,
    :horizon, :location, :target])

ForecastTable(
  models       = InlineStrings.String31["Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE"]
  output_type  = quantile
  task_id_cols = [:reference_date, :target_end_date, :horizon, :location, :target]
  rows         = 345
)

Equal-weight quantile mean (Vincentization)

The simplest combination: at each (location, τ) take the unweighted mean of the three model quantile values.

julia

first(DataFrame(combine(ft, QuantileEnsemble(:mean))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14550.7
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15901.3
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17320.7
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18325.0
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	19261.7

Or the median ensemble the COVID-19 hub used as its default:

julia

first(DataFrame(combine(ft, QuantileEnsemble(:median))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13738.0
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14791.0
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15934.0
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17477.0
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18676.0

Mixture (linear opinion pool)

Average the distributions rather than the quantile values: F = Σᵢ wᵢ Fᵢ. Each model's quantiles are reconstructed into a continuous distribution (PCHIP interior, Normal tails), and the mixture CDF is then inverted at each requested level by bisection. On quantile input this path is deterministic, so n_samples does not apply to it; that field governs the :sample path only.

julia

first(DataFrame(combine(ft, MixtureEnsemble())), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13732.6
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14633.9
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15480.2
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	16558.8
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17435.7

This gives a different answer from Vincentization in general: averaging quantile values is not the same operation as averaging CDFs. On this slice the largest gap between the two is:

julia

vinc = DataFrame(combine(ft, QuantileEnsemble(:mean))).value
mix = DataFrame(combine(ft, MixtureEnsemble())).value
maximum(abs.(vinc .- mix))

6513.631889712095

Geometric (logarithmic) pool

Multiply the member densities instead of averaging them — a product of experts, sharper than the linear pool where the models agree:

julia

first(DataFrame(combine(ft, LogarithmicPool())), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13683.2
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14756.2
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15795.7
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17124.4
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18227.2

Robust mean (trimmed / winsorised)

Drop or clamp the most extreme model at each (location, τ) before averaging, for robustness to an outlier submission:

julia

first(DataFrame(combine(ft, TrimmedMean(; fraction = 0.2))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13738.0
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14791.0
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15934.0
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17477.0
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18676.0

With mode = :winsorise the extremes are clamped to the surviving range rather than dropped:

julia

first(DataFrame(combine(ft, TrimmedMean(; fraction = 0.2, mode = :winsorise))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13738.0
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14791.0
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15934.0
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17477.0
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18676.0

Those two tables are identical here, and both match the median ensemble above. fraction trims round(fraction · n) models from each end, capped so at least one value survives, so with three models nothing is trimmed until fraction rises above about 0.17, and once it does only one value is left: trimming leaves the median, and winsorising clamps both extremes onto it. The two modes diverge only where at least three values survive, so four models agree at every fraction too; five or more, trimmed lightly, is where they part.

Hand-supplied weights

Either combination operation takes an EnsembleWeights:

julia

w = EnsembleWeights(DataFrame(
    model_id = ["Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE"],
    weight = [0.2, 0.4, 0.4]
))

first(DataFrame(combine(ft, QuantileEnsemble(:mean; weights = w))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14147.6
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15304.4
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	16512.6
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17574.6
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18481.2

julia

first(DataFrame(combine(ft, MixtureEnsemble(; weights = w))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13619.4
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14468.0
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15265.6
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	16271.8
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17028.6

A history to learn weights from

The estimators below need past forecasts with matching observations. The bundled slice is a single date, so this page builds a small synthetic history for the same three models. Real use would supply past hub submissions here.

Sample-typed history, for the score-driven estimators:

julia

const MODELS = ["Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE"]

rng = MersenneTwister(20221217)
T, K = 12, 200

# A latent signal the models track, with the observation landing near but not on
# it. Forecasting the realised value exactly would leave nothing for the
# estimators, or for the recalibration further down, to work on.
signal = 100.0 .+ 20.0 .* randn(rng, T)
train_obs = DataFrame(t = 1:T, observed = signal .+ 6.0 .* randn(rng, T))

rows = DataFrame[]
for (mid, sd) in zip(MODELS, (35.0, 15.0, 22.0)), t in 1:T
    push!(rows, DataFrame(model_id = mid, output_type = "sample",
        output_type_id = 1:K, t = t,
        value = signal[t] .+ sd .* randn(rng, K)))
end
train_ft = ForecastTable(reduce(vcat, rows); task_id_cols = [:t])

ForecastTable(
  models       = ["Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE"]
  output_type  = sample
  task_id_cols = [:t]
  rows         = 7200
)

MOBS-GLEAM_FLUH is the sharpest member of this synthetic history, so the estimators below should favour it.

Quantile-typed history, for QRA and BLP, which take quantile input:

julia

levels = sort(unique(flu.output_type_id))

qrows = DataFrame[]
for (mid, sd) in zip(MODELS, (35.0, 15.0, 22.0)), t in 1:T
    push!(qrows, DataFrame(model_id = mid, output_type = "quantile",
        output_type_id = levels, t = t,
        value = signal[t] .+ sd .* quantile.(Normal(), levels)))
end
qtrain_ft = ForecastTable(reduce(vcat, qrows); task_id_cols = [:t])

ForecastTable(
  models       = ["Flusight-baseline", "MOBS-GLEAM_FLUH", "PSI-DICE"]
  output_type  = quantile
  task_id_cols = [:t]
  rows         = 828
)

Weights from CRPS-stacking

fit(CRPSStacking(), ...) optimises a simplex weight vector against CRPS on sample forecasts. The result plugs straight into either combination operation:

julia

stacked = fit(CRPSStacking(), train_ft, train_obs)
DataFrame(weights(stacked))

3×2 DataFrame

Row	model_id	weight
	String	Float64
1	Flusight-baseline	3.7579e-5
2	MOBS-GLEAM_FLUH	0.999899
3	PSI-DICE	6.38643e-5

julia

first(DataFrame(combine(ft, MixtureEnsemble(; weights = stacked))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13738.0
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14790.9
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15933.8
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17476.7
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18675.7

Weights from QRA

QRA fits a quantile regression of past observations on past per-model forecasts. Two configurations matter:

Joint (per_quantile_weights = false): one weight vector across all τ, usable by either operation.
Per-τ (per_quantile_weights = true): a different weight vector at each τ, usable by QuantileEnsemble.

julia

qra = fit(
    QRA(; per_quantile_weights = true, enforce_normalisation = true,
        intercept = false),
    qtrain_ft, train_obs
)
first(DataFrame(combine(ft, QuantileEnsemble(:mean; weights = qra))), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	13738.0
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	14791.0
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	15934.0
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	17477.0
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	18676.0

Some configurations do not reduce to a weight vector at all — fits with an intercept, unconstrained fits, or fits spanning several task groups. For those weights(::FittedQRA) returns nothing:

julia

loose = fit(QRA(; enforce_normalisation = false), qtrain_ft, train_obs)
weights(loose) === nothing

true

Passing such a fit where weights are expected says which of the three cases it is, and how to refit:

julia

try
    MixtureEnsemble(; weights = loose)
catch err
    println(err.msg)
end

FittedQRA does not expose ensemble weights: this QRA fit is unconstrained, so its coefficients are a regression, not weights summing to one; refit with `enforce_normalisation = true` for a weight vector.

You can still apply such a fit directly with combine(ft, loose), which uses the regression coefficients to predict quantiles. Nothing constrains those coefficients, so the predicted quantiles can cross — on this small history they do, and the package warns about it. Fitting with per_quantile_weights = true and noncross = true constrains the ordering at the training points only; on a new forecast table the predictions can still cross, so sort each task's values before submission.

Score-driven estimators

These take any scoring function you supply, with the signature score(samples, y; w) returning a scalar. ScoringRules.jl is the natural companion; the package itself depends on no scoring library, so this page defines a small weighted CRPS to keep the example self-contained and to show the contract explicitly:

julia

function crps(dat::AbstractVector, y::Real; w = nothing)
    n = length(dat)
    ww = w === nothing ? fill(1.0 / n, n) : w ./ sum(w)
    ex = sum(ww[i] * abs(dat[i] - y) for i in 1:n)
    ee = sum(ww[i] * ww[j] * abs(dat[i] - dat[j]) for i in 1:n, j in 1:n)
    return ex - 0.5 * ee
end

crps (generic function with 1 method)

Generic stacking against that score:

julia

DataFrame(weights(fit(Stacking(crps), train_ft, train_obs)))

3×2 DataFrame

Row	model_id	weight
	String	Float64
1	Flusight-baseline	1.17977e-9
2	MOBS-GLEAM_FLUH	1.0
3	PSI-DICE	2.41417e-9

Performance-based weighting, scoring each member independently with no optimiser:

julia

DataFrame(weights(fit(InverseScore(crps), train_ft, train_obs)))

3×2 DataFrame

Row	model_id	weight
	String	Float64
1	Flusight-baseline	0.0102816
2	MOBS-GLEAM_FLUH	0.840436
3	PSI-DICE	0.149282

Adaptive weighting, which walks the time column and updates after each round, so a model forecasting badly this week counts for less next week:

julia

hedged = fit(Hedge(crps; time_col = :t), train_ft, train_obs)
DataFrame(weights(hedged))

3×2 DataFrame

Row	model_id	weight
	String	Float64
1	Flusight-baseline	0.00714136
2	MOBS-GLEAM_FLUH	0.861509
3	PSI-DICE	0.13135

Its trajectory records the weights after every update, and weight_stability summarises how much each model's weight moved:

julia

weight_stability(hedged)

3×2 DataFrame

Row	model_id	total_variation
	String	Float64
1	Flusight-baseline	0.255761
2	MOBS-GLEAM_FLUH	0.481874
3	PSI-DICE	0.238742

Hierarchical weighting fits a vector per stratum, shrunk toward a shared one, for when locations behave differently but not independently. It needs a stratifying column, so this builds a two-location history in which a different member is the sharp one in each:

julia

prows = DataFrame[]
orows = DataFrame[]
for loc in ("CA", "TX")
    sds = loc == "CA" ? (12.0, 30.0, 25.0) : (30.0, 12.0, 25.0)
    for t in 1:T
        push!(orows, DataFrame(location = loc, t = t,
            observed = signal[t] + 6.0 * randn(rng)))
        for (mid, sd) in zip(MODELS, sds)
            push!(prows, DataFrame(model_id = mid, output_type = "sample",
                output_type_id = 1:K, location = loc, t = t,
                value = signal[t] .+ sd .* randn(rng, K)))
        end
    end
end
pool_ft = ForecastTable(reduce(vcat, prows); task_id_cols = [:t, :location])
pool_obs = reduce(vcat, orows)

pooled = fit(PartialPooling(crps; strata = [:location]), pool_ft, pool_obs)
pooled.weights

6×3 DataFrame

Row	model_id	weight	location
	String	Float64	String
1	Flusight-baseline	0.332238	TX
2	MOBS-GLEAM_FLUH	0.667762	TX
3	PSI-DICE	9.44878e-11	TX
4	Flusight-baseline	0.634054	CA
5	MOBS-GLEAM_FLUH	0.365946	CA
6	PSI-DICE	9.663e-11	CA

Each location's weights favour the member that is sharp there, which a single global vector could not express.

Training on a trailing window, and comparing schemes

Windowed restricts any estimator to the most recent times, and backtest compares schemes out of sample by expanding the training window one step at a time.

backtest scores each fold with a function you supply, of the shape (forecast, observations) -> score. Built from the crps above:

julia

function fold_crps(fc, obs)
    d = DataFrames.innerjoin(DataFrame(fc), obs; on = :t)
    per = DataFrames.combine(DataFrames.groupby(d, :t),
        [:value, :observed] =>
            ((v, y) -> crps(Float64.(v), Float64(first(y)))) => :s)
    return sum(per.s) / nrow(per)
end

rolling = Windowed(CRPSStacking(), 6; time_col = :t)

backtest(train_ft, train_obs,
    ["expanding" => CRPSStacking(), "rolling" => rolling];
    time_col = :t, min_train = 4, score_fn = fold_crps)

16×3 DataFrame

Row	scheme	t	score
	String	Int64	Float64
1	expanding	5	3.50134
2	rolling	5	3.54268
3	expanding	6	3.30746
4	rolling	6	3.29721
5	expanding	7	3.61556
6	rolling	7	3.62154
7	expanding	8	3.63011
8	rolling	8	3.65514
9	expanding	9	4.38324
10	rolling	9	4.40804
11	expanding	10	4.51709
12	rolling	10	4.32567
13	expanding	11	3.46529
14	rolling	11	3.40258
15	expanding	12	4.61817
16	rolling	12	4.37887

Recalibrated mixture (beta-transformed linear pool)

BLP corrects a miscalibrated linear pool. It fits a Beta to the pool's PIT values on quantile-typed history, then reweights the pool's quantile levels rather than the models. Its usual motivation is a pool whose tails are too narrow, but the transform works in either direction, and this synthetic history is the other case: the members carry far more spread than the observations warrant, so the fitted Beta is peaked at 0.5 and the recalibrated interval comes out much narrower than the raw pool's.

julia

blp = fit(BLP(), qtrain_ft, train_obs)
first(DataFrame(combine(ft, blp)), 5)

5×9 DataFrame

Row	model_id	output_type	output_type_id	reference_date	target_end_date	horizon	location	target	value
	String	Symbol	Float64	Date	Date	Int64	String	String15	Float64
1	hub-ensemble	quantile	0.01	2022-12-17	2022-12-24	1	US	wk inc flu hosp	21444.7
2	hub-ensemble	quantile	0.025	2022-12-17	2022-12-24	1	US	wk inc flu hosp	21891.7
3	hub-ensemble	quantile	0.05	2022-12-17	2022-12-24	1	US	wk inc flu hosp	22307.7
4	hub-ensemble	quantile	0.1	2022-12-17	2022-12-24	1	US	wk inc flu hosp	22834.5
5	hub-ensemble	quantile	0.15	2022-12-17	2022-12-24	1	US	wk inc flu hosp	23136.0

Because it recalibrates the pooled distribution rather than estimating per-model weights, weights(::FittedBLP) is nothing.

What's where in the data

ForecastEnsembles extends DataFrames.combine, so one combine covers both the ensemble combinations above and this group-and-count:

julia

DataFrames.combine(
    DataFrames.groupby(flu, [:model_id, :location]),
    nrow => :n_quantiles
)

15×3 DataFrame

Row	model_id	location	n_quantiles
	String31	String	Int64
1	Flusight-baseline	US	23
2	Flusight-baseline	06	23
3	Flusight-baseline	12	23
4	Flusight-baseline	36	23
5	Flusight-baseline	48	23
6	MOBS-GLEAM_FLUH	US	23
7	MOBS-GLEAM_FLUH	06	23
8	MOBS-GLEAM_FLUH	12	23
9	MOBS-GLEAM_FLUH	36	23
10	MOBS-GLEAM_FLUH	48	23
11	PSI-DICE	US	23
12	PSI-DICE	06	23
13	PSI-DICE	12	23
14	PSI-DICE	36	23
15	PSI-DICE	48	23

Three models × five locations × 23 quantile levels: 345 rows.

Worked example: combining three flu hospitalisation forecasts ​

Equal-weight quantile mean (Vincentization) ​

Mixture (linear opinion pool) ​

Geometric (logarithmic) pool ​

Robust mean (trimmed / winsorised) ​

Hand-supplied weights ​

A history to learn weights from ​

Weights from CRPS-stacking ​

Weights from QRA ​

Score-driven estimators ​

Training on a trailing window, and comparing schemes ​

Recalibrated mixture (beta-transformed linear pool) ​

What's where in the data ​

Worked example: combining three flu hospitalisation forecasts

Equal-weight quantile mean (Vincentization)

Mixture (linear opinion pool)

Geometric (logarithmic) pool

Robust mean (trimmed / winsorised)

Hand-supplied weights

A history to learn weights from

Weights from CRPS-stacking

Weights from QRA

Score-driven estimators

Training on a trailing window, and comparing schemes

Recalibrated mixture (beta-transformed linear pool)

What's where in the data