meta.shrinkage : An R Package for Meta-Analyses for Simultaneously Estimating Individual Means

: Meta-analysis is an indispensable tool for synthesizing statistical results obtained from individual studies. Recently, non-Bayesian estimators for individual means were proposed by applying three methods: the James–Stein (JS) shrinkage estimator, isotonic regression estimator, and pretest (PT) estimator. In order to make these methods available to users, we develop a new R package meta.shrinkage . Our package can compute seven estimators (named JS, JS+, RML, RJS, RJS+, PT, and GPT). We introduce this R package along with the usage of the R functions and the “average-min-max” steps for the pool-adjacent violators algorithm. We conduct Monte Carlo simulations to validate the proposed R package to ensure that the package can work properly in a variety of scenarios. We also analyze a data example to show the ability of the R package.

Usually, the goal of meta-analyses is to summarize individual studies to find some common effect [5,9]. The idea of estimating the common mean was originated from mathematical statistics and stratified sampling designs under fixed effect models (pp. 55-103 of [10][11][12][13]). In biostatistical methodologies, the estimation method based on random effect models [14] is popular. In either model, the goal of meta-analyses is usually to estimate the common mean by combining the estimators of individual means (Section 2.1).
In some scenarios, however, the goal of estimating the common mean is questionable. In these scenarios, meta-analyses can still be informative by looking at individual studies' means (e.g., by a forest plot). Aside from these simple meta-analyses, Bayesian posterior means provide a more sophisticated summary of individual means in a metaanalysis [15][16][17][18][19].
Recently, Taketomi et al. [20] proposed non-Bayesian estimators of individual means by applying three methods: the James-Stein (JS) shrinkage estimator, isotonic regression estimator, and pretest (PT) estimator. Their frequentist estimators were shown to be superior to the individual studies' estimators via decision theoretic criteria and Monte Carlo simulation experiments. These frequentist estimators also found their suitable applications to real data examples.
In this article, we propose a new R package meta.shrinkage in order to implement the frequentist estimators in [20]. Our package can calculate seven estimators (namely δ JS , δ JS+ , δ RML , δ RJS , δ RJS+ , δ PT and δ GPT ; see Section 3 for their definitions). We introduce this R package along with the usage of the R functions and the "average-min-max" steps for the pool-adjacent violators algorithm. We conduct Monte Carlo simulations to validate the proposed R package to ensure that the package can work properly in a variety of scenarios. We made the package freely available in the Comprehensive R Archive Network (CRAN): Available online: https://CRAN.R-project.org/package=meta.shrinkage (accessed on 14 November 2021).
This article is organized as follows. Section 2 gives the background, including a quick review of meta-analyses. Section 3 introduces the proposed R package. Section 4 conducts simulation studies to validate the proposed R package. Section 5 includes a data example to illustrate the proposed R package. Section 6 concludes the article with discussions. The appendices give the R code to reproduce the numerical results of this article.

Meta-Analysis
This subsection reviews the basic concepts for meta-analysis. To clarify the concepts, we introduce some notations and assumptions for a metaanalysis. Define G as the number of studies, where G stands for groups in the meta-analysis. For each i = 1, 2, . . . , G, let Y i be an estimator for a target estimand µ i that is unknown. Let y i be a realized value of the random variable Y i . We assume that the error is normally distributed so that Y i ∼ N µ i , σ 2 i , where σ 2 i > 0 is a known variance for the error distribution. Thus, µ i is the mean of Y i . The observed data are {y i : i = 1, 2, . . . , G} in this meta-analysis. We will not consider a setting where the variance is unknown [21,22], as this setting does not follow the framework of meta-analyses based on "summary data". Without a loss of generality, we assume that µ i = 0 corresponds to the null value.
Traditionally, the objective of meta-analyses is to estimate the common mean, denoted as µ. It is defined by the fixed effect model assumption of µ ≡ µ 1 = . . . = µ G or the random-effects model assumption of µ i ∼ N µ, τ 2 for i = 1, 2, . . . , G, where τ 2 is the between-study variance. In this article, we assume neither, since the aforementioned models do not always fit the data at hand. For instance, if the studies have ordered means (e.g., µ 1 ≤ . . . ≤ µ G ), the model is neither fixed nor random [20]. Indeed, many real meta-analyses have some covariates to systematically explain the reason for increasing means µ 1 ≤ . . . ≤ µ G or decreasing means µ 1 ≥ . . . ≥ µ G (see Section 5). In such circumstances, there is no general way to define the common mean µ. Below, we discuss what meta-analyses can do in the absence of the common mean.
The inverse variance weights in Equation (1) and (2) make it convenient to apply the classical decision theory [20]. Consequently, Taketomi et al. [20] were able to theoretically verify the (local) improvement of their estimators upon Y in terms of the WMSE. In practice, one may also be interested in the total MSE (TMSE) criterion, defined as [11,25], while the TMSE makes the theoretical analysis complex [11]. Section 4 will employ the TMSE criterion to assess the performance of all the improved estimators that will be introduced in the proposed R package.
Below, we introduce our R package, which can compute several improved estimators suggested by [20]. The goal of our package is to compute δ(y) ≡ (δ 1 (y), . . . , δ G (y)) from y that are the realized values for a random vector Y ≡ (Y 1 , . . . , Y G ).
Before embarking on details, we explain the basic variables used in the package. Let i be an estimate for µ i under a known variance σ 2 i > 0 for i = 1, 2, . . . , G. Thus, one needs to prepare {(y i , σ i ); i = 1, 2, . . . , G} to perform a meta-analysis. Accordingly, the input variables in the R console are two vectors: • y: a vector for y i s; • s: a vector for σ i s.
Below is an example for the input variables for a dataset of G = 14 in the R console.
The inverse variance weights in Equation (1) and (2) make it convenient to apply the classical decision theory [20]. Consequently, Taketomi et al. [20] were able to theoretically verify the (local) improvement of their estimators upon in terms of the WMSE. In practice, one may also be interested in the total MSE (TMSE) criterion, defined as [∑ ( ( ) − ) ] [11,25], while the TMSE makes the theoretical analysis complex [11].
Section 4 will employ the TMSE criterion to assess the performance of all the improved estimators that will be introduced in the proposed R package. Below, we introduce our R package, which can compute several improved estimators suggested by [20]. The goal of our package is to compute ( ) ≡ ( ( ), … , ( )) from that are the realized values for a random vector ≡ ( , … , ).

R Package meta.shrinkage
This section introduces our proposed R package meta.shrinkage, which can compute seven improved estimators for individual means, denoted as , , , , , , and . We will divide our explanations into four sections: Section 3.1 for and , Section 3.2 for , Section 3.3 for and , Section 3.4 for , and . Before embarking on details, we explain the basic variables used in the package. Let ~( , ) be an estimate for under a known variance > 0 for = 1,2, … , . Thus, one needs to prepare {( , ); = 1,2, … , } to perform a meta-analysis. Accordingly, the input variables in the R console are two vectors:  y: a vector for s;  s: a vector for s.
Below is an example for the input variables for a dataset of = 14 in the R console.  [20], in which gastric cancer data [28] was analyzed. The values for s are Cox regression estimates for the effect of chemotherapy on disease-free survival (DFS) for gastric cancer patients, and the values are the SEs.

James-Stein Estimator
The James-Stein (JS) estimator is defined as This dataset comes from the data analysis of [20], in which gastric cancer data [28] was analyzed. The values for y i s are Cox regression estimates for the effect of chemotherapy on disease-free survival (DFS) for gastric cancer patients, and the σ i values are the SEs.

James-Stein Estimator
The James-Stein (JS) estimator is defined as This estimator is a variant from the primitive JS estimator [29], which was derived under homogeneous variances (σ i = 1 for ∀i). The JS estimator reduces the WMSE by shrinking the vector Y toward 0. The degree of shrinkage is determined by the factor i /σ 2 i ) that typically takes its value from 0 (0% shrinkage) to 1 (100% shrinkage). In rare cases, it becomes greater than one (overshrinkage).
It was proven that δ JS has a smaller WMSE than Y when G ≥ 3 [20]. That is, Equations (1) and (2) hold. Thus, δ JS is an improve estimator without any restriction.
The positive-part JS estimator can further reduce the WMSE by avoiding the overshrinkage phenomenon of (G where (.) + ≡ max(0, .). In our R package, the function "js(.)" can compute δ JS and δ JS+ . Below is the usage in the R console. This estimator is a variant from the primitive JS estimator [29], which was derived under homogeneous variances ( = 1 for ∀ ). The JS estimator reduces the WMSE by shrinking the vector toward . The degree of shrinkage is determined by the factor ( − 2)/(∑ / ) that typically takes its value from 0 (0% shrinkage) to 1 (100% shrinkage). In rare cases, it becomes greater than one (overshrinkage).
It was proven that has a smaller WMSE than when ≥ 3 [20]. That is, Equations (1) and (2) hold. Thus, is an improve estimator without any restriction. The positive-part JS estimator can further reduce the WMSE by avoiding the overshrinkage phenomenon of ( − 2)/(∑ / ) > 1: where (. ) ≡ max (0, . ). In our R package, the function "js(.)" can compute and . Below is the usage in the R console.

Restricted Maximum Likelihood Estimators under Ordered Means
We consider the restricted maximum likelihood (RML) estimator when the individual means are ordered. Without loss of generality, we consider the increasing order ≤ ⋯ ≤ . This indicates that belongs to {( , … , ): ≤ ⋯ ≤ }. If there are such parameter constraints, they should be incorporated into the estimators in order to improve the estimation accuracy [20].
The RML estimator satisfying ≤ ⋯ ≤ is calculated by The above formula is called the pool-adjacent violators algorithm (PAVA). The computation requires the "average-min-max" steps ( Figure 1). We developed an R function "rml(.)" in our R package to perform the matrix-based computation of Figure 1. This R function initially computes all the elements of the matrix ( Figure 1). Then, it applies the command "max(apply(z, 1, min))", where "1" indicates that "min" applies to the rows of the matrix "z". This yields easy-to-understand code in the R program.

Restricted Maximum Likelihood Estimators under Ordered Means
We consider the restricted maximum likelihood (RML) estimator when the individual means are ordered. Without loss of generality, we consider the increasing order µ 1 ≤ . . . ≤ µ G . This indicates that µ belongs to {(µ 1 , . . . , µ G ) : µ 1 ≤ . . . ≤ µ G }. If there are such parameter constraints, they should be incorporated into the estimators in order to improve the estimation accuracy [20].
The RML estimator satisfying The above formula is called the pool-adjacent violators algorithm (PAVA). The computation requires the "average-min-max" steps ( Figure 1). We developed an R function "rml(.)" in our R package to perform the matrix-based computation of Figure 1. This R function initially computes all the elements of the matrix ( Figure 1). Then, it applies the command "max(apply(z, 1, min))", where "1" indicates that "min" applies to the rows of the matrix "z". This yields easy-to-understand code in the R program. To understand the reason why the matrix-based computation ( Figure 1) is necessary, we give an example for calculating from a dataset = ( , , ). By setting = 2 and = 3 in Figure 1, the calculation for proceeds as follows: The matrix in the preceding formula keeps four sub-averages of ( , , ), including . Any one of the four components of the matrix can be . Hence, the matrix is necessary as well as sufficient to calculate . That aside, matrices are easy to manipulate in R. The RML estimator ≡ ( , … , ) gives a smaller WMSE than [20]. For theories and applications of the PAVA, we refer to [30][31][32][33][34]. We decided not to use the "pava(.)" function available in the R package "Iso" [33] to ensure the independence of our package from others. Nonetheless, we checked that "rml(.)" and "pava(.)" gave numerically identical results.
So far, we have assumed that the studies (i.e., = 1,2, … , ) are ordered to achieve ≤ ⋯ ≤ . However, usually, the studies in a raw dataset may be arbitrarily ordered, and hence, one needs to find covariates to order the studies. For instance, one can use publication years if the values increase with them. More generally, we assume that there exists an increasing sequence of covariates ( < ⋯ < ) or a decreasing sequence of covariates ( > ⋯ > ) to achieve the order ≤ ⋯ ≤ .
In our R package, the function "rml(.)" can compute . The function allows users to enter covariates when studies are not ordered. For instance, we enter the estimates (y), SEs (s), and the proportion of males (x) from the COVID-19 data with = 11 [7,20].
In this data, the estimates (y) are the log risk ratios (RRs) calculated from two-by-two contingency tables examining the association (mortality vs. hypertension). As found in the previous studies [7,20], there was a decreasing sequence of the proportion of males ( > ⋯ > ) that could achieve the order ≤ ⋯ ≤ . Then, we have the following: To understand the reason why the matrix-based computation ( Figure 1) is necessary, we give an example for calculating δ RML 2 from a dataset Y = (Y 1 , Y 2 , Y 3 ). By setting i = 2 and G = 3 in Figure 1, the calculation for δ RML 2 proceeds as follows: The matrix in the preceding formula keeps four sub-averages of (Y 1 , Y 2 , Y 3 ), including Y 2 . Any one of the four components of the matrix can be δ RML 2 . Hence, the matrix is necessary as well as sufficient to calculate δ RML 2 . That aside, matrices are easy to manipulate in R.
So far, we have assumed that the studies (i.e., i = 1, 2, . . . , G) are ordered to achieve µ 1 ≤ . . . ≤ µ G . However, usually, the studies in a raw dataset may be arbitrarily ordered, and hence, one needs to find covariates to order the studies. For instance, one can use publication years if the µ i values increase with them. More generally, we assume that there exists an increasing sequence of covariates (x 1 < . . . < x G ) or a decreasing sequence of covariates ( In our R package, the function "rml(.)" can compute δ RML . The function allows users to enter covariates when studies are not ordered. For instance, we enter the estimates (y), SEs (s), and the proportion of males (x) from the COVID-19 data with G = 11 [7,20]. To understand the reason why the matrix-based computation ( Figure 1) is necessary, we give an example for calculating from a dataset = ( , , ). By setting = 2 and = 3 in Figure 1, the calculation for proceeds as follows: The matrix in the preceding formula keeps four sub-averages of ( , , ), including . Any one of the four components of the matrix can be . Hence, the matrix is necessary as well as sufficient to calculate . That aside, matrices are easy to manipulate in R. The RML estimator ≡ ( , … , ) gives a smaller WMSE than [20]. For theories and applications of the PAVA, we refer to [30][31][32][33][34]. We decided not to use the "pava(.)" function available in the R package "Iso" [33] to ensure the independence of our package from others. Nonetheless, we checked that "rml(.)" and "pava(.)" gave numerically identical results.
So far, we have assumed that the studies (i.e., = 1,2, … , ) are ordered to achieve ≤ ⋯ ≤ . However, usually, the studies in a raw dataset may be arbitrarily ordered, and hence, one needs to find covariates to order the studies. For instance, one can use publication years if the values increase with them. More generally, we assume that there exists an increasing sequence of covariates ( < ⋯ < ) or a decreasing sequence of covariates ( > ⋯ > ) to achieve the order ≤ ⋯ ≤ .
In our R package, the function "rml(.)" can compute . The function allows users to enter covariates when studies are not ordered. For instance, we enter the estimates (y), SEs (s), and the proportion of males (x) from the COVID-19 data with = 11 [7,20]. - In this data, the estimates (y) are the log risk ratios (RRs) calculated from two-by-two contingency tables examining the association (mortality vs. hypertension). As found in the previous studies [7,20], there was a decreasing sequence of the proportion of males ( > ⋯ > ) that could achieve the order ≤ ⋯ ≤ . Then, we have the following: In this data, the estimates (y) are the log risk ratios (RRs) calculated from two-by-two contingency tables examining the association (mortality vs. hypertension). As found in the previous studies [7,20], there was a decreasing sequence of the proportion of males (x 1 > . . . > x 11 ) that could achieve the order µ 1 ≤ . . . ≤ µ 11 . Then, we have the following: . By the option "test=TRUE", one can test if the values are properly ordered by a sequence of covariates. Figure 2 shows the output including the LOWESS plot and a correlation test based on Kendall's tau via "cor.test(x,y,method="kendall")". The test confirmed that the means were ordered by a "decreasing" sequence ( > ⋯ > ). We suggest the 10% significance level to declare the increasing or decreasing trend. In meta-analyses, 5% is too strict and not realistic since the number of studies is limited.

Shrinkage Estimators under Ordered Means
The estimator introduced above (in Section 3.2) can be improved by the JS type shrinkage. Based on the idea of Chang [35], Taketomi et al. [20] proposed a JS-type estimator under the following order restriction: where "RJS" stands for "restricted JS", (. ) is the indicator function, where ( ) = 1 or ( ) = 0 if is true or false, respectively. Note that has a smaller WMSE than [20,35], meaning that gives more precise estimates than .
We see that the estimates are ordered, which is consistent with the prescribed order µ 1 ≤ . . . ≤ µ 11 . By the option "test=TRUE", one can test if the µ i values are properly ordered by a sequence of covariates. Figure 2 shows the output including the LOWESS plot and a correlation test based on Kendall's tau via "cor.test(x,y,method="kendall")". The test confirmed that the means were ordered by a "decreasing" sequence (x 1 > . . . > x 11 ). We suggest the 10% significance level to declare the increasing or decreasing trend. In meta-analyses, 5% is too strict and not realistic since the number of studies is limited. . By the option "test=TRUE", one can test if the values are properly ordered by a sequence of covariates. Figure 2 shows the output including the LOWESS plot and a correlation test based on Kendall's tau via "cor.test(x,y,method="kendall")". The test confirmed that the means were ordered by a "decreasing" sequence ( > ⋯ > ). We suggest the 10% significance level to declare the increasing or decreasing trend. In meta-analyses, 5% is too strict and not realistic since the number of studies is limited.

Shrinkage Estimators under Ordered Means
The estimator introduced above (in Section 3.2) can be improved by the JS type shrinkage. Based on the idea of Chang [35], Taketomi et al. [20] proposed a JS-type estimator under the following order restriction: where "RJS" stands for "restricted JS", (. ) is the indicator function, where ( ) = 1 or ( ) = 0 if is true or false, respectively. Note that has a smaller WMSE than [20,35], meaning that gives more precise estimates than .

Shrinkage Estimators under Ordered Means
The estimator introduced above (in Section 3.2) can be improved by the JS type shrinkage. Based on the idea of Chang [35], Taketomi et al. [20] proposed a JS-type estimator under the following order restriction: where "RJS" stands for "restricted JS", I(.) is the indicator function, where I(A) = 1 or I(A) = 0 if A is true or false, respectively. Note that δ RJS has a smaller WMSE than δ RML [20,35], meaning that δ RJS gives more precise estimates than δ RML . The RJS estimator can be further corrected by the following positive-part RJS estimator: Consequently, δ RJS+ has a smaller WMSE than δ RJS [20]. In our R package, the function "rjs(.)" can compute δ RJS and δ RJS+ possibly with the aid of covariates.
For instance, one can analyze the COVID-19 data as follows: Algorithms 2022, 10, x FOR PEER REVIEW 7 of 17 The RJS estimator can be further corrected by the following positive-part RJS estimator: Consequently, has a smaller WMSE than [20]. In our R package, the function "rjs(.)" can compute and possibly with the aid of covariates.
For instance, one can analyze the COVID-19 data as follows: - In the above commands, the input includes the optional argument "id" that signifies the leading authors and publication years of the 11 studies. We did so simply to make an informative output. If the input does not include this argument, the output shows the ordered sequence of 1, 2, …, 11 as shown in Section 3.2.

Estimators under Sparse Means
We now consider discrete shrinkage schemes by pre-testing : = 0 vs. : ≠ 0 for = 1,2, … , . The idea was proposed by Bancroft [36], who developed pretest estimators (see also more recent works [20,[37][38][39][40][41][42][43]). In the meta-analytic context, Taketomi et al. [20] adopted the general pretest (GPT) estimator of Shih et al. [41], which is defined as follows: Here 0 ≤ ≤ ≤ 1, 0 < < 1, and is the upper pth quantile of (0,1) for 0 < < 1. To implement the GPT estimator, the values of , , and must be chosen. For any value of and , as well as a function , one can show that ≡ ( , … , ) has smaller WMSE and TMSE values than , provided ≈ [20]. One may choose = 1/2 (50% shrinkage), as suggested by [20,41]. To facilitate the interpretability of the pretests, one may choose = 0.05 (5% level) and The special case of = = = 0.05 leads to the usual pretest (PT) estimator In the above commands, the input includes the optional argument "id" that signifies the leading authors and publication years of the 11 studies. We did so simply to make an informative output. If the input does not include this argument, the output shows the ordered sequence of 1, 2, . . . , 11 as shown in Section 3.2.

Simulation: Validating the R Package
This technical section is devoted to the numerical verification of the proposed package via Monte Carlo simulation experiments. Users of the package may skip this section.
We conducted simulations to investigate the operating performance of the seven estimators implemented in the proposed R package (Section 3). Our simulation design added new scenarios to the original ones in [20]. Hence, the simulation not only added new knowledge on the performance of the seven estimators, but it also validated the proposed R package. The output shows that δ PT and δ GPT lead to 0%, 50%, or 100% shrinkage of Y. The estimates of δ PT yield 0 for 12 studies, and the estimates of δ GPT yield 0 for 10 studies.

Simulation: Validating the R Package
This technical section is devoted to the numerical verification of the proposed package via Monte Carlo simulation experiments. Users of the package may skip this section.
We conducted simulations to investigate the operating performance of the seven estimators implemented in the proposed R package (Section 3). Our simulation design added new scenarios to the original ones in [20]. Hence, the simulation not only added new knowledge on the performance of the seven estimators, but it also validated the proposed R package.
Our simulations were based on 10,000 repetitions using Y (r) , where 1 ≤ r ≤ 10, 000. Let δ Y (r) = δ 1 Y (r) , . . . , δ G Y (r) be one of the seven estimators in the rth repetition. To assess the TMSE, we computed it via the Monte Carlo average: As the TMSE and WMSE gave the same conclusion, we reported on the former. Appendix A provides the R code for the simulations, which can reproduce the results of the following section. Figure 3 compares the estimators Y, δ JS , δ JS+ , δ RML , δ RJS , δ RJS+ , δ PT , and δ GPT . In Scenarios (a) and (e), the smallest TMSE values were achieved by δ RML , δ RJS , and δ RJS+ , which appropriately accounted for the ordered means. Thus, these ordered mean estimators provided some advantages over the standard estimator Y. Here, users needed to specify the option "decreasing = FALSE" (Scenario (a)) or "decreasing = TRUE" (Scenario (e)) to capture the true ordering of the means. On the other hand, δ PT and δ GPT produced unreasonably large TMSE values, since they wrongly imposed the sparse mean assumptions.

Simulation Results
In Scenario (b), the smallest TMSE values were attained by δ PT followed by δ GPT , as they took advantage of the sparse means. The TMSE values for δ RML , δ RJS , and δ RJS+ were also small by accounting for the ordered means. Hence, these pretest and restricted estimators produced significant advantages over the standard estimator Y.
In Scenario (c), δ JS and δ JS+ performed the best, but the advantage over Y was modest. In this scenario, δ RML , δ RJS , and δ RJS+ produced quite large TMSE values and performed the worst, since they wrongly assumed the ordered means. Additionally, δ PT and δ GPT produced large TMSE values since they wrongly assumed the sparse means. This was the only scenario where the standard estimator Y was enough.
In Scenario (d), the smallest TMSE values were attained by δ GPT , as it captured the sparse means. The performance of δ PT , δ JS , and δ JS+ was also good. On the other hand, δ RML , δ RJS , and δ RJS+ gave large TMSE values due to the unordered means.
In summary, our simulations demonstrated that the seven estimators implemented in the proposed R package exhibited desired operating characteristics. If the true means were ordered, the restricted estimators (δ RML , δ RJS , and δ RJS+ ) showed definite advantages over the standard estimator Y. In addition, the pretest estimators (δ PT and δ GPT ) produced the best performance under the sparse means. Finally, the JS estimators (δ JS and δ JS+ ) modestly but uniformly improved upon Y across all the scenarios.
We therefore conclud that there are good reasons to apply the proposed R package to estimate µ in order to improve the accuracy of estimation.

Data Example
This section analyzes a dataset to illustrate the methods in the proposed package, demonstrating their possible advantages over the standard meta-analysis. Appendix B provides the R code, which can reproduce the following results.
We used the blood pressure dataset containing = 10 studies, where each study examined the effect of a treatment to reduce blood pressure. The dataset is available in the R package mvmeta. Available online: https://CRAN.R-project.org/package=mvmeta (accessed on 14 January 2022). Each study provided the treatment's effect estimate on the systolic blood pressure (SBP) and the treatment's effect on the diastolic blood pressure (DBP), as shown in Table 1. In the following analysis, we focus on the treatment's effect estimates for the SBP and regard those for the DBP as covariates. Table 1. The 10 studies from the blood pressure data. Each study provided the treatment's effect on the systolic blood pressure (SBP) and the treatment's effect on the diastolic blood pressure (DBP). We aimed to improve the individual treatment effects on the SBP by the methods in the R package meta.shrinkage. For this purpose, we utilized the covariate information to implement meta-analyses under ordered means (Sections 3.2 and 3.3).

Treatment Effect on SBP
To apply the proposed methods, we changed the order of the 10 studies by the increasing order of the covariates (see the first column of Table 2). Under this order, we defined , where = 1, … ,10 as the treatment effect estimates on the SBP (see in the second column of Table 2). Table 2 shows a good concordance between and the covariates; the smallest covariate (−7.87) yielded the smallest outcome ( = −17.93), and the largest covariate (−2.08) yielded the largest outcome ( = −6.55). However, the values of were not perfectly ordered.

Data Example
This section analyzes a dataset to illustrate the methods in the proposed package, demonstrating their possible advantages over the standard meta-analysis. Appendix B provides the R code, which can reproduce the following results.
We used the blood pressure dataset containing G = 10 studies, where each study examined the effect of a treatment to reduce blood pressure. The dataset is available in the R package mvmeta. Available online: https://CRAN.R-project.org/package=mvmeta (accessed on 14 November 2021). Each study provided the treatment's effect estimate on the systolic blood pressure (SBP) and the treatment's effect on the diastolic blood pressure (DBP), as shown in Table 1. In the following analysis, we focus on the treatment's effect estimates for the SBP and regard those for the DBP as covariates. We aimed to improve the individual treatment effects on the SBP by the methods in the R package meta.shrinkage. For this purpose, we utilized the covariate information to implement meta-analyses under ordered means (Sections 3.2 and 3.3).
To apply the proposed methods, we changed the order of the 10 studies by the increasing order of the covariates (see the first column of Table 2). Under this order, we defined Y i , where i = 1, . . . , 10 as the treatment effect estimates on the SBP (see Y in the second column of Table 2). Table 2 shows a good concordance between Y and the covariates; the smallest covariate (−7.87) yielded the smallest outcome (Y 1 = −17.93), and the largest covariate (−2.08) yielded the largest outcome (Y 11 = −6.55). However, the values of Y i were not perfectly ordered. We therefore considered the order-restricted estimators (δ RML , δ RJS , and δ RJS+ ) by imposing the assumption that the true treatment effects were ordered. Table 2 shows that these restricted estimators satisfied δ 1 ≤ . . . ≤ δ 10 (see the fifth through the seventh columns of Table 2). Since it was reasonable to impose a concordance between the treatment effects on the SBP and DBP, these estimators may have been advantageous over the standard estimates (Y). In this data example, the JS estimators (δ JS and δ JS+ ) were almost identical to the standard estimates Y. Hence, there would be little advantage to shrinkage in the JS estimators.
By the "rml(.)" function, we tested if the µ i values were ordered by a sequence of covariates. Figure 4 shows a correlation-based test and the LOWESS plot. The test confirmed that the means were ordered by an increasing sequence (x 1 < . . . < x 10 ). The P-value was highly significant (p = 0.002). Therefore, we validated the assumption of ordered means.  We therefore considered the order-restricted estimators ( , , and ) by imposing the assumption that the true treatment effects were ordered. Table 2 shows that these restricted estimators satisfied ≤ ⋯ ≤ (see the fifth through the seventh columns of Table 2). Since it was reasonable to impose a concordance between the treatment effects on the SBP and DBP, these estimators may have been advantageous over the standard estimates ( ). In this data example, the JS estimators ( and ) were almost identical to the standard estimates . Hence, there would be little advantage to shrinkage in the JS estimators.
By the "rml(.)" function, we tested if the values were ordered by a sequence of covariates. Figure 4 shows a correlation-based test and the LOWESS plot. The test confirmed that the means were ordered by an increasing sequence ( < ⋯ < ). The P-value was highly significant (p = 0.002). Therefore, we validated the assumption of ordered means.

Conclusions and Future Extensions
This article introduced an R package meta.shrinkage (https://CRAN.R-project.org/package=meta.shrinkage (accessed on 15 November 2021)), which we made freely available on CRAN. It was first released on 19 November 2021 (version 0.1.0), following our original methodological article published on 20 October 2021 [20]. We hope that the timely release of our package facilitates the appropriate use of the proposed methods for interested readers. As the precision and reliability of the developed statistical methods are important, we conducted extensive simulation studies to validate the proposed R package

Conclusions and Future Extensions
This article introduced an R package meta.shrinkage (https://CRAN.R-project.org/ package=meta.shrinkage (accessed on 15 November 2021)), which we made freely available on CRAN. It was first released on 19 November 2021 (version 0.1.0), following our original methodological article published on 20 October 2021 [20]. We hope that the timely release of our package facilitates the appropriate use of the proposed methods for interested readers. As the precision and reliability of the developed statistical methods are important, we conducted extensive simulation studies to validate the proposed R package (Section 4). We also analyzed a data example to show the ability of the R package (Section 5).
To implement isotonic regression in our R package, we proposed a matrix-based algorithm for the PAVA (Figure 1). This algorithm is easy to program in the R environment, where matrices are convenient to manipulate. However, if one tries to implement the PAVA in other programing environments, the matrix-based algorithm may be inefficient, especially for a meta-analysis with a very large number of studies. The number G > 500 makes the proposed algorithm slow to compute the output. However, real meta-analyses rarely have more than 100 studies, so this issue may not arise at the practical level.
The proposed R package can only handle normally distributed data; it cannot handle data that are non-normally distributed, asymmetrically distributed, or discrete-valued. To analyze such data, one should consider extensions of the meta-analysis methods in [20] toward asymmetric distributions for skewed or discrete response variables. Shrinkage and pretest estimators exist, such as [61] for the exponential distribution, [62] for the gamma distribution, and [63,64] for the Poisson distribution. Thus, meta-analytical applications of shrinkage estimators to asymmetric or non-normal models are relevant research directions.