Maximum likelihood estimation and EM algorithm of Copas-like selection model for publication bias correction

Abstract

Publication bias occurs when the published research results are systematically unrepresentative of the population of studies that have been conducted, and is a potential threat to meaningful meta-analysis. The Copas selection model provides a flexible framework for correcting estimates and offers considerable insight into the publication bias. However, maximizing the observed likelihood under the Copas selection model is challenging because the observed data contain very little information on the latent variable. In this article, we study a Copas-like selection model and propose an expectation-maximization (EM) algorithm for estimation based on the full likelihood. Empirical simulation studies show that the EM algorithm and its associated inferential procedure performs well and avoids the non-convergence problem when maximizing the observed likelihood.

1. Introduction

Meta-analysis provides an increasingly important means for assessing the results from several related studies. By synthesizing information from individual studies, meta-analysis usually has higher statistical power than a single study and can answer research questions that cannot be answered by a single study (Jackson and others, 2011). However, the validity of the results of a meta-analysis is threatened by a well-recognized but challenging issue, publication bias. Publication bias occurs when the decision to publish research results depends on the direction of the results (positive versus negative findings), statistical significance of the results, and other factors (Dickersin, 2005). In the presence of publication bias, the studies that can be identified from the published literature are a biased selection of the research that has been done in that particular area, and may result in misleading conclusions. This impacts the inferential results of both treatment effects and between-study variance (Jackson, 2007).

A great deal of effort has been devoted to developing statistical methods to detect and correct for publication bias by using graphical methods or parametric models of the selection mechanism. The most popular graphical method is the funnel plot, which assesses the asymmetry of a scatter plot of the treatment effects estimated from individual studies against a measure of the study size. The presence of asymmetry in a funnel plot can indicate potential publication bias (Egger and others, 1997; Sterne and others, 2000, 2001). Statistical tests based on funnel plot asymmetry are available in the literature (Begg and Mazumdar, 1994; Egger and others, 1997; Macaskill and others, 2001; Sterne and others, 2001). Duval and Tweedie (2000a) further developed the “trim and fill” method for estimating and adjusting for the number and outcomes of missing studies in meta-analysis. Other methods include modeling the publication bias by including impact factors of the publishing journals (Baker and Jackson, 2006). An alternative approach assumes selection models that impose parametric distributions to characterize the underlying publication mechanism by which effect estimates are selected to be observed (Rothstein and others, 2006). In particular, the selection model by Copas and colleagues (Copas and Li, 1997; Copas, 1999; Copas and Shi, 2000, 2001) provides a flexible framework, in which the probability of selection is allowed to depend on both the effect size and its standard error. Mavridis and others (2013) implemented the Bayesian method for model fitting under the Copas model. As stated, the major advantage of using the Bayesian framework is that the prior knowledge on the probability of publication can be utilized in the specification of the prior distributions. In applications, the Copas selection model is very useful. In a recent empirical evaluation, Carpenter and others (2009) analyzed 157 meta-analyses with binary outcomes and suggested that, the overall interpretation of the Copas selection model was clear for approximately 80% of meta-analyses.

As acknowledged by Copas and Shi (Copas and Shi, 2001), the data often contain very little information about the parameters in the selection model. They suggested conducting a sensitivity analysis. Our simulation studies confirm that maximizing the likelihood of the observed data under the Copas model is challenging and often has non-convergence issues. In the existing R package “metasens”, the copas() function utilizes a grid search algorithm for the parameters in the selection model (Halekoh and others, 2006). Although this computational algorithm avoids the non-convergence issue, the variance estimator reported from the R function “copas” underestimates the true variation by ignoring the fact that the subset of parameters is estimated through the grid search. Motivated by the limitations of existing inferential procedures and the challenges in the computing algorithm, we study a Copas-like selection model and propose an expectation-maximization (EM) algorithm for estimating the maximum likelihood under the Copas-like model. A new feature of our method is that we recast the biased-sampling problem as a missing data problem. This method leads to an alternative likelihood based inference on bias corrected estimation of overall effect size and heterogeneity parameters, and serves as an useful complement to the existing Copas approach.

The rest of this article is organized as follows. In Section 2, we first introduce the Copas selection model and then describe the proposed EM algorithm for estimation. In Section 3, we conduct simulation studies to assess the finite sample performance of the EM algorithm, and compare it with existing estimating methods. We apply the proposed method to a meta-analysis of the relationship between second-hand tobacco smoke and lung cancer in Section 4 and provide a brief discussion in Section 5.

2. Model and estimation

2.1. Model

Consider |$N$| studies to be reviewed. For study |$i$| in a meta-analysis, let |$\hat \theta_i$| denote the estimated effect and |$s_i^2$| denote the estimated variance of |$\hat \theta_i$|⁠. We assumed that |$\hat \theta_i$| and |$s_i^2$| are independent. A random effects model is often assumed (DerSimonian and Laird, 1986), 

$${\hat{\theta }}_i=\theta +\tau u_i+s_i{ϵ}_i,$$

(2.1)

where the random effect |$u_i$| and the error |$\epsilon_i$| are independent and follow a standard normal distribution. Here, the parameter |$\theta$| is the underlying effect size and the parameter |$\tau^2$| describes the between-study heterogeneity. Following Copas and Li (1997), a separate selection model that uses a latent variable |$Z_i$| for the publication process is postulated as  

$$Z_i={\gamma }_0+{\gamma }_1/s_i+{\delta }_i,$$

(2.2)

where |$\delta_i\sim N(0,1)$|⁠, |$\textrm{corr}(\epsilon_i, \delta_i)=\rho$| and |$\textrm{corr}(u_i, \delta_i)=0$|⁠. It is further assumed that the study |$i$| is published if the latent variable |$Z_i>0$|⁠. Clearly, equations (2.1) and (2.2) illustrate that this selection model accounts for the impacts of both the estimated effect size and its precision on the probability of being published. The parameter |$\gamma_0$| controls the overall likelihood of a study being published, and the parameter |$\gamma_1$| characterizes how the chance of publication depends on the sample size. In general, |$\gamma_1$| is positive, so that studies with larger sample sizes are more likely to be published. If |$\rho>0$|⁠, smaller studies with larger estimated effect size will tend to be published. Note that the model in equation (2.1) is slightly different from the equation (1) of Copas and Shi (2000). Here, we treat the within-study variance as known and use the observed within-study variance to approximate the true within-study variance. Under this modified model, the number of parameters is fixed and does not increase with the number of studies in the meta-analysis. Hereafter, we refer this modified Copas model as Copas-like model.

Under the selection model, the probability of the |$i$|th study being published is |$\textrm{P}(Z_i>0) = \textrm{P}(\delta_i > -\gamma_0-\gamma_1/s_i)=\Phi(\gamma_0+\gamma_1/s_i)$|⁠, where |$\Phi(\cdot)$| is the cumulative distribution function of the standard normal distribution. The log-likelihood of observing |$\hat \theta_i$| given |$Z_{i}>0$| and |$s_i$|⁠, |$i=1, \cdots, n$|⁠, can be written as  

$$\begin{array}{ccc}& & \hspace{0.17em}\sum _{i=1}^N\log \text{P}({\hat{\theta }}_i|Z_i>0,s_i)=\sum _{i=1}^N\log \left\{\frac{\text{P}(Z_i>0|{\hat{\theta }}_i,s_i)f({\hat{\theta }}_i)}{\text{P}(Z_i>0|s_i)}\right\}\hfill \\ & & \hspace{0.17em}=\sum _{i=1}^N-\frac{1}{\mathrm{2}}\log ({\tau }^2+s_i^2)-\frac{({\hat{\theta }}_i-\theta {)}^2}{2({\tau }^2+s_i^2)}-\log \mathrm{\Phi }({\gamma }_0+{\gamma }_1/s_i)+\log \mathrm{\Phi }(v_i)+C,\end{array}$$

(2.3)

where |$v_{i}= \left\{ \gamma_0+\gamma_1/s_i+\rho s_i(\hat\theta_i-\theta)/(\tau^2+s_i^2) \right\} /\sqrt{1-\rho^2 s_i^2/(\tau^2+s_i^2)}$| and |$ C$| is a constant. When |$\rho=0$|⁠, i.e., |$\hat \theta_i$| and |$Z_i$| are unrelated, the probability of publication is not impacted by the effect size or its precision, the above |$\log{\Phi(v_{i})}$| reduces to |$\log\Phi\left\{ {\gamma_0+\gamma_1/s_i } \right\}$| and the likelihood function is the likelihood obtained from a standard univariate random effects model for meta-analyses.

As acknowledged by Copas and Shi (2001) and confirmed in our simulation studies, maximizing the likelihood (2.3) is computationally challenging, because the likelihood takes its maximum over a very flat plateau and often has non-convergence results. As an alternative, for given |$(\gamma_0, \gamma_1)$|⁠, the estimators of parameters |$(\theta, \tau, \rho )$| can be stably obtained by maximizing the likelihood. Accordingly, the existing R function “copas” for fitting the Copas selection model uses the grid search algorithm for |$(\gamma_0, \gamma_1)$| (Halekoh and others, 2006). Besides the computational intensity of the grid search algorithm, it is challenging to appropriately estimate the variance of estimators for |$(\theta, \tau, \rho )$|⁠. That is because this approach does not easily account for the uncertainty on the estimated |$(\gamma_0, \gamma_1)$|⁠. The R function “copas” ignores such additional variation due to the grid search algorithm and provides underestimated standard errors for parameters |$(\theta, \tau, \rho )$|⁠, resulting in potentially misleading inference.

2.2. EM algorithm

To account for the limitations of existing inferential procedures and the challenges in the computing algorithm, we exploit the underlying feature of the biased-sampling mechanism and derive an EM algorithm to simultaneously obtain the MLE for |$(\theta, \tau, \rho )$| and |$(\gamma_0, \gamma_1)$|⁠. We denote the observed data as |$O=\{O_i, i=1, \cdots, N\}=\{ \widehat\theta_{i}, s_i, i=1,\cdots, N\}$|⁠. For each observed study |$i$|⁠, |$i=1, \cdots, N$|⁠, we assume that there are additional |$m_i$| (⁠|$m_i\geq 0$| and is an integer) studies that are not observed due to publication bias, where the value of |$m_i$| is unknown. Based on the symmetry assumption (2.1), we let |$O_i^*=\{\widehat{\theta}^*_{ij}=-(\widehat{\theta_i}-2\theta), s_i, I(Z_{ij}<0), j=1, \cdots, m_i\}$| be the corresponding unpublished data. By the biased data generating mechanism for the |$i$|th study, the random integer |$m_i$| follows a geometric distribution |$(1-p_i)^{m}p_i$|⁠, |$m=0, 1, \cdots, $| where |$p_i=P(Z_i>0|\widehat{\theta_i}, s_i)={{\Phi(v_i)}}$| is the probability of the |$i$|th study being published specified in equation (2.2).

Following the principle of the EM algorithm, we treat |$\{O_i, O_i^*, i=1, \cdots, N\}$| as the “complete data” and derive the full likelihood, including the published data and the unpublished data. The corresponding log-likelihood based on the complete data is  

$$\begin{array}{c}\log L^{\ast }(\theta ,{\tau }^2,{\gamma }_0,{\gamma }_1,\rho )=\sum _{i=1}^N\left[\log \left\{P(Z_i>0,\widehat{{\theta }_i}|s_i)\right\}+m_i\log \left\{P(Z_i<0,-(\widehat{{\theta }_i}-2\theta )|s_i)\right\}\right].\end{array}$$

We denote the current parameter value in the EM step by |$\psi^*= (\theta^*,\tau^{*},\gamma_0^{*},\gamma_1^{*},\rho^*)$|⁠. In the E-step, given the observed data and |$\psi^*$|⁠, the conditional expectation of the unobserved |$m_i$| is calculated as  

$$E(m_i|O;{\psi }^{\ast })=\frac{1-P(Z_i>0|O;{\psi }^{\ast })}{P(Z_i>0|O;{\psi }^{\ast })}=\frac{1-\mathrm{\Phi }(v_i|{\psi }^{\ast })}{\mathrm{\Phi }(v_i|{\psi }^{\ast })}.$$

It follows that the expected complete-data log-likelihood function, given the current parameter estimate, is  

$$\begin{array}{cc}\hfill \log L_E(\psi |O;{\psi }^{\ast })& \hspace{0.17em}=\sum _{i=1}^N[\log \left\{P(Z_i>0,\widehat{{\theta }_i}|s_i)\right\}\hfill \\ & \hspace{0.17em}+E(m_i|O,{\psi }^{\ast })\log \left\{P(Z_i<0,-(\widehat{{\theta }_i}-2\theta )|s_i)\right\}].\end{array}$$

(2.4)

In the M-step, we maximize the conditional expectation of the log-likelihood of the complete data from the E-step, and update the estimates for the unknown parameters. With a given convergence criterion, we can solve for |$\psi$| iteratively. We use the estimation procedure by Louis (1982) to obtain the variance estimator.

The purpose of introducing the latent variable |$m_i$| for each observed (published) study is to recast the biased sampling issue as a missing data problem. The distribution property of |$m_i$| is fully determined by the publication mechanism specified in equations (2.1) and (2.2). Note that imputing the unpublished studies by using the EM algorithm is different from using the trim and fill method of Duval and Tweedie (2000b). First, the number of unobserved studies |$m_i$| for study |$i$| is a random variable here, rather than a binary variable (⁠|$0$| or |$1$|⁠) as in the trim and fill method. Specifically, for each published study, the proposed method allows for more than one corresponding unpublished study, which provides more model flexibility. Second, our method considers the uncertainty of the estimated number of unpublished studies, which is ignored by the trim and fill method. Third, the key to the trim and fill method is estimating the total number of unpublished studies by a rank approach, while the validity of the proposed method relies on the Assumption (2.2) on the publication mechanism. The comparison between the Copas selection model and the trim and fill method in terms of model assumptions and empirical performance has been well studied in the literature (Schwarzer and others, 2010; Jin and others, 2015). The justification and comparison of these different assumptions for the mechanism of publication bias is beyond the scope of this paper.

3. Simulation study

We conducted simulation studies to evaluate the finite sample performance of the proposed EM algorithm and compared it to that of two methods: the grid search method implemented by the function copas() in the R package metasens, and the direct maximization of the observed likelihood in (2.3). This function fits the Copas selection model over a grid of (by default) 400 points (Halekoh and others, 2006). To directly maximize the observed likelihood, we used the optim() function in R (R Development Core Team, Version 2.14.1).

As illustrated in Bürkner and Doebler (2014), at least 30 studies are required to have reasonable power to detect the publication bias. Hence, we considered three different sizes of meta-analyses, including a moderate sample size (N = 25), and two large sample sizes (N = 50 or 100). We set the overall effect size at 0.4, and generated the study-specific variance from the square of a normal random variable |$N(0.25,0.5)$|⁠. We generated the random errors |$(\epsilon_i,\delta_i)$| from a bivariate standard normal distribution with correlation |$\rho$|⁠. In the supplementary materials available at Biostatistics online, we considered more simulation settings with different values of |$\rho$|⁠. Since the findings were similar, below we only showed the results when |$\rho=0.4$|⁠. We chose different combinations of parameters to cover a range of settings. By letting |$(\gamma_0,\gamma_1)=(-0.6,0.8)$| and |$(\gamma_0,\gamma_1)=(-1,0.6)$|⁠, we obtained 15% and 30% non-publishing rates, respectively. For each non-publishing rate, we considered two scenarios in terms of the between-study heterogeneity: relatively low between-study heterogeneity, with |$\tau=0.5$|⁠, and relatively high heterogeneity, with |$\tau=1$|⁠.

Tables 1 and 2 summarize the empirical biases, empirical standard errors, average estimated standard errors and coverage probabilities of the |$95\%$| confidence intervals based on 1000 simulations. When the between-study heterogeneity was relatively low (⁠|$\tau=0.5$|⁠), the proposed method produced empirically unbiased estimates with very small biases, while the other two methods had relatively larger empirical biases, especially for the effect size parameter |$\theta$|⁠. Furthermore, the coverage probabilities obtained from the proposed method were close to the nominal level. As expected, the Copas R function underestimated the variation of the estimators, leading to poor coverage probabilities. The non-convergence or non-positive definite covariance matrix of the direct maximization was encountered in about 80% |$\sim$| 90% of the simulations, even with a relatively large sample size (N = 100). When summarizing the simulation results, we have excluded simulations in which there was non-convergence. The decision of a case of non-convergence was determined by either the non-convergence indicator in the output of the R function “optim” or the error message “non-finite finite-difference value”. We further calculated the relative efficiency (RE) of the three methods, defined as the mean square error (MSE) from the Copas R function or from the direct maximization, divided by that of the proposed method. The range of RE was from 1.0 to 1.3 when comparing the Copas R function with the proposed method, and the RE when comparing the direct maximization of likelihood function in equation (2.3) with the proposed method ranged from 1.3 to 1.6. The efficiency gain can be explained by the computational efficiency of the EM algorithm. When the between-study heterogeneity increased to |$1$|⁠, we obtained relatively larger but still small empirical biases by using the proposed EM algorithm.

Table 1.

Summary statistics over 1000 runs when |$\tau=0.5$|⁠. All the entries are multiplied by 100

		Proposed Method				Copas R function				Observed Likelihood
n	PAR	BIAS	SD	SE	CP	BIAS	SD	SE	CP	BIAS	SD	SE	CP	NCR
overall non-publishing rate |$P$| = 15%
25	|$\theta$|	–0.4	12.4	12.3	93.1	1.4	12.7	11.8	91.6	–2.8	15.0	15.2	92.5	92.8
	|$\tau$|	–2.0	10.1	9.9	90.5	–2.7	10.1	–	–	–4.6	14.4	9.5	87.5
50	|$\theta$|	–0.1	8.8	8.8	93.8	1.7	9.3	8.5	91.9	–1.4	10.5	10.6	94.5	88.1
	|$\tau$|	–0.8	7.0	7.0	93.2	–1.3	7.0	–	–	–2.1	10.7	7.0	91.4
100	|$\theta$|	0.1	6.4	6.3	94.0	1.6	7.0	6.1	90.9	–1.1	7.9	7.6	93.6	86.4
	|$\tau$|	–0.4	5.1	5.0	93.1	–0.8	5.1	–	–	–0.6	5.1	5.0	94.0
overall non–publishing rate |$P$| = 30%
25	|$\theta$|	0.1	12.4	13.4	94.2	2.7	12.1	11.5	91.6	–2.6	14.5	15.7	94.5	92.1
	|$\tau$|	–2.2	10.0	9.3	90.3	–2.6	9.4	–	–	–4.3	13.1	8.9	88.4
50	|$\theta$|	–0.3	9.4	9.7	93.4	2.3	9.3	8.3	91.0	–1.8	10.6	10.9	94.5	85.9
	|$\tau$|	–0.7	7.3	6.7	90.6	–1.3	6.7	–	–	–1.6	8.3	6.5	92.0
100	|$\theta$|	0.0	6.9	7.1	93.9	1.9	7.2	6.0	88.0	–1.4	7.8	7.8	95.2	85.8
	|$\tau$|	–0.4	5.3	4.8	91.6	–0.8	4.6	–	–	–0.5	4.6	4.7	94.8

|$n$|⁠, sample size; PAR, parameter; BIAS, empirical bias; SD, empirical standard error; SE, average of asymptotic standard error estimates from the first 2 by 2 submatrix from the information matrix; CP, 95% confidence interval coverage probability; NCR, non convergence rate; –, not available.

Table 2.

Summary statistics over 1000 runs when |$\tau=1$|⁠. All the entries are multiplied by 100

		Proposed Method				Copas R function				Observed Likelihood
n	PAR	BIAS	SD	SE	CP	BIAS	SD	SE	CP	BIAS	SD	SE	CP	NCR
overall non-publishing rate |$P$| = 15%
25	|$\theta$|	0.1	22.1	21.8	92.9	2.9	22.4	20.9	90.5	–5.1	28.6	29.3	92.7	96.2
	|$\tau$|	–4.1	16.9	16.1	88.6	–4.6	16.8	–	–	–6.5	18.2	15.4	86.1
50	|$\theta$|	–0.2	15.5	15.8	94.9	2.6	15.8	15.1	93.3	–3.5	19.5	20.1	94.4	94.6
	|$\tau$|	–2.0	11.7	11.6	92.5	–2.5	11.6	–	–	–3.0	11.4	11.3	92.1
100	|$\theta$|	0.0	11.4	11.2	94.3	2.4	12.1	10.9	91.7	–2.6	14.4	14.2	94.4	94.0
	|$\tau$|	–0.7	8.6	8.3	93.0	–1.2	8.4	–	–	–1.1	8.4	8.2	93.5
overall non–publishing rate |$P$| = 30%
25	|$\theta$|	0.2	22.1	25.7	94.9	4.0	22.1	20.6	92.7	–4.4	27.5	30.1	94.1	96.9
	|$\tau$|	–4.9	16.7	16.3	88.8	–4.4	16.2	–	–	–6.6	18.1	14.8	88.0
50	|$\theta$|	0.1	16.3	18.6	95.4	3.7	16.1	15.0	91.5	–4.1	20.5	20.9	94.9	95.5
	|$\tau$|	–2.4	12.1	12.0	90.5	–2.3	11.2	–	–	–2.8	11.1	11.0	90.7
100	|$\theta$|	0.6	12.5	13.5	94.2	3.9	12.4	10.8	89.1	–2.3	14.2	14.4	94.9	94.6
	|$\tau$|	–0.9	8.9	8.6	90.3	–1.0	8.0	–	–	–0.9	8.0	8.0	94.2

|$n$|⁠, sample size; PAR, parameter; BIAS, empirical bias; SD, empirical standard error; SE, average of asymptotic standard error estimates from the first 2 by 2 submatrix from the information matrix; CP, 95% confidence interval coverage probability; NCR, non convergence rate; –, not available.

One potential limitation of the EM algorithm is that the initial values may play an important role in determining the maxima. To check this issue, we compared the estimates from one initial value with the estimates that represent the maximizer over the 10 maximizations that used different initial values. The results were almost identical. Note that the algorithm converged fairly quickly for this model and the use of multiple starting points was not computationally intensive.

4. Data application

We applied the proposed EM algorithm to a meta-analysis of the relationship between second-hand tobacco smoke and lung cancer. Concerns regarding the effects of second-hand smoke have been central to the debate over the regulation of tobacco products. Since the 1980s, the effect of passive smoking on the development of lung cancer has been extensively studied (Jinot and Bayard, 1992). Hackshaw and others analyzed the results from 37 studies that evaluated the risk of developing lung cancer in women who were lifelong non-smokers, comparing women whose husbands currently smoked to those whose husbands had never smoked (Hackshaw and others, 1997). A random effects meta-analysis model was fitted and the pooled odds ratio was reported as 1.24, with a 95% confidence interval of [1.13, 1.36]. The authors concluded that married women who were lifelong non-smokers but were living with husbands who smoked were 24% more likely to develop lung cancer than those who lived with husbands who did not smoke. The data are available in supplementary materials available at Biostatistics online.

The funnel plot in Figure 1 suggested some evidence of publication bias in that meta-analysis. In addition, Egger’s test (Egger and others, 1997) yielded a P-value of |$0.035$|⁠, and the test based on the trim and fill method (Duval and Tweedie, 2000b) yielded a consistent conclusion, with a P-value of |$0.004$|⁠. To account for this publication bias, we revisited this meta-analysis. The direct maximization of the log-likelihood in equation (2.3) using the optim function in R did not lead to a convergent result. We fit the data with the Copas selection model using the proposed EM algorithm, Copas R function, and compared the results with those obtained by using the random effects model without adjusting for publication bias. The results are summarized in Table 3. The estimated log odds ratios from the three methods were qualitatively similar. The estimated odds ratio from the proposed method was 1.22 (⁠|$=\exp(0.202)$|⁠), with a 95% confidence interval of [1.12, 1.34], which was only slightly less than the estimates from the random effects model that did not account for publication bias. One interesting difference between these two methods was that the between-study heterogeneity |$\tau$| was statistically significant under the proposed method (P = 0.004) but was not significant under the standard random effect model (P = 0.095). Furthermore, the Copas R function did not provide variance estimation for between-study heterogeneity |$\tau$|⁠. Based on the proposed method using the Copas model, the lifelong non-smoking women who lived with husbands who smoked were 22% more likely to develop lung cancer than those who lived with husbands who did not smoke [95% CI: 1.12, 1.34], and there was statistically significant heterogeneity across the studies,which suggests further investigations on sources of heterogeneity.

Table 3.

A meta–analysis of 37 studies on second–hand smoke

	Proposed Method			Copas R Function			Standard method
	Est.	SE	|$P$|–value	Est.	SE	|$P$|–value	Est.	SE	|$P$|–value
|$\theta$|	0.202	0.045	|$<$|0.001	0.180	0.050	|$<$|0.001	0.218	0.049	|$<$|0.001
|$\tau$|	0.135	0.051	0.004	0.135	–	–	0.022	0.018	0.095

Standard method, method using the random effects meta–analysis model without accounting for publication bias;

|$\theta$|⁠, log odds ratio of the risk of developing lung cancer in lifelong non–smoking women, comparing those whose husbands currently smoked with those whose husbands had never smoked; |$\tau$|⁠, between–study heterogeneity; –, not available.

Fig. 1.

Funnel plot for the second-hand smoke data based on 37 studies (Hackshaw and others, 1997).

5. Discussion

Publication bias is a well-recognized challenge in meta-analysis. We have proposed a Copas-like selection model, which is slightly different from the model of Copas and Shi (2000). Here, we treat the within-study variance as known and use the observed within-study variance to approximate the true within-study variance. Under this modified model, the number of parameters is fixed and does not increase with the number of studies in the meta-analysis. Under the Copas-like selection model, we proposed an EM algorithm to obtain the maximum full likelihood estimator, which works well in correcting publication bias and obtaining valid inference. The proposed EM algorithm exploits the underlying feature of the biased-sampling mechanism, including the symmetry assumption, and recasts the biased sampling issue as a missing data problem. As expected and demonstrated in the numerical studies, the proposed EM algorithm yields stable estimates, in contrast to the direct maximization of the observed likelihood. Compared with the grid search method implemented in the Copas R function, the proposed method is more computationally efficient. Note that the identifiability of the parameters in the Copas selection model relies on the distributional assumptions on the study-specific effects and the latent variables. In practice, investigators are encouraged to conduct a sensitivity analysis to evaluate the consistency of results under different specifications of distributions.

Multivariate meta-analysis, which jointly analyzes multiple and possibly correlated outcomes, has recently received increasing attention (Jackson and others, 2011). However, the development of methods that account for publication bias has mainly focused on univariate meta-analysis (Begg and Mazumdar, 1994; Egger and others, 1997; Peters and others, 2006). Although conceptually attractive, a direct extension of the Copas model to multivariate outcomes is expected to result in computational challenges and inconsistent inference: applying univariate methods separately to each outcome may lead to different conclusions on a number of unpublished studies. Proposing robust methods to account for these features in the multivariate meta-analysis framework is beyond the scope of this paper and is worthy of future research.

6. Software

Software in the form of R code, together with a sample data set and documentation is available on request from the corresponding author ([e-mail]).

Acknowledgments

Conflict of Interest: None declared.

Supplementary material

supplementary materials is available online at http://biostatistics.oxfordjournals.org.

Funding

This work was supported in part by grants from the National Institutes of Health (CA016672 and CA193878) and the Agency for Health Care Research and Quality.

Author notes

*

The first two authors are equally contributed.