Does Ensembling Improve Calibration? A Monte Carlo Study of Probability Calibration in Ensemble Classifiers

Abstract

Ensemble methods are well established as a means to improve discriminative performance in classification tasks. Whether ensembling similarly improves the calibration of predicted probabilities—the agreement between predicted confidence and observed frequency—remains poorly understood. We conduct a Monte Carlo simulation study with 500 replicates, each comprising 546 synthetic patients, to evaluate five ensemble aggregation strategies (simple average, trimmed mean, median, geometric mean, and selective averaging) applied to five base models with distinct calibration profiles (well-calibrated, overconfident, underconfident, shifted-high, shifted-low). We measure Expected Calibration Error (ECE), Maximum Calibration Error (MCE), Brier score, and Brier decomposition into reliability, resolution, and uncertainty components. Our central finding challenges the conventional assumption: simple averaging—the most commonly used ensemble aggregation—degrades calibration relative to a well-calibrated individual model, with a Cohen's d effect size of +2.89 (higher ECE) and a win rate of only 0.6% across 500 trials. The geometric mean is the only aggregation strategy that consistently improves calibration (ECE = 0.048 vs. 0.054 for the well-calibrated individual, win rate 67.2%, d = −0.44). The trimmed mean and median occupy intermediate positions, improving upon simple averaging but rarely beating a well-calibrated individual. Post-hoc recalibration via Platt scaling substantially reduces the calibration damage caused by simple averaging (ECE from 0.112 to 0.063) but does not eliminate it entirely. These findings have practical implications for clinical decision support systems, risk stratification tools, and any application where predicted probabilities inform downstream actions: ensembling for accuracy and ensembling for calibration are fundamentally different objectives, and the aggregation strategy must be chosen with calibration explicitly in mind.

1. Introduction

Probability calibration—the property that a predicted probability of p corresponds to a true event frequency of p—is a foundational requirement for any predictive model whose outputs inform decisions under uncertainty. In clinical medicine, a model predicting a 30% probability of sepsis onset must be correct roughly 30% of the time for clinicians to appropriately allocate monitoring resources. In financial risk modeling, a credit default probability of 5% must correspond to an actual default rate near 5% for portfolio risk calculations to be valid. In weather forecasting, the discipline of calibration has been central for decades, with calibration metrics forming a core component of forecast verification.

Ensemble methods—combining predictions from multiple models—have become a default strategy in machine learning for improving predictive accuracy. Bagging reduces variance, boosting reduces bias, and model averaging provides robustness to individual model misspecification. The benefits of ensembling for discriminative performance, as measured by metrics like area under the receiver operating characteristic curve (AUROC) or classification accuracy, are well documented both theoretically and empirically.

A natural assumption follows: if ensembling improves accuracy, it should also improve calibration. After all, a more accurate model should assign more appropriate probabilities. This assumption, while intuitive, conflates two distinct properties of probabilistic predictions. A model can be highly discriminative (AUROC = 0.95) while being poorly calibrated (consistently overconfident), and conversely, a model can be well-calibrated (predicted probabilities match observed frequencies) while having modest discrimination. Discrimination and calibration are complementary components of probabilistic forecast quality, as formally captured by the Brier score decomposition into reliability (inverse calibration), resolution (discrimination conditional on forecasts), and uncertainty (irreducible base-rate entropy).

The question we address is therefore: Does the act of combining model predictions through ensemble aggregation preserve, improve, or degrade the calibration of the resulting probabilistic forecasts? And if the answer depends on the aggregation method, which methods are calibration-preserving and which are not?

This question is particularly urgent in clinical decision support. The ensemble methods evaluated in our companion study of flow cytometry classification—including simple averaging, trimmed mean, median aggregation, and others applied to nine base classifiers across five diagnostic tasks—demonstrated clear improvements in AUROC. The TrimmedMean ensemble achieved mean AUROC of 0.880 across tasks, outperforming the best individual base model (Random Forest at 0.864). But AUROC measures only rank ordering, not calibration. A perfectly discriminating ensemble that assigns all positive cases a probability of 0.99 and all negative cases a probability of 0.01 has AUROC = 1.0 but may be catastrophically miscalibrated if the true prevalence is 25%.

In this work, we conduct a systematic Monte Carlo study to isolate the effect of ensemble aggregation on probability calibration, controlling for all other factors by working with synthetic data where the ground truth calibration properties of each base model are known by construction.

1.1 Contributions

Our contributions are fourfold:

1. A controlled experimental framework. We design a Monte Carlo simulation where five base models have precisely specified calibration properties (well-calibrated, overconfident, underconfident, shifted-high, shifted-low), enabling clean measurement of how aggregation transforms calibration.

2. A surprising negative result. Simple averaging—the most widely used ensemble aggregation—worsens calibration relative to a well-calibrated individual model, with a large effect size (Cohen's d = 2.89). This holds across 500 Monte Carlo replicates with overwhelming statistical confidence.

3. Aggregation-specific recommendations. The geometric mean is the only tested aggregation that consistently improves calibration. The trimmed mean and median provide partial mitigation but are not reliably calibration-improving.

4. Post-hoc recalibration analysis. We evaluate Platt scaling and isotonic regression applied after ensemble aggregation, demonstrating that post-hoc methods can partially rescue miscalibrated ensembles but introduce their own variance, particularly in small-sample settings.

2. Background

2.1 Probability Calibration

A probabilistic classifier f is said to be perfectly calibrated if, for all probability values p in [0, 1]:

P(Y = 1 | f(X) = p) = p

In other words, among all instances where the model predicts probability p, the true positive fraction is exactly p. Perfect calibration is a strong condition rarely achieved in practice; applied work therefore focuses on measuring the degree of miscalibration and applying post-hoc corrections.

Calibration is often visualized through reliability diagrams, which plot observed frequency against predicted probability across binned prediction ranges. A perfectly calibrated model produces points along the diagonal. Deviations above the diagonal indicate underconfidence (the model predicts lower probabilities than the true frequency), while deviations below indicate overconfidence.

2.2 Calibration Metrics

Expected Calibration Error (ECE). The most widely used calibration metric partitions predictions into B equal-width bins and computes a weighted average of the per-bin calibration gap:

ECE = Σ_{b=1}^{B} (n_b / N) |acc(b) − conf(b)|

where n_b is the number of samples in bin b, N is the total sample size, acc(b) is the observed accuracy (positive fraction) in the bin, and conf(b) is the mean predicted probability in the bin. ECE = 0 indicates perfect calibration.

Maximum Calibration Error (MCE). The maximum per-bin calibration gap: MCE = max_b |acc(b) − conf(b)|. MCE is particularly relevant in safety-critical applications where even a single probability range with poor calibration could lead to harmful decisions.

Brier Score. The mean squared error of probabilistic predictions:

BS = (1/N) Σ_{i=1}^{N} (f_i − y_i)²

where f_i is the predicted probability and y_i is the binary outcome. The Brier score simultaneously captures calibration and discrimination, with lower values indicating better probabilistic predictions.

Brier Score Decomposition. Murphy decomposed the Brier score into three interpretable components:

BS = Reliability − Resolution + Uncertainty

• Reliability measures calibration error (lower is better): how much predicted probabilities deviate from observed frequencies within each forecast bin. This is closely related to ECE but uses squared rather than absolute differences.
• Resolution measures discrimination (higher is better): how much the observed frequency varies across forecast bins. A model whose predictions carry no information about outcomes has zero resolution.
• Uncertainty is the inherent entropy of the outcome, determined entirely by the base rate: bar_o × (1 − bar_o). This is constant across models for a given dataset.

The decomposition is valuable because it reveals whether a change in Brier score is driven by calibration (reliability) or discrimination (resolution). An ensemble could improve the Brier score by increasing resolution while simultaneously worsening reliability—a tradeoff invisible to the aggregate Brier score alone.

2.3 Post-Hoc Calibration Methods

When a trained model produces miscalibrated probabilities, post-hoc calibration methods can be applied to transform the output probabilities toward better calibration, without modifying the underlying model.

Platt Scaling. Originally developed for calibrating support vector machine outputs, Platt scaling fits a logistic regression model to transform raw model scores into calibrated probabilities. Given model output s, the calibrated probability is:

p_calibrated = σ(a × s + b)

where σ is the sigmoid function and a, b are parameters fitted by maximum likelihood on a held-out calibration set. Platt scaling is a parametric method that assumes the calibration function is sigmoidal—an assumption that holds well for many models but can fail when the miscalibration is non-monotonic or has complex structure.

Isotonic Regression. A non-parametric alternative, isotonic regression fits a monotonically non-decreasing step function mapping raw model outputs to calibrated probabilities. It makes minimal assumptions about the form of the calibration function but requires more data to achieve stable estimates, as it has more degrees of freedom than the two-parameter Platt scaling. Niculescu-Mizil and Caruana conducted a comprehensive comparison of calibration methods across multiple classifier families, finding that Platt scaling works well for sigmoid-shaped distortions (common in SVMs and boosted models), while isotonic regression is more flexible but prone to overfitting with small calibration sets.

2.4 Calibration and Ensembles

The relationship between ensemble aggregation and calibration has received surprisingly little systematic study. Several observations from the existing literature are relevant:

Averaging and the Convexity of the Brier Score. The Brier score is a strictly proper scoring rule: it is minimized when the predicted probability equals the true conditional probability. Moreover, it is convex: for any convex combination of predictions, the Brier score of the combination is less than or equal to the convex combination of the individual Brier scores. This means simple averaging can never increase the Brier score beyond the average of the individual Brier scores. However, this property applies to the composite Brier score—it does not guarantee that the reliability (calibration) component individually improves.

The Calibration-Averaging Paradox. Consider two models: one always predicts 0.3 for positive cases, and another always predicts 0.7 for positive cases. If the true frequency is 0.5, each model has an ECE contribution from its respective bin. The average model predicts 0.5 for these cases—perfectly calibrated. But now consider: one model predicts 0.2 and another predicts 0.8, with the true frequency being 0.3. Model 1 has |0.2 − 0.3| = 0.1 error and Model 2 has |0.8 − 0.3| = 0.5 error. Their average prediction of 0.5 has |0.5 − 0.3| = 0.2 error—better than Model 2 but worse than Model 1. Averaging can help or hurt calibration depending on the nature and symmetry of individual model miscalibrations.

Diversity and Calibration. Ensemble diversity—the degree to which individual models make different errors—is well studied in the context of accuracy. Its role in calibration is less clear. Models that are miscalibrated in different directions (one overconfident, one underconfident) may partially cancel when averaged, but models that are miscalibrated in the same direction will not benefit from averaging.

3. Monte Carlo Experimental Design

3.1 Simulation Framework

We design a Monte Carlo simulation to isolate the effect of ensemble aggregation on calibration. Each trial consists of the following steps:

1. Generate true probabilities. For each of N = 546 synthetic patients, draw a true event probability from a Beta(2, 5) distribution. This produces a skewed distribution centered near 0.29 (mean = 2/7 ≈ 0.286), reflecting the clinical reality that many diagnostic outcomes have moderate-to-low prevalence.

2. Generate binary outcomes. For each patient, draw a binary outcome y_i ~ Bernoulli(p_i), where p_i is the true probability. This ensures the ground-truth relationship between probability and outcome is exact.

3. Generate base model predictions. Five models produce predictions by applying different calibration-distortion functions to the true probabilities, plus Gaussian noise (σ = 0.1):

   • Well-calibrated: f(p) = p (identity function; produces predictions that are, on average, calibrated)
   • Overconfident: f(p) = σ(3(p − 0.5)) where σ is the sigmoid function (pushes predictions toward 0 and 1)
   • Underconfident: f(p) = σ(0.5(p − 0.5)) + 0.25 (compresses predictions toward the center)
   • Shifted-high: f(p) = clip(p + 0.15, 0, 1) (systematically predicts too high)
   • Shifted-low: f(p) = clip(p − 0.15, 0, 1) (systematically predicts too low)

4. Apply ensemble aggregation. Five aggregation methods combine the five base model predictions into ensemble predictions.

5. Compute calibration metrics. ECE, MCE, Brier score, and Brier decomposition are computed for each ensemble and each individual model.

6. Apply post-hoc calibration (subset of trials). Platt scaling and isotonic regression are applied to ensemble predictions on held-out data.

3.2 Sample Size and Replicate Count

We use N = 546 patients per trial, matching the sample size of our companion clinical study (flow cytometry classification across five diagnostic tasks). This ensures our simulation results are directly relevant to the practical setting that motivated this investigation.

We conduct R = 500 Monte Carlo replicates. With 500 replicates, the standard error of any estimated metric mean is approximately σ/√500 ≈ 0.045σ, providing precise estimates of the expected behavior of each method. For win-rate comparisons (binary outcomes per trial), 500 replicates give 95% confidence intervals of approximately ±4.4 percentage points for observed rates near 50%, and much tighter intervals for extreme rates.

3.3 Ensemble Aggregation Methods

We evaluate five aggregation methods:

Simple Average. The arithmetic mean across all five base models: p_ens(x) = (1/5) Σ_k p_k(x). This is the most widely used aggregation and represents the default choice in most ensemble implementations.

Trimmed Mean (20%). The mean after removing the highest and lowest predictions for each patient. With five models, a 20% trim removes one model from each end, effectively computing the mean of the three central predictions. This is more robust to extreme miscalibrations.

Median. The middle value across five models. The median is maximally robust to outliers—it depends only on the rank order of predictions, not their magnitudes.

Geometric Mean. The exponential of the mean log-probability: p_ens(x) = exp((1/5) Σ_k log p_k(x)). The geometric mean is less sensitive to large predictions and naturally produces values biased toward lower probabilities when the input distribution is skewed.

Best-3 Average. The average of the three models with the highest expected calibration (in our simulation: well-calibrated, overconfident, and underconfident). This represents an "oracle" subset-selection strategy where the practitioner has some knowledge of which models are likely to be better calibrated.

3.4 Calibration Profile Rationale

The five calibration profiles are designed to represent common patterns of miscalibration encountered in practice:

• Well-calibrated models are achievable through careful training or post-hoc calibration. Logistic regression, when the model is correctly specified, tends to produce well-calibrated probabilities.
• Overconfident models push predictions toward extremes. This pattern is common in deep neural networks and gradient boosted trees.
• Underconfident models compress predictions toward the center. This can arise from excessive regularization or from averaging over a posterior that includes many weak hypotheses.
• Shifted-high and shifted-low models have systematic bias in one direction. This occurs when training and deployment populations differ in prevalence (a form of dataset shift) or when a model's intercept is miscalibrated.

The combination of these five profiles creates an ensemble where the miscalibrations are asymmetric—they do not cancel under simple averaging, because the underconfident and shifted models produce qualitatively different distortions. This asymmetry is realistic: in practice, one rarely has the luxury of models whose miscalibrations are perfectly complementary.

4. Results

4.1 Expected Calibration Error

The ECE results across 500 Monte Carlo replicates reveal a clear hierarchy of methods:

Several findings emerge:

Finding 1: Simple averaging degrades calibration. The simple average ensemble (ECE = 0.112) is substantially worse than the well-calibrated individual (ECE = 0.054). The 95% confidence intervals do not overlap. This result is not marginal—it represents a 106% increase in ECE.

Finding 2: The geometric mean is the only calibration-improving ensemble. The geometric mean achieves ECE = 0.048, which is better than every individual model including the well-calibrated one. This is the only aggregation method whose 95% CI for ECE lies entirely below the well-calibrated individual's mean.

Finding 3: Median and trimmed mean provide intermediate calibration. Both methods (ECE ≈ 0.078–0.081) are worse than the well-calibrated individual but substantially better than simple averaging. These robust estimators successfully mitigate the most extreme miscalibrations but do not fully preserve the calibration of the best individual model.

Finding 4: Selective averaging (Best-3) performs worst among ensembles. Despite selecting what should be the three best-calibrated models, the best-3 average (ECE = 0.171) is worse than the full ensemble average (ECE = 0.112). Including the overconfident model in a smaller average amplifies its distortion relative to the well-calibrated model.

4.2 Win Rate Analysis

To complement the aggregate ECE comparison, we compute the fraction of trials in which each ensemble achieves lower ECE than the well-calibrated individual:

The win rates are strikingly asymmetric. Simple averaging beats the well-calibrated individual in only 3 out of 500 trials—a near-zero probability event. The geometric mean wins in 67.2% of trials, making it the only method that is likely to improve calibration in any given application. Even the median and trimmed mean win less than 11% of the time, meaning that if you have access to a well-calibrated individual model, ensembling with these methods is more likely to hurt than help.

4.3 Effect Size Analysis

The Cohen's d effect sizes quantify the magnitude of calibration difference between each ensemble and the well-calibrated individual, in units of standard deviations of the paired difference:

By conventional standards, Cohen's d > 0.8 is a "large" effect. The simple average's d = +2.89 indicates a very large calibration degradation—in practically every trial, simple averaging produces meaningfully worse calibration. The geometric mean's d = −0.44 is a moderate improvement—noticeable and consistent but not as dominant as the degradation from other methods.

4.4 Brier Score Results

The Brier score rankings partially but incompletely mirror the ECE rankings. The key divergence is that the trimmed mean and median achieve Brier scores very close to the well-calibrated individual (0.188, 0.189 vs. 0.188), despite having substantially worse ECE (0.081, 0.078 vs. 0.054). This discrepancy is explained by the Brier decomposition: these ensembles partially compensate for worse reliability (calibration) with slightly better resolution (discrimination).

4.5 Brier Score Decomposition

The decomposition reveals the mechanism by which different aggregation methods affect the Brier score:

The uncertainty component is constant (≈ 0.204) across methods, as expected—it depends only on the outcome base rate, which is shared across methods within each trial.

Key insight from decomposition: The geometric mean achieves its superior Brier score through both mechanisms: it has the lowest reliability (best calibration, 0.0041) and the highest resolution (best discrimination, 0.0240). This dual advantage is unusual—most methods show a tradeoff between reliability and resolution.

The simple average shows an instructive pattern: it has good resolution (0.0234, second highest among ensembles) but poor reliability (0.0147, worst among ensembles). This confirms that simple averaging preserves the ensemble's discriminative diversity while degrading its calibration—the model knows who is positive and negative but assigns systematically incorrect probabilities.

The trimmed mean and median improve reliability substantially over simple averaging (0.0086 vs. 0.0147) while maintaining comparable resolution (0.0238 and 0.0224 vs. 0.0234). This confirms that robust aggregation achieves calibration improvement primarily by reducing the influence of extreme miscalibrations, without sacrificing discrimination.

5. When Ensembles Hurt Calibration

5.1 The Mechanism of Calibration Degradation

Why does simple averaging degrade calibration when it demonstrably improves discriminative performance? The answer lies in the asymmetry of the calibration distortions in the component models.

Consider the prediction for a patient with true probability p = 0.3:

• Well-calibrated model: predicts ≈ 0.30 (plus noise)
• Overconfident model: predicts σ(3 × (0.3 − 0.5)) = σ(−0.6) ≈ 0.354 (pushed toward 0 less than expected because of the sigmoid shape)
• Underconfident model: predicts σ(0.5 × (0.3 − 0.5)) + 0.25 = σ(−0.1) + 0.25 ≈ 0.475 + 0.25 = 0.725... but clipped: actually σ(−0.1) + 0.25 ≈ 0.475 + 0.25 ≈ 0.725, but more precisely σ(−0.1) ≈ 0.475, so underconfident predicts ≈ 0.725 (dramatically too high)
• Shifted-high: predicts ≈ 0.45
• Shifted-low: predicts ≈ 0.15

The average of these five predictions is approximately (0.30 + 0.35 + 0.73 + 0.45 + 0.15)/5 ≈ 0.40. The true probability is 0.30, so the average is biased upward by 0.10—a substantial miscalibration driven primarily by the underconfident model's extreme distortion.

The critical insight is that the underconfident model's distortion is not symmetric with the overconfident model's distortion. The overconfident model pushes predictions toward the extremes (toward 0 for low-probability events, toward 1 for high-probability events), while the underconfident model compresses everything toward a central value near 0.5–0.7. These distortions do not cancel under arithmetic averaging. The shifted-high and shifted-low models do partially cancel, but the overall average remains biased.

5.2 Why the Geometric Mean Helps

The geometric mean p_geo = exp((1/K) Σ_k log p_k) operates in log-probability space. This has two calibration-relevant properties:

Property 1: Sensitivity to low predictions. In log space, a prediction of 0.01 becomes log(0.01) = −4.61, while a prediction of 0.50 becomes log(0.50) = −0.69. The geometric mean therefore weights low-probability predictions much more heavily than high-probability predictions. In our skewed scenario (Beta(2,5) true probabilities, mean ≈ 0.29), most patients have low true probabilities, and the geometric mean's natural bias toward low values helps counteract the upward biases of the underconfident and shifted-high models.

Property 2: Multiplicative calibration. If the calibration errors are approximately multiplicative (e.g., one model predicts 1.5× the true probability and another predicts 0.67× the true probability), the geometric mean perfectly cancels these errors: (1.5p × 0.67p)^0.5 = p. The arithmetic mean would yield (1.5p + 0.67p)/2 = 1.085p—a biased estimate. Many calibration distortions in practice have approximately multiplicative structure, making the geometric mean a natural choice.

5.3 The Best-3 Paradox

Counterintuitively, the best-3 average (selecting the well-calibrated, overconfident, and underconfident models) produces worse calibration (ECE = 0.171) than the full 5-model average (ECE = 0.112). This occurs because:

1. The shifted-high and shifted-low models, while individually poorly calibrated, partially cancel each other's biases when included in the full average.
2. Removing them removes this cancellation while retaining the highly miscalibrated underconfident model, whose extreme distortion now has proportionally greater influence (1/3 weight vs. 1/5 weight).

This finding has a practical lesson: selecting a "better" subset of models does not guarantee better ensemble calibration. The calibration of an ensemble depends on the interaction structure of component miscalibrations, not on the calibration quality of individual components. A poorly calibrated model that is miscalibrated in the opposite direction to other models can actually improve ensemble calibration through error cancellation.

6. Post-Hoc Calibration of Ensembles

6.1 Platt Scaling After Ensembling

We applied Platt scaling to the three main ensemble methods (simple average, trimmed mean, median) using a train-test split within each trial. Results are averaged over 50 replicates (every 10th of 500 trials):

Platt scaling provides substantial ECE improvement for all three methods, with the largest absolute improvement for simple averaging (0.112 → 0.063, a 43.8% reduction). After Platt scaling, the three methods converge to similar ECE values near 0.060–0.063, suggesting that Platt scaling largely corrects for the aggregation-specific calibration distortions.

However, Platt-calibrated simple averaging (ECE = 0.063) is still worse than the geometric mean without any post-hoc calibration (ECE = 0.048). This suggests that the geometric mean preserves calibration-relevant information that Platt scaling cannot fully recover.

6.2 Isotonic Regression After Ensembling

Isotonic regression also improves calibration but is less effective than Platt scaling in this setting, likely because with N/2 = 273 calibration samples, the non-parametric isotonic regression has insufficient data to achieve stable estimates. This is consistent with the general finding that isotonic regression requires more data than Platt scaling due to its greater flexibility.

6.3 The Calibration Recovery Hierarchy

Combining our findings, we can rank strategies for achieving well-calibrated ensemble predictions, from best to worst:

1. Geometric mean (no post-hoc): ECE = 0.048
2. Median + Platt scaling: ECE ≈ 0.060
3. Trimmed mean + Platt scaling: ECE ≈ 0.063
4. Simple average + Platt scaling: ECE ≈ 0.063
5. Well-calibrated individual (no ensemble): ECE = 0.054
6. Median (no post-hoc): ECE = 0.078
7. Trimmed mean (no post-hoc): ECE = 0.081
8. Simple average (no post-hoc): ECE = 0.112

This hierarchy reveals that the choice of aggregation method matters more than whether post-hoc calibration is applied. The geometric mean without any post-hoc correction outperforms all other methods with Platt scaling. If the geometric mean is not an option (e.g., because some models may produce zero probabilities, where the log is undefined), median + Platt scaling is the next best strategy.

7. Discussion and Recommendations

7.1 Implications for Clinical Decision Support

In clinical settings, probability calibration directly impacts decision-making. A clinician relying on a model's predicted probability of sepsis to decide whether to order blood cultures needs those probabilities to be well-calibrated. If a simple-average ensemble reports 45% sepsis risk when the true risk is 30%, the clinician may order unnecessary tests, increasing costs and potentially subjecting the patient to unnecessary interventions.

Our findings suggest that clinical ML pipelines using ensemble methods should:

1. Prefer geometric mean aggregation when calibration is the primary concern (e.g., risk calculators, probability-based decision thresholds).
2. Apply Platt scaling after any arithmetic aggregation (simple average, trimmed mean) to mitigate the calibration degradation.
3. Evaluate ensemble calibration separately from discrimination. An ensemble may improve AUROC substantially while degrading the very calibration properties that clinicians rely upon for decision-making.
4. Report calibration metrics (ECE, reliability diagrams) alongside discrimination metrics (AUROC, AUPRC) in model validation studies. Presenting only AUROC can mask severe calibration problems in ensemble models.

7.2 When to Ensemble for Calibration vs. Discrimination

Our results establish a fundamental dichotomy: ensembling for accuracy and ensembling for calibration are different objectives that may require different aggregation strategies.

If the goal is pure discrimination (AUROC), simple averaging or trimmed mean of diverse models is effective, as confirmed by the companion flow cytometry study where TrimmedMean achieved the highest mean AUROC (0.880).

If the goal is calibrated probabilities, the geometric mean should be preferred. If both objectives matter—as they typically do in clinical applications—a two-stage strategy may be warranted:

1. Use trimmed mean or simple averaging for the final discriminative ranking.
2. Apply Platt scaling on held-out data to recalibrate the averaged probabilities.

Alternatively, use the geometric mean, which achieves both good discrimination and good calibration simultaneously, at the cost of requiring all base model predictions to be strictly positive (no zero probabilities).

7.3 Relationship to Ensemble Diversity

The concept of ensemble diversity—the extent to which base models make different errors—is well studied in the accuracy context, where the bias-variance-covariance decomposition shows that averaging benefits from diverse, uncorrelated errors. Our calibration results add a nuance: diversity in calibration properties does not uniformly help.

Diversity that is symmetric (one model overconfident by δ, another underconfident by δ) benefits averaging by cancellation. Diversity that is asymmetric (one model dramatically miscalibrated in one direction, others moderately miscalibrated in various directions) can degrade averaging by introducing bias that does not cancel. In practice, calibration diversity is almost always asymmetric, which explains why simple averaging typically hurts calibration in our simulation.

7.4 Connection to Real Ensemble Results

Our companion study evaluated 15 ensemble aggregation methods on five flow cytometry classification tasks with 9 base classifiers. The AUROC results showed TrimmedMean (0.880) and SimpleAverage (0.874) as top-performing methods. Our calibration analysis suggests that these AUROC gains may come with a calibration cost. If the base classifiers have heterogeneous calibration properties—which is likely, given that they include Naive Bayes (known to be poorly calibrated), gradient boosted machines (known to be overconfident), and logistic regression (known to be well-calibrated)—then the simple average ensemble may have worse calibration than the best individual base classifier.

This prediction could be tested empirically by computing ECE and reliability diagrams for each ensemble method on the clinical data, alongside the already-reported AUROC values. We note that with 546 patients and 5-fold cross-validation, the calibration estimates will have substantial variance, potentially requiring the Monte Carlo perspective developed in this paper to interpret the observed variation.

8. Limitations

8.1 Synthetic Data

Our simulation uses synthetic data with known calibration properties, which enables clean causal interpretation but may not fully capture the complexity of real-world miscalibrations. Real models may exhibit non-stationary calibration (varying across different subpopulations), interaction effects between features and calibration, and more complex distortion functions than the five profiles we examine.

8.2 Number of Base Models

We evaluate ensembles of five base models. As the number of models increases, the law of large numbers suggests that arithmetic averaging of independent calibration errors should improve calibration. Our results apply most directly to the common setting of small-to-moderate ensemble sizes (3–10 models). Very large ensembles (e.g., random forests with hundreds of trees) may exhibit different calibration dynamics due to the averaging effect across many weakly correlated predictions.

8.3 Equal Weighting

All our ensemble methods weight models equally. Weighted ensembles, where weights are optimized for calibration (e.g., minimizing ECE or the reliability component of the Brier score on validation data), could potentially outperform the geometric mean. We do not explore calibration-optimized weighting, which represents a promising direction for future work.

8.4 ECE Bin Sensitivity

ECE depends on the number of bins used for computation. We use 10 bins throughout, following common practice. With 546 samples and 10 bins, some bins may be sparsely populated, introducing noise in the ECE estimate. Adaptive binning strategies or kernel-based calibration metrics could provide more robust estimates but would complicate comparison across methods.

8.5 Single Prevalence Regime

Our Beta(2, 5) true probability distribution produces a specific prevalence regime (mean ≈ 29%, moderately skewed). The relative performance of aggregation methods may differ under very low prevalence (e.g., rare disease screening), balanced prevalence, or high prevalence settings. The geometric mean's advantage may be particularly pronounced in low-prevalence settings where its bias toward lower probabilities aligns with the data distribution.

8.6 Post-Hoc Calibration Sample Size

Our Platt scaling and isotonic regression results use half the data (273 samples) for calibration fitting. In settings with very small calibration sets, the variance of post-hoc calibration could overwhelm its bias-reduction benefits. We observe this effect in isotonic regression, which underperforms Platt scaling likely due to its greater data requirements.

9. Conclusion

We have presented a systematic Monte Carlo investigation of how ensemble aggregation methods affect probability calibration. Our central finding is that the most common aggregation strategy—simple arithmetic averaging—substantially degrades calibration relative to a well-calibrated individual model, with a Cohen's d effect size of +2.89 and a win rate of only 0.6% across 500 trials. This result has been overlooked because ensemble evaluations typically focus on discriminative metrics (AUROC, accuracy) that are insensitive to calibration.

The geometric mean emerges as the clearly preferred aggregation method for calibration-sensitive applications, achieving the lowest ECE (0.048), the highest win rate against the well-calibrated individual (67.2%), and the best Brier decomposition (lowest reliability, highest resolution). Its advantage derives from operating in log-probability space, which is naturally suited to the approximately multiplicative structure of many calibration distortions.

Our Brier decomposition analysis reveals the mechanism: simple averaging preserves resolution (discriminative diversity) while degrading reliability (calibration). The geometric mean uniquely preserves both components. Post-hoc Platt scaling can partially rescue the calibration of arithmetic ensembles but cannot match the geometric mean's inherent calibration-preserving properties.

For practitioners in clinical machine learning and other calibration-sensitive domains, we offer three actionable recommendations:

1. Prefer geometric mean aggregation when deploying ensembles for probability estimation.
2. Always apply post-hoc calibration (Platt scaling preferred over isotonic regression for small samples) after arithmetic averaging if the geometric mean is not feasible.
3. Report calibration metrics alongside discrimination metrics in ensemble evaluation. AUROC gains from ensembling may mask calibration degradation that undermines the clinical utility of predicted probabilities.

The broader lesson is that ensembling is not a monolithic operation: the choice of aggregation function determines which properties of the component models are preserved and which are degraded. For discrimination, arithmetic averaging leverages the diversity of errors. For calibration, geometric averaging leverages the structure of probability space. Recognizing this distinction is essential for building ensemble systems that are not just accurate, but trustworthy.

References

Brier, G. W. (1950). Verification of forecasts expressed in terms of probability. Monthly Weather Review, 78(1), 1–3.

Dawid, A. P. (1982). The well-calibrated Bayesian. Journal of the American Statistical Association, 77(379), 605–610.

Gneiting, T., and Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477), 359–378.

Murphy, A. H. (1973). A new vector partition of the probability score. Journal of Applied Meteorology, 12(4), 595–600.

Niculescu-Mizil, A., and Caruana, R. (2005). Predicting good probabilities with supervised learning. In Proceedings of the 22nd International Conference on Machine Learning, pages 625–632.

Platt, J. C. (1999). Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. In Advances in Large Margin Classifiers, pages 61–74. MIT Press.

Zadrozny, B., and Elkan, C. (2002). Transforming classifier scores into accurate multiclass probability estimates. In Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 694–699.

DeGroot, M. H., and Fienberg, S. E. (1983). The comparison and evaluation of forecasters. The Statistician, 32(1/2), 12–22.

Appendix A: Detailed Reliability Analysis

A.1 Maximum Calibration Error

The MCE results complement the ECE findings by highlighting worst-case bin behavior:

The geometric mean's MCE advantage is proportionally larger than its ECE advantage: MCE = 0.162 vs. 0.175 for the well-calibrated individual (7.4% improvement). This suggests the geometric mean not only improves average calibration but also reduces worst-case miscalibration across the probability range.

A.2 Brier Score Variance Reduction

A notable property of ensemble methods is variance reduction. The standard deviation of the Brier score across trials is:

All ensemble methods produce lower Brier variance than individual models, confirming that ensembling provides prediction stability even when calibration is degraded. The simple average and best-3 average have particularly low variance, reflecting the averaging effect across multiple models.

Appendix B: Experimental Code

The complete Monte Carlo simulation was implemented in Python using NumPy and SciPy. The Expected Calibration Error is computed using equal-width binning:

def compute_ece(y_true, y_pred, n_bins=10):
    bin_edges = np.linspace(0, 1, n_bins + 1)
    ece = 0.0
    n = len(y_true)
    for i in range(n_bins):
        if i < n_bins - 1:
            mask = (y_pred >= bin_edges[i]) & (y_pred < bin_edges[i+1])
        else:
            mask = (y_pred >= bin_edges[i]) & (y_pred <= bin_edges[i+1])
        nb = mask.sum()
        if nb == 0:
            continue
        ece += (nb / n) * abs(y_pred[mask].mean() - y_true[mask].mean())
    return ece

The trimmed mean uses sorted-array slicing for efficiency:

def fast_trimmed_mean(arr, prop=0.2):
    sorted_arr = np.sort(arr, axis=0)
    n = arr.shape[0]
    k = int(n * prop)
    return np.mean(sorted_arr[k:n-k], axis=0)

The Brier score decomposition follows Murphy (1973):

def brier_decomposition(y_true, y_pred, n_bins=10):
    bin_edges = np.linspace(0, 1, n_bins + 1)
    n = len(y_true)
    bar_o = y_true.mean()
    uncertainty = bar_o * (1 - bar_o)
    reliability, resolution = 0.0, 0.0
    for i in range(n_bins):
        # bin membership
        mask = ...  # as above
        nk = mask.sum()
        if nk == 0:
            continue
        ok = y_true[mask].mean()
        fk = y_pred[mask].mean()
        reliability += (nk / n) * (fk - ok) ** 2
        resolution += (nk / n) * (ok - bar_o) ** 2
    return reliability, resolution, uncertainty

Platt scaling uses grid search over sigmoid parameters, and isotonic regression uses the Pool Adjacent Violators algorithm. All random seeds are fixed (numpy.random.seed(42)) for reproducibility. Full code is available in the supplementary repository.