A Meta-Analysis for Simultaneously Estimating Individual Means with Shrinkage, Isotonic Regression and Pretests

Abstract

Meta-analyses combine the estimators of individual means to estimate the common mean of a population. However, the common mean could be undefined or uninformative in some scenarios where individual means are “ordered” or “sparse”. Hence, assessments of individual means become relevant, rather than the common mean. In this article, we propose simultaneous estimation of individual means using the James–Stein shrinkage estimators, which improve upon individual studies’ estimators. We also propose isotonic regression estimators for ordered means, and pretest estimators for sparse means. We provide theoretical explanations and simulation results demonstrating the superiority of the proposed estimators over the individual studies’ estimators. The proposed methods are illustrated by two datasets: one comes from gastric cancer patients and the other from COVID-19 patients.

1. Introduction

Meta-analysis is a statistical method used to summarize results from published studies [1]. While it has been applied to all areas of science, meta-analysis has been a particularly important tool in educational studies [2,3] and medical studies [4,5]. Especially, meta-analysis has played an important role in studies on the impact of COVID-19 [6,7,8].

The most basic method in meta-analyses uses a fixed-effect model, where all of the studies have a common mean [5,9]. This model has long been recognized in mathematical statistics and stratified sampling designs, where the common mean estimator has rigorously been investigated (pp. 55–103 of [10]; [11,12,13]).

Another basic model is the random-effects model [14], where the studies have different means. This model imposes normally distributed random-effects to account for between-study deviations from the common mean [15]. In either model, the goal of meta-analyses is to estimate the common mean by combining the estimators of individual means.

However, in some scenarios, the common mean is undefined or uninformative. In these scenarios, the focus of meta-analyses can be on individual studies’ means. While individual means are estimated by published studies, they could be improved by statistical techniques. For instance, Raudenbush and Bryk [16] pointed out the possibility for improving individual studies’ estimators in meta-analyses by using empirical Bayes estimators under the normal random-effects. Schmid [17] also considered similar Bayes estimators that improve upon individual studies’ estimators. Röver and Friede [18,19] tried to improve an individual estimator by another individual estimator. All these Bayes estimators are defined by posterior means under the normal priors (random-effects). Evidently, the posterior means are related to but different from the individual studies’ estimates. However, there are some cases where the normal priors (random-effects) do not apply (Section 2.2). This motivates us to develop non-Bayesian estimation.

In this article, we consider some improved estimation methods of individual means without normal random-effects. First, we will explain how to employ shrinkage estimators to obtain improved estimators of individual means. Unlike the existing Bayes estimators, our frequentist shrinkage estimators do not involve any prior distribution. Next, we will provide theoretical and numerical results to show the superiority of the proposed estimators over the individual studies’ estimators. Finally, we will analyze two datasets to illustrate the proposed methods.

This article is organized as follows. Section 2 provides a background including a review of meta-analyses. Section 3 proposes the estimators of individual means and studies their theoretical properties. Section 4 conducts simulation studies, and Section 5 analyzes real datasets. Section 6 concludes the article.

2. Background

2.1. Meta-Analysis

This subsection reviews the basic models for meta-analysis under the normal distribution, including the fixed-effect model and random-effects model [1,5].

First, we introduce some mathematical notations and basic assumptions for meta-analyses. Let $G>0$ be the number of studies ($G$ stands for Groups) in a meta-analysis. For $i=1,2,\dots ,G$, the i-th study provides an estimate $Y_i$ for an unknown mean ${\mu }_i$ with a known variance ${\sigma }_i^2>0$. This implies that each study provides ${\sigma }_i$ as the standard error (SE) of the estimate $Y_i$. Hence, in meta-analyses, the variance ${\sigma }_i^2$ is known. We note that there could be a different setting, where the variance is unknown [20,21].

Next, we further assume the normal model $Y_i~N\left({\mu }_i,{\sigma }_i^2\right)$. This model’s special case, the fixed-effect model, imposes the common mean assumption. That is, $\mu \equiv {\mu }_1=\dots ={\mu }_G$. The random-effects model imposes normally distributed means given by ${\mu }_i~N\left(\mu ,{\tau }^2\right)$ for $i=1,2,\dots ,G$. In either case, the main goal of meta-analyses is to estimate$\mu$ by optimally combining $Y_i$s [1,5,9].

In a meta-analysis, observed data consist of $\left\{Y_i;i=1,2,\dots ,G\right\}$. An example is the two-sample t-test for a continuous outcome, where $Y_i$ is the (standardized) mean difference estimator from the i-th study, and ${\sigma }_i^2$ is the pooled variance estimator [5]. Another example is the $2\times 2$ table analysis for association between an event and a risk factor, where $Y_i$ is the logarithm of the risk ratio (or odds ratio), and ${\sigma }_i^2$ is estimated variance [5]. Another example is a survival analysis for censored data, where $Y_i$ can be the logarithm of an estimate of relative risks from the Cox model. Without loss of generality, we assume that ${\mu }_j=0$ corresponds to the null value.

Finally, we describe the main idea of meta-analyses. The idea is to optimally combine the results from the $G$ studies, based on either the fixed-effect model or random-effects model. The first model leads to the fixed-effect meta-analysis under the assumption $\mu \equiv {\mu }_1=\cdots ={\mu }_G$. The second model leads to the random-effects meta-analysis under the assumption ${\mu }_i~N\left(\mu ,{\tau }^2\right)$. In either case, the target is to estimate the common (overall) mean $\mu$ by

$${\widehat{\mu }}^{\mathrm{Fixed}}=\frac{{{\displaystyle \sum }}_{i=1}^G{\left({\sigma }_i^2\right)}^{-1}{Y}_i}{{{\displaystyle \sum }}_{i=1}^G{\left({\sigma }_i^2\right)}^{-1}},$$

$${\widehat{\mu }}^{\mathrm{Rand}}=\frac{{{\displaystyle \sum }}_{i=1}^G{\left({\sigma }_i^2+{\widehat{\tau }}^2\right)}^{-1}{Y}_i}{{{\displaystyle \sum }}_{i=1}^G{\left({\sigma }_i^2+{\widehat{\tau }}^2\right)}^{-1}},$$

where the value of ${\widehat{\tau }}^2$ is an estimator from the data [5,14].

2.2. Does the Estimation of the Common Mean Make Sense?

Given the data of $\left\{Y_i: i=1,2,\dots ,G\right\}$, one usually conducts the fixed-effect meta-analysis and/or random-effects meta-analysis. The Q test for homogeneity can help select between the fixed-effect and random-effects model [1,5]. If the fixed-effect model is rejected (i.e., the random-effects model is selected), one might explore a systematic source of heterogeneity by meta-regression. One might also test the goodness-of-fit of the normal random-effects model [22]. In any case, statistical inference for the common mean depends on some statistical models.

However, the above statistical models do not always fit the data at hand. If the studies have a monotonic-effect, e.g., ${\mu }_1\le \dots \le {\mu }_G$, the model is neither fixed nor random. Thus, there is no general way to estimate the common parameter unless some specific models are imposed. If a covariate is available to explain some structure among ${\mu }_i$s, a meta-regression may be used (Chapter 20 of [1]). However, we do not wish to assume the presence of covariates to conduct meta-regression in this article.

While the common mean estimate can be an informative summary of individual estimates, this single value rarely describes the whole results. Thus, meta-analysis typically displays individual estimates $\left(y_1,\dots ,y_G\right)$, providing an opportunity for learning the data. For instance, forest plots show the individual estimates $\left(y_1,\dots ,y_G\right)$ along with their 95% confidence intervals (CIs): $\left(y_1,\dots ,y_G\right)\pm 1.96\times \left({\sigma }_1,\dots ,{\sigma }_G\right)$. Funnel plots show $\left(y_1,\dots ,y_G\right)$ against$\left({\sigma }_1,\dots ,{\sigma }_G\right)$; see [1,5,23,24] for these plots. In some cases, one is interested in the ${\mu }_i$ that has the largest size [16]. Therefore, reporting individual estimates is an important part of meta-analysis.

2.3. Individual Estimates Could Be Improved

Despite the potential advantages of improving the individual estimators by the shrinkage methods [11,12,25,26], the development of such estimators has rarely been discussed in the literature of meta-analyses. In other words, the goal of meta-analyses has been the estimation of the common mean by disregarding the individual means. This causes a gap between meta-analyses and shrinkage estimation methods. Therefore, there is room for improving the standard individual estimators $\left(y_1,\dots ,y_G\right)$ for $\left({\mu }_1,\dots ,{\mu }_G\right)$ in some scenarios.

We raise scenarios where the common mean is undefined or uninformative.

The first example is the case of the monotonic-effect, ${\mu }_1\le {\mu }_2\le \dots \le {\mu }_G$, where the common mean seems to be not well-defined. In this case, an alternative is to perform meta-regression with the aid of observable covariates. We however do not assume the availability of covariates in this article.

The second example is the case of sparse normal means [27,28], where most of ${\mu }_i$s are equal to zero, e.g., $\left({\mu }_1,\dots ,{\mu }_{10}\right)=\left(-5,0,0,0,0,0,0,0,0,5\right)$. This is relevant to meta-analyses that fail to reject many null hypotheses. In this case, the interest may be in the estimation of the nonzero means rather than the common mean.

3. Proposed Estimators

We propose a variety of shrinkage estimators for simultaneously estimating individual means$\left({\mu }_1,\dots ,{\mu }_G\right)$, and discuss their properties. In particular, we theoretically prove that the proposed estimators have better precision than the individual studies’ estimates $\left(y_1,\dots ,y_G\right)$ in terms of a mean squared error criterion under some conditions.

The problem setting is as follows: Let $Y_i~N\left({\mu }_i,{\sigma }_i^2\right)$ be an estimator for an unknown mean ${\mu }_i$ under a known variance ${\sigma }_i^2>0$ for $i=1,2,\dots ,G$. As mentioned in Section 2, the variance ${\sigma }_i^2$ is known from the $i$-th study. Observed data consist of $\left\{Y_i;i=1,2,\dots ,G\right\}$. We regard $\mathit{Y}\equiv \left(Y_1,\dots ,Y_G\right)$ as a simultaneous estimator of $\mathit{\mu }\equiv \left({\mu }_1,\dots ,{\mu }_G\right)$. That is, we regard $\mathit{Y}$ as an observable random vector to estimate the unknown vector $\mathit{\mu }$.

Let $\mathit{\delta }\left(\mathit{Y}\right)\equiv \left({\delta }_1\left(\mathit{Y}\right),\dots ,{\delta }_G\left(\mathit{Y}\right)\right)$ be a simultaneous estimator for $\mathit{\mu }$. We state that $\mathit{\delta }\left(\mathit{Y}\right)$ improves upon $\mathit{Y}$ if

$$E\left[{\displaystyle \sum }_{i=1}^G\frac{{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right]<E\left[{\displaystyle \sum }_{i=1}^G\frac{{\left(Y_i-{\mu }_i\right)}^2}{{\sigma }_i^2}\right]=G, \exists \left({\mu }_1,\dots ,{\mu }_G\right),$$

(1)

$$E\left[{\displaystyle \sum }_{i=1}^G\frac{{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right]\le E\left[{\displaystyle \sum }_{i=1}^G\frac{{\left(Y_i-{\mu }_i\right)}^2}{{\sigma }_i^2}\right]=G, \forall \left({\mu }_1,\dots ,{\mu }_G\right).$$

(2)

We then call $\mathit{\delta }\left(\mathit{Y}\right)$ an improved estimator of $\mathit{Y}$.

We call $\mathit{\delta }\left(\mathit{Y}\right)=\mathit{Y}$ the standard unbiased estimator because of the unbiasedness $E\left(\mathit{Y}\right)=\mathit{\mu }$. It is a benchmark estimator for comparing different estimators. Recall that $\mathit{Y}$ is simply the collection of the individual studies’ estimates. If “$\forall \left({\mu }_1,\dots ,{\mu }_G\right)$” in the criterion (2) holds only for a restricted parameter space, $\mathit{\delta }\left(\mathit{Y}\right)$ is only locally improved.

The criteria (1) and (2) employ the weighted mean squared errors (WMSEs). The inverse variance weights guarantee the equal contributions of all studies to the estimation error. The weighs also make it convenient to apply the classical decision theory. The special case of ${\sigma }_1=\dots ={\sigma }_G=1$ gives the total MSE (TMSE), defined as $E\left[{{\displaystyle \sum }}_{i=1}^G{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2\right]$ [11,25,29]. We use the WMSEs for theoretical convenience while the TMSE criterion is also a relevant criterion. We will use both the WMSE and TMSE in our numerical assessments.

Indeed, there are infinitely many estimators $\mathit{\delta }\left(\mathit{Y}\right)$ that improve upon $\mathit{Y}$, including very complex and unrealistic ones [25]. In addition, it is quite easy to find an estimator that locally improves upon $\mathit{Y}$, such as $\mathit{\delta }\left(\mathit{Y}\right)=\mathbf{0}$. While the problem of deriving/assessing estimators has been intensively discussed in the statistical decision theory, it has rarely been appreciated in meta-analytical settings. The goal of this article is to facilitate the applications of shrinkage estimators in the context of meta-analyses.

In the subsequent sections, we will introduce three estimators that help reduce the WMSE and TMSE by shrinking $\mathit{Y}$ toward a restricted space of $\left({\mu }_1,\dots ,{\mu }_G\right)$. Section 3.1 will discuss the shrinkage towards the zero vector $\mathbf{0}=\left\{\left({\mu }_1,\dots ,{\mu }_G\right):{\mu }_1=\cdots ={\mu }_G=0\right\}$, the most traditional shrinkage scheme. Section 3.2 will consider the shrinkage toward $\mathbf{0}$ under constraints $\left\{\left({\mu }_1,\dots ,{\mu }_G\right):{\mu }_1\le \dots \le {\mu }_G\right\}$. Section 3.3 will explore the shrinkage towards the sparse space $\{\left({\mu }_1,\dots ,{\mu }_G\right):{{\displaystyle \sum }}_{i=1}^{G}I\left({\mu }_i\ne 0\right)<G\}$.

3.1. Shrinkage Estimation for Means

To estimate $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_G\right)$, we propose the James–Stein (JS) estimator of the form

$${\mathit{\delta }}^{\mathrm{JS}}\equiv \left({\delta }_1^{\mathrm{JS}},\dots ,{\delta }_G^{\mathrm{JS}}\right)\equiv \left(1-\frac{G-2}{{{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2}\right)\mathit{Y}.$$

This estimator is a modification of the original JS estimator [29] that was derived under the unit variances (${\sigma }_i=1$ for $\forall i$). See Appendix A.1 for the details. The JS estimator reduces variance by shrinking the vector $\mathit{Y}$ toward $\mathbf{0}$ while it produces bias. The degree of shrinkage is determined by the factor $\left(G-2\right)/({{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2)$ that usually ranges from 0 (0% shrinkage) to 1 (100% shrinkage), and occasionally becomes greater than 1 (overshrinkage).

It can be shown in Appendix A.1 that ${\mathit{\delta }}^{\mathrm{JS}}$ has the following WMSE

$$E\left[{\displaystyle \sum }_{i=1}^G\frac{{\left({\delta }_i^{\mathrm{JS}}\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right]=G-{\left(G-2\right)}^{2}E\left[\frac{1}{{\chi }_G^2\left(\lambda \right)}\right],$$

where ${\chi }_G^2\left(\lambda \right)$ is a random variable having a noncentral ${\chi }_{}^2$-distribution with the noncentral parameter $\lambda ={{\displaystyle \sum }}_{i=1}^G{\mu }_i^2/{\sigma }_i^2$ and the degrees of freedom $G$. Thus, ${\mathit{\delta }}^{\mathrm{JS}}$ has a smaller WMSE than $\mathit{Y}$ when $G\ge 3$. Indeed, the WMSE is minimized at $\mathit{\mu }=\mathbf{0}$ at which the WMSE is $G-1/\left(G-2\right)$ by the inverse moment of the central ${\chi }_{}^2$-distribution with$\lambda =0$. Thus, the JS estimator gains the greatest advantage if all the individual means are zero. This gain is appealing for meta-analyses for small individual effects (true means close to zero). Even if $\mathit{\mu }\ne \mathbf{0}$, the JS estimator has a smaller WMSE than $\mathit{Y}$. The reduction of the WMSE diminishes as $\lambda$ departs from zero.

One might ask where the special formula of the JS estimator comes from. The JS estimator ${\mathit{\delta }}^{\mathrm{JS}}$ can be derived as an empirical Bayes estimator under the prior ${\mu }_i~N\left(0,{\sigma }_i^2{\tau }^2\right)$:

$${\mathit{\delta }}^{\mathrm{JS}}=\left(1-\widehat{\left(\frac{1}{1+{\tau }^2}\right)}\right)\mathit{Y},$$

where the shrinkage factor $\widehat{\left(1/\left(1+{\tau }^2\right)\right)}\equiv \left(G-2\right)/{{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2$ is the estimator of $1/\left(1+{\tau }^2\right)$. See Appendix A.2 for the detailed derivations. Thus, if ${\mu }_i~N\left(0,{\sigma }_i^2{\tau }^2\right)$, the JS estimator minimizes the Bayes risk, and hence, is the optimal estimator under the prior.

A minor modification to the JS estimator can reduce the WMSE further. The modification is made in order to avoid the effect of an overshrinkage, $\left(G-2\right)/({{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2)>1$, by which all the signs of $\mathit{Y}$ are reverted. The overshrinkage phenomenon occurs with a small probability, and therefore, the modification is minor in the majority of cases. A modified estimator is the positive-part JS estimator

$${\mathit{\delta }}^{\mathrm{JS}+}\equiv \left({\delta }_1^{\mathrm{JS}+},\dots ,{\delta }_G^{\mathrm{JS}+}\right) \equiv {\left(1-\frac{G-2}{{{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2}\right)}^{+}\mathit{Y},$$

where ${(.)}^{+}\equiv \max \left(0,.\right)$. Consequently, ${\mathit{\delta }}^{\mathrm{JS}+}$ has a smaller WMSE than ${\mathit{\delta }}^{\mathrm{JS}}$ (p. 356, Theorem 5.4 of [25].

In summary, this subsection proposes two estimators (${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$) that improve upon the standard unbiased estimator $\mathit{Y}$.

3.2. Estimation under Ordered Means

Our next proposal is a shrinkage estimator under ordered means. We consider the case where the ordering constraints ${\mu }_1\le \dots \le {\mu }_G$ are known by the study design. Thus, the parameter $\mathit{\mu }$ is known to be restricted on the space $\left\{\left({\mu }_1,\dots ,{\mu }_G\right):{\mu }_1\le \dots \le {\mu }_G\right\}$. For instance, suppose that $i \left(=1,2,\dots ,G\right)$ represents the time index ($i=1$ for the oldest study, and $i=G$ for the newest study) at which a treatment effect ${\mu }_i$ is estimated. Then, one may assume a trend ${\mu }_1\le \dots \le {\mu }_G$ due to the improvements of treatments over time.

For instance, the true means may be $\left({\mu }_1,\dots ,{\mu }_5\right)=\left(-2,-1, 0, 1, 2\right)$. This trend could be modeled by a meta-regression with ${\mu }_i=a+bi$, where values $a$ and $b$ are unknown. In practice, one does not know any structure of the means (e.g., linear regression) except for ${\mu }_1\le \dots \le {\mu }_5$. If some knowledge, such as a linear model on covariates, is true, one could use meta-regression. However, we do not adopt any model, permitting various non-linear settings such as $\left(-2,-2, 0, 0,0\right)$ and $\left(-2,-1, 4, 4,5\right)$.

The use of the standard unbiased estimator $\mathit{Y}=\left(Y_1,\dots ,Y_G\right)$ is not desirable under the ordering constraints. Due to random variations, the estimator $\mathit{Y}=\left(Y_1,\dots ,Y_G\right)$ can be outside the parameter space, namely, $\mathit{Y}\notin \left\{\left({\mu }_1,\dots ,{\mu }_G\right):{\mu }_1\le \dots \le {\mu }_G\right\}$. Under this setting, an estimator accounting for the parameter restriction improves upon the unrestricted estimator $\mathit{Y}$, even though the former is a biased estimator [30].

The restricted maximum likelihood (RML) estimator satisfying ${\delta }_1\le \dots \le {\delta }_G$ is calculated by the pool-adjacent-violators algorithm (PAVA)

$${\delta }_i^{\mathrm{RML}}=\underset{s\le i}{\max }\underset{t\ge i}{\min }\frac{{{\displaystyle \sum }}_{j=s}^t{Y}_j}{t-s+1}.$$

This gives the RML estimator ${\mathit{\delta }}^{\mathrm{RML}}\equiv \left({\delta }_1^{\mathrm{RML}},\dots ,{\delta }_G^{\mathrm{RML}}\right)$, which has a smaller WMSE than $\mathit{Y}$. For an example of $G=3$, one has the data of $\left(Y_1,Y_2,Y_3\right)$, and the PAVA results in

$${\delta }_1^{\mathrm{RML}}=\min \left(Y_1,\frac{Y_1+Y_2}{2}, \frac{Y_1+Y_2+Y_3}{3}\right),$$

$${\delta }_2^{\mathrm{RML}}=\max \left\{\min \left(\frac{Y_1+Y_2}{2}, \frac{Y_1+Y_2+Y_3}{3}\right), \min \left( Y_2, \frac{Y_2+Y_3}{2} \right)\right\},$$

$${\delta }_3^{\mathrm{RML}}=\max \left( \frac{Y_1+Y_2+Y_3}{3}, \frac{Y_2+Y_3}{2},Y_3 \right).$$

Hence, ${\delta }_i^{\mathrm{RML}}$ is equal to $Y_i$ itself or an average including $Y_i$. Of course, ${\delta }_i^{\mathrm{RML}}=Y_i \forall i$ if $Y_1\le \dots \le Y_G$. For theories and applications of the PAVA, we refer to [31,32,33]. The max min formula written above can be found in Chapter 8 of [30] or [34].

It is clear that ${\mathit{\delta }}^{\mathrm{RML}}$ is different from order statistics $Y_{\left(1\right)}\le \dots \le Y_{\left(G\right)}$ that are a permutation of $\left(Y_1,\dots ,Y_G\right)$ and also improve the WMSE in some cases [34]. However, the permuted estimator loses the information of individual studies’ identifications, and therefore, is not considered in this article.

Below, we further improve ${\mathit{\delta }}^{\mathrm{RML}}$ with the aid of the JS estimator. Let $I(.)$ be the indicator function; $I\left(A\right)=1$ or $I\left(A\right)=0$ if $A$ is true or false, respectively. We adjust the estimator of Chang [35] who proposed the JS-type estimator under the order restriction as follows:

$${\mathit{\delta }}^{\mathrm{RJS}}\equiv \left(1-\frac{G-2}{{{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2}\right)\mathit{Y}I\left(Y_1\le \dots \le Y_G\right)+{\mathit{\delta }}^{\mathrm{RML}}\left(1-I\left(Y_1\le \dots \le Y_G\right)\right),$$

where “RJS” stands for “Restricted JS”. Note that ${\mathit{\delta }}^{\mathrm{RJS}}$ has a smaller WMSE than ${\mathit{\delta }}^{\mathrm{RML}}$, and hence, the former improves upon the latter [35].

We further improve the RJS estimator by the positive-part RJS estimator given by

$${\mathit{\delta }}^{\mathrm{RJS}+}\equiv {\left(1-\frac{G-2}{{{\displaystyle \sum }}_{i=1}^G{Y}_i^2/{\sigma }_i^2}\right)}^{+}\mathit{Y}I\left(Y_1\le \dots \le Y_G\right)+{\mathit{\delta }}^{\mathrm{RML}}\left(1-I\left(Y_1\le \dots \le Y_G\right)\right)$$

Consequently, ${\mathit{\delta }}^{\mathrm{RJS}+}$ has a smaller WMSE than ${\mathit{\delta }}^{\mathrm{RJS}}$ (Theorem 5.4 of Lehmann and Casella [25]). Note that, if $Y_1\le \dots \le Y_G$ is not satisfied, then ${\mathit{\delta }}^{\mathrm{RML}}={\mathit{\delta }}^{\mathrm{RJS}}={\mathit{\delta }}^{\mathrm{RJS}+}$.

In summary, this subsection proposes three estimators (${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, and ${\mathit{\delta }}^{\mathrm{RJS}+}$) that improve upon the standard unbiased estimator $\mathit{Y}$ under ${\mu }_1\le \dots \le {\mu }_G$.

3.3. Estimation under Sparse Means

Our third proposal is a shrinkage estimator under sparse normal means where most of the ${\mu }_i$s are zero [27,28]. The vector $\left({\mu }_1,\dots ,{\mu }_G\right)$ is called sparse if the number ${\displaystyle \sum }_{i=1}^{G}I\left({\mu }_i\ne 0\right)$ is much smaller than $G$, e.g., $\left({\mu }_1,\dots ,{\mu }_{10}\right)=\left(-5,0,0,0,0,0,0,0,0,5\right)$. In practice, one does not know which components are zeros, and how many components are zero. Nonetheless, one could assume that many of $\left({\mu }_1,\dots ,{\mu }_G\right)$ are zero. However, the elements of $\left(Y_1,\dots ,Y_G\right)$ are almost always nonzero, which disagree with the true values $\left({\mu }_1,\dots ,{\mu }_G\right)$.

Under the sparse means, it is quite reasonable to estimate ${\mu }_i$ as exactly zero if $Y_i$ is close to zero. Accordingly, one can use a thresholding estimator $Y_{i}I(\left|Y_i\right|>c_i)$ for a critical value $c_i>0$. The idea was proposed by Bancroft [36] who formulated pretest estimators that incorporate a preliminary hypothesis test into estimation. Judge and Bock [37] extensively studied pretest estimators with applications to econometrics; see also more recent works [38,39,40,41,42,43]. Among all, we particularly note that Shih et al. [41] proposed the general pretest (GPT) estimator that includes empirical Bayes and Types I-II shrinkage pretest estimators for the univariate normal mean.

We modify the GPT estimator to be adopted to meta-analyses as follows:

$${\delta }_i^{\mathrm{GPT}}=Y_{i}I\left(\left|\frac{Y_i}{{\sigma }_i}\right|>z_{{\alpha }_1}\right)+q\left(\frac{Y_i}{{\sigma }_i}\right)Y_{i}I\left(z_{{\alpha }_2}<\left|\frac{Y_i}{{\sigma }_i}\right|\le z_{{\alpha }_1}\right), i=1,2,\dots ,G.$$

for $0\le {\alpha }_1\le {\alpha }_2\le 1$, $q:\left(-\infty ,\infty \right)\mapsto \left(0,1\right)$, and $z_p$ is the upper p-th quantile of $N\left(0,1\right)$ for $0<p<1$. To implement the GPT estimator, the values of ${\alpha }_1$ and ${\alpha }_2$, and the probability function $q(.)$ must be chosen. They cannot be chosen to minimize the WMSE and TMSE criteria since pretest estimators do not permit tractable forms of MSEs [38,41]. Fortunately, for any value of ${\alpha }_1$ and ${\alpha }_2$, and a function $q$, one can show that ${\mathit{\delta }}^{\mathrm{GPT}}\equiv \left({\delta }_1^{\mathrm{GPT}},\dots ,{\delta }_G^{\mathrm{GPT}}\right)$ has smaller WMSE and TMSE values than $\mathit{Y}$ provided $\mathit{\mu }\approx \mathbf{0}$; see Appendix A.3 for the proof.

For the above reasons, we apply statistically interpretable choices of ${\alpha }_1$, ${\alpha }_2$, and $q(.)$. The special case of ${\alpha }_1={\alpha }_2=\alpha =0.025$ leads to the usual pretest estimator

$${\delta }_i^{\mathrm{PT}}=Y_{i}I\left(\left|\frac{Y_i}{{\sigma }_i}\right|>1.96\right),$$

for which $q(.)$ is arbitrary. Thus, we retain $Y_i$ if $H_0:{\mu }_i=0$ is rejected in favor of $H_1:{\mu }_i\ne 0$ at the 5% level. Otherwise, we discard $Y_i$, and conclude ${\mu }_i=0$.

For the GPT estimator, we set $q\left(z\right)=1/2$ (50% shrinkage) as suggested by Shih et al. [41]. To facilitate the interpretability of pretests, we chose ${\alpha }_1=0.025$ (5% level) and ${\alpha }_2=0.05$ (10% level). Thus, we propose the estimator

$${\delta }_i^{\mathrm{GPT}}=Y_{i}I\left(\left|\frac{Y_i}{{\sigma }_i}\right|>1.96\right)+\frac{Y_i}{2}I\left(1.645<\left|\frac{Y_i}{{\sigma }_i}\right|\le 1.96\right), i=1,2,\dots ,G.$$

Thus, if $H_0:{\mu }_i=0$ is rejected at the 5% level, we set ${\delta }_i^{\mathrm{GPT}}=Y_i$. If $H_0:{\mu }_i=0$ is accepted at the 5% level, but rejected at the 10% level, we set ${\delta }_i^{\mathrm{GPT}}=Y_i/2$. If $H_0:{\mu }_i=0$ is accepted at the 10% level, we set ${\delta }_i^{\mathrm{GPT}}=0$. Thus, ${\delta }_i^{\mathrm{GPT}}$ gives a weaker shrinkage than ${\delta }_i^{\mathrm{PT}}$ does. Obviously, we obtain a relationship $\left|{\delta }_i^{\mathrm{PT}}\right|\le \left|{\delta }_i^{\mathrm{GPT}}\right|\le \left|Y_i\right|$.

The GPT estimator introduced above is not an empirical Bayesian estimator. If ${\alpha }_1=0$, ${\alpha }_2=1$, and $q\left(y\right)=1-{\sigma }_i^2/\max \left\{{\sigma }_i^2,y^2\right\}$ were chosen, the resultant GPT estimator would be an empirical Bayes estimator [41]. However, we do not consider this estimator in our analysis.

In summary, this subsection proposes two estimators (${\mathit{\delta }}^{\mathrm{PT}}$ and ${\mathit{\delta }}^{\mathrm{GPT}}$) that improve upon the standard unbiased estimator $\mathit{Y}$ provided $\mathit{\mu }\approx \mathbf{0}$.

4. Simulation

We conducted Monte Carlo simulations to examine the performance of the proposed estimators. The purpose of the simulations is to clarify how the proposed estimators improve upon the standard unbiased estimator $\mathit{Y}$ in terms of the WMSE and TMSE, defined as

$$\mathrm{WMSE}\equiv E\left[{\displaystyle \sum }_{i=1}^G\frac{{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right],$$

$$\mathrm{TMSE}\equiv E\left[{\displaystyle \sum }_{i=1}^G{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2\right].$$

We also examined the bias defined as,

$$\mathrm{Bias}\equiv E\left[{\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right], i=1,2,\dots ,G.$$

4.1. Simulation Design

We generated normally distributed data $Y_i~N\left({\mu }_i,{\sigma }_i^2\right)$, $i=1,2,\dots ,G$, under the following scenarios:

Scenario (a):

Common means: $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_G\right)=\left(1,1,1,1,1,1,1,1,1,1\right)$,

Scenario (b):

Zero means: $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_G\right)=\left(0,0,0,0,0,0,0,0,0,0\right)$,

Scenario (c):

Ordered means: $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_G\right)=\left(-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5\right)$,

Scenario (d):

Sparse means: $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_G\right)=\left(-5,0,0,0,0,0,0,0,0,5\right)$.

For the known variances, we generated ${\sigma }_i^2~{\chi }_{\mathrm{df}=1}^2/4$, restricted in the interval $\left[0.009,0.6\right]$, as previously considered [15,44], leading to $E\left[{\sigma }_i^2\right]=0.173$. Figure 1 depicts simulated data $\left(Y_1,\dots ,Y_G\right)$ and $\left({\sigma }_1,\dots ,{\sigma }_G\right)$ under Scenarios (a) and (c). Using the data, we computed the proposed estimators ${\mathit{\delta }}^{\mathrm{JS}}$, ${\mathit{\delta }}^{\mathrm{JS}+}$, ${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, ${\mathit{\delta }}^{\mathrm{RJS}+}$, ${\mathit{\delta }}^{\mathrm{PT}}$, and ${\mathit{\delta }}^{\mathrm{GPT}}$ for estimating $\mathit{\mu }$. Our interest is to see how and when these estimators improve upon the standard unbiased estimator $\mathit{Y}$. The WMSE, TMSE and Bias of all the estimators were computed by the Monte Carlo average of 10,000 replications of such data.

4.2. Simulation Results

Table 1. Simulation results under Scenario (a): $\mu =\left({\mu }_1,\dots ,{\mu }_{10}\right)=\left(1,1,1,1,1,1,1,1,1,1\right)$.

	$\mathit{Y}$	${\mathit{\delta }}^{\mathrm{JS}}$	${\mathit{\delta }}^{\mathrm{JS}+}$	${\mathit{\delta }}^{\mathrm{RML}}$	${\mathit{\delta }}^{\mathrm{RJS}}$	${\mathit{\delta }}^{\mathrm{RJS}+}$	${\mathit{\delta }}^{\mathrm{PT}}$	${\mathit{\delta }}^{\mathrm{GPT}}$
TMSE for $\mathit{\mu }$	1.743	1.608	1.608	0.503	0.503	0.503	3.579	3.137
WMSE for $\mathit{\mu }$	10.045	9.635	9.635	5.761	5.761	5.761	17.491	15.586
Bias for ${\mu }_1$	0.01	−0.04	−0.04	−0.21	−0.21	−0.21	−0.15	−0.12
Bias for ${\mu }_2$	0.00	−0.05	−0.05	−0.12	−0.12	−0.12	−0.15	−0.12
Bias for ${\mu }_3$	0.00	−0.05	−0.05	−0.07	−0.07	−0.07	−0.16	−0.13
Bias for ${\mu }_4$	0.00	−0.05	−0.05	−0.04	−0.04	−0.04	−0.16	−0.13
Bias for ${\mu }_5$	0.00	−0.05	−0.05	−0.01	−0.01	−0.01	−0.16	−0.13
Bias for ${\mu }_6$	0.00	−0.05	−0.05	0.01	0.01	0.01	−0.16	−0.13
Bias for ${\mu }_7$	0.00	−0.05	−0.05	0.04	0.04	0.04	−0.16	−0.13
Bias for ${\mu }_8$	0.00	−0.05	−0.05	0.07	0.07	0.07	−0.16	−0.13
Bias for ${\mu }_9$	0.00	−0.05	−0.05	0.12	0.12	0.12	−0.15	−0.13
Bias for ${\mu }_{10}$	0.00	−0.05	−0.05	0.21	0.21	0.21	−0.16	−0.13

Notes:$\mathrm{TMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2\right],$ $\mathrm{WMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G\frac{{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right],$ $\mathrm{Bias}\equiv E\left[{\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right]$.

,

Table 2. Simulation results under Scenario (b): $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_{10}\right)=(0,0,0,0,0,0,0,0,0,0$).

	$\mathit{Y}$	${\mathit{\delta }}^{\mathrm{JS}}$	${\mathit{\delta }}^{\mathrm{JS}+}$	${\mathit{\delta }}^{\mathrm{RML}}$	${\mathit{\delta }}^{\mathrm{RJS}}$	${\mathit{\delta }}^{\mathrm{RJS}+}$	${\mathit{\delta }}^{\mathrm{PT}}$	${\mathit{\delta }}^{\mathrm{GPT}}$
TMSE for $\mathit{\mu }$	1.743	0.353	0.223	0.503	0.503	0.503	0.489	0.558
WMSE for $\mathit{\mu }$	10.045	2.051	1.287	5.761	5.761	5.761	2.837	3.236
Bias for ${\mu }_1$	0.01	0.00	0.00	−0.21	−0.21	−0.21	0.00	0.00
Bias for ${\mu }_2$	0.00	0.00	0.00	−0.12	−0.12	−0.12	0.00	0.00
Bias for ${\mu }_3$	0.00	0.00	0.00	−0.07	−0.07	−0.07	0.00	0.00
Bias for ${\mu }_4$	0.00	0.00	0.00	−0.04	−0.04	−0.04	0.00	0.00
Bias for ${\mu }_5$	0.00	0.00	0.00	−0.01	−0.01	−0.01	0.00	0.00
Bias for ${\mu }_6$	0.00	0.00	0.00	0.01	0.01	0.01	0.00	0.00
Bias for ${\mu }_7$	0.00	0.00	0.00	0.04	0.04	0.04	0.00	0.00
Bias for ${\mu }_8$	0.00	0.00	0.00	0.07	0.07	0.07	0.00	0.00
Bias for ${\mu }_9$	0.00	0.00	0.00	0.12	0.12	0.12	0.00	0.00
Bias for ${\mu }_{10}$	0.00	0.00	0.00	0.21	0.21	0.21	0.00	0.00

Notes:$\mathrm{TMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2\right],$ $\mathrm{WMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G\frac{{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right],$ $\mathrm{Bias}\equiv E\left[{\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right]$.

,

Table 3. Simulation results under Scenario (c): $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_{11}\right)=\left(-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5\right)$.

	$\mathit{Y}$	${\mathit{\delta }}^{\mathrm{JS}}$	${\mathit{\delta }}^{\mathrm{JS}+}$	${\mathit{\delta }}^{\mathrm{RML}}$	${\mathit{\delta }}^{\mathrm{RJS}}$	${\mathit{\delta }}^{\mathrm{RJS}+}$	${\mathit{\delta }}^{\mathrm{PT}}$	${\mathit{\delta }}^{\mathrm{GPT}}$
TMSE for $\mathit{\mu }$	1.914	1.847	1.847	1.314	1.309	1.309	2.826	2.547
WMSE for $\mathit{\mu }$	11.041	10.831	10.831	9.259	9.240	9.240	14.600	13.494
Bias for ${\mu }_1$	0.01	0.07	0.07	−0.03	−0.02	−0.02	0.02	0.02
Bias for ${\mu }_2$	0.00	0.05	0.05	0.00	0.01	0.01	0.03	0.03
Bias for ${\mu }_3$	0.00	0.03	0.03	0.00	0.00	0.00	0.08	0.06
Bias for ${\mu }_4$	0.00	0.02	0.02	0.00	0.00	0.00	0.16	0.13
Bias for ${\mu }_5$	0.00	0.01	0.01	0.00	0.00	0.00	0.20	0.17
Bias for ${\mu }_6$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
Bias for ${\mu }_7$	0.00	−0.01	−0.01	0.00	0.00	0.00	−0.20	−0.17
Bias for ${\mu }_8$	0.00	−0.02	−0.02	0.00	0.00	0.00	−0.16	−0.13
Bias for ${\mu }_9$	0.00	−0.04	−0.04	0.00	0.00	0.00	−0.08	−0.06
Bias for ${\mu }_{10}$	0.00	−0.05	−0.05	0.00	−0.01	−0.01	−0.03	−0.02
Bias for ${\mu }_{11}$	0.00	−0.06	−0.06	0.03	0.03	0.03	−0.01	−0.01

Notes:$\mathrm{TMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2\right],$ $\mathrm{WMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G\frac{{\left({\mathit{\delta }}_i\left(Y\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right],$ $\mathrm{Bias}\equiv E\left[{\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right]$.

and

Table 4. Simulation results under Scenario (d): $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_{10}\right)=\left(-5,0,0,0,0,0,0,0,0,5\right)$.

	$\mathit{Y}$	${\mathit{\delta }}^{\mathrm{JS}}$	${\mathit{\delta }}^{\mathrm{JS}+}$	${\mathit{\delta }}^{\mathrm{RML}}$	${\mathit{\delta }}^{\mathrm{RJS}}$	${\mathit{\delta }}^{\mathrm{RJS}+}$	${\mathit{\delta }}^{\mathrm{PT}}$	${\mathit{\delta }}^{\mathrm{GPT}}$
TMSE for $\mathit{\mu }$	1.743	1.719	1.719	0.823	0.823	0.823	0.739	0.794
WMSE for $\mathit{\mu }$	10.1045	9.897	9.897	7.272	7.272	7.272	4.281	4.599
Bias for ${\mu }_1$	0.01	0.10	0.10	0.01	0.01	0.01	0.01	0.01
Bias for ${\mu }_2$	0.00	0.00	0.00	−0.21	−0.21	−0.21	0.00	0.00
Bias for ${\mu }_3$	0.00	0.00	0.00	−0.11	−0.11	−0.11	0.00	0.00
Bias for ${\mu }_4$	0.00	0.00	0.00	−0.06	−0.06	−0.06	0.00	0.00
Bias for ${\mu }_5$	0.00	0.00	0.00	−0.02	−0.02	−0.02	0.00	0.00
Bias for ${\mu }_6$	0.00	0.00	0.00	−0.02	−0.02	−0.02	0.00	0.00
Bias for ${\mu }_7$	0.00	0.00	0.00	0.02	0.02	0.02	0.00	0.00
Bias for ${\mu }_8$	0.00	0.00	0.00	0.05	0.05	0.05	0.00	0.00
Bias for ${\mu }_9$	0.00	0.00	0.00	0.10	0.10	0.10	0.00	0.00
Bias for ${\mu }_{10}$	0.00	0.00	0.00	0.21	0.21	0.21	0.00	0.00
Bias for ${\mu }_{11}$	0.00	−0.09	−0.09	0.00	0.00	0.00	0.00	0.00

Notes:$\mathrm{TMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2\right],$ $\mathrm{WMSE}\equiv E\left[{{\displaystyle \sum }}_{i=1}^G\frac{{\left({\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right)}^2}{{\sigma }_i^2}\right],$ $\mathrm{Bias}\equiv E\left[{\delta }_i\left(\mathit{Y}\right)-{\mu }_i\right]$.

summarize the simulation results for comparing eight estimators: $\mathit{Y}$, ${\mathit{\delta }}^{\mathrm{JS}}$, ${\mathit{\delta }}^{\mathrm{JS}+}$,${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, ${\mathit{\delta }}^{\mathrm{RJS}+}$, ${\mathit{\delta }}^{\mathrm{PT}}$, and ${\mathit{\delta }}^{\mathrm{GPT}}$.

The JS estimators (${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$) improved upon $\mathit{Y}$ in all the scenarios (see the WMSE and TMSE in Table 1, Table 2, Table 3 and Table 4). Scenario (b) produced the greatest reduction in the values of the WMSE and TMSE, where all the true means are equal to zero (Table 2). This is exactly what the theoretical results imply (Section 3.1). The two estimators ${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$ showed similar performance except for Scenario (b), where ${\mathit{\delta }}^{\mathrm{JS}+}$ exhibited smaller values for the WMSE and TMSE than ${\mathit{\delta }}^{\mathrm{JS}}$ (Table 2). In summary, we observed a clear advantage of the two JS estimators (${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$) over the standard estimator ($\mathit{Y}$) in all the scenarios examined.

The order-restricted estimators (${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$ and ${\mathit{\delta }}^{\mathrm{RJS}+}$) were even more advantageous than ${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$ according to the smallest values of the WMSE and TMSE in Table 1 and Table 3. The substantially reduced values of the WMSE and TMSE were brought by the strict constraints by isotonic regression. We found that ${\mathit{\delta }}^{\mathrm{RJS}+}$ had the smallest value of the WMSE and TMSE under Scenarios (a) and (c) (Table 1 and Table 3). However, we observed large and systematic biases of the order-restricted estimators, which are caused by isotonic regression.

The pretest estimator ${\mathit{\delta }}^{\mathrm{PT}}$ showed the best performance under Scenario (d) (Table 4), followed by the general pretest estimator ${\mathit{\delta }}^{\mathrm{GPT}}$. However, the comparison between ${\mathit{\delta }}^{\mathrm{PT}}$ and ${\mathit{\delta }}^{\mathrm{GPT}}$ depends on the scenarios: ${\mathit{\delta }}^{\mathrm{PT}}$ was better in Scenarios (b) and (d) while ${\mathit{\delta }}^{\mathrm{GPT}}$ was better in Scenarios (a) and (c). Both ${\mathit{\delta }}^{\mathrm{PT}}$ and ${\mathit{\delta }}^{\mathrm{GPT}}$ did not perform well in Scenarios (a) and (c), where the sparsity hypothesis $H_0:{\mu }_i=0$ fails. Naturally, the sparsity of individual means is essential for the pretest estimators to perform well.

In summary, our simulations have shown that the proposed estimators have their own advantageous scenarios: ${\mathit{\delta }}^{\mathrm{JS}+}$ for Scenario (b); ${\mathit{\delta }}^{\mathrm{RJS}+}$ for Scenarios (a) and (c); ${\mathit{\delta }}^{\mathrm{PT}}$ or ${\mathit{\delta }}^{\mathrm{GPT}}$ for Scenario (d). These proposed estimators substantially reduced the values of the WMSE and TMSE over the standard estimator $\mathit{Y}$. Furthermore, the proposed estimators ${\mathit{\delta }}^{\mathrm{JS}}$, ${\mathit{\delta }}^{\mathrm{JS}+}$,${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, and ${\mathit{\delta }}^{\mathrm{RJS}+}$ improved upon $\mathit{Y}$ across all the scenarios. We therefore conclude that there are good reasons to apply the proposed estimators to estimate $\mathit{\mu }$ in order to improve the accuracy of estimation.

5. Data Analyses

We analyze two datasets to illustrate the proposed method. One dataset comes from gastric cancer patients and the other from COVID-19 pneumonia patients.

5.1. Gastric Cancer Data

We illustrate the proposed methods using a dataset from the Global Advanced/Adjuvant Stomach Tumor Research International Collaboration (GASTRIC) [45]. The data provide survival times for gastric cancer patients from 14 studies. One of the goals of this collaboration is to estimate the effects of chemotherapy on disease-free survival (DFS) for gastric cancer patients [46]. We used the data available in the R package surrosurv [47].

Let $Y_i$ be the estimator for a Cox regression coefficient for a treatment indicator (1=treatment vs. 0=control) and let ${\sigma }_i$ be the SE for $i=1,\dots ,14$.

Table 5. The 14 studies on gastric cancer patients for estimating the effects of chemotherapy on DFS.

	Sample Size	$\mathbf{Treatment} \mathbf{Effect}={\mathit{Y}}_{\mathit{i}}$	$\mathbf{SE}={\mathit{\sigma }}_{\mathit{i}}$	p-Value
Study 1	269	−0.18312	0.15372	0.234
Study 2	190	−0.72266	0.28686	0.012 **
Study 3	252	−0.48507	0.33192	0.144
Study 4	536	−0.23961	0.21558	0.266
Study 5	219	−0.13226	0.14691	0.368
Study 6	306	−0.27228	0.14416	0.059 *
Study 7	88	−0.5867	0.24885	0.018 **
Study 8	281	−0.13969	0.14542	0.337
Study 9	271	−0.1004	0.16404	0.541
Study 10	178	−0.31143	0.17038	0.068 *
Study 11	126	−0.04949	0.19818	0.803
Study 12	180	−0.11685	0.16476	0.478
Study 13	206	−0.13044	0.19268	0.498
Study 14	186	0.04391	0.17632	0.803

* 10% significance; ** 5% significance.

summarizes $Y_i$, ${\sigma }_i$, and sample sizes for the 14 studies. It is reasonable to assume the asymptotic normal approximation $Y_i~N\left({\mu }_i,{\sigma }_i^2\right)$ since the sample size of each study is large enough. Table 5 also shows the p-values for testing $H_0: {\mu }_i=0$ through Z-values $Y_i/{\sigma }_i$.

Table 5 shows that only two studies reached 5% significance. Thus, it is reasonable to assume ${\mu }_i\approx 0 \forall i$, and therefore, the JS estimators are suitable. It is also reasonable to assume ${\mu }_i=0$ for most of $i$s, and therefore, the PT and GPT estimators are suitable. We therefore calculated ${\mathit{\delta }}^{\mathrm{JS}}$, ${\mathit{\delta }}^{\mathrm{JS}+}$, ${\mathit{\delta }}^{\mathrm{PT}}$, and ${\mathit{\delta }}^{\mathrm{GPT}}$ to estimate $\mathit{\mu }$ as shown in

Table 6. Estimates of the treatment effects on DFS from the 14 studies on gastric cancer patients.

	$\mathit{Y}$	${\mathit{\delta }}^{\mathrm{JS}}$	${\mathit{\delta }}^{\mathrm{JS}+}$	${\mathit{\delta }}^{\mathrm{PT}}$	${\mathit{\delta }}^{\mathrm{GPT}}$
${\mu }_1$	−0.18	−0.10	−0.10	0	0
${\mu }_2$	−0.72	−0.40	−0.40	−0.72	−0.72
${\mu }_3$	−0.49	−0.27	−0.27	0	0
${\mu }_4$	−0.24	−0.13	−0.13	0	0
${\mu }_5$	−0.13	−0.07	−0.07	0	0
${\mu }_6$	−0.27	−0.15	−0.15	0	−0.14
${\mu }_7$	−0.59	−0.32	−0.32	−0.59	−0.59
${\mu }_8$	−0.14	−0.08	−0.08	0	0
${\mu }_9$	−0.10	−0.06	−0.06	0	0
${\mu }_{10}$	−0.31	−0.17	−0.17	0	−0.16
${\mu }_{11}$	−0.05	−0.03	−0.03	0	0
${\mu }_{12}$	−0.12	−0.06	−0.06	0	0
${\mu }_{13}$	−0.13	−0.07	−0.07	0	0
${\mu }_{14}$	0.04	0.02	0.02	0	0

All values are rounded up to 2 digits.

.

Table 6 shows that ${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$ shrink the values of $\mathit{Y}$ by 45%, namely, ${\mathit{\delta }}^{\mathrm{JS}}={\mathit{\delta }}^{\mathrm{JS}+}=\left(1-0.45\right)\times \mathit{Y}$. We also observe that ${\mathit{\delta }}^{\mathrm{PT}}$ and ${\mathit{\delta }}^{\mathrm{GPT}}$ cause 0%, 50% or 100% shrinkage of $\mathit{Y}$. The values of ${\mathit{\delta }}^{\mathrm{PT}}$ suggested the null effects (${\mu }_i=0$) for 12 studies while the values of ${\mathit{\delta }}^{\mathrm{GPT}}$ suggested the null effects for 10 studies. By shrinkage, the estimators ${\mathit{\delta }}^{\mathrm{JS}}$, ${\mathit{\delta }}^{\mathrm{JS}+}$, ${\mathit{\delta }}^{\mathrm{PT}}$, and ${\mathit{\delta }}^{\mathrm{GPT}}$ more accurately estimate $\mathit{\mu }=\left({\mu }_1,\dots ,{\mu }_{14}\right)$ than $\mathit{Y}$ do. The reason is that the data structure resembles Scenario (b), where ${\mathit{\delta }}^{\mathrm{JS}}$ and ${\mathit{\delta }}^{\mathrm{JS}+}$ performed the best (Table 2) or Scenario (d), where ${\mathit{\delta }}^{\mathrm{PT}}$ and ${\mathit{\delta }}^{\mathrm{GPT}}$ performed the best (Table 4).

5.2. COVID-19 Data

We illustrate the proposed methods using data from Pranata et al., [7] who examined the effect of hypertension on mortality for patients with COVID-19 pneumonia. One of their conclusions was that hypertension increases mortality. We obtained mortality data from Pranata et al. [7], consisting of 11 published studies (

Table 7. Summary of the 11 studies (shown by the author and publication year) on COVID-19 patients for examining the effects of hypertension on mortality.

Study	Sample Size	Male (%)	Log (RR)	SE	p-Value
Akbari 2020	440	56.4	0.6881	0.6732	0.307
Bai 2000	127	63.0	0.5933	0.2754	0.031 **
Cao 2020	102	52.0	1.1756	0.2821	0.000 **
Chen 2020	123	49.0	0.5365	0.2493	0.031 **
Chen T 2020	274	62.0	0.6780	0.1713	0.000 **
Fu 2020	200	49.5	0.5878	0.3302	0.075 *
Grasselli 2020	1591	82.0	0.4637	0.0956	0.000 **
Li 2020	102	58.0	0.5247	0.3272	0.109
Luo 2020	403	47.9	1.2326	0.1489	0.000 **
Yuan 2020	27	45.0	2.8904	1.4263	0.043 **
Zhou 2020	191	62.0	1.1378	0.2097	0.000 **

* 10% significance; ** 5% significance.

).

Let $Y_i$, $i=1,\dots ,11$, be the log of the risk ratio (RR) calculated from two-by-two contingency tables examining the association (mortality vs. hypertension). Let ${\sigma }_i$ be the SE for $i=1,\dots ,11$. Table 7 summarizes the 11 studies. Table 7 also shows the p-values for testing the no association ($H_0: {\mu }_i=0$) by Z-values $Y_i/{\sigma }_i$ assuming $Y_i~N\left({\mu }_i,{\sigma }_i^2\right)$.

Before applying the proposed methods, we changed the order of the 11 studies by the percentages of male patients (

Table 8. Estimated effects of hypertension on mortality from the 11 studies on COVID-19 patients.

Study	Male (%)	Mean	$Y$	${\mathit{\delta }}^{\mathrm{JS}}$	${\mathit{\delta }}^{\mathrm{JS}+}$	${\mathit{\delta }}^{\mathrm{RML}}$	${\mathit{\delta }}^{\mathrm{RJS}}$	${\mathit{\delta }}^{\mathrm{RJS}+}$
Grasselli 2020	82.0	${\mu }_1$	0.46	0.44	0.44	0.46	0.46	0.46
Bai 2000	63.0	${\mu }_2$	0.59	0.56	0.56	0.59	0.59	0.59
Chen T 2020	62.0 *	${\mu }_3$	0.68	0.64	0.64	0.68	0.68	0.68
Zhou 2020	62.0 *	${\mu }_4$	1.14	1.08	1.08	0.78	0.78	0.78
Li 2020	58.0	${\mu }_5$	0.52	0.50	0.50	0.78	0.78	0.78
Akbari 2020	56.4	${\mu }_6$	0.69	0.65	0.65	0.78	0.78	0.78
Cao 2020	52.0	${\mu }_7$	1.18	1.12	1.12	0.78	0.78	0.78
Fu 2020	49.5	${\mu }_8$	0.59	0.56	0.56	0.78	0.78	0.78
Chen 2020	49.0	${\mu }_9$	0.54	0.51	0.51	0.78	0.78	0.78
Luo 2020	47.9	${\mu }_{10}$	1.23	1.17	1.17	1.23	1.23	1.23
Yuan 2020	45.0	${\mu }_{11}$	2.89	2.74	2.74	2.89	2.89	2.89

All values are rounded up to 2 digits. * ties are sorted by the age (62 years-old in Chen T 2020; 56 years-old in Zhou 2020).

). The largest percentage (82%) gave the smallest RR ($Y_1=0.46$) while the smallest percentage (45%) gave the largest RR ($Y_{11}=2.89$). However, the values of $Y_i$s are not monotonically increasing ($Y_1<{\mathrm{Y}}_2<{\mathrm{Y}}_3<{\mathrm{Y}}_4>{\mathrm{Y}}_5<{\mathrm{Y}}_6<{\mathrm{Y}}_7>{\mathrm{Y}}_8>{\mathrm{Y}}_9<{\mathrm{Y}}_{10}<{\mathrm{Y}}_{11}$in Table 8). While Pranata et al. [7] regressed $Y_i$s on the percentages by meta-regression, we prefer not to assume a linear model. Alternatively, we considered the order-restricted estimators (${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, and ${\mathit{\delta }}^{\mathrm{RJS}+}$) by merely assuming that the true means are ordered by the percentages of male patients.

Table 8 shows the proposed estimators ${\mathit{\delta }}^{\mathrm{JS}}$, ${\mathit{\delta }}^{\mathrm{JS}+}$, ${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, and ${\mathit{\delta }}^{\mathrm{RJS}+}$. The JS shrinkage estimators were ${\mathit{\delta }}^{\mathrm{JS}}={\mathit{\delta }}^{\mathrm{JS}+}=\left(1-0.05\right)\times \mathit{Y}$. They reduced only 5% of $\mathit{Y}$, and hence gave very weak shrinkage toward zero. This is because ${\mu }_i$s are unlikely to be zero. The order-restricted estimator ${\mathit{\delta }}^{\mathrm{RML}}$ successfully enforced order restrictions, ${\mathit{\delta }}_1^{\mathrm{RJS}}\le \dots \le {\mathit{\delta }}_{11}^{\mathrm{RJS}}$, which was not achieved by $\mathit{Y}$. Since the order restrictions were not met for $\mathit{Y}$, we had ${\mathit{\delta }}^{\mathrm{RML}}=$ ${\mathit{\delta }}^{\mathrm{RJS}}={\mathit{\delta }}^{\mathrm{RJS}+}$ by definition. Thus, the last three estimators can clearly incorporate the knowledge that hypertension-related mortality is monotonically influenced by gender without imposing a linear model.

In summary, it would be advantageous and informative to add the order restriction to individual studies’ estimators on the association (mortality vs. hypertension). The resultant order-restricted estimators are expected to provide more accurate estimates than individual studies’ estimators. The reason is that the data structure resembles Scenario (c), where ${\mathit{\delta }}^{\mathrm{RML}}$, ${\mathit{\delta }}^{\mathrm{RJS}}$, and ${\mathit{\delta }}^{\mathrm{RJS}+}$ performed the best (Table 3).

6. Conclusions and Future Extensions

We have shown that the proposed shrinkage estimators can be more precise than the standard unbiased estimators for estimating individual means in meta-analyses. Theoretical evidence (Section 3) and simulation studies (Section 4) identified the remarkable advantages of the proposed estimators under certain scenarios of the true individual means. However, no single estimator is the best for all the scenarios. Hence, it is important to choose a suitable estimator for each scenario. We have provided two data examples to demonstrate the scenarios under which the proposed estimators are beneficial (Section 5). However, the choice of estimators depends on the goal of research, and hence, it is difficult to set a purely statistical strategy to choose the estimator. What we have concluded in this article is that there is good potential for the proposed estimators by their careful applications to real meta-analyses.

Recall that the shrinkage effects on the general pretest estimators, $\left({\delta }_1^{\mathrm{GPT}},\dots ,{\delta }_G^{\mathrm{GPT}}\right)$, are component-wise, and hence, do not involves the JS-type shrinkage as in Section 3.1 and Section 3.2. In the presence of pretests, it is unclear how the JS-type shrinkage can be introduced. This issue is an interesting topic, yet it is beyond the scope of the article.

This article focused on the shrinkage toward zero because a meta-analysis is typically a collection of individual studies of weak effects. However, other shrinkage schemes can also be considered, such as the shrinkage toward an arbitrary mean vector, a truncated space, a polyhedral convex, or a positive halfplane [12,25,48,49,50,51,52,53,54]. Different shrinkage schemes could be considered depending on the goals of meta-analyses. Therefore, there is good potential for shrinkage estimation methods to work with meta-analyses.

This article focuses on univariate response. One may consider extensions of this article to multiple (e.g., bivariate) responses. For instance, educational research may involve meta-analyses of bivariate test scores on verbal and mathematics tests [55], and mathematics and statistics tests [44]. In our data example of gastric cancer patients (Section 5.1), bivariate correlated endpoints are of great interest [55]. Correlations between responses should be employed to improve the efficiency and reduce the bias of univariate analyses [56,57,58,59,60,61,62,63,64,65] and to predict a primary outcome by the secondary outcome [66,67,68,69]. Multivariate shrinkage estimators of multivariate normal means, such as [70,71], can be considered for this extension.

One may also consider extensions of this article to non-normal (e.g., gamma distributed) responses. For non-normal distributions, shrinkage estimators exist, such as [72] for gamma distribution. Thus, introduction of shrinkage estimators to multivariate and/or non-normal models, and even their regression models, is a promising future direction.

Finally, we note that all meta-analyses are valid under the questionable assumption of no publication bias. Frequentist and Bayesian models for checking the publication biases are often considered [73,74]. It would be interesting to see if the proposed shrinkage methods can mitigate the effect of publication biases.