We developed separate models for vaccines and drugs to account for differences in triage, confounders, risk factors, and decision categories. Here, we focused on the vaccine models, which used largely the same methodological framework as the drug models; drug-specific methods and results are presented in the Supplementary Methods and Results.

2.1 Data Sources2.1.1 Alert Data

Alerts were retrieved from the company’s signal management system. Alerts for product-event pairs are routinely generated and classified as: (i) disproportionality analysis (DPA), when the Multi-item Gamma Poisson Shrinker lower bound (5th percentile) of the 90% credibility interval (CI) of the Empirical Bayes Geometric Mean (EBGM) (EB05) is ≥ 2 [24]; (ii) time to onset (TTO; vaccine model only), when two-sample Kolmogorov-Smirnov tests [25] yield p-values <0.01 for TTO distributions within a 30-day window; (iii) designated medical event (DME), the list of which is defined and maintained in-house, when ≥ 4 reports of an unlisted DME (or ≥ 1 for newly launched products) are reported; or (iv) aberration, when ≥ 4 cases exceed the forecast upper bound from a negative binomial time-series regression fitted to the preceding six months [26].

Statistical alerts are generated periodically using cumulative data through each alert-period end. Triage is a repeated process: each review considers new information since the previous triage, but records decisions are based on the cumulative data. For each alert, the reviewer makes one of two primary triage decisions: either to close the alert or to initiate further investigation. When an alert is closed, the reviewer should select a single one of 13 predefined decision-rationale categories (Table 1). An “Other” category is available and should be used in situations where none of the specific categories are judged applicable, or when multiple categories are equally applicable and no single rationale predominates; in such cases, the reviewer must use the free-text comment field to clarify the main rationale(s) for closure. For brevity, we refer to these decision rationale categories as “decisions” in the remainder of the manuscript and in the ML modelling. The retrieved dataset therefore represented one row per product-event-alert-period-end-alert-type combination and included the reviewer’s routine triage outcome, selected decision rationale (Table 1), and free-text comment; counts of spontaneous and total reports; new spontaneous reports since 1, 3, 6, 9 and 12 months; counts of serious and fatal reports (if any); EBGM with its 90% CI (EB05, EB95); and flags for DME, label-listed, serious, and medication error events.

Table 1 Predefined decisions to document closure or investigation of statistical alerts

The initial vaccine model used cumulative alerts through December 2022 for 17 marketed vaccines. The refined model was retrained on an expanded dataset covering all marketed vaccines (n = 36) and including cumulative alerts through July 2024 to capture recent trends.

2.1.2 Individual Case Safety Reports Data

ICSR data were retrieved from the company safety database and included attributes selected for their clinical relevance; several were identified using CUSP. Based on CUSP and expert judgement, we included ICSR attributes such as event duration and outcome, product dosage and indication, action taken with the product, presence of dechallenge or rechallenge, and TTO. Patient information (age, sex, medical history) was also included. These variables (Table 2) were organised by their relation to the report, patient, event, or product.

Table 2 ICSR attributes retrieved from the company safety database

The initial vaccine model was trained and tested on spontaneous ICSRs received through 16 December 2022. The refined model’s train/test dataset included all report types (not only spontaneous) received through 16 September 2024.

2.1.3 Product-Level Metadata

Product-specific features for the refined models were compiled outside the safety database using PV expert knowledge and internal product records. These differed by product type and included vaccine target population, target disease, and novelty status.

2.2 Data Pre-processing2.2.1 Alert Data

Because alerts were generated repeatedly, the dataset was restricted to the most recent alert per product-event pair, so each pair contributed a single observation. However, to increase the number of samples for the decision “Merits further investigation”, we kept the most recent alert with that decision.

During pre-processing, fields were added: counts of how often each product-event pair triggered each alert type (DPA/TTO/DME/aberration); the immediately preceding previous decision and period‑end for the retained alert; and a flag for alerts migrated from legacy signal management systems that may have been triggered or triaged under earlier processes. Historical alerts migrated from legacy systems with missing statistical scores and event characteristics were excluded to ensure data completeness.

For the refined models, PV experts reviewed alerts with multiple closing decisions in the train/test dataset and, where possible, selected the single most relevant decision. Alerts remaining with multiple decisions were excluded to avoid many low-frequency combinations that would impair model training. For product-event pairs with at least one “Event being monitored, no change” decision, we selected the most recent alert carrying that decision instead of the most recent alert overall, to increase the number of observations.

The “Other” decision category is non-specific and varies in usage. This category was therefore reviewed by experts and, where appropriate, alerts were reclassified into more specific decision categories. Alerts that had been reclassified were then included under the corresponding category for training/testing. For example, “Other” was sometimes selected to indicate that the cases in the alert remained consistent with the previous review and that the decision was unchanged, although the same decision as previously should instead have been selected. In other instances, “Other” was chosen even though the free-text comment aligned with existing decisions such as “Insufficient information available” or “Event term is uninformative or non-specific”. For this study, only alerts with a confirmed “Other” decision after this review were retained as “Other” for training/testing.

2.2.2 Individual Case Safety Reports Data

Underrepresented values for event outcome and action taken were grouped with related, more frequent categories.

Several attributes were converted to binary flags (e.g., dose_flag = 1 if dose, daily dose, or total regimen dose > 0 with units provided, or if dose number > 0). Because product indications, lot numbers, and reporter countries had many distinct non-null values, we counted distinct values per product-event pair and counted missing values for indication and lot number.

For the refined models, TTO was supplemented with free text timing expressions (e.g., immediately; less than an hour) that could be mapped to the same TTO categories. For vaccines, missing or unknown actions taken were set to “Not applicable”. Dechallenge and rechallenge were retained when explicitly reported as positive or negative and set to “Not applicable” otherwise.

Further pre-processing details and variable definitions are provided in Supplementary Table 1.

2.3 Feature Engineering

TTO and event duration were categorised as: 0–< 24 h (hours), 24 h–< 7 days (days), 7–< 30 days (weeks), 30–< 365 days (months), ≥ 365 days (years); negative or missing values were coded as “Unknown”.

Alerts and ICSRs were linked by product and Medical Dictionary for Regulatory Activities (MedDRA) Preferred Term (PT) and by matching report dates within the alert period. Reports with significant follow-up after the alert-period end were excluded to ensure alignment with the alert timeline. Where multiple ICSRs matched a single alert, they were aggregated at the product-event alert level by counting cases for each category value and expressing these counts as proportions of the total cases underlying the alert [21].

Categorical variables were “one-hot encoded” (i.e., converted into multiple yes/no (1/0) variables, one per category, indicating whether that category applied) to make the data compatible with ML models [27]. A feature correlation matrix was generated, and highly correlated attributes were removed after consultation with PV experts. For the initial model, the DME and event listed in product label flags from the alert data source were removed because they were highly correlated with the corresponding flags from the ICSR data source and were missing for some historical alerts. The number of spontaneous reports (Spont_N), which was strongly correlated with the number of fatal reports (Fatal_N), and the EBGM for spontaneous reports (Spont_EBGM), which was highly correlated with the lower 90% confidence bound of the EBGM for spontaneous reports (Spont_EB05) were also removed. Table 3 summarises the features retained for the initial model after feature selection.

Table 3 Features included in the initial model

For the refined models, we applied targeted feature engineering guided by explainability reviews and PV expert input. Redundant binary features were consolidated to retain a single representative value. Highly correlated features were removed. Table 4 lists the features used in the refined models, comprising a subset of those retained for the initial model plus several newly engineered features identified during expert review.

Table 4 Features included in the refined model2.4 Machine Learning Framework

We evaluated several multiclass algorithms (Decision Tree [28], K-Nearest Neighbours [29], Linear Support Vector Machine [30], Logistic Regression, Random Forest [31], and XGBoost [32, 33]), prioritising interpretability alongside performance given the regulatory context. XGBoost is a supervised tree-based ensemble method that implements an optimised gradient-boosting framework and is well suited to fast, scalable learning from structured (tabular) data.

Models were developed in Python (v3.8.10) using Scikit-learn (v0.24.1) [34] and executed on JupyterHub and Databricks. Performance was assessed on a held-out test set using a combination of class-level precision, recall, and F1-score, with weighted-average F1 as the primary selection metric to account for class imbalance. Both Random Forest and XGBoost achieved acceptable accuracy (> 60%). After hyperparameter tuning and applying class weights, XGBoost performed best on our large, sparse, high-dimensional, imbalanced dataset (data not shown).

Hyperparameters for the XGBoost (v2.1.1.) classifier were optimised via grid search [34] over learning rate, tree depth, number of estimators, and regularisation terms using a 3-fold cross-validation to reduce overfitting, chosen as a compromise between variance reduction and computational cost given the large dataset and model complexity. The weighted F1-score was used as the tuning metric to address the imbalanced multiclass setting. We used the “multi:softprob” objective to produce per-class probabilities for each instance, and used these raw, uncalibrated probabilities to assess prediction confidence and output reliability.

Models were trained on 75% of the data and tested on the remaining 25% (test dataset). To maintain representativeness across all alert triage classes and remove any potential influences of data ordering, random stratified sampling was used during data partitioning. This effectively mitigated sampling bias, ensuring an equitable distribution of examples in both the training and testing sets. To address class imbalance, class weights were adjusted during training to penalise minority-class misclassification more heavily, encouraging the model to focus on underrepresented classes. Model predictions on the test dataset were compared with routine triage decisions made by safety experts.

2.5 Performance Evaluation

Predictive performance was evaluated on the test dataset using classification reports with class-level precision (proportion of predicted positives that are true positives), recall (sensitivity, proportion of true positives correctly identified), and F1-score (harmonic mean of precision and recall). Aggregated metrics, summarising performance across all classes, included macro-average F1-score (equal class weighting) and weighted-average F1-score (accounting for class imbalance). This is particularly important for imbalanced datasets, where accuracy alone can be misleading.

Models produce a probability for every class; the highest probability class becomes the model’s predicted label (the “top” or 1st-ranked prediction) and is assessed by standard metrics (confusion matrix, precision/recall/F1). We also examined the ranked probabilities (2nd/3rd-ranked) to interpret model uncertainty and alternative plausible classifications and referred to these explicitly when discussing ordered class probabilities.

A confusion matrix [35] visualised prediction performance across decision classes, aiding identification of specific error types and supporting explainability analyses and model refinement. Because triage is a repeated, time-dependent process, we also performed retrospective time-based evaluations: models were trained on alerts up to an earlier cut-off and then applied to alerts from subsequent months that had already been triaged but were unseen by the models. Performance on these time-based test sets was similar to that on the main held-out test set (data not shown), supporting the models’ ability to generalise to future alerts without changing the overall conclusions.

2.6 Explainability and Refinement

Explainability analyses using SHAP (v0.44.1) [36] were performed to improve the model. SHAP quantifies feature importance and the direction/magnitude of each feature’s contribution. We assessed importance at three levels. Model-level global SHAP values identified features that most strongly influence predictions and revealed key patterns and interactions. Class-level global explanations clarified decision logic for each decision category. Sample-level local SHAP values explained individual predictions by showing which features pushed the prediction toward or away from specific classes. Misclassified test samples were examined by comparing SHAP values for the true and predicted classes to identify features driving incorrect predictions and systematic error patterns.

Model refinement proceeded iteratively over a two-year period with PV experts: (1) train model; (2) generate SHAP explanations (global summaries and local plots for misclassifications); (3) expert review; (4) apply refinements (feature additions/removals, recoding, procedural changes); and (5) retrain model. This was repeated until stakeholders were satisfied with performance, interpretability and operational fit.

2.7 Prospective Validation

A four-month prospective validation study (November 2024–February 2025) evaluated the refined model’s practical utility. During the study, safety reviewers received the model’s top 3 predictions and key features supporting each prediction; feature importance details were available on request. Reviewers recorded whether they agreed with one of the top 3 predictions or none and noted if any non-selected prediction was nevertheless an acceptable alternative. Reviewers also estimated time saved during triage when using the model versus the standard manual process.

To assess potential bias from model support, a blinded validation was conducted in the first month for one vaccine. Safety reviewers alternated between using model predictions and performing manual triage: half of the PTs were reviewed by Reviewer A with model support, while Reviewer B triaged them manually; the remaining PTs were reviewed with roles reversed.