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

TL;DR

This study evaluates various estimators for between-study variance and overall effect in meta-analysis of standardized mean differences through extensive simulations, highlighting their bias, coverage, and performance.

Contribution

It introduces and compares new and existing estimators for between-study variance and overall effect in meta-analysis, providing comprehensive simulation results.

Findings

KDB estimator shows improved accuracy in variance estimation.

Interval estimators vary in coverage and precision depending on the method.

Simulation results guide best practices for meta-analysis of SMD.

Abstract

Methods for random-effects meta-analysis require an estimate of the between-study variance, $\tau^2$. The performance of estimators of $\tau^2$ (measured by bias and coverage) affects their usefulness in assessing heterogeneity of study-level effects, and also the performance of related estimators of the overall effect. For the effect measure standardized mean difference (SMD), we provide the results from extensive simulations on five point estimators of $\tau^2$ (the popular methods of DerSimonian-Laird, restricted maximum likelihood, and Mandel and Paule (MP); the less-familiar method of Jackson; the new method (KDB) based on the improved approximation to the distribution of the Q statistic by Kulinskaya, Dollinger and Bj{\o}rkest{\o}l (2011) ), five interval estimators for $\tau^2$ (profile likelihood, Q-profile, Biggerstaff and Jackson, Jackson, and the new KDB method), six point estimators of the overall effect (the five related to the point estimators of $\tau^2$ and an estimator whose weights use only study-level sample sizes), and eight interval estimators for the overall effect (five based on the point estimators for $\tau^2$; the Hartung-Knapp-Sidik-Jonkman (HKSJ) interval; a modification of HKSJ; and an interval based on the sample-size-weighted estimator).

Tables

Table 1. Table 1: Simulation studies on meta-analysis of SMD. Estimators of τ 2 superscript 𝜏 2 \tau^{2} : DL - DerSimonian and Laird estimator, ML - Maximum likelihood, REML - Restricted maximum likelihood estimator, HE - Hedges estimator, HS - Hunter-Schmidt estimator, CA - Cochran ANOVA, PM or MP - Mandel-Paule estimator, P ​ M C ​ A 𝑃 subscript 𝑀 𝐶 𝐴 PM_{CA} - two-step Cochran ANOVA, P ​ M D ​ L 𝑃 subscript 𝑀 𝐷 𝐿 PM_{DL} -two-step DerSimonian-Laird, HM - Hartung-Makambi, SJ - Sidik-Jonkman, S ​ J C ​ A 𝑆 subscript 𝐽 𝐶 𝐴 SJ_{CA} - alternative Sidik-Jonkman, BM - Bayes Modal estimator; Estimators of δ 𝛿 \delta : SSW(H) - sample-size-weighted (Hedges 1982), EW - equal weights, IV - inverse-variance-weighted, HS - Hunter and Schmidt (1990) total-sample-size-weighted; Coverage of δ 𝛿 \delta : IV - confidence interval centered at inverse-variance-weighted estimator of δ 𝛿 \delta with z quantiles, IV+t - confidence interval centered at inverse-variance-weighted estimator of δ 𝛿 \delta with t quantiles, IV+HKSJ - Hartung-Knapp-Sidik-Jonkman confidence interval centered at inverse-variance-weighted estimator of δ 𝛿 \delta .

Study	SMD measure	$\delta$	${\tau }^2$	n and/or $\overline{n}$	K	${\hat{\tau }}^2$	Coverage of ${\tau }^2$	$\hat{\delta }$	Coverage of $\delta$
Viechtbauer 2005	Hedges′s d	0, 0.2 ,0.5 ,0.8	0, 0.01 ,0.025, 0.05 ,0.1	$\overline{n}=20,40,80,160,320$	5, 10, 20, 40, 80	DL		IV
				$n_i\sim N(\overline{n},\overline{n}/3)$		ML
				$n_{Ti}=n_{Ci}=n_i/2$		REML
						HE
						HS
Friedrich et al. 2008	Hedges′s d	0.2 ,0.5 ,0.8	0, 0.5	$n_T=n_C=10,100$	5, 10, 30	DL		IV	IV
Petropoulou and Mavridis 2017		0 & 0.5	0,0.01,0.05,0.5	$n_T=n_C\sim U(20,\mathrm{\hspace{0.33em}200})$	5,10,20, 30	20		IV	IV
						estimators
						of ${\tau }^2$
Langan et al. 2018	Hedges′s d	0.5	depends on $I^2$, 0%,15%	$n=40$, $n\sim U(40,400),$	2 ,3, 5,	REML		IV	IV
			30%, 45%, 60%,	$n=400$, $n\sim U(2000,4000)$,	10, 20, 30,	CA			IV+t
			75%, 90%, 95%	$n=40$+$N\sim U(2000,4000)$	50, 100	PM or MP			IV +HKSJ
				$n_T=n_C=n/2$		$P{M}_{CA}$
						$P{M}_{DL}$
						HM
						SJ
						$S{J}_{CA}$
Lin 2018	Cohen′s d	0, 0.2, 0.5, 0.8, 1	0,0.2,0.5	U(5,10), U(10,20),	5, 10, 20, 50	DL		IV	IV
	Hedges′s d			U(20,30), U(30,50)
				U(50,100), U(100,500),
				U(500,1000), $n_T=n_C$
Hamman et al. 2018	Hedges′s d	0,0.1,0.15,0.25	0,0.1,0.5,1,	$\overline{n}=4,6,8,10,12,$	5,10,15,25,35,	REML		IV	IV +HKSJ
		0.35,0.5,0.6,0.75,	2.5,5,10	14,16,20,25	45,55,75,100,125			SSW(H)
		1,1.25,1.5,2.5						EW
Marín-Martínez and		0.2,0.5,0.8	0,0.04,0.08,0.16,0.32	$\overline{n}=30,50,80,100$	5,10,20,40,100	DL		FE
Sánchez-Meca 2010				$n_T=n_C$		ML		IV
								HS

Table 2. Table 2: Combinations of parameters in the simulations for SMD

SMD	Equal study sizes	Unequal study sizes	Results in
			Appendix
$K$ (number of studies)	5, 10, 30	5, 10, 30
$n$ or $\overline{n}$ (average (individual) study size,	20, 40, 100, 250	30 (12,16,18,20,84),
total of the two arms)	30, 50, 60, 70	60 (24,32,36,40,168),
For $K=10$ and $K=30$, the same set of		100 (64,72,76,80,208),
unequal study sizes is used 2 or 6 times, respectively.		160 (124,132,136,140,268)
$q$ (proportion of each study in the control arm)	1/2, 3/4	1/2, 3/4
${\tau }^2$ (variance of random effect)	0(0.5)2.5	0(0.5)2.5	A1, A2
$\delta$ (true value of the SMD)	0, 0.2, 0.5, 1, 2	0, 0.2, 0.5, 1, 2	B1, B2

Equations

$${\delta }_i=\frac{{\mu }_{iT}-{\mu }_{iC}}{{\sigma }_i}.$$

$$g_i=J(m_i)\frac{{\bar{x}}_{iT}-{\bar{x}}_{iC}}{s_i},$$

$$J(m)=\frac{\Gamma \left(\frac{m}{2}\right)}{\sqrt{\frac{m}{2}}\Gamma \left(\frac{m-1}{2}\right)},$$

$$s_i^2=\frac{(n_{iT}-1)s_{iT}^2+(n_{iC}-1)s_{iC}^2}{n_{iT}+n_{iC}-2}.$$

$$v_i^2=\frac{n_{iT}+n_{iC}}{n_{iT}{n}_{iC}}+\left(1-\frac{(m_i-2)}{m_{i}J(m_i{)}^2}\right)g_i^2,$$

$$\frac{\sqrt{n_i{q}_i(1-q_i)}}{J(m_i)}{g}_i\sim t_{m_i}([n_i{q}_i(1-q_i){]}^{1/2}{\delta }_i).$$

$${\hat{\theta }}_i\sim N({\theta }_i,{\sigma }_i^2)\text{and}{\theta }_i\sim N(\theta ,{\tau }^2),$$

$${\hat{\theta }}_{RE}=\frac{\sum _{i=1}^K{\hat{w}}_i({\hat{\tau }}^2){\hat{\theta }}_i}{\sum _{i=1}^K{\hat{w}}_i({\hat{\tau }}^2)},$$

$$Q({\tau }^2)=\sum _{i=1}^K\frac{({\theta }_i-{\hat{\theta }}_{RE}{)}^2}{{\hat{\sigma }}_i^2+{\tau }^2}=E_{KDB}(Q).$$

$$Q({\tau }_L^2)=F_{Q;0.975}Q({\tau }_U^2)=F_{Q;0.025}$$

$$\widehat{Var}({\hat{\theta }}_{SSW})=\frac{\sum {\tilde{n}}_i^2(v_i^2+{\hat{\tau }}^2)}{(\sum {\tilde{n}}_i{)}^2},$$

Full text

Keywords between-study variance, heterogeneity, random-effects model, meta-analysis, mean difference, standardized mean difference

1 Introduction

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.

[FIGURE:]

2 Study-level estimation of standardized mean difference

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$ :

[TABLE]

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$).

3 Standard random-effects model

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,

[TABLE]

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

[TABLE]

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).

4 Point and interval estimation of $\tau^{2}$ by the Kulinskaya-Dollinger-Bjørkestøl method (KDB)

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

[TABLE]

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$:

[TABLE]

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

5 Sample-size-weighted (SSW) point and interval estimators of $\theta$

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

[TABLE]

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

6 Simulation study

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.

7 Methods of estimation of $\tau^{2}$ and $\delta$ used in simulations

Point estimators of $\tau^{2}$

• •

DL - method of DerSimonian and Laird (1986)

• •

J - method of Jackson (2013)

• •

KDB - method based on corrected null moment of $Q$ per Kulinskaya et al. (2011)

• •

MP - method of Mandel and Paule (1970)

• •

REML - restricted maximum-likelihood method

Interval estimators of $\tau^{2}$

• •

BJ - method of Biggerstaff and Jackson (2008)

• •

J - method of Jackson (2013)

• •

KDB - Q-profile method based on corrected null distribution of $Q$ per Kulinskaya et al. (2011)

• •

PL - profile-likelihood confidence interval based on $\hat{\tau}_{REML}^{2}$

• •

QP - Q-profile confidence interval of Viechtbauer (2007)

Point estimators of $\delta$

Inverse-variance-weighted methods with $\tau^{2}$ estimated by:

• •

DL

• •

J

• •

REML

• •

KDB

• •

MP

and

• •

SSW - weighted mean with weights that depend only on studies sample sizes

Interval estimators of $\delta$

Inverse-variance-weighted methods using normal quantiles, with $\tau^{2}$ estimated by:

• •

DL

• •

J

• •

KDB

• •

MP

• •

REML

Inverse-variance-weighted methods with modified variance of $\hat{\delta}$ and t-quantiles as in Hartung and Knapp (2001) and Sidik and Jonkman (2002)

• •

HKSJ (DL) - $\tau^{2}$ estimated by DL

• •

HKSJ KDB - $\tau^{2}$ estimated by KDB

and

• •

SSW KDB - SSW point estimator of $\delta$ with estimated variance given by (5.1) and t-quantiles

8 Results

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.

8.1 Bias in estimation of $\tau^{2}$ (Appendix A1)

The five estimators (DL, REML, J, MP, and KDB) have biases whose traces fan out from the same small positive bias at $\tau^{2}=0$. As $\tau^{2}$ increases, KDB remains positive and increases slowly; MP stays close to 0 (and slightly below); REML stays negative, with a negative slope; J stays negative, with a more-negative slope; and DL becomes increasingly negative, showing noticeable curvature. The value of $\delta$ has little effect on this pattern (except that the bias of DL and J has smaller magnitude when $\delta=2$).

As $n$ or $\bar{n}$ increases, the traces for KDB, MP, and REML flatten, and their bias is essentially 0 when $n$ or $\bar{n}=100$. The traces for J and DL become less steep, but substantial bias remains at $n=250$ (or $\bar{n}=160$).

As $K$ increases, the trace for KDB flattens somewhat, but the traces for the other estimators become steeper.

The traces for J and DL are slightly less steep when $q=.75$ than when $q=.5$.

In summary, the patterns of bias indicate a choice among the five estimators of $\tau^{2}$ (DL, REML, J, MP, and KDB). When $n\leq 40$, MP is closer to unbiased than KDB when $K=5$, the magnitudes of their biases are roughly equal when $K=10$, and KDB is closer to unbiased when $K=30$. When $n\geq 100$, MP, KDB, and REML are nearly unbiased. DL and J seriously underestimate $\tau^{2}$. The average of MP and KDB should be close to unbiased.

8.2 Coverage in estimation of $\tau^{2}$ (Appendix A2)

The five estimators (PL, QP, BJ, J, and KDB) share the feature that their coverage decreases as $\tau^{2}$ increases from 0 to 0.5. At $\tau^{2}=0$ all five have coverage $\geq.95$. KDB is highest (e.g., .99 when $q=.5$ and $n=20$), but it drops below .95 at $\tau^{2}=0.5$ or $\tau^{2}=1.0$ and remains slightly below .95. BJ is next highest (e.g., .98 when $q=.5$ and $n=20$), and it remains above .95 (say, .96 to .97) when $K=5$ and $K=10$. QP is close to .95 for $\tau^{2}\geq 0.5$. PL is between BJ and QP, and it remains above QP when $K=30$. The trace for BJ behaves quite differently when $K=30$ than when $K\leq 10$, decreasing steeply and linearly to around .77 at $\tau^{2}=2.5$ (when $q=.5$ and $n=20$). When $K\leq 10$, J is between BJ and QP; and when $K=30$, it also decreases linearly, but less steeply (e.g., to around .92 at $\tau^{2}=2.5$ when $q=.5$ and $n=20$).

Coverage does not change noticeably as $n$ or $\bar{n}$ increases, when $K\leq 10$. When $K=30$, the slopes of BJ and J become slightly less steep, and the traces of the other estimators move closer together and are closer to .95.

Setting aside the behavior of BJ and J when $K=30$, and of the other estimators when $n=20$ and $K=30$, the traces move closer together as $K$ increases.

The slopes of BJ and J when $K=30$ are less steep when $q=.75$ than when $q=.5$.

As $\delta$ increases, the coverage of QP at $K=30$ increases slightly; it is substantially closer to .95 when $\delta=2$.

In summary, all five interval estimators of $\tau^{2}$ have coverage substantially above .95 when $\tau^{2}=0$. When $\tau^{2}\geq 0.5$, QP is generally closest to .95. The unusual behavior of BJ (and, to a lesser extent, J) when $K=30$ adds to the evidence against it.

8.3 Bias and mean squared error in estimation of $\delta$ (Appendix B1)

When $\delta=0$, the bias of the six estimators (DL, REML, MP, KDB, J, and SSW) follows a single trace, close to 0, for all values of $\tau^{2}$. When $\delta>0$, SSW stays close to 0, and the others shift down, to increasingly negative bias, as $\delta$ increases, and their traces separate. For example, when $\delta=1$, $q=.5$, $n=20$, $K=10$, and $\tau^{2}=2.5$, the bias ranges from $-0.05$ (KDB) to $-0.07$ (DL), and MP, REML, and J (in that order) have intermediate values.

When $\delta\leq 0.5$, bias has little relation to $\tau^{2}$. When $\delta\geq 1$, however, the bias of the estimators other than SSW (especially DL) becomes increasingly negative as $\tau^{2}$ increases. (The plot for $n=20$ and $K=30$ in B1.17 shows an extreme example.)

Where bias is nonzero, increasing $n$ or $\bar{n}$ moves the traces toward (or to) 0, decreasing separation between them. Some plots (e.g.,  B1.37 and B1.39) show slight evidence of greater separation among traces when sample sizes are unequal and $\delta=2$.

For the most part, $K$ has little or no effect on bias. Some plots suggest that, where bias is nonzero, separation among traces increases as $K$ increases, especially from $K=10$ to $K=30$.

Bias does not differ noticeably between $q=.5$ and $q=.75$.

SSW essentially avoids the bias that we found in the inverse-variance-weighted estimators of $\delta$. To provide an additional measure of its performance (besides coverage, discussed below). we estimated the mean squared error of SSW and the best two inverse-variance-weighted estimators, KDB and MP. Appendix E1 includes figures that plot (versus $\tau^{2}$) the ratios MSE(SSW)/MSE(KDB) and MSE(SSW)/MSE(MP) for the five values of $\delta$, the two values of $q$, and $n_{i}=20,\;40,\;100,\;250$. For most situations the two ratios are essentially equal and differ little among values of $K$ and $q$. In most situations the traces are essentially flat as $\tau^{2}$ increases; otherwise, they curve downward as $\tau^{2}$ approaches 0. As $n$ increases, the ratios approach 1. For example, when $\delta=0$ and $K=5$, they decrease from around 1.1 when $n=20$ to nearly 1.0 when $n=250$. As $\delta$ increases ($\geq 0.5$), the ratios at small $\tau^{2}$ decrease. This pattern is first noticeable when $\delta=0.5$ and $n=20$ and $K=30$; and as $\delta$ increases, it becomes more pronounced at that combination of $n$ and $K$ and extends to larger $n$ (with $K=30$) and to $n\leq 40$ and $K=10$. When $\delta=2$, $q=.5$, $n=20$, and $K=30$, the traces for the two ratios are separate and curve up from around 0.55 at $\tau^{2}=0$ to slightly $<1$ when $\tau^{2}=2.5$. The patterns are similar for $q=.75$.

In summary, the bias of SSW is close to 0, and the other five estimators (DL, REML, J, MP, and KDB), which use inverse-variance weights, have greater (and negative) bias, amounting to 5 – 10% when sample sizes are small and $\delta\geq 1$. This bias increases as $\tau^{2}$ increases. SSW usually has slightly greater mean squared error than KDB and MP when $n$ is small, but its MSE can be substantially smaller, especially for small $\tau^{2}$.

8.4 Coverage in estimation of $\delta$ (Appendix B2)

Coverage of the estimators that rely on inverse-variance weights and normal critical values (DL, REML, MP, KDB, and J) is influenced most by $\tau^{2}$ ($=0$ versus $>0$) and $K$. At $\tau^{2}=0$ their coverage is around .97, but at $\tau^{2}=0.5$ (when $q=.5$ and $n=20$) it is mostly below .95: .90 to .91 when $K=5$, .92 to.94 when $K=10$, and .93 to .95+ when $K=30$. As $\tau^{2}$ increases, their coverage either is flat (REML, MP, and KDB) or decreases (DL, J). DL almost always has the lowest coverage, and the gap between it and J widens as $K$ increases.

Except for the effect of $K$ on SSW KDB at $\tau^{2}=0$ (above .99 when $K=5$, decreasing to .97 when $K=30$) and below-nominal coverage of HKSJ and HKSJ KDB in a region whose definition involves mainly $\delta$, $n$ (or $\bar{n}$), $K$, and $\tau^{2}$, the coverage of SSW and the HKSJ-type estimators is close to .95 for all $K$ and $\tau^{2}$. The challenging situations in that region generally involve $\delta=1$ or $\delta=2$, the smaller $n$ or $\bar{n}$, $K=10$ or $K=30$, and the smaller $\tau^{2}$. For example, when $\delta=2$ and $n=20$ or $\bar{n}=30$ and $K=30$, coverage can be as low as .84 when $\tau^{2}=0$.

When $\tau^{2}\geq 0.5$, coverage of DL, REML, MP, KDB, and J increases as $K$ increases, usually staying below .95.

The effect of $\delta$ on coverage is slight except for some situations involving $\delta=2$. When $n=20$ and $K=30$, all of the estimators except SSW KDB have low coverage at $\tau^{2}=0$, ranging from $<.82$ to .86. Their traces rise as $\tau^{2}$ increases; when $\tau^{2}=2.5$, KDB and HKSJ KDB are almost .94, and DL is .86 (up from .84). The pattern is similar when $n=40$ and $K=30$, but much reduced.

When $\tau^{2}\geq 0.5$, coverage decreases slightly as $n$ or $\bar{n}$ increases (except for SSW and the two HKSJ-type estimators); coverage is somewhat lower when sample sizes are unequal.

The plots show at most slight differences between $q=.5$ and $q=.75$.

In summary, except when $\delta=2$ and $K=30$, HKSJ and HKSJ KDB have coverage closest to .95; they differ little, and departures from .95 (toward lower coverage) are seldom serious. SSW KDB is rather conservative when $K=5$ and for other $K$ when $\tau^{2}=0$. Otherwise it provides reliable, albeit slightly conservative, coverage. When $\delta=2$ and $K=30$, SSW KDB is the best alternative. All of the estimators that use inverse-variance weights and critical values from the normal distribution (DL,REML, J, MP, and KDB) often have coverage substantially below .95.

9 Discussion: Practical implications for meta-analysis

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.

Funding

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

Appendices

• •

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

• •

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

10 Appendices

A1. Bias of point estimators of $\hat{\tau}^{2}$ for $\hat{\tau}^{2}=0.0(0.5)2.5$.

For bias 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 point estimators of $\tau^{2}$ 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)

A2. Coverage of of interval estimators of $\hat{\tau}^{2}$ for $\hat{\tau}^{2}=0.0(0.5)2.5$.

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)

B1. Bias and mean squared error of point estimators of $\hat{\delta}$ for $\hat{\tau}^{2}=0.0(0.5)2.5$.

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)

B2. Coverage of interval estimators $\hat{\delta}$ for $\hat{\tau}^{2}=0.0(0.5)2.5$.

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}$.