Aggregate-data estimation of an individual patient data linear random effects meta-analysis with a patient covariate-treatment interaction term

Individual patient-data meta-analysis of randomized controlled trials is the gold standard for investigating how patient factors modify the effectiveness of treatment. Because participant data from primary studies might not be available, reliable alternatives using published data are needed. In this paper, I show that the maximum likelihood estimates of a participant-level linear random effects meta-analysis with a patient covariate-treatment interaction can be determined exactly from aggregate data when the model’s variance components are known. I provide an equivalent aggregate-data EM algorithm and supporting software with the R package ipdmeta for the estimation of the “interaction meta-analysis” when the variance components are unknown. The properties of the methodology are assessed with simulation studies. The usefulness of the methods is illustrated with analyses of the effect modification of cholesterol and age on pravastatin in the multicenter placebo-controlled regression growth evaluation statin study. When a participant-level meta-analysis cannot be performed, aggregate-data interaction meta-analysis is a useful alternative for exploring individual-level sources of treatment effect heterogeneity.


INTRODUCTION
Subgroup analyses are important for identifying possible sources of treatment effect heterogeneity (Wang and others, 2007). A formal way to assess subgroup effects is through the estimate of a treatmentcovariate interaction in a regression analysis, the 'covariate' being the variable that defines the subgroups of interest (Brookes and others, 2001). The power and generalizability of this estimate can be improved using meta-analysis, where data from multiple compatible trials is combined to obtain a pooled estimate of the subgroup effect (Rothwell, 2005). However, available methods for studying subgroup effects with meta-analysis of aggregate data have important practical and inferential limitations.
Two approaches for estimating a treatment-covariate interaction with an aggregate-data meta-analysis have been discussed in the literature (Fisher and others, 2011). In 'meta-analysis of interactions' the metaanalyst obtains study estimates of a given treatment-covariate interaction and pools them across studies using a standard fixed or random effects model (Simmonds and Higgins, 2007). The major drawback 274 S. A. KOVALCHIK of this method is that covariate-treatment interactions are not typically included in primary clinical trial reports. The second approach is meta-regression, which regresses a measure of treatment effect on aggregated covariates, typically the covariate means for the trial sample. Meta-regression is more practical than meta-analysis of interactions because the covariate sample characteristics and estimates of treatment effect are routinely reported in the primary studies, but the ecological association between a covariate and treatment is a potentially biased and inefficient estimate of the individual-level association (Berlin and others, 2002;Lambert and others, 2002).
Compared to meta-regression and meta-analysis of interactions, a treatment-covariate interaction estimated from a regression analysis of complete individual patient data (IPD), where both between-and within-trial information is included, is more powerful and less vulnerable to bias (Schmid and others, 2004;Riley and others, 2007;Koopman and others, 2008). Yet, owing to the difficulties of collecting IPD, a regression analysis of the original trial data will generally not be feasible (Riley and others, 2008). In fact, <5% of current published meta-analyses in the medical sciences are based on patient-level data (Kovalchik, 2012). For this reason, it would be useful to establish whether aggregate statistics are sufficient to obtain the estimates of an IPD meta-analysis with individual-level covariates. Olkin and Sampson (1998) showed that the treatment contrasts in a fixed effects IPD metaanalysis are equivalent to the aggregate-data treatment effect estimates in a fixed effects meta-analysis. Mathew and Nordström (1999) later showed that this result holds when studies have heterogeneous residual variance. More generally, Lin and Zeng (2010) found that the efficiency loss in estimating the regression parameters for standard fixed effects parametric or semiparametric IPD models using meta-regression is small when maximum likelihood estimation is used and there is minimal residual between-trial heterogeneity. No study has considered to what extent these findings can be extended to models with random effects and individual-level associations.
In this paper, I show that the maximum likelihood estimates for an IPD linear random effects model with a patient-covariate-treatment interaction term for a categorical covariate can be determined exactly from published clinical trial data when the model's variance parameters are known. I then provide variance estimators and an optimization algorithm for estimating all of the model parameters when the variance components are unknown. The statistical properties of the estimation methodology are investigated in simulation studies. The methods are demonstrated with an aggregate-data meta-analysis of real clinical trial data to quantify how age and cholesterol modify the effectiveness of the statin pravastatin for the prevention of cardiovascular disease progression.

MODEL
A general form for a linear mixed effects meta-analytic model with continuous outcome is (2.1) Each quantity of (2.1) is a study-specific vector for j = 1, . . . , J studies. The Y j is the n j × 1 column vector of continuous outcomes based on i = 1, . . . , n j subjects, X j is an n j × p matrix of covariates, and Z j is a subset of X j containing q variables with random effects. The study random effects α j follow a shared multivariate normal distribution with unrestricted covariance-variance D, α j ∼ N (0, D). The within-trial errors are iid normal with study-specific variance, j ∼ N (0, σ 2 j I n j ), where I m is the m × m identity matrix.
The specific case of model (2.1) that will be considered in this paper is a random effects interaction model with a categorical covariate. Given t j treated and c j control subjects in the jth study, define the treatment group indicator, G j = 1 t j 0 c j , and the diagonal matrix with G j elements on the diagonal as D(G j ).

275
Suppose that the categorical covariate of interest has K categories and let the dummy vector indicating membership in the kth subgroup be S j (k). In terms of these quantities, the interaction model's fixed and random effects design matrices are The model specified by (2.2) incorporates trial heterogeneity in baseline and treatment effects α j through the bivariate normal random effects having covariance-variance D = τ 0 τ 01 τ 01 τ 1 . The covariate and interaction effects are constant for all trials.
To see how the individual-level and standard aggregate-level approaches differ, consider the metaregression model to investigate the treatment-covariate interactions of model (2.2), Here,θ j is the study estimate of treatment effect. For the linear model with continuous outcome this iŝ θ j =ȳ 1 j −ȳ 0 j , the mean difference in outcome between the treatment and control group. The parameter θ j is a study random treatment effect, , whose variance is considered known. Importantly, meta-regression does not include individual-level information about covariates. Instead, covariates enter as the sample mean of their corresponding variable in the IPD model, in this case, the proportion of individuals in the kth subgroup,s(k) j . Although one could consider using the ecological parameters γ for inference about individual-level treatmentcovariate interactions, these quantities are vulnerable to bias and inefficiency (Berlin and others, 2002;Lambert and others, 2002;Fisher and others, 2011). Consequently, meta-analysis experts and software developers generally recommend that the covariates of meta-regression be limited to study-level variables Viechtbauer, 2007).

Estimation with known variances
Maximum likelihood (Hartley and Rao, 1967) is a common method used to estimate the parameters of the linear mixed effects regression model. Given the variance components σ j and D, the maximum likelihood estimates (MLEs) for the fixed and random effects of model (2.1) are (Laird and Ware, 1982), The variance for this estimator is The empirical Bayes estimator for the random effects iŝ and has the variance In the supplementary material (Section A available at Biostatistics online), I show that aggregate data is sufficient to compute the above MLEs and their variances under the "interaction meta-analysis" (2.2) with known variance components. However, in general, the variance components will not be known. Previous work of Kovalchik and Cumberland (2012) derived an aggregate-data estimate for the error of the patientlevel interaction term in a general linearized random effects meta-analysis. The motivation in that study was to estimate what the power of the IPD test of interaction would be if the patient-level data of the primary studies were available. In the following subsections, I build on this work by providing an EM algorithm to estimate all parameters of the IPD random effects model defined by (2.2) from aggregate data when the model's variance components are unknown.

Estimation with unknown variances
Given the parameters α j and β, maximum-likelihood estimates for the variance components of model (2.2) areD = J j=1α jα j /J andσ 2 j =ˆ jˆ j /n j , withˆ j = Y j − X jβ − Z jα j . We see thatD can be reconstructed exactly from the aggregate-data form of theα j . Because the maximum likelihood estimates in a random effects model have been recognized to underestimate the target parameters, the denominator J − 1 is used in place of J .
The estimator forσ 2 j can be formed from the sample variances s 2 k j of the continuous outomce y i j within each of the K subgroups. Variance decomposition yields the expression whereμ k j is the study-specific mean estimate for the kth treatment-covariate subgroup. The above expressions are equivalent to the E-step of the EM algorithm for the linear mixed model of (2.2) (Laird and Ware, 1982). When the subgroup-within-treatment sample variances s 2 k j are not available, an estimate for the s 2 k j can be obtained from the overall trial sample variance s 2 Aggregate-data estimation of an individual patient data meta-analysis 277

Optimization algorithm
This section outlines a complete aggregate-data EM algorithm for the interaction meta-analysis based on the M-step described in Section 3.1 and the E-step described in Section 3.2. The algorithm estimates the effects and variance components of model (2.2) using the following sufficient aggregate data: the sample sizes (t k j and c k j ) and outcome means (ȳ t . j andȳ c . j ) in each subgroup and treatment arm; the sample variance (s 2 j ) for the outcomes of each trial.
The subscript '.' indicates the index over which the mean is taken.

SIMULATION STUDY
The theoretical results have shown that grouped statistics are sufficient to construct an equivalent EM algorithm for the IPD interaction meta-analysis. Simulation studies were performed to assess the performance of the optimization algorithm for trial number and sample size settings that are typical of current meta-analyses of randomized trials (Kovalchik, 2012;Engels and others, 2000). In these analyses, y i j were generated under model (2.2) for fixed β, σ 2 j = 1 2 + 1/j, and D = . Separate simulations were run with a two-group and three-group categorical covariate. Covariate group assignment was drawn from a multinomial distribution with an equal probability for each group. Meta-analyses of {10, 15, 20} trials, 200 total subjects in each trial and equal allocation to treatment and control groups were considered. The aggregate-data used for estimation were the outcome means and sample sizes in each treatment arm and covariate subgroup, and the trial sample variances for the y i j outcomes. Summaries of the bias, standard error, mean standard error for the aggregate-data estimates of β and D were based on 1000 iterations of each simulation scenario.
The results of the simulation studies for the two-category covariate show that the aggregate-data approach yields unbiased fixed effects for all of the trial number and sample size conditions considered (Table 1). Overall, the mean standard errors (MSEs) were in good agreement with the empirical standard errors (ESE). The estimates for the random effects variance were accurate for all conditions. The MSEs for the variance components are not reported because variance estimates forD are not considered in this work.  The supplementary material (Section B available at Biostatistics online) shows that the statistical properties for the two-category covariate case were unchanged when a three-category covariate was used.

APPLICATION
The Regression Growth Evaluation Statin Study (REGRESS) (Jukema and others, 1995) was a placebocontrolled, multicenter study to assess the effect of treatment with the HMG-CoA reductase inhibitor pravastatin on progression and regression of coronary atherosclerosis in male patients with a serum cholesterol between 4 and 8 mmol/L (155 and 310 mg/dL). The primary endpoint of the study was the change from baseline in the average mean segment diameter of the coronary artery after 2 years of treatment. Less change is more indicative of treatment benefit. The study found that pravastatin was associated with a statistically significant mean difference of 0.04 mm in the mean change in segment diameter.
To demonstrate the paper's methodology, I investigate how baseline cholesterol and age influenced responsiveness to pravastatin using the grouped data reported in Table 2. The data have been permuted within the treatment group and the study in order to be used for purely demonstrative purposes. The optimization algorithm of Section (3.3) was applied to estimate the aggregate-data IPD interaction model with a two-group categorization of baseline cholesterol ( 6, > 6 mmol/L) and a three-group categorization of age (< 50, 50-59, 60 years). To assess the performance of the methods, these estimates were compared Table 2

. Trial and subgroup summary statistics on change in mean segment coronary artery diameter (in mm) used to estimate interaction meta-analysis in the REGRESS study.
Cholesterol ( with the results of complete IPD analyses, which were fit with the R package lme4 (Pinheiro and Bates, 2000). The aggregate-data and IPD meta-analyses yielded very similar results (Table 3). Each found that a baseline cholesterol level less than or equal to 6 mmol/L is associated with a 0.04 mm greater reduction with pravastatin treatment, suggesting that subjects with higher cholesterol tend to have more treatment benefit from the statin. However, this difference was not statistically significant. The results for age, which are provided as supplementary material (Section C available at Biostatistics online), showed similar consistency between methods. Aggregate-data estimation of an individual patient data meta-analysis

281
The random effects of the interaction meta-analysis can be useful for evaluating intertrial heterogeneity. For the cholesterol analysis, the left panel of Figure 1 shows the study-specific deviations from the population intercept; the right panel shows deviations from the population treatment effect.

DISCUSSION
The methods provided in the present paper broaden the findings for linear models discussed by Olkin and Sampson (1998), Mathew and Nordström (1999), and Lin and Zeng (2010). In comparing the properties of patient-level and aggregate-data analyses, these authors focused on study-level effects and fixed-effect conditions. The present paper has shown that the estimates of an IPD random effects metaanalysis with an individual-level interaction term, for any categorical covariate, can be reliably estimated from published aggregate data. These are important extensions of earlier work because between-trial heterogeneity is common  and there is growing interest in exploring to what extent patient-level modifiers can explain trial differences in treatment effect.
There are a number of advantages to gathering participant level data that no current aggregate-based method can substitute. These include the ability to ensure uniformity in variable definitions and eligibility criteria and to diagnose model fit. It is because of these strengths that the presented methodology is not intended to compete with IPD meta-analysis but rather to provide a good aggregate-data approximation to a specific kind of IPD interaction model when participant-level data cannot be obtained. Owing to the heightened concern over individual-data security that has come with the genomic era (Homer and others, 2008), there is likely to be a growing demand for increasingly sophisticiated aggregate-data methods in the future.
The approach introduced here has important advantages over available aggregate-data methods. Unlike meta-regression, it allows the assessment of individual-level associations. And, in contrast to meta-analysis of interactions, there are fewer barriers to obtaining the data necessary to perform the analysis. Metaanalysis of interactions requires that primary investigators perform a regression analysis with an interaction term and make their findings available. The data needed for the method presented here can be obtained by stratifying the standard outcome analysis recommended by the CONsolidated Standards of Reporting Trials (CONSORT) (Schulz and others, 2010) on the covariate of interest.
The IPD random effects model considered in this paper was limited to a continuous outcome and a single categorical covariate with fixed covariate and interaction effects. Without requiring additional data from the primary studies, it is unclear how the paper's methods could be extended to include multiple covariates or to allow random covariate effects. Though it would be a challenge to extend the aggregate-data approach in this direction, the work of Lin and Zeng (2010) suggests that, with some additional assumptions, the single interaction model could be generalized to standard parametric and semiparametric models of binary and time-to-event outcomes.
In clinical trial research, the aims to personalize treatment and avoid spurious findings from subgroup analyses have often seemed irreconcilable (Lagakos, 2006). The approach in this paper offers a practical solution to this dilemma for trials with a continuous endpoint. The methodology and supporting software give clinical researchers a reliable tool to approximate an individual level assessment of effect modification with aggregate data from multiple trials of the same treatment. To facilitate the utility of this tool, clinical trialists should be encouraged to make stratified analyses of the study endpoint available upon request or as supplementary material.