Efficient Integration of Aggregate Data and Individual Participant Data in
One-Way Mixed Models

Neha Agarwala 1 and Junyong Park 2 and Anindya Roy1

1Department of Mathematics and Statistics,
University of Maryland, Baltimore County, Baltimore, MD, USA

2Department of Statistics,
Seoul National University, Seoul, South Korea

Abstract

Often both Aggregate Data (AD) studies and Individual Participant Data (IPD)
studies are available for specific treatments. Combining these two sources of
data could improve the overall meta-analytic estimates of treatment effects.
Moreover, often for some studies with AD, the associated IPD maybe available,
albeit at some extra effort or cost to the analyst. We propose a method for
combining treatment effects across trials when the response is from the
exponential family of distribution and hence a generalized linear model
structure can be used. We consider the case when treatment effects are fixed
and common across studies. Using the proposed combination method, we evaluate
the wisdom of choosing AD when IPD is available by studying the relative
efficiency of analyzing all IPD studies versus combining various percentages
of AD and IPD studies. For many different models, design constraints under
which the AD estimators are the IPD estimators, and hence fully efficient, are
known. For such models we advocate a selection procedure that chooses AD
studies over IPD studies in a manner that force least departure from design
constraints and hence ensures an efficient combined AD and IPD estimator.

Key words: meta-analysis, treatment-control difference, random effect,
efficiency, Individual participant data, design.

##  1 Introduction

Meta-analysis of individual participant data (IPD) is the gold standard
statistical approach in systemic reviews of randomized clinical trials. 1
However, if IPD studies systematically differ from studies without access to
IPD, synthesizing information solely based on IPD studies may lead to data
availability bias or reviewer selection bias. 2 In many cases, IPD may not
always be publicly available and access to IPD may be restricted due to
limited resources. If only summary data are available on quantities of
interest, meta-analysis of aggregate data (AD) approaches are the only
recourse an investigator has for combining information about variables of
interest across different studies. In addition to AD studies, if some IPD
studies are available, combining these two levels of data could improve the
overall meta-analysis estimates, compared to utilizing AD studies alone. 3, 4
Substitution of AD with IPD is generally advocated ‘despite the extra cost,
time and complexity required to obtain and manage raw data’. 5 However, in
certain situations the meta-analysis based exclusively on aggregate data (AD-
MA) may yield estimators that are comparable to those based solely on
individual participant data (IPD-MA) and hence substitution of AD with IPD may
not be worthwhile due to the extra cost.

In this paper, we discuss how to aggregate information from IPD studies and AD
studies for a generic model and obtain a combined estimate for the parameter
of interest that is cost effective yet efficient. We denote the combined
estimator based on individual participant data and aggregate data meta-
analysis as IPD-AD-MA. We explore many such combinations and investigate
trade-offs between efficiency gains from substituting AD with IPD versus the
cost of obtaining IPD studies when it is not easily available. The focus of
the paper is to analyze whether it matters which combination of IPD and AD are
selected in obtaining the IPD-AD-MA estimator. We then propose a novel
selection algorithm for selecting which AD to substitute with IPD for maximum
efficiency gain.

One common application of combining studies is testing the effectiveness of
new treatments or medical interventions in randomized clinical trials (RCTs).
Typically, RCTs tend to differ in study design and are often conducted across
different centers, which may result in conflicting evidence. Thus, meta-
analysis approaches have become the norm to integrate the results by borrowing
power across different trials and provide an all-inclusive conclusion. This
paper aims to exploit the design of the trials. Combination of results from
IPD and AD studies work well when the model across the studies is simple with
a common parameter of interest, e.g. treatment effect, and without disparate
covariates. Due to randomization, RCTs on a common treatment provide ideal
setup for exploring combination designs of IPD-AD-MA.

Consider $k$ independent studies with the $j$th study resulting in the
estimated effect size $\hat{\beta}_{j}$, an estimate of the population effect
size $\beta_{j}$ for $j=1,\dots,k$. Suppose $\widehat{v(\hat{\beta}_{j})}$ is
the estimated variance of $\hat{\beta}_{j}$. When studies are homogeneous i.e,
$\beta_{j}=\beta$ for all $j$, $\beta$ is estimated efficiently by a fixed
effects model, using weighted combination of $\hat{\beta}_{j}$ with the
weights taken to be the reciprocal of $\widehat{v(\hat{\beta}_{j})}$. When
studies are heterogeneous, a random effects model is considered to account for
the between study variances, $\tau^{2}$. The meta-analysis estimate of $\beta$
is a weighted mean of $\widehat{v(\hat{\beta}_{j})}$ where the weights for
each study is given by the reciprocal of
$\widehat{v(\hat{\beta}_{j})}+\tau^{2}$. Some widely used estimators for
$\tau^{2}$ are based on the Cochran’s homogeneity test statistic.6, 7, 8
Several other estimators for the heterogeneity variance are available.9, 10,
11, 12, 13, 14, 15 Many papers like van Houwelingen et al (2002), Ritz et al
(2008) have investigated the multivariate extension of these meta-analysis
models.16, 17

The meta-analysis literature is rich with efficiency comparison of AD-MA
estimators with the estimators obtained from full data based on different
models. For treatment vs. control comparison with continuous outcome, Olkin
and Sampson (1998) showed AD-MA estimator is equivalent to IPD-MA estimator if
there is no study-by-treatment interactions and variances are constant across
trials.18 Mathew and Nordstrom (1999) further showed that this equivalence
holds even if the error variances are different across trials.19 The
performance of IPD-MA estimator has been found to be similar but not identical
to AD-MA estimator empirically (e.g. Whitehead, 2002, Ch. 5).20 For a more
general linear model with fixed treatment and random trial effect, Mathew and
Nordstorm, 2010 provided conditions under which AD analysis and IPD analysis
coincide.21 For all commonly used parametric and semi-parametric models, Lin
and Zeng (2010) showed that IPD-MA estimator has no gain in efficiency over
AD-MA estimator asymptotically in the context of fixed effects model and also
provided the condition for their equality.22 Liu et al (2015) introduced a
meta-analysis approach for heterogeneous studies by combining the confidence
density functions derived from the summary statistics of individual studies.23
Doneal et al (2015) compared the performance of IPD-MA and AD-MA using
different estimation procedures in generalized linear mixed model for binary
outcomes.24

One potential limitation in standard meta-analysis approach is the requirement
of a common set of parameters across studies. Different studies often tend to
include different sets of covariates. For meta-analysis of IPD, Jackson et al
proposed a method to estimate the fully adjusted effect across studies with
different set of confounders.25 Recently, Kundu et al developed an extension
of meta-analysis method for fixed-effects model to combine information from
studies with disparate covariate information.26

An extensive efficacy analysis of the one-stage and two-stage statistical
methods for combining IPD and AD in meta-analysis for continuous outcome was
explored by Riley et al among others.27, 28, 29, 30 In the simple situation of
a fixed effects model with a single continuous outcome and covariate,
Yamaguchi et al proposed a method to reconstruct the missing IPD for AD trials
by a Bayesian sampling procedure and use the mixture of simulated IPD and
collected IPD for an IPD meta-analysis.31 Over the past decade, meta-analysis
methods for mixture data have also been developed for dichotomous outcomes and
time-to-event data, some based on reconstruction of IPD.32, 33, 29 Other
popular methods for integrating binary data is random-effects mixed treatment
comparison (MTC) models and likelihood based approaches.34, 35, 36, 37, 38

In section 2, we start with the case of a continuous response following a
linear model and provide the combined treatment effects across trials when the
treatment effect is fixed and common across trials while the trial effect is
random. Assuming the observations within and between the studies are
independent, we investigate the loss of efficiency from using combined
estimator with various percentages of AD and IPD studies. When treatments are
fixed and trial effects are random, Mathew and Nordstorm (2010) derived the
necessary and sufficient condition for the IPD-MA estimator to coincide with
AD-MA estimator for a general within trial covariance matrix. 21 The condition
for equality requires that the fraction of observations corresponding to any
given treatment to be same across trials. In practice, it is more likely to
have studies with differential allocation to treatments. For such models, we
study the relative efficiency of an estimator based on combining IPD and AD
studies, denoted by IPD-AD-MA to the IPD-MA estimator under systematic
departures from the same allocation proportion condition. We further propose a
method to select the IPDs among the available studies so as to get the maximum
efficiency in terms of the combined estimator.

In section 3, we propose a method of combining information across IPD and AD
studies for a multidimensional parameter in a generalized linear mixed model
(GLMM) and study the performance of the combined estimator empirically. This
is a more general setup where the covariate or response may be categorical or
continuous and the common parameter of interest can be multidimensional. In
addition, the random effects may not necessarily be independent. As a special
case, we consider a logistic model with a similar setup to the LMM framework.
For this case, we derive a relative efficiency expression and use the
expression to propose efficient selection of IPD when combining IPD and AD
studies for the IPD-AD-MA estimator.

We use a real data example to illustrate efficient selection of IPD when
synthesizing information on the common parameter of interest across IPD and AD
studies for a linear mixed model and a logistic model with random study
effect. For all our analyses, we assume that all studies or trials are
independent, which is a general and common assumption in the meta-analysis
literature.

##  2 Efficient Aggregation in Linear Mixed Model (LMM)

Consider that there is one continuous outcome of interest and assume that
there are two groups, namely treatment(T) group and control(C) group for all
$k$ independent studies. Let $y_{ji}$ be the response of the $i$th participant
in study $j$ and $x_{ji}$ is the corresponding treatment indicator. Suppose
further that, $n_{jT}$ and $n_{jC}$ be the number of persons in the treatment
group and the control group, respectively, for $j$th study with
$n_{jT}+n_{jC}=n_{j}$. Suppose, for $k_{1}$ studies, IPD are available and for
the remaining $k-k_{1}=k_{2}$ studies, we have access to only AD. Let $S_{1}$
and $S_{2}$ denote the set of IPD studies and AD studies, respectively where
$S_{1}\cup S_{2}=S$, $S$ being the set with all studies. The model is

$\displaystyle y_{ji}=\alpha_{j}+\beta x_{ji}+\epsilon_{ji},\quad j=1,...,k,$
(1) $\displaystyle\epsilon_{ji}\sim N(0,\sigma_{j}^{2}),\quad\alpha_{j}\sim
N(\alpha,\sigma_{\alpha}^{2}),$

where $\alpha_{j}$ and $\epsilon_{ji}$ are assumed independent. Our parameter
of interest is the common treatment effect $\beta$.

###  2.1 Aggregation

For model (1), $\alpha$ acts as a nuisance parameter common across studies.
Thus, the AD-MA estimator and the IPD-MA estimator does not necessarily
coincide although the parameter of interest is common across studies (Mathew
and Nordstorm 2010).21 In fact, the IPD-MA estimator, which is the Best Linear
Unbiased Estimator (BLUE), is more efficient than the AD-MA estimator for any
finite and fixed $k$ and $n_{j}$. 21 The two estimators coincide if and only
if the vectors $(n_{jT}/n_{j},n_{j2}/n_{j})$ are all equal for $j=1,...,k$.
Asymptotically, the AD-MA estimator has the same efficiency as the IPD-MA
estimator.

For the $k_{1}$ studies with access to IPD, $\textbf{y}_{j}$ is normal with
mean and covariance given by

$\displaystyle E(\boldsymbol{y_{j}}|X_{j})$
$\displaystyle=\alpha\boldsymbol{1_{n_{j}}}+\beta X_{j}$ $\displaystyle
Cov(\boldsymbol{y_{j}}|X_{j})$
$\displaystyle=H_{j}=\sigma_{\alpha}^{2}\boldsymbol{1_{n_{j}}}\boldsymbol{1_{n_{j}}^{T}}+\sigma_{j}^{2}I_{n_{j}}.$

For the $k_{2}$ AD studies, the maximum likelihood estimates (MLE),
$\hat{\beta_{j}}$, and their estimated variances,
$\widehat{v(\hat{\beta}_{j})}$ are available. The model for AD study is

$\displaystyle\hat{\beta}_{j}\sim N(\beta,\widehat{v(\hat{\beta}_{j})}),\quad
j=1,\dots,k_{2}$ $\displaystyle
v(\hat{\beta}_{j})={\big{(}n_{j}\pi_{j}(1-\pi_{j})\big{)}}^{-1}{\sigma_{j}^{2}}$

where $\pi_{j}=n_{jT}/n_{j}$ is the proportion of treatment in study $j$.

Integration of IPD and AD uses the standard weighted combination approach with
weights being inversely proportional to the variance, where the variance for
AD part is simply the variance of MLE whereas for the IPD, it is the variance-
covariance matrix of the marginal distribution of the data. With the above
model, the combined estimator of $\boldsymbol{\theta}=(\alpha,\beta)^{\prime}$
and the variance of the combined estimator are

$\displaystyle\hat{\boldsymbol{\theta}}_{IPD-AD-MA}$ $\displaystyle=$
$\displaystyle(U^{T}\Sigma^{-1}U)^{-1}U^{T}\Sigma^{-1}Y^{*},$ $\displaystyle
Cov(\hat{\boldsymbol{\theta}}_{IPD-AD-MA})$ $\displaystyle=$
$\displaystyle(U^{T}\Sigma^{-1}U)^{-1},$ (2)

where
$Y^{*}=(\textbf{y}_{1}^{\prime},\dots,\textbf{y}_{k_{1}}^{\prime},\hat{\beta}_{1},\dots,\hat{\beta}_{k_{2}})^{\prime}$
and

$\displaystyle U=\begin{pmatrix}\boldsymbol{1_{n_{1}}}&X_{1}\\\ \dots&\dots\\\
\boldsymbol{1_{n_{k_{1}}}}&X_{k_{1}}\\\ 0&\boldsymbol{1_{k_{2}}}\\\
\end{pmatrix},\;\Sigma=\begin{pmatrix}\Sigma_{1}&0\\\ 0&\Sigma_{2}\\\
\end{pmatrix},\;\Sigma_{1}$ $\displaystyle=\begin{pmatrix}H_{1}\\\ &\ddots\\\
&&H_{k_{1}}\\\ \end{pmatrix},\;\Sigma_{2}=\begin{pmatrix}v(\hat{\beta}_{1})\\\
&\ddots\\\ &&v(\hat{\beta}_{k_{2}})\\\ \end{pmatrix}.$

Since the estimator in (2) is unbiased, for analyzing efficiency of the
estimator, we obtain the expression for its variance. To derive the variance
of the combined estimator for the treatment effect, we assume $\sigma_{j}^{2}$
and $\sigma_{\alpha}^{2}$ are known. However, both the variance components are
estimated from the data in section 2.3.4 to reflect practical considerations.
The variance is given by

$\displaystyle v(\hat{\beta}_{IPD-AD-MA})={\bigg{[}\sum_{j\in
S_{1}}\frac{n_{j}\pi_{j}(1+n_{j}(1-\pi_{j})b_{j})}{\sigma_{j}^{2}a_{j}}+\sum_{j\in
S_{2}}\frac{n_{j}\pi_{j}(1-\pi_{j})}{\sigma_{j}^{2}}-\frac{(\sum_{j\in
S_{1}}\frac{n_{j}\pi_{j}}{\sigma_{j}^{2}a_{j}})^{2}}{\sum_{j\in
S_{1}}\frac{n_{j}}{\sigma_{j}^{2}a_{j}}}\bigg{]}}^{-1}$ (3)

where
$a_{j}=1+n_{j}b_{j},\;b_{j}=\frac{\sigma^{2}_{\alpha}}{\sigma_{j}^{2}}\;\&\;\pi_{j}=\frac{n_{jT}}{n_{j}}$.

The gold standard for the combined analysis is having access to all IPD
studies. Thus, we compare the efficiency of $\hat{\beta}_{IPD-AD-MA}$ with the
“best” estimator, which is the maximum likelihood estimate (MLE) based on IPD
from all studies. The variance of the “best” estimator of $\beta$, also known
as the minimum variance unbiased estimator is

$\displaystyle v(\hat{\beta}_{IPD-
MA})={\bigg{[}\sum_{j=1}^{k}\frac{n_{j}\pi_{j}(1+n_{j}(1-\pi_{j})b_{j})}{\sigma_{j}^{2}a_{j}}-\frac{\big{(}\sum_{j=1}^{k}\frac{n_{j}\pi_{j}}{\sigma_{j}^{2}a_{j}}\big{)}^{2}}{\sum_{j=1}^{k}\frac{n_{j}}{\sigma_{j}^{2}a_{j}}}\bigg{]}}^{-1}.$

The variance of the combined estimator of $\beta$, assuming only AD results
are available for all $k$ studies, is

$\displaystyle v(\hat{\beta}_{AD-
MA})={\bigg{[}\sum_{j=1}^{k}\frac{n_{j}\pi_{j}(1-\pi_{j})}{\sigma_{j}^{2}}\bigg{]}}^{-1}.$

For the balanced homoscedastic case, when $n_{j}=n$ and
$\sigma_{j}^{2}=\sigma^{2}$, the variance expression for the combined
estimator simplifies to

$\displaystyle v(\hat{\beta}_{IPD-AD-
MA})=\frac{\sigma^{2}}{n}{\bigg{[}\frac{nb}{a}\sum_{j\in
S_{1}}\pi_{j}(1-\pi_{j})+\sum_{j\in
S_{2}}\pi_{i}(1-\pi_{i})+\frac{1}{ak_{1}}\Big{(}\sum_{j\in
S_{1}}\pi_{j}\Big{)}\Big{(}\sum_{j\in S_{1}}(1-\pi_{j})\Big{)}\bigg{]}}^{-1}$

where $a=1+nb,\;b=\frac{\sigma^{2}_{\alpha}}{\sigma^{2}}$. Similarly, the
variance expression for the all IPD estimator reduces to

$\displaystyle v(\hat{\beta}_{IPD-MA})$
$\displaystyle=\sigma^{2}{\bigg{[}\frac{n^{2}b}{a}\sum_{j=1}^{k}\pi_{j}(1-\pi_{j})+\frac{n}{ak}\Big{(}\sum_{j=1}^{k}\pi_{j}\Big{)}\Big{(}\sum_{j=1}^{k}(1-\pi_{j})\Big{)}\bigg{]}}^{-1}.$

The relative efficiency of the combined estimator with respect to the
estimator with all IPD studies is $RE(\hat{\beta}_{IPD-AD-
MA})=\frac{v(\hat{\beta}_{IPD-MA})}{v(\hat{\beta}_{IPD-AD-MA})}$ which for the
simple case when $n_{j}=n$ and $\sigma_{j}^{2}=\sigma^{2}$ is given by

$\frac{\bigg{[}\frac{n^{2}b}{a}\sum_{j\in S_{1}}\pi_{j}(1-\pi_{j})+n\sum_{j\in
S_{2}}\pi_{j}(1-\pi_{j})+\frac{n}{ak_{1}}\Big{(}\sum_{j\in
S_{1}}\pi_{j}\Big{)}\Big{(}\sum_{j\in
S_{1}}(1-\pi_{j})\Big{)}\bigg{]}}{{\bigg{[}\frac{n^{2}b}{a}\sum_{j=1}^{k}\pi_{j}(1-\pi_{j})+\frac{n}{ak}\Big{(}\sum_{j=1}^{k}\pi_{j}\Big{)}\Big{(}\sum_{j=1}^{k}(1-\pi_{j})\Big{)}\bigg{]}}}.$
(4)

We will use the above expression for relative efficiency for simulation
purposes in section 2.3.

###  2.2 Selection

The main question that we address in this paper is whether, while substituting
IPD studies with the corresponding AD results or vise versa, it matters which
studies are selected. If yes, how can we minimize the loss in efficiency
occurring from using AD results instead of the IPD results for a study?

The variance of the combined estimator in equation (3) can be simplified to

$\displaystyle v(\hat{\beta}_{IPD-AD-MA})={\bigg{[}\sum_{j\in
S_{1}}\frac{n_{j}}{\sigma_{j}^{2}a_{j}}(\pi_{j}-\tilde{\pi}_{S_{1}})^{2}+\sum_{j=1}^{k}\frac{n_{j}\pi_{j}(1-\pi_{j})}{\sigma_{j}^{2}}\bigg{]}}^{-1}$
(5)

where $\tilde{\pi}_{S_{1}}=\frac{\sum_{j\in
S_{1}}{\frac{\pi_{j}n_{j}}{\sigma_{j}^{2}a_{j}}}}{\sum_{j\in
S_{1}}{\frac{n_{j}}{\sigma_{j}^{2}a_{j}}}}$.

The expression for variance in equation (5) does not involve the response
$\boldsymbol{y}$ and depends only on $\pi_{j}$, $n_{j}$ and $\sigma_{j}$. This
gives us a way of selecting $k_{1}$ IPD studies among the $k$ studies to
optimize the efficiency of combined estimator. In other words, selection of
IPD studies depends on choosing the best subset $S_{1}$ that maximizes
$\sum_{j\in
S_{1}}\frac{n_{j}}{\sigma_{j}^{2}a_{j}}(\pi_{j}-\tilde{\pi}_{S_{1}})^{2}$.
Essentially this is a combinatorial optimization problem which can be
formulated as follows:

Selection Problem: Given $k>1$ and $1\leq k_{1}\leq k$, and $\mathcal{P}_{k}$,
the power set of $\\{1,\dots,k\\}$, find

$A_{opt}=\underset{A\in\mathcal{P}_{k}:|A|=k_{1}}{\makebox{\rm arg
max}}\,\sum_{j\in A}v_{j}(u_{j}-\bar{u}_{A})^{2}$ (6)

where $(u_{1},v_{1}),\dots,(u_{k},v_{k})$ with $u_{j}>0$, $v_{j}>0$, are
prespecified constants with $\bar{u}_{A}=\frac{\sum_{j\in
A}v_{j}u_{j}}{\sum_{j\in A}u_{j}}$ for any $A\in\mathcal{P}_{k}$.

For the balanced homoscedastic case, the variance further reduces to

$\displaystyle v(\hat{\beta}_{IPD-AD-
MA})={\frac{\sigma^{2}}{n}\bigg{[}\sum_{j\in
S_{1}}\frac{1}{a}(\pi_{j}-\tilde{\pi}_{S_{1}})^{2}+\sum_{j\in
S}\pi_{j}(1-\pi_{j})\bigg{]}}^{-1}.$

For this special case, an exact algorithm exists for finding the optimum set,
$A_{opt}$ of IPD studies. Firstly we arrange the studies in increasing order
of their $\pi_{j}$ values. We include studies with extreme $\pi_{j}$ values in
the IPD set, alternating between the two ends to choose a total of $k_{1}$ IPD
studies.

This is consistent with the result from Mathew and Nordstorm (2010) where the
IPD-MA estimator coincides with AD-MA estimator when the fraction of
observations corresponding to any given treatment is same across trials in a
linear model with fixed treatments and random trial effects. When combining
IPD and AD studies, we propose allocating studies with similar proportion of
treatments to the AD set and studies with widely varying proportion of
treatments to the IPD set.

Algorithm 1 Sequential Selection Algorithm: SSA

$k_{1}\geq 2$

$\mathcal{S}\leftarrow 1:k$

$(m,n)\leftarrow\underset{p,q\in\mathcal{S}}{\arg\max}\;(u[p]-u[q])^{2}/[1/v[p]+1/v[q])$

$A_{SSA}\leftarrow\\{m,n\\}$

$B\leftarrow\mathcal{S}\backslash A_{SSA}$

$D\leftarrow v[m]+v[n]$

$M\leftarrow(v[m]*u[m]+v[n]*u[n])/D$

while $|A_{SSA}|\neq k_{1}$ do

$l\leftarrow\underset{r\in B}{\arg\max}\;v[r]*(u[r]-M)^{2}/(D+v[r])$

$A_{SSA}\leftarrow A_{SSA}\bigcup\\{l\\}$

$B\leftarrow B\backslash\\{l\\}$

$M\leftarrow(D*M+v[l]*u[l])/(D+v[l])$

$D\leftarrow D+v[l]$

end while

For a more general case when $n_{j}$ or $\sigma_{j}$ may not be equal for all
$j$, we propose an approximate sequential algorithm for the selection problem
stated in (6). The sequential selection algorithm, SSA is a forward selection
algorithm that involves choosing the first two studies so that the objective
function is maximized. Then, we select the next $k_{1}-2$ studies sequentially
using the weighted mean based on the preceding selected studies. It sidesteps
the computationally intensive task of calculating the weighted mean for every
subset $S_{1}$ with cardinality $k_{1}$. Two vectors, $u$ and $v$ are fed as
an input to the SSA algorithm and the output is the set $A_{SSA}$ of $k_{1}$
studies. The performance of the algorithm is assessed using simulation studies
and the results are presented in section 2.3.3.

###  2.3 Simulation

Assuming the variance components are known without uncertainty, the variance
expression in (5) does not involve any data and hence the results in the
following sections do not involve any simulation. The discussion on estimation
of the variance components and its sensitivity to the selection algorithm is
covered in section 2.3.4.

For linear mixed effects model, we present two scenarios to show the
importance of selecting appropriate studies for constructing the combined
estimator. For simplicity, we consider both $k$ and $n$ to be small and equal
to $10$ and $\sigma^{2}=2.5,\sigma_{\alpha}^{2}=0.025$.

####  2.3.1 Uniform distribution of proportion of treatment

For this scenario, we use the uniform distribution to randomly generate $k=10$
proportion of treatment over the interval $[0,1]$ as provided in Table 1:

Table 1: Uniform distribution of proportion of treatment across studies study | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
proportion | 0.1 | 0.2 | 0.3 | 0.3 | 0.3 | 0.5 | 0.6 | 0.6 | 0.8 | 0.8

Figure 1 shows that if the desired relative efficiency is 0.9, we could
achieve that with many possible combinations of 60% IPD and 40% AD studies.
However 60% IPD may not be cost effective. But if we choose the right
combination of 40% IPD and 60% AD, we can have a RE of 0.956 in this case. One
such combination which maximizes the RE for 40% IPD is study no 1, 2, 9 and 10
for IPD with treatment proportion 0.1, 0.2, 0.8 and 0.8, respectively, and AD
for the rest of the studies.

Figure 2 shows the maximum and minimum relative efficiency that can be
achieved for different combinations of the IPD and AD studies. In the figure,
the numbers at the top and bottom of the line plot indicate the IPD studies
which achieve the maximum and minimum relative efficiency, respectively. For
example, among all combinations with 20% IPD and 80% AD, the combination which
yields maximum relative efficiency includes study 1 and 10 for the IPD part
and rest of the studies for the AD part. The minimum efficiency is achieved
with the combination of study 3 and 4 for the IPD part and the rest for the AD
part.

####  2.3.2 Bathtub distribution of proportion of treatment

An example of unbalanced proportion of treatment is given in Table 2. The
importance of selecting the best subset of IPD studies for the combined
estimate is more clearly seen in this case. The AD-MA estimator has a relative
efficiency of $0.44$ compared to the IPD-MA estimator. From Figure 3, we see
that it is hard to obtain a relative efficiency of 0.9 even with 80% IPD and
20% AD.

Figure 1: The relative efficiency of all 45 possible combinations for each of 20% IPD and 80% IPD, 210 possible combinations for each of 40% IPD and 60% IPD are plotted. Red dashed line represents the RE for IPD-MA estimator which is 1, blue dashed line is RE for AD-MA estimator which is 0.79 and the green dashed line represents the desired RE, say 0.9, for example. Figure 2: Plot showing the maximum and minimum relative efficiency for each percentage of IPD studies and the numbers at the top and bottom of the line plot indicate the IPD studies which achieves the maximum and minimum relative efficiency, respectively. Table 2: Unbalanced distribution of proportion of treatment across studies study | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
proportion | 0.1 | 0.1 | 0.1 | 0.1 | 0.2 | 0.8 | 0.9 | 0.9 | 0.9 | 0.9
Figure 3: The relative efficiency of all 45 possible combinations for each of
20% IPD and 80% IPD, 210 possible combinations for each of 40% IPD and 60% IPD
are plotted. Red dashed line represents the RE for IPD-MA estimator which is
1, blue dashed line is RE for AD-MA estimator which is 0.79 and the green
dashed line represents the desired RE, say 0.9, for example. Figure 4: Plot
showing the maximum and minimum relative efficiency for each percentage of IPD
studies and the numbers at the top and bottom of the line plot indicate the
IPD studies which achieves the maximum and minimum relative efficiency,
respectively.

The situation can worsen for severely unbalanced distribution of proportion of
treatment. However, the loss in efficiency is an issue only with fixed small
sample size studies. Asymptotically, the AD-MA combined estimator are quite
efficient compared to IPD-MA combined estimator. Hence, the selection of IPD
in data integration is useful for when combining treatment effect across
randomized clinical trials (RCTs) with small cohort size. In such cases, one
can obtain maximum efficiency of the combined estimator for a given number of
IPD studies using the selection algorithm.

####  2.3.3 Selection of IPD studies through sequential algorithm

To assess the performance of sequential algorithm when $n_{j}$ or $\sigma_{j}$
are not equal, we consider $k=30$ studies each with $n=10$ and
$\sigma_{\alpha}^{2}=0.025$. We simulate the $\sigma_{j}$’s and $\pi_{j}$’s as
following:

$\displaystyle\pi$ $\displaystyle\sim 0.5\ \text{Beta}(\alpha=2,\beta=9)+0.5\
\text{Beta}(\alpha=9,\beta=2),$ $\displaystyle\sigma^{2}$
$\displaystyle\sim\text{Inv-Gamma}(\alpha=2,\beta=5).$

For $k_{1}=2,\dots,10$, we first obtain the optimal set with $k_{1}$ IPD
studies for which the variance is minimized and then record the number of
times the optimal set matches with the set obtained through the sequential
algorithm. Table 3 presents the count of the matches for 100 simulations and
the mean ratio of variance for the set from sequential algorithm to the
variance of the optimal set. For larger $k_{1}$, the number of combinations is
large and hence finding the optimal set can be time consuming and
computationally expensive. The sequential algorithm, on the other hand, is
quite fast and performs relatively well.

Table 3: Count of matches for a set of size $k_{1}$ out of $k=10$ studies | $k_{1}$
---|---
Count Match | $10$ | $9$ | $8$ | $7$ | $6$ | $5$ | $4$ | $3$ | $2$
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0
2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 8 | 100
3 | 0 | 0 | 0 | 0 | 0 | 3 | 15 | 90 |
4 | 0 | 0 | 0 | 0 | 0 | 14 | 83 |  |
5 | 0 | 0 | 0 | 0 | 7 | 83 |  |  |
6 | 0 | 0 | 0 | 3 | 93 |  |  |  |
7 | 0 | 0 | 6 | 97 |  |  |  |  |
8 | 0 | 6 | 94 |  |  |  |  |  |
9 | 9 | 94 |  |  |  |  |  |  |
10 | 91 |  |  |  |  |  |  |  |
Mean Ratio | 1 | 1 | 1 | 1 | 1 | 1.001 | 1.001 | 1.001 | 1

####  2.3.4 Variance component estimation

In this section, we discuss about the estimation of the variance components
$\sigma^{2}_{j}$ and $\sigma^{2}_{\alpha}$ and its sensitivity to the
sequential selection algorithm. When $\sigma^{2}_{j}\neq\sigma^{2}_{l}$ for
some $j,l$, $\sigma^{2}_{j}$ is estimated using the estimated variance of the
treatment effect and the number of participants for the two groups in each
study. One can also test for homogeneity of study variances and obtain a
pooled estimate of $\sigma^{2}$ with the assumption
$\sigma^{2}_{j}=\sigma^{2}$ $\forall j$. Furthermore, $\sigma_{\alpha}^{2}$ is
estimated using some pilot IPD studies, that are accessible without any
difficulty or effort.

To study the sensitivity of the variance component estimation to the selection
algorithm, two different scenarios are considered: $\sigma^{2}_{j}=\sigma^{2}$
$\forall j$ and $\sigma^{2}_{j}\neq\sigma^{2}_{l}$ for some $j,l$. For each
scenario, we generate $10000$ data sets from the linear mixed effects model in
(1) where $\alpha=0.5,\beta=1.5,\sigma^{2}_{\alpha}=0.025$ with $k=10$ studies
with $n_{j}=50$ participants and proportion of treatment as provided in table
2. The $\sigma_{j}$’s are simulated as:

$\sigma^{2}\sim\text{Inv-Gamma}(\alpha=2,\beta=5).$

For each simulated data set, we find the optimal set of $k_{1}=5$ IPD studies
when the $\sigma_{j}$’s and $\sigma_{\alpha}$’s are known without uncertainty.
Next, we derive the estimates of the variance components and obtain the set of
$k_{1}=5$ IPD studies using the sequential algorithm based on these estimates.
Table 4 reports the count of exact matches between the two sets. For both the
cases, we found that the selection algorithm is robust to the estimation of
these parameters.

Table 4: Count of matches for a set of size $k_{1}=5$ out of $k=10$ studies Count Match | $\sigma_{j}^{2}=\sigma^{2}$ | $\sigma_{j}^{2}\neq\sigma^{2}$
---|---|---
0 | 0 | 0
1 | 0 | 0
2 | 0 | 0
3 | 0 | 386
4 | 607 | 3636
5 | 9393 | 5974

###  2.4 World Values Survey: Efficient Recombination

We illustrate the significance of selection of IPD studies when combining
information across studies using data from the World Values Survey (WVS)
(https://www.worldvaluessurvey.org/wvs.jsp). The WVS is a large consortium of
social survey data from around 100 regions across 6 waves (1981-1984,
1990-1994, 1995-1998, 1999-2004, 2005-2009 and 2010-2017). The data consist of
scores based on questionnaire covering a broad range of topics such as
economic development, democratization, religion, gender equality, social
capital, and subjective well-being.

For our analysis, we consider a simple linear mixed effects model based on the
data for wave 6 in USA which contains 2232 participants. The dependent
variable is life satisfaction score which ranges from 1 (completely
dissatisfied) to 10 (completely satisfied). Based on financial satisfaction
score, the predictor is a binary variable taking the value 1 for scores more
than 5 and 0 otherwise. We consider category 1 as treatment and 0 as control
for the predictor. We focus on individuals who are single and also account for
heterogeneity due to age, sex, ethnicity and education in the study population
through stratification. The variable age is categorized into intervals of 20
years as 0-20, 20-40, 40-60, 60-80 and more than 80. Ethnicity has 5
categories: Non-Hispanic white, Non-Hispanic Black, Hispanic, Non-Hispanic
more than 2 races and other races while sex is categorized into male and
female groups. The variable education has many categories starting with no
formal education to doctorate degree. The strata are formed by considering
different combinations of age, sex, ethnicity and education. We remove strata
that have either only 1 or only 0 as response, resulting in treatment
proportions to be 1 or 0. We also remove strata with number of participants
(stratum size) less than 3. The total number of strata after exclusion is 30
with strata sizes ranging from 3 to 24 and the proportion of treatment ranging
from 0.1 to 0.89. Each strata will be considered as a study and meta analysis
will be in combining the strata specific results to obtain overall population
results.

We assume that we have access to the AD for each study, i.e., the estimate of
the treatment effect and the estimated standard error along with the number of
participants in the treatment group and the control group are available for
each study. The objective is to assess whether getting IPD is useful and
determine which IPDs should be combined with AD to obtain maximum efficiency.
We randomly sample 5 pilot studies among the 30 total studies and estimate
$\sigma_{\alpha}^{2}=0.144$ using model (1). First, we consider the case when
$\sigma^{2}_{j}$’s are not same. The estimated variance for an all-IPD
estimator of the treatment effect, $\hat{\beta}_{IPA-MA}$ is 0.033 whereas for
an estimator based on all-AD estimator, $\hat{\beta}_{AD-MA}$ the estimated
variance is 0.038 resulting in an approximately $11\%$ loss in efficiency.
However, if we include 5 best subset IPD studies selected using the algorithm
we proposed, the estimated variance for the combined estimator,
$\hat{\beta}_{IPD-AD-MA}$ is 0.035 and hence we have a $7\%$ gain in
efficiency compared to $\hat{\beta}_{AD-MA}$.

Under the assumption $\sigma^{2}_{j}=\sigma^{2}$, there is an approximate
$12\%$ loss in efficiency for the AD-MA estimator compared to the IPD-MA
estimator. If we include the 5 best chosen IPD studies, there is a $9\%$ gain
in efficiency for the IPD-AD-MA estimator relative to the AD-MA estimator.
Figure 5 shows which 5 studies are selected using the sequential algorithm.
Both assumptions yield nearly the same set of best 5 IPD studies which
reinforces that the selection algorithm is not sensitive to variance component
estimation. Table 5 summarizes the pooled estimates and the standard errors
for various scenarios.

Table 5: Pooled estimate and standard error of the treatment effect for IPD-MA estimator, IPD-AD-MA estimator with 5 best subset IPD and the AD-MA estimator | $\sigma_{j}^{2}=\sigma^{2}$ |  | $\sigma_{j}^{2}\neq\sigma^{2}$
---|---|---|---
IPD Trials | $\boldsymbol{\hat{\beta}}$ | se($\boldsymbol{\hat{\beta}}$) | IPD Trials | $\boldsymbol{\hat{\beta}}$ | se($\boldsymbol{\hat{\beta}}$)
None | 1.624 | 0.221 | None | 1.603 | 0.195
1:3, 29:30 | 1.762 | 0.207 | 1:2, 28:30 | 1.746 | 0.182
1:30 | 1.784 | 0.205 | 1:30 | 1.774 | 0.179
Figure 5: Plot showing the studies that are selected against the proportion of
treatments when the study variances are homogeneous versus when they are
heterogeneous.

##  3 Aggregation in Generalized Linear Mixed Model

For studies with non-Gaussian data, generalized linear models (GLM) form a
large class, including popular models for binary response and count response.
To aggregate several GLM studies with the same parameters of interest, one
could use a meta analysis approach similar to that in the linear model case.
We consider a generalized linear mixed model (GLMM) setup which allows for
random effect terms for the study effect. We use a simple aggregation
framework which will allow us to derive an efficiency expression for the
treatment effect estimator similar to that in the LMM. We could then use the
efficiency expression, written in terms of the treatment allocations to decide
which AD could be replaced by IPD for maximum gain in efficiency.

We present the setup in a slightly more general form with several fixed
effects and random effects terms. Suppose, for $j=1,\dots,k_{1}$ studies in
set $S_{1}$, we have access to the IPD data whereas for $k_{2}$ studies in set
$S_{2}$, only the summary statistics for the parameter of interest are
available where $S_{1}\cup S_{2}=S$ is the total set of studies with $|S|=k$.

The IPD data $(Y_{ji},X_{ji}|\beta_{j},\alpha_{j})\sim
f_{j}(Y_{ji},X_{ji};\beta_{j},\alpha_{j})$ where $f$ is a probability density
belonging to the regular exponential family. We assume a random effect model
for $\beta_{j}$ and $\alpha_{j}$, where both $\beta_{j}$ and $\alpha_{j}$ may
be multidimensional. We further assume a joint multivariate gaussian model for
the random effects allowing for non-zero correlation between the $\beta_{j}$
and $\alpha_{j}$ where $\beta$ is the parameter vector of interest and the
vector $\alpha$ includes all nuisance parameters:

$\displaystyle\begin{pmatrix}\beta_{j}\\\ \alpha_{j}\\\
\end{pmatrix}\bigg{\rvert}\begin{pmatrix}\beta\\\ \alpha\\\ \end{pmatrix}\sim
MVN\Bigg{(}\begin{pmatrix}\beta\\\ \alpha\\\
\end{pmatrix},\begin{pmatrix}\Sigma_{\beta\beta}&\Sigma_{\beta\alpha}\\\
\Sigma_{\alpha\beta}&\Sigma_{\alpha\alpha}\\\ \end{pmatrix}\Bigg{)}$

with $MVN$ denoting the multivariate normal density. This is a general model
and can be reduced to specific mixed effects models or fixed effects model by
constraining the corresponding components in the covariance matrix to be zero.
This setup allows having additional covariates and the parameter corresponding
to the additional covariates can be included in $\alpha_{j}$.

For the $k_{2}$ AD studies, we have the MLE estimates of the main parameter
vector with their estimated covariance matrices,
$(\hat{\beta}_{j},\widehat{V(\hat{\beta}_{j})})$. The model for random effects
is given by

$\displaystyle\hat{\beta}_{j}|\beta_{j}$ $\displaystyle\sim
N(\beta_{j},\widehat{V(\hat{\beta}_{j})})$ $\displaystyle\beta_{j}|\beta$
$\displaystyle\sim N(\beta,\Sigma_{\beta}).$

We will use a composite likelihood where the likelihood from the studies with
IPD and the likelihood from the studies with AD are simply multiplied
together.39 For the studies with only summary statistics available, we
consider an approximate likelihood based on the asymptotic normal model for
the maximum likelihood estimators (MLE). For the studies with IPD, we can
write the full GLMM likelihood. But given the complexity of the GLMM
likelihood involving integrals we use an approximate likelihood based on the
Laplace approximation. We estimate the common parameter of interest using
profile likelihood approach applied to the composite likelihood.

The full (log) likelihood function is obtained by combining IPD studies and AD
studies as:

$\displaystyle L(\beta,\alpha)$ $\displaystyle=\prod_{j\in
S_{1}}{L_{j}(\beta,\alpha)}\prod_{j\in S_{2}}{L_{j}(\beta)},$
$\displaystyle\log L(\beta,\alpha)$ $\displaystyle=l(\beta,\alpha)=\sum_{j\in
S_{1}}{l_{j}(\beta,\alpha)}+\sum_{j\in S_{2}}{l_{j}(\beta)}.$

For each $k\in S_{1}$, the IPD part of log-likelihood, $l_{j}(\beta,\alpha)$
is constructed using the Laplace approximation (where the family $f$ satisfies
the standard smoothness assumptions for the Laplace approximation to work) to
expand $l_{j}(\beta,\alpha)$ around the MLE,
$(\hat{\beta}_{j},\hat{\alpha}_{j})$ from study $j$ and ignoring the higher
order terms $o(||\beta_{j}-\hat{\beta}_{j}||^{2})$. Then,
$(\beta_{j},\alpha_{j})$ is integrated out in the Gaussian integral to obtain

$l_{j}(\beta,\alpha)=l_{j}(\hat{\beta}_{j},\hat{\alpha}_{j})-\frac{1}{2}\Bigg{[}\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\ \end{pmatrix}-\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\
\end{pmatrix}\Bigg{]}^{T}\Delta_{j}^{-1}\Bigg{[}\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\ \end{pmatrix}-\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\ \end{pmatrix}\Bigg{]}-\frac{1}{2}\log|\Delta_{j}|$ (7)

where
$\Delta_{j}=\Bigg{[}\begin{matrix}\begin{pmatrix}\Sigma_{\beta\beta}&\Sigma_{\beta\alpha}\\\
\Sigma_{\alpha\beta}&\Sigma_{\alpha\alpha}\\\
\end{pmatrix}+\Bigg{[}\begin{pmatrix}I_{\beta_{j}\beta_{j}}&I_{\beta_{j}\alpha_{j}}\\\
I_{\alpha_{j}\beta_{j}}&I_{\alpha_{j}\alpha_{j}}\\\
\end{pmatrix}\bigg{\rvert}_{\hat{\beta}_{j},\hat{\alpha}_{j}}\Bigg{]}^{-1}\end{matrix}\Bigg{]}^{-1}$.

The AD part of likelihood for each $k\in S_{2}$ is given by

$l_{j}(\beta)={-\frac{1}{2}\begin{pmatrix}\hat{\beta_{j}}-\beta\\\
\end{pmatrix}^{T}\big{(}\widehat{V(\hat{\beta}_{j})})+\Sigma_{\beta\beta}\big{)}^{-1}\begin{pmatrix}\hat{\beta_{j}}-\beta\\\
\end{pmatrix}}\\\
-\frac{1}{2}\log|\widehat{V(\hat{\beta}_{j})})+\Sigma_{\beta\beta}|.$ (8)

Putting (7) and (8) together, we have the full log likelihood for
$(\beta,\alpha)$ as

$\displaystyle l(\beta,\alpha)$ $\displaystyle=$ $\displaystyle\sum_{j\in
S_{1}}l_{j}(\hat{\beta}_{j},\hat{\alpha}_{j})-\frac{1}{2}\sum_{j\in
S_{1}}\Bigg{[}\begin{pmatrix}\hat{\beta_{j}}\\\ \hat{\alpha_{j}}\\\
\end{pmatrix}-\begin{pmatrix}\hat{\beta_{j}}\\\ \hat{\alpha_{j}}\\\
\end{pmatrix}\Bigg{]}^{T}\Delta_{j}^{-1}\Bigg{[}\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\ \end{pmatrix}-\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\ \end{pmatrix}\Bigg{]}-\frac{1}{2}\sum_{j\in
S_{1}}\log|\Delta_{j}|$ (9) $\displaystyle-\frac{1}{2}\sum_{j\in
S_{2}}(\hat{\beta_{j}}-\beta)^{T}\big{(}\widehat{V(\hat{\beta}_{j})})+\Sigma_{\beta\beta}\big{)}^{-1}(\hat{\beta_{j}}-\beta)-\frac{1}{2}\sum_{j\in
S_{2}}\log|\widehat{V(\hat{\beta}_{j})}+\Sigma_{\beta\beta}|.$

which then can be maximized to obtain the estimates of mean parameters and the
variance components. To reduce the computational burden, we lower the
dimension by expressing the mean parameter as a function of the variance
components in (10).

$\displaystyle\begin{pmatrix}\hat{\beta}(\Sigma)\\\ \hat{\alpha}(\Sigma)\\\
\end{pmatrix}$ $\displaystyle=\Bigg{[}\sum_{j\in S_{1}}\Delta_{j}+\sum_{j\in
S_{2}}\Omega_{j}\Bigg{]}^{-1}\ \Bigg{[}\sum_{j\in
S_{1}}\Delta_{j}\begin{pmatrix}\hat{\beta_{j}}\\\ \hat{\alpha_{j}}\\\
\end{pmatrix}+\sum_{j\in S_{2}}\Omega_{j}\begin{pmatrix}\hat{\beta_{j}}\\\
\hat{\alpha_{j}}\\\ \end{pmatrix}\Bigg{]}$ (10)

where

$\displaystyle\Delta_{j}$
$\displaystyle=\Delta_{j}(\Sigma)=\Bigg{[}\begin{pmatrix}\Sigma_{\beta\beta}&\Sigma_{\beta\alpha}\\\
\Sigma_{\alpha\beta}&\Sigma_{\alpha\alpha}\\\
\end{pmatrix}+\Bigg{[}\begin{pmatrix}I_{\beta_{j}\beta_{j}}&I_{\beta_{j}\alpha_{j}}\\\
I_{\alpha_{j}\beta_{j}}&I_{\alpha_{j}\alpha_{j}}\\\
\end{pmatrix}\bigg{\rvert}_{\hat{\beta}_{j},\hat{\alpha}_{j}}\Bigg{]}^{-1}\Bigg{]}^{-1}$
$\displaystyle\Omega_{j}$
$\displaystyle=\Omega_{j}(\Sigma_{\beta\beta})=\Bigg{[}\begin{matrix}(\Sigma_{\beta\beta}+\widehat{V(\hat{\beta}_{j})})^{-1}&0\\\
0&0\\\ \end{matrix}\Bigg{]}.$

We then maximize (9) with respect of $\Sigma$ to obtain $\hat{\Sigma}$ and use
the expression from (10) to obtain the combined estimator,
$\hat{\beta}(\hat{\Sigma})$ denoted by $\hat{\beta}_{IPD-AD-MA}$. The
estimated variance of $\beta$ is then given by

$\displaystyle\widehat{V(\hat{\beta}_{IPD-AD-MA})}$
$\displaystyle=\Bigg{[}\sum_{j\in S_{1}}\Delta_{j}(\hat{\Sigma})+\sum_{j\in
S_{2}}\Omega_{j}(\hat{\Sigma}_{\beta\beta})\Bigg{]}^{-1}_{\ \ [1,1]}.$ (11)

In the GLMM framework, for the IPD-MA estimator, we essentially combine the
MLE estimates and the estimated variances for all the parameters in the model,
including the parameter of interest and the nuisance parameter whereas for the
AD-MA estimator, we integrate MLE estimates and the estimated variances for
only the parameter of interest. Thus, there may be a loss of information for
the AD-MA estimator when compared to the IPD-MA estimator depending on the
correlation of the parameter of interest with the nuisance parameters.

In GLMM, we use the data to get the MLE estimates and integrate them to derive
the IPD-MA estimator. In LMM, however, one can directly use the data for the
IPD studies instead of using the MLE estimates. While the general methodology
of aggregation of GLMM studies can be applied to linear mixed models with the
identity link, explicit computation of the likelihood and the estimators in
the LMM allows more efficient strategies for data aggregation.

###  3.1 Selection in the Logistic Model

For a special case of generalized linear mixed model, we provide the variance
expression for IPD selection. We consider a logistic model with similar setup
as (1), where the treatment effect, $\beta$ is fixed and the trial effect is
random with with $\alpha_{j}\sim N(\alpha,\sigma^{2}_{\alpha})$. Suppose
further that, $n_{jT}$ and $n_{jC}$ be the number of persons allocated for
treatment and control, respectively, for $j$th study where
$n_{jT}+n_{jC}=n_{j}$. Let $x_{ji}$ denote the treatment allocation for
individual $i$ in study $j$, with proportion of treatment $\pi_{j}$ for
$j=1,\dots,k$ and $i=1,\dots,n_{j}$. The response $y_{ji}$ is binary with
probability of success,

$\displaystyle P(y_{ji}=1|x_{ji}=1)$ $\displaystyle=$ $\displaystyle p_{j1}$
$\displaystyle=\frac{\exp(\alpha_{j}+\beta)}{1+\exp(\alpha_{j}+\beta)},$ (12)
$\displaystyle P(y_{ji}=1|x_{ji}=0)$ $\displaystyle=$ $\displaystyle p_{j0}$
$\displaystyle=\frac{\exp(\alpha_{j})}{1+\exp(\alpha_{j})}.$

Using the variance expression in (11), the variance of the combined estimator
can be simplified to

$\displaystyle V(\hat{\beta}_{IPD-AD-MA})$ $\displaystyle={\bigg{[}\sum_{j\in
S_{1}}\frac{h_{j}}{c_{j}}(g_{j}-\tilde{g})^{2}+\sum_{j=1}^{k}h_{j}^{-1}\bigg{]}}^{-1}$
(13)

where

$\displaystyle h_{j}$
$\displaystyle=a_{j}^{-1}+b_{j}^{-1}=V(\hat{\beta}_{j}),$ $\displaystyle
g_{j}$ $\displaystyle=\frac{1}{{(a_{j}/b_{j})}^{-1}+1},$ $\displaystyle c_{j}$
$\displaystyle=\sigma^{2}_{\alpha}h_{j}+a_{j}^{-1}b_{j}^{-1},$ $\displaystyle
a_{j}$ $\displaystyle=n_{jT}\ p_{j1}(1-p_{j1}),$ $\displaystyle b_{j}$
$\displaystyle=n_{jC}\ p_{j0}(1-p_{j0}).$

The problem of using (13) to select the “best” combination of IPD and AD is
identical to the selection problem stated in (6). In order to use the
sequential selection algorithm proposed in section 2.2, the unknowns $g_{j}$
and $h_{j}/c_{j}$ need to be estimated using AD. There are two possible
solutions.

First, with the assumption of rare disease,
$a_{j}/b_{j}=({n_{jT}}/{n_{jC}})*(\text{odds ratio}_{j})$ where the
$\text{odds ratio}_{j}$ can be estimated from $\hat{\beta}_{j}$. Therefore,
with this assumption $h_{j}/c_{j}$ and $g_{j}$ can be computed as

$\displaystyle g_{j}$
$\displaystyle=\frac{1}{({n_{jT}}/{n_{jC}})*\exp{(\hat{\beta}_{j})}+1}$ (14)
$\displaystyle\frac{h_{j}}{c_{j}}$
$\displaystyle=\frac{1/\widehat{V(\hat{\beta}_{j})}}{\sigma^{2}_{\alpha}/\widehat{V(\hat{\beta}_{j})})+g_{j}(1-g_{j})}$

provided we have AD for all studies along with $n_{jT}$ and $n_{jC}$. We can
then compute the variance expression (13) to choose the IPD studies that
‘maximize’ the efficiency of the combined estimator using the sequential
algorithm (1) .

The second solution requires knowledge of the number of cases and controls for
each study in addition to treatment and control group totals. In fact, with
this information we can reconstruct the individual $2\times 2$ contingency
tables and hence estimate $g_{j}$ and $h_{j}/c_{j}$ even without the rare
disease assumption. However, efficient selection of IPD studies is of no
additional advantage for binary data in RCTs in case the $2\times 2$
contingency tables are available since one can construct the IPD from the
table.

###  3.2 Simulation results

We illustrate the advantages of combining estimates in GLMM through a limited
simulation experiment using binary data and the logistic link. Suppose
$X_{ji}$ denotes the treatment allocation for individual $i$ in study $j$,
with proportion of treatment, $\pi_{j}$ where $j=1,\dots,k$ and
$i=1,\dots,n_{j}$. The response $y_{ji}$ is binary with probability of success
as in model (12) where both $\beta_{j}$ and $\alpha_{j}$ are assumed to be
random with $\beta_{j}|\beta\sim N(\beta,\Sigma_{\beta})$ and
$\alpha_{j}|\alpha\sim N(\alpha,\Sigma_{\alpha})$, independently, where the
true parameters are taken to be $\beta=0.5,\alpha=0.5,\sigma^{2}_{\beta}=0.5$
and $\sigma^{2}_{\alpha}=0.5$. We simulated the true $\pi_{j}$’s from a
uniform distribution over the interval $[0,1]$. The results are shown in Table
6. We present the estimates, estimated bias and standard error for different
choices of $n=n_{j}$ and $k$ along with the relative efficiency with respect
to the IPD-MA estimator. Similar to the linear model, we could have different
combinations of IPD and AD for the IPD-AD-MA estimator, even when the
percentage of IPD is fixed. However, since there can be many such combinations
for $k=50$, we only report the results for one randomly chosen combination for
each possibility of no IPD, $20\%$ IPD, $40\%$ IPD, $60\%$ IPD, $80\%$ IPD and
all IPD.

The mean square error (MSE) for the all AD estimator is 0.021 for
$k=50,n=100$, 0.013 for $k=50,n=500$, 0.012 for $k=100,n=100$ and 0.006 for
$k=100,n=500$. The MSE for the all IPD estimator is 0.019 when $k=50,n=100$
and 0.010 for $k=100,n=100$ and stays the same for other choices of $k$ and
$n$. While the bias doesn’t change much when $n$ is fixed and $k$ is varying,
the MSE decreases with increase in either sample size or number of studies or
both. On the other hand, the relative efficiency approaches 1 as percentage of
IPD studies is increased for each choice of $k$ and $n$.

Table 6: Model parameter estimates for different choices of $n$ and $k$ Scenario | (% IPD, % AD) | Estimate | Bias | Std. Error | Rel. eff.
---|---|---|---|---|---
k=50, n=100 | (0, 100) | 0.43 | -0.07 | 0.127 | 0.926
| (20, 80) | 0.43 | -0.07 | 0.127 | 0.925
| (40, 60) | 0.431 | -0.069 | 0.126 | 0.94
| (60, 40) | 0.433 | -0.067 | 0.125 | 0.961
| (80, 20) | 0.436 | -0.064 | 0.125 | 0.982
| (100, 0) | 0.44 | -0.06 | 0.126 | 1.000
k=50, n=500 | (0, 100) | 0.482 | -0.018 | 0.113 | 0.975
| (20, 80) | 0.482 | -0.018 | 0.113 | 0.976
| (40, 60) | 0.483 | -0.017 | 0.113 | 0.979
| (60, 40) | 0.483 | -0.017 | 0.113 | 0.983
| (80, 20) | 0.485 | -0.015 | 0.113 | 0.986
| (100, 0) | 0.488 | -0.012 | 0.112 | 1.000
k=100, n=100 | (0, 100) | 0.438 | -0.062 | 0.09 | 0.875
| (20, 80) | 0.438 | -0.062 | 0.09 | 0.873
| (40, 60) | 0.439 | -0.061 | 0.09 | 0.89
| (60, 40) | 0.441 | -0.059 | 0.089 | 0.911
| (80, 20) | 0.444 | -0.056 | 0.089 | 0.942
| (100, 0) | 0.448 | -0.052 | 0.088 | 1.000
k=100, n=500 | (0, 100) | 0.482 | -0.018 | 0.074 | 0.973
| (20, 80) | 0.482 | -0.018 | 0.074 | 0.975
| (40, 60) | 0.483 | -0.017 | 0.074 | 0.979
| (60, 40) | 0.484 | -0.016 | 0.074 | 0.986
| (80, 20) | 0.485 | -0.015 | 0.074 | 0.994
| (100, 0) | 0.488 | -0.012 | 0.074 | 1.000

###  3.3 Real data analysis

We consider the dataset in Yusuf et al. (1985) which includes results from the
long-term trials of oral beta blockers on its effectiveness for reducing
mortality.40 The data consists of $2\times 2$ tables from 22 clinical trials.
We illustrate the aggregation of (log) odds ratio under the logistic model
(12) with independent random effects for the study effect and treatment
effect. Although the full data is available, for the purpose of illustration
we assume access to IPD for some trials and access to only meta-analysis
results for the remaining trials. The cohort size for the smallest trial is 77
whereas the largest trial has a cohort of 3837 individuals. To show the
significance of selection for the logistic model, we randomly sampled
$n_{j}=50$ individuals for each study $j$.

We assume we have access to the AD for each study,i.e. the MLE estimates of
the log-odds ratio and its standard error along with the treatment and control
group totals. We apply the sequential algorithm using the estimates in (14)
where the $\sigma^{2}_{\alpha}$ is estimated using a randomly selected pilot
IPD studies. In practice, it can be estimated using IPD studies that are
conveniently accessible since the $\sigma^{2}_{\alpha}$ seems to have little
or no effect in choosing the final set of selected IPD studies. The estimated
variance for the IPD-MA estimator and AD-MA estimator are 0.033 and 0.044,
respectively, which implies a $26\%$ loss in efficiency. If we include the
best 5 IPD studies chosen using the sequential algorithm, the estimated
variance for the IPD-AD-MA estimator is 0.035, resulting in a mere $6\%$ loss
in efficiency relative to the IPD-MA estimator, which is a gain of $26\%$
efficiency relative to the AD-MA estimator. However, the IPD-AD-MA estimator
with the worst 5 IPD studies results in a estimated variance of 0.044, which
is no any efficiency gain from the AD-MA estimator.

The pooled estimate (log) odds ratio with its standard error for the IPD-MA
estimator, the best combination of IPD-AD studies with $k_{1}=2,5$ and $8$,
respectively and the AD-MA estimator are reported in Table 7. The IPD trials
in the table are ordered with respect to the $g_{j}$, analogous to the
proportion of treatment, $\pi_{j}$ in the linear mixed effects model. The
standard error for the best combined estimator decreases as number of IPD
studies is increased. In conclusion, a combined estimator with just 8 selected
IPD studies out of the total 22 studies can achieve efficiency very close to
the IPD-MA estimator.

Table 7: Pooled estimate and standard error of the (log) odds ratio for the IPD-MA estimator, the best IPD-AD-MA estimator with 2 IPDs, 5 IPDs and 8 IPDs, respectively and the AD-MA estimator IPD Trials | $\boldsymbol{\hat{\beta}}$ | se($\boldsymbol{\hat{\beta}}$)
---|---|---
None | 0.001 | 0.211
1,22 | 0.006 | 0.206
1:3,21:22 | 0.017 | 0.204
1:5,20:22 | 0.001 | 0.200
1:22 | 0.057 | 0.199

##  4 Discussion

In this paper, we provide a method of combining information across independent
studies for a generalized mixed effects model in a multivariate setup. We show
that the combined estimator is efficient compared to the all IPD estimator
asymptotically through simulation. For a much simpler linear mixed effects
model with heterogeneous studies, we investigate the performance of the
combined estimator for various distribution of treatment proportion across
studies. We advocate a method for selection of AD studies over IPD studies to
ensure a fully efficient combined estimator. When the number of studies is
large, we propose an approximate sequential algorithm to select the best
combination of IPD and AD studies.

We assume the same set of covariates for each study. In the future, we plan to
expand our method to include disparate covariate information for multivariate
mixed effects models. Another possible direction would be to quantify the
degree of unbalancedness of proportion of treatments among studies. That would
help in investigating the impact of selecting AD studies for multiple
treatment problems as well. It would be interesting to extend the sequential
algorithm to a more general GLM framework, which requires a result on
optimality of the AD estimators in terms of the design parameters. Once such a
result has been established then an analogous algorithm could be devised by
optimizing a measure of departure from the optimality condition.

We also look at an interesting application of our novel selection algorithm in
the context of the recently popular ‘Split and Conquer’ approach to analysis
of ‘Big Data’. With the increasing need of resources to store and analyze
large data, meta-analysis methods are becoming very popular. Due to the
availability and accessibility of massive amount of data in several fields, a
new research paradigm is focused towards divide and recombine (D&R) approach
for big data analysis.41, 42, 43 The main goal of our method is to aid in the
analysis of Big Data like Genome Wide Association Studies, Electronic Health
Records, large socio-economic survey data, using the split, analyze and
aggregate approach in order to reduce computational and storage limitations
but still provide estimates with efficiency close to that from the IPD
analysis.

##  5 Data Availability Statement

The data that support the findings of this study are available in github
repository ’meta-analysis’ at https://github.com/agneha2010/meta-analysis.
These data were derived from the following resources available in the public
domain: Yusuf et al. (1985) and https://www.worldvaluessurvey.org/wvs.jsp.40

## References

  * 1 Chalmers I. The Cochrane collaboration: preparing, maintaining, and disseminating systematic reviews of the effects of health care. Ann N Y Acad Sci 1993; 703: 156–163; discussion 163–165.
  * 2 Ahmed I, Sutton AJ, Riley RD. Assessment of publication bias, selection bias, and unavailable data in meta-analyses using individual participant data: a database survey. BMJ 2012; 344.
  * 3 Jones AP, Riley RD, Williamson PR, Whitehead A. Meta-analysis of individual patient data versus aggregate data from longitudinal clinical trials. Clinical Trials 2009; 6(1): 16–27.
  * 4 Cooper H, Patall EA. The relative benefits of meta-analysis conducted with individual participant data versus aggregated data.. Psychological Methods 2009; 14(2): 165.
  * 5 Riley RD. Commentary: like it and lump it? Meta-analysis using individual participant data. International Journal of Epidemiology 2010; 39(5): 1359–1361.
  * 6 DerSimonian R, Laird N. Meta-analysis in Clinical Trials. Controlled clinical trials 1986; 7(3): 177–188.
  * 7 DerSimonian R, Laird N. Meta-analysis in clinical trials revisited. Contemporary Clinical Trials 2015; 45: 139–145.
  * 8 Viechtbauer W. Conducting meta-analyses in R with the metafor package. Journal of Statistical Software 2010; 36(3): 1–48.
  * 9 Veroniki AA, Jackson D, Viechtbauer W, et al. Methods to estimate the between-study variance and its uncertainty in meta-analysis. Research Synthesis Methods 2016; 7(1): 55–79.
  * 10 Panityakul T, Bumrungsup C, Knapp G. On estimating residual heterogeneity in random-effects meta-regression: a comparative study. Journal of Statistical Theory and Applications 2013; 12(3): 253–265.
  * 11 Sidik K, Jonkman JN. A comparison of heterogeneity variance estimators in combining results of studies. Statistics in Medicine 2007; 26(9): 1964–1981.
  * 12 Viechtbauer W. Bias and efficiency of meta-analytic variance estimators in the random-effects model. Journal of Educational and Behavioral Statistics 2005; 30(3): 261–293.
  * 13 Hedges LV, Olkin I. Statistical methods for meta-analysis.  Academic Press .  2014\.
  * 14 Hartung J, Knapp G, Sinha BK. Statistical meta-analysis with applications. 738.  John Wiley & Sons .  2011\.
  * 15 McCulloch CE, Neuhaus JM. Generalized linear mixed models. Encyclopedia of Biostatistics 2005; 4.
  * 16 Van Houwelingen HC, Arends LR, Stijnen T. Advanced methods in meta-analysis: multivariate approach and meta-regression. Statistics in Medicine 2002; 21(4): 589–624.
  * 17 Ritz J, Demidenko E, Spiegelman D. Multivariate meta-analysis for data consortia, individual patient meta-analysis, and pooling projects. Journal of Statistical Planning and Inference 2008; 138(7): 1919–1933.
  * 18 Olkin I, Sampson A. Comparison of meta-analysis versus analysis of variance of individual patient data. Biometrics 1998: 317–322.
  * 19 Mathew T, Nordström K. On the equivalence of meta-analysis using literature and using individual patient data. Biometrics 1999; 55(4): 1221–1223.
  * 20 Whitehead A. Meta-analysis of controlled clinical trials. 7.  John Wiley & Sons .  2002\.
  * 21 Mathew T, Nordström K. Comparison of one-step and two-step meta-analysis models using individual patient data. Biometrical Journal: Journal of Mathematical Methods in Biosciences 2010; 52(2): 271–287.
  * 22 Lin DY, Zeng D. On the relative efficiency of using summary statistics versus individual-level data in meta-analysis. Biometrika 2010; 97(2): 321–332.
  * 23 Liu D, Liu RY, Xie M. Multivariate meta-analysis of heterogeneous studies using only summary statistics: efficiency and robustness. Journal of the American Statistical Association 2015; 110(509): 326–340.
  * 24 Thomas D, Platt R, Benedetti A. A comparison of analytic approaches for individual patient data meta-analyses with binary outcomes. BMC Medical Research Methodology 2017; 17(1): 1–12.
  * 25 Collaboration FS. Systematically missing confounders in individual participant data meta-analysis of observational cohort studies. Statistics in Medicine 2009; 28(8): 1218–1237.
  * 26 Kundu P, Tang R, Chatterjee N. Generalized meta-analysis for multiple regression models across studies with disparate covariate information. Biometrika 2019; 106(3): 567–585.
  * 27 Riley RD, Dodd SR, Craig JV, Thompson JR, Williamson PR. Meta-analysis of diagnostic test studies using individual patient data and aggregate data. Statistics in Medicine 2008; 27(29): 6111–6136.
  * 28 Riley RD, Lambert PC, Staessen JA, et al. Meta-analysis of continuous outcomes combining individual patient data and aggregate data. Statistics in Medicine 2008; 27(11): 1870–1893.
  * 29 Riley RD, Simmonds MC, Look MP. Evidence synthesis combining individual patient data and aggregate data: a systematic review identified current practice and possible methods. Journal of Clinical Epidemiology 2007; 60(5): 431–e1.
  * 30 Idris NRN, Abdullah MH. A study on the effects of different levels of data on the overall meta-analysis estimates. Far East Journal of Mathematical Sciences 2015; 96(1): 73.
  * 31 Yamaguchi Y, Sakamoto W, Goto M, et al. Meta-analysis of a continuous outcome combining individual patient data and aggregate data: a method based on simulated individual patient data. Research Synthesis Methods 2014; 5(4): 322–351.
  * 32 Guyot P, Ades A, Ouwens MJ, Welton NJ. Enhanced secondary analysis of survival data: reconstructing the data from published Kaplan-Meier survival curves. BMC Medical Research Methodology 2012; 12(1): 1–13.
  * 33 Riley RD, Steyerberg EW. Meta-analysis of a binary outcome using individual participant data and aggregate data. Research Synthesis Methods 2010; 1(1): 2–19.
  * 34 Donegan S, Williamson P, D’Alessandro U, Garner P, Smith CT. Combining individual patient data and aggregate data in mixed treatment comparison meta-analysis: individual patient data may be beneficial if only for a subset of trials. Statistics in Medicine 2013; 32(6): 914–930.
  * 35 Saramago P, Sutton AJ, Cooper NJ, Manca A. Mixed treatment comparisons using aggregate and individual participant level data. Statistics in Medicine 2012; 31(28): 3516–3536.
  * 36 Ravva P, Karlsson MO, French JL. A linearization approach for the model-based analysis of combined aggregate and individual patient data. Statistics in Medicine 2014; 33(9): 1460–1476.
  * 37 Sutton AJ, Kendrick D, Coupland CA. Meta-analysis of individual-and aggregate-level data. Statistics in Medicine 2008; 27(5): 651–669.
  * 38 Jackson C, Best N, Richardson S. Improving ecological inference using individual-level data. Statistics in Medicine 2006; 25(12): 2136–2159.
  * 39 Varin C, Reid N, Firth D. An overview of composite likelihood methods. Statistica Sinica 2011: 5–42.
  * 40 Yusuf S, Peto R, Lewis J, Collins R, Sleight P. Beta blockade during and after myocardial infarction: an overview of the randomized trials. Progress in Cardiovascular Diseases 1985; 27(5): 335–371.
  * 41 Lee JY, Brown JJ, Ryan LM. Sufficiency revisited: Rethinking statistical algorithms in the big data era. The American Statistician 2017; 71(3): 202–208.
  * 42 Chen X, Xie Mg. A split-and-conquer approach for analysis of extraordinarily large data. Statistica Sinica 2014: 1655–1684.
  * 43 Cheung MWL, Jak S. Analyzing big data in psychology: A split/analyze/meta-analyze approach. Frontiers in Psychology 2016; 7: 738.