mmeta : An R Package for Multivariate Meta-Analysis

This paper describes the core features of the R package mmeta , which implements the exact posterior inference of odds ratio, relative risk, and risk diﬀerence given either a single 2 × 2 table or multiple 2 × 2 tables when the risks within the same study are independent or correlated.


Introduction
Epidemiological studies often involve comparisons between two populations with binary outcomes.Data from these studies are usually summarized by a single or multiple 2 × 2 tables.To quantify the association between an exposure and a certain disease, comparative measures between two risks, e.g., odds ratio (OR), relative risk (RR), and risk difference (RD), are frequently used.A Bayesian approach has been widely applied to obtain the posterior distributions of these comparative measures that reflect evidence from the data and available prior knowledge.Bayesian inference on a single study based on a 2 × 2 table has been investigated by several researchers.Specifically, Nurminen and Mutanen (1987) derived the exact posterior distributions of OR, RR, and RD under independent beta prior distributions with integer hyperparameters.Marshall (1988) extended the results of OR by using hypergeometric func-tions (Gauss 1812) to allow the hyperparameters being any positive numbers.Nadarajah and Kotz (2007) gave a formula for RD using the Appell hypergeometric function.Chen and Luo (2011) corrected the formula by Nadarajah and Kotz (2007) and further simplified the formula to avoid divergence of the Appell hypergeometric function.Hora and Kelley (1983) and Hashemi, Nandram, and Goldberg (1997) extended the results of Nurminen and Mutanen (1987) on RR to beta prior distributions with any positive hyperparameters.
Multiple 2 × 2 tables often arise in meta-analysis which combines statistical evidence from multiple studies.Two risks within the same study are possibly correlated because they share some common factors such as environment and population structure.For example, in genetic association studies, people in the same study are likely to live in the same community sharing similar environmental factors or similar ancestors (Lee 1996).Riley (2009) has showed via simulation studies that separate meta-analysis of correlated outcomes can lead to biased estimates of variances of the summary effect sizes.In contrast, multivariate meta-analysis summarizes simultaneously all outcomes of interest instead of conducting many separate univariate meta-analysis.Multivariate meta-analysis has recently received lots of attention (e.g., Reitsma, Glas, Rutjes, Scholten, Bossuyt, and Zwinderman 2005;Chu and Cole 2006;Riley, Abrams, Sutton, Lambert, and Thompson 2007;Riley, Thompson, and Abrams 2008;Hamza, Reitsma, and Stijnen 2008).An excellent overview of multivariate meta-analysis can be found in Jackson, Riley, and White (2011) and Mavridis and Salanti (2012).In the multivariate meta-analysis with a binary outcome and a categorical exposure, two modeling strategies have been commonly used: A bivariate general linear mixed-effects model on the transformed proportions (Reitsma et al. 2005;Arends, Hamza, Houwelingen, Heijenbrok-Kal, Hunink, and Stijnen 2008) and a bivariate generalized linear mixed-effects model on the transformed risks (e.g., logit or probit transformations; Houwelingen, Zwinderman, and Stijnen 1993;Houwelingen, Arends, and Stijnen 2002;Chu and Cole 2006;Chu, Guo, and Zhou 2010).However, these two methods are based on the transformed proportions or the transformed risks and thus the interpretation is transformation dependent.Multivariate meta-analysis can be conducted using various software packages including Stata (StataCorp.2011), SAS (SAS Institute Inc. 2011), R (R Core Team 2013).Specifically, the mvmeta command in Stata performs fixed-and random-effects multivariate meta-regression analysis.The SAS PROC MIXED routine was the first that popularized multivariate meta-analysis (Houwelingen et al. 2002).More recently, the SAS macro METADAS was made available to fit bivariate meta-analysis models for diagnostic test accuracy studies (Takwoingi, Guo, and Deeks 2008).R package metaSEM (Cheung 2012) can be used to conduct univariate and multivariate meta-analysis using structural equation modeling (SEM) via the OpenMx package.In addition, a new R package mvmeta (Gasparrini 2012) can perform fixed-and random-effects multivariate meta-analysis and meta-regression.
Instead of modeling the transformed proportions or transformed risks, we use a Sarmanov family of correlated beta prior distributions (referred to as Sarmanov beta prior distributions; Sarmanov 1966) to model the risks directly; see for example, Chen, Chu, Luo, Nie, and Chen (2014a).The correlation parameter can be intuitively interpreted as the correlation coefficient between risks.In addition, the Sarmanov beta prior distribution has the following advantages in modeling.Firstly, it allows for both positive and negative correlations; secondly, it only needs specification of marginal distributions and the correlation parameter, which has important advantages in Bayesian inference because it is often easier to specify and interpret an univariate prior than a bivariate prior; thirdly, it is pseudo-conjugate to binomial distributions, i.e., the Sarmanov beta prior distribution can be expressed as a linear combination of independent bivariate beta distributions (Lee 1996), which enables us to derive closedform expressions of the exact posterior distributions for study-specific comparative measures.Such closed-form expressions offer computational convenience when the exact posterior distributions of the study-specific comparative measures are also of interest.We have used the Sarmanov beta prior distribution to make exact posterior inference of some comparative measures (e.g., OR, RR, and RD; Chen et al. 2014a;Chen, Luo, Chu, Su, and Nie 2014b;Chen, Luo, Chu, and Wei 2013).This paper describes the mmeta package as a collection of a new family of models different from those in the aforementioned packages.Specifically, the inference of the overall and study-specific comparative measures (i.e., OR, RR, and RD) are inferred under the Sarmanov beta prior distributions.The functions of the mmeta package have been written in the R language, with some Fortran 77 and C routines which are interfaced through R. The package is built following the S3 formulation of R methods with dependencies on R packages HI (Petris, Tardella, and Gilks 2013) and aod (Lesnoff and Lancelot 2012).The mmeta package (currently version 2.2) is available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=mmeta.
The paper is organized as follows.In Section 2 we outline the exact Bayesian posterior inference approach.We describe the features of two main functions in the mmeta package and the analysis of two real datasets in Section 3. In Section 4, we provide a brief discussion.

Models and inference on overall comparative measures
For the ith study (i = 1, . . ., I; I is the number of studies), let n ji , y ji and p ji be the number of subjects, number of subjects experienced a certain event, and the risk of experiencing the event in the jth group (j = 1, . . ., J; J is the number of groups), respectively.For simplicity, we consider only the settings with two groups under comparison (i.e., J = 2) and the extension to cases with more than 2 groups is straightforward.We assume the following Bayesian hierarchical model.At the first stage, we assume that given the study-specific risks (p 1i , p 2i ), y 1i and y 2i are independently distributed binomial variables, i.e., (y 1i , y 2i )|(n 1i , n 2i , p 1i , p 2i ) ∼ Binomial(y 1i |n 1i , p 1i ) × Binomial(y 2i |n 2i , p 2i ). (1) This conditional independence assumption is reasonable because y 1i and y 2i are calculated using subjects from different groups.To complete the Bayesian hierarchical model, we need to impose a parametric prior distribution on the study-specific risks (p 1i , p 2i ).Here we consider a family of distributions first proposed by Sarmanov (1966), and later studied extensively by Cole, Lee, Whitmore, and Zaslavsky (1995), Lee (1996), Shubina and Lee (2004), Danaher and Hardie (2005), and Chen et al. (2014a).The Sarmanov beta prior distribution is constructed such that the marginal distribution of the random-effects in the jth group p ji is equal to a beta distribution with shape parameters (a j , b j ) and the correlation coefficient between p 1i and p 2i is ρ (Sarmanov 1966;Lee 1996).Specifically, we denote beta(p; , and δ 2 j = µ j (1 − µ j )/(a j + b j + 1).The joint prior distribution of the study-specific risks (p 1i , p 2i ), referred to as Sarmanov beta prior distribution, is where g(p 1 , p With the Bayesian hierarchical model specified in equations ( 1) and (2), the log marginalized likelihood function for the unknown hyperparameters (a where P BB (y; n, a, b) is the probability mass function of a beta-binomial distribution, i.e., The last expression in equation ( 3) has been derived by Danaher and Hardie (2005) and an outline of its derivation is provided in the Appendix for interested readers.We refer to equation (3) as Sarmanov beta-binomial model.As a benefit of using Sarmanov beta prior distributions, the log marginalized likelihood function has a closed-form expression, which avoids numerical approximation of integrals.Hence the Bayesian hierarchical model specified in equations ( 1) and (2) has great computational advantage over commonly used multivariate generalized linear mixed-effects models.When ρ = 0, the Sarmanov beta-binomial model reduces to the independent beta-binomial model, i.e., the product of two beta-binomial distributions.
The hyperparameters (a 1 , b 1 , a 2 , b 2 , ρ) can be estimated by maximizing the log likelihood log L(a 1 , b 1 , a 2 , b 2 , ρ).We implement it through R (R Core Team 2013) with the optim function, which uses a quasi-Newton method with box constraints on the ranges of parameters.Denote (â 1 , b1 , â2 , b2 , ρ) the maximum likelihood estimates based on the log likelihood function in equation (3).We use the delta method to obtain the variance of the overall comparative measures, namely the overall odds ratio estimate, and the overall risk difference estimate,

Inference on study-specific comparative measures
Denote the study-specific comparative measures The statistical evidence of these comparative measures from the ith study can be quantified by the posterior distributions, i.e., Pr(θ where θ i = OR i , RR i , or RD i and data i = (y 1i , n 1i , y 2i , n 2i ).Note that the true values of the hyperparameters (a 1 , b 1 , a 2 , b 2 , ρ) are often unknown.One solution is to simply replace the hyperparameters by their estimates.Such an approach is called empirical Bayes method (Efron andMorris 1973, 1975;Gelman, Carlin, Stern, and Rubin 2004;Carlin and Louis 2009).The coverage property of the credible intervals using the empirical Bayes method has been investigated via simulation studies in Chen et al. (2014b).The conclusion is that the credible interval without accounting for the uncertainty on the hyperparameters still performs reasonably well, when the number of studies is moderate (Chen et al. 2014b).
An important property of the Sarmanov beta prior distribution for p 1 and p 2 is that it can be written as a linear combination of independent bivariate beta distributions (Lee 1996), where v k (k = 1, . . ., 4) are weights defined by After some algebra, the posterior distribution of p 1 and p 2 given data is also a linear combination of independent bivariate beta distributions, where α j = a j + y ji , β j = b j + n ji − y ji (j = 1, 2) and the weights ω k (k = 1, . . ., 4) are defined as and the normalizing constant C is calculated as The exact posterior distributions of the comparative measures (i.e., OR, RR, and RD) under the Sarmanov beta prior distribution take the following generic form where F (•, •; •; •) denotes the hypergeometric function Gauss (1812) defined by where denotes the Appell function of the first kind defined by denotes the ascending factorial.
The function multipletables() is to conduct inference based on multiple 2×2 tables.Specifically, the hyperparameters' maximum likelihood estimates ( a 1 , b 1 , a 2 , b 2 , ρ) and the inference on the overall comparative measures are obtained as described in Section 2.1.The posterior distributions of the study-specific comparative measures can be obtained either by the exact method as stated in Section 2.2 or by the samples obtained from adaptive rejection Metropolis sampling (ARMS, Gilks, Best, and Tan 1995) implemented in the R package HI, which is an interface to the C code originally developed by Wally Gilks.The argument method can be either "exact" or "sampling" to control if either the exact posterior distributions or the ARMS samples of the posterior distributions are used.The sampling method is set as default in this package because the current version of Gauss hypergeometric and Appell functions may diverge for some studies with extremely large numbers of subjects.The posterior means of the study-specific comparative measures, the corresponding 95% equal tail credible intervals (the interval between the 2.5% and 97.5% quantiles, referred to as 95% ET CI), and the 95% highest posterior density regions (referred to as 95% HDR) are obtained based on ARMS samples of the posterior distribution.To ensure reproducibility when the sampling method is used, a random seed can be set (using set.seedcommand) before calling multipletables() and singletable() functions.Various plots can be generated by multipletables(), which will be illustrated in the rest of this section.In contrast, the function singletable() is to conduct exact posterior inference based on a single 2 × 2 table for the given prior distributions of risks.This function can be used as a sensitivity analysis tool to investigate the posterior distributions of the comparative measures under various pre-specified prior distributions.The details of the function singletable() will be given in Section 3.5.
data: A data frame that contains y1, n1, y2, n2, and studynames.The details of the data structure is described in Section 3.2.
model: A character string specifying the model.Options are "Independent" and "Sarmanov" (default)."Independent" is the independent beta-binomial model; "Sarmanov" is the Sarmanov beta-binomial model.
method: A character string specifying the method.Options are "exact" and "sampling"."exact" denotes the exact method; "sampling" (default) is the method based on ARMS samples of the posterior distribution obtained with the R package HI.
alpha: A numeric value specifying the significance level.Default value is set to 0.05.
nsam: A numeric value specifying the number of samples if method = "sampling".Default value is set to 10000.

Data structure
The structure of data in multipletables() requires the input of a data frame with five columns, y1, n1, y2, n2, and studynames.The meanings of y1, n1, y2, and n2 can vary for different study designs.Users can define their own data frame to be used in multipletables().
For example, a data frame named Bellamy based on a meta-analysis of the association between gestational diabetes mellitus and type 2 diabetes mellitus (Bellamy, Casas, Hingorani, and Williams 2009) can be defined as follows.

Example: colorectal dataset
The dataset colorectal consists of data from twenty published case-control studies of the N-acetyltransferase 2 (NAT2) acetylation status and colorectal cancer risk.NAT2 is a lowpenetrance gene that regulates metabolizing enzymes.The activity of the enzymes is classified as rapid and slow acetylators.Ye and Parry (2002) investigated the association between rapid NAT2 acetylator status (event) and colorectal cancer (case) by conducting a meta-analysis based on twenty published case-control studies from January 1985 to October 2001.The data are summarized in Table 1.We define the odds ratio as the ratio of odds of having rapid NAT2 acetylator status comparing those with colorectal cancer to those without.The colorectal dataset example takes around 3 minutes and 1 minute to run using "exact" and "sampling" methods (10,000 samples), respectively.
To start analyzing the dataset, we first load the mmeta package and the colorectal dataset.
R> library("mmeta") The likelihood ratio test of H 0 : ρ = 0 yields a p value of 0.08 with χ 2 test statistic being 3.152.The estimates of the hyperparameters, the estimated mean and the 95% confidence interval (CI) of the overall odds ratio are provided.In addition, the posterior means and the 95% credible intervals (CI) of all study-specific odds ratios are given.If the argument model = "Independent", the independent Beta-Binomial model is fitted to the dataset.If the argument method = "sampling", adaptive rejection Metropolis sampling implemented in R package HI is used to obtain the posterior inference.
The forest plot with the 95% CI of the overall odds ratio and the 95% CIs of the study-specific odds ratios as shown in Figure 1 can be obtained using the plot function with the argument type = "forest".
R> plot(multiple.OR, type = "forest", addline = 1, xlabel=c(0.3,0.5,2,4)) The argument addline is to add a blue dotted reference line to the plot.If the argument file is specified, (e.g., file = "multiple_OR_forest"), the plot will be saved as Figure 1: Forest plot of 20 study-specific and the overall odds ratios with 95% CIs.
"./mmeta/multiple_OR_forest.eps",where "./" denotes the current working directory and the directory mmeta is created automatically if it does not exist.The argument select (e.g., select = 1:4) can be set to select multiple target studies to be displayed.If the argument ciShow = TRUE (by default), the numbers of the posterior means and the CIs will be displayed at the right side of the forest plot.Many standard R plotting arguments (e.g., xlabel, ylabel, ylim, xlim) can be set in the plot function.See the help file for more details.
The posterior density functions of some target studies can be overlaid as shown in Figure 2 with the argument type = "overlap".
The posterior density functions of these target studies can be viewed in a side-by-side manner as in Figure 3 if the argument type = "sidebyside", where both the prior and posterior distributions are displayed.
R> plot(multiple.OR, type = "sidebyside", select = c(4, 14, 16, 20), + ylim = c(0, 2.7), xlim = c(0.5,1.5)) Figure 3 displays the prior and posterior distributions of the study-specific odds ratios for these four target studies.Such a plot is useful to investigate the difference between the prior and posterior distributions, hence the strength of contribution from individual studies.

Example: withdrawal dataset
Tricyclic antidepressants are effective in preventing headaches and have become a standard modality in headache prevention.To investigate the efficacy and related adverse effects of tricyclic antidepressants in the treatment of headaches, Jackson et al. (2010) reported a metaanalysis based on multiple clinical trials from year 1964 to year 2009.Among several outcomes of interest, proportion of withdrawal during a trial is paid special attention because it is a very important measure of adverse effects and it plays a critical role in drug safety.One question of interest is whether the probability of withdrawing due to adverse effects is increased by the tricyclic treatment compared with the placebo.This can be measured by relative risk (defined as the ratio of risks of withdrawal comparing those in the tricyclic treatment group to those in the placebo group).The numbers of withdrawals due to adverse effects in sixteen clinical trials are summarized in Table 2.The withdrawal dataset example takes around 2 minutes and 1 minute to run using the "exact" and the "sampling" method (with 10,000 samples), respectively.
To start analyzing the dataset, we first load the withdrawal dataset.
R> data("withdrawal", package = "mmeta") The available data have the following structure R> str(withdrawal)  The likelihood ratio test of H 0 : ρ = 0 yields a p value of 0.65 with χ 2 test statistic being 0.207.The estimates of the hyperparameters, the estimated mean and the 95% CI of the overall relative risk are provided.In addition, the posterior mean and 95% CI of each studyspecific relative risk are given.The forest plot with the confidence interval of the overall relative risk and the credible intervals of the study-specific relative risks as shown in Figure 4 can be obtained using the plot function with the argument type = "forest".
Moreover, the posterior density function of each target study can be viewed in a side-by-side manner as in Figure 6 if the argument type = "sidebyside", where both the prior and posterior distributions are displayed.
Because the estimated relative risk for the study of Loldrup, Langemark, Hansen, Olesen, and Bech (1989) (mean: 6.325, 95% CI: [3.749, 10.782]) is much larger than those for the other studies, it can be potentially influential to the analysis results.To evaluate the influence of this study, we remove it and reanalyze the dataset by calling the function multipletables().
The results can be investigated using summary() (not shown).
R> multiple.RR.sens <-multipletables(data = withdrawal [-11, ], + measure = "RR", model = "Sarmanov") R> summary(multiple.RR.sens) The likelihood ratio test of zero correlation coefficient results in a p value of 0.40 with the χ 2 test statistics being equal to 0.707.Although the study of Loldrup et al. (1989) slightly changes the overall and study-specific relative risk estimates, it is not influential because the overall relative risk estimates are not significant regardless of it being included or not.
If the risk difference is of interest, we define it as the difference of risks of withdrawal comparing those in the treatment group to those in the placebo group.The function multipletables() is called to conduct posterior inference of the risk differences.The forest plot, the overlaid, and side-by-side plots of the posterior density functions for risk difference can also be obtained by using the plot function with the argument type = "forest", type = "overlap", and type = "sidebyside", respectively.

R>
Note that the study of Loldrup et al. (1989) can be influential on the analysis results because the estimated risk difference (mean: 0.593, 95% CI: [0.509, 0.668]) is much larger than those in other studies.To evaluate the sensitivity of the inference on this study, we remove it and reanalyze the dataset by calling the function multipletables().The results can be investigated using summary() (not shown).
R> multiple.RD.sens <-multipletables(data = withdrawal [-11, ], + measure = "RD", model = "Sarmanov") R> summary(multiple.RD.sens) The likelihood ratio test of zero correlation coefficient results in a p value of 0.40.Although the study of Loldrup et al. (1989) slightly changes the overall and study-specific risk difference estimates, it is not influential because the overall risk difference estimates are not significant regardless of it being included or not.

Using the singletable() function
When the inference on a specific study based on a single 2 × 2 table is of interest, the function singletable() can be used as a sensitivity analysis tool to conduct the exact posterior inference on comparative measures under various prior distributions.The arguments used in a call to the function singletable() are singletable <-function(y1, n1, y2, n2, measure, model = "Sarmanov", method = "sampling", alpha = 0.05, nsam = 10000) In the following we summarize the main arguments of singletable().
y1, n1: Integers indicating the number of events and the total number of subjects in group 1.
y2, n2: Integers indicating the number of events and the total number of subjects in group 2.
model: A character string specifying the model.Options are "Independent" and "Sarmanov" (default)."Independent" is the independent beta-binomial model; "Sarmanov" is the Sarmanov beta-binomial model.
method: A character string specifying the method.Options are "exact" and "sampling"."exact" denotes the exact method; "sampling" (default) is the method based on ARMS samples of the posterior distribution obtained with the R package HI.
rho: A numeric value specifying the correlation coefficient for the Sarmanov beta prior distribution.Default value is set to 0.
alpha: A numeric value specifying the significance level.Default value is set to 0.05.
nsam: A numeric value specifying the number of samples if method is sampling.Default value is set to 10,000.
To illustrate the use of the function singletable(), we consider the study by Ladero et al. (1991)  single study of two populations with binary outcomes.The theory used for model fitting is summarized briefly, and the two major functions (multipletables() and singletable()) of the package are described in details.Practical use of the mmeta package is illustrated with two examples of meta-analysis based on multiple 2 × 2 tables and one example of a single 2 × 2 table.As a future research direction, we would like to expand the functionality of this package to conduct meta-regression analysis using the Sarmanov beta prior distributions as illustrated in Chen et al. (2014a) and Chen et al. (2014b).Moreover, we have investigated many available non-commercial algorithms and software packages to compute the hypergeometric function and the Appell function.To the best of our knowledge, we cannot find one that provides stable computation for studies with extremely large numbers of subjects.Developing a robust algorithm for computation of these two functions are also part of our future research.

Figure 2 :
Figure2: Posterior distributions of study-specific odds ratios for four selected studies in one plot.

Figure 3 :
Figure 3: Posterior distributions of study-specific odds ratios for four studies in four separate plots.

Figure 5 :
Figure 5: Posterior distributions of study-specific relative risk for four studies in one plot.

Figure 6 :
Figure 6: Posterior distributions of study-specific relative risks for four studies in four separate plots.

Table 2 :
Jackson et al. (2010)tables() is called to conduct exact posterior inference of relative risks.Data from a meta-analysis of sixteen studies on the association between withdrawal due to the adverse effects and the tricyclic treatment inJackson et al. (2010).No.events: Number of individuals who withdrew from the study.No.individuals: Number of individuals who started the study.Figure 4: Forest plot of 16 study-specific and the overall relative risks with 95% CIs.