An overdispersion model in meta-analysis

Abstract

We propose a general approach based on the concept of overdispersion for specification of a random effects model (REM) in meta-analysis. This approach is similar to that used in generalized linear models, and includes the traditional REM as a particular case. A key feature of the model is the interpretation of the multiplicative factor as an intra-class correlation parameter. We provide several motivating examples, discuss statistical inference, and compare the new and standard methods on two examples of published meta-analyses. Estimation of the overdispersion parameter in the proposed model is compared in simulations to that of the traditional between-studies variance in the case of normal means. For small values of heterogeneity, the coverage of the confidence intervals for the overdispersion parameter is more stable.

Keywords

Intra-class correlation between-study variability random effects model heterogeneity

1 Introduction

In a traditional random effects meta-analysis the variance of the estimated effect size parameter ${\hat{θ}}_{i}$ in the ith study (i = 1, …, K) is comprised of two components: a unique within-study variance $σ_{i}^{2}$ (generally assumed known) and a random effects variance τ². Assuming an underlying normal distribution for the ith study effect size estimate ${\hat{θ}}_{i}$ , the composite model is:

{\hat{θ}}_{i} \sim N (θ, σ_{i}^{2} + τ^{2}), i = 1, ..., K,

(1.1)

where θ is the overall effect size. More generally, other distributional assumptions can be made for the ith study effect size estimate:

{\hat{θ}}_{i} \sim F (\cdot; θ_{i}, σ_{i}^{2}) with θ_{i} \sim G (\cdot; θ, τ^{2}) .

(1.2)

Examples of non-normal distributions F and G are F binomial and G normal, or F normal and G a skewed normal or skewed t (see Lee and Thompson, 2008). Note, that the variance $σ_{i}^{2}$ may include a sample size normalization.

The model (1.1) is an additive random effects model (REM). However, in some circumstances, a multiplicative model

{\hat{θ}}_{i} \sim N (θ, ϕ σ_{i}^{2}), i = 1, ..., K

(1.3)

might be more appropriate. Note that the $σ_{i}^{2}$ in (1.1) and (1.3) are not the same. The multiplicative model (1.3) was introduced by Thompson and Sharp (1999). The parameter φ is a multiplicative random effects parameter that is constant over all studies. When τ² = 0 in (1.1) and φ = 1 in (1.3), the REM reduces to a fixed effects model, and both (1.1) and (1.3) are identical.

The basic REM (1.1) is based on the assumption that the K studies having particular characteristics are randomly chosen from a population of similar studies. The assumption of randomness, though generally accepted, is not subject to any testing procedure. For example, Peto (1987) noted that non-randomness occurs in designs with chronologically consecutive studies.

The proposed model is constructed by replacing φ in (1.3) by

ϕ_{i} = 1 + a (n_{i}) γ, i = 1, ..., K,

(1.4)

where a(n_i) are known linear functions of the sample sizes n_i, i = 1, …, K. This model has an advantage of permitting inflation or deflation (over- or underdispersion) of the variances, in contrast to the additive model in (1.1) which permits only inflation of the variances. A positive feature of this model is that γ can be interpreted as an intra-class correlation parameter or a transformation of it; see Section 3 for details.

The proposed overdispersion model (ODM) is motivated by overdispersion due to intra-class correlation in the normal model or in the binomial model. Some background on overdispersion is provided in Section 2. The proposed model is formally introduced in Section 3. Method of moments and maximum likelihood estimators of the correlation parameter γ are given in Section 4; two well-known examples are provided in Section 5; a simulation study for the case of normal means is described in Section 6; final comments are in Section 7.

2 Background on overdispersion

The idea of overdispersion arises naturally for data from a Poisson distribution. The mean and variance of a Poisson variable are equal, but in some cases the data may have a sample variance that is smaller or larger than the mean. The mean–variance specification has been used as a definition for the terms under- or overdispersion. The paper by Bliss (1953) is one of the earliest to discuss this concept. Cox (1983) describes a model in which the variance has an inflation or deflation factor. An extensive examination of overdispered binomial data is provided by Anderson (1988).

In the case of a one-parameter family, the notion of overdispersion is relatively clear because the variance is a known function of the mean (see Cox, 1983). When the postulated mean–variance relationship does not hold, overdispersion means inflation of the variance above the nominal value. For a beta-binomial model (see Section 3.2), the variance is inflated relative to the binomial variance. In general, the inflation may be by a constant multiplier or it may depend on covariates.

For a location-scale family, such as the normal distribution, the scale can take on any positive value. So there is no notion of data being overdispersed relative to such a model. However, in the case that the data is modelled by a distribution with a particular scale, it can be over- or underdispersed relative to that assumption. In the general meta-analytic context overdispersion means variance inflation relative to that of the fixed effects model. Thus, the REM (1.1) is overdispersed whenever the variance component τ² is non-zero. Similarly, the REM (1.3) is overdispersed when ϕ > 1.

In the more general context, overdispersion can arise from omitted unobserved variables. Regardless of origin of overdispersion, the main interest lies in the resulting marginal models, i.e., the composite models in which the mean and correlation functions are modelled directly and separately (see Pryseley et al., 2011). In meta-analysis, two possible mechanisms are clustering in the population and mixing resulting in a compound distribution. Different modelling approaches can lead to the same marginal model, so in general it is difficult to infer the precise cause, or underlying process, leading to the overdispersion. See Kim (2006), Hinde and Demétrio (1998) and references therein for further details.

3 Overdispersed models in meta-analysis

The proposed ODM was motivated by overdispersed normal data and by overdispersed binomial data when overdispersion is due to intra-class correlation.

The following examples make use of what has been called the intra-class correlation (ICC) matrix, namely, R = (r_jj_′) with r_jj= 1, r_jj_′ = ρ (j ≠ j′), j, j′ = 1, …, n, where ρ is the common correlation between two distinct observations. This n-dimensional matrix has one eigenvalue equal to 1 + (n – 1)ρ and n – 1 eigenvalues equal to 1 – ρ. Because R is positive semi-definite, ρ is constrained to the interval: –1/(n – 1) ≤ ρ ≤ 1. For even moderately sized samples the intra-class correlation is predominantly positive.

3.1 Overdispersed normal data

For participant j in study i, denote the outcome by Y_ij = θ + b_i + ε_ij, j = 1, …, K_i, i = 1, …, K. Further, let the random effect variables b_i be normally and independently distributed, $b_{i} \sim N (0, τ_{i}^{2})$ , and the errors ε_ij be normally and independently distributed with a common distribution within a study, i.e., $ε_{i j} \sim N (0, λ_{i}^{2}) .$ The random effects b_i and the errors ε_ij are mutually independent with $var (Y_{i j}) = σ_{i}^{2} = τ_{i}^{2} + λ_{i}^{2}$ and $cov (Y_{i j}, Y_{i j^{'}}) = τ_{i}^{2}$ for j ≠ j′. Thus, a key feature of this model is that any two participants j and j′from the same study have a correlation ρ_i depending only on the study:

ρ_{i} = corr (Y_{i j}, Y_{i j^{'}}) = τ_{i}^{2} / (τ_{i}^{2} + λ_{i}^{2}),

where ρ_i is the intra-class correlation coefficient (ICC). The n_i × n_i covariance matrix V_i for vector $Y_{i} = (Y_{i 1}, ..., Y_{i n_{i}})$ is

V_{i} = τ_{i}^{2} J_{n_{i}} + λ_{i}^{2} I_{n_{i}} = σ_{i}^{2} [(1 - ρ_{i}) I_{n_{i}} + ρ_{i} J_{n_{i}}],

(3.1)

where $J_{n_{i}}$ is the n_i × n_i matrix of 1’s and $I_{n_{i}}$ is the identity matrix of dimension n_i × n_i.

In this parameterization,

λ_{i}^{2} = σ_{i}^{2} (1 - ρ_{i}), τ_{i}^{2} = σ_{i}^{2} ρ_{i} .

(3.2)

For ρ_i > 0, the hierarchical model

Y_{i j} = θ_{i} + ε_{i j}, θ_{i} \sim N (θ, τ_{i}^{2})

(3.3)

with parameters $τ_{i}^{2}$ and $λ_{i}^{2},$ and the marginal model

Y_{i j} \sim N (θ, V_{i})

(3.4)

with parameters $σ_{i}^{2}$ and ρ_i, are equivalent. However, in the marginal model negative values of $τ_{i}^{2}$ are permissible because they correspond to the negative values of the ρ_i in the interval –1/(n_i – 1) < ρ_i < 0. The only restrictions on the ρ_i arise from the requirements of $τ_{i}^{2} + λ_{i}^{2} > 0$ and V_i to be positive definite. The study mean ${\bar{Y}}_{i}$ has the variance

var ({\bar{Y}}_{i}) = \frac{λ_{i}^{2}}{n_{i}} + τ_{i}^{2} = \frac{σ_{i}^{2}}{n_{i}} [1 + (n_{i} - 1) ρ_{i}] for ρ_{i} > - 1 / (n_{i} - 1) .

Assuming homogeneous correlations, i.e., ρ_i = ρ for all ICCs, we propose the meta-analytic REM of the one-sample study level statistics:

{\hat{θ}}_{i} \sim N (θ, \frac{σ_{i}^{2}}{n_{i}} [1 + (n_{i} - 1) ρ]) for ρ > - 1 / (\max {n_{i}} - 1)

(3.5a)

equivalently,

{\hat{θ}}_{i} \sim N (θ, \frac{λ_{i}^{2}}{n_{i}} [1 + n_{i} γ]) for γ = \frac{ρ}{1 - ρ} > - \frac{1}{\max {n_{i}}} .

(3.5b)

This model can also be called an overdispersed (for ρ > 0) or underdispersed (for ρ < 0) normal model, in which the estimate ${\hat{θ}}_{i} \sim N (θ, ϕ_{i} σ_{i}^{2} / n_{i})$ with φ_i = [1 + (n_i – 1)ρ] 0. Note that φ_i is linear in n_i. The form (3.5b) is more useful for estimation because the sample variances $s_{i}^{2}$ are estimates of $λ_{i}^{2} = σ_{i}^{2} (1 - ρ) .$

3.2 Overdispersed binomial data and the beta-binomial model

Consider the hierarchical model in which the group-level sums X_i of Bernoulli random variables have parameter π_i, i.e., X_i ∼ Bin(n_i, π_i); the probabilities π_i have a common beta distribution with parameters α, β, 0 < α, β < 1; that is π_i ∼ Beta(α, β). Then X_i has a beta-binomial distribution with the density

f (x_{i} | n_{i}, α, β) = (\begin{matrix} n_{i} \\ x_{i} \end{matrix}) \frac{B (α + x_{i}, β + n_{i} - x_{i})}{B (α, β)}, 0 \leq x_{i} \leq n_{i},

(3.6)

where B(u, v) is the ordinary beta function (Euler integral). Note that

E (X_{i}) = \frac{n_{i} α}{α + β}, var (X_{i}) = \frac{n_{i} α β (n_{i} + α + β)}{{(α + β)}^{2} (α + β + 1)} .

Reparametrize (3.6) with π = α/(α + β) and θ = 1/(α + β+ 1), 0 < π, θ < 1, or equivalently, α = π (1 – θ)/θ, β= (1 – π)(1 – θ)/θ. Then E(X_i) = n_i π and it follows that

var (X_{i}) = n_{i} π (1 - π) [1 + (n_{i} - 1) θ] .

Note the linear factor in the variance.

Alternatively, a marginal ICC-based model can be obtained as follows. Let Bernoulli variables Y_ij ∼ Bin(1, π) and corr(Y_ij, Y_ij_′) = ρ Then for $X_{i} = \sum_{j = 1}^{n_{i}} Y_{i j}$ , var(X_i) = n_iπ(1 – π)[1 + (n_i – 1)ρ]. This variance expression is the same as for the beta-binomial model with ρ instead of θ, but now ρ can be negative though restricted to the interval ρ > –1/[max{n_i} – 1]. Interestingly, the reparametrized density (3.6) is still valid for negative values of ρ as long as

\begin{array}{l} ρ > - \frac{π}{n_{\max} - π - 1} for 0 < π < \frac{1}{2}, \\ ρ > - \frac{1 - π}{n_{\max} - (1 - π) - 1} for 1 > π > \frac{1}{2} \end{array}

(3.7)

(Prentice, 1986). The beta-binomial model is widely used in both frequentist and Bayesian meta-analysis of binomial data (see e.g., Young-Xu and Chan, 2008). For a general development of hierarchical and marginal approaches to the generalized linear models see Pryseley et al. (2011).

3.3 Comparative studies

In a two-arm study with a treatment (T) and a control (C), hierarchical models can be used to compare directly two-sample effect measures θ_i, i = 1, …, K. Examples of these measures are standardized mean differences, risk differences or differences of log-odds ratios. As we shall see, a variety of marginal arm-level models are possible, and there is no simple way to choose among them.

Comparing treatment with control for normal data. Suppose the effect size of interest is the difference θ = µ_T – µ_C of population means in two arms of a study. Denote the sample sizes by n_T and n_C, the overall sample size n = n_T + n_C, and the variances by $σ_{T}^{2}$ and $σ_{C}^{2}$ for the treatment and control arms of a trial, respectively. Let the observations

Y_{i j k} \sim N (μ_{i C} + θ z_{k}, V_{i k}), i = 1, ..., K; j = 1, ..., n_{i k}; k = T or C,

(3.8)

where µ_iC is a control arm mean, z_k is an indicator variable for T or C, and the arm-level covariance matrices V_ik are given by (3.1) with an additional subscript k to denote an arm of the ith study. Denote correlations within the respective arms by ρ_iC and ρ_iT. Additionally, let ρ_iCT be the correlation across the arms of the ith study. The studies are assumed to be independent.

Suppose the difference of means µ_T – µ_C is estimated by the difference ${\bar{Y}}_{T} - {\bar{Y}}_{C}$ of the sample means. Suppressing the subscript i, the variance of the difference of the means is given by

var ({\bar{Y}}_{T} - {\bar{Y}}_{C}) = \frac{σ_{T}^{2}}{n_{T}} [1 - ρ_{T}] + \frac{σ_{C}^{2}}{n_{C}} [1 - ρ_{C}] + σ_{T}^{2} ρ_{T} + σ_{C}^{2} ρ_{C} - 2 σ_{T} σ_{C} ρ_{C T}

(3.9a)

= \frac{σ_{T}^{2}}{n_{T}} [1 + (n_{T} - 1) ρ_{T}] + \frac{σ_{C}^{2}}{n_{C}} [1 + (n_{C} - 1) ρ_{C}] - 2 σ_{T} σ_{C} ρ_{C T} .

(3.9b)

Both forms (3.9a) and (3.9b) provide different insights into the model. The standard REM has a random effect only in the treatment arm (see Thompson and Sharp (1999), models 6a,b), so that ρ_C = ρ_CT = 0, and (3.9a), (3.9b) reduce to

var ({\bar{Y}}_{T} - {\bar{Y}}_{C}) = {\frac{σ_{T}^{2}}{n_{T}} + \frac{σ_{C}^{2}}{n_{C}}} + \frac{σ_{T}^{2}}{n_{T}} [n_{T} - 1] ρ_{T}

(3.10a)

= \frac{σ_{T}^{2}}{n_{T}} [1 + (n_{T} - 1) ρ_{T}] + \frac{σ_{C}^{2}}{n_{C}}

(3.10b)

Another option is to assume independence between the arms (ρ_CT = 0), and equal correlations ρ_C = ρ_T = ρ within each arm, resulting in

var ({\bar{Y}}_{T} - {\bar{Y}}_{C}) = {\frac{σ_{T}^{2}}{n_{T}} + \frac{σ_{C}^{2}}{n_{C}}} (1 - ρ) + ρ (σ_{T}^{2} + σ_{C}^{2})

(3.11a)

= \frac{σ_{T}^{2}}{n_{T}} [1 + (n_{T} - 1) ρ] + \frac{σ_{C}^{2}}{n_{C}} [1 + (n_{C} - 1) ρ] .

(3.11b)

Define R = n_T/n_C to be the allocation ratio, and

v (R) = σ_{T}^{2} + σ_{C}^{2} + \frac{σ_{T}^{2}}{R} + R σ_{C}^{2},

(3.12)

The term v(R) is a precision measure of the difference between sample means. It is a function of the allocation ratio $R, σ_{T}^{2}$ and $σ_{C}^{2},$ but for given variances we are only concerned with v(R) as a function of R. Note that v(R) depends only on the allocation ratio and not on the study size, n = n_T + n_C.

For $\hat{θ} = {\bar{Y}}_{T} - {\bar{Y}}_{C}$ , equations (3.10a) and (3.10b) can be written as

var (\hat{θ}) = \frac{v (R)}{n} {1 + \frac{σ_{T}^{2}}{v (R)} (n_{T} - 1) (1 + R^{- 1}) ρ} .

(3.13)

From (3.12), $v (1) = 2 (σ_{T}^{2} + σ_{C}^{2})$ , so that equations (3.11a) and (3.11b) can be written as

var (\hat{θ}) = \frac{v (R)}{n} {1 + [\frac{n v (1)}{2 v (R)} - 1] ρ}

(3.14a)

= \frac{v (R) (1 - ρ)}{n} {1 + [\frac{n v (1)}{2 v (R)}] [\frac{ρ}{1 - ρ}]} .

(3.14b)

The form (3.14b) is more useful for estimation because the sample variances $s_{i k}^{2}$ estimate the error variances $λ_{i k}^{2} = σ_{i k}^{2} (1 - ρ)$ . A key point is that the variances are of the form

\frac{v (R_{i})}{n_{i}} (1 + a_{i} ρ), i = 1, ..., K,

(3.15)

where a_i = bn_i + c are linear functions of the sample sizes n_i. When $σ_{T}^{2} = σ_{C}^{2}$ , the values a_i = a(n_i, R_i) are constants depending on the study design and sample sizes, and not on the variances. This conclusion is also true for the case of balanced design (R = 1) in equations (3.14a) and (3.14b). Therefore for a moderately balanced data a choice of a_i = n_i/2 in (3.15) constitutes a simple and robust ODM in lieu of (3.11).

Formulation (3.15) also holds for other standard two-sample effect measures used in meta-analysis, such as, the standardized mean difference, the log mean ratio and the binary data based effect measures, discussed in the following.

Comparing treatment with control for binary data. Two-sample effect measures for binary outcomes compare the probabilities (risks), or functions of them, such as log-odds of an event in the treatment and the control arms of a study. The following four effect measures are widely used.

For binary outcomes with the probabilities p_T and p_C in the treatment and the control arm, denote the difference of risks by RD = p_T – p_C, the log relative risk by LRR = log(p_T/p_C), the log-odds ratio by LOR = log([p_T/(1 – p_T)]/[p_C/(1 – p_C)]), and the difference of arcsine-transformed risks by $ASD = 2 \arcsin (\sqrt{p_{T}}) - 2 \arcsin (\sqrt{p_{C}})$ . Then, assuming binomial distributions and independence within and between the arms,

var (\hat{RD}) \approx \frac{1}{n} (R + 1) {p_{C} (1 - p_{C}) + R^{- 1} p_{T} (1 - p_{T})},

(3.16a)

var (\hat{LRR}) \approx \frac{1}{n} (R + 1) {p_{C}^{- 1} - 1 + R^{- 1} (p_{T}^{- 1} - 1)},

(3.16b)

var (\hat{LOR}) \approx \frac{1}{n} (R + 1) {p_{C}^{- 1} + {(1 - p_{C})}^{- 1} + R^{- 1} (p_{T}^{- 1} + {(1 - p_{T})}^{- 1})},

(3.16c)

var (\hat{ASD}) \approx \frac{1}{n} {(R + 1)}^{2} / R .

(3.16d)

For overdispersed binomial distributions in one or both arms, still assuming independence between the two arms, equations (3.13) and (3.14a) with v(R) for each effect measure given by equations (3.16) can be obtained by the standard delta-method.

In general, in the fixed effects meta-analysis summary effects ${\hat{θ}}_{i} \sim N (θ, v_{i} (R_{i}) / n_{i}),$ i = 1, …, K. The variance v_i(R_i) in the ith study often depends on the effect of interest θ_i and the potential (vector) nuisance parameters which in the normal case are the within-arms variances. Assuming that the correlations within and between the arms are known multiples of the same parameter ρ, the general expression for the variances v(R_i)(1 + a_iρ)/n_i with a_i = b_in_i + c_i is applicable to all standard effect measures.

The proposed REM of meta-analysis for the two-sample study level statistics is

{\hat{θ}}_{i} \sim N (θ, \frac{v_{i} (R_{i})}{n_{i}} (1 + a_{i} γ)) for γ > - 1 / \max {a_{i}} .

(3.17)

As before, the values a_i = b_in_i + c_i are linear in the n_i. In the nonnormal case there may be additional nuisance parameters. The values b_i = b(R_i, θ_i, ψ_i) and c_i = c(R_i, θ_i, ψ_i) depend on the allocation ratios R_i, parameters of interest θ_i and (vector) nuisance parameters ψ_i. The model

{\hat{θ}}_{i} \sim N (θ, ϕ_{i} \frac{v_{i} (R_{i})}{n_{i}}) with ϕ_{i} = 1 + a_{i} γ

(3.18)

can also be called an over- or underdispersed normal model. The choice of a_i = n_i/v_i(R_i) combined with the restriction γ ≥ 0 results in the standard REM.

In the model (3.17) γ is a nuisance parameter independent of study i and the effect of interest θ, and it appears in the variance function as a multiplicative factor. Because γ > –1/max(a_i), underdispersion is also possible though usually not likely when sample sizes are moderate or large. For particular choices of coefficients b_i and c_i, the overdispersion in a meta-analysis can be interpreted through a marginal model involving the ICC. For binary effect measures, when the value of γ is less than 1, γ can be directly interpreted as the ICC ρ. For various measures based on the sample means for normal data, the parametrization γ = ρ/(1 – ρ) is more appropriate, cf (3.5b, 3.14b). Alternatively, for a simple ODM a_i = n_i can be used.

In (3.17), note that the variance of ${\hat{θ}}_{i}$ is v_i(R_i) times 1/n_i + a_iγ/n_i. When all the studies are large, the precision does not depend on the sample sizes n_i and is proportional to the precision v_i(R_i) of the original measurements. In contrast, in the standard REM, the variances $v_{i}^{2} (R_{i}) / n_{i} + τ^{2}$ result in relatively equal weights for large studies.

4 Estimation of γ

For known γ, the weights for estimating θ are $w_{i}^{*} = w_{i} / (1 + a_{i} γ),$ and the ‘corrected’ Q-statistic is given by $Q * = \sum w_{j}^{*} {({\hat{θ}}_{j} - {\bar{θ}}_{^{w *}})}^{2},$ where ${\bar{θ}}_{^{w *}} = \sum w_{i}^{*} θ_{i} / \sum w_{i}^{*}$ is a linear estimate of the effect size. The Q-statistic has a chi-square distribution with K – 1 degrees of freedom. In practice, γ needs to be estimated, and an estimated value $\hat{γ}$ (a consistent estimate) is substituted for γ in the weighted mean ${\bar{θ}}_{w *} .$ Rukhin (2009) obtains confidence intervals for the weighted mean ${\bar{θ}}_{w *} .$

The parameter γ can be estimated by a variety of methods, no one of which is uniformly the best, regardless of the criterion. We discuss here several alternatives. In this section the variances $σ_{i}^{2},$ and therefore the values a_i, are assumed known as is customary in meta-analysis methodology.

4.1 Method of moments estimation of γ

For Cochran’s statistic $Q = \sum w_{i} {({\hat{θ}}_{i} - {\bar{θ}}_{w})}^{2},$ the weights are inversely proportional to the variances, namely, $w_{i} = n_{i} / σ_{i}^{2}$ and ${\bar{θ}}_{w} = \sum w_{i} {\hat{θ}}_{i} / \sum w_{i} .$ Under the null hypothesis of no over- or underdispersion γ= 0, the Q-statistic has a chi-square distribution with K – 1 degrees of freedom, so that E(Q) = K – 1.

Under alternatives (3.17), $var ({\hat{θ}}_{i}) = w_{i}^{- 1} (1 + a_{i} γ),$ and for equal a_i = a, the Q-statistic has a scaled chi-square distribution: ${(1 + a γ)}^{- 1} Q \sim χ_{K - 1}^{2} .$ In the general case of unequal a_i, the distribution of Q is the same as that of a linear combination $\sum_{1}^{m} c_{i} Z_{i}^{2}$ , where each Z_i has a standard normal distribution, similar to the distribution in the standard REM (Biggerstaff and Jackson, 2008).

The exact non-null distribution of the Q-statistic is complicated. A number of approximations have been developed. One scheme is to approximate the distribution of $Z = \sum_{1}^{m} c_{i} Z_{i}^{2}$ by the distribution of a random variable cU^r, where c is a function of the coefficients c_i, for example, the geometric mean $c = Π c_{i}^{1 / m} .$ Other types of power means ${[\sum c_{i}^{r} / m]}^{1 / r}$ could also be used. The distribution of U is usually chosen to depend on one or two parameters, as in the gamma distribution. The choice of c, r and the two parameters then permits matching several moments. This idea is due to Tukey and Wilks (1946). For more detailed discussion and further references, see Section 3 of Olkin and Saner (2001).

To find the moments of Q under alternatives (3.17), one needs the covariances

cov ({\hat{θ}}_{i}, \bar{θ}) = (1 + a_{i} γ) / W,

(4.1)

cov ({\hat{θ}}_{j} - \bar{θ}, {\hat{θ}}_{i} - \bar{θ}) = δ_{i j} \frac{1 + a_{j} γ}{w_{j}} - \frac{1 + a_{j} γ}{W} - \frac{1 + a_{i} γ}{W} + \frac{1 + {\bar{a}}_{w} γ}{W},

(4.2)

where $W = \sum w_{i},$ ${\bar{a}}_{w} = \sum w_{i} a_{i} / W$ and δ_ij is the Kronecker delta. If the variances $σ_{i}^{2}$ are known,

E (Q) = K - 1 + (K \bar{a} - {\bar{a}}_{w}) γ,

(4.3)

where $\bar{a} = \sum a_{i} / K .$ The estimate of γ from (4.3) should satisfy the condition $\hat{γ} > - 1 / \max (a_{i}) .$ This condition is always met for the case of equal values a_i ≡ a, in which case E(Q) = (K – 1)(1 + aγ) > 0. For the sample equivalent, Q – (K – 1) ≥ $- (K \bar{a} - {\bar{a}}_{w}) /$ max(a_i), which may or may not be satisfied; in the latter case define $\hat{γ} = - 1 / a_{\max} .$ If only positive values of γ are acceptable, then $\hat{γ}$ can be truncated at zero.

The estimator in the spirit of random effects estimates of variance (e.g., DerSimonian and Laird, 1986) is

{\hat{γ}}_{M} = \max (\frac{Q - (K - 1)}{K \bar{a} - {\bar{a}}_{w}}, - \frac{1}{a_{\max}});

(4.4)

underdispersion is present for Q < K – 1.

Alternatively, we can use a two-sided test of γ= 0 based on Q. Overall this model appears to be more flexible than the standard REM because Q < K – 1 in a large proportion of meta-analyses, which is the case of negative γ.

4.2 Confidence intervals for γ

Several approaches to obtain a confidence interval for γ_M are described in the following, though others are available. The first is based on the gamma distribution with two parameters r and λ. Denote the gamma distribution function by G(r, λ) with density

g (z) = \frac{z^{r - 1} e^{- z / λ}}{λ^{r} Γ (r)}; r, λ, z > 0.

Following Biggerstaff and Tweedie (1997), the variance of Q under alternatives (3.17) is

var (Q) = 2 γ^{2} [\sum a_{i}^{2} - 2 \sum w_{i} a_{i}^{2} / W + {\bar{a}}_{w}^{2}] + 4 γ (K \bar{a} - {\bar{a}}_{w}) + 2 (K - 1) .

(4.5)

For equal a_i = a the variance (4.5) simplifies to var(Q) = 2(K – 1)(aγ + 1)² in accordance with the scaled $χ_{K - 1}^{2}$ distribution of Q in that case.

One approximation for the distribution of Q when the a_i differ is to match the first two moments of the gamma distribution:

r (γ) = \frac{{E (Q)}^{2}}{var (Q)} and λ (γ) = \frac{E (Q)}{var (Q)} .

(4.6)

For simplicity of notation we write r and λ for r(γ) and λ(γ). The estimator ${\hat{γ}}_{M}$ is distributed as a linear function of a gamma variable. Define $c = K \bar{a} - {\bar{a}}_{w}$ . Then the density of ${\hat{γ}}_{M}$ is

f (t; γ) = c \frac{λ^{r}}{Γ (r)} {[c t + K - 1]}^{r - 1} e^{- λ [c t + K - 1]}, λ > 0, r > 0, t > - (K - 1) / c .

(4.7)

This density function can be used to find a confidence interval for γ(Hedges and Olkin, 1985, pp. 54–56). Define the lower and upper functions L(γ) and U(γ) by

\begin{array}{l} L (γ) & = \int_{\hat{γ}}^{\infty} f (t; γ) d t & = \int_{λ Q}^{\infty} \frac{1}{Γ (r)} u^{r - 1} e^{- u} d u, \\ U (γ) & = \int_{- \frac{K - 1}{c}}^{\hat{γ}} f (t; γ) d t & = \int_{0}^{λ Q} \frac{1}{Γ (r)} u^{r - 1} e^{- u} d u . \end{array}

(4.8)

A 100(1 – α)% confidence interval for γ based on this approximate distribution is now obtained as a solution to L(γ) = α₁ and U(γ) = α_for α₁ + α₂ = α.

A second approach is to invert the Q test. When the correct weights $w_{i}^{*} = w_{i}^{*}$ (γ) = w_i/(1 + a_iγ) are used, the corrected Q-statistic given by $Q * (γ) = \sum w_{j}^{*} {({\hat{θ}}_{j} - {\bar{θ}}_{w *})}^{2}$ has the $χ_{K - 1}^{2}$ distribution. The confidence interval is constructed as

{γ > - 1 / a_{\max} : χ_{K - 1; α / 2}^{2} \leq Q * (γ) \leq χ_{K - 1; 1 - α / 2}^{2}} .

(4.9)

Viechtbauer (2007) shows that the standard REM confidence intervals for τ² based on this approach, named Q-profile, perform very well, better than the intervals using the method of moments and better than the restricted maximum likelihood confidence intervals described next.

4.3 Maximum likelihood estimation of γ

The likelihood function for the normal distribution with mean µ (omitting constants) is

l (γ, μ) = - \frac{1}{2} \sum \log (\frac{σ_{i}^{2}}{n_{i}} (1 + a_{i} γ)) - \frac{1}{2} \sum \frac{{(θ_{i} - μ)}^{2}}{σ_{i}^{2} (1 + a_{i} γ) / n_{i}} .

(4.10)

The ML estimates of γ and µ are given by the solution of the equations

\sum \frac{a_{i} {({\hat{θ}}_{i} - μ)}^{2}}{σ_{i}^{2} {(1 + a_{i} γ)}^{2} / n_{i}} = \sum \frac{a_{i}}{1 + a_{i} γ},

(4.11a)

\sum \frac{({\hat{θ}}_{i} - μ)}{σ_{i}^{2} (1 + a_{i} γ) / n_{i}} = 0.

(4.11b)

From (4.11b), the estimate of the mean is $μ_{ML} = \sum w_{i}^{*} {\hat{θ}}_{i} / W *,$ where $w_{i}^{*} = w_{i} /$ (1 + a_iγ) for $w_{i} = n_{i} / σ_{i}^{2}$ , and $W * = \sum w_{i}^{*} .$ For an initial choice of γ calculate µ_ML, which, when inserted in (4.11a), permits the solution in γ. This new value of γ is then used to get a new value of µ_ML. The process is then repeated. Convergence should be fairly rapid because the left-hand side and the right-hand side of (4.11a) are decreasing and convex functions of γ, so points of crossing will surface readily. Note that 1 + a_iγis constrained to be positive.

For equal a_i = a, (4.11b) simplifies and yields ${\hat{μ}}_{ML} = \sum w_{i} {\hat{θ}}_{i} / \sum w_{i}$ . Because ${\hat{μ}}_{ML}$ does not depend on γ, equation (4.11a) is simplified to $\sum w_{i} {({\hat{θ}}_{i} - μ)}^{2} = K (1 + a γ)$ . After substitution of ${\hat{μ}}_{ML}$ for µ the ML estimator of γ is $\hat{γ} = (Q - K) / K a$ . This estimator is biased because under the hypothesis γ = 0, EQ = K − 1 and the expectation $E (\hat{γ} | γ = 0) = - 1 / K a .$

4.4 Restricted maximum likelihood estimation (REML)

The restricted likelihood is

l_{R} (γ, μ) = - \frac{1}{2} \log (\sum w_{i}^{*}) - \frac{1}{2} \sum w_{i}^{*} {(θ_{i} - μ)}^{2} + \frac{1}{2} \sum \log (w_{i}^{*}),

(4.12)

for $w_{i}^{*} = w_{i} / (1 + a_{i} γ) .$ The derivative of $w_{i}^{*}$ with respect to γ is ${(w_{i}^{*})^{'}}_{γ} = - w_{i}^{*} a_{i} / (1 + a_{i} γ)$ and the REML equation for γ reduces to

{(W *)}^{- 1} \sum w_{i}^{*} \frac{a_{i}}{1 + a_{i} γ} + \sum w_{i}^{*} {(θ_{i} - μ)}^{2} \frac{a_{i}}{1 + a_{i} γ} = \sum \frac{a_{i}}{1 + a_{i} γ},

(4.13)

where $W * = \sum w_{i}^{*} .$ The mean µ_REML is obtained as $μ_{REML} = \sum w_{i}^{*} {\hat{θ}}_{i} / W * .$ The left- and right-hand sides of (4.13) are decreasing convex functions of γ, so an iterative procedure should readily yield a solution. For equal values, a_i = a,

{\hat{γ}}_{REML} = [Q - (K - 1)] / [(K - 1) a] = {\hat{γ}}_{M} .

For K = 2 studies,

{\hat{γ}}_{REML} = \frac{{(θ_{1} - θ_{2})}^{2} w_{1} w_{2} - (w_{1} + w_{2})}{w_{1} a_{2} + w_{2} a_{1}} = \frac{{(θ_{1} - θ_{2})}^{2} - σ_{1}^{2} / n_{1} - σ_{2}^{2} / n_{2}}{a_{1} σ_{1}^{2} / n_{1} + a_{2} σ_{2}^{2} / n_{2}} .

The REML confidence intervals are given by all values of γ which satisfy

l_{R} (γ) \geq l_{R} ({\hat{γ}}_{REML}) - χ_{1; 1 - α}^{2} / 2,

where $χ_{1; 1 - α}^{2}$ is the (1 − α) percentage point of chi-square distribution with 1 degree of freedom (Hardy and Thompson, 1996; Hedges and Olkin, 1985).

4.5 Mandel–Paule (MP) method

Under alternatives γ ≠ 0, the Q-statistic with the new weights $w_{i}^{*} = w_{i} / (1 + a_{i} γ)$ is

Q * (γ) = \sum \frac{w_{i}}{(1 + a_{i} γ)} {({\hat{θ}}_{i} - {\bar{θ}}_{w * (γ)})}^{2} .

(4.14)

This is a strictly decreasing function of $γ > - 1 / a_{\max},$ so a unique estimate of γ can be found from the estimating equation Q*(γ) = K − 1 given that a solution exists. This may not be the case due to the restricted range of γ; see the next subsection for details. This method was introduced originally by Mandel and Paule (1970) and we refer to it as the MP method. The motivation comes from the fact that, for the true weights, Q*(γ) has a chi-square distribution with K − 1 degrees of freedom and therefore the first moment is K − 1.

The MP algorithm was initially introduced in the context of inter-laboratory studies, but was then adopted for meta-analysis by Rukhin (2003) and DerSimonian and Kacker (2007). In the standard REM, the MP estimate of τ² coincides with the REML estimate given that the variances are equal to sample variances. As noted by Rukhin, Biggerstaff and Vangel (2000), the MP algorithm is also a generalized Bayes procedure for the common parameter θ. Desirable statistical properties and computational simplicity make this method of estimation our favourite.

Confidence intervals are constructed by the Q-profile method (4.9). A further improvement is possible by replacing K − 1 by an improved approximation for E₀(Q) correct to order 1/n. In general, improved approximations are effect-measure specific. The correction of Welch (1951) should be used for normal means; approximation by Kulinskaya, Dollinger and Bjørkestøl (2011) for standardized mean differences. Similarly, the Q-profile confidence intervals for γ should use improved approximations instead of the chi-square distribution.

Table 1

The Mandel–Paule estimates of the heterogeneity parameter γ, the Q-profile confidence intervals (4.9) for γ and the resulting estimates of the combined effect $\hat{θ}$

$Q_{\lim} = \lim_{γ \to - 1 / a_{\max}} Q * (γ)$	$\hat{γ}$	Q-profile CI for γ	$\hat{θ}$
$Q_{\lim} < χ_{K - 1; α / 2}^{2}$	–1/a_max	(–1/a_max, –1/a_max)	${\hat{θ}}_{1}$
$χ_{K - 1; α / 2}^{2} < Q_{\lim} < K - 1$	–1/a_max	$(- 1 / a_{\max}, {\hat{γ}}_{U})$	${\hat{θ}}_{1}$
$K - 1 < Q_{\lim} < χ_{K - 1; 1 - α / 2}^{2}$	${\hat{γ}}_{MP}$	$(- 1 / a_{\max}, {\hat{γ}}_{U})$	${\hat{θ}}_{ω *}$
$χ_{K - 1; 1 - α / 2}^{2} < Q_{\lim}$	${\hat{γ}}_{MP}$	$({\hat{γ}}_{L}, {\hat{γ}}_{U})$	${\hat{θ}}_{ω *}$

Source: Authors’ own

4.6 Estimating and testing γ in the case of underdispersion

As noted, the range of γ values is constrained from below by −1/a_max. As mentioned previously, the Q*(γ) statistic is a strictly decreasing function of γ given by (4.14). For simplicity, let a_max = a₁. When γ= −1/a₁, the corrected variance of the effect θ₁ given by $σ_{1}^{2} (1 + a_{1} γ)$ is zero and the weight $ω_{1}^{*}$ is undefined. Therefore, the value of the weighted mean ${\bar{θ}}_{w *}$ and the value of Q*(γ) at the point γ= −1/a₁ are both undefined. This can easily be remedied by defining these values as the limit values when γ → −1/a_max.

The limiting value of the weighted mean when γ → −1/a₁ is θ₁, and the limiting value of Q is Q_lim ≈ Q(−1/a₁ + ε) for a small value of ε. The value of Q_lim determines the estimates of θ and γ. A flowchart for different scenarios is provided in Table 1. In the table, a_max = max(a_i), ${\hat{γ}}_{MP}$ is the Mandel−Paule estimate of $γ, {\hat{γ}}_{L}$ and ${\hat{γ}}_{U}$ are the lower and upper confidence limits of the Q-profile interval, ${\hat{θ}}_{1}$ is the effect estimate corres¬ponding to the study with a_i = a_max and ${\hat{θ}}_{w *}$ is the weighted mean with the weights $w_{i}^{*} = w_{i} / (1 + a_{i} \hat{γ}) .$

The upper confidence limit ${\hat{γ}}_{U}$ is positive if and only if the standard Q-statistic value $Q = Q * (0) > χ_{K - 1; α / 2}^{2} .$ Thus a significant underdispersion occurs only when $Q < χ_{K - 1; α / 2}^{2}$ . This is a subcase of line 2 in Table 1.

When an improved approximation by a distribution F(⋅) is used to approximate the distribution of Q, the first moment of this distribution should be substituted for K − 1, and critical values of the distribution F(⋅) for the chi-square critical values in the first column of Table 1.

4.7 Confidence intervals for θ

For a weighted mean ${\bar{θ}}_{I V} = \sum w_{i} θ_{i} / \sum w_{i}$ with inverse variance weights w_i = n_i/v_i, the variance is $var ({\bar{θ}}_{I V}) = W^{- 1}$ , where W = ∑w_i. Because the true variances are unknown, θ is estimated by ${\hat{θ}}_{\hat{w}} = \sum_{i} {\hat{w}}_{i} {\hat{θ}}_{i} / \hat{W}$ , where ${\hat{w}}_{i}$ is a consistent estimator of w_i in the ith study. However, the variance of ${\hat{θ}}_{\hat{w}}$ is no longer W⁻¹, because the weights are estimated. In fact, W⁻¹ underestimates this variance, as noted by Li, Shi and Roth (1994) and by Rukhin (2009), so that the coverage of the conventional interval ${\hat{θ}}_{\hat{w}} \pm z {\hat{W}}^{- 1 / 2}$ is lower than the nominal level.

An alternative suggested by Rukhin (2009) is to estimate the variance of ${\hat{θ}}_{\hat{w}}$ by

ν = \sum \frac{w_{i}^{2}}{W (W - w_{i})} {({\hat{θ}}_{i} - {\bar{θ}}_{I V})}^{2} .

(4.15)

He shows that confidence intervals of the form ${\bar{θ}}_{w} \pm t_{α / 2} \sqrt{ν}$ perform well under normality. These confidence intervals need to be studied in the context of other effect measures used in meta-analysis.

5 Examples

5.1 Effects of diuretics on pre-eclampsia

A well-known meta-analysis of nine trials evaluated an effect of diuretics on pre-eclampsia (Collins, Yusuf and Peto, 1985). These data have been studied repeatedly, as for example in Hardy and Thompson (1996), Biggerstaff and Tweedie (1997) and Viechtbauer (2007). The basic data with odds ratios, and their logs are provided in Table 2a. The effect sizes indicate considerable heterogeneity in the populations in these trials. This heterogeneity can be seen from the differences in the overall risk $\tilde{p}$ of pre-eclampsia. Note that studies 5, 7 and 8 have the largest variances, due to very low overall risks in these studies, and studies 3, 4 and 9 have the smallest variances. There is a significant heterogeneity of risks of pre-eclampsia both in the treatment group (Q = 239.0 at 8 df) and in the control group(Q = 209.2); the slopes of the QQ plots (not shown) also suggest clear overdispersion in both arms, suggesting the model (3.11) as a potential fit. The data is fairly balanced, which suggests choice a_i = n_i/2.

Table 2a

Odds ratios (OR) of effect of diuretics on pre-eclampsia patients; n_C and n_T are the numbers of patients in the control and treatment (diuretics) groups, and x_C and x_T are the numbers of patients who developed pre-eclampsia in these groups. LOR is the log-odds ratio, and the pooled overall risk $\tilde{p}$ = (x_T + x_C)/(n_T + n_C).

Trial	x_C	n_C	x_T	n_T	OR	LOR	var(LOR)	$\tilde{p}$
1	14	136	14	131	1.043	0.042	42.613	0.105
2	17	134	21	385	0.397	–0.924	61.105	0.073
3	24	48	14	57	0.326	–1.122	18.692	0.362
4	18	40	6	38	0.229	–1.473	23.316	0.308
5	35	760	12	1011	0.249	–1.391	202.399	0.114
6	175	1336	138	1370	0.743	–0.297	39.599	0.015
7	20	524	15	506	0.770	–0.262	124.308	0.034
8	2	103	6	108	2.971	1.089	144.824	0.038
9	40	102	65	153	1.145	0.135	17.309	0.412
Total	345	3183	291	3759

Source: Author’s own. Columns 1–5 have been taken from Table 1 of Hardy and Thompson (1996).

Table 2b

Confidence intervals for log-odds ratios and for odds ratios for diuretics on pre-eclampsia example; ODM is the overdispersion model with a_i = n_i/2

Model	Heterogeneity	LOR	L	U	OR	L	U
FEM	τ ² = γ = 0	–0.398	–0.573	–0.223	0.672	0.564	0.800
REM	τ ²_DL = 0.23	–0.517	–0.916	–0.117	0.596	0.400	0.889
REM	τ ²_REML = 0.30	–0.518	–0.956	–0.080	0.596	0.384	0.923
ODM	γ_M= 0.0074	–0.428	–0.785	–0.071	0.652	0.456	0.932
ODM	γ_MP= 0.0158	–0.460	–0.912	–0.008	0.631	0.402	0.992

Source: Authors’ own

The Q-statistic value is Q = 27.265, the total sample size N = 6942 and ${\bar{a}}_{w} = \sum w_{i} n_{i} / (2 W) = 877.8$ , where w_i are the fixed effects model (FEM) weights listed in Table 2c. From formulae (4.4) and (4.14), with $\bar{a} = \bar{n} / 2$ , we obtain ${\hat{γ}}_{M} = 0.00743$ , and ${\hat{γ}}_{MP} = 0.01580$ (CI: 0.00227, 0.09213). For comparison, were we to use the variances given in (3.14), the values would be ${\hat{γ}}_{M} = 0.00736$ and ${\hat{γ}}_{MP} = 0.01619$ (CI: 0.00226, 0.09479). These values are directly interpretable as the estimated ICCs and their confidence limits. The conventional DerSimonian–Laird estimate of the variance component is $τ_{D L}^{2} = 0.23$ and $τ_{REML}^{2} = 0.30.$

To see the effect of these estimates of heterogeneity on the inference about the log-odds ratio, we compare the corresponding estimates for LOR and OR, and confidence intervals for LOR and OR in Table 2b. Note that all models yield confidence intervals that do not include zero for LOR or one for OR. The FEM yields the shortest confidence interval, as expected. The four REMs yield confidence intervals with widths ranging from 0.47 to 0.59 for OR and from 0.71 to 0.90 for LOR. The standard REM (with DL) and the ODM (with γ_M) yield the shortest confidence intervals.

Table 2c shows the weights attributed to each study using the five methods: FEM, REM (with DL), REM (with REML), ODM (with M) and ODM (with MP). As expected, the FEM shows a discrepancy with the four REMs. The ODM yields the largest estimated weights for studies 3, 4, 6 and 9, whereas for the REM studies 2, 5–7 and 9 yield the largest values. However, these differences are generally small, except for studies 5, 7 and 9.

The differences in weights are due to the fact that the precision of the studies 5 and 7 as measured by the inverse of the LOR variance is very low, whereas for study 9 it is rather high. These differences in precision have direct bearing on the ODM weights of the studies, whereas the REM weights of the studies are to large extent equalized.

Table 2c

of inverse variance weights (%) for diuretics on pre-eclampsia example. We compare three models, the FEM, the standard REM and the ODM. The estimators of heterogeneity used are given in the brackets: DerSimonian–Laird (DL), restricted maximum likelihood (REML), method of moments (M) and the Mandel–Paule (MP).

Trial	FEM	REM(DL)	REM(REML)	ODM(M)	ODM(MP)
1	5.0	10.7	10.9	10.4	10.7
2	6.8	11.9	11.9	9.6	8.9
3	4.5	10.2	10.4	13.4	16.3
4	2.7	7.9	8.3	8.6	11.0
5	7.0	12.1	12.0	3.8	3.1
6	54.5	17.0	15.9	20.5	16.2
7	6.6	11.8	11.9	5.7	4.8
8	1.2	4.5	5.1	2.7	2.9
9	11.8	13.9	13.6	25.1	26.0

Source: Authors’ own

5.2 Cochrane logo data

This example is a meta-analysis of seven randomized controlled clinical trials of a corticosteroids treatment to reduce mortality in premature babies (Crowley, Chalmers and Keirse, 1990). The data is used in the Cochrane Collaboration logo.¹ This example illustrates the case of line 2 in Table 1. The effect measure is the difference in risks Δ = p_T − p_C. Summary data is found in Table 3. For this data, Cochran’s Q-statistic is Q = 3.82 at 6 degrees of freedom, resulting in the p-value of 0.701. The DL and REML methods estimate the random variance component τ² by zero. The FEM indicates that the treatment significantly reduces mortality. The inverse variance estimate of the risk difference is −0.0547  0.0245.

There is significant heterogeneity of risks of mortality in the treatment (Q = 23.361 at 6 df), but not in the control (Q = 4.358) arms, suggesting the model (3.10). The a_i values are calculated from (3.13) with $σ_{T}^{2} = p_{T} (1 - p_{T})$ and v(R) given by (3.16). The first trial by Auckland has the largest a value, a₁ = 326.66, and therefore the range of γ is γ > −1/a₁ = −0.0031. The moment estimate of γ is ${\hat{γ}}_{M}$ = −0.0049. This value is less than −0.0031, and we chose the ${\hat{γ}}_{M} = - 0.0031$ . The value of Q_lim is approximated as Q(−1/a₁ + 0.00000001) = 5.47, which is less than K − 1 = 6 and greater than $χ_{K - 1; 0.025}^{2} = 1.237$ . The Mandel–Paule estimate of γ is not defined because Q − (K − 1) is negative across the whole range of the γ values. The lower confidence limit of the 95% Q-profile confidence interval is also equal to −0.0031, and the upper limit is 0.0376. The value of γ is not statistically different from zero. If we were to accept the underdispersion, and to use the value of ICC γ = −0.0031, our combined risk difference would be equal to the risk difference of −0.0439 in the Auckland trial.

Table 3

Risk difference for the effect of corticosteroids on the mortality in premature babies (Cochrane logo data); n_T and n_C are the sample sizes in treatment (T) and control (C), and x_T and x_C are the numbers of events in the arms of each study; d is the estimated difference in proportions of events between the two arms of each study.

Study	n_T	x_T	n_C	x_C	d
Auckland	532	36	538	60	–0.0439
Block	69	1	61	5	–0.0675
Doran	81	4	63	11	–0.1252
Gamsu	131	14	137	20	–0.0391
Morrison	67	3	59	7	–0.0739
Papageorgiou	71	1	75	7	–0.0792
Tauesch	56	8	71	10	0.0020
Totals	1007	67	1004	120

Source: http://www.cochrane.org/about-us/history/our-logo

6 Simulation results

In this section we provide a simulation study of the properties of the Q-profile confidence intervals for the ICC parameter γ in comparison to those for the random effects variance τ². Because the Q-profile confidence intervals for τ² were found by Viechtbauer (2007) to perform the best in the standard REM we limit our study to these intervals.

In addition we present a simulation study of confidence intervals for the combined effect θ. We include both conventional and improved intervals based on the variance (equation 4.15) combined with critical values from the t-distribution (Rukhin, 2009) for both REM and ODM, when γ and τ² are estimated by the moment and the MP methods. Also included are confidence intervals for γ obtained by the improved MP method based on the Welch correction (Welch, 1951) to the distribution of Q, denoted by ${\hat{γ}}_{W M P} .$

The size n_i of the studies was generated from a normal distribution with mean n and variance n/4 rounded to the nearest integer and was left-truncated by 3. To make the two models comparable, we generate estimates ${\hat{θ}}_{i}$ from equation (2.1b) with variances $λ_{i}^{2} = 1.$ Without loss of generality we have chosen the true effect θ = 0. Thus ${\hat{θ}}_{i}$ ~ (0, 1/n_i + γ). Estimated variances ${\hat{σ}}_{i}^{2}$ are generated from a scaled chi-square distribution with n_i − 1 degrees of freedom, that is, $(n_{i} - 1) {\hat{σ}}_{i}^{2} \sim χ_{n_{i} - 1}^{2} .$

The following configurations from Viechtbauer (2007) are included in the simulations: K = (10; 20; 30; 50; 80), n = (10; 20; 40; 80; 160), γ between 0 and 1 in steps of 0.1, and θ is set to 0. A total of 10 000 iterations were run for each combination. Therefore the standard error of the empirical coverage probabilities is at most 0.005.

Achieved confidence levels of confidence intervals for γ and τ² at the nominal confidence level of 0.95 for the number of studies K = 10, 30 and 80 are presented in Figure 1. The Q-profile confidence intervals for τ² behave well for τ² ≥ 0.1. Their actual coverage can be considerably lower (for small n) or higher (for large n) than the nominal confidence level under the null.

Figure 1

Achieved confidence levels of the γ and τ² confidence intervals at the nominal confidence level of 0.95. The methods used are represented as follows: squares (Q-profile interval for τ² based on χ² distribution) triangles (Q-profile interval for γ based on χ² distribution) and circles (Q-profile interval for γ with Welch correction). Here K is the number of studies and N is the total sample size of each study.

Confidence levels of the Q-profile confidence intervals for γ based on the χ² approximation to the distribution of Q are considerably lower than the nominal level. The coverage improves for large sample sizes n.

The Q-profile confidence intervals for γ based on the Welch approximation to the distribution of Q are reasonably good across all values of γ even for small sample sizes, and are recommended. The main advantage of these intervals in the comparison to the Q-profile confidence intervals for τ² is their stability when γ is small.

Achieved confidence levels of the confidence intervals for the combined effect θ at the nominal confidence level of 0.95 for the number of studies K = 10, 30 and 80 are presented in Figure 2. Coverage of the traditional confidence intervals with weights incorporating ${\hat{γ}}_{M}$ or ${\hat{γ}}_{M P}$ is considerably lower than nominal. The same is true, to a lesser extent, for traditional confidence intervals with weights incorporating ${\hat{τ}}_{D L}^{2}$ and ${\hat{τ}}_{M P}^{2} .$

Coverage is greatly improved for the confidence intervals obtained by the method of Rukhin (2009). These intervals, centered at the weighted means with the weights incorporating ${\hat{τ}}_{D L}^{2}$ or ${\hat{γ}}_{W M P},$ with the variances calculated by (4.15) and the critical values from the t distribution with K − 1 degrees of freedom, are recommended.

7 Discussion

The tradition of reporting meta-analytic results is in terms of FEM or REMs. However, whereas the FEM is well defined, there can be alternative representation of the REMs. We have pursued a model that uses an inflationary or deflationary factor applied to the FEM variance, and demonstrated its effect on the estimation of several parameters. The main rationale of our model is the interpretation of the multiplicative factor through the intra-class correlation parameter. To examine more visibly the effects of using different models, we give results for two well-known data examples.

The multiplicative model (1.3) introduced by Thompson and Sharp (1999) was used in a number of publications (Moreno et al. 2009, 2012) to adjust for a possible small study bias. Moreno et al. (2012) argue that this is a FEM because the studies with larger sample sizes have larger weights. When all the studies are large, the traditional REM results in approximately equal weights. In contrast, our ODM apportions larger weights to studies with higher precision.

An important advantage of the ODM is the possibility of interpretation and modelling of the heterogeneity in meta-analysis through the ICC parameter as it relates to the studies and their ‘populations’. We assumed that the ICC parameter was constant between the studies. When the number of studies is limited, as often is the case in meta-analysis, this assumption is necessary for estimation. For a large number of studies, meta-regression can be used for modelling the disparities among the homogeneous in respect to ICC subgroups of studies. Such disparities constitute one of the possible but so far neglected reasons for heterogeneity between studies.

Figure 2

Achieved confidence levels of the combined effect θ confidence intervals at the nominal confidence level of 0.95. The variance estimators used in weights are represented as follows: open squares rhomboids open triangles reverse triangles filled triangles ( with the Rukhin (2009) correction) and filled squares with the Rukhin (2009) correction). Here K is the number of studies and N is the total sample size of the two arms of each study.

We have performed simulations to compare estimation quality for the overdispersion parameter in the new model versus the traditional variance component in the case of normal means. The results are encouraging: the estimation quality in our simulations is considerably better for small values of the overdispersion parameter. Our results agree with other simulations, such as, Viechtbauer (2007), with respect to the poor quality of traditional estimators of the between-studies variance τ² for small values of heterogeneity. It is worth noting that we used an improved Welch approximation (Welch, 1951) in order to achieve satisfactory coverage of the confidence intervals for γ. Additional simulations need to be carried out for other popular effect measures, such as, the standardized mean difference and the odds ratio.

Footnotes

Acknowledgements

The authors thank two anonymous referees, an associate editor and the editor for helpful suggestions on the presentation of the contents of this article.

Note

References

Anderson

(1988) Some models for overdispersed binomial data. Australian Journal of Statistics, 30, 125–48.

Biggerstaff

Jackson

(2008) The exact distribution of Cochrans’ heterogeneity statistic in one-way random effects meta-analysis. Statistics in Medicine, 27, 6093–110.

Biggerstaff

Tweedie

(1997) Incorporating variability in estimates of heterogeneity in the random effects model in meta-analysis. Statistics in Medicine, 16, 753–68.

Bliss

(1953) Fitting the negative binomial distribution to biological data. Biometrics, 9, 176–96.

Collins

Yusuf

Peto

(1985) Overview of randomised trials of diuretics in pregnancy. British Medical Journal, 291, 97–104.

Cox

(1983) Some remarks on overdispersion. Biometrika, 70, 269–74.

Crowley

Chalmers

Keirse

(1990) The effects of corticosteroid administration before preterm delivery: an overview of the evidence from controlled trials. Br J Obstet Gynaecol, 97, 11–25.

DerSimonian

Kacker

(2007) Random-effects model for meta-analysis of clinical trials: an update. Contemporary Clinical Trials, 28, 105–14.

DerSimonian

Laird

(1986) Meta-analysis in clinical trials. Controlled Clinical Trials, 7, 177–88.

10.

Hardy

Thompson

(1996) A likelihood approach to meta-analysis. Statistics in Medicine, 15, 619–29.

11.

Hedges

Olkin

(1985) Statistical methods for meta-analysis. Orlando: Academic Press.

12.

Hinde

Demétrio

CGB

(1998) Overdispersion: models and estimation. Computational Statistics & Data Analysis, 27, 151–70.

13.

Kim

(2006) Overdispersion. In Encyclopedia of Statistical Sciences. New York: John Wiley & Sons, DOI: 10.1002/0471667196.ess0701.pub2.

14.

Kulinskaya

Dollinger

Bjørkestøl

(2011) Testing for homogeneity in meta-analysis I. The one parameter case: Standardized mean difference. Biometrics, 67, 203–12.

15.

Lee

Thompson

(2008) Flexible parametric models for random-effects distributions. Statistics in Medicine, 27, 418–34.

16.

Shi

Roth

(1994) The bias of the commonly-used estimate of variance in meta-analysis. Communications in Statistics - Theory and Methods, 23, 1063–85.

17.

Mandel

Paule

(1970) Interlaboratory evaluation of a material with unequal number of replicates. Analytical Chemistry, 42, 1194–97.

18.

Moreno

Sutton

Ades

Stanley

Abrams

Peters

Cooper

(2009) Assessment of regression-based methods to adjust for publication bias through a comprehensive simulation study. BMC Medical Research Methodology, 9, 2.

19.

Moreno

Sutton

Thompson

Ades

Abrams

Cooper

(2012) A generalised weighting regression-derived meta-analysis estimator robust to small study effects and heterogeneity. Statistics in Medicine, 31, 1407–17.

20.

Olkin

Saner

(2001) Approximation for trimmed Fisher procedures in research synthesis. Statistical Methods in Medical Research, 10, 267–76.

21.

Peto

(1987) Why do we need systematic overviews of randomized trials? Statistics in Medicine, 6, 233–40.

22.

Prentice

(1986) Binary regression using an extended beta-binomial distribution, with discussion of correlation induced by covariate measurement errors. Journal of the American Statistical Association, 81, 321–27.

23.

Pryseley

Tchonlafi

Verbeke

Molenberghs

(2011) Estimating negative variance components from gaussian and non-gaussian data: a mixed models approach. Computational Statistics & Data Analysis, 55, 1071–85.

24.

Rukhin

(2003) Two procedures of meta-analysis in clinical trials and interlaboratory studies. Tatra Mountains Mathematical Publications, 26, 155–68.

25.

Rukhin

(2009) Weighted means statistics in interlaboratory studies. Metrologia, 46, 323–31.

26.

Rukhin

Biggerstaff

Vangel

(2000) Restricted maximum likelihood estimation of a common mean and the Mandel-Paule algorithm. Journal of Statistical Planning and Inference, 83, 319–30.

27.

Thompson

Sharp

(1999) Explaining heter¬ogeneity in meta-analysis: a comparison of methods. Statistics in Medicine, 18, 2693–708.

28.

Tukey

Wilks

(1946) Approximation of the distribution of the product of beta variables by a single beta variable. Annals of Mathematical Statistics, 17, 318–24.

29.

Viechtbauer

(2007) Confidence intervals for the amount of heterogeneity in meta-analysis. Statistics in Medicine, 26, 37–52.

30.

Welch

(1951) On the comparison of several mean values: an alternative approach. Biometrika, 38, 330–36.

31.

Young-Xu

Chan

(2008) Pooling overdispersed binomial data to estimate event rate. BMC Medical Research Methodology, 8, 358.