Compliance with Ethics Guidelines

The original trials used in this analysis obtained patient consent to participate and were performed in accordance with the Helsinki Declaration of 1964. All ethical approvals and permission to access and use non-public data have been granted for this study. All remaining data are publicly available.

Data Sources and Extraction

RWD from patients with FL was obtained from the Haematological Malignancy Research Network (HMRN) [17], International Organization for Medical Oncology (iOMEDICO) [18], SCHOLAR-5 [19], Lymphoma Epidemiology of Outcomes (LEO) [20], Retrospective Cohort Study of Treatment Outcomes of Adult Patients with Relapsed or Refractory Low-Grade Follicular Lymphoma (ReCORD-FL) [21], Batlevi et al. [3], and National LymphoCare Study (NLCS) [22]. Data from LEO (number of patients [n] starting 3L = 441), SCHOLAR-5 (n = 143) and NLCS (n = 401) were available from the start of index therapy. Data were available by line (L) of therapy but considered as the same pool of patients (e.g., patients receiving 4L were repeated from the pool of patients receiving 3L) from the Tumor Registry of Lymphatic Neoplasia (TLNext; ClinicalTrials.gov, NCT00889798) conducted by iOMEDICO (3L-8L, n = 69), ReCORD-FL (3L-6L, n = 187), and Batlevi et al. (3L-6L, n = 299). HMRN data (n = 82) were only reported for 3L. Data were obtained from the UK (HMRN), USA (LEO, NLCS, and Batlevi et al.), Germany (TLNext), or a combination of these countries (ReCORD-FL and SCHOLAR-5).

Individual patient data for OS from patients receiving mosunetuzumab were obtained from the expansion cohort of the GO27981 trial (n = 90) for patients with RR FL receiving 3L + therapy. A total of 8 OS events occurred among 90 evaluable patients in the August 2021 data cut used in this analysis.

Data from all sources were digitized using standard software and the Guyot method [23], except for NLCS, from which individual patient data were available (Fig. 1). Patients in these datasets who started 3L therapy could start subsequent lines of therapy and appear in multiple lines. Therefore, as patients could not be individually identified, they could be counted multiple times. If left unadjusted, this would artificially inflate the number of unique patients. As a result, a resampling with replacement approach was taken to reduce the impact of this artifact by mixing different lines of therapy for the digitized data sets. For example, if 30 unique patients were tracked from 3L in the database (and no new patients entered at later lines), and there were 30 observations in the 3L, 20 in the 4L, and 10 in the 5L in the data, resampling was performed from the pooled 60 observations to have a total of 30 patients while keeping 50% of patients in the 3L, 33% in the 4L, and 17% in the 5L. This approach allows for adjustment of the number of observations for an artificial inflation in precision while still taking into account the evolution of the patients in subsequent lines of therapy. We implicitly assumed that the incidence of patients entering at each line in the data would be the same as the incidence we observe when we track patients from 3L progressing to subsequent lines over time. No adjustments in the population characteristics of the data have been made because no individual patient data were available from the majority of the RWD sources. This could be considered as a relevant limitation, but the absence of heterogeneity adjustments (e.g., exclusion criteria, covariate adjustment, etc.) to the sources of RWD allows a wide range of survival estimates to be considered in this setting. This also allowed the patient characteristics of a given trial to be covered by the umbrella of all the studies selected. As a result, this approach can be useful to inform the extrapolations even with dispersed and heterogeneous datasets.

Extracting Survival

All analyses were performed using the R statistical programming language. Some publications showed the number of patients starting each line of therapy but displayed the survival curves by mixing patients starting a new line with patients already tracked from previous lines. Similar to the resampling approach described in the previous subsection, if the number of patients starting each line of therapy was available, the number of patients starting at each line of therapy was matched to avoid artificial inflation of observations for the time to event data. The analyses to obtain smoothed hazards and estimate disease-wide OS relied on the same resampling procedure. KM curves were estimated during the resampling process (1000 resamples), and pointwise OS was obtained together with the corresponding uncertainty at 3-month intervals for ≤ 8 years for each RWD source.

For robustness, results were also obtained for patients receiving only 3L therapy. Because counting patients multiple times does not affect 3L therapy patients, the 3L-only group was used to test whether patterns observed in the 3L + group were due to patients being included multiple times in the pooled sample.

Smoothed Hazard

Smoothed hazards were obtained through the muhaz function from the muhaz package [24] in R by using kernel-based methods (R Foundation, Indianapolis, IN). The Epanechnikov boundary kernel function [25] was applied, with local bandwidths defined as by Mueller and Wang (1994) [26]. Smoothed hazards were obtained for KM curves generated for each bootstrapped sample. The same resampling method used to obtain KM curves without artificially increasing precision (as mentioned in the previous subsection) was used to calculate the smoothed hazards of each RWD source. Observations were resampled following this resampling method 1000 times. At each resample, the smoothed hazard was estimated for each RWD source. The final smoothed hazard for each RWD source was constructed from the 1000 iterations by averaging at each time point for which the smoothed hazard value was calculated, and the confidence intervals were obtained from the 95% quantiles. Because resampling produced some samples with few events at longer follow-up time points and no hazards could be estimated at those time points, smoothed hazards are shown through 6 years. Hazard trends were assessed for each average smoothed hazard to identify hazard patterns compatible with well-known survival distributions such as exponential, Gompertz, Weibull, log-normal, and log-logistic distributions (Fig. 2).

Fig. 2

Smoothed nonparametric hazard estimates using Epanechnikov boundary kernel function for patients with FL receiving a 3L + therapy or b 3L therapy. 3L third line; 3L + , third line or later; FL follicular lymphoma; HMRN Haemetological Malignancy Research Network; LEO Lymphoma Epidemiology of Outcomes; NLCS National LymphoCare Study; ReCORD-FL Retrospective Cohort Study of Treatment Outcomes of Adult Patients with Relapsed or Refractory Low-Grade Follicular Lymphoma; TLNext Tumor Registry of Lymphatic Neoplasia

Bayesian Random-Effects Meta-Analysis

Several different approaches could be taken to integrate the previously extracted information in a Bayesian context, like assuming the data comes from a specific distribution (e.g., Weibull) and meta-analyzing those parameters [27]. Since we would like to compare specific distributions against the traditional extrapolations that would be performed in a HTA submission, we used an alternative approach already used to elicit beliefs from experts [28, 29], and we applied it to RWD through pointwise survival. To do so, pointwise survival was predicted every 3 months over ≤ 8 years from the KM curves generated at each resample. Because some samples did not have sufficient numbers of events at longer follow-up periods, survival was not estimated past 8 years. Survival at each 3-month time point was averaged across 1000 resamples for each RWD source to obtain average survival and 95% CIs.

Pointwise survival was synthesized using a Bayesian random-effects meta-analysis with flat priors, adapting the code from the National Institute for Health and Care Excellence Decision Support Unit Technical Support Document 20 [30], Appendices A.1.2 and B.5, for a univariate case. Flat priors were used at this stage as no external information could be used to inform the parameters. As pointwise survival at any point in time must be contained between 0 and 1, the flat prior for the random effects must be set so that the outcomes verify this restriction, in this case as a uniform (0, 0.5). The true outcomes were predicted with random effects included, as in the Technical Support Document 20. True outcomes were used because the focus of this analysis is the value for the overall population and not for a specific study. That is, the aim is to obtain long-term extrapolations in a generic population informed by different RWD studies. Pointwise survival at any time point (done here at 5 and 8 years) was assumed to follow a log-normal distribution, following the same principles applied when calculating the pointwise confidence intervals for a KM estimator [31]. Overall disease-wide pointwise survival and the corresponding uncertainty were estimated at 5 and 8 years. For robustness, the same analysis was again performed for patients receiving 3L therapy.

Bayesian Extrapolations to Mosunetuzumab

Current data from the GO29781 trial present challenges in terms of long-term extrapolation due to few OS events and short follow-up periods. To increase the long-term credibility of these data, the results from the meta-analyzed pointwise survival at 5 and 8 years were used as informative priors in a Bayesian extrapolation. The survival function for each statistical distribution (e.g., exponential, Weibull, etc.) were used at each time point by setting the survival function with the parameter for time and survival fixed at specific values. The equation is then solved for the scale or rate parameter (depending on the distribution) as a function of the specific survival value (e.g., 0.5), time point (e.g., 8 years), and remaining parameters if any (e.g., in the case of the exponential distribution there would be none as it only has 1 parameter). The derived analysis is shown in Table S1 in the electronic Supplementary material. It is assumed that the log of these parameters is distributed as a normal distribution, as is commonly done in widely used statistical packages [32]. As such, the lower and upper 95% survival values for OS at 8 years obtained from the meta-analysis can be used to calculate the mean and variance of these parameters. Note that the informative prior is not set over the statistical distributions themselves but rather on one of their parameters (rate/scale). An example applied to the exponential distribution is shown in Table S1 in the electronic Supplementary material. The remaining parameters (if any) were set as a normal distribution with large uncertainty (precision = 0.0001) while ensuring that they stayed within valid ranges. An exception to this method was the Gompertz distribution, where the shape parameter takes the frequentist coefficient as a prior to ensure convergence.

Comparison of Bayesian and Frequentist Extrapolations

Frequentist extrapolations using the individual patient data from the GO29781 trial alone were computed as a benchmark with which Bayesian estimates were compared in terms of outcomes and uncertainty. Both frequentist and Bayesian estimates were adjusted by assuming that weekly hazard could not be smaller than general-population hazard. This method was chosen because it could be applied to both frequentist and Bayesian estimates and because other preferred methods like internal additive hazards [33] did not converge properly with the limited number of events in the data. The method described by Guyot et al. [34] was not used here because it constrains distributions with decreasing hazards like the log-normal distribution to have very low survival curves in order to conform with the long-term survival in the general population and, therefore, is better suited for more flexible hazards such as those obtained from extrapolation methods like splines. These methods were not used here as the focus was on parametric distributions, which are also more commonly used by different HTA evaluations. Weekly mortality hazard for the general population was computed using a heterogeneous cohort (n = 90) matched by age and sex to the cohort from the GO29781 trial. Because age varied across the cohort of patients, the average age of the heterogeneous cohort grew more slowly over time than that of the homogeneous cohort due to older patients dying at a higher rate than younger patients, which translated into a longer tail in the survival distribution. Estimation was performed in the rjags package [35] using 5000 iterations, 3 chains, and the n.adapt parameter set to 500. Extrapolations were performed using exponential, Weibull, log-normal, log-logistic, and Gompertz distributions (Table S1 in the electronic Supplementary material). Apart from the smoothed hazards, no other method, such as the marginal likelihood method, was used to evaluate the best-fitting distribution in the extrapolations because of the limited number of events observed for mosunetuzumab.

Bayesian extrapolations were compared with frequentist extrapolations and with evidence from the RWD source using smoothed hazards. An overall assessment of mean and median survival, long-term survival, reduction of area under the curve, and reduction in the coefficient of variation (CV) of estimated total life years was performed. Reduction in uncertainty was measured by area under the curve between 0 and 30 years, as well as by a reduction in CV.