 Research
 Open Access
 Published:
Performance of several types of betabinomial models in comparison to standard approaches for metaanalyses with very few studies
BMC Medical Research Methodology volume 22, Article number: 319 (2022)
Abstract
Background
Metaanalyses are used to summarise the results of several studies on a specific research question. Standard methods for metaanalyses, namely inverse variance random effects models, have unfavourable properties if only very few (2 – 4) studies are available. Therefore, alternative metaanalytic methods are needed. In the case of binary data, the “commonrho” betabinomial model has shown good results in situations with sparse data or few studies. The major concern of this model is that it ignores the fact that each treatment arm is paired with a respective control arm from the same study. Thus, the randomisation to a study arm of a specific study is disrespected, which may lead to compromised estimates of the treatment effect. Therefore, we extended this model to a version that respects randomisation.
The aim of this simulation study was to compare the “commonrho” betabinomial model and several other betabinomial models with standard metaanalyses models, including generalised linear mixed models and several inverse variance random effects models.
Methods
We conducted a simulation study comparing betabinomial models and various standard metaanalysis methods. The design of the simulation aimed to consider metaanalytic situations occurring in practice.
Results
No method performed well in scenarios with only 2 studies in the random effects scenario. In this situation, a fixed effect model or a qualitative summary of the study results may be preferable. In scenarios with 3 or 4 studies, most methods satisfied the nominal coverage probability. The “commonrho” betabinomial model showed the highest power under the alternative hypothesis. The betabinomial model respecting randomisation did not improve performance.
Conclusion
The “commonrho” betabinomial appears to be a good option for metaanalyses of very few studies. As residual concerns about the consequences of disrespecting randomisation may still exist, we recommend a sensitivity analysis with a standard metaanalysis method that respects randomisation.
Introduction
Metaanalyses (MAs) are used to summarise the results of studies on a specific research question. If the number of studies is large and the sample sizes within these studies are not too small, standard inverse variance random effects models (IVREMs) can provide valid estimates. However, if only a few (≤ 10) studies are included in the MA, the IVREMs perform poorly [1,2,3]. The DerSimonianLaird (DSL) method leads to too narrow 95% confidence intervals (CIs) with poor coverage probabilities below 95%, especially in the case of few studies. The HartungKnappSidikJonkman (HKSJ) method generally holds the type I error, but frequently results in extremely wide 95% CIs in the case of very few (2 – 4) studies.
The worse performance in the case of few studies is a particular challenge, because such MAs are frequently performed in systematic reviews of interventions. For example, in an analysis of 14,886 MAs from the Cochrane Library, the median number of studies in MAs was 3 and the 3rd quartile was 6 [4].
In the case of binary data, alternatives to IVREM methods have been proposed. Because outcome (success, failure) and study arm (treatment, control) for each patient can be reconstructed from studies’ fourfold tables, the generalised linear mixed model (GLMM) framework (and generally speaking all logistic regression models accounting for dependent data) can be applied to MAs [5]. The “commonrho” betabinomial model (BBM) showed good results when pooling data of randomised controlled trials (RCTs), especially in the case of very few studies and/or rare events in the MA [6, 7]. However, there are some concerns about the model because it ignores the fact that each treatment arm is paired with a respective control arm, both originating from the same study (disrespecting the randomisation to a study arm of a specific study).
Therefore, we conducted a simulation study to compare existing BBMs and extensions with established models, such as GLMMs, HKSJ and DSL, especially in situations with very few (2 – 4) studies and for a wide range of risks, including rare events. Our focus was on BBM extensions that accounted for the pairing of a treatment arm with a control arm of the same study by implementing a random effect for the study or conditioning on the study in the maximum likelihood estimation.
The outline of the paper is as follows. In the 2nd chapter, we describe the statistical models for MAs that were included in the comparison (Models section) and explain how the simulation study was conducted (Simulation study section). In the 3rd chapter, we present the results of the simulation study. In the 4th chapter, we discuss the results, and in the 5th chapter, we conclude with final remarks and recommendations for practice.
Methods
We consider situations where K studies compare a binary outcome between two study arms (i = 1 [or T for treatment] and i = 0 [or C for control]). For each study k (k = 1, …, K), n_{kT} and n_{kC} denote the sample size for the treatment and control arm, y_{kT} and y_{kC} the number of events in the treatment and control arm, and θ_{k} the treatment effect with a specific withinstudy variance \({\sigma}_k^2\). We are interested in the estimation of the overall treatment effect θ and use the effect measures odds ratio (OR) and relative risk (RR) to quantify the effect between the treatment and control arm.
Models
We compared the following metaanalytic models and methods in our simulation study:

Betabinomial models (BBMs)

◦ Standard (“commonrho”) betabinomial model (BBST)

◦ Standard betabinomial model with an additional random treatment effect (BBFR)

◦ Two “commonbeta” betabinomial models (BBCB1 and BBCB2)


Generalised linear mixed models (GLMMs)

◦ Generalised linear mixed model with a fixed intercept and random treatment effect (GLFR)

◦ Generalised linear mixed model with a random intercept and random treatment effect (GLRRI)


Inverse variance random effects model (IVREM)

◦ DerSimonianLaird (DSL) method

◦ HartungKnappSidikJonkman (HKSJ) method


MantelHaenszel (MH) method

Peto odds ratio (POR) method

Collapsed table (COLL)
and describe them in the following sections.
Betabinomial models
Standard betabinomial model
The standard betabinomial model (BBST) [6, 7] assumes that the observed number of events in the control arm y_{kC} (k = 1, …, K) follows a binomial distribution Bin(π_{kC}, n_{kC}), where the event probability π_{kC} is not fixed, but beta distributed with parameters α_{C} and β_{C}. These parameters are assumed to be constant over all control arms of the studies. The individual binary event z_{kCj} (j = 1, …, n_{kC}; \({y}_{kC}=\sum_j^{n_{kC}}{z}_{kC j}\)) is sampled with a different π_{kC}. The expected value and variance of π_{kC} are:
and
with
and
and y_{kC} is betabinomially distributed with the expected value
and variance
Because the probabilities for two individual binary events in the control arm are sampled from the same beta distribution, these events are correlated. The intraclass correlation \({\rho}_C= corr\left({z}_{{kC j}_1},{z}_{kC{j}_2}\right)\) between two individual binary events in the control arm k (k = 1, …, K; j_{1}, j_{2} = 1, …, n_{kC}; j_{1} ≠ j_{2}) is
and is assumed to be equal over all y_{kC} (k = 1, …, K). Further, it is assumed that individual binary events from different control (and treatment) arms are uncorrelated, \(corr\left({z}_{k_1{Cj}_1},{z}_{k_2C{j}_2}\right)=0\) for k_{1} ≠ k_{2}.
The log likelihood of the betabinomial distribution of all control arms can be written in closed form as
with
where Γ denotes the gamma function,
and
The same formulas hold true for the number of events in the treatment arm y_{kT} (k = 1, …, K) with n_{kT}, π_{T}, α_{T}, β_{T}, μ_{T}, ν_{T} and ρ_{T}. The log likelihood for the treatment arms ℓ_{T}(α_{T}, β_{T}) is given accordingly.
Importantly, in the BBST it is assumed that ρ_{C} = ρ_{T} = ρ, which is equivalent to α_{C} + β_{C} = α_{T} + β_{T}. In other words, all individual binary events within a study arm are correlated with the same ρ, regardless of therapy.
The treatment effect θ = b_{T} = g^{−1}(μ_{T})/g^{−1}(μ_{C}) is modelled via the link function
where b_{0} denotes the risk of an event in the control arm and i the study arm (1 = treatment; 0 = control). In our simulation study, the link functions are the logit and the natural log to measure the treatment effect as log OR and log RR.
Because g(μ_{C}) = b_{0}, g(μ_{T}) = b_{0} + b_{T} and α_{C} + β_{C} = α_{T} + β_{T}, one can write
Therefore, only three parameters (α_{T}, α_{C}, β_{C}) have to be estimated in this model.
One advantage of the BBST is that no continuity correction has to be used if there are studies with no events in one study arm (singlezero studies). Furthermore, studies without any events in both study arms (doublezero studies) are not ignored in the analysis and contribute to the overall effect estimation [6]. The only situations where the BBST cannot estimate the OR and RR (but can estimate the risk difference) are situations where no events occur in one study arm (e.g. the treatment arm) over all studies.
In the BBST, the event probability in the control arm π_{kC} is random but the treatment effect is considered to be fixed across all studies. Thus, although the BBST is a true random effects model, from a metaanalytic point of view, it is a model with a fixed treatment effect.
Furthermore, the BBST estimates the treatment effect via μ_{T} and μ_{C}. Therefore, the fact that the treatment and control arm originate from the same study is ignored in the process of parameter estimation. Thus, the BBST disrespects the randomisation to a study arm of a specific study. According to Senn [8] and Piepho et al. [9] it is unlikely that this property is detrimental in situations where the same treatments are evaluated across trials, because the effects are comparable between the studies. This was indicated by recent simulation results, where the BBST performed well [6, 7].
Standard betabinomial model with additional random effect
To deal with the aforementioned properties of BBST as a fixed effect model that disrespects randomisation, we implemented another BBM (BBFR) where the treatment effect θ = b_{T} = g^{−1}(μ_{T})/g^{−1}(μ_{C}) is modelled as
with \({\gamma}_k\sim N\left(0,{\sigma}_{BBFR}^2\right)\). By adding a random effect to the treatment effect, this model respects the randomisation to a study arm of a specific study.
Like the BBST, the BBFR takes the information of all studies into account and therefore needs no continuity correction when single or doublezero studies are included in the MA.
When constructing the 95% CI for the treatment effect b_{T}, Mathes and Kuss [7] showed that using the tdistribution rather than the normal distribution led to better performance of the BBST. Therefore, for both models (BBST, BBFR) we calculated 95% CIs for the treatment effect b_{T} using the tdistribution
where \({\hat{b}}_{T(BB)}\) is the estimated treatment effect of BBST or BBFR, \({\hat{\sigma}}_{BB}\) the corresponding estimated standard error and t_{df; 0.975} the 0.975 quantile of the tdistribution with df degrees of freedom. We chose two different numbers of degrees of freedom: df = K − 1, which are the degrees of freedom for the HKSJ method, and df = 2K − 2, which is a reasonable compromise between the 2K and 2K − 3 degrees of freedom that were used by Mathes and Kuss [7] and that led to too narrow and too wide intervals, respectively.
Commonbeta BBM
In the BBST it is assumed that the intraclass correlation ρ is equal for all treatment and control arms implying that α_{C} + β_{C} = α_{T} + β_{T} holds true. Therefore, this model is sometimes called the “commonrho” model. Another possibility is to assume that β is equal in all groups (β_{C} = β_{T} = β). We call this model the “commonbeta” BBM. Similar to the BBST, the likelihood function of the “commonbeta” BBM can be written in closed form. Guimaraes [10], and more recently in the metaanalytic context Mathes and Kuss [11], were able to show that a negative binomial regression model for count panel data can be interpreted as this “commonbeta” BBM.
We considered two versions of the “commonbeta” BBM in our simulation. BBCB1, which disrespects the randomisation to a study arm of a specific study by conditioning on study group, while BBCB2 respects the randomisation to a study arm of a specific study by conditioning on the study.
As for BBST and BBFR, we constructed 95% CIs using the tdistribution for the treatment effect b_{T}
where \({\hat{b}}_{T(BBCB)}\) is the estimated treatment effect of BBCB1 or BBCB2, \({\hat{\sigma}}_{BBCB}\) the corresponding estimated standard error, and t_{df; 0.975} the 0.975 quantile of the tdistribution with df = K − 1 or 2K − 2 degrees of freedom.
Generalised linear mixed models
Generalised linear mixed models (GLMMs) [5] are probably the most common alternatives to standard MA methods (Inverse variance random effects models section) with binary data because of their flexibility. A GLMM with random treatment effect θ_{k} = θ + ϵ_{k} can be expressed as
where g(∙) is the link function for the OR (logit) or RR (log), π_{ki} the probability of an event in study arm i (i = 1: treatment; i = 0: control) of study k (k = 1, …, K), γ_{k} the intercept (baseline risk of an event) of study k and ϵ_{k}~N(0, τ^{2}).
We included two GLMMs in our simulation. The first model has a fixed intercept and a random treatment effect (GLFR) (similar to model 2 in Jackson et al. [5] originally suggested by Simmonds and Higgins [12]). The second GLMM included has a random intercept (\({\gamma}_k\sim N\left(\gamma, {\sigma}_{GLRRI}^2\right)\)) and an uncorrelated random treatment effect (GLRRI), similar to model 3 in Jackson et al. [5].
Like the BBST, the GLMM is a random effects model. But as the treatment effect is random it is more comparable to metaanalytic REMs than the BBST.
Like BBMs, GLMMs take the information of all studies into account and therefore do not need a continuity correction for single or doublezero studies.
We calculated 95% CIs for log OR and log RR using normal approximation, therefore
where \({\hat{\theta}}_{GLMM}\) is the estimated overall effect (log OR or log RR) in the analysed model (GLFR or GLRRI) and \({\hat{\sigma}}_{GLMM}\) the corresponding standard error.
Inverse variance random effects models
The metaanalytic random effects model (REM) assumes no common effect for all studies but instead assumes that the mean of all study effects is the mean of the distribution of the true effect [13]. The study effects are usually assumed to follow a normal distribution. The treatment effect in study k can be expressed as \({\hat{\theta}}_k={\theta}_k+{\varepsilon}_k\) with θ_{k} = θ_{REM} + δ_{k}, δ_{k}~N(0, τ^{2}) and \({\varepsilon}_k\sim N\left(0,{\sigma}_k^2\right)\).
The overall effect θ_{REM} of this REM can be estimated by the inverse variance approach
where \({w}_{k(REM)}=1/\left({\sigma}_k^2+{\tau}^2\right)\) are the studyspecific weights, \({\sigma}_k^2\) is the withinstudy variance, and τ^{2} is the betweenstudy variance (heterogeneity).
In the case of binary data and for OR and RR as the effect sizes, \({\hat{\theta}}_k\) and \({\hat{\theta}}_{REM}\) are generally analysed on the logarithmic scale, thus representing the log OR and the log RR. A continuity correction is applied to singlezero studies by adding 0.5 to the number of events in both groups. Doublezero studies are ignored for parameter estimation.
DerSimonianLaird method
The DerSimonianLaird (DSL) method [14] was regarded as the gold standard for performing MAs until it was criticized in recent years for being anticonservative (i.e., producing too narrow CIs) [15].
The DSL estimator \({\hat{\theta}}_{DSL}\) uses weights \({w}_{k(DSL)}=1/\left({\sigma}_k^2+{\tau}_{DSL}^2\right)\) where
is estimated using the methodofmoments principle [14, 16]. Q is the heterogeneity statistic \(Q=\sum_{k=1}^K{w}_{k(FEM)}\times {\left({\hat{\theta}}_k{\hat{\theta}}_{FEM}\right)}^2\) and \({\hat{\theta}}_{FEM}\) is the pooled effect of a fixed effect model with weights \({w}_{k(FEM)}=1/{\sigma}_k^2\).
The 95% CI for θ_{REM} is given by
The standard error is given by
We included the DSL method in our simulation because it is still in use and is important, at least for sensitivity analyses [17].
HartungKnappSidikJonkman method using the PauleMandel heterogeneity variance estimator
For the method of HartungKnappSidikJonkman (HKSJ) [18, 19] different weights w_{k(HKSJ)} can be applied to calculate the overall estimator \({\hat{\theta}}_{HKSJ}\), depending on what estimator for the betweenstudy variance (heterogeneity) is used. For the analysis presented here, we used \({w}_{k(HKSJ)}={w}_{k(PM)}=1/\left({\sigma}_k^2+{\tau}_{PM}^2\right)\), where \({\tau}_{PM}^2\) is the PauleMandel (PM) estimator of τ^{2} [20,21,22]. The PM estimator of τ^{2} is derived from the generalised Q statistic
by setting \(Q\left({\tau}_{PM}^2\right)\) to its expected value K − 1 with \({w}_{k(PM)}=1/\left({\sigma}_k^2+{\tau}_{PM}^2\right)\) and \(\hat{\theta}\left({\tau}_{PM}^2\right)=\left(\sum_{k=1}^K{w}_{k(PM)}\times {\hat{\theta}}_k\right)/\left(\sum_{k=1}^K{w}_{k(PM)}\right)\). The equation is solved by iterating \({\tau}_{PM}^2\) until convergence is reached. If no solution exists, \({\hat{\tau}}_{PM}^2\) is set to 0.
The 95% CI for θ_{REM} is given by
where t_{K − 1; 0.975} is the 0.975 quantile of the tdistribution with K − 1 degrees of freedom,
and
In general, this CI tends to be wider than the 95% CI of the DSL method. However, in very homogeneous data situations, this is not always the case. Therefore, Knapp and Hartung [23] suggested an ad hoc modification of q, q^{∗} = max {1, q}. If the ad hoc modification is used, the 95% CI of the HKSJ method is always wider than the 95% CI of the DSL method. In the simulation study, we followed the recommendations of the literature to carry out sensitivity analyses using a fixed effect model or the DSL method to decide on whether the modification is needed [17, 24]. If the 95% CI of HKSJ was narrower than the 95% CI of DSL, the ad hoc modification was used.
We included the HKSJ method in our simulation because it is wellestablished that it performs better than the DSL method [17] and is recommended as the standard approach in the Cochrane Handbook in combination with the PM estimator for τ^{2} [25]. Furthermore, it is the IQWiG standard method for MAs with five or more studies [26]. It is wellknown that for MAs with less than five studies, 95% CIs of the HKSJ method tend to be too wide but in general, the method does not lead to increased type I errors.
Other models
MantelHaenszel method
The MantelHaenszel (MH) method [27] assumes a fixed effect model which is based on the assumption that all studies in the MA have a common effect θ_{FEM}. The resulting estimate is a weighted average of the studyspecific risk ratios or risk differences.
The MH odds ratio of the overall effect is given by
where \({\hat{OR}}_{k(MH)}=\frac{y_{kT}\times \left({n}_{kC}{y}_{kC}\right)}{\left({n}_{kT}{y}_{kT}\right)\times {y}_{kC}}\) is the estimated odds ratio, \({w}_{k(MH)}=\frac{\left({n}_{kT}{y}_{kT}\right)\times {y}_{kC}}{n_k}\) the weight, and n_{k} = n_{kT} + n_{kC} the sample size of study k (k = 1, …, K).
No continuity correction was applied. Therefore, single and doublezero studies were ignored during the analysis. We estimated 95% CIs using normal approximation.
We included the MH method in our analysis because it is the standard fixed effect model for binary data in Cochrane Reviews and performs better than the standard (inverse variance) fixed effect model in the case of sparse data [25].
Peto odds ratio method
The Peto odds ratio (POR) method [28] was introduced as the effect estimator for the real underlying OR in the data situation of rare events. Later it was shown that POR is an independent effect measure and cannot be used as approximation of the true OR in all (rare event) data situations [29].
The pooled log POR is estimated as
where O_{k} is the observed number of events in the treatment arm, E_{k} the expected number of events in the treatment arm under the assumption of no treatment effect, and V_{k} the variance of their difference.
No continuity correction was applied. Singlezero studies are included in the analysis by definition, whereas doublezero studies are ignored. We estimated 95% CIs using the normal approximation.
Although this method has major limitations [29, 30], we considered the POR in our analysis because according to the Cochrane Handbook, it is the standard MA method for small intervention effects or very rare events [25].
Collapsed table
This method (COLL) ignores the fact that data were collected from different studies and aggregates them in one single fourfold table. The effect is estimated using standard methods for 2 × 2 tables [31]. The procedure assumes a constant underlying risk of an event across all studies, which is rarely given, and therefore the method is vulnerable to Simpson’s paradox [32, 33].
Because of its simplicity and because the differences in event rates across studies might be negligible in scenarios with few events and studies, we included the method in our analysis.
We applied a continuity correction if necessary and estimated 95% CIs using normal approximation.
Simulation study
We performed a simulation study using SAS/STAT® software Version 9.4 (SAS Institute Inc., Cary, NC, USA) for Microsoft Windows 10 to compare the statistical properties of the different metaanalytic methods. In an attempt to investigate the methods in realistic scenarios, the values of the design factors in the simulation study were chosen from MAs actually performed. For this purpose, we used the review by Turner et al. [4] The review included 1991 systematic reviews from the Cochrane Database of Systematic Reviews and analysed 14,886 MAs of dichotomous outcomes from 77,237 single studies.
Design of simulations
In our simulation study, we varied the following parameters: the number of studies in the MA, the sample size of a single study, the event probability in the control arm, the heterogeneity between studies in the MA, and the effect size in the MA.
We chose 2, 3, 4, 5 and 10 as the number of studies in the MA. Ten studies were chosen in addition to gain an impression of how the models perform in rather uncritical scenarios. For each number of studies in the MA, we simulated 10,000 data sets each under the null (H_{0}) and under the alternative hypothesis (H_{1}). In each data set, we sampled the values of the other parameters from distributional functions based on the data from Turner et al. [4] For example, the sample size n_{k} for a single study in the MA was sampled using a lognormal distribution with μ_{ND} = 4.615 and σ_{ND} = 1.1, resulting in a mean overall sample size of 185.7 (median: 102.0; 1st quartile: 50.0; 3rd quartile: 213.0). Table 1 shows the distributional assumptions and the resulting data values. The data were simulated according to the REM from the Inverse variance random effects models section. The simulation process was as follows:

1.
For each MA, we sampled the true risk π_{C, true} in the control arm and the heterogeneity τ^{2} between the studies in the MA from the distributional functions given in Table 1. Under H_{1}, we did the same for the effect size θ = log OR and log RR. Under H_{0}, the effect size was set to zero (i.e. OR = 1 and RR = 1).

2.
For the kth study in the MA, we

a.
sampled the study size n_{k} and the size of the treatment (n_{kT}) and control arm (n_{kC}) using the distributional functions given in Table 1

b.
sampled the number of events in the control arm y_{kC} using a binomial distribution with π_{C, true} as event probability and n_{kC} as number of experiments

c.
sampled an individual heterogeneity variance \({\tau}_k^2\) using the sampled true heterogeneity and assuming that it follows a normal distribution within the kth MA

d.
calculated the true risk in the treatment arm π_{kT, true} using π_{C, true}, θ, \({\tau}_k^2\) and the following formula: \({\pi}_{kT, true}=\exp \left(\textrm{logit}\left({\pi}_{C, true}\right)+\theta +{\tau}_k^2\right)/\left(1+\exp \left(\textrm{logit}\left({\pi}_{C, true}\right)+\theta +{\tau}_k^2\right)\right)\)

e.
sampled the number of events in the treatment arm y_{kT} using a binomial distribution with π_{kT, true} as event probability and n_{kT} as number of experiments.

a.
Overall, we simulated 5 (number of studies in the MA 2, 3, 4, 5, 10) × 2 (under H_{0} and under H_{1}) × 10,000 data sets = 100,000 MAs each for the OR and for the RR.
We performed a sensitivity analysis to assess the robustness of the results regarding heterogeneity. For each MA, we calculated Cochran’s Q test [13] in order to gain an impression of whether the results depend on homogeneity of the data situations. Although we are aware that this dichotomisation is somewhat arbitrary, we used the Cochran’s Q test for the purpose of sensitivity analyses because in practical applications, MAs will frequently not be performed, at least when the test for heterogeneity is statistically significant.
Parameter estimation in the models
For parameter estimation, we used the SAS/STAT® software procedure NLMIXED for BBST and BBFR, COUNTREG for BBCB1 and BBCB2, GLIMMIX for GLFR and GLRRI, and FREQ for MH and COLL. For HKSJ, DSL and POR, we programmed our own syntax that was validated using R 3.3 [36] and the metafor package [37].
Because we used the COUNTREG procedure for parameter estimation, we were only able to estimate the OR but not the RR. In the GLIMMIX procedure, we used maximum likelihood estimation based on adaptive quadrature (METHOD = QUAD) with 1 quadrature point (QPOINTS = 1), which is equivalent to the Laplace approximation. We decided to use the Laplace approximation because we assumed that this would be most robust [38].
Performance measures
To assess the performance of the methods we calculated the following measures:

Number of converged simulation runs with estimated effect and standard error (R): Sometimes the procedures converged and an effect was estimated but no standard error was given (most notably when using the NLMIXED procedure). Because such results would cause problems of interpretation, we counted these runs as nonconverged. All other measures were based on R, the number of converged simulation runs with an estimated effect and standard error.

(Absolute) bias \({\hat{\theta}}_r{\theta}_r\): Difference between the estimated (\({\hat{\theta}}_r\)) and true effect (θ_{r}); r = 1, …, R.

Percentage bias under H_{1}\(\left(100\times \left({\hat{\theta}}_r{\theta}_r\right)\right)/{\theta}_r\): Ratio of the bias (\({\hat{\theta}}_r{\theta}_r\)) and true effect (θ_{r}); r = 1, …, R.

Coverage probability: Proportion of converged simulation runs where the 95% CI included the true effect θ_{r}; r = 1, …, R.

Length of 95% CI: Difference (CI_{U, r} − CI_{L, r}) of upper (CI_{U, r}) and lower (CI_{L, r}) confidence limit of the 95% CI for θ_{r}; r = 1, …, R.

Power under H_{1}: Proportion of converged simulation runs under H_{1} where the 95% CI excluded the null effect.
Bias, percentage bias and length of 95% CI were calculated on the corresponding log scale, i.e. log OR or log RR. For these measures, the median as well as the 1st and 3rd quartile are presented.
The simulation code containing the data generation, parameter estimation, and the calculation of the performances measures is available in the Supporting Information.
Results
In the following sections, we describe the results for the OR of all methods under the null hypothesis (Results for the odds ratio under the null hypothesis section), and alternative hypothesis (Results for the odds ratio under the alternative hypothesis section). In the Direct comparison of results for betabinomial models section, we compare the results of the BBMs, especially BBST and BBFR. The results of the RR are discussed in the Results for the relative risk section. The results of the sensitivity analysis are presented in the Sensitivity analysis section and the main results are summarized in the Summary of results section.
Although the results of all methods are described in this chapter, in tables and figures we focused on the BBST, BBCB1, BBCB2, GLFR, HKSJ and DSL. As BBST and BBFR yielded almost identical results (Direct comparison of results for betabinomial models section and Discussion), we refrained from showing them both in tables and figures. The results of all methods can be found in the Supporting Information.
For the BBMs, the results for coverage probability, length of 95% CI and power, with CIs using 2K − 2 degrees of freedom, are presented and discussed. Results for the CIs with K − 1 degrees of freedom, which were generally worse, can be found in the Supporting Information.
Results for the odds ratio under the null hypothesis
Number of converged simulation runs
The methods HKSJ, DSL, COLL and POR were only marginally affected by convergence problems (< 0.5% for all scenarios) due to their construction. The same held true for GLFR (R ≥ 9988 for all scenarios). The BBMs (BBST, BBFR, BBCB1, BBCB2) converged in more than 95% of the simulation runs. The number of converged runs was lower for MH in scenarios with 2 studies (R = 9394) but increased up to 10,000 in scenarios with 10 studies. The GLRRI had the lowest number of converged runs, with 8469 (2 studies) to 6896 runs (10 studies) (see Fig. 1 and Table S1 in the Supporting Information).
Bias
For all methods, the median bias was similar, mainly positive, increased with an increasing number of studies, and was low (< 0.04 on log OR scale). Because bias was calculated on the log OR scale (bias \(=\log \left(\hat{OR}\right)\log (OR)=\log \left(\hat{OR}/ OR\right)\)), this could be interpreted as relative effect of the ORs. Therefore, the estimated OR increased up to 4% in median if compared to the true OR (see Table 2 and Table S2 in the Supporting Information).
Coverage probability
Coverage probabilities were at or above 95% for the BBMs (BBST, BBFR, BBCB1, BBCB2) and HKSJ for all scenarios. Coverage probability for GLFR was at or above 95% for 2, 3 and 4 studies and fell below 95% for 5 (93.9%) and 10 studies (90.4%). All other methods had coverage probabilities below 95% for all scenarios (GLRRI: ≤ 85.7%; DSL: ≤ 93.3%; MH: ≤ 81.6%; POR: ≤ 80.6%; COLL: ≤ 82.2%) (see Fig. 2 and Table S3 in the Supporting Information).
BBMs (BBST, BBFR, BBCB1, BBCB2), HKSJ and GLFR had coverage probabilities ≥97.3% in scenarios with 2 studies. In scenarios with 3, 4, 5, and 10 studies, BBMs (BBST, BBFR, BBCB1, BBCB2) were closer to 95% than HKSJ (≥ 96.2% in all scenarios). The same applied to the GLFR in scenarios with 3 and 4 studies (see Fig. 2 and Table S3 in the Supporting Information).
Coverage probabilities of the BBMs were closer to 95%, i.e., less conservative, if confidence intervals with 2K − 2 degrees of freedom were used compared to the use of K – 1 degrees of freedom (see Table S3 in the Supporting Information).
Length of 95% CI
The length of the 95% CI for log OR was largest and approximately the same for HKSJ and GLFR in the scenario with 2 studies. The length was far shorter for BBMs. With an increasing number of studies in the MA, the length of the 95% CIs converged between the methods, but HKSJ always had the widest CIs (see Table 3 and Table S4 in the Supporting Information).
Results for the odds ratio under the alternative hypothesis
The results under the alternative hypothesis were quite similar to the results under the null hypothesis for all performance measures investigated. Therefore, for the number of converged simulation runs, bias, squared error, coverage probability and length of 95% CI, only important differences between null and alternative hypothesis are mentioned.
Number of converged simulation runs
The number of converged runs for GLRRI increased by about 100 runs up to between 8590 (2 studies) and 7083 runs (10 studies). GLRRI still had the lowest convergence rate (see Table S5 in the Supporting Information).
Bias
The median absolute bias of log OR was higher than under the null hypothesis and mainly positive, but still small (< 0.05), for most methods. The most notable exception was BBCB2 with a median bias > 0.10 for all scenarios (see Table 4 and Table S6 in the Supporting Information).
Percentage bias
The percentage bias of log OR was defined as \(\left(100\times \left(\log \left(\hat{OR}\right)\log (OR)\right)\right)/\log (OR)\) and because log OR < 0, a negative percentage bias means an overestimation of the log OR. The median percentage bias was low (about − 6 and 0%) for the BBMs (BBST, BBFR, BBCB1), except for BBCB2, and the GLMMs (GLFR, GLRRI). The median values for HKSJ, DSL, MH, POR and COLL were higher, with values between − 14 and − 7%. BBCB2 had much worse median values, about − 30% (see Table 5 and Table S7 in the Supporting Information).
Coverage probability
Only HKSJ had coverage probabilities at or above 95% for all scenarios. BBCB2 had coverage probabilities < 95% for scenarios of 3 or more studies. The other BBMs (BBST, BBFR, BBCB1) had coverage probabilities at, above, or marginally below 95% for all scenarios. Coverage probabilities for GLFR were at or above 95% only for scenarios with 2, 3 or 4 studies. All other methods had coverage probabilities below 95% for all scenarios (GLRRI: ≤ 86.6%; DSL: ≤ 92.5%; MH: ≤ 82.2%; POR: ≤ 82.1%; COLL: ≤ 83.3%) (see Fig. 3 and Table S8 in the Supporting Information).
Length of 95% CI
Results for the length of the 95% CI were similar to results under the null hypothesis, with HKSJ and GLFR having the broadest intervals in all scenarios (see Table S9 in the Supporting Information).
Power
Power for methods with coverage probability of ≥95% under the null hypothesis (BBST, BBFR, BBCB1, BBCB2, and HKSJ in all scenarios and GLFR in scenarios with 2, 3 and 4 studies) was quite low (still < 40% in scenarios with 10 studies). Power for BBMs (BBST, BBFR, BBCB1, BBCB2) was higher than for HKSJ in all scenarios. In scenarios with 2 studies, none of these methods showed a power > 5%, that is, no method yielded satisfactory results. In scenarios with 3 and 4 studies, power was still low (maximum 21.0% for BBCB1) but the differences in the methods became visible. Power was highest for BBST, BBFR and BBCB1, followed by BBCB2 and HKSJ with the lowest power (see Fig. 4 and Table S10 in the Supporting Information).
Methods with coverage probabilities < 95% under the null hypothesis (GLRRI, DSL, MH, POR, COLL) had higher power up to 55% (see Table S10 in the Supporting Information).
The small power values were to be expected due to the fact that the true ORs were near 1 (> 0.83) for about 25% of the simulations, the moderate sample sizes (mean: 185.7; median: 102.0), and only few studies in the MAs.
Direct comparison of results for betabinomial models
BBST and BBFR showed almost identical results. We assumed that one reason could be the maximum likelihood approximation method. Therefore, we tried to use another approximation method, the GaussHermite quadrature with two quadratic points (QPOINTS = 2). However, in the case of a higher number of quadrature points, a floating point exception error occurred at one point, stopping the whole simulation. Therefore, we could not run the complete simulation using more quadrature points. We tried to reanalyse the simulated data using the GaussHermite quadrature with 2 quadrature points for the BBFR. From the few existing results, it appears that BBST and BBFR vary the most if the ORs in the studies are very heterogeneous. Of notice, these were mainly data situations with strong heterogeneity (P value of Q test < 0.05).
The coverage probability of the BBCB2 was below the nominal level under the alternative hypothesis in scenarios with 3 or more studies. In contrast, BBCB1 performed well and very similarly to BBST considering all performance measures.
Results for the relative risk
For the RR, we only considered BBST, GLFR, HKSJ, DSL and COLL in our simulation study. Because the results of BBST and BBFR were almost identical, we focussed on BBST. As mentioned earlier, it was not possible to compute RRs for BBCB1 and BBCB2. GLRRI and MH performed quite poorly considering the OR. Thus, we saw no need to reconsider them again. DSL and COLL also performed quite poorly, but we considered them due to their simplicity and as “benchmark” methods.
Results for RR were comparable to OR, except for the fact that BBST and GLFR struggled with convergence problems. BBST converged only in about 82% of each scenario (about 98% for OR) and GLFR between 90 and 95% (> 99% for OR). Noticeably, coverage probability for GLFR was higher and nearer the 95% level than that for OR (see Tables S11 – S20 in the Supporting Information).
Sensitivity analysis
The sensitivity analysis considering only data situations with no statistically significant heterogeneity (P value of Q test > 0.05) did not fundamentally alter the simulation results of the performance measures. The biggest influence was seen in coverage probabilities. Coverage probabilities below 95% increased and approached 95%. Coverage probabilities above 95% remained more or less stable. The biggest improvement was seen in MH, POR and COLL, but their coverage probabilities were still far below 95% (e.g. OR under H_{0} in scenarios with 4 studies: MH: 88.9%, POR: 88.4%, COLL: 88.9%). The most relevant improvement was shown in the GLFR, where coverage probabilities were far more near 95% in scenarios of 5 and 10 studies than in data situations where the appropriateness of pooling was not considered (see Fig. 5 as an example).
Summary of results
The betabinomial models BBST and BBCB1 (BBCB1 only for OR) performed well in scenarios with 3, 4, 5 and 10 studies. HKSJ was the only standard MA method that had adequate performance measures in these scenarios, although the coverage probability was very high (> 98%) for scenarios < 5 studies. In scenarios with only 2 studies, no method showed coverage probabilities near 95%; especially HKSJ and GLFR were far too conservative. Power was very low in this scenario; therefore, for 2 studies no method appeared appropriate. In scenarios with 3 and 4 studies, BBST, BBCB1 (BBCB1 only for OR), GLFR and HKSJ performed best, with the first two methods having higher power than the last two. In scenarios with 5 and 10 studies, BBST, BBCB1 (BBCB1 only for OR) and HKSJ performed best, with the first two methods having higher power than the last one.
Discussion
We conducted a simulation study to compare BBMs with various standard metaanalytic methods in the case of very few studies. The BBST and the BBCB1 (BBCB1 only for OR) showed good results in the given data situations that were based on realistic data situations of Cochrane Reviews. The only standard MA method that showed acceptable results was HKSJ [39].
The attempt to extend the BBST by a random treatment effect term for the study, due to concerns about disrespecting the randomisation, failed. From the few results we obtained from the BBFR GaussHermite quadrature, it could be seen that the 2 models varied the most in situations where a lot of heterogeneity between the study effects was present. As BBFR uses a random effect attached to the treatment effect in addition to the random intercept, this behaviour is to be expected, because only then “enough” additional heterogeneity remains to be estimated, i.e. not all heterogeneity goes into the random intercept. In practice, pooling will often not be appropriate in situations where even in the case of few studies there is such strong heterogeneity in the data. Thus, there is probably little benefit in using other approximation methods than Laplace in the case of sparse data, because either it has no impact on results, or it has only an impact in situations where there is a lot of heterogeneity. Our finding is in agreement with a recent study that assessed different approximations methods to perform metaanalyses using GLMMs in the case of rare events. This study found that the GaussHermite quadrature is not superior to the Laplace approximation [38]. Thus, there seems to be no benefit in using BBFR in practical applications for metaanalyses of few studies. However, further research is necessary to obtain findings that are more conclusive.
Disrespecting the randomisation to a study arm of a specific study had no strong influence on the results of BBST and BBCB1. This was already pointed out more than 10 years ago [8] and is in line with a recent publication on armbased (disrespecting randomisation) and contrastbased models (respecting randomisation) in network MA, where the authors conclude that both models are suitable for network MAs [40]. One reason for this is presumably that, especially in the case of few studies, heterogeneity cannot be estimated properly. As our simulation mirrors real MA situations, the problem may not be important in practice, as probably only a few data situations occur where this critical aspect of BBST and BBCB1 actually has negative consequences.
A very important point to keep in mind is the fact that we used adequate data situations based on RCT data. Our conclusions could be flawed when there is doubt about randomisation. Therefore, using another method respecting the randomisation to a study arm of a specific study as sensitivity analysis seems to be a reasonable approach. One could try a BBFR with more quadrature points than 1, but this might not work because of convergence problems. Thus, we recommend the use of standard procedures such as HKSJ in these situations. If the results and especially the point estimates are quite different, we would refrain from using BBST as a final method, because there probably is a problem with disrespecting the randomisation to a study arm of a specific study.
Surprisingly, BBCB2, which strictly respects randomisation, performed very badly regarding coverage probability under the alternative hypothesis. The narrow 95% CIs suggest that the standard errors are underestimated. This may be because only half of observations, concrete each study instead of each single group, contribute to the estimation of the parameters compared to BBCB1, where each arm of a study contributes to the estimation. Moreover, the assumption of homogeneity within one study may lead to overdispersion.
Although GLMMs are an intuitive alternative to the standard MA methods, GLFR and GLRRI performed quite poorly in many of our investigated scenarios. In Jackson et al. [5] both models (labelled as model 2 and 3) and their reparametrized versions were examined. In contrast to our results, the performance of GLFR was worse than that of GLRRI. Apart from this, the coverage probabilities for both models were below 95% for almost all 15 investigated scenarios, including 2 scenarios with 3 and 5 studies. The authors used maximum likelihood estimation based on adaptive quadrature with 7 (GLFR) and 1 quadrature point (GLRRI), whereas we used 1 quadrature point for each model. We do not think that the different number of quadrature points for GLFR can fully explain the differences in results. The noticeable difference in the simulation studies might be the difference in heterogeneity. Jackson et al. [5] simulated data with a median true τ for log OR of 0.024, whereas our median τ was \(\sqrt{0.273}=0.52\), that is, considerably larger. As seen in setting 15 in Jackson et al. [5], both models performed worse if a lot of heterogeneity was present (coverage probabilities below 90%). Therefore, in scenarios with higher heterogeneity, the GLFR with 1 quadrature point can be more robust than with 7 points, leading to better performance measures. The opposite is true for the GLRRI. The greater the heterogeneity, the lesser the assumption of one random effect for the intercept may be justified, resulting in worse performance measures. In agreement with this assumption, the GLRRI often did not converge in the case of large heterogeneity, i.e. large values of τ^{2}. As heterogeneity in the data is better detectable with an increasing number of studies in the MAs, this could also explain why, counterintuitively, the convergence of GLRRI decreases with an increasing number of studies. The partly differing results of these models in different simulation studies indicate how important the parameter settings in simulation studies could be, and thus stress the necessity of (independent) replications of results from simulation studies [41].
Some problems exist when using the RR in MAs [42, 43]. The problems arise from using the log link for the RR while the event probability π is bound to [0, 1] and are more pronounced if the event probabilities are large. Apart from problems with convergence, BBST and GLFR showed quite satisfactory results. A posthoc analysis revealed that the main reason for the worse convergence were large baseline probabilities, namely values of π_{kC} > 0.5. In addition, in the case of π_{kC} < 0.5, high heterogeneity between the studies and small RRs led to nonconverged simulations. Therefore, the log link affects the usability of these methods and further research is needed to improve the performance of these models.
An exact method for combining the effect sizes of the studies has recently been proposed [44]. This method was originally proposed for continuous data but could be easily implemented for binary data when using the logit or log link. Unfortunately, this method does not solve the problem of studies with no events in one or both study arms, because the effect size estimates of the studies are used to construct the 95% CI of the overall effect. Therefore, the same drawbacks exist with a continuous correction as with standard MA methods (HKSJ, DSL, MH for OR and RR).
Günhan et al. [45] investigated a Bayesian BBST. In a simulation study this model showed inappropriate coverage probabilities for very low OR (log OR < − 2), suggesting that BBMs are not suitable in situations where very large effects are expected. Because we tried to investigate realistic scenarios, such extreme values occurred in only a few simulations (≤ 10 in 10,000 simulations), which might be an explanation why no problems regarding coverage probability of the BBST were observed in our simulation study.
Our study has some limitations. As for all simulation studies, our results are only valid for data constellations we investigated. Because we based our simulation data on MAs actually performed in Cochrane Reviews, i.e. rather realistic scenarios, this is probably not a big problem. Another limitation is that we only investigated studies with balanced (1:1) randomisation schemes in the MA. In the case of a mixture of different randomisation schemes (1:1, 2:1, 3:1) the result of the BBST could be affected by the fact that the model disrespects randomisation. This highlights the importance of a sensitivity analysis using a method that respects randomisation.
Concerns could exist about the interpretability of the results with varying effect sizes under the alternative hypothesis. Therefore, we reanalysed the data by classifying the true effect sizes into four groups (< 1st quartile, between 1st quartile and median, between median and 3rd quartile, > 3rd quartile; data not shown). As expected, this had no impact on coverage probability. Likewise, for all other performance measures the values were bigger (smaller) for high effect size categories and smaller (bigger) for lower effect size categories. Thus, the given results can be interpreted as the mean/median of the different simulated effect sizes.
By chance, there were no doublezero studies in our MAs, but only singlezero studies. The number of MAs with at least one singlezero study varied from about 36% (2 studies) to 58% (10 studies). In an additional analysis (data not shown), the exclusion of MAs with singlezero studies led to similar results. However, this could be different if the number of single and doublezero studies increases, which would require further investigation.
Conclusion
In the case of very few (2 – 4) studies, the betabinomial models BBST and BBCB1 (BBCB1 only for odds ratio) are valuable alternatives to standard random effects metaanalytic models, if the corresponding 95% confidence intervals for BBST and BBCB1 are constructed using the t distribution with 2K − 2 degrees of freedom.
For metaanalyses with 2 studies, no general recommendation for a specific model can be given due to very conservative coverage probabilities and very low power of all investigated methods. The application of a fixed effect model, if appropriate, or a qualitative summary of the study results could be a solution. For metaanalyses with 3 and 4 studies, the BBST and BBCB1 can be recommended in conjunction with a sensitivity analysis using HKSJ or another adequate method for a random effects model. For metaanalyses with 5 or more studies, the use of HKSJ is recommended. BBST and BBCB1 are useful methods for sensitivity analyses in this case.
Availability of data and materials
All data generated in this study are available as supplementary material.
Abbreviations
 BBM:

Betabinomial model
 BBST:

Standard (“commonrho”) betabinomial model
 BBFR:

Standard betabinomial model with an additional random treatment effect
 BBCB:

“commonbeta” betabinomial model
 COLL:

Collapsed table
 DSL:

DerSimonianLaird
 GLMM:

Generalised linear mixed model
 GLFR:

Generalised linear mixed model with a fixed intercept and random treatment effect
 GLRR:

Generalised linear mixed model with a random intercept and random treatment effect
 HKSJ:

HartungKnappSidikJonkman
 IVREM:

Inverse variance random effects model
 MA:

Metaanalysis
 MH:

MantelHaenszel
 OR:

Odds ratio
 POR:

Peto odds ratio
 RCT:

Randomised controlled trial
 RR:

Relative risk
References
Bender R, Friede T, Koch A, Kuss O, Schlattmann P, Schwarzer G, et al. Methods for evidence synthesis in the case of very few studies. Res Synth Methods. 2018;9(3):382–92.
IntHout J, Ioannidis JP, Borm GF. The HartungKnappSidikJonkman method for random effects metaanalysis is straightforward and considerably outperforms the standard DerSimonianLaird method. BMC Med Res Methodol. 2014;14:25.
Veroniki AA, Jackson D, Bender R, Kuss O, Langan D, Higgins JPT, et al. Methods to calculate uncertainty in the estimated overall effect size from a randomeffects metaanalysis. Res Synth Methods. 2019;10(1):23–43.
Turner RM, Davey J, Clarke MJ, Thompson SG, Higgins JP. Predicting the extent of heterogeneity in metaanalysis, using empirical data from the Cochrane database of systematic reviews. Int J Epidemiol. 2012;41(3):818–27.
Jackson D, Law M, Stijnen T, Viechtbauer W, White IR. A comparison of seven randomeffects models for metaanalyses that estimate the summary odds ratio. Stat Med. 2018;37(7):1059–85.
Kuss O. Statistical methods for metaanalyses including information from studies without any eventsadd nothing to nothing and succeed nevertheless. Stat Med. 2015;34(7):1097–116.
Mathes T, Kuss O. A comparison of methods for metaanalysis of a small number of studies with binary outcomes. Res Synth Methods. 2018;9(3):366–81.
Senn S. Hans van Houwelingen and the art of summing up. Biom J. 2010;52(1):85–94.
Piepho HP, Williams ER, Madden LV. The use of twoway linear mixed models in multitreatment metaanalysis. Biometrics. 2012;68(4):1269–77.
Guimarães P. A simple approach to fit the betabinomial model. Stata J. 2005;5(3):385–94.
Mathes T, Kuss O. Betabinomial models for metaanalysis with binary outcomes: variations, extensions, and additional insights from econometrics. Res Methods Med Health Sci. 2021;2(2):82–9.
Simmonds MC, Higgins JPT. A general framework for the use of logistic regression models in metaanalysis. Stat Methods Med Res. 2016;25(6):2858–77.
Borenstein M, Hedges LV, Higgins JPT, Rothstein HR. Introduction to Metaanalysis. Chichester: Wiley; 2009.
DerSimonian R, Laird N. Metaanalysis in clinical trials. Control Clin Trials. 1986;7(3):177–88.
Higgins JP, Thompson SG, Spiegelhalter DJ. A reevaluation of randomeffects metaanalysis. J Royal Stat Soc Ser. 2009;172(1):137–59.
DerSimonian R, Kacker R. Randomeffects model for metaanalysis of clinical trials: an update. Contemp Clin Trials. 2007;28(2):105–14.
Jackson D, Law M, Rucker G, Schwarzer G. The HartungKnapp modification for randomeffects metaanalysis: a useful refinement but are there any residual concerns? Stat Med. 2017;36(25):3923–34.
Hartung J, Knapp G. A refined method for the metaanalysis of controlled clinical trials with binary outcome. Stat Med. 2001;20(24):3875–89.
Sidik K, Jonkman JN. A simple confidence interval for metaanalysis. Stat Med. 2002;21(21):3153–9.
Paule RC, Mandel J. Consensus values and weighting factors. J Res Natl Bur Stand. 1982;87(5):377–85.
Veroniki AA, Jackson D, Viechtbauer W, Bender R, Knapp G, Kuss O, et al. Recommendations for quantifying the uncertainty in the summary intervention effect and estimating the betweenstudy heterogeneity variance in randomeffectsmetaanalysis. Cochrane Database Syst Rev. 2015;1:1–72.
Langan D, Higgins JPT, Jackson D, Bowden J, Veroniki AA, Kontopantelis E, et al. A comparison of heterogeneity variance estimators in simulated randomeffects metaanalyses. Res Synth Methods. 2019;10(1):83–98.
Knapp G, Hartung J. Improved tests for a random effects metaregression with a single covariate. Stat Med. 2003;22(17):2693–710.
Wiksten A, Rucker G, Schwarzer G. HartungKnapp method is not always conservative compared with fixedeffect metaanalysis. Stat Med. 2016;35(15):2503–15.
Cochrane Handbook for Systematic Reviews of Interventions version 6.0 (updated July 2019) [www.training.cochrane.org/handbook].
Allgemeine Methoden; Version 6.1 [https://www.iqwig.de/methoden/allgemeinemethodenv61.pdf].
Mantel N, Haenszel W. Statistical aspects of the analysis of data from retrospective studies of disease. J Natl Cancer Inst. 1959;22(4):719–48.
Yusuf S, Peto R, Lewis J, Collins R, Sleight P. Beta blockade during and after myocardial infarction: an overview of the randomized trials. Prog Cardiovasc Dis. 1985;27(5):335–71.
Brockhaus AC, Bender R, Skipka G. The Peto odds ratio viewed as a new effect measure. Stat Med. 2014;33(28):4861–74.
Brockhaus AC, Grouven U, Bender R. Performance of the Peto odds ratio compared to the usual odds ratio estimator in the case of rare events. Biom J. 2016;58(6):1428–44.
Bradburn MJ, Deeks JJ, Berlin JA, Russell Localio A. Much ado about nothing: a comparison of the performance of metaanalytical methods with rare events. Stat Med. 2007;26(1):53–77.
Altman DG, Deeks JJ. Metaanalysis, Simpson's paradox, and the number needed to treat. BMC Med Res Methodol. 2002;2:3.
Lievre M, Cucherat M, Leizorovicz A. Pooling, metaanalysis, and the evaluation of drug safety. Curr Control Trials in Cardiovasc Med. 2002;3(1):6.
Fleishman AI. A method for simulating nonnormal distributions. Psychometrika. 1978;43(4):521–32.
Fan X, Felsövályi A, Sivo SA, Kennan SC. SAS for Monte Carlo studies: a guide for quantitative researchers. Cary: SAS Institute; 2002.
R: A language and environment for statistical computing [https://www.Rproject.org/].
Viechtbauer W. Conducting metaanalyses in R with the metafor package. J Stat Softw. 2010;36(3):1–48.
Ju K, Lin L, Chu H, Cheng LL, Xu C. Laplace approximation, penalized quasilikelihood, and adaptive gaussHermite quadrature for generalized linear mixed models: towards metaanalysis of binary outcome with sparse data. BMC Med Res Methodol. 2020;20(1):152.
van Aert RCM, Jackson D. A new justification of the HartungKnapp method for randomeffects metaanalysis based on weighted least squares regression. Res Synth Methods. 2019;10(4):515–27.
White IR, Turner RM, Karahalios A, Salanti G. A comparison of armbased and contrastbased models for network metaanalysis. Stat Med. 2019;38(27):5197–213.
Morris TP, White IR, Crowther MJ. Using simulation studies to evaluate statistical methods. Stat Med. 2019;38(11):2074–102.
Beisemann M, Doebler P, Holling H. Comparison of randomeffects metaanalysis models for the relative risk in the case of rare events: a simulation study. Biom J. 2020;62(7):1597–630.
Bakbergenuly I, Hoaglin DC, Kulinskaya E. Pitfalls of using the risk ratio in metaanalysis. Res Synth Methods. 2019;10(3):398–419.
Michael H, Thornton S, Xie M, Tian L. Exact inference on the randomeffects model for metaanalyses with few studies. Biometrics. 2019;75(2):485–93.
Günhan BK, Rover C, Friede T. Randomeffects metaanalysis of few studies involving rare events. Res Synth Methods. 2020;11(1):74–90.
Acknowledgements
None.
Funding
Open Access funding enabled and organized by Projekt DEAL. The Institute for Quality and Efficiency in Health Care (IQWiG) supported this work. No external funding was received.
Author information
Authors and Affiliations
Contributions
MF: design of study and simulation, programming and performing simulation, interpretation of results, writing manuscript. LB: design of study and simulation, interpretation of results, revision of manuscript. RB: design of study and simulation, interpretation of results, revision of manuscript. OK: development of methods, interpretation of results, revision of manuscript. GS: design of study and simulation, interpretation of results, revision of manuscript. TM: development of methods, design of study and simulation, programming and performing simulation, writing manuscript. The authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no financial competing interests. OK and TM were involved in previous studies on betabinomial models for metaanalyses.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1:
Data simulation and parameter estimation.
Additional file 2:
Table S1. Number of converged simulation runs (out of 10,000) for the odds ratio under the null hypothesis. Table S2. (Absolute) bias for the log odds ratio under the null hypothesis. Table S3. Coverage probability (%) for the odds ratio under the null hypothesis. Table S4. Length of 95% confidence interval for the log odds ratio under the null hypothesis. Table S5. Number of converged simulation runs (out of 10,000) for the odds ratio under the alternative hypothesis. Table S6. (Absolute) bias for the log odds ratio under the alternative hypothesis. Table S7. Percentage bias (%) for the log odds ratio under the alternative hypothesis. Table S8. Coverage probability (%) for the odds ratio under the alternative hypothesis. Table S9. Length of 95% confidence interval for the log odds ratio under the alternative hypothesis. Table S10. Power (%) for the odds ratio under the alternative hypothesis. Table S11. Number of converged simulation runs (out of 10,000) for the relative risk under the null hypothesis. Table S12. (Absolute) bias for the log relative risk under the null hypothesis. Table S13. Coverage probability (%) for the relative risk under the null hypothesis. Table S14. Length of 95% confidence interval for the log relative risk under the null hypothesis. Table S15. Number of converged simulation runs (out of 10,000) for the relative risk under the alternative hypothesis. Table S16. (Absolute) bias of the log relative risk under the alternative hypothesis. Table S17. Percentage bias (%) for the log relative risk under the alternative hypothesis. Table S18. Coverage probability (%) for the relative risk under the alternative hypothesis. Table S19. Length of 95% confidence interval for the log relative risk under the alternative hypothesis. Table S20. Power (%) for the relative risk under the alternative hypothesis.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Felsch, M., Beckmann, L., Bender, R. et al. Performance of several types of betabinomial models in comparison to standard approaches for metaanalyses with very few studies. BMC Med Res Methodol 22, 319 (2022). https://doi.org/10.1186/s12874022017793
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12874022017793
Keywords
 Betabinomial model
 Generalised linear mixed models
 Metaanalyses
 Simulation study
 Few studies