Replicating a Published Multinomial NMA: Response Rates in Relapsed/Refractory Myeloma

A parameter-level replication of van Beurden-Tan et al. 2022 — and what it took to make the published numbers appear

Author

Xiaoge Zhang, PhD

Published

August 26, 2026

WarningWhat this is

A replication and appraisal of a published analysis, not a new clinical finding. Nothing here should be read as evidence about how these treatments compare; everything here is evidence about how one published comparison behaves when you rebuild it.

Every input comes from one openly licensed paper:

van Beurden-Tan CHY, Sonneveld P, Uyl-de Groot CA. Multinomial network meta-analysis using response rates: relapsed/refractory multiple myeloma treatment rankings differ depending on the choice of outcome. BMC Cancer. 2022;22:591. doi:10.1186/s12885-022-09571-8 — CC-BY 4.0.

The model is a fixed-effect Bayesian multinomial logistic NMA over 17 phase III randomised trials and 16 treatments, ported from the paper’s WinBUGS code to Stan, and cross-checked by running the original code in JAGS.

What this replication found

Finding How it was established
1 The published results reproduce exactly — 16 response rates to within 0.27 percentage points, SUCRA to within 0.2 Independent Stan implementation, validated against every direct comparison in the network
2 But not from the data file the paper supplies. Appendix A’s data block holds the uncorrected counts; running it gives the leading treatment 4.0% against a published 42% The Methods state a zero-correction the appendix data does not contain; applying it makes the results appear
3 The headline is a 47% claim, not a ranking. The treatment reported as best holds 47.1% of the posterior probability of being best; the runner-up holds 38.5% Rank probabilities computed under the paper’s own priors and correction
4 The winner changes identity under an undiscussed constant and under a conventional prior Sensitivity sweep over the zero-correction constant and the prior width
NoteVerdict

The analysis is sound and it replicates. Its conclusion that treatment rankings depend on the choice of response outcome is real, and this rebuild reproduces it.

What does not survive is the strength of the ranking. Under the paper’s own settings, the three leading treatments hold 93% of the probability of being best on complete response between them, without separating from each other — and the identity of the leader flips when a continuity-correction constant nobody justified moves from 1 to 3, or when the prior is narrowed to a width that most current practice would call conservative.

The network is precise about which treatments are worst and vague about which is best. Decisions are made at the end it cannot resolve.

The question

A conventional network meta-analysis of response splits patients in two: responders and everyone else. But myeloma response is not binary. It is graded — complete response, partial response, less than partial — and collapsing that grading throws away the distinction clinicians actually use.

A multinomial NMA keeps all three categories in one likelihood. The paper’s claim is that this matters: rank the treatments by complete response and one regimen leads; rank them by objective response, which pools complete and partial together, and a different one does.

That claim is correct. The question this page asks is a different one: how much of the ranking is carried by the data, and how much by choices made along the way?

The evidence network

Figure 1: The 17 trials as a network of 16 treatments. Dexamethasone, the reference treatment, is highlighted. Edge width counts trials; the doubled edge is the two lenalidomide trials, MM-009 and MM-010. Dashed edges carry no information about complete response, because both arms report zero complete responders.

Two structural facts matter later.

Twelve of the sixteen treatments appear in exactly one trial. Their estimates are anchored through a single comparison, and everything else about them is inferred indirectly.

The network contains one closed loop — dexamethasone to bortezomib to thalidomide and back — so inconsistency between direct and indirect evidence could be tested in exactly one place. The paper does not test it. That is not unusual, but a loop that exists and is never examined is worth naming.

What the ranking actually looks like

Figure 2: Posterior probability that each treatment occupies each rank, under the paper’s own priors and its own zero-correction. Treatments are in one order in both panels so movement between outcomes is visible; the boxed column is rank 1. SUCRA scores are printed on the right.

Read the boxed column, not the SUCRA scores.

On complete response, the leading treatment holds 47% of the probability of being best and the runner-up 39%. The two are separated by one SUCRA point — 93 against 92 — which is not a distinction any decision should rest on. The third holds 7%, and between them the top three account for 93% of the probability of being best without resolving which of the three it is.

On objective response the same three reorder, and elotuzumab plus lenalidomide and dexamethasone moves from SUCRA 31 to SUCRA 74 — twelfth to fourth. That movement has a mechanism: in ELOQUENT-2 elotuzumab lowered complete response (24/325 against 14/321, odds ratio 0.57) while raising overall response. The same drug ranks near the bottom or near the top depending only on which response definition the committee is looking at.

Meanwhile the bottom of the network is sharp. Dexamethasone, oblimersen and thalidomide each place around half their probability mass on a single rank. The model knows what does not work.

The replication

Every one of the sixteen complete response rates, both credible interval bounds included, against the published Figure 2.

Mean absolute error 0.27 percentage points on the rates, 0.2 points on SUCRA. The fourth column is what the same model produces when fed the data file the paper publishes for replication: mean absolute error 19.6 points.
Rank Treatment Published This rebuild From the Appendix A data SUCRA ours / published
1 PomBorDex 42% [17, 71] 41.9% [16.9, 70.6] 4.0% 93 / 93
2 CarLenDex 40% [21, 61] 40.4% [20.6, 61.0] 3.2% 92 / 92
3 DaraLenDex 34% [17, 55] 34.5% [16.6, 54.9] 2.6% 85 / 84
4 DaraBorDex 27% [10, 53] 26.9% [9.7, 52.2] 2.2% 75 / 76
5 CarDex 24% [9, 48] 24.0% [8.7, 47.2] 1.9% 70 / 70
6 IxaLenDex 21% [8, 39] 20.9% [8.4, 38.1] 1.3% 64 / 64
7 PanoBorDex 20% [7, 40] 19.5% [7.2, 39.2] 1.5% 59 / 60
8 PLDBor 19% [5, 45] 18.7% [5.1, 43.6] 1.4% 57 / 57
9 LenDex 11% [5, 21] 11.5% [5.0, 20.5] 0.7% 42 / 42
10 BorThalDex 11% [2, 33] 10.6% [1.8, 33.3] 0.5% 39 / 39
11 Bor/BorDex 11% [4, 24] 10.8% [3.8, 22.8] 0.7% 37 / 37
12 EloLenDex 8% [3, 18] 8.1% [2.6, 17.9] 0.5% 31 / 30
13 PomDex 6% [0, 29] 5.8% [0.4, 31.7] 58.8% 22 / 22
14 Thal/ThalDex 3% [1, 11] 3.2% [0.6, 10.8] 0.1% 14 / 14
15 OblDex 4% [0, 23] 3.6% [0.0, 22.5] 13.6% 13 / 13
16 Dex 1% [1, 2] 1.3% [0.6, 2.2] 0.1% 6 / 6

Why the appendix data does not work

Appendix A is titled WinBUGS code, init and data files. Its data block holds the GMY302 row as 0, 19, 95 against n = 114. Those sum exactly, so no correction has been applied.

The Methods say one was:

“In case there were zero responders in at least one category within an RCT, a zero-correction factor of k = 1 was added to all the fields in the data table of that specific trial to properly run the NMA.”

Three trials qualify. Add the correction and the published numbers appear. Leave it out and the model does something specific and instructive rather than merely wrong: the study baselines for those trials become unidentified — the log-odds of partial versus complete response is infinite when there are no complete responses — and the pooled reference baseline, which the code takes as a plain average over the six trials containing the reference treatment, follows the two runaway values. The reference response rate collapses to 0.00003% and carries the whole network down with it.

The same data block lists GMY302 twice, with different counts — 10, 19, 95 in the first row, which exceeds its own n by ten patients, and 0, 19, 95 in the last. Neither is the corrected value the analysis used.

This was confirmed two ways: an independent Stan implementation, and the paper’s own WinBUGS code run in JAGS with its stated three chains, 25,000 burn-in and 80,000 iterations. Both engines agree with each other and disagree with the published figures by an order of magnitude. In the JAGS run the Gelman-Rubin statistic — the diagnostic the paper reports as showing convergence — comes back at 1.74.

What the ranking rests on

The zero-correction constant sets the winner. The paper adds 1. Nothing in the paper explains why 1, and no sensitivity analysis varies it.

k = 1 (published) k = 2 k = 3
PomBorDex, probability best 47.1% 45.1% 43.8%
CarLenDex, probability best 38.5% 42.4% 44.2%
Leading response rate 41.9% 46.8% 49.2%
BorThalDex SUCRA 39 45 48

NICE DSU TSD 2 recommends a different convention — add 1 to the denominator and 0.5 to the numerator — and separately notes that trials with zero cells in both arms contribute no evidence about the treatment effect and can be excluded. The half-patient version cannot be expressed here at all: a multinomial likelihood takes integer counts. That is a real gap in the method’s inherited practice, not an error by these authors.

The prior matters more than the correction. The paper’s only statement about its priors is the phrase # vague priors in a code comment. Replacing dnorm(0, 0.001) with normal(0, 2.5) — still wide enough to admit odds ratios of 130, far beyond anything in this network — moves the leader at every value of k tested:

vague (published) normal(0, 2.5)
PomBorDex, probability best 47.1% 21.7%
CarLenDex, probability best 38.5% 63.5%

One treatment’s estimate is generated entirely by the correction. Both arms of GMY302 report zero complete responses. Oblimersen’s published complete response rate of 4% [0, 23] therefore rests on the single patient the correction invents — and that rate occupies a place in a sixteen-treatment ranking, which shifts every other treatment’s rank distribution.

How the replication was done

Six steps, in order. The same sequence works on any published Bayesian analysis.

Step What it catches
1 Parse the trial data from the source HTML, not the XML or the PDF Table cells that hold two arms’ values; the XML concatenates “5” and “41” into “541”
2 Cross-check every count against a second published rendering Transcription error — here, the appendix’s own data block, which also reveals the treatment node mapping
3 Write down the direct evidence before fitting anything A model that disagrees with the trials it contains, whatever the paper reports
4 Report posterior contraction, not only convergence Parameters where the posterior is still the prior, which R-hat cannot see
5 Run the original code in its own dialect Confusing a porting error with a finding
6 Vary every undocumented constant, not only the ones the paper varied Conclusions that rest on an arbitrary choice

How it is put together

The data is parsed from the paper’s HTML into arm-level counts, cross-checked against the WinBUGS data block, and never edited by hand. The model is a single Stan file of about thirty lines whose model block is four; the multinomial likelihood over three ordered categories is a softmax with the reference category’s logit pinned at zero, which is exactly what the four WinBUGS lines it replaces compute. The zero-correction, the prior widths and the set of trials entering the pooled baseline are all inputs to the model rather than constants inside it, so every table above is the same code with different arguments.

Reproduction is Rscript over numbered files: parse, cross-check, direct evidence, network, fit, diagnose, replicate, sweep, plot. The source PDF and appendices are not redistributed here; the paper is open access under CC-BY 4.0 and everything needed to run this is in it.

Before the multinomial model, the same network was fitted as a conventional binomial NMA on complete response alone — collapsing partial and less-than-partial into one category. It reproduces every direct comparison in the network: APEX 9.30 crude against 10.38 estimated, ASPIRE 4.53 against 4.57, the two lenalidomide trials pooling from 18.50 and 5.96 to 9.49.

It also breaks in a way the multinomial model does not. Posterior contraction — one minus the ratio of posterior to prior standard deviation — separates the parameters the data determines from the ones it does not:

binomial multinomial
Oblimersen + Dex 0.13 0.37
Pomalidomide + Dex 0.13 0.53
Well-identified treatments 0.99 0.72 – 0.89

The two treatments whose trials report no complete responses are unidentified in the binomial model and identified in the multinomial one, because the partial and less-than-partial counts carry information the complete-response count alone does not. This is the strongest argument for the paper’s method, and the paper does not make it. The trade is visible in the third row: every other treatment loses precision, because each now spends two parameters instead of one.

A note on diagnostics from the broken version, which is the most transferable thing in this exercise. Under vague priors the binomial fit returned R-hat 1.002 and a smallest effective sample size of 3,714 — a textbook pass — while assigning 98.7% of the probability of being best to two treatments whose posteriors were their priors. Four chains wandering the same flat ridge agree with each other perfectly. R-hat tests whether chains agree, not whether the posterior means anything.

Narrowing the prior does not fix that; it hides it. Oblimersen’s interval goes from [0, 100] to [0.2, 7.0] as the prior tightens, while its contraction falls from 0.19 to 0.07. The interval reads better precisely as it carries less data.

Limitations

JAGS is not WinBUGS. The original was run in WinBUGS, which is no longer practical to run here. JAGS shares the BUGS language and a Gibbs-family sampler and is the closest available substitute, but a residual difference in sampler behaviour cannot be ruled out. The conclusion that the appendix data does not produce the published results rests on two independent engines agreeing, not on either one being WinBUGS.

The zero-correction reading is inferred, not confirmed. The Methods sentence is unambiguous about the constant but not about the fields. Adding 1 to all three categories in both arms of the three affected trials reproduces the published figures to 0.27 percentage points, which is strong evidence that this is what was done, but it is evidence rather than documentation.

The appraisal is of one analysis, not of a literature. Zero cells, undeclared continuity corrections and rankings reported without their rank probabilities are common. Nothing here establishes how common.

Nothing here is a clinical conclusion. The relative effects estimated in this page are the paper’s, recovered; they inherit every limitation of a fixed-effect model over 17 heterogeneous trials, including the pooling of bortezomib with bortezomib plus dexamethasone into one node, the pooling of three thalidomide dose arms from a four-arm trial into one, and the absence of any between-study variance term.

Scope held to the source

Fixed effect, seventeen trials, sixteen treatments, dexamethasone as reference, three ordered response categories. No random-effects comparison, no node-splitting, no meta-regression — the paper runs none of these, and adding them would answer a different question than “does this replicate”.

Reference

van Beurden-Tan CHY, Sonneveld P, Uyl-de Groot CA. Multinomial network meta-analysis using response rates: relapsed/refractory multiple myeloma treatment rankings differ depending on the choice of outcome. BMC Cancer. 2022;22:591. doi:10.1186/s12885-022-09571-8. Published under CC-BY 4.0; trial counts and published estimates are reproduced here under that licence.

Methodological standards referenced: NICE DSU Technical Support Document 2 (generalised linear models for network meta-analysis), from which the paper’s competing-risk multinomial model is taken.