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

TL;DR

This study evaluates the performance of various estimators for between-study variance and overall effect in meta-analyses of mean difference through extensive simulations and an example.

Contribution

It introduces and compares multiple new and existing estimators for between-study variance and overall effect in meta-analyses, providing comprehensive empirical evidence.

Findings

WT method shows improved bias and coverage for $ au^2$ estimation.

Interval estimators vary in accuracy; some provide better coverage.

New estimators outperform traditional methods in certain scenarios.

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 mean difference (MD), we review 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; and a new method (WT) based on the improved approximation to the distribution of the $Q$ statistic by \cite{kulinskaya2004welch}), five interval estimators for $\tau^2$ (profile likelihood, Q-profile, Biggerstaff and Jackson, Jackson, and the new WT 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). We obtain empirical evidence from extensive simulations and an example.

Tables

Table 1. Table 1: Simulation studies on meta-analysis of MD. Estimators of τ 2 superscript 𝜏 2 \tau^{2} : DL - DerSimonian and Laird, ML - Maximum likelihood, REML - Restricted maximum likelihood, HE - Hedges, HS - Hunter and Schmidt; Intervals for τ 2 superscript 𝜏 2 \tau^{2} : QP - Q-profile, BJ - Biggerstaff and Jackson, J - Jackson, KBH - Knapp et al. (2006). Estimators of θ 𝜃 \theta : IV - inverse-variance-weighted; Intervals for θ 𝜃 \theta : all confidence intervals are centered at IV estimator of μ 𝜇 \mu . IV with z quantiles: HKSJ, HK2, SJ, SJ2 and KR. Confidence intervals using t quantiles with variance estimators by: HKSJ: Hartung and Knapp (2001), Sidik and Jonkman(2002);  HK2: Röver et al. (2015);  SJ and SJ2: Sidik and Jonkman (2006); KR: Kenward and Roger (1997).

Study	(${\sigma }_T^2$, ${\sigma }_C^2$)	$\mu$	${\tau }^2$	$n$ and/or $\overline{n}$	$K$	${\hat{\tau }}^2$	${\tau }^2$ intervals	$\hat{\mu }$	$\mu$ intervals
Viechtbauer 2005	$(10,10)$	0,1,2,4	0,0.125,0.25,0.5,1	$\overline{n}=20,40,80,160,320$	5, 10,	DL
	pooled variance			$n_i\sim N(\overline{n},{\left(\overline{n}/3\right)}^2)$	20, 40, 80	ML
				$n_{iT}=n_{iC}=n_i$		REML
						HE
						HS
Knapp et al. 2006	100, 10 across studies	0	1, 2.5, 5, 10, 20	20, 40	5, 10, 20, 50, 100
	pooled variance						KBH
Friedrich et al 2008	${\sigma }_T={\sigma }_C=\sigma =10,\mathrm{\hspace{0.33em}40},\mathrm{\hspace{0.33em}70}$	$0.2\sigma ,\mathrm{\hspace{0.33em}0.5}\sigma ,\mathrm{\hspace{0.33em}0.8}\sigma$	$0,\mathrm{\hspace{0.33em}0.5}\sigma$	$n=10,\mathrm{\hspace{0.33em}100}$, $n_T=n_C$;	5,  10,  30	DL		IV	IV
				$n_T=n_C/2$,  $n_T=2n_C$
				(with $n_T+n_C=2n$)
Lin 2018	${\sigma }_{iT}={\sigma }_{iC}\sim U(1,5)$	0, 0.5, 1, 2, 5	0, 0.25, 1	U(5,10), U(10,20),	5,  10,  20, 50	DL		IV	IV
	pooled variance			U(20,30),U(30,50)
				U(50,100),U(100,500),
				U(500,1000); $n_{Ti}=n_{Ci}=n_i/2$
IntHout et al. 2014	(1,1)	0	depends on $I^2$ and the $n_i$,	$n_T=n_C=n=25,\mathrm{\hspace{0.33em}50},$	2(1)10(5)20	DL		IV	IV
	pooled variance		$I^2=$0, 25%,  50%, 75%, 90%	100,  250,  500,  1000					HKSJ
				$\overline{n}$=100,  250, 500, 1000;
				25%, 50%, 75% large (=$10\times$ small)
Parlett and Riley 2017	one normal mean	1	0.01, 0.05, 0.1, 1	30	3, 5, 7, 10, 100	REML		IV	IV
	${\sigma }^2=0.1n$								KR
									HKSJ,  HK2
									SJ,  SJ2
Petropoulou and Mavridis 2017	(1,1)	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	17		IV	IV
	pooled variance		${\mu }_i\sim N(\mu ,{\tau }^2)$			estimators
			${\mu }_i\sim \mathrm{Laplace}(0,1)$			of ${\tau }^2$
Jackson 2013	one normal mean	0	0, 0.029, 0.069, 0.206, 1.302	implicit in ${\sigma }_i^2$	5; 10, 20, 40	DL	BJ
	$({\sigma }_1^2,\dots ,{\sigma }_5^2)=$						J
	(0.009, 0.046, 0.122, 0.265, 0.600)						QP

Table 2. Table 2: Data patterns in the simulations for MD

Parameter	Equal study sizes	Unequal study sizes	Full 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)
For $K=10$ and $K=30$, the same set of		60 (24,32,36,40,168),
unequal study sizes is used twice or six times,		100 (64,72,76,80,208),
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.01)0.1; 0(0.1)1	0(0.01)0.1; 0(0.1)1	A1, A2; A3, A4
${\sigma }_C^2,{\sigma }_T^2$(within-study variances)	(1,1), (1,2)	(1,1), (1,2)
$\mu$	0	0	B1, B2; B3, B4
${\tau }^2$ (variance of random effect)	0(0.1)1	0(0.1)1	A5, A6
${\sigma }_C^2,{\sigma }_T^2$ (within-study variances)	(10,10), (10,20)	(10,10), (10,20)
$\mu$	0	0	B5, B6

Equations

$${\mu }_i={\mu }_{iT}-{\mu }_{iC},estimated\; by{y}_i={\bar{x}}_{iT}-{\bar{x}}_{iC},$$

$$v_i^2={\hat{\sigma }}_i^2=s_{iT}^2/n_{iT}+s_{iC}^2/n_{iC}.$$

$${\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)},$$

$$g_i=\left(\frac{{\sigma }_{iT}^4}{n_{iT}^2{f}_{iT}}+\frac{{\sigma }_{iC}^4}{n_{iC}^2{f}_{iC}}\right);$$

$${\kappa }_1\approx I-1+2\sum _i{w}_i^2{g}_i{p}_i^2;{\kappa }_2\approx 2(I-1)+14\sum _i{w}_i^2{g}_i{p}_i^2.$$

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

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

$$\hbox{E}(Q)\approx I-1+2\sum_{i}w_{i}^{2}g_{i}p_{i}^{2}+\tau^{2}\big{(}W-W_{(2)}/W\big{)},$$

$${\hat{\tau }}_{CDL}^2=max\left(\frac{Q-(I-1)-2{\sum }_i{w}_i^2{g}_i{p}_i^2}{\hat{W}-{\hat{W}}_{(2)}/\hat{W}},0.\right)$$

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

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 mean difference as the effect measure.

If the studies can be assumed to have the same true effect, a meta-analysis can use 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 MD, 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 varies widely among effect measures. Veroniki et al. (2016) mention this only in passing (in Section 6.1) as a hypothetical possibility. Second, any review on the topic, such as Veroniki et al. (2016), currently can draw on only limited empirical information on the comparative performance of the methods. Table 1 gives the (short) list of previous simulation studies for MD. It shows that only two studies (Viechtbauer (2005) and Petropoulou and Mavridis (2017)) considered and compared several point estimators of $\tau^{2}$. Unfortunately, these and other studies (wrongly) used the pooled variance for MD. The pooled variance is equivalent to the unpooled variance, given by equation (2.1), only when the true variances and sample sizes are equal within the study. Also, simulations with one normally distributed effect measure are informative only under this scenario. Petropoulou and Mavridis (2017) use simulation to study 20 estimators of heterogeneity variance and their impact on coverage and length of 95% confidence intervals for the overall effect; but to assess bias of the estimators, they use mean absolute error, which is not a measure of bias; it is the linear counterpart of mean squared error.

Several studies have considered the quality of estimation of MD, but not estimation of heterogeneity variance. Friedrich et al. (2008) report extensive simulations for MD, but they use only the DL method to estimate $\tau^{2}$ and do not report on its quality. Lin (2018) provides similar simulations for MD. IntHout et al. (2014) study coverage for MD and combinations of different-sized studies, using the standard method and the Hartung-Knapp-Sidik-Jonkman (HKSJ) method based on the DL estimator of $\tau^{2}$ (Hartung and Knapp (2001), Sidik and Jonkman (2002)). Partlett and Riley (2017) provide another in-depth study of coverage based on restricted-maximum-likelihood (REML) estimation of $\tau^{2}$. They also recommend HKSJ confidence intervals.

Only two studies (Knapp et al. (2006) and Jackson (2013)) considered the quality of interval estimation of $\tau^{2}$ for MD.

To address this gap in information on methods of estimating the heterogeneity variance for MD, 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 an improved approximation to the distribution of the $Q$ statistic for MD (Kulinskaya et al., 2004) based on the Welch (1951) test. We also study coverage of confidence intervals for $\tau^{2}$ achieved by five methods, including the Q-profile method of Viechtbauer (2007), Q-profile methods based on the improved Welch-type 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. For comparison, we add an estimator (SSW) whose weights depend only on the sample sizes of the Treatment and Control arms. We study 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 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 mean difference effect measure is

[TABLE]

with variance $\sigma^{2}_{i}=\sigma_{iT}^{2}/n_{iT}+\sigma_{iC}^{2}/n_{iC}$, estimated by

[TABLE]

$s_{iT}^{2}$ and $s_{iC}^{2}$ do not depend on $\mu_{iT}$ and $\mu_{iC}$, so $\hat{\sigma}_{i}^{2}$ does not involve $\mu_{i}$. In the best-case scenario for traditional meta-analysis methods, for normal data, the sample means are independent of the sample variances (and therefore of inverse-variance-based weights). However, the relation of the between-study variance $\tau^{2}$ and the within-study variances $\sigma_{i}^{2}$ may affect quality of estimation. Sometimes the pooled variance is used instead of $v_{i}^{2}$ given by Equation (2.1). Unequal variances in the Treatment and Control arms, however, can adversely affect estimation (Kulinskaya et al. (2004)).

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 Welch-type and corrected DerSimonian-Laird methods

Because the $\hat{w}_{i}(\hat{\tau}^{2})$ in (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. However, this approximation is used in all moment methods for estimating $\tau^{2}$. As an alternative one can use an improved, effect-measure-specific approximation to the expected value of $Q$. Corrected Mandel-Paule-type moment methods for estimating $\tau^{2}$ equate the $Q$ statistic, with weights $\hat{w}_{i}(\tau^{2})$, to the first moment of an improved approximate null distribution. The corrected DerSimonian-Laird method estimates $\tau^{2}$ from the first moment of $Q$ under alternatives.

More-realistic approximations to the null distribution of $Q$ are available for several effect measures. These approximations do not treat the estimates $\hat{\sigma}_{i}^{2}$ as equal to $\sigma_{i}^{2}$. For MD, Kulinskaya et al. (2004) proposed an approximation based on the method of Welch (1951). This method calculates the first two corrected moments of $Q$, $\kappa_{1}=\hbox{E}[{Q}]$ and $\kappa_{2}=\hbox{Var}[{Q}]$, under the null hypothesis of homogeneity and then approximates the null distribution of $Q$ by an F distribution: $\hat{c}F_{K-1,\hat{f}_{2}}$ with matched moments. The estimated degrees of freedom $\hat{f}_{2}$ and the scale factor $\hat{c}$ are functions of $K$, the $n_{iT}$ and $n_{iC}$, and the $\hat{\sigma}_{iT}^{2}$ and $\hat{\sigma}_{iC}^{2}$.

To simplify notation, let $W=\sum{w_{i}}$, $W_{(k)}=\sum{w_{i}}^{k}$, and $p_{i}=1-w_{i}/W$, and let

[TABLE]

where $f_{ik}=n_{ik}-1$ is the number of degrees of freedom of for group $k$ of study $i$, $k=T,\;C$. Then the null moments of $Q$ for MD are (Kulinskaya et al., 2004),

[TABLE]

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

[TABLE]

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

We also propose a new WT 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. (2004) (i.e., the scaled F distribution with $K-1$ and $\hat{f}_{2}$ degrees of freedom based on the corrected first two moments 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.

We also propose another new method of estimating $\tau^{2}$ based on the improved first moment of $Q$ under alternatives, in the spirit of DerSimonian and Laird (1986). Under the alternatives,

[TABLE]

and the corrected DerSimonian-Laird (CDL) estimator is given by

[TABLE]

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

For comparison with the inverse-variance-weights estimators, we include 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 effective-sample-size-based weights were suggested in (Hedges and Olkin, 1985, p.110).

The interval estimators corresponding to SSW (SSW WT and SSW CDL) use the SSW point estimator as their center, and their 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{WT}}^{2}$ or $\hat{\tau}_{\mathit{CDL}}^{2}$, respectively.

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 MD, but gaps in evidence remain. For example, the possibility that the variances in the two arms may differ is rarely, if ever, reflected in simulations. Usually, pooled variances are used (Viechtbauer, 2005), or, equivalently, only the $y_{i}$ are simulated, as in IntHout et al. (2014) and Partlett and Riley (2017).

The range of values of $\mu_{i}$ in previous simulation studies varies from [math] in (Viechtbauer, 2005) to $56$ in Friedrich et al. (2008), but the value of $\mu_{i}$ is unimportant because $\hat{\sigma}_{i}^{2}$ does not involve $\mu_{i}$. In Viechtbauer (2005) $\tau^{2}$ and the $\sigma_{i}^{2}$ are commeasurable and vary from $0.0625$ to $1$. IntHout et al. (2014) set $\tau^{2}=C\bar{\sigma}^{2}$, where $\bar{\sigma}^{2}$ is the average within-study variance, and vary $C$ from $1/3$ to $9$. In Friedrich et al. (2008) the $\tau^{2}$ values, if non-zero, are much larger at $\tau^{2}=n\sigma^{2}/4$ for $n=10$ and $100$. Partlett and Riley (2017) also use comparatively large values of $\tau^{2}$, varying from $0.01$ to $1$, whereas $\sigma^{2}_{i}=0.1/30$ for $n=30$.

6.1 Design of the simulations

Our simulation study assesses the performance of six methods for point estimation of between-study variance $\tau^{2}$ (DL, REML, J, MP, WT and CDL) and five methods of interval estimation of $\tau^{2}$ (Q-profile-based methods corresponding to DerSimonian-Laird and Welch, 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 $\mu$ in the random-effects models for MD.

We vary six parameters: the between-study variance ($\tau^{2}$) and the within-study variances ($\sigma_{T}^{2}$ and $\sigma_{C}^{2}$), in addition to the number of studies ($K$), the total sample size ($n$), and the proportion of observations in the Control arm ($q$). We set the overall true MD $\mu=0$ because the estimators of $\tau^{2}$ do not involve $\mu$ and the estimators of $\mu$ are equivariant. Configurations we study are listed in Table 2.

To cover both small and large values of the ratio of within-study to between-studies variance, separately from the value of $\tau^{2}$, we use two series of within-study variances. We generate the within-study sample variances $s_{ij}^{2}$ ($j=T,\;C$) from chi-squared distributions as $\sigma_{ij}^{2}\chi_{n_{ij}-1}^{2}/(n_{ij}-1)$. We generate the estimated mean differences $y_{i}$ from a normal distribution with mean $\mu$ and variance $\sigma_{iT}^{2}/n_{iT}+\sigma_{iC}^{2}/n_{iC}+\tau^{2}$. We obtain the estimated within-study variances as $v_{i}^{2}=s_{iT}^{2}/n_{iT}+s_{iC}^{2}/n_{iC}$.

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 in the estimation of $\mu$ is $0.01$ (for $n=20$) or less for the first series, and $0.02$ or less for the second series of simulations. The simulation standard error for estimated coverage of $\tau^{2}$ or $\mu$ at the $95\%$ confidence level is roughly $\sqrt{0.95\times 0.05/10,000}=0.00218$.

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.

6.2 Analysis of the simulation data

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$, $\sigma_{C}^{2}$, and $\sigma_{T}^{2}$. 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$. The values of $\sigma_{C}^{2}$, and $\sigma_{T}^{2}$ consist of a cross of two factors, equal/unequal and small/large ($\sigma_{C}^{2}=1$ and $\sigma_{T}^{2}=1$, $\sigma_{C}^{2}=10$ and $\sigma_{T}^{2}=10$, $\sigma_{C}^{2}=1$ and $\sigma_{T}^{2}=2$, and $\sigma_{C}^{2}=10$ and $\sigma_{T}^{2}=20$). 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 $\mu$. Our summaries of results in Section 8 are based on examination of the figures in the corresponding Appendices.

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

Point estimators of $\tau^{2}$

• •

DL - method of DerSimonian and Laird (1986)

• •

J - method of Jackson (2013)

• •

MP - method of Mandel and Paule (1970)

• •

REML - restricted maximum-likelihood method

• •

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

• •

CDL - method based on corrected first moment of $Q$ under alternatives per Kulinskaya et al. (2004) and DerSimonian and Laird (1986).

Interval estimators of $\tau^{2}$

• •

BJ - method of Biggerstaff and Jackson (2008)

• •

J - method of Jackson (2013)

• •

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

• •

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

• •

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

Point estimators of $\mu$

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

• •

DL

• •

J

• •

REML

• •

MP

• •

WT

• •

CDL

and

• •

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

Interval estimators of $\mu$

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

• •

DL

• •

J

• •

MP

• •

REML

• •

WT

• •

CDL

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

• •

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

• •

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

and

• •

SSW WT - SSW point estimator of $\mu$ with estimated variance given by (5.1) with $\hat{\tau}^{2}=\hat{\tau}_{\mathit{WT}}^{2}$ and t-quantiles

• •

SSW CDL - SSW point estimator of $\mu$ with estimated variance given by (5.1) with $\hat{\tau}^{2}=\hat{\tau}_{\mathit{CDL}}^{2}$ 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}$ (Appendices A1, A3, and A5)

All of the estimators (DL, REML, J, MP, CDL and WT) have positive bias when $\tau^{2}=0$. In favorable situations (e.g., equal $n$’s, $q=.5$, and $\sigma_{C}^{2}=\sigma_{T}^{2}=1$, as in Figure A1.1), the bias of the estimators other than WT is slightly greater than 0 when $\tau^{2}=0$, and it decreases to near 0 as $\tau^{2}\rightarrow 1$ and becomes closer to 0 as $n$ increases. In some situations, however, the bias is much greater across the range of $\tau^{2}$; in the most extreme case in our simulations ($n=20$, $K=5$, $\sigma_{C}^{2}=10$, and $\sigma_{T}^{2}=20$; Figure A5.7), the bias at $\tau^{2}=0$ ranges from 1.6 (MP, J) to 1.8 (DL, REML) and decreases only to 1.4 to 1.7 at $\tau^{2}=1$. Generally, the bias decreases slightly as $K$ increases. Unbalanced arms ($q=.75$) magnify the bias, especially for the smaller values of $n$ and $\bar{n}$ and unequal variances, i.e. when $\sigma_{T}^{2}=2$ or $\sigma_{T}^{2}=20$. Among the estimators other than CDL and WT, DL and REML generally have the most bias, followed by J and MP.

The trace in the bias of CDL parallels that of DL, but it is considerably less biased than all the standard estimators. This favorable difference between CDL and the standard estimators is especially pronounced for small and unequal sample sizes. CDL is practically unbiased in the case of equal within-arm variances, and its bias is considerably less than that of DL in the most extreme cases of small and unequal sample sizes combined with unequal variances, as in Figure A5.7.

The bias of WT follows a strikingly different pattern. For small $\tau^{2}$ it is positive and smaller than (or equal to) that of the other estimators. As $\tau^{2}$ increases, the bias of WT becomes negative and takes increasingly more negative values, roughly linearly in $\tau^{2}$. The crossover point decreases as $K$ increases (e.g., 0.15 at $K=5$, 0.08 at $K=10$, and 0.05 at $K=30$, when $n=20$, $q=.5$, $\sigma_{C}^{2}=1$, and $\sigma_{T}^{2}=1$), but is substantially larger when $\sigma_{C}^{2}$ and $\sigma_{T}^{2}$ are large (Figures A5.1–A5.8). The slope against $\tau^{2}$ flattens substantially as $n$ (or $\bar{n}$) increases.

In summary, except for CDL and WT, the estimators of $\tau^{2}$ (DL, REML, J, and MP) have non-negligible positive bias, especially for small sample sizes ($n\leq 40$) and small values of $\tau^{2}$. Overall, CDL is the least biased, and WT is the least biased when the values of $\tau^{2}$ are considerably smaller than the within-study variances $\sigma^{2}_{i}$ (say, when $\tau^{2}\leq 0.2$ for $\sigma^{2}=1$, and when $\tau^{2}\leq 1$ for $\sigma^{2}=10$). All other estimators become acceptable for larger sample sizes $n\geq 100$.

8.2 Coverage in estimation of $\tau^{2}$ (Appendices A2, A4, and A6)

The relation between coverage of the interval estimators (PL, QP, BJ, J, and WT) and $\tau^{2}$ involves other variables. When $n$ or $\bar{n}$ is $\geq 100$, most of the interval estimators have coverage close to .95 when $\tau^{2}\geq 0.1$, but noticeably above .95 when $\tau^{2}=0$. (When $\sigma_{C}^{2}$ and $\sigma_{T}^{2}$ are large, PL often has coverage around .97.) When $K=30$ and $n=20$ (or $\bar{n}=30$) or $n=40$ (or $\bar{n}=60$), QP, PL, and BJ have quite low coverage at the smallest values of $\tau^{2}$ (e.g., .67 to .77 when $\tau^{2}=0$ and .83 to .88 when $\tau^{2}=0.1$, and $n=20$, $K=30$, $q=.75$, $\sigma_{C}^{2}=1$, and $\sigma_{T}^{2}=1$; Figures A4.3 and A2.3); the low coverage sometimes extends to most $\tau^{2}\in[0,1]$, and the departures are substantially greater for $q=.75$ than for $q=.5$. In those situations the coverage of WT is usually close to .95 and is seldom below .90. Otherwise, when $\tau^{2}>0$, the coverage of WT is usually between .95 and .96 and sometimes slightly greater. The impact of $\sigma_{C}^{2}$ and $\sigma_{T}^{2}$ (large vs. small, unequal vs. equal) is generally small.

In summary, none of the interval estimators of $\tau^{2}$ (PL, QP, BJ, J, and WT) consistently achieve coverage close to .95 (i.e., between .94 and .96). All have difficulty at $\tau^{2}=0$, usually overcoverage; the departures of PL extend to other small $\tau^{2}$, and its coverage is often greater than .96 but sometimes less than .94. Meta-analyses in which the studies have small sample sizes are challenging for PL, QP, BJ, and J, which in some situations have coverage well below nominal for all $\tau^{2}\in[0,1]$, especially when the number of studies is larger ($K=30$ vs. $K=5$ and $K=10$). Overall, WT comes closest to providing nominal coverage of $\tau^{2}$. (The contrast in behavior between the WT interval and point estimators is surprising, but the two are defined in different ways.)

8.3 Bias in estimation of $\mu$ (Appendices B1, B3, and B5)

Because the estimated MD and its estimated variance are independent, all the estimators of $\mu$ are practically unbiased in all situations.

8.4 Coverage in estimation of $\mu$ (Appendices B2, B4, and B6)

When $\tau^{2}\geq 0.1$ and $K\leq 10$, the methods that use critical values from the normal distribution (DL, REML, J, MP, CDL and WT) have coverage substantially below .95 ($<.90$ when $K=5$ and around .92 when $K=10$), and those coverages generally change only slightly as $\tau^{2}$ increases to 1. WT is lowest, in part because it underestimates $\tau^{2}$ (Figure B1.1). HKSJ and HKSJ WT have coverage close to .95. Coverage of SSW WT exceeds .95 and decreases toward .95 as $\tau^{2}$ increases. The estimators other than HKSJ- and SSW-type have coverage slightly $>.95$ when $\tau^{2}=0$. HKSJ and HKSJ WT are around .94 when $\tau^{2}=0$, and SSW WT is $>.99$ when $K=5$, decreasing to around .965 when $K=30$. For $0<\tau^{2}<0.1$, all estimators other than SSW-type have coverage below nominal for small sample sizes, considerably below for unequal sample sizes $\bar{n}\leq 60$ and $K\leq 10$ (Appendix B4).

When the studies’ $n$’s are equal, coverage of most profile-type estimators decreases as $n$ increases; this pattern is absent when the $n$’s are unequal. The coverage of WT increases to the level of the other profile-type estimators.

As $K$ increases, the coverage of the profile-type estimators approaches .95 (from below).

When $n=20$ or $\bar{n}=30$, the traces for the various estimators show more separation when $q=.75$ than when $q=.5$. This pattern is most noticeable when the $\sigma^{2}$’s are large and $K=30$. It is in these circumstances that SSW CDL is considerably better than SSW WT, achieving nominal coverage for larger values of $\tau^{2}$ or $K$, as in Figure B2.7.

Undercoverage is somewhat less when $\sigma_{C}^{2}=\sigma_{T}^{2}=10$ than when $\sigma_{C}^{2}=\sigma_{T}^{2}=1$. Coverage differs little between $\sigma_{T}^{2}=2$ and $\sigma_{T}^{2}=1$ and between $\sigma_{T}^{2}=20$ and $\sigma_{T}^{2}=10$.

In summary, HKSJ and HKSJ WT generally (but not uniformly) have the best coverage. Their coverage is not always within $\pm.01$ of .95; it may be considerably below nominal for $\tau^{2}<0.1$ when sample sizes are small; but in situations where clear differences separate the interval estimators, HKSJ and HKSJ WT are much closer to .95. DL, WT, MP, REML, and J exhibit very serious undercoverage when $K=5$ and nontrivial undercoverage when $K=10$.

9 Discussion: Practical implications for meta-analysis

The results of our simulations give a rather disappointing picture of the current state of meta-analysis. In brief:

Small sample sizes are rather problematic even for such a well-behaved effect measure as the mean difference, and meta-analyses that involve numerous small studies are especially challenging.

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 $50$ or less in $25\%$ of meta-analyses that used MD as an effect measure, and less than $110$ in $50\%$ of them Kulinskaya et al. (2014).

For MD, the between-study variance, $\tau^{2}$, is usually overestimated near zero, but the Welch-type method provides better point estimation for $\tau^{2}<0.1$, and reliable interval estimation across all values of $\tau^{2}$, $n$, and $K$. The estimates of $\mu$ are unbiased, and HKSJ intervals provide good coverage.

Arendacká (2012) and Liu et al. (2017) propose new confidence intervals for $\tau^{2}$ in the one-way heteroscedastic random-effects model. These intervals can be used directly in meta-analysis of means in noncomparative studies. Both publications include extensive simulations and compare their intervals with those of Knapp et al. (2006). Both proposals seem to do very well for normal distributions and very small sample sizes. It should be possible to extend these methods to MD in comparative two-arm designs; this extension will be pursued elsewhere.

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

This results in the unbiased estimation of $\theta$ for the mean differences because the estimated variances are independent of the estimated effects. An alternative approach 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. We prefer to use weights proportional to an effective sample size, $\tilde{n}_{i}=n_{iT}n_{iC}/n_{i}$. 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). A good estimator of $\tau^{2}$, such as MP or CDL can be used as $\hat{\tau}^{2}$. Further, confidence intervals for $\theta$ centered at $\hat{\theta}_{\mathit{SSW}}$ with $\hat{\tau}_{\mathit{WT}}^{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. MD: Plots for bias and coverage of $\tau^{2}$.

• •

Appendix B. MD: Plots for bias, mean squared error, and coverage of estimators of the mean difference $\mu$.

A1. Bias of $\hat{\tau}^{2}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For bias of $\hat{\tau}^{2}$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$ and a set of values of $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 $\tau^{2}$ are

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

A2. Coverage of $\hat{\tau}^{2}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For coverage of $\hat{\tau}^{2}$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$ , a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ and a set of values of $n$ (= 20, 40, 100, 250) 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)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jacksons interval)

A3. Bias of $\hat{\tau}^{2}$ for $\tau^{2}=0.0(0.01)0.1$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For bias of $\hat{\tau}^{2}$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.01)0.1$, a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ , and a set of values of $n$ (= 20, 40, 100, 250) 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)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

A4. Coverage of $\hat{\tau}^{2}$ for $\tau^{2}=0.0(0.01)0.1$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For coverage of $\hat{\tau}^{2}$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.01)0.1$ , a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ and a set of values of $n$ (= 20, 40, 100, 250) 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)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jacksons interval)

A5. Bias of $\hat{\tau}^{2}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=10$, $\sigma_{T}^{2}=10,\;20$.

For bias of $\hat{\tau}^{2}$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$, a value of $\sigma_{C}^{2}=10$, a value of $\sigma_{T}^{2}=10,\;20$ , and a set of values of $n$ (= 20, 40, 100, 250) 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)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

A6. Coverage of $\hat{\tau}^{2}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=10$, $\sigma_{T}^{2}=10,\;20$.

For coverage of $\hat{\tau}^{2}$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$ , a value of $\sigma_{C}^{2}=10$, a value of $\sigma_{T}^{2}=10,\;20$ and a set of values of $n$ (= 20, 40, 100, 250) 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)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jacksons interval)

B1. Bias and mean squared error of point estimators $\hat{\mu}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For bias of $\mu$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$, a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ , and a set of values of $n$ (= 20, 40, 100, 250) 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 WT estimator of $\tau^{2}$) use the above values of $\mu$ 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 $\mu$ are

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

• •

SSW (sample-size weighted)

B2. Coverage of $\hat{\mu}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For coverage of $\mu$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$, a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ , and a set of values of $n$ (= 20, 40, 100, 250) 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 $mu$ are the companions to the inverse-variance-weighted point estimators

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

and

• •

HKSJ (Hartung-Knapp-Sidik-Jonkman)

• •

HKSJ WT (HKSJ with WT 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}_{WT}$

B3. Bias and mean squared error of point estimators $\hat{\mu}$ for $\tau^{2}=0.0(0.01)0.1$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For bias of $\mu$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.01)0.1$, a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ , and a set of values of $n$ (= 20, 40, 100, 250) 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 WT estimator of $\tau^{2}$) use the above values of $\mu$ 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 $\mu$ are

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

• •

SSW (sample-size weighted)

B4. Coverage of $\hat{\mu}$ for $\tau^{2}=0.0(0.01)0.1$, $\sigma_{C}^{2}=1$, $\sigma_{T}^{2}=1,\;2$.

For coverage of $\mu$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.01)0.1$, a value of $\sigma_{C}^{2}=1$, a value of $\sigma_{T}^{2}=1,\;2$ , and a set of values of $n$ (= 20, 40, 100, 250) 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 $\mu$ are the companions to the inverse-variance-weighted point estimators

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

and

• •

HKSJ (Hartung-Knapp-Sidik-Jonkman)

• •

HKSJ WT (HKSJ with WT 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}_{WT}$

B5. Bias and mean squared error of point estimators $\hat{\mu}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=10$, $\sigma_{T}^{2}=10,\;20$.

For bias of $\mu$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.01)0.1$, a value of $\sigma_{C}^{2}=10$, a value of $\sigma_{T}^{2}=10,\;20$ , and a set of values of $n$ (= 20, 40, 100, 250) 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 WT estimator of $\tau^{2}$) use the above values of $\mu$ 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 $\mu$ are

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

• •

SSW (sample-size weighted)

B6. Coverage of $\hat{\mu}$ for $\tau^{2}=0.0(0.1)1.0$, $\sigma_{C}^{2}=10$, $\sigma_{T}^{2}=10,\;20$.

For coverage of $\mu$, each figure corresponds to a value of $\mu(=0,0.2,0.5,1,2)$, a value of $q(=.5,.75)$, a value of $\tau^{2}=0.0(0.1)1.0$, a value of $\sigma_{C}^{2}=10$, a value of $\sigma_{T}^{2}=10,\;20$ , and a set of values of $n$ (= 20, 40, 100, 250) 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 $mu$ are the companions to the inverse-variance-weighted point estimators

• •

DL (DerSimonian-Laird)

• •

REML (restricted maximum likelihood)

• •

MP (Mandel-Paule)

• •

WT (Corrected Mandel-Paule moment estimator based on Welch-type approximation for Q distribution)

• •

J (Jackson)

• •

CDL (Corrected DerSimonian-Laird)

and

• •

HKSJ (Hartung-Knapp-Sidik-Jonkman)

• •

HKSJ WT (HKSJ with WT 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}_{WT}$