Bayesian multivariate meta-analysis of multiple factors

1 INTRODUCTION

As a disease condition is typically associated with many risk and protective factors in medical sciences, many randomized controlled trials and observational studies considered multiple factors.1-4 Reliable summary of association between each factor and the disease condition is crucial for the design of a multi-factor intervention program. The growth of interest in evidence-based medicine has led to a dramatic increase in attention paid to systematic reviews and meta-analyses. In prevention studies, it has become increasingly popular to perform meta-analyses on multiple risk and protective factors to summarize existing evidence; however, currently, nearly all meta-analyses are performed on each factor separately.5-7 Different studies usually focus on different subsets of all risk and protective factors and may only selectively report some significant factors in peer-reviewed articles; some factors may be reported by only a few studies. Hence, if we organize the collected data in a matrix with rows and columns indexing studies and factors, respectively, then the data matrix is expected to contain many missing entries; see the example in Table 1. The conventional meta-analysis separately estimates each factor's association with the disease condition, so it cannot borrow information across the factors. Moreover, results from separate meta-analyses may not be directly comparable because they may be based on different subpopulations. This limits medical investigators to select most important factors for the design of a multi-factor intervention program.

Recently, Serghiou et al8 introduced field-wide systematic review and meta-analysis to report and assess the entire field of putative factors for a disease condition. Based on this concept, researchers can learn the selective availability and different adjustments of factors and the patterns of modeling. Although multiple factors were collected, the authors pooled the results for each factor separately; this is not efficient to analyze the multivariate data from a field-wide systematic review.

This article proposes multivariate meta-analysis of multiple factors (MVMA-MF) to jointly synthesize all risk and protective factors in the field-wide systematic review. Compared with separate meta-analyses that use different subpopulations, the results produced by jointly modeling of multiple factors have better interpretation and generalizability because they are based on the whole population in the multivariate meta-analysis. Multivariate meta-analysis methods have gained much attention in the recent literature,9-12 and they improve effect estimates by borrowing information from the correlations between multiple endpoints.13-15 The multivariate meta-analysis models have been applied to several areas, such as meta-analysis of diagnostic tests,16-18 meta-analysis of multiple outcomes,19, 20 and network meta-analysis of mixed treatment comparisons.21-26 Mixed treatment comparisons use both direct and indirect evidence of treatment contrasts to synthesize the comparisons between multiple treatments; its focus differs from MVMA-MF, because MVMA-MF is interested in estimating the effect of multiple factors, not the contrasts between them.

Multivariate random-effects model generally requires correlations within each collected study. In some situations, within-study correlations are known to be 0. For example, in meta-analysis of diagnostic tests, the study-specific sensitivity and specificity are uncorrelated within studies because they are calculated from the true-negative and true-positive patients, respectively.10 However, in MVMA-MF, the factors can be correlated within each study because they may be measured on same patients. Such within-study correlations are unknown unless individual patient data are available. Ignoring within-study correlations in the standard multivariate random-effects model may have great impact on the estimated overall effect sizes.27

To deal with unknown within-study correlations, this article considers an alternative Bayesian model for MVMA-MF. The conventional multivariate model partitions the overall variance-covariance matrix into 2 parts: the within-study level that is due to sampling error and the between-study level that is due to heterogeneity between the collected studies. Instead of partitioning the overall correlations into the 2 levels, the alternative model directly specifies 1 single overall correlation matrix; hence, it may be viewed as a hybrid approach. This model is the Bayesian version of the model introduced by Riley et al.28 Currently, Riley's model is implemented in a frequentist way, such as using the restricted maximum likelihood method.29 However, the data of MVMA-MF are usually fairly sparse (Table 1), and our simulations in Section 5 show that the frequentist method generated poor 95% confidence intervals (CIs) for sparse data; also, the algorithm for maximizing the (restricted) likelihood did not converge for many simulated data. Instead of using the frequentist method, we present a fully Bayesian method to perform MVMA-MF. Although Bayesian methods have been introduced to deal with multivariate meta-analyses by Lu and Ades30 and Wei and Higgins,31 the former paper considered mixed treatment comparisons, which focused on treatment contrasts, and the latter assumed known within-study covariance matrices. Hence, both Bayesian methods may not be directly applied to MVMA-MF. The Bayesian version of Riley's model will be detailed in Section 4. We will use simulations and reanalyze the real data in Serghiou et al8 to demonstrate the benefit of the joint Bayesian model for MVMA-MF, and we will conclude this article with brief discussion.

2 THE MOTIVATING PTERYGIUM DATA

Pterygium is one of the most common eye conditions, and it may be related to many risk factors.32 Instead of reporting only 1 risk factor at a time, Serghiou et al8 collected the odds ratios of all putative risk factors for pterygium. Specifically, they identified 60 eligible studies reporting on a total of 65 risk factors. Since most risk factors were only reported in less than 3 studies, we focus on the following 8 risk factors, each of which was reported in at least 4 studies: (1) occupation type (outdoor vs indoor); (2) smoking status (yes vs no); (3) education attainment (low vs high); (4) use of hat (yes vs no); (5) use of spectacles (yes vs no); (6) area of residence (rural vs urban); (7) use of sunglasses (yes vs no); and (8) latitude of residence (low vs high). These risk factors are sorted according to their frequencies reported in the collected studies (from high to low). Also, we cleaned the data by removing the log odds ratios that were obtained using multivariate regression model, because they were adjusted for different risk factors in different studies. Table 1 presents the cleaned data, and Figure 1 shows the network plot of the 8 risk factors. The network indicates that most pairs of risk factors are simultaneously reported in some studies, but several pairs, such as area of residence and use of hat, are not. From Table 1, the risk factor latitude of residence was reported in only 4 studies, while 23 studies reported occupation type. Most studies reported different subsets of the 8 risk factors, and many entries in Table 1 are missing. The estimated overall effect sizes produced by univariate models may be poor because some risk factors have few samples. Also, many factors (eg, area of residence and education attainment) are expected to be correlated, so a multivariate model may be more appropriate for this dataset than univariate models.

3 CONVENTIONAL META-ANALYSIS MODELS

This section reviews some existing models for general multivariate meta-analysis. Suppose that n independent studies are collected; each study reports a p-dimensional vector of effect sizes, denoted as y_i = (y_i1, …, y_ip)^T, and denote its within-study variance-covariance matrix as S_i (i = 1, …, n). The conventional univariate meta-analysis pools the results for each j = 1, …, p separately; a fixed- or random-effects model is applied to the data , where v_ij is the within-study variance, ie, the jth diagonal element in S_i.33, 34 Since most studies are conducted by different research teams in different places using different methods, the studies are usually expected to be heterogeneous.35 Also, the random-effects model may produce more conservative results than the fixed-effects model,36, 37 so this article focuses on the random-effects setting that accounts for the heterogeneity between studies. We denote the univariate model as Model U, which ignores both within- and between-study correlations.

Multivariate meta-analysis recently gains much attention to simultaneously synthesize the p-dimensional effect sizes.9-11 Given that the within-study variance-covariance matrices S_i are known, the commonly used random-effects model is specified as follows to analyze the multivariate data :

where μ_i represents study i's underlying true effect sizes, μ = (μ₁, …, μ_p)^T contains the overall effect sizes, and T = (τ_ij) is the p × p between-study variance-covariance matrix. We denote this multivariate model as Model M. Note that within-study variances are routinely reported in published articles, but within-study correlations are usually unavailable. Let D_i =  diag (S_i) be the diagonal matrix consisting of the within-study variances. To analyze the data when the within-study correlations are unknown, a naïve multivariate method is to simply ignore these correlations by setting them to 0 but still account for the between-study correlations; we denote this model as Model M₀. Nevertheless, ignoring within-study correlations could lead to poor estimated effect sizes, especially when the within-study correlations are comparable with or greater than the between-study correlations.27 The following section introduces an alternative model that can incorporate both within- and between-study correlations for MVMA-MF.

4 MULTIVARIATE META-ANALYSIS OF MULTIPLE FACTORS

4.1 Multivariate hybrid meta-analysis model

Model M may be deemed ideal to perform MVMA-MF: It uses the within-study correlations between multiple factors that are usually unknown, and it provides a benchmark for the performance of other potential models. Note that study i's marginal variance-covariance matrix in Model M is M_i = S_i + T, and it can be written as M_i = (D_i + Δ)^1/2R_i(D_i + Δ)^1/2, where contains between-study variances (ie, diagonal elements in T), and D_i is the diagonal matrix containing within-study variances (ie, diagonal elements in S_i). The marginal correlation matrix of study i, R_i, is determined by both S_i and T, and thus, it needs to be estimated if the within-study correlations are unknown. It may be inefficient to use the data from merely n studies to simultaneously estimate all R_i's, which involve too many parameters.

Alternatively, extending the bivariate model in Riley et al,28 we consider a multivariate model that does not require within-study correlations to perform MVMA-MF:

where is a diagonal matrix that consists of additional variances beyond sampling error because of between-study heterogeneity for the p effect sizes. In this model, all collected studies are assumed to share a common marginal correlation matrix R. This assumption effectively reduces the number of parameters to be estimated and accounts for both within- and between-study correlations. Like Model M, the alternative model partitions the marginal variances of y_i, D_i + Ψ, into the within- and between-study levels; however, it directly uses the matrix R to model the overall correlations, instead of partitioning the correlations into the previous 2 levels. Therefore, the alternative model may be deemed hybrid, and we denote it as Model H.

4.2 Missing data

So far, the models for complete data have been discussed. In MVMA-MF, each collected study often reports a small subset of the complete factors, and many factors are missing as in Table 1. It is straightforward to extend the 4 methods (ie, Models U, M, M₀, and H) to deal with missing data under the assumption of missing at random. Model H for missing data will be detailed here. Suppose that contains the complete p factors in study i; however, we only observe t_i factors and their within-study variances (1 ≤ t_i ≤ p). Denote the effect sizes of the t_i factors as , which is a t_i-dimensional subvector of , and let D_i be the t_i × t_i diagonal matrix containing the within-study variances. We write , where is a t_i × p matrix indicating missingness.31 Specifically, for each j = 1, …, t_i, e_ij = (e_ij1, …, e_ijp)^T with e_ijk = 1 if the observed y_ij is the effect size of factor k and e_ijk = 0 otherwise. For example, for study 3 in the pterygium data (Table 1),

Recall that μ, Ψ, and R represent the overall mean effect sizes, the between-study variances, and the marginal correlation matrix for the complete p factors, respectively. Hence, for the observed t_i-dimensional vector y_i, its overall mean is X_iμ, its marginal variances form the diagonal matrix , and its marginal correlation matrix is . Consequently, the hybrid model for missing data can be specified as

Although the derivation here is based on the assumption of missing at random,38 the simulations in Section 5 will compare the performance of Model H with Models M, M₀, and U when some factors are missing under various mechanisms (including missing not at random).

5 SIMULATIONS

We conducted simulations under various settings to compare the performance of the hybrid model (Model H) with the ideal model that uses within-study correlations (Model M), the model that ignores within-study correlations but accounts for between-study correlations (Model M₀), and the univariate model that ignores both types of correlations (Model U). Bias and root mean squared error (RMSE) of point estimate and 95% CI/CrI coverage probability were used to evaluate the models' performance. We set the number of studies in each simulated MVMA-MF dataset to 30 and considered 5 factors in total. The true overall effect sizes of the 5 factors were 0, ie, μ = (0, 0, 0, 0, 0)^T. Also, the between-study standard deviation τ was 1 for each factor; the within-study standard deviation of each factor σ was set to 0.5, 1, or 2. These choices for σ represented different extents of heterogeneity between studies. Moreover, we considered the exchangeable correlation structure for both the between- and within-study correlation matrices, R_B = (r_Bij) and R_W = (r_Wij), which were determined by the correlation parameters ρ_B and ρ_W, respectively; that is, r_Bij = ρ_B and r_Wij = ρ_W for 1 ≤ i ≠ j ≤ 5. The between-study correlations ρ_B were 0.2, 0.5, or 0.8, and the within-study correlations ρ_W were drawn from U(0, 0.3), U(0.3, 0.6), or U(0.6, 0.9). For each setting, 1000 MVMA-MF datasets were simulated using the ideal Model M with S_i = σ²R_W and T = τ²R_B. Three missingness scenarios were considered: (I) all 5 factors were observed in all studies, ie, the data were complete; (II) the data of factors 1, 3, and 5 in 10 studies were missing completely at random; and (III) the smallest 10 effect sizes of factors 1, 3, and 5 were missing. The missingness that was not at random in scenario (III) may be viewed as the effect of reporting bias. We also considered a missingness scenario that was similar to (III) but contained more missing data: (III′) the smallest 25 effect sizes of factors 1, 3, and 5 were missing; in this case, the 3 factors were only available from 5 studies, so the simulated MVMA-MF data were much sparser than the previous scenarios. Both the restricted maximum likelihood (REML) and Bayesian methods were applied to implement the 4 models. For the Bayesian analyses, we used the priors specified in Section S1.2. The Markov chain Monte Carlo (MCMC) algorithm was used for estimation with 1 chain, having a run of 10 000 updates after a 10 000-run burn-in period.

Table 2 presents the results when ρ_B was 0.5, and Tables S1 and S2 present those when ρ_B was 0.2 and 0.8, respectively. First, since the data of factors 2 and 4 were complete under all scenarios, their corresponding results produced by the 4 models were almost the same; all models led to nearly unbiased estimates and proper 95% CI/CrI coverage probabilities for these 2 factors. Second, the results of factors 1, 3, and 5 produced by the 4 models also differed little when the data were complete or missing completely at random. This is expected from the perspective of missing data analysis.38 Third, if the missingness was not at random, the results produced by the 4 models were noticeably different. The univariate model led to the largest bias and RMSE and the lowest 95% CI/CrI coverage probability. Since Model M₀ still accounted for the between-study correlations, its performance was similar to Model H when ρ_W was small compared with ρ_B. However, Model H outperformed Model M₀ when ρ_W was larger than ρ_B, because the within-study level dominated the estimation of the overall effect sizes in this situation, but Model M₀ ignored correlations at this level. Finally, note that Model M was ideal as it used the within-study correlations that are usually unavailable in real world. Although Model H did not use the within-study correlations, the biases and RMSEs produced by Model H were fairly close to the ideal Model M under various settings. Also, the 95% CI/CrI coverage probability produced by Model H was generally higher than those produced by Models M₀ and U.

In most situations, the biases and RMSEs of point estimates obtained using the REML method were close to those obtained using the Bayesian method. However, under the missingness scenarios (III) and (III′), the 95% CrI coverage probabilities for factors 1, 3, and 5 obtained using the Bayesian method were generally higher than those obtained using the REML method. As noted in Section 4.3, this may be due to that the estimated variance-covariance matrix was inconsistent when many observations were missing and the dimension was close to the sample size (only 5 observations for each of factors 1, 3, and 5). Also, under scenario (III′), the biases and RMSEs produced by Model H using the REML method were noticeably different from those using the Bayesian method for factors 1, 3, and 5; again, these differences may be caused by the sparsity of the simulated data. In addition, when Model H was used under scenario (III′), the optimization algorithm for the REML estimates did not converge for many simulated replicates. We ran around 2500 iterations to obtain 1000 datasets that produced converged REML estimates, and these 1000 datasets were used to produce the results for scenario (III′). Hence, the Bayesian method is preferred to implement the multivariate hybrid model when the dimension of factors is high but the observations are limited.

Moreover, we conducted an additional simulation study driven by the pterygium data. The conclusions are similar to those presented here. For more details, see Section S2.

6 REAL-DATA ANALYSIS

The pterygium dataset has been introduced in details as our motivating example. Since the within-study correlations were unknown, we used Models H, M₀, and U but not Model M to estimate the overall log odds ratios of the 8 risk factors. Because of the sparsity of the dataset, the models were implemented using Bayesian analyses. The estimates were based on 3 MCMC chains, each containing a run of 100 000 updates after a 100 000-run burn-in period. The convergence of the chains was checked using their trace plots. Table 3 presents the median overall log odds ratios with 95% CrIs. Figure 2 shows the posterior density plots of the 8 risk factors.

For risk factors that were reported in a relatively large number of studies (eg, occupation type, smoking status, and education attainment), the 3 models yielded similar estimated overall log odds ratios; their density curves were also fairly similar. For risk factors that were only reported in a few studies (eg, use of spectacles, area of residence, and latitude of residence), the peaks of the density curves produced by Model H were narrower and higher compared with those produced by Models M₀ and U, indicating that Model H produced narrower 95% CrIs. Also, for the risk factor use of sunglasses, the location of its density curve produced by Model H was noticeably different from those produced by the other 2 models. Figure 3 depicts the estimated overall correlations between the 8 risk factors produced by Model H. It shows that many factors were correlated and some correlations were fairly high. Hence, ignoring the within-study correlations could lead to fairly different estimated overall log odds ratios; the unknown within-study correlations need to be carefully considered in MVMA-MF.

The priors in the above real-data analyses were specified in Section S1.2; specifically, uniform priors U(0, 10) were assigned for the between-study standard deviations. We examined the results' sensitivity to prior specification by changing U(0, 10) priors to U(0, 3), U(0, 100), and U(0, 1000). In addition, we altered the order of the 8 risk factors in the data analysis; this basically changed the priors assigned to the correlation matrices in Models H and M₀ (see Section S1). Table S3 presents the results. In sum, changing priors had little impact on the risk factors occupation type, smoking status, and education attainment, which were reported in at least 10 studies. When the priors for the between-study standard deviations varied from U(0, 3) to U(0, 100), the point estimates of the overall log odds ratios of all 8 risk factors remained nearly unchanged. The CrIs for risk factors use of sunglasses and latitude of residence noticeably became wider. Compared with other risk factors, since the data for these 2 factors were limited, the priors had larger influence on their posteriors. When the priors for the between-study standard deviations varied from U(0, 100) to U(0, 1000), most results did not change. Changing the priors for the correlation matrices in Models H and M₀ had some impact on the CrIs for the risk factors that were reported in no more than 6 studies. Nevertheless, the CrIs produced by Model H were still narrower than those produced by Models M₀ and U. For more details, see Section S3.

We also performed additional sensitivity analyses to investigate the impact of risk factor selection on the estimated overall log odds ratios. We considered 2 scenarios for the set of risk factors to be included in MVMA-MF: (i) subdataset consisted of the k most frequently reported risk factors and (ii) subdataset consisted of the k least frequently reported risk factors (k = 2, …, 8). For example, when k = 2, the subdataset under scenario (i) contained the risk factors (1) occupation type and (2) smoking status; the subdataset under scenario (ii) contained the risk factors (7) use of sunglasses and (8) latitude of residence. Model H was implemented to analyze these subdatasets.

Figure 4 shows the 95% CrIs of the overall log odds ratios under both scenarios; the labels of risk factors used in this figure can be found in Table 3's first column. In scenario (i), starting from the 2 most frequently reported risk factors, infrequently reported risk factors were iteratively added to the MVMA-MF. Figure 4A shows that the estimated overall log odds ratios of risk factors 4 and 5 had some changes as new factors were added to the subdatasets. The 95% CrIs of the 3 most frequently reported risk factors 1–3 changed little. This might be explained by 2 reasons. First, the correlations between these 3 factors were weak (Figure 3), so the addition of risk factor 3 had little impact on estimating the effect sizes of factors 1 and 2. Second, the later added factors 4–8 provided much fewer information (no more than 6 studies) compared with factors 1–3, which were reported in more than 10 studies, so the correlations contributed little to the effect size estimation for factors 1–3.

Compared with scenario (i), Figure 4B shows larger changes of estimated overall log odds ratios in scenario (ii). Under this scenario, starting from the 2 least frequently reported risk factors, more frequently reported factors were iteratively added to the MVMA-MF. The 95% CrIs of infrequently reported risk factors, such as 5 and 6, became narrower as more reported risk factors were included in the MVMA-MF. This illustrates the benefit of jointly modeling multiple factors: The inference on infrequently reported factors can be strengthened by borrowing information from frequently reported factors through the correlations between them.

7 DISCUSSION

This article has introduced MVMA-MF with application to the pterygium data. In contrast to the tradition of meta-analyzing each single factor separately, this article encourages researchers to collect all possible factors and analyze them jointly to enhance the estimation of overall effect sizes. A Bayesian multivariate hybrid model has been introduced to implement MVMA-MF in which within-study correlations are usually unknown. The simulations have shown that, when many factors are missing not at random, Bayesian MVMA-MF can use the information across multiple factors to impute missing data and thus can produce better estimates than separate meta-analyses. Our results are consistent with those produced by Kirkham et al14 using frequentist bivariate and univariate meta-analyses; that is, “borrowing of strength” from the correlations between factors can reduce bias under the setting of missing not at random. For example, consider 2 highly correlated factors with a correlation coefficient close to 1, and the results for factor 1 are fully observed. Even if the results for factor 2 are missing not at random, the missing results can be largely recovered through the high correlation between the 2 factors, so the bias in the estimate of factor 2 may be greatly reduced.

The hybrid model effectively reduces model complexity by specifying a common marginal correlation matrix for all studies. This model could be a proper tradeoff between goodness-of-fit and model complexity in some situations; however, if the collected studies' marginal correlation matrices differ a lot, the hybrid model may fit the data poorly. Future work is needed to generalize the hybrid model by allowing more degrees of freedom and select the best model using deviance information criterion.40

An important issue of MVMA-MF is to incorporate the effect size of a certain factor that has been adjusted for other factors. For example, in the original pterygium data presented in Serghiou et al,8 many collected studies only reported log odds ratios that were obtained using multivariate regression after adjusting different factors (eg, age and gender), while log odds ratios without any adjustments were unavailable from these studies. We did not include such data in Table 1 because of the inconsistent adjustments. How to incorporate such data with different adjustments is of great interest to enrich the MVMA-MF and enhance the robustness and precision of the results.41, 42 We leave this to future studies.

Another interesting but challenging problem is to identify the factors' missingness mechanism and impute the missing factors more properly when the missingness is not at random. The factors could be missing completely at random, missing at random, or missing not at random.43 If the missingness of a factor does not depend on any observed or unobserved data, then this factor is missing completely at random. It is also possible that this factor is missing because of some observed data in some earlier studies; this is the case of missing at random. Bayesian inference under these 2 missingness mechanisms is valid when restricted to the observed data.38 However, the missingness not at random closely relates to publication and reporting bias, when some studies or factors were intended to be not published or reported because of their negative results.44, 45 Although our simulations showed that the proposed hybrid model performed better than Models M₀ and U when missingness was not at random, both its bias and 95% CrI coverage probability were expected to be further improved. Approaches to correcting publication or reporting bias have been introduced and widely used in univariate meta-analysis46; similar methods are highly needed for multivariate meta-analysis.

CONFLICTS OF INTEREST

None.