 Research article
 Open Access
 Published:
Efficient twostep multivariate random effects metaanalysis of individual participant data for longitudinal clinical trials using mixed effects models
BMC Medical Research Methodology volume 19, Article number: 33 (2019)
Abstract
Background
Mixed effects models have been widely applied in clinical trials that involve longitudinal repeated measurements, which possibly contain missing outcome data. In metaanalysis of individual participant data (IPD) based on these longitudinal studies, joint synthesis of the regression coefficient parameters can improve efficiency, especially for explorations of effect modifiers that are useful to predict the response or lack of response to particular treatments.
Methods
In this article, we provide a valid and efficient twostep method for IPD metaanalyses using the mixed effects models that adequately addresses the betweenstudies heterogeneity using random effects models. The twostep method overcomes the practical difficulties of computations and modellings of the heterogeneity in the onestep method, and enables valid inference without loss of efficiency. We also show the twostep method can effectively circumvent the modellings of the betweenstudies heterogeneity of the variancecovariance parameters and provide valid and efficient estimators for the regression coefficient parameters, which are the primary objects of interests in the longitudinal studies. In addition, these methods can be easily implemented using standard statistical packages, and enable synthesis of IPD from different sources (e.g., from different platforms of clinical trial data sharing systems).
Results
To assess the proposed method, we conducted simulation studies and also applied the method to an IPD metaanalysis of clinical trials for new generation antidepressants. Through the numerical studies, the validity and efficiency of the proposed method were demonstrated.
Conclusions
The twostep approach is an effective method for IPD metaanalyses of longitudinal clinical trials using mixed effects models. It can also effectively circumvent the modellings of the betweenstudies heterogeneity of the variancecovariance parameters, and enable efficient inferences for the regression coefficient parameters.
Background
In clinical trials for drug developments, the outcome variables often involve longitudinal repeated measurements. In these trials, the primary statistical analyses usually compare the responses from experimental treatments with a control treatment at the end of a followup period. However, during the followup period, dropouts or missing data often arise, and may seriously influence the validity and precision of the statistical inference. Conventionally, simple imputation methods such as LOCF (last observation carried forward) had been adopted to handle these problems, but these adhoc approaches have been known to provide biased results and possibly produce erroneous conclusions [1, 2]. Recently, regulatory guidelines concerning missing data problems in clinical trials have been issued [1, 3, 4], and these invalid practices are rapidly substituted by more sophisticated methods. The mixed effects model [5, 6], especially known as the MMRM (mixed effects model for repeated measures) [7], is one of these methods and has been widely adopted for primary analyses of recent clinical trials in drug developments, which allows for valid statistical inference based on the direct likelihood approach under the missing at random (MAR) assumptions [5, 6].
The mixed effects model are increasingly used in the context of evidence synthesis as well. Recent advances of clinical trials data sharing systems (e.g., ClinicalStudyDataRequest.com, the YODA project) have enabled evidence synthesis of large individual participant level datasets and explorations of effect modifiers that are useful to predict the response or lack of response to particular treatments, because the powers of treatmentbycovariate interaction tests increasingly gain using rich statistical information from multiple clinical trials [8, 9]. However, there are additional complexities in the individual participant data (IPD) metaanalysis to synthesize longitudinal trial datasets based on the mixed effects models, due to the longitudinal nature of the datasets. The most important issue is about addressing correlations among different time points in the repeated measurements and the betweenstudies heterogeneity. In the IPD metaanalyses of ordinary univariate regression models such as analysis of covariance (ANCOVA) models, the synthesizing analyses can be simply implemented by multivariate metaanalysis [10, 11], but for the mixed effects models, the betweenstudies heterogeneity should be adequately considered not only for the fixed effects coefficients, but also for the variancecovariance parameters, simultaneously. In addition, to use the statistical information among multiple trials maximally, we should consider the correlations among the estimators of model parameters to gain efficiency. In the IPD metaanalyses using mixed effects models, Ishak et al. [12] considered to aggregate only summarized data of each time point, but their methods do not use the correlation information. Also, Jones et al. [13] proposed to implement IPD metaanalyses of longitudinal data using mixed effects models, but their methods only assumed fixed effects models for the regression coefficient parameters, and the heterogeneity among the studies was not considered. In existing IPD metaanalyses, these problems have also not been addressed and fixed effects models that assume no heterogeneity have been usually used.
In this article, we provide a valid and efficient method for IPD metaanalyses using the mixed effects models that adequately address the betweenstudies heterogeneity using random effects models. We especially show the ordinary onestep methods have practical difficulties on computations and modellings of the heterogeneity, and we provide a twostep approach as an alternative effective procedure to overcome these problems. We will show that the twostep method can effectively circumvent the modellings of the betweenstudies heterogeneity of the variancecovariance parameters and provides valid and efficient estimators for the regression coefficient parameters, which are the primary objects of interests. In addition, these methods can be easily implemented using standard statistical software, and have advantages to synthesize IPD from different sources of clinical trial datasets (e.g., from different platforms of clinical trial data sharing systems). We also demonstrate the effectiveness of these methods via simulation studies and a real data application to an IPD metaanalysis of clinical trials for new generation antidepressants [14].
Methods
Mixed effects models for longitudinal data
We consider a metaanalysis of N longitudinal clinical trials with n_{i} individual participants within the ith trial (i = 1, … , N). Let Y_{ijt} be continuous repeated measurement outcome of tth time point of participant j in trial i, (j = 1, … , n_{i}; t = 1, … , T), possibly to be missing at some points. Note that the time points planned to be measured the outcome variables are uncommon across the N trials, we can formally regard these unmeasured time point data as ‘missing’ without loss of generality.
At first, we discuss linear mixed effect models [5, 6] for the analysis of individual trials separately,
where Y_{ij} = (Y_{ij1}, … , Y_{ijT})^{T} is a Tvector of repeated measurements for participant j in trial i, X_{ij} = (x_{ij1}, … , x_{ijT})^{T} is a T × p matrix of fixed effect explanatory covariates, β_{i} is a pvector of fixed regression coefficients, Z_{ij} = (z_{ij1}, … , z_{ijT})^{T} is a T × q matrix of randomeffect covariates of participant j in trial i (usually a subset of X_{ij}), and b_{ij} is a qvector of random effects, assumed to follow a multivariate normal distribution b_{ij}~MVN(0, D_{i}), where D_{i} is the q × q covariance matrix of the random effects distribution. ε_{ij} is a Tvector of random error, assumed to follow a multivariate normal distribution ε_{ij}~MVN(0, Σ_{i}), where Σ_{i} is the T × T covariance matrix of the random error distribution. The random effects b_{ij} and the error ε_{ij} are assumed to be independent, and thus the outcome variables among different participants are also independent. Y_{ij} marginally follows MVN(X_{ij}β_{i}, V_{ij}), where \( {\boldsymbol{V}}_{ij}={\boldsymbol{Z}}_{ij}\boldsymbol{D}{\boldsymbol{Z}}_{ij}^T+{\boldsymbol{\Sigma}}_i \). The model parameters can be estimated by maximum likelihood (ML) or restricted maximum likelihood (REML) methods [6]; we suppose to use the REML estimator here owing to its theoretical property [6, 15]. If missing outcome data are involved, these likelihoodbased inference methods provide valid statistical inference results under MAR assumptions [5, 6], and thus the linear mixed models have been widely adopted in longitudinal data analyses. In addition, recently, owing to the developments of regulatory guidelines concerning missing data treatments in clinical trials [3, 4], the mixed effects models have been increasingly popular. Especially, an alternative parameterized model named as MMRM [7] has been adopted for primary analyses of many clinical trials, which directly models the marginal distribution of Y_{ij}’s and parameterize the covariance matrix V_{ij} not specifying D and Σ_{i}. The primary objects of interests are usually the fixed effects parameters, and we have rather little interests in the variancecovariance parameters in these clinical trials [16, 17]. Thus, through the reparameterization to the marginal covariance matrix V_{ij}, possible model misspecifications concerning the random effects structures might be effectively circumvented in the statistical inference. The procedures of statistical analyses and their theoretical properties are identical between the two parameterized models, so we conduct the following discussions using the former model (1) in this article.
IPD metaanalyses based on the mixed effects models
Based on the mixed effects model (1), IPD metaanalysis can be implemented in a straightforward manner. Jones et al. [13] proposed a fixed effects model with an assumption the regression coefficients β_{i}’s are common across the N trials using both of onestep and twostep approaches. However, the betweenstudies heterogeneity is an essential issue to be addressed in modern evidence synthesis researches [9], so the random effects models are preferred in general [9, 18]. While, there are several crucial problems to be addressed in applying them to the mixed effects models for IPD metaanalysis:
1. When we consider the heterogeneity across studies, we should adequately address the heterogeneities for both of the fixed effects and variancecovariance parameters in the synthesizing analyses. A straightforward approach would be to model these parameters as random effects, the joint random effects model for b_{ij} and the components of V_{ij}. This strategy would also require complicated computations that cannot be implemented in standard statistical packages, and we should completely specify the statistical model without any misspecifications in order to assure the consistency of parameter estimation for fixed effects. Besides, we have little interests in the heterogeneities of variancecovariance parameters, i.e., the components of V_{ij} in the longitudinal clinical trials [16, 17].
2. In most of metaanalyses of longitudinal trials, the timepoints to be planned to collect outcome measurements are not common across studies. In the mixed effects model analyses, the time variables are often treated as factor variables (especially for the MMRM analyses [19]), and then all of the outcome variables at several time points are possibly missing (unmeasured) in some studies. In these cases, when we conduct a onestep random metaanalysis using standard statistical packages (e.g., PROC MIXED in SAS, lme4 [20] for R), the statistical model can be partially not identifiable, i.e., we cannot obtain the model estimates.
Given these problems, there are practical difficulties in implementing ordinary onestep IPD metaanalyses. Note that for the first problem, Jones et al. [13] proposed another approach to modelling the variancecovariance parameters as different parameters across studies in their fixed effects model framework, but these approaches are not implementable using standard packages. In addition, this problem has not been usually considered in many existing IPD metaanalyses, possibly due to its practical difficulties, and common variancecovariance structures have usually been assumed. Besides, the second issue can be resolved by making an original program that maximises the (restricted) likelihood function that assumes different regression functions across the studies (and possibly different random effects structures), but it requires usually large efforts in making the complicated original programs in casebycase practices.
Twostep marginal multivariate random effects analyses
To circumvent the complicated modelling problems and computational difficulties, we provide an alternative twostep approach for IPD metaanalyses based on the mixed effects models. Our proposition is to divide the model parameters into two components, fixed effects parameters and variancecovariance parameters, and to conduct synthesis analyses after obtaining studyspecific aggregated statistics separately. Because we have little interests in the random effects in these clinical trials and the mixed effects models are typically designed to assess fixed effects [16, 17], here we specifically focus on the synthesis analyses for the fixed effects parameters. As noted later, the same method can be adapted to the synthesis analyses of the variancecovariance parameters. The computations are easily implementable, but we can show the resultant inference is valid and efficient, and that its theoretical properties do not depend on any heterogeneous structures of withinstudies variancecovariance parameters.
We describe the procedure as follows.

Step1. Analyse the IPD from N longitudinal trials separately, using the mixed effects models (1), and obtain the REML estimates of the model parameters (\( {\widehat{\boldsymbol{\beta}}}_i,{\widehat{\boldsymbol{D}}}_i,{\widehat{\boldsymbol{\Sigma}}}_i \)), i = 1, 2, … , N. In addition, obtain the variancecovariance matrix estimates of \( {\widehat{\boldsymbol{\beta}}}_i \), \( {\widehat{\boldsymbol{S}}}_i\left({\widehat{\boldsymbol{\beta}}}_i\right) \) for i = 1, 2, … , N.

Step2. Conduct synthesizing analyses for the regression coefficient parameters \( \Big({\widehat{\boldsymbol{\beta}}}_i \), \( {\widehat{\boldsymbol{S}}}_i\left({\widehat{\boldsymbol{\beta}}}_i\right)\Big) \) that are the primary objects of interests in these metaanalyses, dropping statistical information about the withinstudies covariance parameters (\( {\widehat{\boldsymbol{D}}}_i,{\widehat{\boldsymbol{\Sigma}}}_i \)),
Inference concerning the model parameters (μ, Ω) can be conducted using the log likelihood or the restricted log likelihood functions [21, 22],
where \( {\boldsymbol{W}}_i={\left(\boldsymbol{\Omega} +{\widehat{\boldsymbol{S}}}_i\left({\widehat{\boldsymbol{\beta}}}_i\right)\right)}^{1} \).
The computation of Step1 is straightforwardly conducted using standard packages for linear mixed models, e.g., PROC MIXED in SAS, lme4 [20] for R. Also, Step2 can be also simply implemented using existing packages for multivariate metaanalyses, e.g., mvmeta for Stata [23, 24] and R [25].
The resultant estimators of (μ, Ω) have the following properties.
Remark A
There are no efficiency losses regardless of dropping the withinstudies variancecovariance information of (\( {\widehat{\boldsymbol{D}}}_i,{\widehat{\boldsymbol{\Sigma}}}_i \)) in the multivariate metaanalysis models (2) for inferences of (μ, Ω), i.e., the ML and REML estimators of (μ, Ω) on the multivariate metaanalysis models (2) are the most efficient estimators, and their efficiency bound is equivalent to that of the onestep estimators for the completely specified correct model.
This property is assured because the ML and REML estimators of β_{i} and (D_{i}, Σ_{i}) are orthogonal [6] in the sense of Cox and Reid [26], i.e., the nondiagonal submatrices of the information matrix of mixed model (1) are zeromatrices. Therefore, the statistical information of (\( {\widehat{\boldsymbol{D}}}_i,{\widehat{\boldsymbol{\Sigma}}}_i\Big) \) does not contribute to synthesizing β_{i}’s when we conduct onestep IPD metaanalyses under the completely specified correct model. Thus, fully efficient estimator of (μ, Ω) can be obtained via the twostep method above, even if the statistical information of (D_{i}, Σ_{i}) are dropped. This property is confirmed in the simulation studies in Section 3.
Remark B
When the heterogeneity for the variancecovariance parameters of (D_{i}, Σ_{i}) is modelled by a random effects model G(D, Σ), we can obtain the consistent and fully efficient estimator of (μ, Ω) using the twostep method under any distribution form of G(D, Σ).
Note this property is important that it assures the validity and optimality of the twostep method without modelling the random effects model G(D, Σ) explicitly, even if there exists substantial heterogeneity for the variancecovariance parameters. This property is assured because the ML and REML estimators of β_{i} and (D_{i}, Σ_{i}) are orthogonal as the same arguments with Remark A. Besides, if we adopt the onestep approach, we should completely specify the joint random effects distributions of β_{i} and (D_{i}, Σ_{i}) correctly to assure the consistency. These parametric assumptions are also not substantial objects of interests in the evaluations of treatment effects. Through the twostep approach, we can effectively circumvent the unnecessary additional complicated assumptions as well as the computational efforts in practices.
In addition, in IPD metaanalyses of longitudinal clinical trials, the timepoints of outcome measurements are possibly uncommon across studies. As mentioned regarding the second issue in the previous section, when the time variables are often treated as factor variables (e.g., for the typical MMRM analyses [19]), the onestep approach has practical and computational difficulties.
Remark C
For synthesis of longitudinal clinical trials in which timepoints of outcome measurements are possibly uncommon, we can apply the twostep approach formally via dropping the corresponding timepoints variables in Step1. Then, in Step2, the corresponding components of \( {\widehat{\boldsymbol{\beta}}}_i \) should be treated as missing outcomes variables.
The missing problems are common not only in twostep IPD metaanalysis, but in various multivariate metaanalysis, and have been extensively discussed in previous papers [21, 22]. Accordingly, in many software packages for multivariate metaanalysis [23,24,25], the standard modules can treat these partially missing outcomes adequately. Usually, these packages conduct available data analyses based on direct likelihood methods [23, 24].
In addition, in some cases, we can obtain the IPD for multiple clinical trials from different platforms of clinical trial databases. Recently, the clinical trial data sharing has been increasingly advanced [27, 28], and researchers can access the IPD of clinical trials for scientific research purposes, e.g., via ClinicalStudyDataRequest.com. However, in many platforms, the users cannot download the datasets and can only analyse on the remote accesses on the platform, in which case the IPD can only be analysed separately. In such cases, the onestep IPD metaanalyses cannot be conducted, but the twostep approaches can still be applied as valid and efficient procedure. In analysing longitudinal trial datasets, the twostep marginal synthesis method would be an effective strategy for evaluating treatment effects and interactions of the candidates of effect modifiers.
Results
Simulation studies
We conducted simulation studies to evaluate the validity and efficiency of the proposed method under practical IPD metaanalysis scenarios. In this context, validity corresponds to the accuracy of point estimation, standard error (SE) evaluation, and the coverage rate of the confidence interval (CI). In addition, the efficiency corresponds to the precision of estimation and the power of statistical test.
Simulation 1
First, we conducted simulations to evaluate the efficiency of twostep marginal multivariate metaanalysis approach. In particular, we would evaluate the equivalence of the twostep marginal method and the fully efficient onestep IPD metaanalysis. Simulation data were generated as three timepoints longitudinal data, following the mixed effects model involving twoway interactions,
where b_{ij} is a random intercept for the jth participant of ith trial following N(0,4.5^{2}) and ε_{ijt} is the random error for the tth visit of jth participant of ith trial following N(0,3.2^{2}), treat_{ij} is a binary variable indicating treatment group. time_{2, ijt} and time_{3, ijt} are dummy variables indicating timepoints of 2nd and 3rd visits, respectively (reference is the 1st visit). The parameter settings were mimicked to the Furukawa et al. [14]‘s example in Section 4. z_{ij} is a covariate variable that is a candidate to be an effect modifier, which is a measurement that is associated with response or lack of response to a particular therapy. The twoway interaction terms represent the potential usefulness of z_{ij} as an effect modifier. Following the mixed effects model (3), we generated N longitudinal clinical trial datasets, for N = 5, 10, 15, and 20. In addition, the number of participants of each trial was generated from a discrete uniform distribution on [50, 500].
In these first simulation studies, the primary purpose is to assess the validity and efficiency of the twostep approach provided in Section 2. Also, we would evaluate the equivalence of the twostep marginal method and the fully efficient onestep IPD metaanalysis. At first, to circumvent the discussions to be complicated, we assumed there was no heterogeneity for β_{i} and (D_{i}, Σ_{i}) among the N trials, which is a special case of the fixed effects models considered by Jones et al. [13]. In other words, all β_{i} (i = 1, … , N) are equal to the grand mean μ. The parameter settings of μ were also mimicked to the Furukawa et al. [14], μ_{1} = − 0.04, μ_{2} = − 0.14, μ_{3} = 1.56, μ_{4} = 0.83, μ_{5} = − 0.46, μ_{6} = 0.57, μ_{7} = 0.14, μ_{8} = 0.15, μ_{9} = 0.05, μ_{10} = − 0.06. Under these settings, we conducted 10,000 simulations, and we analyzed the simulated datasets using the onestep method by the correctly specified mixed effects model and the twostep multivariate metaanalysis approach in Section 2.
The results are presented in Table 1. We presented empirical means and SE of the regression coefficient estimates, means of the SE estimates (\( \widehat{\mathrm{SE}} \)), empirical coverage rates of the 95% CIs, and empirical powers for the tests of individual regression coefficients. The contexts of the results are too rich, so we present the ones for selected coefficients, μ_{2}, μ_{4}, μ_{5}, μ_{7}, μ_{9}, μ_{10}. The results for the remaining parameters were similar. These results consistently showed that the fully efficient onestep approach was equivalent to the marginal twostep approach under all of the scenarios. Both the mean and SE of 10,000 estimates were mostly the same, and the SE estimate could accurately estimate the actual SE of the estimates. In addition, coverage rates of 95% CIs were accurate and the power of the two approaches were equivalent. As a whole, the twostep approach was shown to provide accurate results without loss of efficiency, regardless of dropping the variancecovariance information. The efficiency was also retained for the interaction tests.
Simulation 2
Second, we conducted other simulations to evaluate the validity of the twostep random effects multivariate metaanalysis approaches under substantial heterogeneity of β_{i} as well as (D_{i}, Σ_{i}), among N longitudinal trials. We considered the mixed effects models (3) for individual trials, and additionally considered heterogeneity of β_{i} and (D_{i}, Σ_{i}). We considered a compound symmetry structure for Ω, which assumes a common variance among all of the components of β_{i}, β_{ik}~N(μ_{k}, τ^{2}), k = 1, … , 10, where τ^{2} = 0.10, 0.20. We assumed no correlations among β_{ik}s, and also assumed heterogeneity among the random effects distributions of b_{ij}~N(0, ζ^{2}), where ζ is generated from a continuous uniform distribution on [2.5, 9.0]. We considered two settings: (a) outcome variables are measured on all three timepoints, and (b) outcome variables are measured on two timepoints for 40% of trials and on three timepoints for the remaining 60% trials.
Under these conditions, the twostep approach can be easily implemented using standard statistical packages. Again, note that the twostep method does not require modelling of the heterogeneity structures of (D_{i}, Σ_{i}), and the validity and efficiency of the inference are assured under any heterogeneous structures of (D_{i}, Σ_{i}). Besides, to implement valid inference using onestep approach, we should adopt the correct statistical model involving the random effects distributions for both of β_{i} and (D_{i}, Σ_{i}). One adhoc approach uses a Bayesian hierarchical modelling involving random effects models, but it requires complicated computations using special software (e.g., OpenBUGS [29]).
The simulation results are presented in Tables 2 and 3, which correspond to the settings (a) and (b) respectively. We conducted 10,000 simulations and provided the same outputs as Table 1. We considered 6 scenarios, setting N = 5, 10, or 15 and τ^{2} = 0.10 or 0.20. The results also showed unbiasedness of the twostep estimates under all of the scenarios. In addition, the validities of SE estimation and CIs (inversely, the sizes of the corresponding statistical tests) could be confirmed. Again, note that we only took parametric assumptions for the random effects distribution of β_{i} and no restrictive assumptions concerning the heterogeneity of (D_{i}, Σ_{i}). The validity of the twostep method was demonstrated under the flexible and weak assumptions. The efficiency is not compared with those of the onestep methods here, but the equivalence is theoretically assured (Remark B).
Applications to IPD metaanalysis for new generation antidepressants
To illustrate our method, we analyzed the dataset from an IPD metaanalysis for new generation antidepressants by Furukawa et al. [14] The authors conducted IPD metaanalysis from four placebocontrolled, doubleblinded randomized clinical trials (1482 participants) for patients undergoing acute phase treatment for major depression. The outcome variable was the change score on the 17item Hamilton Rating Scale for Depression (HRSD) or the MontgomeryAsberg Depression Rating Scale (MADRS). In their analyses, the latter was converted into the former using the conversion algorithm based on the item response theory [30]. In particular, the authors investigated whether the baseline depression severity modifies the efficacy of the antidepressants using IPD metaanalysis. The outcome measurements were repeatedly measured at several timepoints in the four trials, but the timepoints at which the outcome variables were planned to be measured were not common across the four trials.
We considered the following mixed effects model on 1st, 2nd, 4th and 6th weeks involving twoway interaction terms based on the model (1),
where week_{1, ij}, week_{2, ij} and week_{4, ij} are dummy variables for the time points (reference: 6th week). We first attempted a onestep random effects metaanalysis that modelled all of the regression coefficients as random effects among the 4 trials, but the REML estimate was not computable using standard statistical packages. There might also be betweenstudies heterogeneity for the variancecovariance parameters, which was not modelled in this analysis. Thus, we provided the results of a onestep fixed effects metaanalysis only consider correlations of repeated measured outcomes among different time points within individual participants. The REML estimates, 95% CI, and Pvalues are presented in Table 4.
Next, we conducted IPD metaanalyses using the twostep approach based on the REML estimates of the fixed effects coefficients of the 4 trials. In the multivariate random effects metaanalysis models, we considered several assumptions for betweenstudies variancecovariance structures Ω, and the fitting of the candidate models was evaluated using the Bayesian information criterion (BIC) [31]. The best fitting model was the compound symmetry model (BIC = 125.73), whereas for the fixed effects model, BIC = 316.89.
In Table 3, we also provided the results of twostep multivariate metaanalysis of the fixed effects model and the random effects model using the compound symmetry model. In the random effects model, the betweenstudies SD estimate was 0.047 and the betweenstudies correlation coefficient was 0.064. The results for the fixed effects model analyses using onestep and twostep approaches were consistent, which was expected as described in Remark A. However, note that even if there was betweenstudies heterogeneity for the variancecovariance parameters, the twostep approach provides valid inference results. Besides, comparing the fixed and random effects twostep analyses, the coefficient estimates were somewhat different. In addition, the SE estimates of the fixed effects coefficients were larger for the random effects model, indicating that uncertainty might be larger when using this model. Considering the simulation results, even when there is substantial heterogeneity in β_{i}s, the twostep random effects approach provides valid inference results. The interaction tests of β_{i13} were not significant by significance level of 5% for any of the three models, consistent with the original analyses of Furukawa et al. [14], but these tests should be the most powerful under the assumption that the statistical models were correct.
Discussion
IPD metaanalysis is becoming increasingly popular in recent medical studies, especially as an effective approach for exploring effect modifiers for precision medicine [8]. Given recent trends in clinical trial data sharing [27, 28], the practical values of these methodologies are bound to increase in future studies. The efficiencies of the onestep and twostep approaches have also been discussed for IPD metaanalyses based on the conventional ANCOVA, logistic and Cox regression models [32]. Besides, as noted above, IPD metaanalyses using the mixed effects models involve more complicated problems, which should address the heterogeneity of variancecovariance structures of random effects distributions. In this paper, we presented a flexible twostep method that can effectively circumvent these difficulties, and demonstrated its validity and efficiency.
In the simulation studies, the coverage rates of 95% CIs for the random effects model retained nearly at the nominal level, but previous studies showed that SE are possibly underestimated under small N settings in random effects metaanalyses, and the standard CIs can have serious undercoverage properties [33,34,35,36]. In the simulation studies, the proper coverage properties are possibly caused by the adoption of a compound symmetry structure for the heterogeneity covariance matrix Ω. In multivariate metaanalyses, when the heterogeneity variances are assumed to be common across different outcomes (this assumption is often adopted in network metaanalysis [37, 38]), the statistical information across the different outcomes can be shared for the common parameter estimation, so the accuracy is generally gained. When we adopted different variance assumptions for Ω, the CIs had undercoverage property under small N setting (data not shown). Note that several improved methods for addressing the undercoverage property have been developed, e.g., the KenwardRoger method [39] and the Bartletttype corrections [40].
In this paper, we discussed the linear mixed effects models for continuous Gaussian outcome variables, but the twostep approach can be similarly adapted to the generalized linear mixed effects models for nonGaussian outcome variables [41] and the frailty models for survival outcome variables [42]. For these outcome variables, it should be noted that there are possibly correlations between the estimators of fixed effects coefficients and variancecovariance parameters. Thus, the twostep analyses may cause losses of efficiency on the inferences. However, the twostep approach does not require parametric assumptions for heterogeneities for the withinstudies variancecovariance parameters, and validity of the inference is generally assured. In addition, this method can be easily implemented using standard statistical packages. Therefore, the twostep marginal approach would also be an effective solution for these settings. Further research is warranted to evaluate the expansion of the twostep model to these settings.
Conclusions
The twostep approach is an effective method for synthesizing the regression coefficient parameters efficiently for IPD metaanalyses of longitudinal clinical trials. Using this method, modelling of heterogeneities of the withinstudies variancecovariance parameters can be effectively circumvented without additional parametric assumptions, which are not of primary interests in the analyses of these clinical trials [16, 17, 19]. If we adopt the onestep approach to modelling the heterogeneity, complicated modelling is required and the validity is possibly violated if the adopted model is misspecified. Based on the relaxed assumptions of the twostep approach, we can implement efficient inferences for the treatment effects and interactions. In addition, the twostep approach has advantages for implementations because it can be easily performed using standard statistical packages [23,24,25]. It can also be effectively applied to IPD metaanalyses using clinical trials data sharing systems in which individual trials can only be analyzed separately.
Abbreviations
 ANCOVA:

analysis of covariance
 BIC:

Bayesian information criterion
 CI:

confidence interval
 HRSD:

Hamilton Rating Scale for Depression
 IPD:

individual participant data
 MADRS:

MontgomeryAsberg Depression Rating Scale
 MAR:

missing at random
 ML:

maximum likelihood
 MMRM:

mixed effects model for repeated measures
 REML:

restricted maximum likelihood
 SE:

standard error
References
 1.
Little RJ, D'Agostino R, Cohen ML, Dickersin K, Emerson SS, Farrar JT, Frangakis C, Hogan JW, Molenberghs G, Murphy SA, et al. The prevention and treatment of missing data in clinical trials. N Engl J Med. 2012;367(3):1355–60.
 2.
O'Neill RT, Temple R. The prevention and treatment of missing data in clinical trials: an FDA perspective on the importance of dealing with it. Clin Pharm Ther. 2012;91(3):550–4.
 3.
National Research Council. The prevention and treatment of missing data in clinical trials. Washington, DC: National Academies Press; 2010.
 4.
European Medicines Agency. Guideline on missing data in confirmatory clinical trials. London: European Medicine Agency; 2010.
 5.
Laird NM, Ware JH. Randomeffects models for longitudinal data. Biometrics. 1982;38(4):963–74.
 6.
Verbeke G, Molenberghs G. Linear mixed models for longitudinal data. New York: Springer; 2000.
 7.
Mallinckrodt CH, Clark WS, David SR. Accounting for dropout bias using mixedeffects models. J Biopharm Stat. 2001;11(1–2):9–21.
 8.
Riley RD, Lambert PC, AboZaid G. Metaanalysis of individual participant data: rationale, conduct, and reporting. BMJ. 2010;340:c221.
 9.
Higgins JPT, Green S. Cochrane handbook for systematic reviews of interventions. Chichester: WileyBlackwell; 2008.
 10.
Riley RD, Kauser I, Bland M, Thijs L, Staessen JA, Wang J, Gueyffier F, Deeks JJ. Metaanalysis of randomised trials with a continuous outcome according to baseline imbalance and availability of individual participant data. Stat Med. 2013;32(16):2747–66.
 11.
Riley RD, Price MJ, Jackson D, Wardle M, Gueyffier F, Wang J, Staessen JA, White IR. Multivariate metaanalysis using individual participant data. Res Synth Methods. 2015;6(2):157–74.
 12.
Ishak KJ, Platt RW, Joseph L, Hanley JA, Caro JJ. Metaanalysis of longitudinal studies. Clin Trials. 2007;4(5):525–39.
 13.
Jones AP, Riley RD, Williamson PR, Whitehead A. Metaanalysis of individual patient data versus aggregate data from longitudinal clinical trials. Clin Trials. 2009;6(1):16–27.
 14.
Furukawa TA, Maruo K, Noma H, Tanaka S, Imai H, Shinohara K, Ikeda K, Yamawaki S, Levine SZ, Goldberg Y, et al. Initial severity of major depression and efficacy of antidepressants: individualparticipant data metaanalysis. Acta Psychiatr Scand. 2018;137(6):450–8.
 15.
McCulloch CE, Searle SR, Generalized NJM. Linear, and mixed models. 2nd ed. New York: Wiley; 2008.
 16.
Mallinckrodt CH, Kaiser CJ, Watkin JG, Molenberghs G, Carroll RJ. The effect of correlation structure on treatment contrasts estimated from incomplete clinical trial data with likelihoodbased repeated measures compared with last observation carried forward ANOVA. Clin Trials. 2004;1(6):477–89.
 17.
Mallinckrodt CH, Lane PW, Schnell D, Peng Y, Mancuso JP. Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Info J. 2008;42(4):303–19.
 18.
Higgins JPT, Thompson SG, Spiegelhalter DJ. A reevaluation of randomeffects metaanalysis. J R Stat Soc A. 2009;172(1):137–59.
 19.
Gosho M, Hirakawa A, Noma H, Maruo K, Sato Y. Comparison of biascorrected covariance estimators for MMRM analysis in longitudinal data with dropouts. Stat Methods Med Res. 2017;26(5):2389–406.
 20.
Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixedeffects models using lme4. J Stat Soft. 2015;67(1).
 21.
Jackson D, Riley R, White IR. Multivariate metaanalysis: potential and promise. Stat Med. 2011;30(20):2481–98.
 22.
Mavridis D, Salanti G. A practical introduction to multivariate metaanalysis. Stat Methods Med Res. 2013;22(2):133–58.
 23.
White IR. Multivariate Radomeffects metaanalysis. Stata J. 2009;9:40–56.
 24.
White IR. Multivariate randomeffects metaregression: updates to mvmeta. Stata J. 2011;11:255–70.
 25.
Gasparrini A, Armstrong B, Kenward MG. Multivariate metaanalysis for nonlinear and other multiparameter associations. Stat Med. 2012;31(29):3821–39.
 26.
Cox DR, Reid N. Parameter orthogonality and approximate conditional inference. J R Stat Soc B. 1986;49(1):1–39.
 27.
Taichman DB, Sahni P, Pinborg A, Peiperl L, Laine C, James A, Hong ST, Haileamlak A, Gollogly L, Godlee F, et al. Data sharing statements for clinical trials: a requirement of the International Committee of Medical Journal. Ann Intern Med. 2017;167(1):63–5.
 28.
Taichman DB, Sahni P, Pinborg A, Peiperl L, Laine C, James A, Hong ST, Haileamlak A, Gollogly L, Godlee F, et al. Data sharing statements for clinical trials. BMJ. 2017;357:j2372.
 29.
Lunn D, Spiegelhalter D, Thomas A, Best N. The BUGS project: Evolution, critique and future directions. Stat Med. 2009;28(25):3049–67.
 30.
Carmody TJ, Rush AJ, Bernstein I, Warden D, Brannan S, Burnham D, Woo A, Trivedi MH. The Montgomery Asberg and the Hamilton ratings of depression: a comparison of measures. Eur Neuropsychopharmacol. 2006;16(8):601–11.
 31.
Schwarz G. Estimating the dimension of a model. Ann Stat. 1978;6(2):461–4.
 32.
Burke DL, Ensor J, Riley RD. Metaanalysis using individual participant data: onestage and twostage approaches, and why they may differ. Stat Med. 2017;36(5):855–75.
 33.
Brockwell SE, Gordon IR. A comparison of statistical methods for metaanalysis. Stat Med. 2001;20(6):825–40.
 34.
Jackson D, Riley RD. A refined method for multivariate metaanalysis and metaregression. Stat Med. 2014;33(4):541–54.
 35.
Noma H. Confidence intervals for a randomeffects metaanalysis based on Bartletttype corrections. Stat Med. 2011;30(28):3304–12.
 36.
Bender R, Friede T, Koch A, Kuss O, Schlattmann P, Schwarzer G, Skipka G. Methods for evidence synthesis in the case of very few studies. Res Synth Methods. 2018. https://doi.org/10.1002/jrsm.1297.
 37.
Salanti G, Higgins JP, Ades AE, Ioannidis JP. Evaluation of networks of randomized trials. Stat Methods Me Res. 2008;17(3):279–301.
 38.
White IR. Network metaanalysis. Stata J. 2015;15:951–85.
 39.
Morris TP, Fisher DJ, Kenward MG, Carpenter JR. Metaanalysis of Gaussian individual patient data: twostage or not twostage? Stat Med. 2018;37(9):1419–38.
 40.
Noma H, Nagashima K, Maruo K, Gosho M, Furukawa TA. Bartletttype corrections and bootstrap adjustments of likelihoodbased inference methods for network metaanalysis. Stat Med. 2018;37(7):1178–90.
 41.
Fitzmaurice GM, Laird NM, Ware J. Applied longitudinal analysis. Hoboken, NJ: John Wiley & Sons; 2011.
 42.
Duchateau L, Janssen P. The frailty model. New York: Springer; 2008.
Acknowledgements
Not applicable.
Funding
This study was supported by GrantinAid from the Japan Agency for Medical Research and Development (Grant number: JP17dk0307072) and GrantsinAid for Scientific Research from the Japan Society for the Promotion of Science (Grant numbers: JP15K15954, JP17K19808) in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
The original datasets were available upon agreements between this research group and the sponsor companies, so are not publicly available.
Author information
Affiliations
Contributions
HN, KM, MG contributed to the developments of methods. HN, KM, TAF contributed to the acquisition of data. HN, KM, MG, SZL, YG, SL, TAF contributed to the analysis and interpretation of data. HN drafted the manuscript. All authors reviewed, read and approved the final version to be submitted.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
HN has received lecture fees from Boehringer Ingelheim and Kyowa Hakko Kirin, and research support from Kyowa Hakko Kirin and GSK. He is an associate editor of this journal; MG received lecture and/or consultant fees from Daiichi Sankyo, Ferring Pharmaceuticals, Novartis, and Tiho, received travel fees from Takeda, and received manuscript fees from Kowa; SZL received honoraria or research support or consultancy fees or travel support from Shire Pharmaceuticals, F. HoffmannLa Roche and Eli Lilly; SL has received honoraria for consulting from LB Pharma, Lundbeck, Otsuka, TEVA, Geodon Richter, Recordati, LTS Lohmann, Boehringer Ingelheim, and for lectures from Janssen, Lilly, Lundbeck, Otsuka, SanofiAventis and Servier; TAF has received lecture fees from Janssen, Meiji, MitsubishiTanabe, MSD and Pfizer. He has received research support from MitsubishiTanabe.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Noma, H., Maruo, K., Gosho, M. et al. Efficient twostep multivariate random effects metaanalysis of individual participant data for longitudinal clinical trials using mixed effects models. BMC Med Res Methodol 19, 33 (2019). https://doi.org/10.1186/s1287401906761
Received:
Accepted:
Published:
Keywords
 Individual participant data metaanalysis
 Mixed effects models
 Multivariate random effects model
 Longitudinal studies
 Heterogeneity