Testing for publication bias in meta-analysis under Copas selection model

In this subsection we restate the Copas selection model proposed by copas1997inference, and describe a slight modification. Let $Y_{i}$ denote the reported effect size in the $i$-th study. We assume $Y_{i}$ follows a normal distribution with mean $\mu_{i}$ and variance $\sigma_{i}^{2}$, where $\sigma_{i}^{2}$ is the true study variance which is often not observed, and we observe only an estimator $s_{i}^{2}$ of the true variance. Suppose the underlying effect size of the study is $\mu_{i}$ and it is assumed to follow a normal distribution with mean $\mu$, the population averaged effect size, and variance $\tau^{2}$. The parameter $\tau^{2}$ describes the between study heterogeneity, while $\sigma_{i}^{2}$ denotes the within-study sampling variance.

We model the observed outcome of study $i$ using the usual random effects model as

Model (2.1) can also be written as

According to the Copas selection model, the study $i$ is published if and only if a postulated latent variable $Z_{i}>0$, where this $Z_{i}$ can be written as

with $\gamma_{0}$ and $\gamma_{1}$ being unknown parameters which control the marginal probability of observing a study and $\delta_{i}$ is a random variable following a standard normal distribution. The parameter $\gamma_{0}$ controls the overall proportion of selection and $\gamma_{1}$ controls how the selection depends on the size of the study. The sign of $\gamma_{1}$ is expected to be positive, since larger studies are more likely to be accepted for publication. Publication bias is modeled by assuming a correlation $\rho$ between $\epsilon_{i}$ and $\delta_{i}$, such that

The correlation $\rho$ controls how the probability of observing a study is influenced by the effect size of a study. From equation (2.2), we can calculate the marginal probability of the $i$-th study with a standard error $s_{i}$ being published as

where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution.

To make inference on the population averaged effect size $\mu$, we can construct the log-likelihood of the observed studies,

where

In the above likelihood, the unknown parameters $\sigma^{2}_{i}$ can be approximated using the observed variables to reduce the number of unknown parameters. In copas2000meta, a function of the reported sampling variance $s_{i}^{2}$ is used to approximate $\sigma^{2}_{i}$, based on the derivation from the conditional probability. However, we propose to directly use the reported sampling variance $s_{i}^{2}$ to approximate the unknown within-study variance $\sigma_{i}^{2}$. In the publication process, the study is first conducted and then selected. Regardless of whether the study is suppressed or observed, the value of $s_{i}^{2}$ is obtained before the publication process, which should provide a fairly good estimation of the true variance $\sigma_{i}^{2}$. By replacing $\sigma_{i}^{2}$ with $s_{i}^{2}$, the unknown parameters in the likelihood function are reduced to $\mu,\tau,\rho,\gamma_{0}$ and $\gamma_{1}$.

When the correlation parameter $\rho=0$, it indicates that $Y_{i}$ and $Z_{i}$ are uncorrelated. In this case, the quantity $\Phi(v_{i})$ in the likelihood function (2.4) reduces to the marginal probability $\Phi(\gamma_{0}+\gamma_{1}/s_{i})$. The likelihood then reduces to the likelihood of the univariate random-effects model. Therefore, the univariate random-effects meta-analysis leads to correct inference for the overall effect $\mu$ and between-study variance $\tau^{2}$ when the correlation $\rho$ is zero. On the other hand, when $\rho\neq 0$, the standard meta-analysis leads to biased estimation of the parameters $\mu$ and $\tau^{2}$. Thus, under the model specified by equations (2.1) and (2.2), testing for publication bias is equivalent to testing for $H_{0}:\rho=0$.

However, the standard results for asymptotic tests may not apply because of the following challenges:
(N1) Parameters $\gamma_{0}$ and $\gamma_{1}$ are not identifiable under the null hypothesis;
(N2) The Fisher information matrix is singular under the null hypothesis.
The asymptotic distribution of the likelihood ratio test under the above conditions is often not a $\chi^{2}$ distribution. Hypothesis testing procedures are non-standard, and have to be considered on a case-by-case basis (davies1977hypothesis; davies1987hypothesis; rotnitzky2000likelihood). For example, the test of homogeneity in mixture models or mixture regression models, genetic admixture models, case-control studies with contaminated controls and testing for partial differential gene expression in microarray studies (qin2011hypothesis; ning2014class; HongJASA2017). Naive tests ignoring the above non-regularity conditions can produce misleading results. For example, the Wald test based on $\hat{\rho}/\hat{se}(\hat{\rho})$ is invalid because the negative Hessian matrix of the log-likelihood is singular under the null hypothesis.

Motivated by the method in Rathouz1999Biometrika, when $\gamma_{0}$ and $\gamma_{1}$ are fixed at some known numbers, the submatrix of the Fisher information matrix is non-singular, and consequently the non-regular conditions (N1) and (N2) are avoided. To account for these facts, we propose a score test constructed as:

where $(\hat{\mu},\hat{\tau}^{2})$ is the constrained MLE of $(\mu,\tau^{2})$ under $\rho=0$ for a given pair of values for $(\gamma_{0},\gamma_{1})$,

and $\hat{I}_{\rho|\mu,\tau}(\gamma_{0},\gamma_{1},\mu,\tau,\rho)$ is the estimated partial information matrix where its explicit form is defined in the Supplementary Material.

This test statistic is to calculate the score test statistic $Z(\gamma_{0},\gamma_{1},\hat{\mu},\hat{\tau})$ for a given pair of values for $\gamma_{0}$ and $\gamma_{1}$, and then construct $T$ by taking the supreme value of $Z(\gamma_{0},\gamma_{1},\hat{\mu},\hat{\tau})$ over a grid of $\gamma_{0}$ and $\gamma_{1}$. The range of $\gamma_{0}$ and $\gamma_{1}$ could be decided similarly to the sensitivity analysis in the Copas model. Usually, the range for $\gamma_{0}$ could be $[-2,2]$, and the range for $\gamma_{1}$ depends on the range of $s_{i}$ to make sure that the probability of the largest study being published is greater than a certain probability threshold (e.g. 0.9) and the probability of the smallest study being published is less than a certain probability threshold (e.g. 0.1). Additional discussion regarding how to choose the grid can be found in Appendix C of the Supplementary Materials. Deriving the asymptotic distribution of $T$ can be challenging and it may not be of a simple form for practical use. In such a case, a bootstrap procedure can be used if the asymptotic distribution is well behaved.

In this subsection, we provide the asymptotic distribution of the proposed test, followed by a parametric bootstrap procedure for practical application.

Theorem 1 provides the limiting distribution of the proposed test statistic $T$, and also provides the theoretical foundation for the following parametric bootstrap procedure which offers a simpler way to obtain a p-value for the test.

1. 1.

   For $n$ observed studies of a meta-analysis, calculate the test statistic in equation (2.5), denoted by $T$.

2. 2.

   Estimate $(\hat{\mu}_{0},\hat{\tau}_{0})$, which is the MLE of $(\mu,\tau)$ under the null hypothesis of $\rho=0$.

3. 3.

   Generate $B_{boot}$ independent bootstrap samples (e.g. $B_{boot}$=200), each consisting of $n$ observations. Specifically, for $i\in\{1,\dots n\}$ we generate $\epsilon_{i}\sim N(0,1)$ independently, and $s_{i,Boot}$ are drawn from the within-study standard errors of the original data with replacement. Then the bootstrap sample is constructed by sampling $\mu_{i}\sim N(\hat{\mu}_{0},\hat{\tau}_{0}^{2})$, and $y_{i,Boot}=\mu_{i}+s_{i,Boot}\epsilon_{i}.$

4. 4.

   Evaluate the score test statistic for each bootstrap sample. Denote $T^{*}=(T_{1}^{*},\dots,T_{B_{boot}}^{*})$.

5. 5.

   The p-value is computed by $\textup{p-value}=\frac{\#\{T^{*}>T\}}{B_{boot}}$.

In meta-analysis, the number of studies is often small or moderate. Therefore the computation time for the above parametric bootstrap is reasonable.

We first studied the influence of number of grid points $p$ on the performance of the test in terms of both type I error and power. Figure 1 presented the empirical rejection rates (i.e., type I error and power) under various numbers of grid points and sample size settings. In order to evaluate the loss of power due to the lack of knowledge on the parameters $(\gamma_{0},\gamma_{1})$ that characterize the publication process, we included a test statistic calculated from equation (2.6) with $(\gamma_{0},\gamma_{1})$ fixed at the true value. This test, although not practical in real application due to the unknown truth of $(\gamma_{0},\gamma_{1})$, can serve as a gold standard that informs the best possible statistical power under the Copas model.

Figure 1 revealed the following interesting findings. First, under the null hypothesis of no publication bias, i.e. $\rho=0$, for all scenarios with different numbers of grid points, the empirical type I errors were close to the nominal level. This suggested that the parametric bootstrap described in Section 2.3 performed well. Secondly, by selecting more grid points, the statistical power of the proposed test increased. For example, when the sample size $n=40$ and $\rho=0.6$, the power increased from 85.1% with $p=1$ to 93.6% with $p=2500$. For most cases, the increase of power was rather limited when $p$ changed from $9$ to $2500$. Thirdly, by comparing the empirical power of the proposed test with the power of the test with $(\gamma_{0},\gamma_{1})$ fixed at the true value (i.e., the last column), we found that the proposed test only suffered a minimal loss of power due to the unknown $(\gamma_{0},\gamma_{1})$. These findings suggested that in practice, the proposed test can be performed by taking several random points from the grid (e.g. $n_{g}=10$), and it can still ensure a close statistical power to the ideal test with known $\gamma_{0}$ and $\gamma_{1}$.

We compared the proposed test with the Trim and Fill test, Egger’s test, and the test proposed in copas2000meta. From the results presented in Table 1, it was reasonable to consider the proposed test with only a few grid points. Thus we considered the proposed test with 9 points selected from the grid.

Figure 2 showed the comparison of type I error and power. For all scenarios, the type I errors were well controlled for the proposed test, but were slightly inflated for Egger’s test. For the power comparison, the proposed test had the highest power among all scenarios. The Trim and Fill method outperformed the Egger’s test and Copas naive test when sample size was small, while the Copas naive test outperformed the Trim and Fill method and Egger’s test for larger sample sizes.

In the Copas model, Equation (2.1) is the random-effects assumption commonly used in the area of meta-analysis, which captures both the normality of estimates from each study, and the heterogeneity across studies. Equation (2.2) is a working model describing the publication process, which captures several main features. First, there is stronger suppression on smaller studies compared to larger studies. They therefore set the latent variable to be positively correlated to the precision of the estimate, i.e., $1/s_{i}$, which can be considered as a proxy of the size of the study. Second, if there is publication bias, studies with larger effect sizes ($Y_{i}$) has higher chance to be published given the precision ($s_{i}$). In the Copas model, this is captured by the correlation $\rho$, where a larger $\rho$ leads to stronger correlation between $Y_{i}$ and $Z_{i}$.

It is interesting to investigate the sensitivity of the proposed method, which is based on the Copas model, when the Copas model is not the correct model. To come up with some alternatives of the Copas selection model, we retain the random-effects model in Equation (2.1), and alter the latent variable $Z_{i}$. We also maintain the two features where stronger suppression on smaller studies compared to larger studies, and studies with larger effect sizes ($Y_{i}$) have higher chances to be published given the precision ($s_{i}$) in the existence of publication bias. We considered the following two alternative models.

The first model modifies the functional relationship between $Z_{i}$ and the term $1/s_{i}$, by changing equation (2) to

so that $Z_{i}$ is a linear function of $1/s_{i}^{2}$, instead of $1/s_{i}$. This is possible since the variance $s_{i}^{2}$ is roughly proportional to sample size of study $i$, and could be a better proxy of the size of the study.

Secondly, we change the way of characterizing publication bias. Instead of assuming a correlation $\rho$, we let $Z_{i}$ to be directly related to $y_{i}/s_{i}$ in the existence of publication bias, by assuming

where $c$ is some positive scaling constant. This model links $Z_{i}$ with the $z$-score $y_{i}/s_{i}$, as the selection process could be depending on the $z$-score directly. Under Equation (3.2), it still holds that $\rho=0$ means no publication bias and studies with larger effect sizes ($Y_{i}$) have higher chances to be published given the precision ($s_{i}$) in the existence of publication bias. In our simulation, we let $\rho$ vary from $0$ to $0.8$ to match the previous setting (in fact, in this alternative model $\rho$ does not need to be within $0$ and $1$), and set $c=0.5$ to obtain a proper range for the power curve.

To investigate the power and type I error of the proposed test under misspecified models, we simulated data from (3.1) and (3.2), and compared the performance of our method with existing methods. The results are presented in Figure 3. Under these two settings where the Copas model is the misspecified model, the power of the proposed test with $9$ points from the grid had similar power as the Trim and Fill method when sample size was small, and outperformed other methods with the increase of sample size. The type I error of the proposed method was still well-controlled around the nominal levels.

In summary, our simulation studies demonstrated that the proposed score test under Copas selection model behaves well with controlled Type I error and competitive statistical power relative to the existing tests. Furthermore, the implementation of the proposed test requires a small number of grid points and is computationally straightforward using parametric bootstrap.

Antidepressants are among the world’s most widely prescribed drugs, where many meta-analyses were conducted to synthesize evidence from existing trials to obtain a more reliable and generalizable results regarding treatment efficacy and safety (barbui2011efficacy; cipriani2018comparative). By comparing trials of antidepressants that were published in the literature and unpublished trials obtained from the Food and Drug Administration (FDA), Turner et al. turner2008selective identified strong evidence for selective bias due to publication towards results favoring active interventions.

Although collecting information for unpublished trials is normally not feasible for meta-analyses, this unique dataset from turner2008selective provide an opportunity to validate our method in an ideal situation, where the unpublished studies were actually observed. Among the 74 FDA-registered studies, 31% were not published. Figure 4 shows the contour enhanced funnel plots for all the published trials and all the registered trials.

By comparing the two funnel plots, we observed clear evidence that several non-significant results were suppressed for publication. Applied on all the published studies, our test yielded a p-value less than $0.01$, showing strong evidence of publication bias. Egger’s test also captured the evidence for bias with a p-value less than $0.01$. The Trim and Fill method, however, has a p-value of $0.25$, possibly driven by the left most trial with relatively large precision and small effect, which suggests its lack of power in the existence of outlying studies.

A more interesting investigation is to apply the above methods on all the FDA-registered trials. By including all the unpublished studies, our test had a p-value of $0.24$ which suggested insufficient evidence for publication bias. Egger’s test, however, had a p-value of $0.07$, still showing marginal evidence for SSE which was likely not caused by publication bias. Trim and Fill test maintained a non-significant p-value of $0.50$, as the symmetry of the funnel plot is improved by adding all the unpublished results.

This case study demonstrated the capability our method for testing publication bias, which on the other than also showed the different focus of our test compared to existing methods. Based on the Copas model, our test is to detect the evidence for publication bias, where majority of the current methods, including Egger’s test and the Trim and Fill method, is for testing SSE. Combing our test with existing methods, we are able to provide a more comprehensive understanding of the potential sources of bias.

When patients have surgery on their abdomen, they are at risk to develop ileus, which is the inability of the intestine (bowel) to contract normally and move waste out of the body. Chewing gum could be one possible way to prevent ileus, since it tricks the body into thinking it is eating, and may help the digestive system to start working again. To investigate the effect of gum-chewing, Short et al. short2015chewing did a systematic review, which included relevant trials focussed on people having bowel surgery, caesarean section, or other abdominal surgery types. Using time to first flatus (TFF) as the primary outcome, slight improvement was identified comparing the gum-chewing groups to the control groups with an overall decrease of TFF ranging from 7.92 to 12.64 hours across different subgroups of surgery types.

To invesigate the potential risk for publication bias, we applied our method to the subgroup of trials for abdominal surgery types other than colorectal surgery and Caesarean, which was referred to the“Other Surgery” subgroup in short2015chewing. This subgroup included in total 43 studies and 4224 participants. Figure 5 presented the funnel plot of the mean different of TFF comparing the control group to the gum-chewing group. From Figure 5, we observed a possible trend of selection by sign, where more studies with positive mean differences were observed. Many larger studies had mean difference near $0$, but only one study among the 43 studies had negative mean difference.

Applying our method with a $9$-point grid, we obtained a p-value of $0.03$, suggesting significance evidence for publication bias. The Trim and Fill method had a p-value of $0.50$ and the Egger’s test had a p-value of $0.12$. The non-significant p-values from the Trim-and-Fill method and the Egger’s test were likely due to the lack of funnel shape of the studies, where larger study were not observed to be more centered compared to the smaller studies, and the fact that the funnel plot was relatively symmetric.

In this study, our test was able to capture the evidence for potential publication bias, which was consistent with our observation from the funnel plot. We then conducted the Copas sensitivity analysis and obtained the adjusted overall mean difference to be $9.92(7.82,12.02)$, which was smaller compared to the overall mean difference $10.57(8.47,12.68)$ obtained by short2015chewing. Although both results suggested chewing gum has significant effect in terms of decreasing the TFF, our testing approach revealed the risk for publication bias, where sensitivity analyses are recommended in order to explore the possible overall mean difference after adjusting for publication bias. We believe that our method can be used as an complement of the existing methods for evaluating the risk of publication bias, which is crucial in assessing the reliability of conclusions from meta-analysis.