Testing homogeneity of effect sizes in pooling 2x2 contingency tables from multiple studies: a comparison of methods

Abstract

Meta-analysis is a statistical methodology that combines the outcomes of several independent studies. Either a fixed-effects or a random-effects model is used to combine the results of the independent studies. The choice depends on whether the outcomes can be considered homogeneous or not. Several methods have been developed to test the homogeneity of effect sizes. Although many of them have been compared already, they have not been studied in more realistic practical settings where individuals all have their own risk of an event and a complete overview is lacking. In this article, we investigate the performance of 10 statistical methods to test homogeneity of a binary exposure on a binary clinical outcome, using an extensive simulation study and a real life meta-analysis. We evaluated the Type I error and the statistical power. The fixed-effects regression model for treatment and study interaction was found to be a slightly liberal test, while the Q-statistic and the Bliss statistic can be considered conservative tests. The random-effects regression model, the Peto statistic and the I² perform rather poorly compared to the other methods. All chi-square tests that are based on the calculation of a specific odds ratio and the fixed-effects logistic regression analysis perform best. Among these chi-square tests, we recommend the Breslow–Day test, but only for convenience purposes, since it is available in most statistical software packages and the other tests are not.

Keywords:

• Bliss statistic
• Breslow-Day statistic
• Cochran’s Q-statistic
• generalized linear mixed model
• logistic regression
• multi-center trials
• Peto statistic
• Woolf statistic

PUBLIC INTEREST STATEMENT

Testing homogeneity of effect sizes (treatment effects, risk factor associations) is an essential part of meta-analysis since it helps us to interpret the consistency of these effects across groups of participants. The testing can be done using various statistical tools, which can result in different conclusions. For this reason we thought it was important to study the performance of the existing statistical tests for homogeneity. Specific comparisons have been done earlier but we also tried to explore the performance of statistical tests that have never been applied to this specific problem, and to adjust existing statistical tests. We studied the performance of 10 different methods to test the homogeneity of effect sizes using a case study and an extensive simulation study. Based on our study, we recommend one specific statistical test for use for meta-analysts when testing the homogeneity since it outperformed the others in most of the settings we studied.

1. Introduction

Pooling (analysis) results (like odds ratios) from multiple studies is common practice in medical, social and biological sciences. It is the main goal in an aggregate data meta-analysis or in an experimental study where some kind of stratification is applied (like multi-center randomized controlled trails). Researchers would typically be interested mainly in the overall or pooled estimate of a binary exposure or treatment effect, regardless of study differences, since it may better represent an estimate for the full target population than the individual studies. For a binary outcome, this means pooling 2 × 2 contingency tables from several studies. A traditional approach for pooling 2 × 2 contingency tables is the Cochran–Mantel–Haenszel statistic (Breslow, 1994; Breslow & Day, 1980; Mantel & Haenszel, 1959). However, this approach requires that the odds ratios across studies are homogeneous. This emphasizes the need for a test on homogeneity of the odds ratios or more generally the homogeneity of treatment effects on the binary clinical outcome. This topic is not new and many different approaches have been developed.

A classic test for homogeneity of effect sizes is the $Q$-statistic (Cochran, 1954). It creates a weighted effect size using the squared inverse within-study standard error and evaluates if the variation around this pooled effect is larger than that can be expected from the within-study standard errors. It was initially developed for pooling potency estimates from bioassay experiments (Mzolo, Hendriks, & van den Heuvel, 2013). A similar approach used in this area is the method of Bliss (1952), which uses the $Q$-statistic but in a different way. Related measures that are used in systematic reviews are the $H^2$, $I^2$ and the $R^2$ statistics (Higgins & Thompson, 2002). The $I^2$ statistic is also fully determined by the $Q$-statistic and it was initially developed to quantify the amount of heterogeneity instead of testing for homogeneity. The $H^2$ statistic is just an alternative measure for $I^2$, but the $R^2$ statistic describes the inflation from the random effects to the fixed effects, and it is viewed more as a measure than as a test statistic. However, the three measures of heterogeneity can be used as test statistics using appropriate confidence intervals.

Alternatively, fixed- and random-effects logistic regression models can be applied to test homogeneity. This is performed by including an interaction term between study and treatment, and then applying either a Wald test, a score test, or the Likelihood ratio test for its significance (Agresti, 2002).

A method that studies the heterogeneity of the odds ratio directly is the Breslow–Day test adjusted by Tarone (1985). It is some form of Pearson’s chi-square statistic. Other types of chi-square test statistics, which have been developed for specific situations, are Peto statistic (Yusuf, Peto, Lewis, Collins, & Sleight, 1986), Liang and Self statistics (Liang & Self, 1989) and Zelen statistic (Zelen, 1971). Yusuf et al. (1986) originally presented the Peto method for pooling odds ratios, and accompanied it with a homogeneity test (Yusuf et al., 1986). Liang and Self (1989) focused on developing test statistics in case data are sparse, i.e. some cells in the study-specific contingency tables contain only a few participants. Zelen developed the Zelen statistic to test the homogeneity of the relative risks (Zelen, 1971).

Several of these statistics have been discussed, evaluated and compared on their performance of detecting heterogeneity but, as Sutton and Higgins (2008) point out, there is no overview of the results (Sutton & Higgins, 2008). The $Q$-statistic has been shown to have low statistical power when the number of studies is small (Higgins & Thompson, 2002; Higgins, Thompson, Deeks, & Altman, 2003; Huedo-Medina, Sánchez-Meca, Marín-Martínez, & Botella, 2006) and to be conservative for a large number of studies (Higgins & Thompson, 2002). Moreover, the application of the $Q$-statistic to test the homogeneity of treatment effects has come under criticism recently (Hoaglin, 2016). Hoaglin (2016) points out that Cochran did not use the $Q$-statistic itself as a test statistic and that the $Q$-statistic only approaches a chi-square distribution under the null hypothesis asymptotically when within-study sample size increases (since the individual terms will become approximately normally distributed). In addition, the assumptions of normality of effect sizes and that the within-study variances are known do not apply to many meta-analyses. This is especially the case when the effect sizes are risk differences, the logarithm of the risk differences or the logarithm of the odds ratios. Hoaglin’s concerns also apply directly to the widely used $I^2$ as well (Hoaglin, 2016). Jackson (2006) derived a formula to calculate the power of the homogeneity test when using the $Q$-statistic. The author only considered the case when the number of studies is odd and the studies have the same within-variance, and concluded that the homogeneity test based on the $Q$-statistic indeed had low power. Kulinskaya, Dollinger and Bjørkestøl (2011) improved the performance of the $Q$-statistic using two approximations—a Gamma distribution based on expanded first two moments of the $Q$-statistic and a chi-square distribution with fractional degrees of freedom matching the first moment—and recommended using the latter approximation (Kulinskaya et al., 2011). Sánchez-Meca and Marín-Martínez (1997) compared the performance of the Q-statistic with the Schmidt and Hunter rule for homogeneity on a number of continuous outcomes (Sánchez-Meca & Marín-Martínez, 1997). The outcomes that were investigated were the standardized mean difference, the point-biserial correlation coefficient and the standardized correlation coefficient. The Schmidt and Hunter rule rejects the homogeneity of treatment effects if the sampling error variance is less than a specific threshold of the total variance (e.g. 75% or 90%). Sanchez-Meca et al. concluded that all tests had low statistical power, the Q-statistic resulted in nominal Type 1 errors while the Schmidt and Hunter rule overestimated the Type 1 errors (Sánchez-Meca & Marín-Martínez, 1997). Jones, O’Gorman, Lemke, Woolson and Monte Carlo (1989) examined the performances of the following tests: likelihood ratio test, Pearson’s chi-square test, Breslow–Day test (Jones et al., 1989) and Tarone’s adjustment to it (Tarone, 1985), a conditional score test and Liang and Self’s normal approximation to the score test (Liang & Self, 1989). They found that the likelihood ratio test was too liberal. They recommended using the Breslow–Day statistic when the number of studies is large, and Liang and Self’s normal approximation to the score test when the data are sparse, i.e. when there are few participants in some studies. Similar results were also reached by Paul and Donner (1992). Jones et al. (1989) concluded that all tests considered had low power when the number of studies is small (Jones et al., 1989). This was supported by Gavaghan, Moore, and McQuay (2000) who examined the Peto statistic (Yusuf et al., 1986), the Woolf statistic (Woolf, 1955), the Q-statistic (applied to the estimates of the risk difference) (DerSimonian & Laird, 1986), Liang and Self’s normal approximation to the score test, and the Breslow–Day test. These tests were also found to have incorrect Type 1 error rates in case of truly homogeneous sets. The authors concluded that the Q-statistic was too liberal and that the Woolf statistic was conservative, and they recommended the Breslow–Day test for routine use (Gavaghan et al., 2000). Recently, the Cochrane Reviews have started including $I^2$ in its reviews (Breslow, 1994). However, the $I^2$ statistic has been shown to be conservative and has a low power for continuous outcomes when the number of studies is small (Huedo-Medina et al., 2006). Bagheri, Ayatollahi and Jafari (2011) compared the performance of the random-effects logistic regression model (Bagheri et al., 2011) to Cochran’s Q-statistic and the Breslow–Day test, and found that the Breslow–Day test had the highest power. The authors recommended using the random-effects logistic regression model citing practical advantages, though these models are known for their computational difficulties (Wolfinger & O’ Connell, 1993).

This overview demonstrates that a lot of work has already been done. However, some measures have never been investigated (e.g. Bliss statistic), and not all measures have been investigated in all settings. More importantly, all studies assumed that the underlying binomial probabilities for the treatments in the stratified contingency tables would be satisfied and used this in their simulations, while in practical situations this would typically be violated. The aggregate data in a contingency table would be generated by individuals all having a different probability of obtaining the clinical outcome, due to other variables like age and gender. It is unclear whether the above-mentioned statistical methods will remain robust in more realistic practical settings. We use an extensive simulation study trying to account for hidden variations that can occur in practice, and we study the performance of 10 methods, including statistical methods that have not been applied before for this purpose.

This article is organized as follows. In Section 2 the statistical approaches for testing the homogeneity of odds ratios are presented, assuming that there does not exist individual risk. Section 3 describes the simulation model used where individual risk scores are being simulated. In Section 4 we present the simulation results and a case study on tuberculosis to illustrate the methods. The final section provides a discussion of the results and recommendation for use of the methods and possible further research.

2. Statistical models

Before we introduce the methods, we first introduce notation that will be used throughout this section. The binary exposure or treatment will be referred to as just treatment and the clinical outcome is considered binary. Then let $X_{hi}$ be the number of successes for treatment $h$ ($=0,1$) in study $i$ ($=1,2,...,m$) out of $n_{hi}$ trials. The control group is indexed by $h=0$ while the treatment group is indexed by $h=1$. The random variables $X_{h1},X_{h2},\dots ,X_{hm}$ are assumed to be mutually independent (independent studies). For the methods that test homogeneity of effect sizes, it is assumed that the expected value of $X_{hi}$ is equal to $E\left(X_{hi}\right)=n_{hi}\cdot p_{hi}$, with $p_{hi}$ being the proportion of a success for treatment $h$ in study $i$. All methods for testing homogeneity will assume that the proportion $p_{hi}$ satisfy the following form $p_{hi}=F\left({\alpha }_i+\left(\beta +{\gamma }_i\right)\cdot t_{hi}\right)$, with $F$ being some kind of a function (typically Logistic), ${\alpha }_i$ an intercept for study $i$, $\beta$ the (mean) treatment effect, $t_{hi}$ a treatment indicator variable for study $i$ with value 1 when $h=1$ and 0 otherwise, and ${\gamma }_i$ an interaction effect for study $i$. Then homogeneity means that

(1) $H_0:{\gamma }_1={\gamma }_2=\dots ={\gamma }_m=0$(1)

when the parameters ${\gamma }_{i}s$ are considered fixed or

(2) $H_0:\mathrm{v}\mathrm{a}\mathrm{r}\left({\gamma }_i\right)=0$(2)

when they are being assumed random.

Note that in practice the analysis of meta-analysis could also include one or more available covariates, e.g. average age per study or percentage males or females. In this setting, however, we assume that the only available information is the number of successes and the number of units in the treatment and control groups. This way the statistical methods discussed below can be applied to data with the least available information. Many of these methods can be easily generalized to accommodate extra covariates.

2.1. Logistic regression approaches

For fixed-effects logistic regression analysis, the function $F$ is a distribution function of the form $F(x)=exp(x)/\left[1+exp(x)\right]$. We also need to impose a constraint to obtain identifiability of the model. A natural choice is to take ${\gamma }_1=0$ or ${\gamma }_m=0$. Additionally, logistic regression implies that $X_{hi}$ is binomial and that $X_{0i}$ and $X_{1i}$ are independent (this is often not true due to different risks for individuals, but without additional information $X_{hi}$ is assumed binomial).

To test (1) we apply the likelihood ratio test. The likelihood of the full model for $F$ being the logistic distribution function is

$l\left({\alpha }_i,\beta ,{\gamma }_i\right)=\sum _{i=1}^m\left(X_i{\alpha }_i+X_{1i}{\beta }_i+log\left(1-F\left({\alpha }_i\right)\right)+log\left(1-F\left({\alpha }_i+{\beta }_i\right)\right)\right)$

with X_i = X_0i + X_1i being the total number of events and ${\beta }_i=\beta +{\gamma }_i$. The maximum likelihood estimates are denoted by $\hat{\alpha }{}_i^F$, $\hat{\beta }{}^F$ and $\hat{\gamma }{}_i^F$ for the full model and $\hat{\alpha }{}_i^0$ and $\hat{\beta }{}^0$ under the null hypothesis. The likelihood ratio test is now defined by

$\mathrm{L}\mathrm{R}\mathrm{T}=-2\left(l\left({\hat{\alpha }}_i^0,{\hat{\beta }}^0,0\right)-l\left({\hat{\alpha }}_i^F,{\hat{\beta }}^F,{\hat{\gamma }}_i^F\right)\right)$

The LRT has an approximate ${\chi }^2$ distribution with $m-1$ degrees of freedom [28]. The LRT statistic can be determined with standard software packages (we applied the GENMOD procedure of SAS).

Instead of a fixed-effects model, the interaction effects ${\gamma }_i$ can be taken random for the logistic link function. Under this assumption, it is common to also assume that the intercepts ${\alpha }_i$ are random. Thus, we will assume that ${\left({\alpha }_i,{\gamma }_i\right)}^T$ are bivariate normal with mean ${\left(\alpha ,0\right)}^T$ and variance–covariance matrix $\mathrm{\Sigma }$ given by

$\mathrm{\Sigma }=\left[\begin{array}{cc}{\sigma }_C^2& \rho {\sigma }_C{\sigma }_{CT}\\ \rho {\sigma }_C{\sigma }_{CT}& \sigma {}_{CT}^2\end{array}\right],$

with $\sigma {}_C^2$ being the standard deviation for the random intercepts, $\sigma {}_{CT}^2$ the standard deviation for the random interaction term and $\rho$ the correlation between ${\alpha }_i$ and ${\gamma }_i$.

The number of events $X_{hi}$, conditionally on ${\alpha }_i$ and ${\gamma }_i$, is considered binomial with parameters $n_{hi}$ and $p_{hi}=F\left({\alpha }_i+\left(\beta +{\gamma }_i\right)t_{hi}\right)$ and $X_{0i}$ and $X_{1i}$ are being independent given ${\alpha }_i$ and ${\gamma }_i$. The log likelihood of the full model can be written as:

$l\left({\alpha }_i,\beta ,{\gamma }_i\right)=\sum _{i=1}^m{log}_{\hspace{0.17em}}\underset{\mathbb{R}}{\int }\prod _{h=1}^2\left(\begin{array}{c}{n}_{hi}\\ X_{hi}\end{array}\right){\left[F\left(z_0+z_1{x}_{hi}\right)\right]}^{X_{hi}}{\left[1-F\left(z_0+z_1{x}_{hi}\right)\right]}^{n_{hi}-X_{hi}}{\mathrm{\Phi }}_2\left(z_0-\alpha ,z_1-\beta |\mathrm{\Sigma }\right)d{z}_{0}d{z}_1$

with ${\mathrm{\Phi }}_2\left(z_0,z_1|\sum \right)$ being the bivariate normal density given by

$\frac{1}{2\pi {\sigma }_C{\sigma }_{CT}\sqrt{1-{\rho }^2}}exp\left\{-\frac{1}{2\left(1-{\rho }^2\right)}\left[\frac{z{}_0^2}{\sigma {}_C^2}-2\rho \frac{z_0{z}_1}{{\sigma }_C{\sigma }_{CT}}+\frac{z{}_1^2}{\sigma {}_{CT}^2}\right]\right\}$

Homogeneity hypothesis in (2) can be evaluated using a likelihood ratio test again. The LRT becomes

$\mathrm{L}\mathrm{R}\mathrm{T}=-2\left(l\left(\hat{\alpha }{}_R^0,\hat{\beta }{}_R^0,\hat{\mathrm{\Sigma }}{}_R^0\right)-l\left(\hat{\alpha }{}_R^F,\hat{\beta }{}_R^F,\hat{\mathrm{\Sigma }}{}_R^F\right)\right)$

with $\hat{\alpha }{}_R^F,\hat{\beta }{}_R^F$ and $\hat{\mathrm{\Sigma }}{}_R^F$ being the maximum likelihood estimators under the full model and $\hat{\alpha }{}_R^0,\hat{\beta }{}_R^0$ and $\hat{\mathrm{\Sigma }}{}_R^0$ being the maximum likelihood estimators under the null hypothesis. Note that the variance–covariance matrix $\mathrm{\Sigma }$ becomes

$\mathrm{\Sigma }=\left[\begin{array}{cc}\sigma {}_C^2& 0\\ 0& 0\end{array}\right]$

when the null hypothesis would be true. The approximate asymptotic distribution of the LRT is a mixture of a chi-square distribution with two degrees of freedom and a degenerate distribution in zero. This means that we should compare the LRT statistic with a chi-square value at significance level $2\alpha$ (McCulloch, Searle, & Neuhaus, 2008; Verbeke & Molenberghs, 2003). The LRT statistic can be determined with standard software packages (we applied the NLMIXED procedure of SAS).

2.2. Approaches based on the Q-statistic

The Q-statistic was developed by Cochran (1954) and it can be applied to test the homogeneity of any type of effect size. The Q-statistic is given by $Q={\sum }_{i=1}^m{\left({\hat{\beta }}_i-\bar{\beta }\right)}^2/\hat{s}{e}_i^2$, with ${\hat{\beta }}_i$ being the estimated effect size, $\bar{\beta }$ a weighted average, and $\hat{s}e{}_i^2$ an estimated standard error for ${\hat{\beta }}_i$. For the log odds ratio of study $i$, we would have ${\hat{\beta }}_i=log\left(X_{1i}\left(n_{0i}-X_{0i}\right)/\left[\left(n_{1i}-X_{1i}\right)X_{0i}\right]\right)$; $\bar{\beta }$ is the weighted average which is given by

$\bar{\beta }=\sum _{i=1}^m\left({\hat{\beta }}_i/\hat{s}e{}_i^2\right)/\sum _{i=1}^m\left(1/\hat{s}e{}_i^2\right)$

with $\hat{s}e{}_i^2$ being given by

$\hat{s}e{}_i^2=\frac{1}{X_{1i}}+\frac{1}{n_{1i}-X_{1i}}+\frac{1}{X_{0i}}+\frac{1}{n_{0i}-X_{0i}}$

Cochran showed that the Q-statistic is approximately ${\chi }^2$ distributed with $m-1$ degrees of freedom, although he did not apply or discuss log odds ratios as effect sizes. This was done by Woolf (1955).

Bliss (1952) used the Q-statistic in a slightly different way to test the homogeneity hypothesis. Bliss’ test statistic is given by

$B=\left(m-1\right)+\sqrt{\left(\bar{n}-4)\right)/\left(\bar{n}-1\right)}\left\{\left(\bar{n}-2\right)\cdot Q/\bar{n}-\left(m-1\right)\right\}$

with $\bar{n}=\left({\sum }_{i=1}^m\left({\sum }_{h=0}^1{n}_{hi}-2\right)\right)/m$.

Note that Bliss’ method was originally developed for testing the equality of potencies obtained from multiple bioassays similar to that of Cochran. In such cases, $\bar{n}$ is the average of the corresponding degrees of freedom over all bioassays. The corresponding degrees of freedom per bioassay is obtained as a result of applying a parallel line bioassay model, i.e. a linear regression model (Mzolo et al., 2013). In our case we apply a fixed-effects logistic regression per study, so two parameters are estimated per study. Hence, we have chosen the degrees of freedom per study equal to ${\sum }_{h=0}^1{n}_{hi}-2$. Under the null hypothesis, the B-statistic may be approximated by a ${\chi }^2$ distributed random variable with $m-1$ degrees of freedom.

As mentioned in Section 1, Higgins and Thompson (2002) introduced a few measures for quantifying the heterogeneity of treatment effects that are based on the $Q$-statistic (Sutton & Higgins, 2008). The authors suggested a homogeneity test based on a confidence interval of the $I^2$ statistic, given by $I^2=100\cdot \left(H^2-1\right)/H^2$ with $H^2=Q/\left(m-1\right)$. A confidence interval for $H$ is given by $exp\left(log\left(H\right)\pm z_{\alpha /2}\cdot \hat{s}{e}_{log(H)}\right)$ with $z_{\alpha /2}$ being the $\alpha /2$ quantile of the standard normal distribution and

$\hat{s}{e}_{log(H)}=\left\{\begin{array}{cc}\frac{1}{2}\cdot \left[log(Q)-log(m-1)\right]/\left[\sqrt{2Q}-\sqrt{\left(2m-3\right)}\right]& \mathrm{i}\mathrm{f}Q>m\\ \sqrt{{\left[2(m-2)\right]}^{-1}\left[3{\left(m-2\right)}^2-1\right]/\left[3{\left(m-2\right)}^2\right]}& \mathrm{i}\mathrm{f}Q\le m\end{array}\right.$

This confidence interval for $H$ can be transformed to a confidence interval for $I^2$, and the transformed confidence interval can be used to test hypothesis in (1) by verifying whether the value zero lies in the interval or not.

2.3. Chi-square approaches

One of the oldest approaches for testing homogeneity of odds ratio is the Breslow–Day test (Breslow & Day, 1980). They used the Cochran–Mantel–Haenszel pooled odds ratio $\mathrm{O}{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}$ calculated by

$O^{\wedge }{R}_{CMH}=\frac{\sum _{i=1}^m\left(X_{1i}\cdot \left(n_{0i}-X_{0i}\right)/n_i\right)}{\sum _{i=1}^m\left(X_{0i}\cdot \left(n_{1i}-X_{1i}\right)/n_i\right)}$

and assumed the same distributional assumptions for $X_{hi}$ as for the fixed-effects logistic regression analysis. Now define $E_{1i}=E\left(X_{1i}|X_i,{\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}\right)$ as the expected value of $X_{1i}$ given $X_i=X_{0i}+X_{1i}$ and ${\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}$ under the assumption that the null hypothesis in (1) is correct. $E_{1i}$ can be obtained by solving the following quadratic equation:

(3) $\left({\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}-1\right)E{}_{1i}^2-\left(\left(X_i+n_{1i}\right){\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}+\left(n_{0i}-X_i\right)\right)E_{1i}+X_i{n}_{1i}{\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}=0$(3)

For testing the homogeneity, the Tarone test statistic (Tarone, 1985) is used. This statistic is given by

(4) $\sum _m^{i=1}\frac{{\left(X_{1i}-E_{1i}\right)}^2}{\mathrm{v}\mathrm{a}\mathrm{r}\left(X_{1i}|X_i,{\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}\right)}-\frac{{\left(X_1-E_1\right)}^2}{\sum _m^{i=1}\mathrm{v}\mathrm{a}\mathrm{r}\left(X_{1i}|X_i,{\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}\right)}$(4)

with $X_1=\sum _{i=1}^m{X}_{1i}$, $E_1=\sum _{i=1}^m{E}_{1i}$ and $var\left(X_{1i}|X_i,\mathrm{O}{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}\right)$ is given by

(5) $var\left(X_{1i}|X_i,{\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}\right)={\left(\frac{1}{E_{1i}}+\frac{1}{X_i-E_{1i}}+\frac{1}{n_{1i}-E_{1i}}+\frac{1}{n_{0i}-X_i+E_{1i}}\right)}^{-1}$(5)

Note that we applied the Tarone statistic instead of Breslow–Day statistic. The Breslow–Day statistic is based on the Cochran–Mantal–Haenszel odds ratio estimator. This estimator has been shown to be consistent but not efficient. For this reason Tarone modified the Breslow–Day statistic (Tarone, 1985). If (1) is true, the statistic in (4) approximately follows a ${\chi }^2$ distribution with $m-1$ degrees of freedom. The Tarone-adjusted Breslow–Day test can be carried out using standard software packages (we applied the procedure FREQ of SAS).

Zelen (1971) also presented a statistic to test the homogeneity of odds ratios (Zelen, 1971). However, it was found to be biased and inconsistent (Halperin et al., 1977). Instead, Halperin et al. (1977) presented a corrected version of the Zelen statistic, which we will use. It also uses an expected value $E_{1i}$, like Breslow–Day, but now with an alternative estimate of the pooled odds ratio. The odds ratio $O^{\wedge }{R}_{CZ}$ is given by $O^{\wedge }{R}_{CZ}=exp\left({\hat{\beta }}^0\right)$ with ${\hat{\beta }}^0$ being the maximum likelihood estimator of the Fixed-effects logistic regression analysis under the null hypothesis of homogeneity. The expected value $E_{1i}=E\left(X_{1i}|X_i,O^{\wedge }{R}_{CZ}\right)$ can be obtained by solving Equation (3) with ${\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}$ being replaced by $O^{\wedge }{R}_{CZ}$.

Then the corrected Zelen statistic for testing homogeneity is given by

$Z=\sum _{i=1}^m{\left(X_{1i}-E_{1i}\right)}^2/var\left(X_{1i}|X_i,O^{\wedge }{R}_{CZ}\right)$

with the variance given in (5) using the solution $E_{1i}=E\left(X_{1i}|X_i,O^{\wedge }{R}_{CZ}\right)$. If the null hypothesis in (1) is true, the corrected $Z$ statistic will follow an asymptotic ${\chi }^2$ distribution with $m-1$ degrees of freedom (Halperin et al., 1977).

Liang and Self (1989) presented test statistics for sparse data (Liang & Self, 1989). (Recall that sparse data have low cell counts, thus not necessarily low number of studies.) Again, an expected value $E_{1i}$ is used given $X_i$ and a pooled odds ratio $O^{\wedge }{R}_{CL}$. Here $\hat{O}{R}_{CL}$ is given by $O^{\wedge }{R}_{CL}=exp\left(\hat{\beta }{}_{CL}^0\right)$ with $\hat{\beta }{}_{CL}^0$ being the maximum likelihood estimator of the conditional likelihood function given $X_i$ and ${\gamma }_1={\gamma }_2=\dots ={\gamma }_m=0$. The conditional logistic regression can be obtained using standard software packages. (We applied the procedure LOGISTIC of SAS.) The expected value $E_{1i}=E\left(X_{1i}|X_i,O^{\wedge }{R}_{CL}\right)$ can be obtained by solving Equation (3) with $O^{\wedge }{R}_{CMH}$ being replaced by $O^{\wedge }{R}_{CL}$. Liang and Self test statistic $T_2$ is defined by

$T_2=\sum _{i=1}^m{\left(X_{1i}-E_{1i}\right)}^2/var\left(X_{1i}|X_i,O^{\wedge }{R}_{CL}\right)$

with $var\left(X_{1i}|X_i,O^{\wedge }{R}_{CL}\right)$ being given by (5) and $E_{1i}$ being the described solution of (3) using $O^{\wedge }{R}_{CL}$. The statistic $T_2$ is approximately chi-squared distributed with $m-1$ degrees of freedom (Liang & Self, 1989).

Following Liang and Self’s approach, we derived an alternative statistic. Instead of using the conditional logistic regression, we apply a random-effects logistic regression model, assuming that $var\left({\gamma }_i\right)=0$. The odds ratio we apply is $O^{\wedge }{R}_R=exp\left(\hat{\beta }{}_R^0\right)$ with $\hat{\beta }{}_R^0$ being the maximum likelihood estimator described in Section 2.1. The test statistic $T_R$ is equal to $T_2$ but with $E_{1i}$ the solution of (3) with ${\mathrm{O}}^{\wedge }{\mathrm{R}}_{\mathrm{C}\mathrm{M}\mathrm{H}}$ replaced by $O^{\wedge }{R}_R$ and the variance given by (4), with $E_{1i}$ the described solution of (3) using $O^{\wedge }{R}_R$.

Finally, a complete alternative ${\chi }^2$ statistic is the Peto statistic. Yusuf et al. (1986) presented the Peto statistic to test the homogeneity of the odds ratios [10]. The Peto statistic $P$ is given by

$P=\sum _{i=1}^m\frac{{\left(O_i-E_i\right)}^2}{V_i}-\frac{{\left({\sum }_{i=1}^m\left(O_i-E_i\right)\right)}^2}{{\sum }_{i=1}^m{V}_i}$

where $O_i=X_{1i},$ $E_i=X_i{n}_{1i}/n_i$, $n_i=n_{1i}+n_{0i}$ and $V_i$ defined by

$V_i=\frac{X_i\cdot \left(n_i-X_i\right)\cdot n_{1i}\cdot n_{0i}}{n_i^2\cdot \left(n_i-1\right)}$

When the null hypothesis in (1) is true, the Peto statistic will have an approximate ${\chi }^2$ distribution with $m-1$ degrees of freedom (Yusuf et al., 1986).

3. Simulation

3.1. Introduction

Although all approaches discussed above assume equal risk probabilities for all participants in one study for one treatment group, this assumption is hardly true in practice. In reality, risk probabilities of individuals will be influenced by individual’s characteristics. Therefore, the data are simulated so that each individual has its unique risk probability. This is done by creating covariates for each individual that affect the risk probability in addition to the influence of the treatment. Since populations across studies vary, the distributions of the covariates are affected by studies, causing already one form of treatment–study interaction. We also added treatment–study interactions independent of the covariates.

3.2. Simulation model and settings

The simulation model to generate $m$ studies can be described as follows. First, the total number of participants $n_i$ for study $i$ is drawn. We used an overdispersed Poisson distribution with parameter $\lambda$. The value ${\gamma }_i\sim \mathrm{\Gamma }\left(a_0,b_0\right)$ is drawn to make a study-specific parameter ${\lambda }_i=\lambda \cdot exp\left(0.5\cdot {\gamma }_i\right)$. Then $n_i$ is drawn from a Poisson distribution $\left({\lambda }_i\right)$. As mentioned in Section 1, to make the simulation somewhat more realistic we generated covariates, in this case age and gender. The covariate age $a_{ij}$ for the jth participant in study $i$ is created by

$a_{ij}={\alpha }_0+{\eta }_i+{\epsilon }_{ij},$

with ${\alpha }_0$ being the average age over all participants in all studies, ${\eta }_i\sim N\left(0;\tau {}_s^2\right)$ a study effect for age and ${\epsilon }_{ij}\sim N\left(0,\tau {}_{a(s)}^2\right)$ the age effect within study $i$. The covariate gender $g_{ij}$ for the $j{}^{th}$ participant in study i is generated with a uniform variable $z_{ij}$. The jth participant will be assigned the value $g_{ij}=1$ if $z_{ij}>0.5$ and zero otherwise.

Then we need to generate study-specific treatment effects independent of covariates. We use a bivariate normally distributed variable $\left(u_{0i},u_{1i}\right)$ such that

$\left[\begin{array}{c}{u}_{0i}\\ u_{1i}\end{array}\right]\sim N\left(\left[\begin{array}{c}\alpha \\ \alpha +\beta \end{array}\right];\left[\begin{array}{cc}\sigma {}_0^2& \rho {\sigma }_0{\sigma }_1\\ \rho {\sigma }_0{\sigma }_1& \sigma {}_1^2\end{array}\right]\right)$

The parameter $\alpha$ represents an intercept to generate the risk of an event for the control group and $\beta$ is the direct treatment effect. Additionally, the random variable $v_{ij}$ is drawn from a normal distribution, $v_{ij}\sim N(0,\sigma {}^2)$, to generate heterogeneity across participants in the study due to variables other than age and gender.

A treatment indicator ${\delta }_{ij}$ is generated by a uniform $\left(0,1\right)$ distribution. If $w_{ij}$ is uniformaly distributed, then ${\delta }_{ij}=1$ if $w_{ij}>p$ and zero otherwise, with $p$ being the ratio of allocation of control and treatment. The event $X_{ij}$ is Bernoulli $\left(p_{ij}\right)$ distributed with $p_{ij}$ satisfying

(6) $log\left[p_{ij}/\left(1-p_{ij}\right)\right]=\left({\alpha }_a+{\beta }_a\cdot {\delta }_{ij}\right)\cdot a_{ij}+{\alpha }_g\cdot g_{ij}+u_{ij}+v_{ij}$(6)

with ${\alpha }_a$ being the effect of age on the risk of an event irrespective of treatment, ${\beta }_a$ a treatment age interaction effect, ${\alpha }_g$ the effect of gender and $u_{ij}=u_{01}\left(1-{\delta }_{ij}\right)+u_{11}{\delta }_{ij}$ the study treatment effect. Finally, the dichotomous response variable $X_{ij}$ is drawn with a uniform $\left(0,1\right)$ distribution, i.e. $X_{ij}=1$ if $y_{ij}>p_{ij}$ and zero otherwise, with $y_{ij}\sim U\left(0,1\right)$. This procedure is repeated 1,000 times resulting in 1,000 simulation runs for any set of parameters. As per simulation run, it is determined whether the null hypothesis of homogeneity is rejected or not using the approaches in Section 2 which all ignore the individual heterogeneities. The percentage of rejected null hypotheses over the number of simulations is being reported.

Note that the final model in $\left(6\right)$ has a study-specific latent variable $Z_{0i}={\alpha }_a{\eta }_i+u_{0i}$ for the control group, and a study-specific latent variable $Z_{1i}=\left({\alpha }_a+{\beta }_a\right){\eta }_i+u_{1i}$ for the treatment group. The variance–covariance matrix for these latent variables is given by

$\left[\begin{array}{cc}\alpha {}_a^2\tau {}_s^2+\sigma {}_0^2& {\alpha }_a\left({\alpha }_a+{\beta }_a\right)\tau {}_s^2+\rho {\sigma }_0{\sigma }_1\\ {\alpha }_a\left({\alpha }_a+{\beta }_a\right)\tau {}_s^2+\rho {\sigma }_0{\sigma }_1& {\left({\alpha }_a+{\beta }_a\right)}^2\tau {}_s^2+\sigma {}_1^2\end{array}\right]$

when ${\beta }_a=0$, $\rho =1$ and ${\sigma }_0={\sigma }_1$; the correlation between $Z_{0i}$ and $Z_{1i}$ is equal to 1 and there is no treatment–study interaction. In case ${\beta }_a\ne 0$ or $\rho \ne 1$, there is a treatment–study interaction and thus heterogeneity in treatment effects. The size of the heterogeneity is given by $\beta {}_a^2\tau {}_s^2+2\sigma {}_0^2\left(1-\rho \right),$ when ${\sigma }_0={\sigma }_1$. The larger this value, the higher the statistical power of the methods. Thus, the heterogeneity reduces when $\rho$ increases to 1. For the simulation, a choice was made to take $\sigma {}_0^2=\sigma {}_1^2$. Type I error was simulated assuming $\rho =1$ and ${\beta }_a=0$. The different simulation settings that we studied are presented in Table 1.

4. Results

4.1. Meta-analysis of antituberculosis therapy

We applied the statistical methods presented in Section 2 to a meta-analysis on antituberculosis therapy (Pasipanodya, Srivastava, & Gumbo, 2012). This meta-analysis consisted of 26 different regimens (plans or regulated courses). In each regimen, slow acetylators (control) and rapid acetylators (treatment) were applied to two groups of patients. The objective was to test whether microbiological failure after therapy, relapse after therapy and acquired drug resistance (ADR) were as likely to occur among slow and rapid acetylators. The authors limited their selection of studies to randomized, controlled clinical trials in which outcomes were reported in table, graph or text format. All patients who were considered microbiological failures, relapses or ADR during treatment for slow and rapid acetylators were recorded. For each of the three different outcomes, the authors then calculated risk ratios (RR) with 95% confidence intervals using regimen as unit of analysis. Random-effects models were used to combine regimens when $I^2>50\mathrm{\%}$ and the corresponding P-value <0.05; otherwise fixed-effects models were used to combine the RR.

For our study we only consider the results of one outcome, namely microbiological failure. The authors performed a test of homogeneity of the relative risks using the $I^2$ and found $I^2=39\mathrm{\%}$ with a corresponding $P=0.035$ (Pasipanodya et al., 2012). The authors then applied mixed-effects modeling for this outcome, resulting in a combined RR of 2, with a confidence interval [1.5 ; 2.7] for the RR.

However, four of the 26 regimens contained zero cells, and the selected methods cannot deal with this unless a correction to the number of cells is made (Gart, 1966). We chose to remove these four regimens from our data, and apply all the selected approaches to the remaining datasets.

The test statistic, degrees of freedom and the P-value for rejecting the null hypothesis are presented in Table 2. Using 5% level of significance, the null hypothesis is only rejected by the random-effects logistic regression model. Additionally, the P-values range from 0.035 to 0.124. Thus, the proposed methods do not all seem to agree.

4.2. Simulation results

The full set of results is presented in the supplementary material. A subset of these results is shown in Tables 3–6. Tables 3–5 present a subset of the simulation results where $p_{hi}$ depends on age, gender and other heterogeneous factors $\left(v_{ij}\right)$. They represent the results for the most realistic practical settings to contrast these results and to be able to compare our results with literature. Table 6 presents a selection of the Type I error and the statistical power results of the test methods under the assumption of the binomial distribution. In such case, the proportion of a success, $p_{hi}$, does not depend on age, gender or any other participant variables $\left(v_{ij}\right)$. The random-effects logistic regression analysis occasionally showed numerical problems, in particular when heterogeneity between studies was small or not present at all. Note that in around 13.7% of all simulation runs the full random-effects logistic regression model could not converge and the results of these problematic simulation runs are eliminated from the results.

Table 1. Parameter values used for a simulation model

Number of studies	$m$	5 ; 10 ; 15 ; 20
Poisson distribution parameter for a number of participants per study	$\lambda$	100 ; 200
Overdispersion on a number of participants	$a_0=b_0$	1
Mean age	${\alpha }_0$	65
Between-study variance in age	${\tau }_s^2$	10
Within-study variance in age	${\tau }_{a(s)}^2$	15
Heterogeneity in treatment effect	$\rho$	1 ; 0.3 ; 0.6
Allocation proportion of treatment	$p$	0.3 ; 0.5
Overall risk intercept	$\alpha$	−6
Age effect on risk	${\alpha }_a$	0.05
Gender effect on risk	${\alpha }_g$	0.1
Treatment effect	$\beta$	0 ; 2
Age-treatment interaction effect	${\beta }_a$	0 ; 0.01
Between-study variance on risk for treatment $h$	${\sigma }_h^2$	0.2 ; 0.5

Table 2. Results of meta-analysis on the pharmacokinetic variability hypothesis for ADR and failure of antituberculosis therapy

	Test statistic	DF	P-value
FE	33.31	21	0.124
RE	5.33	2	0.035
Q	29.92	21	0.094
B	29.32	21	0.107
$I^2$	[0 ; +58.30]	—	>0.05
BD	32.23	21	0.055
Z	32.27	21	0.055
$T_2$	32.24	21	0.055
$T_R$	32.28	21	0.055
P	31.73	21	0.062

Note that we could not reject the heterogeneity in treatment effects with $I^2$, as the authors suggested. This did not change when we included the four studies and changed the cell counts with 0.5.

Table 3. Type I error rate for $m=5$, $a_a=0.05$, $a_g=0.1$

$\lambda$	$\beta$	${\sigma }_{\mathit{h}}^{\mathbf{2}}$	$p$	FE	RE	Q	B	$I^2$	BD	Z	${\mathit{T}}_{\mathbf{2}}$	${\mathit{T}}_{\mathit{R}}$	P
100	0	0.2	0.3	0.060	0.024	0.028	0.027	0.020	0.045	0.044	0.045	0.045	0.042
			0.5	0.060	0.033	0.039	0.037	0.025	0.054	0.054	0.054	0.054	0.050
		0.5	0.3	0.060	0.059	0.020	0.015	0.010	0.040	0.040	0.040	0.041	0.036
			0.5	0.064	0.054	0.030	0.024	0.014	0.050	0.050	0.050	0.050	0.046
	2	0.2	0.3	0.056	0.035	0.022	0.021	0.014	0.044	0.044	0.040	0.044	0.005
			0.5	0.045	0.032	0.032	0.028	0.021	0.037	0.038	0.037	0.039	0.019
		0.5	0.3	0.059	0.082	0.014	0.012	0.007	0.046	0.044	0.044	0.043	0.003
			0.5	0.052	0.086	0.031	0.028	0.019	0.044	0.045	0.043	0.044	0.015
200	0	0.2	0.3	0.057	0.023	0.044	0.043	0.024	0.050	0.050	0.050	0.050	0.050
			0.5	0.045	0.025	0.034	0.031	0.022	0.040	0.040	0.040	0.040	0.039
		0.5	0.3	0.051	0.070	0.032	0.031	0.020	0.043	0.043	0.043	0.043	0.042
			0.5	0.049	0.099	0.034	0.033	0.023	0.048	0.048	0.048	0.048	0.047
	2	0.2	0.3	0.055	0.036	0.037	0.032	0.014	0.045	0.045	0.045	0.045	0.004
			0.5	0.036	0.040	0.028	0.028	0.017	0.031	0.034	0.033	0.033	0.017
		0.5	0.3	0.052	0.106	0.031	0.030	0.015	0.046	0.050	0.049	0.050	0.008
			0.5	0.050	0.111	0.037	0.032	0.020	0.043	0.043	0.043	0.043	0.013

Table 4. Statistical power for $\rho =0.3$, $m=5$, $l=100$, $a_a=0.05$, $a_g=0.1$

$\beta$	${\beta }_{\mathit{a}}$	${\sigma }_{\mathit{h}}^2$	$\mathit{p}$	FE	RE	Q	B	$I^2$	BD	Z	${\mathit{T}}_{\mathbf{2}}$	${\mathit{T}}_{\mathit{R}}$	P
0	0	0.2	0.3	0.234	0.154	0.174	0.160	0.135	0.223	0.224	0.223	0.225	0.222
			0.5	0.269	0.153	0.225	0.210	0.173	0.258	0.258	0.258	0.258	0.250
		0.5	0.3	0.501	0.367	0.359	0.344	0.299	0.481	0.481	0.481	0.480	0.460
			0.5	0.556	0.429	0.469	0.444	0.402	0.547	0.547	0.547	0.547	0.532
	0.01	0.2	0.3	0.252	0.158	0.192	0.172	0.142	0.237	0.242	0.241	0.242	0.186
			0.5	0.315	0.193	0.272	0.255	0.206	0.301	0.302	0.304	0.302	0.274
		0.5	0.3	0.520	0.423	0.417	0.401	0.347	0.498	0.500	0.500	0.500	0.400
			0.5	0.589	0.465	0.522	0.512	0.465	0.578	0.577	0.577	0.577	0.548
2	0	0.2	0.3	0.264	0.167	0.203	0.188	0.158	0.258	0.267	0.263	0.255	0.116
			0.5	0.323	0.210	0.284	0.275	0.234	0.315	0.318	0.315	0.316	0.212
		0.5	0.3	0.540	0.451	0.419	0.407	0.357	0.510	0.520	0.517	0.520	0.334
			0.5	0.626	0.511	0.570	0.550	0.514	0.617	0.623	0.621	0.621	0.499
	0.01	0.2	0.3	0.255	0.161	0.205	0.188	0.153	0.260	0.266	0.264	0.264	0.110
			0.5	0.303	0.190	0.276	0.260	0.221	0.301	0.306	0.304	0.307	0.173
		0.5	0.3	0.536	0.405	0.418	0.403	0.363	0.507	0.521	0.518	0.520	0.343
			0.5	0.628	0.503	0.575	0.560	0.514	0.623	0.629	0.652	0.624	0.467

Table 5. Statistical power for $r=0.6$, $m=5$, $l=100,$$a_a=0.05$, $a_g=0.1$

$\beta$	${\beta }_{\mathit{a}}$	${\sigma }_{\mathit{h}}^{\mathbf{2}}$	$\mathit{p}$	FE	RE	Q	B	$I^2$	BD	Z	${\mathit{T}}_{\mathbf{2}}$	${\mathit{T}}_{\mathit{R}}$	P
0	0	0.2	0.3	0.166	0.092	0.118	0.111	0.087	0.159	0.158	0.158	0.159	0.158
			0.5	0.181	0.097	0.140	0.128	0.106	0.165	0.165	0.165	0.165	0.161
		0.5	0.3	0.352	0.245	0.229	0.213	0.175	0.337	0.337	0.337	0.336	0.330
			0.5	0.383	0.265	0.290	0.273	0.232	0.356	0.356	0.355	0.357	0.344
	0.01	0.2	0.3	0.176	0.097	0.121	0.113	0.086	0.159	0159	0.159	0.159	0.117
			0.5	0.203	0.129	0.169	0.157	0.127	0.193	0.194	0.193	0.193	0.174
		0.5	0.3	0.351	0.288	0.246	0.235	0.195	0.336	0.337	0.337	0.337	0.255
			0.5	0.422	0.309	0.350	0.331	0.292	0.401	0.402	0.402	0.402	0.362
2	0	0.2	0.3	0.179	0.093	0.118	0.106	0.082	0.164	0.167	0.166	0.165	0.047
			0.5	0.209	0.123	0.173	0.162	0.127	0.199	0.200	0.200	0.199	0.115
		0.5	0.3	0.369	0.300	0.254	0.243	0.208	0.339	0.347	0.346	0.347	0.183
			0.5	0.458	0.328	0.391	0.376	0.318	0.442	0.446	0.444	0.443	0.301
	0.01	0.2	0.3	0.180	0.093	0.121	0.110	0.083	0.174	0.178	0.177	0.177	0.046
			0.5	0.182	0.107	0.157	0.145	0.121	0.172	0.175	0.175	0.174	0.089
		0.5	0.3	0.370	0.276	0.250	0.239	0.208	0.330	0.342	0.340	0.341	0.191
			0.5	0.442	0.312	0.380	0.367	0.324	0.427	0.435	0.432	0.432	0.255

Table 6. Type I error and statistical power for $a_a=a_g={\sigma }^2=0$, $m=5$, $l=100$

$\rho$	$\beta$	${\sigma }_{\mathit{h}}^{\mathbf{2}}$	$\mathit{p}$	FE	RE	Q	B	${\mathit{I}}^{\mathbf{2}}$	BD	Z	${\mathit{T}}_{\mathbf{2}}$	${\mathit{T}}_{\mathit{R}}$	P
1	0	0.2	0.3	0.047	0.072	0.041	0.035	0.024	0.048	0.048	0.048	0.048	0.047
			0.5	0.050	0.084	0.045	0.041	0.028	0.049	0.049	0.049	0.049	0.046
		0.5	0.3	0.048	0.135	0.035	0.034	0.018	0.041	0.041	0.041	0.041	0.039
			0.5	0.046	0.142	0.038	0.035	0.024	0.043	0.043	0.043	0.043	0.042
	2	0.2	0.3	0.054	0.035	0.030	0.028	0.018	0.045	0.047	0.046	0.046	0.056
			0.5	0.058	0.041	0.026	0.023	0.011	0.049	0.052	0.050	0.051	0.004
		0.5	0.3	0.058	0.073	0.030	0.025	0.019	0.049	0.048	0.048	0.048	0.059
			0.5	0.058	0.084	0.027	0.021	0.017	0.042	0.044	0.044	0.044	0.011
0.3	0	0.2	0.3	0.596	0.471	0.583	0.570	0.531	0.595	0.595	0.595	0.595	0.588
			0.5	0.656	0.510	0.646	0.632	0.585	0.651	0.651	0.651	0.651	0.645
		0.5	0.3	0.828	0.774	0.821	0.807	0.785	0.828	0.829	0.829	0.828	0.821
			0.5	0.876	0.791	0.869	0.864	0.839	0.874	0.874	0.874	0.874	0.869
	2	0.2	0.3	0.443	0.252	0.390	0.370	0.320	0.432	0.437	0.437	0.435	0.422
			0.5	0.425	0.299	0.341	0.330	0.272	0.401	0.411	0.408	0.410	0.211
		0.5	0.3	0.726	0.587	0.659	0.641	0.596	0.716	0.718	0.716	0.717	0.699
			0.5	0.727	0.602	0.626	0.607	0.568	0.706	0.712	0.710	0.710	0.534
0.6	0	0.2	0.3	0.428	0.314	0.408	0.392	0.346	0.426	0.426	0.426	0.426	0.416
			0.5	0.479	0.336	0.461	0.438	0.393	0.471	0.471	0.471	0.471	0.464
		0.5	0.3	0.697	0.623	0.680	0.663	0.630	0.695	0.695	0.695	0.695	0.695
			0.5	0.761	0.659	0.752	0.741	0.711	0.758	0.758	0.758	0.758	0.750
	2	0.2	0.3	0.301	0.162	0.248	0.235	0.193	0.281	0.284	0.284	0.284	0.276
			0.5	0.277	0.187	0.210	0.200	0.164	0.254	0.262	0.260	0.260	0.111
		0.5	0.3	0.553	0.404	0.475	0.456	0.408	0.534	0.539	0.536	0.535	0.514
			0.5	0.553	0.420	0.457	0.441	0.396	0.535	0.540	0.536	0.537	0.359

For the realistic practical settings (Table 3), the fixed-effects logistic regression analysis has a slight increased Type I error, when the number of participants within studies is relatively small, irrespective of the number of studies, but this is still acceptable. When sample sizes within studies increase, the Type I error is closer to nominal. The Q-statistic, the Bliss statistic and the $I^2$ statistic constantly underestimate the Type I error (Table 3). This is most apparent when only a third of the participants are exposed to the treatment. This gets slightly better when the number of studies increases (supplementary material). The $I^2$ statistic has the smallest Type I error rates of all the methods, except for a few cases where the Peto statistic gives even smaller values. The random-effects logistic regression model mostly underestimates the Type I error rate when ${\sigma }_0^2$ and ${\sigma }_1^2$ are small. However as ${\sigma }_0^2$ and ${\sigma }_1^2$ increase, it overestimates the Type I error rate (Table 2). This gets even worse when the number of studies increases (supplementary material). The Type I error rates of the chi-square methods using estimates of the odds ratio, i.e. the Breslow–Day statistic, the corrected Zelen statistic, Liang and Self’s $T_2$ and its adjusted version $T_R$, are almost at the nominal value. However, when one-third of the participants is exposed to treatment, the chi-square methods stay just below the nominal level, but this disappears when the number of studies increases. When participants no longer have their individual risk, the Type I error rates (Table 5) are mostly similar to the Type I error rates for the realistic practical settings, except for two cases. The fixed-effects logistic regression no longer gives liberal Type I error rates when there is no treatment effect $\left(\beta =0\right)$.

Investigating the power for the realistic practical settings (Tables 4 and 5), it obviously shows that the statistical power increases as the number of studies increases. This improved power for testing the interaction effect can be explained by the reduced sampling error due to increased sample size. Furthermore, the increase in ${\sigma }_0^2$ and ${\sigma }_1^2$ results in higher power since they increase the treatment interaction effects. This is due to the term $2{\sigma }_0^2\left(1-\rho \right)$ in the size of heterogeneity. The statistical power is negatively correlated with $\rho$, with no interaction effect when $\rho =1$ and ${\beta }_a=0$ (as mentioned earlier). The highest statistical power is obtained by the fixed-effects logistic regression model, followed by the ${\chi }^2$ statistics (Breslow–Day test, the corrected Zelen test, $T_2$ and $T_R$). The lowest statistical power is obtained by the $I^2$ and the random-effects logistic regression model, respectively. These results hold true irrespective of the number of studies. Although Bliss uses the Q-statistic differently from Cochran, the two tests behave almost identically but they are less powerful than the best test statistics. It should be noted that heterogeneity in treatment effects from differences in age effects across studies $\left({\beta }_a>0\right)$ is not increasing the power a lot, a result that holds even when the number of studies increases. This is due to the selected settings. The heterogeneity effect in the simulation was only increased with ${\beta }_a^2{\tau }_a^2={\left(0.01\right)}^2\cdot 10=0.001$, which is really small.

Comparing the power results from the setting that participants have a common risk (Table 6) with the power results from the setting that individuals have their own individual risks (Tables 4 and 5) shows a few interesting differences. First, the power is substantially smaller when individuals are heterogeneous in their risk rates, which makes it more difficult to detect treatment heterogeneity across studies. The power values for testing treatment interactions increase slightly when there would be a treatment effect in case individuals have heterogeneous risk rates. However, this is opposite when individuals are not heterogeneous, since the power values for treatment heterogeneity is substantially smaller for $\beta =2$ than for $\beta =0$ (Table 6 and supplementary material). With respect to choosing the best methods for detecting treatment effect heterogeneity, no real difference could be detected.

5. Discussion

The purpose of this article was to investigate the performance of statistical methods that may test for the homogeneity of odds ratios in meta-analysis. Ten methods were applied to a real-life meta-analysis study and to practically realistic simulated datasets. The simulated datasets were created using different scenarios based on different numbers of studies, participants within studies, unbalanced number of observations for treatment and control groups, whether or not there is an actual treatment effect, individual risks through covariates and the level of heterogeneity.

The Peto statistic is somewhat unreliable with respect to its Type I error rate, a conclusion that has been reached earlier (Gavaghan et al., 2000). The fixed-effects logistic regression analysis is a slightly liberal test, although the inflated Type I error rate is not too large, and it reduces with larger number of studies. This conclusion has been reached earlier (Paul & Donner, 1992), although we could not demonstrate it for settings where all individuals have the same constant risk rates, the settings used in the literature. In our study the inflation happened when the proportion of success depends on individual covariates and with small numbers of studies. This interesting result merits further investigation to find out the exact analytical cause, but it is beyond the scope of this article.

The fact that the fixed-effects logistic regression has the highest power among all other methods is partly due to this inflated Type I error rate. The $I^2$ statistic can be considered a very conservative test. This conclusion had been reached earlier (Huedo-Medina et al., 2006), but only for continuous responses. It makes the $I^2$ an unsuitable test statistic for homogeneity of odds ratios. The random-effects logistic regression analysis is also conservative, albeit less conservative than the $I^2$ for small numbers of studies, but it is less powerful than $I^2$ for a larger number of studies. This low power combined with the unreliable Type I error rate makes the random-effects model unsuitable for testing homogeneity; see also Bagheri et al. (2011). The Q-statistic and the Bliss statistic can also be considered conservative tests, a conclusion reached earlier (Higgins & Thompson, 2002; Higgins et al., 2003; Huedo-Medina et al., 2006). However, heterogeneity in the individual’s risk rates makes it more apparent. The Type 1 error rates of the chi-square methods are close to the nominal value, but Paul and Donner (1989) found $T_2$ to have liberal Type 1 error rates in the setting of a small number of centers and a large number of units. In our study, however, the $T_2$ (and all the ${\chi }^2$ methods) hardly overestimates the Type I error in such setting, but we did not study large numbers of participants within studies. The Type I error rate of $T_2$ is not affected by heterogeneity of individual’s risk rates. Since the allocation of participants to treatment and control groups seems to affect the statistical power, it should be taken into consideration. Additionally, the heterogeneity in individual’s risk rates reduces the power, in particular when there would be no treatment effects present.

Overall, the chi-square tests based on odds ratios and the fixed-effects logistic regression analysis perform best. We recommend the Breslow–Day test adjusted by Tarone for convenience purposes, since it is available in most statistical software packages and the other tests are not. The Breslow–Day test had been recommended before by Jones et al. (1989), but only when the number of studies is large. We also recommend the Breslow–Day test when the number of studies is small. However, when a treatment effect is suspected to be present and we have only a small number of studies ($m\le 10)$, we recommend using the fixed-effects logistic regression over the ${\chi }^2$ tests since it has a slightly higher power than the chi-square test statistics, with an almost nominal Type I error rate. We advise against using the Q-statistic, the Bliss statistic and the $I^2$ statistic because of their poor performance, and due to the recent theoretical worries expressed by Hoaglin (Hoaglin, 2016). The Peto statistic and the random-effects logistic regression are not recommended either due to their poor performance.

This article has its limitations. First, the conclusions presented only hold for dichotomous outcomes. Second, although the simulation method creates data on a unit (subject) level, hence adding more variation, it uses normally distributed latent variables. This assumption may not always be realistic in practice. Finite mixture models for the latent variables can be a potentially useful method. More specifically, the beta-binomial distribution (Skellam, 1948) could be contemplated to perform the interaction test. Another suggestion is to use the methods that estimate between-study variance; see e.g. Brockwell and Gordon (2001) and DerSimonian and Laird (1986) for an overview. The estimated confidence intervals of the between-study variance can be used to test the homogeneity of effect sizes. Finally, Bayesian approaches to test the interaction effects between study and treatment effect are other options for future research.

Acknowledgments

The authors are very grateful to the reviewer whose comments helped improve the results section and thereby expanding the relevance thereof.

Supplementary material

Supplemental material for this article can be accessed at https://doi.org/10.1080/25742558.2018.1478698.

Additional information

Funding

This work was funded by the Netherlands Organization for Scientific Research [grant number: 023.005.087].

Notes on contributors

Osama Almalik

Edwin van den Heuvel is a full-time professor at the Department of Mathematics and Computer Sciences in the Eindhoven University of Technology and he is the head of the Statistics group within the Stochastics section. The Statistics group develops and compares data-analytical methods for analyzing and sampling complex structured correlated datasets. It includes parameter estimation, model fitting, latent variable models, mixed models, missing data, statistical process control, survival and reliability theory, time series analysis, and statistical learning methods. One of the central themes is the analysis of high-dimensional temporal datasets and other large datasets. The group actively explores new research lines in data science. This current research article is part of a PhD project of Osama Almalik (the principal author) on meta-analysis for non-standard or complex situations under supervision of van den Heuvel. It is considered beneficial to medical researchers, as well as statisticians, when carrying out meta-analysis for combining effect sizes from cohort studies and clinical trials.