Permutation inference methods for multivariate meta-analysis

Multivariate meta-analysis is gaining prominence in evidence synthesis research because it enables simultaneous synthesis of multiple correlated outcome data, and random-effects models have generally been used for addressing between-studies heterogeneities. However, coverage probabilities of confidence regions or intervals for standard inference methods for random-effects models (e.g., restricted maximum likelihood estimation) cannot retain their nominal confidence levels in general, especially when the number of synthesized studies is small because their validities depend on large sample approximations. In this article, we provide permutation-based inference methods that enable exact joint inferences for average outcome measures without large sample approximations. We also provide accurate marginal inference methods under general settings of multivariate meta-analyses. We propose effective approaches for permutation inferences using optimal weighting based on the efficient score statistic. The effectiveness of the proposed methods is illustrated via applications to bivariate meta-analyses of diagnostic accuracy studies for airway eosinophilia in asthma and a network meta-analysis for antihypertensive drugs on incident diabetes, as well as through simulation experiments. In numerical evaluations performed via simulations, our methods generally provided accurate confidence regions or intervals under a broad range of settings, whereas the current standard inference methods exhibited serious undercoverage properties.


Introduction
Multivariate meta-analysis is gaining prominence in evidence synthesis research in clinical epidemiology and health technology assessments because it enables simultaneous synthesis of multiple correlated outcome data, and can thus borrow strength in statistical inference from different correlated outcomes (Jackson, Riley and White, 2011). Multivariate metaanalysis methods have been applied to various types of meta-analyses in specific clinical contexts, e.g., meta-analysis for diagnostic test accuracy (Deeks, 2001), network metaanalysis (Caldwell, Ades and Higgins, 2005), and individual participant data meta-analysis (Riley, Lambert and Abo-Zaid, 2010).
In multivariate meta-analyses, heterogeneity in effect sizes from different studies commonly arises, and random-effects models are generally adopted to account for such heterogeneity (Jackson et al., 2011). However, standard inference methods for these randomeffects models depend on large sample approximations for the number of trials synthesized, e.g., the extended DerSimonian-Laird (EDL) methods (Chen, Manning and Dupuis, 2012;Jackson, White and Riley, 2013;Jackson, White and Thompson, 2010) and restricted maximum likelihood (REML) estimation (Jackson et al., 2011). However, the numbers of trials are often small for the large sample approximations to hold. Consequently, inference methods are not valid, i.e., the coverage probabilities of the confidence regions or intervals cannot retain their nominal confidence levels. Moreover, the type-I error probabilities of the corresponding tests cannot be retained, as shown in the simulation studies in Section 5. The same problem with random-effects models is widely recognized in the context of conventional univariate meta-analysis, even when the models are completely specified (Noma, 2011;Veroniki et al., 2019). This invalidity issue may seriously influence the overall conclusions of these evidence synthesis studies. Although several refined methods have been proposed for overcoming this issue (Jackson and Riley, 2014;Noma et al., 2018), they are still based on large sample approximations for the number of studies synthesized.
In this study, we developed alternative effective inference methods that use permutation approaches to resolve the invalidity problem. Exact permutation-based inference methods have been extensively studied in conventional univariate meta-analysis (Follmann and Proschan, 1999;Liu et al., 2018), and simulation-based evidence has demonstrated that these methods perform well even in the context of small studies. First, we provide generalized methods for applications of these exact permutation methods to multivariate meta-analysis for joint inferences of average outcome measure parameters of random effects models.
Through permutation approaches, exact tests and confidence regions of these parameters can be constructed; moreover, these methods enable accurate inferences without large sample approximations. However, in marginal inferences for individual components of the average treatment effect parameters, other components formally become nuisance parameters. Exact permutation methods cannot be formally constructed under such circumstances. Therefore, we propose an alternative effective approach using the local Monte Carlo method . The proposed marginal inference methods are not exact, and are based on asymptotic approximations, but in the simulation-based evaluations in Section 5, we show that they generally retain validity and accuracy under a wide range of settings for multivariate meta-analyses. We also demonstrate the effectiveness of these methods via real data applications to a meta-analysis of diagnostic accuracy studies for airway eosinophilia in asthma (Korevaar et al., 2015) and a network meta-analysis for antihypertensive drugs on incident diabetes (Elliott and Meyer, 2007). The numbers of synthesized studies for these meta-analyses were 12 and 22 respectively, and the large sample approximations may have been violated. Because underestimation of statistical errors might influence the overall conclusions of these studies, we re-evaluated them using the methods proposed in Section 6.

Multivariate random effects model for multivariate meta-analysis
We consider the general random-effects model for multivariate meta-analysis to address between-studies heterogeneity for multiple outcomes, which is a generalization of the DerSimonian-Laird-type random-effects model for conventional univariate meta-analyses (DerSimonian and Laird, 1986). Let denote an estimator of the rth outcome measure in the th study ( 1, 2, … , ; 1,2, … , ). Commonly used effect measures include mean difference, standardized mean difference, risk difference, risk ratio, odds ratio, and hazard ratio. Typically, the ratio measures are log-transformed to allow approximations based on normal distributions. Further, in meta-analyses for diagnostic test accuracy, the logittransformed sensitivity and specificity (or false positive rate) are typically used (Reitsma et al., 2005). Here we consider the multivariate random-effects model for the outcome vector , , … , and its within-study variance-covariance matrix (a matrix), which is assumed to be known and substituted by its valid estimate (Jackson et al., 2011), covariance parameters can be assumed to be equal or different across treatments, and the model of covariance structure is appropriately selected for applications (Jackson et al., 2011).

Exact joint inference methods
For conventional univariate meta-analyses, Follmann and Proschan (1999) and Liu et al.
(2018) developed permutation methods for exact inference of the random-effects model. In synthesizing randomized clinical trials, their idea is based on the re-randomization argument considering that the active treatment and control groups within a trial form a "pair." Further, their idea can be more generally interpreted as an adaptation of the one-sample permutation test based on the symmetric assumptions of test statistics around the null value (Good, 2000).
First, we developed an extended exact permutation method for joint inference of the grand mean vector of the multivariate random-effects model (*). We consider the null hypothesis H0: , and construct an exact statistical test. Follmann and Proschan (1999) and Liu et al. (2018) directly used the grand mean estimator of the univariate random effects model as the test statistic, but in the multivariate random-effects model, the corresponding estimator is vector-valued, and therefore cannot be straightforwardly adapted to the testing problem. Alternatively, we consider appropriate test statistics for the joint null hypothesis. An effective and analytically tractable statistic is the efficient score test (Cox and Hinkley, 1974 The two-sided p-value is calculated by comparing quantiles of ; and . In addition, the corresponding 100 1 % confidence region of can be constructed by a set of that satisfies

;
The confidence limits cannot be expressed in closed form, but can be computed by numerical techniques, e.g., the bisectional method (Burden and Faires, 2010). Further, the confidence regions can be graphically presented on a multidimensional space by plotting the null values that fulfill the above criterion, as shown in Section 6.
Although permutation inference based on the efficient score statistic is an efficient method when using the optimal weights, computations of the CML estimate of for each permutation require iterative calculations (e.g., the Newton-Raphson method) resulting in large overall computational burdens. To circumvent this computational problem, Liu et al.
(2018) proposed using another sign-invariant method-of-moments estimator in the univariate setting. We can provide a multivariate extended estimator,

1
For a study that reports a subset of outcomes, , , and are reduced to the corresponding subvectors and submatrices. Under the permutation scheme, the sign of is independent of the magnitudes of its individual components, and thus is invariant for all combinations of signs of in the permutation inferences. Note that the diagonal elements of with non-positive values should be replaced with 0, as in a conventional DerSimonian-Laird-type estimator (DerSimonian and Laird, 1986). The corresponding non-diagonal elements should also be replaced with 0, which is naturally induced by the fundamental property of covariance. In addition, in that case, is semi-positive definite, and this truncation is equivalent to performing an eigendecomposition on the estimated matrix and setting the negative eigenvalues to 0. Plugging in the methodof-moments estimator into the efficient score statistic instead of the CML estimate, we can construct a "pseudo"-efficient score statistic that is computationally efficient in the permutation inferences,

, , ,
Note that cannot be defined for an incomplete dataset in general, and the validity of inferences might be violated because it is not founded on likelihood-based methods.
However, in certain applications of multivariate meta-analyses, complete outcomes are usually available, e.g., in meta-analyses for diagnostic accuracy studies; sensitivity and specificity are commonly available for all studies, as shown in Section 6.1. Accordingly, this is a computationally efficient alternative for the permutation inference based on the efficient score statistic.

Marginal inference methods
Permutation methods can also be applied to marginal inferences for individual components of , i.e., to construct confidence intervals for individual components of . The permutation algorithm is concretely constructed as follows: Algorithm 2 (Permutation test for marginal inference of using the efficient score statistic).
1. For the multivariate random-effects model (*), compute the efficient score statistic for the test of the composite null hypothesis H0: , where , , . Note that the test statistic , corresponds to the ordinary efficient score statistic for testing H0:

Obtain the permutation null distribution
, . Therefore, the corresponding test is the most powerful test asymptotically. A note for the validity conditions of this inference method is provided in e-

Appendix-B.
The two-sided p-value is also calculated by comparing a quantile of ; , and , . In addition, the corresponding 100 1 % confidence intervals of can be constructed by a set of , that satisfies , 1 ; , The confidence limits cannot be expressed in closed form generally, and therefore must be computed by numerical methods, e.g., the bisectional method (Burden and Faires, 2010).

Simulations
We conducted simulation studies to evaluate the performances of the proposed permutation methods. The simulation settings mimicked the bivariate meta-analyses for diagnostic accuracy studies in Section 6 (Korevaar et al., 2015). We generated two binomial variables that correspond to the numbers of true positives and true negatives in the th study, , and the within-study correlation was 0. The random effects logit and logit were generated from a bivariate normal distribution of (*).
We set the grand mean parameters as  Under all of the scenarios investigated, the coverage probabilities of the EDL, ML, and REML methods were largely below the nominal level (95%), and seriously underestimated the statistical uncertainties, especially for joint inferences. In general, the EDL and REML methods exhibited favorable performances relative to the ML method, as expected, but the coverage rates were in general less than 0.95. In marginal inferences, the undercoverage properties were less extreme, but the coverage probabilities were still generally lower than 0.95. In particular, under small N and large heterogeneity settings, undercoverage properties were especially serious.
Note that for the permutation methods, the coverage probabilities were generally larger Additional simulation studies considering other various settings (bivariate and trivariate settings, and those involving missing outcomes) are presented e-Appendices D, E, and F in Supporting Information.

Meta-analysis of diagnostic accuracy of airway eosinophilia in asthma
To illustrate our method, we analyzed datasets from a meta-analysis of diagnostic accuracy of airway eosinophilia in asthma, performed by Korevaar et al. (2015). Although eosinophilic airway inflammation is associated with elevated corticosteroid responsiveness in asthma, direct airway sampling methods are invasive or laborious. Therefore, minimally invasive markers for diagnosis were investigated. Korevaar et al. (2015) conducted meta-analyses for diagnostic studies of several markers using multivariate meta-analysis methods. In particular, in their meta-analyses of the fraction of exhaled nitric oxide (FeNO) and blood eosinophils, they conducted bivariate random effects meta-analyses for sensitivity and specificity of 12 studies (total participants: 1720 and 1967).
We analyzed these datasets using the permutation methods. We fitted the bivariate random effects model (*) to the logit-transformed sensitivity and FPR of individual studies.
The diagnostic studies are not randomized trials, but the validities of the proposed permutation methods hold by the symmetry assumption concerning the null values as in conventional one-sample permutation tests; thus, they can be applied to these situations based on the symmetric assumptions. The concrete procedure to adapt the permutation methods to binomial variables for the results of individual diagnostic studies is the same as that in Section 5. The one-sample permutations assume that the summary statistics (sensitivity and specificity) distribute around a null vector (a fixed point corresponding to the null hypothesis) in accordance with the algorithms provided in Sections 3 and 4.
In Figure 2, we show the 95% confidence regions for joint inferences of sensitivity and FPR, obtained by the two proposed permutation methods. We exhaustively conducted joint tests on the two-dimensional space of sensitivity and FPR, and the non-significant null values by two-sided 5% significant level tests are depicted by gray points in this space. Estimates of sensitivity and FPR of individual studies varied widely on the two-dimensional plots, and there were substantial heterogeneities among the studies (as noted later). The number of permutations was consistently set to 2400. Although the point estimates of the proposed inference methods were not formally defined, natural choices would be their medianunbiased estimates, and we present numerical median-unbiased estimates in Figure 2. The dashed line represents the 95% confidence regions provided by the REML-based Wald-type method. For both the FeNO and blood eosinophils datasets, the confidence regions of the two permutation-based methods were wider than those of the REML method. These results may reflect the undercoverage properties of the REML method as shown in the simulation studies.
Besides, although the two permutation methods provided exact confidence regions, the obtained regions were slightly different, and the T1 statistic provided narrower confidence regions. These trends were possibly due to the efficiency of the efficient score statistic. Note that the permutation methods did not provide symmetric confidence regions, as they were nonparametric approaches for obtaining the reference distributions of the corresponding test statistics.
In Table 1, we present 95% confidence intervals based on marginal inferences of sensitivity and FPR. We present the results of EDL (Chen et al., 2012;Jackson et al., 2013), ML, REML and permutation methods. REML estimates of the between-studies standard deviation of logit-transformed sensitivity and FPR were 0.558 and 0.687, respectively, for the FeNO dataset, and 0.298 and 0.455, respectively, for the blood eosinophils dataset; accordingly, substantial heterogeneities were observed. The number of permutations was also set consistently to 2400. The proposed permutation method provided wider confidence intervals relative to the EDL, ML, and REML methods, and it may reflect the undercoverage property of EDL, ML, and REML methods and the valid coverage property of the proposed method. Moreover, the confidence intervals were generally narrower than the 95% confidence regions of joint inferences, and the proposed methods provided asymmetric confidence intervals concerning the point estimates. For the random permutations, we conducted sensitivity analyses; the results are presented in e-Appendix G in Supporting Information.
These results suggest that there might have been greater statistical uncertainty for these results when this systematic review was published. However, conventional methods might not accurately evaluate these statistical errors. Considering the simulation results, the permutation-based methods would correct the undercoverage properties of the conventional methods, and thus provide statistically accurate results.

Network meta-analysis of antihypertensive drugs for incident diabetes
As a second illustrative example, we analyzed a dataset from a network meta-analysis to assess the effects of antihypertensive drugs on incident diabetes (Elliott and Meyer, 2007).
The network meta-analysis can be implemented using the multivariate random effects model (*), regarding the multiple comparative effect measure estimators (e.g., log-odds ratio [OR]) as multiple outcome variables (Salanti, 2012). Here, we regarded the diuretic as the reference treatment, and applied the contrast-based network meta-analysis model for the comparative OR. Note that the outcome variables comprise five comparative log OR estimates and involve partially unobserved outcomes. Here, we considered a sensitivity analysis conducted by Elliott and Meyer (2007) excluding three trials (DREAM, HOPE, PEACE; N = 19).
In Table 2, we present the results of network meta-analyses using the ML, REML and proposed permutation method. The number of permutations was set to 2400. The comparative OR estimates and 95% confidence intervals were obtained by the marginal inference methods.
For the between-studies variance-covariance matrix , we adopted a standard compound symmetry structure in which the correlation coefficient was fixed to 0.50 (Higgins and Whitehead, 1996). The estimates of the between-studies standard deviation were 0.081 and 0.114 for the ML and REML estimations, respectively. The results showed that the proposed permutation method provided wider confidence intervals than the ML and REML methods.
They also reflected the undercoverage property of the ML and REML methods, as well as the valid coverage property of the proposed method. In addition, the results of the test of the overall null hypothesis, which correspond to the incidence rates of all six treatments, are equivalent ( ⋯ 0 . The overall test corresponds to the joint inferences of the grand mean parameters. We provide the p-values of the overall test using the ML, REML, and permutation methods in Table 1; all were less than 0.001. The results were similar, but the permutation method could adequately quantify the statistical errors in these cases, based on the numerical evidence from the simulations. A league table that presents comparative OR for all possible pairwise comparisons is provided in e-Appendix H.

Discussion
In this paper, we presented permutation-based inference methods for multivariate meta-    for joint inferences. In general, the EDL and REML methods exhibited superior performances to ML method, as expected, but the coverage rates were in general less than 0.95. In marginal inferences, the undercoverage properties were less extreme, but the coverage probabilities were still generally less than 0.95. In particular, under small N and serious.
It is notable that for the permutation methods, the coverage probabilities were generally higher than the nominal level (

e-Appendix E: Simulations for trivariate meta-analysis settings
We also conducted simulation studies to evaluate the performances of the proposed permutation methods under trivariate meta-analysis settings. The simulation settings were generalized to trivariate meta-analysis settings based on those for bivariate cases by Jackson et al. (2010) and Jackson and Riley (2014). We generated the trivariate outcome variables , , of N studies (=8, 12, 16) from the random-effects model (*). We set the grand mean parameters 0 following Jackson and Riley We analyzed the generated datasets using the ML estimator and the Wald confidence region (interval), the REML estimator and the Wald-type confidence region (interval), and the proposed methods for joint inferences ( , ) and marginal inferences ( ) of the grand mean parameters. The EDL method (Chen et al., 2012;Jackson et al., 2013) was not involved in these simulations, because the EDL estimate could not be calculated in to circumvent the "singularity" of marginal covariance matrices of the outcome vectors, we adopted the Moore-Penrose generalized inverse matrix for calculating inverse matrices. Thus, in all the simulations, we could calculate the efficient score statistics. In some cases, they might not be an "inverse matrix" formally, but the resultant accurate coverage properties (corresponding to approximately 95%, accurately) were achieved by this computational procedure. As a whole, the simulation results clearly demonstrated the validity and effectiveness of the proposed methods.

e-Appendix F: Simulations for trivariate meta-analysis involving missing outcomes
In addition, we conducted simulations under trivariate meta-analysis settings that involve partially missing outcomes. The simulation settings are the same as those in e-Appendix E. Besides, in these simulations, we set 25% of the first and second outcomes , and 50% of the third outcome as not reported (measured). For the incomplete outcomes, the joint inference method using the sign-invariant moment estimator ( ) cannot be applied. Therefore, we did not involve this method in these simulations. Note that the ML and REML computations (R mvmeta package) did not converge in some cases, and we excluded them in calculating the coverage probabilities.
The results are presented in e- Figure 4. Because the statistical information became smaller owing to the missing outcomes, the coverage probabilities of the ML and REML e- Figure 3. Results of simulation for trivariate meta-analyses: Monte Carlo estimates of coverage probabilities of 95% confidence regions (joint inference) and confidence intervals (marginal inference).   Figure 4. Results of simulation for trivariate meta-analyses involving missing outcomes: Monte Carlo estimates of coverage probabilities of 95% confidence regions (joint inference) and confidence intervals (marginal inference).

e-Appendix G: Sensitivity analyses for the number of permutations
To assess the validity of the number of permutations, we conducted sensitivity analyses for the meta-analysis of diagnostic accuracy studies for airway eosinophilia in asthma.
For this example, the 95% confidence interval of sensitivity was (0.642, 0.759). These confidence limits were calculated for 2400 random permutations, and we confirmed that they agreed exactly with those calculated for all possible permutations (2 4096). The permutation number was set to 2400 to control the Monte Carlo standard error at 0.0044 when the estimated p-value is 0.05. To assess the influence of the number of random permutations, we conducted sensitivity analyses, setting the permutation number (B) to 100, 200, 400, 600, 800, 1000, 1500, 2000, 2500, 3000, 3500, or (Edgington and Onghena, 2007;Keller-McNulty and Higgins, 1987), they would be controlled at least to a certain degree by 2000 in these applications.