Let Y j , j = 1 , … , J, denote a safety response of interest in the jth study, which are collected in individual trials but may not be fully reported. For censored outcomes, the censoring mechanism can be defined by bounding variables ( A , B ), with semi-closed boundaries ( A j , B j ] for response variable Y j . Here, both bounding variables could be covariate-dependent. Often, we assume that the outcome model and the censoring mechanism are independent, so the censoring mechanism does not contribute to the inference and model estimation. It is a fundamental assumption for censored data behind most statistical methodologies [58]. In survival analysis, it is also known as noninformative censoring. Assume the random variable Y has the right continuous cumulative distribution function F Y ( y ) = P [ Y ≤ y ]. Denote f Y ( y ) the probability density/mass function of F Y ( y ).

In a meta-analysis of safety data, the frequency of adverse events may not always be reported. For example, left censoring occurs when some severe (grade 3 or higher) events are not observed due to low incidence. In this case, the cutoff boundaries are not random but fixed and study-specific, which spontaneously satisfies the noninformative censoring assumption. Denote the fixed cutoff by c j for each study. The number of subjects having a specific event in the jth study Y j follows a binomial distribution with study-level sample size n j and AE incidence probability θ j Y j ∼ B i n ( n j , θ j ) . Therefore, the likelihood function for both fully observed and censored events can be written by: (2.1) L = ∏ o = 1 O f Y ( y o ) ∏ l = 1 L [ F Y ( c l ) − 0 ] ∏ r = 1 R [ F Y ( n r ) − F Y ( c r ) ] = ∏ o = 1 O f Y ( y o ) ∏ l = 1 L F Y ( c l ) ∏ r = 1 R [ 1 − F Y ( c r ) ] = ∏ o = 1 O f Y ( y o ) ∏ l = 1 L ∑ k l = 0 c l f Y ( k l ) ∏ r = 1 R [ 1 − ∑ k r = 0 c r f Y ( k r ) ] , where O is the number of fully-observed AE outcomes, L the number of left-censored outcomes, and R the number of right-censored outcomes. c l and c r are cutoff values for left-censored and right-censored data, respectively.

The likelihood (2.1) can simultaneously handle both cases of left- and right-censoring for meta-analysis of safety data. In meta-analysis of the grade 3-5 AEs, we often only observed left censoring due to reporting cutoffs. In meta-analysis of the all-grade AEs, right censoring may also occur when some studies report grade 2 or higher AEs only rather than all-grade AEs [4]. If Y j is left-censored, then Y j lies in the semi-closed interval ( A l = 0 − , B l = c l ], where c l is the corresponding cutoff value for left-censored data. If Y j is right-censored, then Y j lies in the semi-closed interval ( A r = c r , B r = n r ], where c r is the corresponding cutoff value for right-censored data and n r is the total number of subjects in the rth study. This is specified to be consistent with the JAGS model implementation in Section 2.2. Compared with multiple imputation approaches that often requires the assumption of normality and involves intensive computation for Bayesian modeling [26], the proposed approaches are computationally efficient to model incomplete data.

The incidence probability can be decomposed through the link function g ( θ ), (2.2) g ( θ j ) = logit ( θ j ) = μ + α j + X j β , where X j is a design matrix for study-level covariates from the j th study, and we consider the logit link as the default link function. Given the overall mean μ and the study covariate effects β associated with X j , we assume the study-specific random effect α j is conditionally independent with mean 0 and variance σ α 2 .

For nonhierarchical intercept μ or slope parameters in β in the linear term of logistic model, a symmetric Cauchy distribution with center 0 and scale 2.5 was suggested [16]. For random effects in the linear term, which include the effects of various study-level factors and study-specific effects α j , we assume they follow mean-zero normal distribution with different variances, with the prior distributions of standard deviation parameters following half-Cauchy distribution and a finite scale parameter of 10 or 25 [14]. In practice, the normality assumption can also be replaced by using heavy-tailed distributions such as t-distribution. In the illustrative anti-PD-1/PD-L1 example, the meta-analysis included 125 studies with a total of 20,218 patients. To identify possible sources of heterogeneity between studies, the following pieces of study-level information were also extracted: trial name, number of treated patients, selected immunotherapy drug, dosing schedule, cancer type, AE frequency (which could be missing), and the study-specified censoring cutoff for AE reporting. In particular, the study-level factors could include the drug type or drug/dose combinations, and/or cancer type, such that it is possible to estimate the AE incidence of a new drug or a new cancer type in a future trial with Bayesian credible interval.

Lastly, the model specification of the proposed Bayesian modeling framework only involves discrete data in the likelihood function, while the prior distribution of parameters is proper. As a result, it is easy to verify the property of Bayesian posterior inference. Specifically, we can have the following remark:

Overall, the proposed Bayesian framework can help to streamline the modeling process of drug safety data and enhance the reporting reliability of meta-analysis.

To carry out posterior inference and assess the performance of the proposed Bayesian model, we apply Just Another Gibbs Sampling (JAGS) to generate samples from the posterior distribution. JAGS makes Bayesian hierarchical models easy to implement using Markov Chain Monte Carlo (MCMC) simulation [40] in R and other computational software. In the presence of censored data in the response variable, an existing function, known as dinterval distribution, is commonly used to model censored data [27, 43]. However, such model specification for censored data in JAGS intrinsically yields a mis-specified deviance function [42]. More technical details are provided in Qi et al. (2022) [45].

Alternatively, we apply a simple but effective approach to censored data specification. To facilitate model implementation for censored observations (when δ j = 0) and avoid the miscalculation of deviance via the dinterval() function in JAGS. Here, we utilize the idea of data augmentation by introducing ancillary indicator variables W j . Each W j separates the left-censored from right-censored observations ( W j = 1 if left-censored, 0 if right-censored). By assuming W j follows a Bernoulli distribution, we have the density function (2.3) f W ( w j ) = F Y ( c l ) w j [ 1 − F Y ( c r ) ] 1 − w j = ∑ k l = 0 c l f Y ( k l ) w j 1 − ∑ k r = 0 c r f Y ( k r ) 1 − w j , which exactly matches the second and third terms for censored observations in (2.1), with the cumulative binomial distribution of incidence probability of AE in the jth study, restricted by a pre-determined study-level cutoff value c j . This approach can also be extended to interval-censored data [45].

A JAGS model specification for the application is provided in the Supplementary Material for illustrative purposes. Together with fully observed data that follow a binomial distribution, the full likelihood implemented in a JAGS model is, in fact, identical to the exact likelihood of observed and censored cases in (2.1). This creates the correct focus of model parameters and produces the proper posterior samples, as well as simultaneously computes the correct deviance for model selection. For example, by calling the deviance module in JAGS, correct deviance information criterion (DIC) [49], penalized deviance [41], posterior averaging information criterion (PAIC) [60], or widely applicable information criterion (WAIC) [55, 53] can be conveniently derived to assess candidate models for Bayesian model selection. It is important and beneficial to select the best model, especially in the presence of complicated model features. Also, from the point of view of convenience, the sampling procedure can smoothly start with default initial values for the censored data. The standard diagnostic tools can be applied to assess the posterior convergences of the model [15].

The total number of studies for each probability level of AE incidence is fixed at 10 ( J = 10) to reflect the typical number of studies in a meta-analysis for subgroups. We also consider J = 25 and 50 as comparisons. The outcome of interest, the number of AEs for each study, is generated from a binomial distribution with the number of patients ( n = 100) and probability of events ( p = 0.2 , 0.05, and 0.01, respectively). The probability of incidence is designed to cover the typical range of the probability levels of AE incidence. The data generation incorporates heterogeneity between studies into the incidence probability of AE. To be specific, the coefficient of study effect was generated from a normal distribution at mean 0 and standard deviation of 0.2.

To assess the performance of the proposed model that incorporates both observed and censored data, we consider six scenarios: (1) no censoring; (2) a low percentage (20%) of left-censoring for all three incidence levels; (3) a moderate percentage (40%) of left-censoring; (4) a high percentage (60%) of left-censoring; (5) a low/moderate percentage (20%/40%) of right-censoring; and (6) a moderate percentage (40%) of left-censoring with a low percentage (20%) of right censoring. In Scenario 1, the number of events for all studies is fully observed. In Scenarios 2, 3 and 4, which include observations that could only be informatively left-censored for zero and low incidence counts to mimic real-world cases. Therefore, in Scenario 2, we treat the 20% of AE data with the lowest observed incidence as censored data and the 80% of AE data that have a higher incidence as observed data. Similarly, in Scenario 3, in order to stress test the robustness of estimation in a more extreme case of censoring, 40% of AE data with low incidence are treated as censored and the remaining 60% are treated as observed. Lastly, in Scenario 4, the top 40% of studies are treated as observed data, and the remaining 60% are treated as censored data. Accordingly, the study-specific left-censored cutoff value is selected from the sorted outcomes in a descending order based on the percentage of left-censoring after data generation. For example, in the scenario of 60% left-censoring, the cutoff value corresponds to the ( ( 1 − 60 % ) J + 1 ) th outcome among sorted values. However, the right-censoring in AE reporting is a different censoring mechanism (e.g., the reporting of grade 2 or higher AEs rather than all-grade AEs) and likely from the random selection of studies. Therefore, in Scenario 5, 20% of studies are randomly treated as right-censored; in Scenario 6, 40% of AE data with low incidence are treated as left-censored data first and then randomly selecting one third of the remaining 60% (e.g., 20% among all) of data as right-censored. The study-specific right-censored cutoff value is generated from a binomial distribution by assuming the probability of grade 2 or higher AE is 50% of all-grade AE incidence.

We compare the proposed model, MAGEC, with four other methods: the pooled estimation method after continuity correction (PEM) [52, 24], the normal approximation method (NAM) [38], the logistic regression method (LRM), and the normal approximation method with robust variance estimator (RVE) [34]. In MAGEC, following the recommendation of the weakly-informative prior distribution for logistic regression models [16], we assign a half-Cauchy prior distribution on the standard deviation of study effects. In PEM [52], we pool observations by incidence level and the 95% nominal confidence intervals (CIs) for AE incidence levels are calculated by the exact binomial test. In NAM, a standard method in practice [38], we use a normal likelihood procedure to estimate the incidence rate by taking the inverse logit of the observed logit incidence [18] of each incidence level weighted by its within-level variance. In LRM, we estimate the AE incidences by an exact method through fitting a generalized linear model with a logit link. In addition, we compare the performance of NAM with and without robust variance estimators [34]. Therefore, in RVE, instead of Fisher information, we implement the sandwich estimator of variance into NAM to improve the robustness of the statistical inference on incidence rates.

We assess model performance in terms of accuracy in the estimation of the AE incidences. More specifically, we compare the five methods using the following metrics: mean absolute deviation (MAD), root mean squared error (RMSE) of the point estimates (for MAGEC, we use the posterior median), and coverage probability (CP) of the 95% nominal confidence interval (for MAGEC, we use the 95% posterior credible intervals), of the incidence rate parameters of interest in each scenario.

The results are based on 10,000 simulated data sets. For each method, we repeated the same data generation procedure in order to be able to compare results across methods. Figure 1 gives boxplots for the point estimates of incidence probabilities by scenario, method and the total number of studies. CPs of three parameters of interest by scenario, method, and the total number of studies are displayed in bar charts in Figure 2. In Table 2, performance in terms of both MAD and RMSE of incidence rates based on the five methods are shown for the four data censoring scenarios. The table that presents two additional data censoring scenarios is available in the Supplementary Material.

When there is no censoring (Scenario 1), the proposed method (MAGEC) has CPs (Figure 2), MADs, and RMSEs (Table 2) for the incidence probabilities that are almost identical to those of PEM and LRM. Of the five methods compared, PEM can be considered the gold standard/benchmark for both interval and point estimation. Our results indicate that MAGEC is not inferior to PEM. They also indicate that the CP for each probability level of AE incidence obtained from NAM and RVE appeared to be unstable in estimating the incidence probabilities of less frequent events compared with the other methods. The point estimates of the incidence probabilities in both NAM and RVE are overestimated in Scenario 1. This finding is consistent with arguments mentioned in the normal approximation for rare events [7] and biased results of estimation for rare events using normal approximation [33].

When 20% of data are censored (Scenario 2), the proposed method (MAGEC) performs better than the others in estimating incidence probabilities; with an increasing total number of studies, its performance in Scenario 2 is comparable to that in Scenario 1. Because censored observations are ignored under the other four methods (PEM, NAM, LRM, and RVE), it is unsurprising that their point estimates for the incidence probabilities are overestimated and that their CPs in Scenario 2 are much lower than those in Scenario 1. In contrast, the performance of MAGEC in Scenario 2 is almost identical to its performance in Scenario 1 for both interval and point estimates, demonstrating its robustness of censored data.

In a more extreme scenario where 40% of data are censored (Scenario 3), the proposed method (MAGEC) performs well, with little information loss compared with Scenarios 1 and 2. However, all other estimators of AE incidence lead to inferior CP due to an increased percentage of censoring. The point estimations obtained from PEM, LRM, NAM, and RVE in Scenario 3 are more biased than those obtained from these methods in Scenario 2. Based on the MAD and RMSE, there are larger deviations from the true values of the incidence probabilities compared with those in Scenario 2. With an increasing total number of studies, their CPs in Scenario 3 are much lower that those in Scenario 1 and 2. When the percentage of censoring reaches 60%, the performance of MAGEC is still the best among all methods. When we incorporate right-censored AE data, MAGEC consistently provides the lowest errors across all settings. In Scenario 6 with censoring possibly on both sides, the results from four other methods are unreliable. For example, the observed errors sometimes are smaller in Scenario 6 (e.g., MAD = 0.015 in PEM) than in Scenario 3 or 5 for one-sided censored data alone (e.g., MAD = 0.033 or 0.020 in PEM), when combining two wrongs makes a right because the analysis of neglecting the left-censored or right-censored data tends to be biased in the opposite direction. However, as illustrated in Table B.1 in the Supplementary Material, the smaller error cannot be warranted across different sample size settings. Overall, the MAGEC model yields not only more stable and superior coverage, but also unbiased estimates of the incidence probabilities in all six scenarios especially with censored event data. The proposed method, handling both observed and censored data took, around 20 seconds in the absence of censoring and 50-70 seconds with an increasing percentage of censoring. In contrast, the other four methods, which only consider observe data, completed the task in less than one second.

Across all scenarios considered above, MAGEC outperforms the other four methods in estimating incidence probabilities when AEs have low incidence and when a high proportion of AEs are censored. Furthermore, the quality of an estimator can be measured by its efficiency, which is defined as the asymptotic variance of an estimator [8]. The larger the variance, the lower the efficiency of an estimator. Here, the asymptotic relative efficiency (RE), defined by the reciprocal of the ratio of the asymptotic variances of two unbiased estimators, is given to examine the amount of information loss from informative censoring. There is a 10%-20% inflation in the RMSE from 0% to 40% censored data, suggesting limited efficiency loss of MAGEC in incidence estimation due to censoring. In summary, the proposed method consistently achieves a reasonable performance in estimating the incidence rates at different probability levels of AE incidence in the presence of censored data.