Get real in individual participant data (IPD) meta-analysis: a review of the methodology.

Individual participant data (IPD) meta-analysis is an increasingly used approach for synthesizing and investigating treatment effect estimates. Over the past few years, numerous methods for conducting an IPD meta-analysis (IPD-MA) have been proposed, often making different assumptions and modeling choices while addressing a similar research question. We conducted a literature review to provide an overview of methods for performing an IPD-MA using evidence from clinical trials or non-randomized studies when investigating treatment efficacy. With this review, we aim to assist researchers in choosing the appropriate methods and provide recommendations on their implementation when planning and conducting an IPD-MA.

Introduction

The evaluation of a novel drug or intervention typically involves a series of randomized clinical trials (RCTs) where its safety and efficacy are extensively tested. Because trials are often relatively small and typically exhibit differences in study design, selection of subjects, studied outcome(s), dosage, choice of comparator intervention, and quality of the conducted research, conflicting evidence occasionally arises. As a consequence, systematic reviews have become an important tool to summarize the evidence from these trials and to generalize their conclusions beyond their specific settings.

Over the past few decades, several methods have been developed to quantify the results of a systematic review. Most of these methods adopt a meta‐analytical rationale and pool the results from individual studies by accounting for various forms of uncertainty (Sutton et al., 2009; Sutton and Higgins, 2008). Hereby, group‐level summary statistics (aggregate data (AD)) that quantify the treatments' relative efficacy or safety are retrieved from the published literature or from study authors and are subsequently synthesized into a weighted average.

Unfortunately, when synthesizing published AD, even rigorously conducted meta‐analyses can be of limited value. In particular, when there is substantial heterogeneity in estimates of relative treatment effect, a weighted average may no longer be informative in medical care. In such situations, it is important to identify whether treatment effects vary across clinical subgroups because of effect modification. Although published AD can be used for exploring modifiers of treatment effect, such approach lacks power when published summary statistics (e.g., mean age) do not vary much across studies. More importantly, the use of published AD to investigate effect modification is prone to bias because it cannot properly take subject‐level characteristics into account. This bias is also known as the “ecological fallacy” (Berlin et al., 2002). Additional problems arise when AD are not available, poorly reported, derived and presented differently across studies (for example, odds ratio versus relative risk), and more likely to be reported (and in greater detail) when statistically or clinically significant (Riley et al., 2010). For this reason, investigators increasingly embark into an individual participant data meta‐analysis (IPD‐MA) (Riley et al., 2010; Thompson, 2009; Stewart and Tierney, 2002). These meta‐analyses include the raw data from each relevant study that is (ideally) identified through a systematic review. By securing IPD of individual trials, it becomes possible to disentangle subject‐level and study‐level sources of heterogeneity in treatment effect (van Walraven, 2010; Michiels et al., 2005; Lyman and Kuderer, 2005). This, in turn, may help to explore effect modification (Simmonds and Higgins, 2007; Stewart and Tierney, 2002) and to consistently adjust for confounding variables (i.e., for differences in baseline characteristics across different treatment arms). Access to IPD may also help to improve data quality (Lyman and Kuderer, 2005; Stewart and Tierney, 2002; Lau et al., 1998), to standardize definitions and analyses (van Walraven, 2010; Lyman and Kuderer, 2005), to obtain complete follow‐up data on all randomized participants (Clarke and Stewart, 1994), to combine studies with different follow‐up times (Sud and Douketis, 2009; Lyman and Kuderer, 2005; Stewart and Tierney, 2002; Lau et al., 1998), to analyze multiple outcomes (Bujkiewicz et al., 2014; Li and Meredith, 2003), to investigate long‐term outcomes (Stewart and Tierney, 2002), and to investigate rare exposures.

Because limited guidance currently exists on dealing with challenges that are specific to IPD‐MA, we embarked on a literature review to identify state‐of‐the‐art methods and to provide recommendations on their implementation. This review aims to help systematic reviewers with a limited background in medical statistics in identifying relevant methods for their IPD‐MA.

Methods

Information sources

We conducted a literature review to identify scientific, peer‐reviewed articles that focus on methods for IPD‐MA to investigate treatment efficacy using evidence from clinical trials or non‐randomized studies. Medline and Embase databases were queried using a combination of relevant keywords in the Pubmed and Ovid platform (Supporting Information 2). The corresponding search strategy was independently reviewed by two librarians. They concluded that the search strategy did not meaningfully improve when medical subject headings terms were included or search terms were truncated. A selected set of key journals was also hand‐searched: the Journal of Research Synthesis Methods and the Journal of the Royal Statistical Society (series A, B, and C). Databases were searched from the earliest date until January 27, 2014. Finally, a last attempt to identify missed key articles was made by contacting key researchers in the field and performing cross‐reference checks. Here, relevant hits were added until March 1, 2015.

Because the search for IPD methods was limited to medical and natural science databases, this review is not intended to be comprehensive. It is, for instance, well known that several IPD‐MA methods have been developed in the field of social sciences. In general, the concepts of meta‐analysis are universally applicable across different scientific domains, whereas the nature of the data may require tailoring of statistical methods. As a consequence, a more extensive search strategy is unlikely to add relevant variations in the identified methods and may fail to yield any new insights. This effect is known as theoretical saturation (Lilford et al., 2001).

Eligibility criteria

Articles were considered eligible for inclusion if they were written in English and addressed issues related to IPD‐MA in intervention research. Key IPD‐MA issues were defined as statistical models, software routines, simulation studies, empirical comparisons, didactic examples, and best practice guidelines. Articles were excluded when IPD‐MA methods were applied to address solely a substantive clinical question, without a methodological focus. Articles for which no full text was available were also excluded.

Search results

A flowchart of the articles identified from our search is presented in Figure 1. A total of 3360 unique records were identified using the search strategy and were deemed eligible for title screening. Out of these, 3206 records were excluded because they did not have a methodological focus or did not correspond to a peer‐reviewed article (e.g., conference poster). From the remaining 154 articles, another 16 records were excluded because of missing full text or non‐English language. Finally, when assessing full text, 16 records were added after cross‐reference checks, and one article was removed because it did not adhere to inclusion criteria. A total of 153 articles were included in the review (Supporting Information 4); a complete list of the included articles can be accessed online at https://www.zotero.org/groups/wp4_‐_ipd_meta‐analysis.

Figure 1

Flow diagram of selected studies.

Results

Conceptual issues

The IPD from a single RCT are often analyzed using regression models where the observed outcome is modeled as a function of the subjects' undergone treatment (and, in many cases, other covariates of interest). The intercept term of these models then represents the study effect (e.g., baseline risk), whereas their regression coefficient represents the treatment effect. Regression models have two major advantages: First, their implementation is possible for numerous outcome types (e.g., continuous, binary, and time‐to‐event), and second, their estimated coefficients can be interpreted in a fairly straightforward manner. Unfortunately, the implementation of traditional regression models becomes challenging when a trial consists of multiple centers or when multiple trials need to be analyzed (as in an IPD‐MA). In those situations, regression models need to account for clustering of subjects within trials (or within centers) by allowing each trial (or center) to have their “own” study and treatment effect. This approach becomes infeasible when individual trials (or centers) are relatively small. Moreover, corresponding regression models no longer yield an overall estimate of relative treatment effect, casting doubt about whether and under what circumstances a certain intervention is advantageous.

Meta‐analysis attempts to overcome these issues by pooling relative treatment effects from multiple trials or centers (henceforth studies) and by borrowing information across these studies during the process. Over the past few decades, two alternate approaches have been proposed for conducting an IPD‐MA (Stewart et al., 2012; Riley et al., 2010). The most common and conceptually least complicated approach is the so‐called two‐stage approach. In this approach, the IPD are first analyzed separately in each study to produce study‐specific estimates of relative treatment effect (and possible treatment‐covariate interactions). A combined estimate is then obtained in the second step by calculating a weighted average (e.g., inverse error‐variance‐based) of the individual estimates. Hereby, one may assume that all studies share a common treatment effect and that differences between estimates from the included studies solely arise because of sampling variation (fixed effects meta‐analysis). Alternatively, it is possible to allow for heterogeneous treatment effects by incorporating variation between studies and performing a random effects meta‐analysis.

In the so‐called one‐stage approach, the IPD from all studies are analyzed simultaneously by adopting a single statistical model that fully accounts (random effects meta‐analysis), partially accounts (mixed effects meta‐analysis), or does not account (fixed effects meta‐analysis) for heterogeneity across studies. These single‐stage models are more commonly known as multilevel or hierarchical models (Simmonds et al., 2005). The implementation of some common one‐stage and two‐stage meta‐analysis model is illustrated in Supporting Information 3.

It is conventionally believed that the one‐stage and two‐stage approaches yield similar estimates of treatment effects (Stewart et al., 2012; Koopman et al., 2008b; Tudur Smith and Williamson, 2007; Steinberg et al., 1997). Indeed, it has been shown that under certain conditions, the one‐stage and two‐stage approaches may lead to equivalent results (Lyman and Kuderer, 2005; Mathew and Nordström, 2010), particularly when interest lies in estimating a single treatment effect estimate. The two‐stage approach offers several advantages favoring its implementation in the medical literature. In particular, the two‐stage approach is conceptually more intuitive and therefore requires less statistical expertise when conducting or interpreting an IPD‐MA (Stewart et al., 2012; Koopman et al., 2007). Furthermore, because the two‐stage approach analyzes each study separately, study‐specific estimates are not influenced by external information. This discrepancy may, for instance, become relevant when effect modification varies across studies and corresponding associations do not follow a well‐defined distribution (Piepho et al., 2012).

Unfortunately, the two‐stage approach is known to have little power for detecting nonlinear associations between continuous exposures and the outcome(s) of interest and for detecting treatment‐covariate interactions (Simmonds and Higgins, 2007; Simmonds et al., 2005; Tudur Smith et al., 2005b; Schmid et al., 2004). In addition, the two‐stage approach may lead to bias in pooled effects, standard errors, between‐study heterogeneity, and correlation between random effects when few studies or few participants (or events) per study are available (Debray et al., 2013; Stijnen et al., 2010; Tobías et al., 2004), when statistical models cannot fully account for follow‐up times (Poppe et al., 2011; Lyman and Kuderer, 2005; Duchateau et al., 2001; Buyse and Piedbois, 1996) or for the time between recurrent events (Haines and Hill, 2011; Jones et al., 2009). Although methodology to overcome some of these limitations has been described (Stijnen et al., 2010), the one‐stage approach offers the highest degree of flexibility for making necessary assumptions and is therefore often considered to be superior. The two‐stage approach can, however, still be useful to explore the available data, to present intermediate results or to identify key challenges when designing a one‐stage meta‐analytical model (Stewart et al., 2012).

Although IPD‐MA is often considered as gold standard, their implementation is no panacea against the limitations of meta‐analyses that are solely based on published AD. In particular, meta‐analysis of raw data may still be difficult when randomization or follow‐up procedures vary, when studies have measured different covariates or when relevant descriptions of samples, settings, and treatments are not available for every study included in the meta‐analysis (Burdett and Stewart, 2002). Furthermore, IPD‐MA may become prone to bias when IPD are only sought for a specific subset of studies (e.g., studies conducted by contacts or friends in the research field) or when the availability of IPD is related to study results (Ahmed et al., 2012).

Several recommendations have been made to accommodate for these potential sources of bias (Ahmed et al., 2012; Stewart and Tierney, 2002). First of all, meta‐analyses should be based on studies that are identified through a systematic review. Secondly, clinicians from the field should actively be involved as they may help to identify whether important studies may have been missed (Clarke and Godwin, 1998), to collect additional follow‐up data where relevant and to resolve inconsistencies between different data sources. Finally, researchers should seek IPD for all relevant studies identified and collect AD for studies where IPD could not be obtained. Methods for combining IPD and AD are discussed in Section 3.5.

Statistical models to estimate an overall summary of treatment effect

In the succeeding texts, we describe the common methodology of one‐stage IPD‐MA models. These hierarchical models can be viewed as a direct extension of traditional regression models and are also known as generalized linear mixed models (GLMMs). In contrast to traditional regression models, hierarchical regression models allow coefficients to follow a certain distribution. For instance, GLMMs may specify that the relative treatment effect between two specific interventions varies across studies according to a normal distribution. The mean of this distribution then simply represents the “average” treatment effect between the two interventions, and its variance indicates the degree of between‐study heterogeneity in treatment effect.

Consider an IPD‐MA of i = 1, …, M independent studies with k = 1, …, N subjects each. Let x be a dummy variable that indicates treatment group (treatment or control) of subject k in study i. We denote the corresponding outcome as y . This outcome can be continuous, binary, ordinal, count, and so on. The GLMM can then be stated as follows:

In this model, E(y ) denotes the expected value of y , the parameter α represents the study effect (e.g., baseline risk), and the parameter β represents the treatment effect. The link function g may take different forms depending on the type of outcome data. A detailed overview of possible implementations is discussed in Section 3.4 and also illustrated in Table 1.

Table 1

Basic statistical models for estimating overall treatment effect

Outcome type	Model type	Basic statistical model	–
Continuous	GLMM	y_ik∼Nμ_ikσi2orNμ_ikσ^2μ_ik=α_i+β_ix_ik	(L1)
Binary	GLMM	y_ik∼Bernoullip_iklogitp_ik=α_i+β_ix_ik	(L2)
Ordinal	GLMM	y_ijk∼Bernoulliq_ijklogitq_ijk=α_i+ζ_j+β_ix_ik	(L3a)
	GLMM	y_ijk∼Bernoulliq_ijklogitq_ijk=α_ij+β_ix_ik	(L3b)
Count	GLMM	y_ik∼Poissonμ_iklnμ_ik=α_i+β_ix_ik	(L4)
Time‐to‐event	Cox PH	h_ikt=h_i0texpβ_ix_ik	(L5a)
	Cox PH	h_ikt=h_0texpα_i+β_ix_ikwhereα_i∼N0σ^2orexpα_i∼σ^2Γσ^−2	(L5b)
	Cox PH	h_ikt=h_i0texpα_i+β_ix_ikwith∑α_i=0	(L5c)
	GLMM	y_ijk∼Poissonμ_ijklnμ_ijk=α_i+β_ix_ik+λ_j+lnt_ijk	(L5d)

GLMM, generalized linear mixed model; PH, proportional hazards.

Overview of statistical models for different outcome types. Each model is discussed in Section 3.4.

When specifying the study effects in model (1), each α may be taken as fixed effects (estimated separately in each study), as a common effect (so α  = α for all studies) or as random effects (α is drawn from a certain distribution). Researchers typically allow for heterogeneous study effects by estimating a separate intercept α in each study (Abo‐Zaid et al., 2013; Higgins et al., 2001; Thompson et al., 2001; Turner et al., 2000) or by specifying a random‐effects distribution on α (Broström and Holmberg, 2011; Simmonds and Higgins, 2007; Goldstein et al., 2002; Thompson et al., 2001; Turner et al., 2000). A random‐effects distribution may also be specified for β when heterogeneity in treatment effects across studies is plausible (Abo‐Zaid et al., 2013; Bowden et al., 2011; Katsahian et al., 2008; Goldstein et al., 2002; Higgins et al., 2001; Thompson et al., 2001). It is conventional and computationally efficient to assume normal distributions for the random effects (Sutton and Higgins, 2008; Thompson et al., 2001). Other distributions are also possible (Broström and Holmberg, 2011; Thompson et al., 2001) but often require more computational power and occasionally lead to convergence issues.

For instance, consider an IPD‐MA with a binary outcome where study effects and treatment effects differ across studies. The following GLMM could then be specified to estimate an overall treatment effect that is adjusted for the presence of between‐study heterogeneity:

In this model, the pooled treatment effect (i.e., the “average” treatment effect across studies) is given by μ, and the degree of between‐study heterogeneity in treatment effect is given by τ. Because the study effects are estimated separately in each study, there is no pooled study effect.

Statistical models to investigate heterogeneity in treatment effect across and within studies

In the previous section, we described common methods to produce an overall summary of treatment effect. These methods make little or no attempt to explain or investigate differences between study results (heterogeneity) across studies or to identify subgroup effects (Fisher et al., 2011; Sutton and Higgins, 2008; Simmonds and Higgins, 2007). There has been increasing interest in using meta‐analysis to go beyond estimating only the overall treatment effect. As a consequence, meta‐analyses increasingly attempt to identify modifiers of treatment effect by exploring the presence of interaction effects.

Traditional approaches such as meta‐regression and subgroup analysis investigate the presence of trial‐level interaction, that is, interaction between treatment status and a specific study‐level covariate. This covariate may represent a certain study characteristic (such as level of blinding) or a summarized subject‐level characteristic (such as mean age) (Fisher et al., 2011; Sutton et al., 2008; Tudur Smith et al., 2005b; Li and Meredith, 2003; Yamaguchi et al., 2002). For instance, Berlin et al. (2002) used meta‐regression to identify whether anti‐lymphocyte antibody induction therapy is more beneficial in patients with elevated panel reactive antibodies (PRA). Hereto, they determined the percent of patients with elevated PRA within each study and estimated their association with the corresponding treatment's effect size (Berlin et al., 2002).

Although trial‐level interactions may indicate the presence of effect modification, it has been demonstrated that such interactions have low statistical power for identifying modifiers of treatment effect and may lead to ecological (aggregation) bias. In particular, associations between aggregated values may not be representative for individual subjects (Lambert et al., 2002). In the example of induction, meta‐regression failed to identify effect modification by PRA because the level of benefit within PRA category changed across studies. A remedy against this pitfall is to use IPD and to investigate the presence of subject‐level interaction (rather than trial‐level interaction). This can be achieved by specifying an interaction term between treatment status and subject‐level covariate in model (1). Details on how to do this in one‐stage and two‐stage approaches are provided in Supporting Information 3. The inclusion of interaction terms may also help to adjust for baseline imbalances despite randomization (Higgins et al., 2001) and to increase the precision of treatment effect estimates (by including a strong prognostic factor). Estimated interaction terms must, however, be interpreted with care as they describe a mixture of trial‐level and subject‐level interactions (Higgins et al., 2001). For this reason, when a modifier of interest has sufficient variation within and across studies, researchers sometimes specify two (instead of one) interaction terms: one interaction between treatment status and the study‐specific mean of the covariate (i.e., the values that would be used in meta‐regression) and one interaction between treatment status and the study‐centered covariate values (Donegan et al., 2012; Fisher et al., 2011). In this manner, it becomes possible to quantify the presence of ecological bias.

Individual participant data meta‐analysis often assume common effects for interaction terms to restrict model complexity and facilitate estimation (Crowther et al., 2012; Teramukai et al., 2004). It is, however, possible that interaction terms are actually heterogeneous across studies. For this reason, some researchers propose to estimate subject‐level interactions in isolation (“meta‐analysis of interaction estimates”). This can be achieved by adopting a two‐stage approach where interaction effects are estimated separately within each study and then pooled in the second stage of the IPD‐MA (Fisher et al., 2011; Simmonds and Higgins, 2007). Alternatively, it is possible to adopt a one‐stage approach by assuming random‐effects distributions for the interaction effects (Simmonds and Higgins, 2007; Schmid et al., 2004). The latter strategy has been demonstrated to have the greatest power and flexibility, even when few studies are available, and is therefore generally recommended (Katsahian et al., 2008; Simmonds and Higgins, 2007). Finally, researchers may also attempt to model nonlinear interactions, for example, by using fractional polynomials (Royston and Sauerbrei, 2004), to model two‐way interactions between treatment and study‐level covariates or even to model three‐way interactions between treatment, subject‐level covariates, and study‐level covariates.

It is important to realize that practical implementation of IPD‐MA with interaction terms requires careful thought, as the range of possible analyses (i.e., options for treating the unknown parameters) could easily lead to over‐fitting and data‐dredging. This problem becomes even more apparent when multiple study‐level and subject‐level covariates are available and IPD are limited. For this reason, it is recommended to pre‐specify statistical models in a study protocol before collecting or analyzing the IPD (Higgins et al., 2001).

Modeling of specific types of outcomes

In the succeeding texts, we describe common statistical models for analyzing specific types of outcomes and describe typical estimation procedures. An overview is provided in Table 1.

Continuous

Continuous outcome data are typically modeled using an identity link function (L1). In this model, β indicates the absolute change in outcome due to treatment for patients from the i th study. If the outcome of interest represents a change from baseline (e.g., a change in depression score since inclusion by the trial), it is appropriate to adjust for baseline responses. Further discussion can be found in Supporting Information 3 and (Riley et al., 2013; Goldstein et al., 2002; Higgins et al., 2001; Thompson et al., 2001; Goldstein et al., 2000).

Binary

When modeling a binary outcome, it is common to use a logit link function (L2). In this model, exp(β ) represents the odds ratio for the treatment effect in the i th study. If β is assumed as a common effect (β  = β for all studies) or as a random effect (e.g., β  ∼ N(β, τ 2)), then exp(β) represents the odds ratio for the overall treatment effect. Examples are provided by (Thomas et al., 2014; Debray et al., 2013; Stijnen et al., 2010; Sutton et al., 2008; Thompson et al., 2001; Wakefield and Salway, 2001; Turner et al., 2000).

Ordinal

The logistic model can easily be extended to account for ordinal outcome‐type data by relating to the proportional odds model (Thompson et al., 2001; Whitehead et al., 2001). Let q denote the probability of subject k in study i having a response in category j or below. Hereby, it is assumed that the categories are ordered in terms of desirability, thus, lower categories are better. The proportional odds model is described in (L3a), where the parameter ζ represents the log odds ratios for each category cutoff. This model assumes that there is a common treatment effect β and a common study effect α within each treatment group. It is possible to relax these assumptions by allowing study effects to vary across different treatment groups (L3b). When the assumption of proportional odds is violated, generalized ordered models or partial proportional models may be used instead of the proportional odds model.

Count

Denote y as the number of events for subject i in study k. A straightforward one‐stage approach may adopt a Poisson distribution with a log‐link (L4). Several extensions have been described that for instance, allow for data with many zeros by introducing a discrete point mass (hierarchical zero‐inflated Poisson regression) (Lee et al., 2006; Yau and Lee, 2001; Hall, 2000).

Time‐to‐event

Let t denote the observed time (either censoring time or event time) for subject k in study i. The outcome y then represents an indicator that the time corresponds to an event (y  = 1) or a censoring time (y  = 0). Common multilevel survival models distinguish between heterogeneity in baseline hazard (e.g., due to differences in incidence of the event between studies given the same treatment) and heterogeneity in relative treatment effect (e.g., due to differences in treatment effect size) (Michiels et al., 2005).

The hierarchical Cox proportional hazards regression model is a popular approach for analyzing survival data from multiple studies as it requires no assumptions regarding the distribution of the baseline hazard rate (Michiels et al., 2005; Tudur Smith et al., 2005a,2005b; Sargent, 1998). This model may account for heterogeneity in baseline hazard (that is, the presence of study effects) by stratifying the baseline hazard by study (L5a). Examples of this strategy are provided in (Crowther et al., 2012; Bowden et al., 2011; Thompson et al., 2010; Katsahian et al., 2008; Tudur Smith and Williamson, 2007; Michiels et al., 2005; Tudur Smith et al., 2005a,2005b; Trikalinos and Ioannidis, 2001). Because L5a does not yield a direct estimate of the trial effect, more advanced approaches attempt to investigate this heterogeneity by introducing a frailty term α . It is then assumed that the hazards within each study are proportional to the same common baseline hazard function (L5b). Examples of this second strategy are provided in (Crowther et al., 2012; Bowden et al., 2011; Katsahian et al., 2008; Rondeau et al., 2008; Tudur Smith et al., 2005a; Yamaguchi et al., 2002; Vaida and Xu, 2000; Yamaguchi and Ohashi, 1999; Sargent, 1998; Yau and McGilchrist, 1998). Alternatively, it may be assumed that the baseline hazard functions for each study have the same shape but may have different magnitudes (L5c) (Michiels et al., 2005). Further extensions of these models may consider time‐dependent frailties (Yau and McGilchrist, 1998).

Similar to other outcome‐type data, heterogeneity in treatment effect can be investigated by specifying a random effects distribution, for example, β  ∼ N(β, τ 2) or β  = β + b where b  ∼ N(0, τ 2) (Rondeau et al., 2008; Michiels et al., 2005; Tudur Smith et al., 2005a; Yamaguchi et al., 2002; Vaida and Xu, 2000; Yamaguchi and Ohashi, 1999). In general, the use of stratified models with random treatment effects (that is, L5a or L5c using normally distributed frailty terms) has been recommended as these models maintain the within‐study structure while allowing for heterogeneity in study and treatment effects (Tudur Smith et al., 2005a; Yamaguchi et al., 2002).

Unfortunately, there are significant computational challenges in the fitting of random treatment effects to time‐to‐event outcome data (Crowther et al., 2012; Fisher et al., 2011; Thompson et al., 2010; Sargent, 1998). For this reason, researchers may choose to approximate the baseline hazard function by defining a piece‐wise constant or a spline function (Royston and Parmar, 2002; Vaida and Xu, 2000; Sargent, 1998; Yau and McGilchrist, 1998). Alternatively, it has been shown that the Cox proportional hazards model can be approximated with a Poisson GLMM by splitting follow‐up time t (and event indicator y ) into j narrow time intervals of fixed length (L5d) (Crowther et al., 2012). The parameter λ then represents the baseline hazard rate during the j th time interval, and β is once again the log hazard ratio for the treatment effect. An additional advantage of this approach is that the proportional hazards assumption can be relaxed, for example, by replacing the Poisson with a log‐gamma distribution or by pre‐specifying the baseline function h 0(t) in L5b (Barrett et al., 2012; Siannis et al., 2010).

Other

Several one‐stage models have been proposed for IPD‐MA where the outcome belongs to another type of non‐normally distributed data (Thompson et al., 2001), recurrent events (Haines and Hill, 2011), repeated measurements (Jones et al., 2009), or multivariate responses (Kim et al., 2013; Riley et al., 2008; Yamaguchi et al., 2002).

Meta‐analysis of individual participant data combined with aggregate data

Subject‐level data are often unavailable for all relevant studies (Ravva et al., 2014; Donegan et al., 2013; Saramago et al., 2012; Jansen, 2012; Riley et al., 2008; Sutton et al., 2008; Riley et al., 2007). Therefore, published AD may represent an additional source of evidence. Including published AD in an IPD‐MA is strongly recommended when studies providing IPD systematically differ from studies for which IPD were not obtained, as this may lead to data availability bias or reviewer selection bias (Ahmed et al., 2012). But even when there is no indication of such bias, inclusion of AD in an IPD‐MA for studies of which IPD are lacking may be considered to increase the statistical power for detecting treatment effects or treatment‐covariate interactions (Donegan et al., 2013; Jansen, 2012; Saramago et al., 2012). In general, there are three alternate approaches for combining IPD and AD in a meta‐analysis:

Meta‐analysis of AD

In this approach, the available IPD are first reduced to AD (Ravva et al., 2014; Riley et al., 2007) and then pooled into a weighted average as in two‐stage meta‐analysis. When investigating the presence of effect modification, caution is warranted to avoid the introduction of ecological bias (Signorovitch et al., 2012; Jackson et al., 2008; Riley et al., 2008, 2007; Tudur Smith et al., 2005b; Wakefield and Salway, 2001).

Meta‐analysis of reconstructed IPD

It is possible to reconstruct IPD from published aggregate information for binary (based on 2 × 2 tables), ordinal (based on the number of responses within each treatment category), and survival type data (based on Kaplan–Meier survival curves) (Guyot et al., 2012; Riley and Steyerberg, 2010; Riley et al., 2007). In particular, this summary information can be used to restore subject‐level information on treatment group and outcome status (or time until event). Afterwards, (traditional) one‐stage models can be used for evidence synthesis. Unfortunately, reconstructed IPD do not have subject‐level information on covariates and may therefore lead to ecological bias when investigating the presence of effect modification.

Hierarchical‐related regression

This approach directly combines the likelihood from each data source (IPD or AD), for example, using Bernoulli (IPD) and binomial (AD) distributions for modeling the individual responses and trial‐specific summaries of binary outcome data (Ravva et al., 2014; Saramago et al., 2012; Jansen, 2012; Sutton et al., 2008; Jackson et al., 2006, 2008; Goldstein et al., 2002). In this manner, it ensures that all studies contribute toward estimation of the overall treatment effect and study‐level covariates and that only IPD are used to estimate subject‐level covariate effects. Hierarchical related regression is also known as shared parameters models (Donegan et al., 2013) and can, for instance, be achieved by fitting a multilevel model that includes a dummy variable indicating whether a data source relates to IPD or AD (Riley et al., 2008, 2007; Goldstein et al., 2000).

Unfortunately, the incorporation of AD in an IPD‐MA does not necessarily improve the validity of estimated treatment effects and interactions. In particular, amongst other issues, the exclusion of subjects post randomization (Tierney and Stewart, 2005) as well as the presence of limited follow‐up data and informative drop‐out may introduce bias in published AD. Because researchers may often find themselves between the Scylla and Charybdis when deciding to include published AD in an IPD‐MA, it has been recommended to assess the sensitivity of meta‐analysis results to the inclusion of AD (Riley and Steyerberg, 2010; Riley et al., 2008, 2007). This can, for instance, be achieved by performing an IPD‐only meta‐analysis and meta‐analysis that combines IPD and AD. An additional advantage of sensitivity analyses is that they may help readers without (advanced) statistical expertise to interpret results in function of the underlying assumptions.

Multiple treatment comparisons

In a different article of this series, we discuss in depth the use of network meta‐analysis (NMA) models for synthesizing information from multiple treatments across different studies (Efthimiou et al., 2015). These models allow the estimation of the relative treatment effect between many (rather than two) specific interventions, even when these interventions have never been compared head‐to‐head. The implementation of NMA models is, however, less straightforward than traditional meta‐analysis models and requires studies to have similar covariate distributions for important effect modifiers (Signorovitch et al., 2010). For this reason, it is important to correct for baseline imbalance between studies (Ali et al., 2013; Jansen, 2012) and explore the presence of treatment‐covariate interactions (Donegan et al., 2012; Jansen, 2012; Signorovitch et al., 2012). This is ideally achieved using IPD.

The specification of IPD‐NMA models can be viewed as an extension of regression model (1) where a dummy variable is included for each (except one) treatment (Higgins et al., 2001). It is, however, more common to describe NMA models as a function of consistency equations, which simplifies the specification of interaction effects (Donegan et al., 2013; Jansen, 2012; Saramago et al., 2012). The use of IPD in NMA has been advocated even if IPD are only available for a fraction of the studies (Jansen, 2012). Recently, several NMA models have been proposed that allow the combined use of IPD and AD in a joint analysis (Donegan et al., 2013; Jansen, 2012; Saramago et al., 2012). These models were shown to have improved network consistency when compared to NMA models that were solely based on AD. Similar to AD‐NMA, the quality of evidence resulting from an IPD‐NMA strongly depends upon the quality and choice of included studies, the complexity and credibility of statistical models, and the transparency of reporting. A questionnaire for assessing the relevance and credibility of NMA studies has recently been reported by the International Society for Pharmacoeconomics and Outcomes Research and may help to improve the quality of IPD‐NMA studies (Jansen et al., 2014).

Cross‐design synthesis

Most trials are conducted to establish efficacy and safety of a single treatment in a specific study setting, and therefore it is often difficult to translate this evidence to “real‐world” effectiveness where different concerns may be at play (Prevost et al., 2000; Stuart et al., 2011). For this reason, there is a growing interest to overcome the weaknesses of individual study designs by carefully combining different types of evidence (Droitcour et al., 1993). In particular, intervention research might benefit by the inclusion of non‐randomized study data such as from (case‐)cohort studies, case‐control studies, or patient registries (Higgins et al., 2013; Kaizar, 2011; Peters et al., 2005; Prevost et al., 2000; Li and Begg, 1994). This strategy is more commonly known as cross‐design synthesis or generalized synthesis of evidence (Sutton and Higgins, 2008; Sutton et al., 1998).

Because non‐randomized studies of interventions (NRSI) are inherently prone to confounding, there is substantial controversy about the magnitude and interpretation of potential differences between the results from RCTs and NRSI (Reeves et al., 2003; Valentine and Thompson, 2013; Schmidt et al., 2013; Golder et al., 2011; Ioannidis et al., 2001; Concato et al., 2000). Whether or not it is valid to combine RCT and NRSI in IPD‐MA and under what conditions is an ongoing debate in the literature (Higgins et al., 2013; Reeves et al., 2003). So far, cross‐design synthesis has only been advocated when the likelihood of conducting randomized trials is low (e.g., due to ethical concerns or practical limitations) or when research questions relate to unintended, rare (often harms), or long‐term outcomes (benefits or harms). The justification of including NRSI may further increase when the derived estimates of relative treatment effects are deemed to be resistant against bias. This becomes more plausible when there are no systematic differences between the groups being compared (e.g., in terms of care provided, exposure to other factors, baseline characteristics, or in terms of how outcomes are determined).

To improve the validity of meta‐analyses incorporating evidence resulting from NRSI, several suggestions have been made. First, researchers should critically appraise the risk of bias in available studies prior to cross‐design synthesis and carefully decide whether inclusion is appropriate (Berger et al., 2014; Higgins et al., 2013; Reeves et al., 2003; Wells et al., 2013). Subsequently, researchers should routinely account for confounding factors within studies. Finally, researchers should explore potential sources of between‐study heterogeneity and compare the results between randomized studies and NRSI (Valentine and Thompson, 2013).

Dealing with missing data

Missing data are a common problem when analyzing IPD from any (single) study and are typically classified into missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR) (Rubin, 1987). These types correspond to situations where the probability of data being missing is completely independent of observed and unobserved factors (MCAR), depends only on the observed data (MAR), or also depends on unobserved data (MNAR) (Donders et al., 2006; Schafer, 1999; Rubin, 1987). When data are missing within a single study, it is generally recommended to apply multiple imputation strategies by adopting a statistical model that assumes MCAR or MNAR. We can distinguish the following missing data scenarios in an IPD‐MA (Sutton et al., 1998):

missing data for some subjects in one or more studies,

missing data for all subjects of one or more studies (e.g., variables that have not been measured, outcomes that are missing, or missing study‐level covariates), and

entire study missing (e.g., when study authors are unable to share IPD).

Because the association between observed and missing variables may differ across studies, it has been recommended to allow for between‐study heterogeneity during imputation. This can, for instance, be achieved by imputing IPD studies separately (two‐stage imputation) (Burgess et al., 2013; Koopman et al., 2008a). Two‐stage imputation is obviously not feasible if data are missing for all subjects within one or more studies. In such scenarios, it is desirable to use the so‐called one‐stage imputation models that borrow information across the available studies (Jolani et al., 2015; Burgess et al., 2013; Resche‐Rigon et al., 2013; White et al., 2008; Yucel, 2008, 2011).

Software and estimation techniques for individual participant data meta‐analysis

Several software packages exist for conducting an IPD‐MA and adopt a frequentist or a Bayesian estimation framework (Tables 2 and 3). Frequentist methods typically maximize the likelihood function of the meta‐analysis model using expectation–maximization (Vaida and Xu, 2000), Newton–Raphson, or Fisher scoring algorithms. Unfortunately, maximization of the (log‐)likelihood function becomes problematic when there are few studies in the meta‐analysis or included studies are small (Crowther et al., 2012; Li et al., 2011; Thompson et al., 2001; Higgins et al., 2001). For this reason, it is often recommended to adopt more advanced estimation techniques such as restricted (or residual) maximum likelihood estimation (Bowden et al., 2011; Higgins et al., 2001), penalized maximum likelihood estimation (Rondeau et al., 2008), penalized quasi‐likelihood (Thompson et al., 2001), extended quasi‐likelihood estimation, penalized partial likelihood estimation, or semi‐parametric‐penalized likelihood estimation. Unfortunately, these techniques are not widely implemented to investigate nonlinear outcomes (e.g., binary, ordinal, or count data).

Table 2

Software packages for fitting meta‐analysis models

Software	Package	Characteristics	Used in
WinBUGS, JAGS, Stan, OpenBUGS	–	Fitting of one‐stage and two‐stage meta‐analysis models using Bayesian Markov chain Monte Carlo (MCMC)	[26, 32, 48, 52, 67, 81, 86, 95, 98, 118, 126, 146, 147]
NONNEM	–	Fitting of GLMM	[87]
R, S‐plus*	–	Unspecified	[9, 103, 104]
–	ecoreg	Estimation of individual‐level covariate‐outcome associations using AD (“ecological inference”) or a combination of AD and individual participant data (IPD) (“hierarchical‐related regression”)	[52]
–	glmmML	Fitting of GLMM using ML. Allows for non‐normal distributions in the specification of random intercepts.	[11]
–	hglm	Fitting of GLMM where the random effect may come from a conjugate exponential‐family distribution	–
–	lme4	Fitting of GLMM using ML or REML (for mixed linear models only)	[11, 30, 56, 65, 110, 153], SI3
–	MASS	Fitting of GLMM using penalized quasi‐likelihood (PQL)	–
–	mvmeta	Meta‐analysis and meta‐regression of AD (two‐stage meta‐analysis)	[30], SI3
–	nlme	Fitting of linear mixed‐effects models using ML or REML	[89, 122]
–	survival	Fitting of Cox PH and mixed effect survival models using penalized partial likelihood estimation (PPL)	[44, 122, 132, 133]
–	frailtypack	Fitting of frailty models using semi‐parametric‐penalized likelihood (SPL)	[97]
–	coxme	Fitting mixed effects Cox PH models	–

GLMM, generalized linear mixed model; AD, aggregate data; PH, proportional hazards; REML, restricted (or residual) maximum likelihood; ML, maximum likelihood.

References are provided in Supporting Information (SI) 4.

An overview of software packages for two‐stage meta‐analysis can be found on http://cran.r‐project.org/web/views/MetaAnalysis.html.

Table 3

Software packages for fitting meta‐analysis models

Software	Package	Characteristics	Used in
SAS	–	Unspecified	[106, 122, 132, 133]
–	PROC GLIMMIX	Fitting of GLMM using penalized quasi‐likelihood (PQL)	[82, 129]
–	PROC GLM	Fitting of GLMM using MOM	[148]
–	PROC LOGISTIC	Fitting of mixed nonlinear models (binary/ordinal/nominal responses) using ML	[147]
–	PROC MIXED	Fitting of mixed linear models using ML, REML, or MOM. Can also perform two‐stage meta‐analysis.	[48, 57, 93]
–	PROC NLMIXED	Fitting of mixed nonlinear models using (approximated) ML	[95, 114, 148]
Stata	–	Unspecified	[3, 4, 45]
–	gllamm	Fitting of GLMM using ML	–
–	mvmeta	Meta‐analysis and meta‐regression of aggregate data (AD) (two‐stage meta‐analysis)	[5]
–	REGOPROB2	Fitting of random effects generalized‐ordered probit models	[26, 95, 122]
–	stmixed	Fitting of flexible parametric survival models with mixed effects	–
–	XT	Fitting of GLMM	–
MLwiN, MLn	–	Fitting of GLMM and survival models using ML, REML, and EM	[41, 48, 122, 126, 136, 147]
FORTRAN, C	–	Unspecified	[60, 97, 137]

GLMM, generalized linear mixed model; REML, restricted (or residual) maximum likelihood; ML, maximum likelihood; MOM, method of moments; EM, expectation maximization.

References are provided in Supporting Information 4.

The Bayesian framework combines the likelihood function with prior information (Yamaguchi et al., 2002; Higgins et al., 2001; Sargent, 1998; Larose and Dey, 1997) and is typically implemented within WinBUGS, JAGS, OpenBUGS, or Stan. The use of prior information may help to overcome convergence issues when studies are small or few studies are available (Table 4). Semi‐informative prior distributions may, for instance, be used when model parameters are unidentifiable due to a lack of data (Jackson et al., 2008). Alternatively, informative prior distributions may be used, for example, when relevant information can be borrowed from historical data (Viele et al., 2014). An additional advantage of the Bayesian framework is that it allows more flexibility in choosing and combining likelihood functions that are appropriate for the data that are being analyzed. As such, the use of Bayesian meta‐analyses may sometimes help to improve the credibility of estimated treatment effects and between‐study heterogeneity. The Bayesian framework, however, requires careful specification of prior distributions, investigation of convergence, and the drawing of inferences from a large simulated sample. Because it tends to be cumbersome and requires substantial statistical expertise, Bayesian IPD‐MA is not pursued by many researchers.

Table 4

Overview of prior distributions used within the Bayesian framework

–	Prior distribution	95% prior belief	Used in
Parameters for fixed effects (regression coefficients)
Non‐informative	β ∼ N(0, 106)	[0.00; + ∞) for the odds ratio	[98, 118]
–	β ∼ N(0, 105)	[0.00; + ∞) for the odds ratio	[26]
–	β ∼ N(0, 104)	[0.00; + ∞) for the odds ratio	[32, 48, 70, 86, 95, 126, 147]
–	β ∼ N(0, 10)	[0.00; 494] for the odds ratio	[86]
Weakly informative	β ∼ N(0, 1.47)	[0.20; 5.03] for the odds ratio	[52]
Parameters for between‐study standard deviation
Non‐informative	τ 2 ∼ Γ − 1(10− 3, 103)	[1.2 × 107, + ∞) for τ	[67]
–	τ 2 ∼ Γ − 1(10− 3, 10− 3)	[3.2 × 104; + ∞) for τ	[48, 81, 147]
–	τ 2 ∼ U(0.01, 100)	[1.62; 9.87] for τ	[70]
–	τ ∼ N(0, 10)I(0)	[0.30; 22.48] for τ	[146]
–	τ ∼ U(0, 10)	[0.25; 9.76] for τ	[32, 118]
Weakly informative	τ ∼ U(0, 5)	[0.12; 4.88] for τ	[86]
–	τ ∼ U(0, 2)	[0.05; 1.96] for τ	[98]
–	τ ∼ N(0, 1)I(0)	[0.03; 2.24] for τ	[26, 95]
–	τ ∼ U(0, 1)	[0.02; 0.98] for τ	[126]
–	τ ∼ N(0, 0.125)I(0)	[0.01; 0.80] for τ	[86]
–	τ 2 ∼ Γ − 1(1, 10− 2)	[0.05; 0.62] for τ	[52]
–	τ ∼ N(0, 0.033)I(0)	[0.01; 0.41] for τ	[86]

References are provided in Supporting Information 4.

Concluding remarks

Individual participant data meta‐analyses offer numerous advantages over meta‐analyses that are solely based on published AD and are therefore considered as a gold standard in evidence synthesis. Nevertheless, the presence of publication bias, data accessibility bias, or reviewer selection bias may still hamper or invalidate the results of an IPD‐MA (Ahmed et al., 2012), and synchronizing between different sources of IPD is not always straightforward. For this reason, it is not recommended to conduct IPD‐MA on an ad hoc basis (e.g., without systematic review). Because the implementation of an IPD‐MA requires additional efforts and statistical expertise, researchers should carefully assess whether the potential advantages (e.g., increased power, reduced bias, and investigating interaction and subgroup effects) outweigh the extra efforts involved. Finally, researchers should be aware that IPD‐MA may have similar issues to meta‐analyses of published AD and are no panacea against poorly designed and conducted primary research.

Funding

The work leading to these results has received support from the Innovative Medicines Initiative Joint Undertaking under grant agreement number 115546, resources of which are composed of financial contribution from the European Union's Seventh Framework Programme (FP7/20072013) and EFPIA companies' in‐kind contribution

Supporting info item

Click here for additional data file.