Copula-Based Estimation Methods for a Common Mean Vector for Bivariate Meta-Analyses

Traditional bivariate meta-analyses adopt the bivariate normal model. As the bivariate normal distribution produces symmetric dependence, it is not flexible enough to describe the true dependence structure of real meta-analyses. As an alternative to the bivariate normal model, recent papers have adopted “copula” models for bivariate meta-analyses. Copulas consist of both symmetric copulas (e.g., the normal copula) and asymmetric copulas (e.g., the Clayton copula). While copula models are promising, there are only a few studies on copula-based bivariate meta-analysis. Therefore, the goal of this article is to fully develop the methodologies and theories of the copula-based bivariate meta-analysis, specifically for estimating the common mean vector. This work is regarded as a generalization of our previous methodological/theoretical studies under the FGM copula to a broad class of copulas. In addition, we develop a new R package, “CommonMean.Copula”, to implement the proposed methods. Simulations are performed to check the proposed methods. Two real dataset are analyzed for illustration, demonstrating the insufficiency of the bivariate normal model.


Introduction
Bivariate outcomes often arise in meta-analyses on scientific studies, such as education and medicine. Educational researchers may analyze bivariate exam scores on verbal and mathematics [1,2], or on mathematics and statistics [3]. Medical experts may analyze bivariate risk scores on myocardial infection and cardiovascular death for diabetes patients [4,5]. Bivariate meta-analyses are statistical methods designed for these metaanalytical studies [6]. Dependence between two outcomes should be considered while performing bivariate meta-analyses. If one simply considers univariate (marginal) analysis for each outcome separately, any possible dependence between the outcomes is ignored. Riley [2] and Copas et al. [7] showed that ignoring the dependence between two outcomes increases the error for estimating parameters due to the loss of information. In medical research, dependence itself can be of clinical importance, e.g., dependence between two survival outcomes in meta-analysis [8][9][10][11].
In the traditional bivariate meta-analyses, the parameters of interest are the means of a bivariate normal model [6]. However, the bivariate normal model is not flexible enough to describe the true dependence structure of real meta-analyses. It will be shown that the bivariate normal mode fits poorly to the dependence structure of real bivariate meta-analyses (Section 8). This has motivated researchers to consider alternative models.
As an alternative to the bivariate normal model, recent papers have adopted "copula" models for bivariate meta-analyses [3,5,[12][13][14][15]. Copula models are flexible as they allow a variety of dependence structures. Copulas consist of both symmetric copulas (e.g., the normal Symmetry 2022, 14, 186 2 of 27 copula) and asymmetric copulas (e.g., the Clayton copula). Copula models have become very popular in all areas of science by replacing the traditional multivariate normal models. In astronomy, Takeuchi [16] constructed the bivariate luminosity density functions using the FGM copula; see reference [17] for the application of the FGM copula to engineering. In ecology, Ghosh et al. [18] applied copulas to model the dependence structure in environmental and biological variables. In environmental science, Alidoost et al. [19] used bivariate copulas in the analysis of temperature. See the survey of [20] for applications to energy, forestry, and environmental sciences. The books of [21,22] are devoted to the applications of copulas in survival analysis; see also references [11,[23][24][25].
While bivariate copula models for meta-analyses are promising, there are only a few methodologically and theoretically solid studies on copula-based bivariate meta-analysis. For instance, the detailed theoretical studies of [3] are limited to the FGM copula. Other copula-based meta-analyses published in biostatistical journals, such as [5,[12][13][14][15], are proposed without theoretical details. Furthermore, copula-based bivariate meta-analyses have not been implemented in a free software environment.
Therefore, the goal of this article is to fully develop the methodologies and theories of the copula-based bivariate meta-analysis for estimating the common mean vector. This work is regarded as a large generalization of our previous methodological/theoretical studies under the FGM copula model [3] to a broad class of copula models. In this article, we obtain theoretical results, including the formula of the information matrix and large sample theories. Our theoretical results guarantee the applications of many copulas, such as the Clayton, Gumbel, Frank, and normal copulas, in addition to the FGM copula. In addition, we developed a new R package, "CommonMean.Copula" [26], to implement the proposed methods under the five copulas. Therefore, the aim of the article is to make a solid development of the methodologies, theories, and practical implementations of copula-based bivariate meta-analysis for the common mean, which are not yet available in the literature.
The article is organized as follows. Section 2 reviews the background of this research. Section 3 introduces the proposed model and estimator. Section 4 provides the asymptotic theory and Section 5 gives confidence sets. Section 6 introduces our new R package. Section 7 conducts simulations to check the accuracy of the proposed methods. Section 8 analyzes two real datasets for illustration. Section 9 extends the proposed methods to non-normal data. Finally, Section 10 concludes with a discussion.

Background
This section reviews the literature on bivariate meta-analyses and the concept of copulas.

Bivariate Meta-Analysis
We review the bivariate meta-analysis method for bivariate continuous outcomes [6,27]. For each study i, let the bivariate outcomes, Y i1 and Y i2 , follow a bivariate normal distribution where ρ i ∈ (−1, 1) is the within-study correlation for each i. In Equation (1), all the responses (Y i s) share the common mean vector (µ). The covariance matrix Ω i is assumed to be known (from the i-th study) in usual bivariate meta-analyses. We do not consider a setting where the covariance is unknown [28,29]. Then, the MLE of the common mean vector is quite easily computed aŝ µ Normal n = μ Normal n,1 µ Normal One could use the R package mvmeta [30], although the above computation is easy. The bivariate normal model (1) does not allow for a different dependence structure between the two outcomes. In practice, the bivariate normal model (1) can be too restrictive, as there are various dependence patterns between two outcomes. For example, to model the luminosity function of galaxies, Takeuchi [16] pointed out that the FGM copula model offers a more ideal shape than the normal copula model from a physical point of view. Such a limitation motivates us to construct a general copula model that can describe various dependence structures.

Copulas
This subsection prepares the basic terms on copulas that will subsequently be used. A copula is a bivariate distribution function whose margins are uniformly distributed on the unit interval [31,32]. Copulas are indispensable tools when modelling a dependence structure between two random variables. We specifically consider the following parametric copulas.
The normal copula: The copula function is where Φ ρ (·, ·) is the cumulative distribution function (CDF) of the bivariate standard normal distribution with correlation ρ and Φ −1 is the inverse of the standard normal CDF Φ. While this copula is easy to understand, it has a complex form involving two implicit functions Φ ρ and Φ −1 . The following two copulas provide simpler forms than the normal copula. The Farlie-Gumbel-Morgenstern (FGM) copula [33]: The copula function is The FGM copula has a very simple form, and is a fundamental copula, which has been extended to a variety of copulas, called the generalized FGM copulas [34][35][36][37][38].
The Clayton copula [39]: The copula function is The Clayton copula is one of the simplest and most frequently used copulas in applications. The Clayton copula is derived from the gamma frailty model, leading to its remarkable popularity in survival data analysis [22,40]. It has a lower tail dependence [31], but is not tractable for modeling negative dependence. The Gumbel copula [41]: The copula function is The Gumbel copula is a popular copula with upper tail dependence [31]. The Gumbel copula does not offer a negative dependence, as in the Clayton copula. The Frank copula [42]: The copula function is The Frank copula does not have tail dependence [31]. Unlike the Clayton and Gumbel copulas, it can model both positive and negative dependences as the normal copula.
Under the null parameter (e.g., θ = 0), all the above copulas reduce to the independence copula Π(u, v) = uv. As the parameter departs from the null, the dependence gets stronger. We define the notations for partial derivatives (if they exist) as For instance, where C [1,1] is called the copula density. The copula is symmetric if C [1,1] . This means that the normal and FGM copulas are symmetric while the Clayton and Gumbel copulas are asymmetric. This symmetry should not be confused with the exchangeability C(u, v) = C(v, u). All the aforementioned parametric copulas are exchangeable.

Proposed Methods
This section proposes a general copula-based approach for estimating a bivariate common mean vector. We first define the bivariate copula model and provide sufficient conditions for the copula parameter to be identifiable. We then develop a maximum likelihood estimator (MLE) for the common mean vector. In addition, we derive the expression for the information matrix.

General Copula Model for the Common Mean
This subsection proposes a new model for estimating the common mean in bivariate meta-analyses. For Here, we call µ = (µ 1 , µ 2 ) the 'common mean vector' since it is common across i = 1, 2, . . . , n. Our target is the estimation of µ when Ω i , i = 1, 2, . . . , n are known. In general, Ω i = Ω j for some i = j, and, therefore, the random vectors Y i , i = 1, 2, . . . , n are independent but not identically distributed (i.n.i.d.). While the marginal normality is specified, the bivariate normality is unspecified. We only specify the equation We now specify a bivariate distribution for Y i . According to Sklar's Theorem [43], for copulas C θ i , i = 1, 2, . . . , n, we define the bivariate CDFs However, since ρ i is known, the copula can be restricted. To see the problem clearly, we define the correlation function ρ C : Θ → R C as where R C ≡ {ρ C (θ) : θ ∈ Θ} denotes the range of ρ C that depends on the choice of C θ . The correlation function ρ C does not depend on µ. For the copula to be useful in real meta-analyses, θ i has to be identifiable from ρ i . This means that one has to be able to solve the equation ρ C (θ) = ρ. Now, we define our general copula model for a bivariate common mean vector.
To explain the flexibility and generality of our model, we give examples for C θ i .

Example 1.
(the normal copula): Under the normal copula, the model in Equation (2) becomes Under this model, the correlation function is the identity function ρ C Normal (ρ) = ρ. In addition, one has the copula parameter space Θ C Normal = (−1, 1) , and the range of correlations R C Normal = (−1, 1). Without doubt, for any ρ i ∈ (−1, 1) , the copula parameter can be identified.

Example 2.
(the FGM copula): Under the FGM copula, the model in Equation (2) becomes Under this model, the correlation function is ρ C FGM (θ) = θ/π for −1 ≤ θ ≤ 1 [44]. Thus, the copula parameter is identified by θ i = πρ i , as long as Hence, θ i can still be identified by ρ i . This boundary enforcement is illustrated in Figure 1.
To explain the flexibility and generality of our model, we give examples for .

Example 2. (the FGM copula): Under the FGM copula, the model in Equation
Under this model, the correlation function is FGM ( ) = / for −1 ≤ ≤ 1 [44]. Thus, the copula parameter is identified by = , as long as Hence, can still be identified by . This boundary enforcement is illustrated in Figure 1.  for > 0. The correlation function does not have a closed-form, and is written as for α i > 0. The correlation function does not have a closed-form, and is written as It is known that lim α→0 ρ C Clayton (α) = 0 and lim α→∞ ρ C Clayton (α) = 1. In addition, if α 2 ≥ α 1 then [45]. Then, we conclude that the range of the correlation is R C Clayton = (0, 1). Thus, one can identify by solving ρ C Clayton (α i ) = ρ i numerically if ρ i > 0. If ρ i ≤ 0 , we suggest the independence model ( Figure 1) Similar to the Clayton copula, the correlation function does not have a closed-form, and is not displayed here. It is known that ρ C Gumbel (1) = 0 and lim β→∞ ρ C Gumbel (β) = 1. If ρ i < 0 , we suggest the independence model as in the Clayton copula.
Again, the correlation function does not have a closed-form, and is not displayed here. It is known that lim γ→−∞ ρ C Frank (γ) = −1 and lim γ→∞ ρ C Frank (γ) = 1. Thus, the Frank copula parameter does not require boundary correction.

Statistical Inference Methods
This subsection develops statistical inference methods under the proposed model. We propose the MLE for µ under the general copula model (Definition 1) in Equation (2).
Suppose that the copula density C exists. Then, the joint density of Y i is where y = (y 1 , y 2 ) and ϕ(·) is the density of N(0, 1). Given the samples, the log-likelihood function is The MLE of the common mean vector is defined aŝ where R = (−∞, ∞) is a real line. The MLE does not have a closed-form expression except for the normal copula. Thus, the MLE can also be obtained by the Newton-Raphson algorithm or some software functions (e.g., the R functions optim or nlm). One may also apply our R package CommonMean.Copula [26], which will be explained in Section 6. where Theorem 1 can be proved by straightforward calculations as Lemma 1 (Appendix A.1.). Theorem 1 helps us interpret the role of the copula C θ i on the information matrix.

Theorem 2. The determinant of I i can be expressed as
In addition, det(I i ) > 0 and I i is positive definite. Proof of Theorem 2. The expression of det(I i ) is obtained by straightforward calculations. Clearly, we have |ρ C (θ i )| < 1. Then, by the Cauchy-Schwarz inequality, Furthermore, by the arithmetic-geometric mean inequality, we have Then we obtain det( Then, by Theorem 1, the information matrix in Equation (3) becomes and its determinant is det(I Normal Example 7. (The FGM copula): Under the FGM copula Then, by Theorem 1, the information matrix in Equation (3) becomes By Theorem 2, its determinant is This result agrees with [3] who considered the FGM model.
Then, by Theorems 1 and 2, we obtain I ) accordingly.

Asymptotic Theory
To assess the sampling variability ofμ n , its asymptotic distribution is presented in this section.
A technical burden comes from the fact that our samples Y i , i = 1, 2, . . . , n are independent and non-identically distributed (i.n.i.d.) owing to heterogeneous variances (Ω i = Ω j , i = j). The existence of the asymptotic distribution requires the stabilization of the information matrix [3,46,47] in large samples. For the asymptotic variance ofμ n , to be defined, we assume the existence of a 2 × 2 positive definite matrix I ≡ lim n→∞ ∑ n i=1 I i /n. We further assume that the copula's derivatives C [4,1] θ , C [3,2] θ , C [2,3] θ , and C [1,4] With these conditions and many other technical conditions given in [48], we establish the consistency and asymptotic normality ofμ n : Theorem 3. Under the copula model (Definition 1), if some regularity conditions hold, then (a) Existence and consistency: With probability tending to one, there exists the MLÊ µ n = (μ n,1 ,μ n,2 ) such thatμ n → p µ , as n → ∞ ; (b) Asymptotic normality: The proof of Theorem 3 and the required regularity conditions are given in the Ph.D dissertation of [48]. The proof approximates n 1/2 (μ n − µ) by the sum of independent random variables, and then applies the weak law of large numbers for i.n.i.d. random variables from Theorem 1.14 in [49] and the Lindeberg-Feller multivariate central limit theorem from Proposition 2.27 in [50]. The proof is fairly technical, but similar to those of Theorem 6.5.1 in [51], Theorem 1 in [47], and Theorem 5.1 in [3].

SE and Confidence Sets
As Section 4 has established the asymptotic theory to evaluate the variability of the proposed MLE, we can derive the SE, confidence interval (CI), and confidence ellipse (CE).
Let g : R 2 → R be a differentiable function, and g(µ) be the parameter of interest. For instance, g(µ) = µ 1 and g(µ) = µ 2 − µ 1 can be considered. The SE of g(μ n ) is This formula is based on the delta method and the large sample approximation Moreover, based on Theorem 3, we construct a 95% CE for µ: where χ 2 df=2,0.95 is be the 95% point of the χ 2 -distribution with two degrees of freedom.

R Package
We implement the proposed methods in an R package CommonMean.Copula [26]. R users can easily compute the MLE with its SE and 95% CI under the FGM, Clayton, Gumbel, Frank, and normal copulas. In this package, the log-likelihood is maximized by the R optim function, where the initial values are set as the univariate estimators For illustration, we fitted the Clayton copula by the following R codes: Symmetry 2022, 13, x FOR PEER REVIEW 10 of

R Package
We implement the proposed methods in an R package CommonMean.Copula [26]. users can easily compute the MLE with its SE and 95% CI under the FGM, Clayton, Gum bel, Frank, and normal copulas. In this package, the log-likelihood is maximized by the optim function, where the initial values are set as the univariate estimators For illustration, we fitted the Clayton copula by the following R codes: ulas by changing "Clayton" to "FGM", "Gumbel", "Frank", or "normal".

Simulation Studies
We conducted Monte Carlo simulations to examine the accuracy of the propos methods. We report the results for the Clayton copula; more results are available fro [48].

R Package
We implement the proposed methods in an R package CommonMean.Copula [ users can easily compute the MLE with its SE and 95% CI under the FGM, Clayton, bel, Frank, and normal copulas. In this package, the log-likelihood is maximized by optim function, where the initial values are set as the univariate estimators For illustration, we fitted the Clayton copula by the following R codes: One can fit othe ulas by changing "Clayton" to "FGM", "Gumbel", "Frank", or "normal".

Simulation Studies
We conducted Monte Carlo simulations to examine the accuracy of the pro ) = −285.65. One can fit other copulas by changing "Clayton" to "FGM", "Gumbel", "Frank", or "normal".

Simulation Studies
We conducted Monte Carlo simulations to examine the accuracy of the proposed methods. We report the results for the Clayton copula; more results are available from [48].
We generated Y i , i = 1, 2, . . . , n, under the Clayton copula with α i ∼ Gamma(64, 1/8), In all three cases, we have Var[α i ] = 1. Without loss of generality, we set µ = (0, 0). To set σ 2 i1 and σ 2 i2 , we followed the simulation setting of [52]. That is, , andμ Clayton n , and their SEs and 95% CIs (CEs) by using the R function CommonMean.Copula (Section 6). We then evaluated the coverage probability (CP) of the 95% CI (CE) to see how the confidence set can cover the true value. We consider a small sample size n ∈ {5, 10, 15} and a large sample size n ∈ {50, 100, 300}. Our simulations are based on 1000 repetitions. , the SDs of the estimates decrease when n increases from n = 5 to n = 300. We report the boxplots summarizing the 1000 repetitions forμ Clayton n,1 in Figure 2. This clearly visualizes how the variability of the estimates vanishes as the sample sizes increase. Table 1 also shows that the SDs are close to the average SEs, except for n = 5 (due to the very small samples). Consequently, the CPs are close enough to the nominal level of 0.95, especially when sample sizes are large, which is consistent with our asymptotic theories. Forμ Clayton n , the CPs of the 95% CEs are also reasonably close to the nominal level. In summary, the proposed estimators and the asymptotic theory work fairly well in finite samples.

Data Analysis
We analyze two real datasets to illustrate the usefulness of the proposed methods.

The Entrance Exam Data
The first dataset we analyzed was the entrance exam scores on mathematics and statistics, which was introduced by [3]. The data come from undergrad students who took written exams from 2013 to 2017 to enter the Graduate Institute of Statistics, National Central University, Taiwan. The possible score range is from 0 to 100 for both subjects. Let i = 1, 2, . . . , 5 be indices for years 2013, 2014, . . . , 2017. Table 2 provides the data, including the values of mathematics (Y i1 = mean math score) and statistics (Y i2 = mean stat score), and their covariance matrix (Ω i ). We fitted the data to the proposed model using the R function CommonMean.Copula(.) in our R package (Section 6). Table 3 summarizes the fitted results for the FGM, Clayton, Gumbel, Frank, and normal copulas. According to the values of the log-likelihood, the Gumbel copula produces the best fit, the Frank copula the second best, and the bivariate normal model the worst fit. The FGM copula failed to capture the dependence and fitted at the boundary θ i = 1 for all i ( Table 2). Since the number of unknown parameters across different copulas is the same, model selection by the Akaike information criterion (AIC) is equivalent to model selection by the log-likelihood value. An alternative way of selecting a copula is based on a leave-one-out cross validation (CV), defined as n,2 are the MLE obtained without the ith sample. Here, CV measures how a sample is predicted by the others under a copula model. A smaller CV corresponds to a better performance of the model. Table 3 reports the values of CV for each copula. It shows that the Clayton copula has the best performance while the Gumbel copula has the worst. The normal copula has the second worst performance. Overall, our analysis clearly shows the insufficiency of the bivariate normal model. Figure 3 shows the 95% CEs for the mean vector µ. This visualizes how the resultant estimates vary from the choice of copulas. Interestingly, the CE under the Clayton copula is far away from the other four, although it has a larger log-likelihood value than the normal copula. The normal and Clayton copulas produce the rotated oval shape of the CEs, representing a positive dependence between math and stat scores. The FGM and Gumbel copulas produce similar shapes for their CE. We adopt the 95% CE given by the Gumbel copula because it has the largest log-likelihood value. Figure 4 gives the 3D plot of the log-likelihood surface under the Gumbel copula model. The plot shows that the estimate of the common meanμ Gumbel n = (37.67, 42.56) attains the global maximum of the log-likelihood function. Figure 4 gives the 3D plot of the log-likelihood surface model. The plot shows that the estimate of the common m attains the global maximum of the log-likelihood function.   attains the global maximum of the log-likelihood function.

The Blood Pressure Data
The second dataset we used contains 10 studies that examined the effectiveness of hypertension treatment for lowering blood pressure. Each study provides complete data on two treatment effects, the difference in systolic blood pressure (SBP) and diastolic blood pressure (DBP) between the treatment and the control groups, where these differences are adjusted for the participants' baseline blood pressures. The within-study correlations of the two outcomes range from ρ i = 0.45 to ρ i = 0.78, exhibiting positive dependence. This dataset is available in R package mvmeta [30] and was previously analyzed by [53].
We fitted the data to the proposed copula models using the R function Common-Mean.Copula(.) in our R package (Section 6). Table 4 summarizes the fitted results for all the copulas. Based on the log-likelihood values, the Frank copula produces the best fit, the Gumbel copula the second best, and the Clayton copula produces the worst fit. The FGM copula failed to capture the dependence and fitted at the boundary θ i = 1 for all i. Again, our analysis reveals the insufficiency of the bivariate normal model; the Frank copula best captured the correlations in the blood pressure data. We also compared CV across all the copulas ( Table 4). The results show that the Clayton copula has the best performance while the normal copula has the worst. Again, our analysis shows the insufficiency of the normal model.  Figure 5 shows the 95% CEs for the mean vector µ. The CE under the Clayton and normal copula are far away from the other three. The CE under the FGM copula was almost fully covered by the CE under the Frank copula. We adopt the 95% CE given by the Frank copula, since it has the largest log-likelihood value ( Table 4).
The second dataset we used contains 10 studies that examined the effectiv hypertension treatment for lowering blood pressure. Each study provides compl on two treatment effects, the difference in systolic blood pressure (SBP) and d blood pressure (DBP) between the treatment and the control groups, where thes ences are adjusted for the participants' baseline blood pressures. The within-stud lations of the two outcomes range from = 0.45 to = 0.78, exhibiting posi pendence. This dataset is available in R package mvmeta [30] and was previously a by [53]. We fitted the data to the proposed copula models using the R function C Mean.Copula(.) in our R package (Section 6). Table 4 summarizes the fitted result the copulas. Based on the log-likelihood values, the Frank copula produces the bes Gumbel copula the second best, and the Clayton copula produces the worst fit. T copula failed to capture the dependence and fitted at the boundary = 1 fo Again, our analysis reveals the insufficiency of the bivariate normal model; the Fra ula best captured the correlations in the blood pressure data. We also compa across all the copulas ( Table 4). The results show that the Clayton copula has performance while the normal copula has the worst. Again, our analysis shows th ficiency of the normal model.  Figure 5 shows the 95% CEs for the mean vector . The CE under the Clay normal copula are far away from the other three. The CE under the FGM copula most fully covered by the CE under the Frank copula. We adopt the 95% CE give Frank copula, since it has the largest log-likelihood value (Table 4).

Extension to Non-Normal Models
So far, we have considered a common mean model under the marginal nor This section explains how the proposed methods can be extended to non-normal m For this reason, we specifically consider a common mean model under the margi ponential distributions. Thus, the common mean vector is = (1/ , 1/ ). We consider the Clayton copula to specify the bivariate distribution becaus simple derivatives with respect to the copula parameter [54]. Therefore, we propo variate common mean Clayton copula model with exponential margins as follows where is known for = 1,2, … , . Note that copula Clayton is a survival cop ( , ) as the usual way to model a survival function [22]. Using similar argum [55], the information matrix with respect to = ( , ) can be decomposed as See Appendix A.3. for detailed derivations. The expression of ( ) is an ex of Theorem 3 to the exponential model. With the information matrix, the propertie

Extension to Non-Normal Models
So far, we have considered a common mean model under the marginal normality. This section explains how the proposed methods can be extended to non-normal models. For this reason, we specifically consider a common mean model under the marginal exponential distributions. _ F j (y) ≡ Pr(Y ij > y) = exp(−λ j y), y > 0, λ j > 0, j = 1, 2.
Thus, the common mean vector is µ = (1/λ 1 , 1/λ 2 ). We consider the Clayton copula to specify the bivariate distribution because it has simple derivatives with respect to the copula parameter [54]. Therefore, we propose a bivariate common mean Clayton copula model with exponential margins as follows: where α i is known for i = 1, 2, . . . , n. Note that copula C Clayton α i is a survival copula for (Y i1 , Y i2 ) as the usual way to model a survival function [22]. Using similar arguments to [55], the information matrix with respect to λ = (λ 1 , λ 2 ) can be decomposed as where φ(α) = 1 3α + 1 See Appendix A.3. for detailed derivations. The expression of I i (λ) is an extension of Theorem 3 to the exponential model. With the information matrix, the properties of the MLE and the asymptotic theory are similar to the normal models.
We conducted Monte Carlo simulations to examine the correctness of Equation (5) by comparing it with their empirical version. We set λ 1 = λ 2 = 1 and α i = 1 for all i. We generated data (Y i1 , Y i2 ), i = 1, . . . , n from the model in Equation (4) and computed the empirical versions of I i,11 (λ) and I i,12 (λ) as The formulas for the derivatives of the log-density are found in Equations (A1) and (A2) in Appendix A.3. Our simulations were based on 1000 repetitions with n ∈ {100, 200, 300, 400, 500}.  We generated data ( , ), = 1, … , from the model in Equation (4)

Conclusions
At present, copula models are very popular in all areas of science. Bivariate metaanalyses are among those research areas that require sophisticated copula-based methods and theories. Nonetheless, there are only a few studies on copula-based bivariate metaanalysis from a methodological/theoretical perspective. This article fully develops the methodologies and theories of the copula-based bivariate meta-analysis, specifically for estimating the common mean vector. These developments will provide solid methodological/theoretical bases that are not available to date.
In this article, we emphasize the flexibility of the proposed copula models that allow for a variety of dependence structures. In the two real data examples, we employed the log-likelihood value as a criterion for model selection (Section 8). Even if the best copula is selected, it still raises the issue of goodness-of-fit, which is difficult to assess under the meta-analysis setting. The classical methods, such as Kolmogorov-Smirnov or Cramérvon Mises type statistics, cannot be directly applied to the non-identically distributed samples for which the empirical distribution function is difficult to interpret. Therefore, the development of goodness-of-fit tests is a possible research direction.
The fundamental assumption made in the proposed model is the common mean model, with known within-study correlations. The common mean assumption, although convenient for summarizing the data for a small number of studies [56], may not always

Conclusions
At present, copula models are very popular in all areas of science. Bivariate metaanalyses are among those research areas that require sophisticated copula-based methods and theories. Nonetheless, there are only a few studies on copula-based bivariate metaanalysis from a methodological/theoretical perspective. This article fully develops the methodologies and theories of the copula-based bivariate meta-analysis, specifically for estimating the common mean vector. These developments will provide solid methodological/theoretical bases that are not available to date.
In this article, we emphasize the flexibility of the proposed copula models that allow for a variety of dependence structures. In the two real data examples, we employed the log-likelihood value as a criterion for model selection (Section 8). Even if the best copula is selected, it still raises the issue of goodness-of-fit, which is difficult to assess under the meta-analysis setting. The classical methods, such as Kolmogorov-Smirnov or Cramér-von Mises type statistics, cannot be directly applied to the non-identically distributed samples for which the empirical distribution function is difficult to interpret. Therefore, the development of goodness-of-fit tests is a possible research direction.
The fundamental assumption made in the proposed model is the common mean model, with known within-study correlations. The common mean assumption, although convenient for summarizing the data for a small number of studies [56], may not always hold in real meta-analyses [6]. Therefore, the extension of the proposed estimator to random means (random-effects models) or ordered means [57,58] is an important direction for future research. To model the random effects, we need another bivariate copula. The estimation problem for these hierarchical copula-based models is beyond the scope of the paper. Nonetheless, the results presented in this paper serve as fundamental knowledge before the exploration of more advanced models. We first prepare a lemma: Lemma A1. Under the general copula model (Definition 1), if C [2,2] θ exists in (0, 1) 2 , the correlation function has alternative expressions where Z i1 and Z i2 have the joint density Proof of Lemma A1. We only prove the first identity for illustration. If C [2,2] θ exists, then where the last equality follows from Stein's identity. Thus, we obtain The proof completes.
Lemma A1 is a generalization of Lemma 3.2 in [3]. Now, we prove Lemma 1 for j = 1 and k = 2. If C [2,2] θ exists, by straightforward calculations, On the other hand, Based on the above results, it suffices to show which is asserted by Lemma A1. Hence, the proof is completed.

Appendix A.2. Derivatives for Copulas
The normal copula: The FGM copula:
The first-order partial derivative of log f y j e α i λ j y j e α i λ 1 y 1 + e α i λ 2 y 2 − 1 , j = 1, 2.
For the inner integral, ∞ 0 x 2 e 2α i x 2 (e α i x 1 + e α i x 2 − 1) 1/α i +4 dx 1 where the second equality follows from integration by parts. For the first integral, (e α i x 1 + e α i x 2 − 1) 1/α i +3 dx 2 = − 1 (3α i + 1)(2α i + 1) 1 (e α i x 1 + e α i x 2 − 1) 1/α i +2 ∞ 0 = 1 (3α i + 1)(2α i + 1) e −2(α i +1)x 1 . The second integral is where the second equality follows from integration by parts. Now, the expectation becomes For the integral in the above expression, we consider its inner integral, (e α i x 1 + e α i x 2 − 1) 1/α i +1 d(x 1 e α i x 1 ) where the second last equality follows from integration by parts. We compute the above two integrals separately. We have On the other hand, we have Then, Combine all the results, one has Let α i x 1 = s and α i x 2 = t, according to [52], we have Hence, we obtain Finally, combining the above results, we have