Simulation study of estimating between-study variance and overall effect in meta-analysis of standardized mean difference

Meta-analysis is a statistical methodology for combining estimated effects from several studies in order to assess their heterogeneity and obtain an overall estimate. In this paper we focus on meta-analysis of standardized mean difference. The data and, often, existing tradition determine the choice of outcome measure. In a comparative study with continuous subject-level data for a treatment arm (T) and a control arm (C), the customary outcome measures are the mean difference (MD) and the standardized mean difference (SMD). The Cochrane Handbook (Higgins and Green, 2011, Part 2, Chapter 9) points out that the choice of MD over SMD depends on whether “outcome measurements in all studies are made on the same scale.” However, fields of application have established preferences: MD in medicine and SMD in social sciences.

If the studies can be assumed to have the same true effect, a meta-analysis uses a fixed-effect (FE) model (common-effect model) to combine the estimates. Otherwise, the studies’ true effects can depart from homogeneity in a variety of ways. Most commonly, a random-effects (RE) model regards those effects as a sample from a distribution and summarizes their heterogeneity via its variance, usually denoted by $\tau^{2}$.The between-studies variance, $\tau^{2}$, has a key role in estimates of the mean of the distribution of random effects; but it is also important as a quantitative indication of heterogeneity (Higgins et al., 2009). In studying estimation for meta-analysis of SMD, we focus first on $\tau^{2}$ and then proceed to the overall effect.

Veroniki et al. (2016) provide a comprehensive overview and recommendations on methods of estimating $\tau^{2}$ and its uncertainty. Their review, however, has two important limitations. First, the authors study only “methods that can be applied for any type of outcome data.” However, the performance of the methods that we study varies widely among effect measures. (Veroniki et al., 2016, Section 6.1) mention this only in passing as a hypothetical possibility. Second, any review on the topic, such as Veroniki et al. (2016), currently can draw on limited empirical information on the comparative performance of the methods. The (short) list of previous simulation studies for SMD is given in Table 1. From this table, it is clear that only three studies (Viechtbauer (2005) Petropoulou and Mavridis (2017) and Langan et al. (2018)) considered and compared several point estimators of $\tau^{2}$.

However, to assess bias of the estimators of heterogeneity variance, Petropoulou and Mavridis (2017) use mean absolute error (the only performance measure reported for SMD in their Table S8), which is not a measure of bias; it is the linear counterpart of mean squared error.

Langan et al. (2018) studied bias and mean squared error of estimators of $\tau^{2}$, as well as coverage of confidence intervals for the overall effect, but they used only one value of SMD=0.5 because they believed that the value of SMD does not matter. We show in our simulations that it does matter. Additionally, both Viechtbauer (2005) and Langan et al. (2018)) consider a very restricted range of $\tau^{2}$ values. There appear to be no studies at all on coverage of $\tau^{2}$. Most of the simulation studies consider only inverse-variance-based estimation of the overall effect, and only two (Langan et al. (2018) and Hamman et al. (2018)) consider t-based confidence intervals for it.

Additionally, all moment-based point estimators of $\tau^{2}$ (i.e., the vast majority of the estimators listed in Table 1), use the inferior $\chi^{2}_{K-1}$ approximation to the distribution of Cochran’s $Q$ statistic, which does not perform well for small to medium sample sizes (Kulinskaya et al., 2011).

To address this gap in information on methods of estimating the heterogeneity variance for SMD, we use simulation to study four methods recommended by Veroniki et al. (2016). These are the well-established methods of DerSimonian and Laird (1986), restricted maximum likelihood, and Mandel and Paule (1970), and the less-familiar method of Jackson (2013). We also include a new method based on improved approximations to the distribution of the $Q$ statistic for SMD (Kulinskaya et al., 2011). We also study coverage of confidence intervals for $\tau^{2}$ achieved by five methods, including the Q-profile method of Viechtbauer (2007), a Q-profile method based on an improved approximation to the distribution of Cochran’s $Q$, and profile-likelihood-based intervals.

For each estimator of $\tau^{2}$, we also study bias of the corresponding inverse-variance-weighted estimator of the overall effect. As our work progressed, it became clear that those inverse-variance-weighted estimators generally had unacceptable bias for SMD. Therefore, we added an estimator (SSW) whose weights depend only on the sample sizes of the Treatment and Control arms. We studied the coverage of the confidence intervals associated with the inverse-variance-weighted estimators, and also the HKSJ interval, the HKSJ interval using the improved estimator of $\tau^{2}$, and the interval centered at SSW and using the improved $\hat{\tau}^{2}$ in estimating its variance.

We assume that each of the $K$ studies in the meta-analysis consists of two arms, Treatment and Control, with sample sizes $n_{iT}$ and $n_{iC}$. The total sample size in Study $i$ is $n_{i}=n_{iT}+n_{iC}$. We denote the ratio of the control sample size to the total by $q_{i}=n_{iC}/n_{i}$. The subject-level data in each arm are assumed to be normally distributed with means $\mu_{iT}$ and $\mu_{iC}$ and variances $\sigma_{iT}^{2}$ and $\sigma_{iC}^{2}$. The sample means are $\bar{x}_{ij}$, and the sample variances are $s^{2}_{ij}$, for $i=1,\ldots,K$ and $j=C$ or $T$.

The standardized mean difference effect measure is

The plug-in estimator $d_{i}=(\bar{x}_{iT}-\bar{x}_{iC})/s_{i}$, known as Cohen’s $d$, is biased in small samples, and we do not consider it further. Instead, we study the unbiased estimator

where $m_{i}=n_{iT}+n_{iC}-2$, and the factor

often approximated by $1-3/(4m-1)$, corrects for bias (Hedges, 1983). This estimator of $\delta$ is sometimes called Hedges’s $g$. The variances in the Treatment and Control arms are usually assumed to be equal. Therefore, $\sigma_{i}$ is estimated by the square root of the pooled sample variance

For the variance of ${g}_{i}$ we use the unbiased estimator

Hedges (1983). The sample SMD ${g}_{i}$ has a scaled non-central $t$-distribution with non-centrality parameter $[n_{i}q_{i}(1-q_{i})]^{1/2}\delta_{i}$ :

Cohen (1988) categorized values of $\delta=0.2,\;0.5,\;0.8$ as small, medium, and large effect sizes. Four studies ( Viechtbauer (2005), Friedrich et al. (2008), Mar´ın-Mart´ınez and Sánchez-Meca (2010); Sánchez-Meca and Mar´ın-Mart´ınez (2010) and Lin (2018)) use these values of SMD in their simulations. However, these definitions of “small,” “medium,” and “large” may not be appropriate outside the behavioral sciences. Ferguson (2009) proposed the values $0.41,\;1.15,\;2.70$ as benchmarks in the social sciences. In an empirical study of 21 ecological meta-analyses by Møller and Jennions (2002), 136 observed values of SMD varied in magnitude from $0.005$ to $3.416$, with mean $0.721$ and $95\%$ confidence interval ($0.622-0.820$).

In meta-analysis, the standard random-effects model assumes that within- and between-study variabilities are accounted for by approximately normal distributions of within- and between-study effects. For a generic measure of effect,

resulting in the marginal distribution $\hat{\theta}_{i}\sim N(\theta,\sigma_{i}^{2}+\tau^{2})$. $\hat{\theta}_{i}$ is the estimate of the effect in Study $i$, and its within-study variance is $\sigma_{i}^{2}$, estimated by $\hat{\sigma}_{i}^{2}$, $i=1,\ldots,K$. $\tau^{2}$ is the between-study variance, which is estimated by $\hat{\tau}^{2}$. The overall effect $\theta$ can be estimated by the weighted mean

where the $\hat{w}_{i}(\hat{\tau}^{2})=(\hat{\sigma}_{i}^{2}+\hat{\tau}^{2})^{-1}$ are inverse-variance weights. The FE estimate $\hat{\theta}$ uses weights $\hat{w}_{i}=\hat{w}_{i}(0)$.

If $w_{i}=1/\hbox{Var}(\hat{\theta}_{i})$, the variance of the weighted mean of the $\hat{\theta}_{i}$ is $1/\sum w_{i}$. Thus, many authors estimate the variance of $\hat{\theta}_{\mathit{RE}}$ by $\left[\sum_{i=1}^{K}\hat{w}_{i}(\hat{\tau}^{2})\right]^{-1}$. In practice, however, this estimate may not be satisfactory Sidik and Jonkman (2006); Li et al. (1994); Rukhin (2009).

Because the $\hat{w}_{i}(\tau^{2})$ in Equation (3.2) involve the $\hat{\sigma}_{i}^{2}$, $K-1$ is an adequate approximation for the expected value of Cochran’s $Q$ statistic only for very large sample sizes. As an alternative one can use one of the improved approximations to the expected value of Cochran’s $Q$. Corrected Mandel-Paule methods for estimating $\tau^{2}$ equate Cochran’s $Q$ statistic with the weights $\hat{w}_{i}(\tau^{2})$ to the first moment of an improved approximate null distribution.

More-realistic approximations to the distribution of $Q$ are available for several effect measures. In these approximations the estimates $\hat{\sigma}_{i}^{2}$ are not treated as equal to the $\sigma_{i}^{2}$. For SMD, Kulinskaya et al. (2011) derived $O(1/n)$ corrections to moments of $Q$ and suggested using the chi-squared distribution with degrees of freedom equal to the estimate of the corrected first moment to approximate the distribution of $Q$. Kulinskaya et al. (2011) give expressions from which it can be calculated, along with a computer program in R.

We propose a new method of estimating $\tau^{2}$ based on this improved approximation. Let $E_{KDB}({Q})$ denote the corrected expected value of $Q$. Then one obtains the KBD estimate of $\tau^{2}$ by iteratively solving

We denote the resulting estimator of $\tau^{2}$ by $\hat{\tau}_{KDB}^{2}$.

We also propose a new KDB confidence interval for the between-study variance. This interval for $\tau^{2}$ combines the Q-profile approach and the improved approximation by Kulinskaya et al. (2011) (i.e., the chi-squared distribution with fractional degrees of freedom based on the corrected first moment of $Q$).

This corrected Q-profile confidence interval can be estimated from the lower and upper quantiles of $F_{Q}$, the cumulative distribution function for the improved approximation to the distribution of $Q$:

The upper and lower confidence limits for $\tau^{2}$ can be calculated iteratively.

In an attempt to avoid the bias in the inverse-variance-weighted estimators, we included a point estimator whose weights depend only on the studies’ sample sizes. For this estimator (SSW), $w_{i}=\tilde{n}_{i}=n_{iT}n_{iC}/(n_{iT}+n_{iC})$; $\tilde{n}_{i}$ is the effective sample size in Study $i$. These weights would coincide with the inverse-variance weights when $\delta=0$. These effective-sample-size-based weights were suggested in (Hedges and Olkin, 1985, p.110).

The interval estimator corresponding to SSW (SSW KDB) uses the SSW point estimator as its center, and its half-width equals the estimated standard deviation of SSW under the random-effects model times the critical value from the $t$ distribution on $K-1$ degrees of freedom. The estimator of the variance of SSW is

in which $v_{i}^{2}$ comes from Equation (2.1) and $\hat{\tau}^{2}=\hat{\tau}_{\mathit{KDB}}^{2}$.

As mentioned in Section 1, other studies have used simulation to examine estimators of $\tau^{2}$ or of the overall effect for SMD, but gaps in evidence remain.

Our simulation study for SMD uses $0\leq\delta\leq 2$ and $0\leq\tau^{2}\leq 2.5$ as realistic for a range of applications.

Our simulation study assesses the performance of five methods for point estimation of between-study variance $\tau^{2}$ (DL, REML, J, MP, and KDB) and five methods of interval estimation of $\tau^{2}$ (Q-profile-based methods corresponding to DerSimonian-Laird and KDB, the generalized Q-profile intervals of Biggerstaff and Jackson (2008) and Jackson (2013), and the profile-likelihood confidence interval based on REML).

We also assess the performance of the point and interval estimators of $\delta$ in the random-effects model for SMD.

We vary five parameters: the overall true SMD ($\delta$), the between-studies variance ($\tau^{2}$), the number of studies ($K$), the studies’ total sample size ($n$ and $\bar{n}$), and the proportion of observations in the Control arm ($q$). The combinations of parameters are listed in Table 2.

All simulations use the same numbers of studies $K=5,\;10,\;30$ and, for each combination of parameters, the same vector of total sample sizes $n=(n_{1},\ldots,n_{K})$ and the same proportions of observations in the Control arm $q_{i}=n_{iC}/n_{i}=.5,\;.75$ for all $i$. The values of $q$ reflect two situations for the two arms of each study: approximately equal (1:1) and quite unbalanced (1:3). The sample sizes in the Treatment and Control arms are $n_{iT}=\lceil{(1-q_{i})n_{i}}\rceil$ and $n_{iC}=n_{i}-n_{iT}$, $i=1,\ldots,K$.

We study equal and unequal study sizes. For equal study sizes $n_{i}$ is as small as 20, and for unequal study sizes $n_{i}$ is as small as 12, in order to examine how the methods perform for the extremely small sample sizes that arise in some areas of application. In choosing unequal study sizes, we follow a suggestion of Sánchez-Meca and Mar´ın-Mart´ınez (2000), who selected study sizes having skewness of 1.464, which they considered typical in behavioral and health sciences. Table 2 gives the details.

The patterns of sample sizes are illustrative; they do not attempt to represent all patterns seen in practice. By using the same patterns of sample sizes for each combination of the other parameters, we avoid the additional variability in the results that would arise from choosing sample sizes at random (e.g., uniformly between 20 and 200).

We use a total of $10,000$ repetitions for each combination of parameters. Thus, the simulation standard error for estimated coverage of $\tau^{2}$ or $\delta$ at the $95\%$ confidence level is roughly $\sqrt{0.95\times 0.05/10,000}=0.00218$.

We generate the true effect sizes $\delta_{i}$ from a normal distribution: $\delta_{i}\sim N(\delta,\tau^{2})$. We generate the values of Hedges’s estimator ${g}_{i}$ directly from the appropriately scaled non-central $t$-distribution, given by Equation (2.2), and obtain their estimated within-study variances from Equation (2.1).

The simulations were programmed in R version 3.3.2 using the University of East Anglia 140-computer-node High Performance Computing (HPC) Cluster, providing a total of 2560 CPU cores, including parallel processing and large memory resources. For each configuration, we divided the 10,000 replications into 10 parallel sets of 1000 replications.

The structure of the simulations invites an analysis of the results along the lines of a designed experiment, in which the variables are $\tau^{2}$, $n$, $K$, $q$, and $\delta$. Most of the variables are crossed, but two have additional structure. Within the two levels of $n$, equal and unequal, the values are nested: $n=20,\;40,\;100,\;250$ and $\bar{n}=30,\;60,\;100,\;160$. We approach the analysis of the data from the simulations qualitatively, to identify the variables that substantially affect (or do not affect) the performance of the estimators as a whole and the variables that reveal important differences in performance. We might hope to describe the estimators’ performance one variable at a time, but such “main effects” often do not provide an adequate summary: important differences are related to certain combinations of two or more variables.

We use this approach to examine bias and coverage in estimation of $\tau^{2}$ and bias and coverage in estimation of $\delta$. Our summaries of results are based on examination of the figures in the corresponding Appendices. Section 8 gives brief summaries, and Appendix A contains more detail.

Our full simulation results, comprising $130$ figures, each presenting $12$ combinations of the 4 values of $n$ or $\bar{n}$ and the 3 values of $K$, are provided in Appendices A and B. A summary is given below.

The results of our simulations for SMD give a rather disappointing picture. In brief:
Because the study-level effects and their variances are related (cf. Equation (2.1) for SMD), the performance of all statistical methods depends on the effect measures, estimates of overall effects are biased, and coverage of confidence intervals is too low, especially for small sample sizes.

The conventional wisdom is that these deficiencies do not matter, as meta-analysis usually deals with studies that are “large,” so all these little problems are automatically resolved. Unfortunately, this is not true, even in medical meta-analyses; in Issue 4 of the Cochrane Database 2004, the maximum study size was $63$ or less in $25\%$ of meta-analyses with $K\geq 3$ that used SMD as an effect measure, and less than $110$ in $50\%$ of them (our own analysis). We have not surveyed typical study sizes in psychology, but Sánchez-Meca and Mar´ın-Mart´ınez (2010), promoting MA in psychological research, use an example with 24 studies in which the smallest study size is $12$ and the largest is $121$. In ecology, typical sample sizes are between 4 and 25 (Hamman et al., 2018). An effect-measure-specific estimator of $\tau^{2}$, such as KDB for SMD, can reduce inherent biases.

Arguably, the main purpose of a meta-analysis is to provide point and interval estimates of an overall effect.

Usually, after estimating the between-study variance $\tau^{2}$, inverse-variance weights are used in estimating the overall effect (and, often, its variance). This approach relies on the theoretical result that, for known variances, and given unbiased estimates $\hat{\theta}_{i}$, it yields a Uniformly Minimum-Variance Unbiased Estimate (UMVUE) of $\theta$. In practice, however, the true within-study variances are unknown, and use of the estimated variances makes the inverse-variance-weighted estimate of the overall effect biased.

Consumers routinely expect point estimates to have no (or small) bias and CIs to have (close to) nominal coverage. Thus, the IV-weighted approach is unsatisfactory because, in general, it cannot produce an unbiased estimate of an overall effect.

A pragmatic approach to unbiased estimation of $\theta$ uses weights that do not involve estimated variances of study-level estimates, for example, weights proportional to the study sizes $n_{i}$. Hunter and Schmidt (1990) and Shuster (2010), among others, have proposed such weights, and Mar´ın-Mart´ınez and Sánchez-Meca (2010) and Hamman et al. (2018) have studied the method’s performance by simulation for SMD. We prefer to use weights proportional to an effective sample size, $\tilde{n}_{i}=n_{iT}n_{iC}/n_{i}$; these are the optimal inverse-variance weights for SMD when $\delta=0$ and $\tau^{2}=0$. Thus, the overall effect is estimated by $\hat{\theta}_{\mathit{SSW}}=\sum\tilde{n}_{i}\hat{\theta}_{i}/\sum\tilde{n}_{i}$, and its variance is estimated by Equation (5.1). Hamman et al. (2018) use weights proposed by Hedges (1982), which differ slightly for very small sample sizes.

A good estimator of $\tau^{2}$, such as MP or KDB, can be used as $\hat{\tau}^{2}$. Further, confidence intervals for $\delta$ centered at $\hat{\delta}_{\mathit{SSW}}$ with $\hat{\tau}_{\mathit{KDB}}^{2}$ in Equation (5.1) can be used.

This approach based on SSW requires further study. For example, in the confidence intervals we have used critical values from the $t$-distribution on $K-1$ degrees of freedom, but we have not yet examined the actual sampling distribution of SSW. The raw material for such an examination is readily available: For each situation in our simulations, each of the $10,000$ replications yields an observation on the sampling distribution of SSW.

The work by E. Kulinskaya was supported by the Economic and Social Research Council [grant number ES/L011859/1].

• •

  Appendix A. SMD: Plots for bias and coverage of $\tau^{2}$.

• •

  Appendix B. SMD: Plots for bias and coverage of standardized mean difference $\delta$.

For coverage of $\hat{\tau}^{2}$, each figure corresponds to a value of $\delta(=0,0.5,1,1.5,2,2.5)$, a value of $q(=.5,.75)$, and a set of values of $n$ (= 20, 40, 100, 250 or 30, 50, 60, 70) or $\bar{n}(=30,60,100,160)$.
Each figure contains a panel (with $\tau^{2}$ on the horizontal axis) for each combination of n (or $\bar{n}$) and $K(=5,10,30)$.
The interval estimators of $\tau^{2}$ are

• •

  QP (Q-profile confidence interval)

• •

  BJ (Biggerstaff and Jackson interval )

• •

  PL (Profile likelihood interval)

• •

  KDB ( KDB - improved Q-profile method based on Kulinskaya, Dollinger and Bjørkestøl (2011))

• •

  J (Jackson’s interval)

For bias of $\hat{\delta}$, each figure corresponds to a value of $\delta(=0,0.5,1,1.5,2,2.5)$, a value of $q(=.5,.75)$, and a set of values of $n$ (= 20, 40, 100, 250 or 30, 50, 60, 70) or $\bar{n}(=30,60,100,160)$.
Figures for mean squared error (expressed as the ratio of the MSE of SSW to the MSEs of the inverse-variance-weighted estimators that use the MP or KDB estimator of $\tau^{2}$) use the above values of $\delta$ and q but only n = 20, 40, 100, 250.
Each figure contains a panel (with $\tau^{2}$ on the horizontal axis) for each combination of n (or $\bar{n}$) and $K(=5,10,30)$.
The point estimators of $\delta$ are

• •

  DL (DerSimonian-Laird)

• •

  REML (restricted maximum likelihood)

• •

  MP (Mandel-Paule)

• •

  KDB (improved moment estimator based on Kulinskaya, Dollinger and Bjørkestøl (2011))

• •

  J (Jackson)

• •

  SSW (sample-size weighted)

For coverage of $\hat{\delta}$, each figure corresponds to a value of $\delta$ (= 0, 0.5, 1, 1.5, 2 , 2.5), a value of q (= .5, .75), and a set of values of n (= 20, 40, 100, 250 or 30, 50, 60, 70) or $\bar{n}$ (= 30, 60, 100, 160).
Each figure contains a panel (with $\tau^{2}$ on the horizontal axis) for each combination of n (or $\bar{n}$) and $K(=5,10,30)$.
The interval estimators of $\delta$ are the companions to the inverse-variance-weighted point estimators

• •

  DL (DerSimonian-Laird)

• •

  REML (restricted maximum likelihood)

• •

  MP (Mandel-Paule)

• •

  KDB (improved moment estimator based on Kulinskaya, Dollinger and Bjørkestøl (2011))

• •

  J (Jackson)

and

• •

  HKSJ (Hartung-Knapp-Sidik-Jonkman)

• •

  HKSJ KDB (HKSJ with KDB estimator of $\tau^{2}$)

• •

  SSW (SSW as center and half-width equal to critical value from $t_{K-1}$

times estimated standard deviation of SSW with $\hat{\tau}^{2}$ = $\hat{\tau}^{2}_{KDB}$.