Flexible random-effects distribution models for meta-analysis

In meta-analysis, the random-effects models are standard tools to address between-study heterogeneity in evidence synthesis analyses. For the random-effects distribution models, the normal distribution model has been adopted in most systematic reviews due to its computational and conceptual simplicity. However, the restrictive model assumption might have serious influences on the overall conclusions in practices. In this article, we first provide two examples of real-world evidence that clearly show that the normal distribution assumption is unsuitable. To address the model restriction problem, we propose alternative flexible random-effects models that can flexibly regulate skewness, kurtosis and tailweight: skew normal distribution, skew t-distribution, asymmetric Subbotin distribution, Jones-Faddy distribution, and sinh-arcsinh distribution. We also developed a R package, flexmeta, that can easily perform these methods. Using the flexible random-effects distribution models, the results of the two meta-analyses were markedly altered, potentially influencing the overall conclusions of these systematic reviews. The flexible methods and computational tools can provide more precise evidence, and these methods would be recommended at least as sensitivity analysis tools to assess the influence of the normal distribution assumption of the random-effects model.


Introduction
In meta-analysis in medical studies, random-effects models have been the primary statistical tools for quantitative evaluation of treatment effects that account for betweenstudies heterogeneity 1,2 . Conventionally, the normal distribution assumption has been adopted in most systematic reviews due to its computational and conceptual simplicity 2,3 . However, the shape of the random-effects distribution reflects how the treatment effects parameters (e.g., mean difference, log relative risk) are distributed in the target population, and are directly associated with the fundamental heterogeneity of treatment effects. If the normal distribution assumption diverges drastically from the true heterogeneous structure, the overall results of the meta-analyses may be misleading. In addition, in recent studies, prediction intervals have been gaining prominence in meta-analyses as a means to quantify heterogeneity and effectiveness in real-world uses of the treatment 4,5 . Because the prediction interval is constructed by the estimated random-effects distribution, it should be directly influenced by the form of the distribution assumptions. Several papers have discussed the flawed uses of the normal distribution model in meta-analyses 6-8 , but to date only a few effective methods have been developed to address this issue, and there are no useful statistical packages that can be handled by non-statisticians. Also, there is a lack of real-world evidence that clearly demonstrates the relevance of this issue.
In this article, we propose random-effects meta-analysis methods with flexible distribution models that can flexibly express skewness, kurtosis, and tailweight: (1) skew normal distribution 9, 10 , (2) skew t-distribution 10,11 , (3) asymmetric Subbotin distribution 10,12 , (4) Jones-Faddy distribution 13 , and (5) sinh-arcsinh distribution 14 . Via application of these five flexible random-effects distribution models to two recently published systematic reviews 15,16 , we will demonstrate that the overall conclusions and interpretations of meta-analyses can be dramatically altered if the normal distribution assumption is not suitable. In addition, we provide a new R package, flexmeta, that can perform meta-analysis with simple code using the flexible random-effects distributions.
We will explicitly show that the implicit uses of the normal distribution assumption might yield misleading results, and that our flexible alternative distributions may provide more valid conclusions for health technology assessments and policy making.

Descriptions of two motivating meta-analyses
We searched for recently published systematic reviews in leading medical journals (e.g., BMJ, JAMA), and found two examples 15,16 that clearly demonstrated the unsuitability of the normal distribution assumption. The first example is a meta-analysis by Rubinstein et al. 15 assessing the benefits and harms of spinal manipulative therapy (SMT) for the treatment of chronic lower back pain. In Figure 1(a), we present a forest plot of their meta-analysis of 23 randomized controlled trials assessing the reduction of pain at 1 month (0-100; 0 = no pain, 100 = maximum pain) for SMT (N = 1629) vs. recommended therapies (N = 1526). The effect measure was the mean difference (MD). Using the ordinary random-effects meta-analysis method based on the normal random-effects distribution, we identified a substantial heterogeneity of the treatment effects, I 2 = 92%, τ 2 = 112.20 (P < 0.01; Cochrane's Q-test). The between-studies heterogeneity should be addressed in synthesis analysis. However, most of the MD estimates fell within relatively narrow range around the mean, although a small number exhibited larger effect sizes. This might imply that the true MD distribution is a skewed, heavy tailed, and sharp distribution.
Although the average MD was estimated as −3.17 (95%CI: −7.85, 1.51) by the DerSimonian-Laird method 3 , the point and interval estimates depend on the normal distribution assumption.
The second example is a meta-analysis by Koutoukidis et al. 16 aimed at estimating the association of weight loss interventions with biomarkers of liver disease in nonalcoholic fatty liver disease. In Figure 1(b), we also present a forest plot of their metaanalysis of 25 randomized controlled trials that assess the weight loss (kg) for moreintensive weight loss interventions (N = 1496) vs. no or lower-intensity weight loss interventions (N = 1062). Again, the effect measure was the MD, and we identified a substantial heterogeneity of the treatment effects: I 2 = 95%, τ 2 = 12.45 (P < 0.01; Cochrane's Q-test). In this case, the MD estimates were not symmetrically distributed, and a certain number of trials exhibited a larger intervention effect than the average MD of −3.51 (95%CI: −5.03, −2.00). Thus, the true MD distribution would be a skewed, heavy tailed distribution. In particular, in predicting the intervention effect of a future trial, the normal distribution model would not suitably fit this dataset. Although the ordinary 95% Higgins-Thompson-Spiegelhalter (HTS) prediction interval 5 was (−11.02, 3.99), it might not express the true nature of the intervention effects in the target population.

The flexible random-effects distribution models
To address the restriction problem of the normal distribution, we propose random-effects meta-analysis methods using five flexible random-effects distribution models. For the notation, we consider that there are studies to be synthesized, and that 1,2, … , is the estimated treatment effect measure in the th study, e.g., mean difference, odds ratio, and hazard ratio; the ratio measures are typically transformed to logarithmic scales. The random-effects models considered here 2, 3 are then defined as where is the true effect size of the th study, and is the within-studies variance, which is usually assumed to be known and fixed to valid estimates. Also, corresponds to the random-effects distribution that expresses the heterogeneous probability distribution of . For the conventional normal-normal random-effects model, corresponds to a normal distribution. The predictive interval for future study 4 is substantially constructed based on the random-effects distribution . To overcome the limitations on the expressive ability of the normal distribution, our proposal is to adopt alternative flexible probability distributions.
Currently, due to developments in statistical distribution theory, various flexible probability distributions are available 10  The skew normal distribution 9, 10 is a generalized version of the conventional normal distribution that allows for skewness. is the location parameter that regulates the center location, and ω is the scale parameter that regulates the dispersion of the distribution; we use these notations similarly for the following four distributions as well; in the graphical displays in Figure 2, we set these parameters to 0 and 1, consistently. Also, is the skewness parameter that adjusts the skewness, and the distribution is positively (negatively) skewed for 0 ( 0 ). When 0 , it accords to a normal distribution N( , ). In Figure 2 The skew t-distribution 11 is also a generalized version of the conventional Student tdistribution that allows for skewness. The t-distribution can express heavy tailweight and kurtosis relative to the normal distribution via varying the degree of freedom ν (> 0). In this case, is the skewness parameter; the distribution is positively (negatively) skewed for 0 ( 0). Also, for 0, it accords to the ordinary t-distribution. In Figure   2

Methods for the treatment effect estimation and prediction
For the random-effects model (*), we can adopt the flexible distribution models for the random-effects distribution . The average treatment effect can be addressed as the mean μ of . As in the conventional DerSimonian-Laird-type normal-normal model, the parameters of can be estimated by frequentist methods (e.g., the maximum likelihood estimation), but in many cases, they require complex numerical integrations;

Applications to the two meta-analyses
We applied the flexible random-effects models to the two meta-analysis datasets described above. As reference methods, we also conducted the same analyses using the normal and t-distribution models. We used R ver. 3 and predictive intervals (PI) were calculated using the posterior samples of the mean of and the predictive distribution of the effect of a future study ~ from MCMC. To evaluate the impacts of adopting the flexible distribution models rather than the ordinary normal distribution, we present graphical displays of the posterior and predictive distributions. In addition, we assessed model adequacies by DIC.

Results
In Table 1  respectively, whereas that of the normal random-effects distribution model was 0.90. In the original paper by Rubinstein et al. 15 , the overall MD test was not statistically significant at the 5% level. However, the overall results were clearly altered by adopting the skewed flexible distribution models, which strongly indicated that the true effect sizes would lie within a narrower range and would be skew distributed. The overall conclusion for the overall MD could be changed using the flexible models. Also, we present a summary of the predictive distribution of this example in Table 2 (a). These results also indicated that the predictive distribution would be strongly skewed.
For the second example, the meta-analysis of non-alcoholic fatty liver disease, we present summaries of the posterior distribution of μ and the predictive distribution in Table 1 respectively. Hence, the overall conclusions might be altered.

Discussion
Conclusions obtained from meta-analyses are widely applied to public health, clinical practice, health technology assessments, and policy-making. If misleading results have been produced by inadequate methods, the impact might be enormous. In this article, we proposed effective methods for meta-analysis using flexible random-effects distribution models, and provided an easily implementable statistical package for these methods.
Through illustrative examples, we clearly showed the restrictions of using the conventional normal random-effects distribution model, which may yield misleading conclusions. The flexible random-effects distribution models represent effective tools for preventing such an outcome. Conventionally, these MCMC computations require special software and high-performance computers; to address these obstacles, we developed a user-friendly package, flexmeta, which was designed to be easily handled by non-statisticians and is freely available online. The proposed methods and the developed tools would help us to provide precise evidence. At a minimum, we recommend using these methods in sensitivity analyses.
In this study, we adopted five flexible distributions, but other probability distributions exist in statistical theory, e.g., see the comprehensive textbook by Azzalini and Capitanio 10 . Other choices might also be considered, but the five distribution models discussed here have sufficient expressive abilities, and significantly different results would not be obtained by adopting other existing distributions. Another choice would be to adopt nonparametric methods 18,19 . However, in meta-analysis in medical studies, the number of studies K is usually not large 22,23 ; consequently, nonparametric methods would be unstable in many applications because they require much larger statistical information (parallel to K) to conduct valid estimation and prediction.
As shown by the real data applications, existing meta-analyses may have reached misleading conclusions due to the straightforward uses of the normal random-effects distribution model. Our proposed methods might change the overall conclusions of these meta-analyses, and systematic re-evaluation of existing meta-analyses would be an interesting topic for future studies. In addition, for future systematic reviews, the flexible methods might be used as standard methods to provide accurate conclusions, at least in sensitivity analyses.

A.1 Skew normal distribution
The skew normal distribution ) is a generalized version of the conventional normal distribution that allows for skewness. The probability density function is where ( ) and Φ( ) are, respectively, the probability density and cumulative distribution functions of the standard normal distribution N(0, 1). is the location parameter that regulates the center location, and ω (> 0) is the scale parameter that regulates the dispersion of the distribution. is the skewness parameter that adjusts the skewness; the distribution is positively (negatively) skewed for > 0 ( < 0). When = 0, the skew normal distribution accords with the normal distribution N( , ). The mean and variance of this distribution are

A.2 Skew t-distribution
The skew t-distribution Capitanio, 2003, 2014) is also a generalized version of the conventional Student t-distribution that allows for skewness. The probability density function is where ( | ) and ( | ) are, respectively, the probability density and cumulative distribution functions of the Student t-distribution with (> 0) degrees of freedom.
is the location parameter that regulates the center location, and ω (> 0) is the scale parameter that regulates the dispersion of the distribution. is the skewness parameter that adjusts the skewness and the distribution is positively (negatively) skewed for > 0 ( < 0). When = 0, this distribution accords with the Student t-distribution with degrees of freedom. The mean and variance of this distribution are The skew t-distribution can express flexible shapes by controlling the degree of freedom, compared with skew normal distribution, especially for the kurtosis and tailweight.

A.3 Asymmetric Subbotin distribution
Subbotin (1923) proposed a symmetric probability distribution that can regulate the kurtosis and tail thickness flexibly. The probability density function of the Subbotin distribution with (> 0) degrees of freedom is Based on the form of the probability density function, this distribution involves a double exponential and trapezoidal-shaped distributions as special cases. The location and scale can be regulated by linear transmission, and the distribution can flexibly express heavy and light tailweight. The asymmetric Subbotin distribution of type II (AS2) ) is an extended version of this distribution, and the probability density function is is the location parameter that regulates the center location, and ω is the scale parameter that regulates the dispersion of the distribution. The distribution is positively (negatively) skewed for > 0 ( < 0), and more kurtosed for smaller ν. The mean and variance of this distribution are As shown in Figure 2, the AS2 distribution can express a sharp skew distribution, which can be seen as an asymmetric double exponential distribution. Also, it can express a more rounded shape, like the skew t-distribution.

A.4 Jones-Faddy distribution
Jones and Faddy (2003) proposed another skewed version of t-distribution, whose probability density function is expressed as

e-Appendix B: Methods for the Bayesian modelling
For the Bayesian random-effects model, the prior distributions for the model parameters of ( ) directly influence to the posterior inferences and predictions. We adopted non-informative priors in the applications, and the default settings of flexmeta also adopt the same prior distributions.
Here, we explain which priors were adopted; the source R and Stan codes are available at our GitHub site (https://github.com/nomahi/flexmeta), and the users can freely customize the prior distribution settings.
For the location and scale parameters and , we consistently adopted the following vague prior distributions for the seven models (involving the ordinary normal and t-distribution models): ~(0, 100 2 ) ~(0, 20) For the degree-of-freedom parameter of the t-distribution, skew t-distribution, and AS2 distribution, we adopted an exponential (0.1) prior that was restricted to k > 2.5 to assure the existence of the second moment (k ≥ 2), as in Fernandez and Steel (1998) and . For the skewness parameter of the skew normal distribution, skew-t distribution, and AS2 distribution, we adopted a proper vague normal prior (0, 5 2 ). For the Jones-Faddy distribution, we assumed uniform priors for the two model parameters and , ,~(1.5, 200) . The lower bound of the uniform distribution is determined to assure the existences of the first, second and third moments . For the SAS distribution, we also adopted vague priors for the skewness and kurtosis parameters, ~(0,100 2 ) and ~(0,100).