MGPS vs RGPS — how the scores are calculated

This page walks through the two Bayesian disproportionality algorithms used on every heatmap page — MGPS (score = EBGM) and RGPS (score = ERAM) — as implemented in the site’s Python pipeline (build_jp_heatmap_pages.py / article_algorithms.py), following Harpaz et al. 2022. If you prefer numbers to symbols, jump to the worked Myocarditis example.

The one-sentence difference

Both start from the same contingency-table counts. The divergence is almost entirely in how they form the expected count E, and in which prior they use to shrink n/E.

Shared inputs: observed count and expected count

For a given cumulative month × brand page (e.g. Pfizer through Dec 2022), every preferred term (PT) contributes:

• a / n — distinct reports of that brand that mention the PT
• N — all reports in the country scope (Pfizer + Moderna + Other)
• n_drug — all reports of the brand
• E_MGPS — Mantel–Haenszel expected count, stratified by the article’s 27 sex × age groups

The stratified expected count is computed as:

In code (score_scope):

# a = brand reports of this PT; evS[s] = all-brand event counts in stratum s
e_strat = sum(drug_s[d][s] * evS[s] / N_s[d][s]
              for s in range(N_STRATA) if N_s[d][s])

PRR on the heatmaps is still the unstratified ratio [a/(a+b)]/[c/(c+d)]. MGPS and RGPS both use stratified (or regression-adjusted) expecteds for their Bayesian scores.

MGPS → EBGM

Model

For each PT i, the relative reporting ratio λ_i has a two-component gamma mixture prior (DuMouchel 1999):

The five hyperparameters (α₁, β₁, α₂, β₂, P) are fitted by empirical-Bayes maximum likelihood over all PTs of the brand in that cumulative run:

def fit_mgps_prior(n, E, seed=0):
    """Fit the 5-parameter two-gamma mixture prior from (n, E) pairs."""
    # ... maximise mixture negative-binomial marginal likelihood ...
    a1, b1, a2, b2 = np.exp(best[:4])
    p = float(sigmoid(best[4]))
    if a1 / b1 > a2 / b2:  # order components by prior mean
        a1, b1, a2, b2, p = a2, b2, a1, b1, 1 - p
    return (a1, b1, a2, b2, p)

Score: posterior geometric mean

After the prior is fit, each PT gets posterior component weights w_i1, w_i2 from the negative-binomial marginals, and

where ψ is the digamma function. EB05 / EB95 are the 5th / 95th percentiles of the same posterior mixture (found by bisection on the mixture CDF).

def mgps_ebgm(n, E, prior):
    a1, b1, a2, b2, p = prior
    w1, w2 = _posterior_weights(n, E, prior)
    return np.exp(w1 * (digamma(a1 + n) - np.log(b1 + E))
                  + w2 * (digamma(a2 + n) - np.log(b2 + E)))

How it is called for one brand in one cumulative month:

n_arr = np.array([rows[p]["a"] for p in pts_d], float)
e_arr = np.array([max(rows[p]["expected"], 1e-9) for p in pts_d], float)
prior = fit_mgps_prior(n_arr, e_arr)
ebgm = mgps_ebgm(n_arr, e_arr, prior)
eb05 = mgps_quantile(n_arr, e_arr, prior, 0.05)

RGPS → ERAM (three steps)

RGPS is run once per preferred term, across all vaccine products (type × manufacturer) in the country’s reports. Pfizer and Moderna are then scored as target products of that fit.

Step 1a — group age/sex strata by event rate

Stratum event rates are shrunk with a mean-1 gamma prior Gamma(g₀, g₀), ordered, and pooled until each group has at least max(10, √N₊/2) events:

def grouped_strata(stratum_events, stratum_reports):
    p0 = total_events / ns.sum()
    g0 = estimate_mean1_gamma(Ns, p0 * ns)
    lam = (Ns + g0) / (es + g0)          # shrunk stratum rates
    n_min = max(10.0, np.sqrt(total_events) / 2.0)
    # walk strata from highest lam → lowest, open a new group when n_min is hit

Step 1b — select masking product predictors

Other products whose grouped-strata-adjusted rates look extreme are selected as regression predictors. With f_j the adjusted expected count and δ₀ from another mean-1 gamma fit:

selected = np.where(
    (np.abs(Nj - fj) > 2.5)
    & ((q01 > mu0) | (q99 < mu0))
)[0]

Step 1c — penalised logistic regression

An empirical-Bayes logistic regression is fit with grouped-stratum intercepts and coefficients for the selected products (other products collapsed into an OTHER category). The site uses the ordinary logistic link (α = 1/2 in the white paper’s extended-logistic family — Empirica’s α-optimisation is unpublished):

beta0 = _elr_newton(cell_g, cell_cat, cell_n, cell_y,
                    n_groups, len(selected), p_group, g0, delta0, p0)
baseline = sigmoid(beta0)   # per grouped stratum

Step 2 — regression-adjusted expected count

For the target brand, the RGPS expected count is the sum of the fitted baseline probabilities over that brand’s reports in each grouped stratum:

t_grouped = np.asarray(t_counts, float) @ Gmat   # target reports → groups
E_target = float(t_grouped @ baseline)

This E_RGPS is the “Expected (RGPS)” column on the Top signals table. When other products soak up much of the event, E_RGPS falls below E_MGPS and ERAM rises.

Step 3 — single-gamma posterior mean (ERAM)

A single Gamma(γ, δ) prior is fit to the (N_j, E_j) pairs of all products for this event. The score is the posterior arithmetic mean:

gamma_a, gamma_b = fit_single_gamma(Nj, Ej_products)  # across products
shape = N_target + gamma_a
rate  = E_target + gamma_b
ERAM  = shape / rate
ER05  = gamma_dist.ppf(0.05, shape, scale=1.0 / rate)

Side-by-side comparison

Why ERAM can exceed EBGM (“masking”)

The heatmaps define

so +100% means ERAM is twice EBGM. That typically happens when Step 2 of RGPS produces a much smaller expected count than MGPS’s stratified E: other products (often other COVID brands, or strongly associated non-COVID vaccines) were selected as predictors and “claimed” part of the background rate. The same observed count a over a smaller E_RGPS yields a larger posterior mean.

Conversely, a negative masking effect means the regression raised the expected count for the target brand (the event looks more common in that brand’s strata / product neighbourhood than the plain stratified rate suggested).

Click any tile on a heatmap page to open the masking explainer — it shows E_MGPS, E_RGPS, the fold change between them, and how many predictors RGPS selected for that PT.

Worked example: Myocarditis · Pfizer · Nov 2022

Everything below is taken from the live US pipeline for the page Pfizer — Nov 2022 (cumulative from Jan 2021), preferred term Myocarditis. No symbols only — each step shows the actual counts that produced EBGM ≈ 1.30 and ERAM ≈ 2.55 on that page.

Step 0 — contingency table (US reports, Jan 2021 → Nov 2022)

Crude (unstratified) expected count would be n_drug × n_event / N = 433,449 × 3,096 / 1,005,816 = 1,334.2. MGPS and RGPS do not use that crude figure.

Step 1 — MGPS expected count (age × sex stratified)

Across the 27 Harpaz sex × age strata, each stratum contributes drug_s × evS / N_s. The largest contributions for Myocarditis are young males:

Step 2 — MGPS prior + EBGM

A two-gamma mixture prior is fit by empirical Bayes over all Pfizer PTs in this cumulative run (not just Myocarditis). The fitted hyperparameters are:

For Myocarditis specifically, the posterior almost entirely prefers the second (near-1) component:

Step 3 — RGPS: which products are “maskers”?

RGPS rebuilds the expected count from a product × stratum table (every vaccine type × manufacturer in the US data). Strata are collapsed into G = 12 rate groups (shrinkage strength g₀ = 0.43). Products whose grouped-strata-adjusted counts look extreme are selected as regression predictors. For Myocarditis, 22 products are selected. The biggest |N − f| gaps:

Janssen is a classic masker here: it has far fewer Myocarditis reports than its age/sex mix would predict (105 vs 277), so the regression attributes part of the background rate to other products’ absence of signal — which lowers Pfizer’s baseline.

Step 4 — RGPS regression-adjusted expected count

After fitting the penalised logistic regression (grouped-stratum intercepts + the 22 selected product coefficients), Pfizer’s expected count is the sum of the fitted baseline probabilities over Pfizer’s reports:

Same observed count (1,932), but the “null” background that Pfizer is compared against has almost halved once masking products are accounted for.

Step 5 — RGPS single-gamma score (ERAM)

A single Gamma(γ, δ) prior is fit across all products’ (N_j, E_j) pairs for Myocarditis. The fit gives γ = 0.270, δ = 0.297. Then:

Step 6 — put them next to each other

So for this PT, PRR says “about 2×”, MGPS says “after age/sex stratification and shrinkage, only about 1.3×”, and RGPS says “after also removing the masking of other products (especially Janssen’s deficit and Moderna’s share), about 2.5×”. That +96% jump is exactly the masking effect shown when you click the Myocarditis tile on the Pfizer Nov 2022 heatmap.

References

• Harpaz R, DuMouchel W, et al. (2022). Signaling COVID-19 Vaccine Adverse Events. Drug Safety 45:765–780.
• DuMouchel W (1999). Bayesian Data Mining in Large Frequency Tables. The American Statistician 53:177–190. (MGPS)
• DuMouchel W & Harpaz R (2012). Regression-Adjusted GPS Algorithm (RGPS). Oracle Empirica white paper.

Adaptations on this site vs. the article’s Empirica runs: each score is a cumulative country-only database run (not all-VAERS); products are VAERS vaccine type × manufacturer; each report is assigned its priority vaccine (Pfizer > Moderna > other); the RGPS link is ordinary logistic (α = 1/2). Scores are capped at 5 for colour on the treemaps only — the stored parquet values and table columns are uncapped.