Keywords: meta-analysis, multivariate analysis, multivariate meta-analysis, non-linear, splines

2.1. Temperature–health associations

The health effects of exposure to extreme temperatures have been frequently investigated, especially regarding the association with high temperatures and heat waves [14]. Analyses are typically based on a time-series design, comparing averaged daily temperature with aggregated daily outcomes, often mortality or morbidity counts occurring in a specific city. Most studies have reported an increased risk for both cold and hot temperatures and have used a variety of functions to describe such non-linear association in regression models [15]. The association has shown a strong geographical pattern, linked with climatological, socio-economic and demographic factors of each city [16]. This issue has prompted the use of multi-city analyses, where the results from different cities are pooled together, and explanation for differences is sought through reference to study-specific characteristics such as those listed previously. Individual data, corresponding to series of daily observations from different cities, are usually available, and the analytical setting corresponds to the two-stage framework detailed in Section 1. Here the termsstudy andcity may be used interchangeably to refer to second-stage units.

Different analytical approaches have been proposed for modelling the health effects of temperature in multi-study settings [13]. Traditional methods rely on the simplification of the dependency, assuming a linear association beyond a threshold [17]. An alternative uses non-linear representations such as splines or polynomials but reduces the estimand to a simple comparison between two specific temperatures [18]. Both these approaches are motivated by the need to summarize the study-specific associations in a single parameter for meta-analysis. Such simplifications are, however, based on strong assumptions about the exposure–response shape or, alternatively, only provide a partial picture of possibly complex dependencies. This supports the development of methods allowing the meta-analysis of truly non-linear multi-parameter relationships.

2.2. Data

The multi-city time-series data used in this analysis were collected as part of the National Morbidity, Mortality and Air Pollution Study (NMMAPS) (http://www.ihapss.jhsph.edu). This publicly available database contains, among other information, daily series of mortality counts and weather and pollution measurements totalling 5114 observations for the period 1987–2000 in each of 108 cities in the USA. For the purpose of this illustrative example, we restrict the analysis to the 20 cities illustrated inFigure 1, situated in the industrial Midwest region, as defined in the original study.

In addition, the database includes city-level measures of several variables on geographical, climatological, demographic and socio-economic characteristics. A summary of average non-accidental mortality, temperature and other study-specific variables across the 20 cities is reported inTable I.

2.3. Modelling strategies

The analysis is performed in two stages, as described in the framework summarized previously. In the first stage, a regression model is applied to individual data for each of the 20 cities included in the analysis in order to derive study-specific estimates of the non-linear exposure–response association, as described in detail in Section 3. The estimated study-specific relationship is entirely defined by the parameters of the function chosen for representing the association, namely, a set of regression coefficients. These coefficients are then used as outcomes for the multivariate meta-analytical model in the second stage, following the methodology described in Section 4. As in conventional univariate meta-analysis, the main aim is to derive a set of regression coefficients defining an average exposure–response association across the studies. In a second step, we aim to test and quantify the amount of heterogeneity and to assess the extent to which this heterogeneity is related to study-specific characteristics, such as those listed inTable I, through multivariate meta-regression.

3.1. First-stage model used in our example

In our example, the first-stage analysis is carried out using a generalized linear model with quasi-Poisson family for overdispersed data, following a standard analytical approach for time-series environmental data [22,15]. The general algebraical definition is given by the following:

(1)

whereg is a log link function of the expectationμ_ti ≡E(Y_ti), withY_ti as the series of 5114 non-accidental daily mortality counts. The exposurez_ti corresponds to outdoor mean daily temperature over lag 0–5, computed as the moving average on dayt and on the previous 5 days, in order to account for lagged effects. Here, the exposure–response functions is chosen as a quadratic B-spline, entirely defined byk − 2 internal knots and 2 boundary knots, wherek corresponds to the dimension of the spline basis and the number of parameters. This is equal to the degrees of freedom (df) spent in the first-stage model to estimate the relationship. Number and location of knots are chosen by the selection method detailed in Section 3.3.

The model also includes two potential confounding variablesc_tji, included a priori, represented by elapsed time and day of the week, modelled through the functionsh_ji and related parametersγ_ji. Elapsed time is used to control for seasonal and long-term trends and specified through a natural cubic B-spline with 7df per year. Day of the week is modelled with six indicator variables through a dummy parameterization.

3.2. Definition of the exposure–response function

In Equation(1), while defining the first-stage regression model, we intentionally omitted the indexi for the functions to indicate that the same function needs to be applied in all the study-specific models. This requirement assures that the estimated coefficients retain the same interpretation in all them studies, so that their meta-analysis can provide meaningful and interpretable results. In our example, this condition is met by placing the internal and boundary knots at the same temperatures in all them studies. Functions for non-linear associations are usually centred on a specific value, so that the results may be interpreted as the effect of the exposure versus a reference [23]. The choice of the centering value depends on interpretational issues and does not affect the fit of the model in(1).

3.3. Control for confounders and model selection

The specification of the study-specific models in(1) can follow two different approaches regarding the control of confounders. In this example, the confounder model, although involving several terms, is conceptually simple, and, consistently with previous time-series environmental studies on the same dataset [22], an identical set of pre-specified predictors is included in all studies. An alternative approach is based on study-specific confounder models. The choice depends, among other issues, on the research setting and the quality of the data. Whatever the strategy adopted to control for confounding in the first stage, the aim is to obtain a valid estimate of the study-specific parametersθ_i, and this choice does not affect the development of the second-stage multivariate meta-analysis.

Different selection methods are available for the definition of the functions. In this illustration, the choice reduces to the identification of the optimal number and location of knots for the quadratic B-spline. In the context of time-series environmental studies, several criteria have been proposed, although consensus as to an optimal approach has not been reached yet. Here we rely on the Q-AIC, a modification of the Akaike information criterion (AIC) for quasi-likelihood models [24]. We first define a list of model candidates, specified by an increasing number of internal knots placed at temperatures corresponding to optimal average percentiles across cities [25], Section 2.4. Boundary knots are placed at temperatures corresponding to the average minimum and maximum temperature. The selected model has the minimum value of the sum of the Q-AIC in all the 20 studies. As with other issues discussed previously concerning the first-stage model, different selection criteria may be preferred in different settings.

3.4. Absolute and relative scales

The meta-analysis of exposure–response relationships requires that the exposure is measured on the same scale in all the studies. Also, in the presence of different exposure ranges, common knots placement may leave some parameters inestimable in some studies. The meta-analysis of their parameters therefore presents additional complications.

The modelling strategy described previously defines the exposure on anabsolute scale (°C). The functions, common to all studies, is coherently defined with knots placed at the same absolute exposure values, so that the estimated coefficients have the same interpretation across studies. Their synthesis through meta-analysis is therefore meaningful. In our example, selected cities have largely overlapping temperature ranges, as showed inTable I, and the specification of common interior and boundary knots does not raise important issues. However, in analyses including studies with more varied temperature ranges, such as of all the 108 US cities, the definition of common knots would be problematic or impractical.

An alternative is to adopt arelative scale, standardizing the study-specific distributions by transforming temperatures to the related study-specific percentiles. This method allows the comparison of studies showing even non-overlapping ranges. The standardization is carried out by selecting the internal knots at the same percentiles in all the studies, which indeed correspond to different absolute temperatures. The interpretation of the estimated coefficients changes, and the comparison must be made on the relative scale of percentiles, as illustrated while commenting on the results in Section 5.4.

4.1. The model

The framework we use is nested within that of the multivariate normal linear mixed model and so follows well-developed lines [27]. It will be presented in the specific context of multi-parameter associations. As anticipated previously, we assume that a first-stage model has been fitted to the data from each of thei = 1, …,m studies, obtainingk-length vectors of regression coefficients, and associatedk ×k estimated (co)variance matricesS_i. These regression coefficients are then used as outcomes for the second stage and are termed from here on as outcome parameters in order to distinguish them from the coefficients of the second-stage meta-analytic model.

Following Jackson and colleagues [11], a model for random-effect multivariate meta-analysis can be written as follows:

(2)

withS_i +Ψ =Σ_i. The marginal model in (2) has independent within-study and between-study components. In the within-study component, the estimated is assumed to be sampled with error from N_k(θ_i,S_i), a multivariate normal distribution of dimensionk, whereθ_i is the vector of true unknown outcome parameters for studyi. In the between-study component,θ_i is assumed sampled from N_k(θ,Ψ), whereΨ is the unknown between-study (co)variance matrix. Hereθ can be interpreted as the population-average outcome parameters, namely the coefficients of the functions defining the average exposure–response association. Model (2) can be extended to multivariate meta-regression, where thek outcomes are modelled in terms of study-level variables defining a set ofp meta-predictorsx_i = [x_1i,x_2i, …,x_pi]^T associated with theith study, where usuallyx_1i = 1 specifies intercept terms. Algebraically,

(3)

Thek ×kp block-diagonal matrixX_i, of rankkp, is derived by the Kronecker product between an identity matrixI_(k) of dimensionk and the vectorx_i, following

(4)

Thekp-dimensional coefficient vectorβ defines the association of thek outcomes withp meta-predictors. They commonly includek intercepts with a similar interpretation ofθ in (2) but related to the population average of studies characterized by a zero value of meta-predictors. The otherk(p − 1) coefficients inβ express how the outcome parameters change in respect to the meta-predictors. The problem can also be re-expressed in the form of a conventional linear mixed model, defining random effectsu_i ∼ N_k(0,Ψ), which represent study-specific deviations from the average. The model in (3) is then written as follows:

(5)

The matrixΨ is completely defined by a set of parametersξ, dependent on the chosen structure and parameterization. If no a priori structure is assumed,k(k + 1) ∕ 2 terms are needed. Optionally, under the assumption that each outcome parameter is explained only in terms of a subset of thep variables, the related columns ofX and entries ofβ can be excluded, defining different linear predictors. Fixed-effect meta-analytic models presuppose that no heterogeneity exists in the distribution of the outcome parameters and that the random variability is explained only by sampling error, assumingΣ_i ≡S_i. As for the univariate case, estimation procedures treatS as known. From here on, the model in (2) for multivariate meta-analysis will be considered a special case of (3)–(5) whereX ≡I andβ ≡θ. The unknown parameters are thereforeβ and, for random-effect meta-analytic models,ξ.

4.2. Estimation

Different estimation methods have been proposed for random-effect multivariate meta-analysis: likelihood-based methods9,28, estimating equations [29], variants of iterative generalized least squares [8,30], Bayesian approaches [31] and multivariate extensions of the method of moments [32]. Here we will concentrate on maximum likelihood (ML) and restricted maximum likelihood (REML), following an extensive literature within the framework of linear mixed models [27,33,34]. These methods are implemented in the R package mvmeta and applied to perform the analysis in Section 5.

The marginal log-likelihood functionℓ(β,ξ) for model (5) may be written as follows [27]:

(6)

withΣ_i(ξ) written here asΣ_i for ease of notation andn as the total number of observations (usually equal tokm where there are no missing values). Assuming thatξ, and thereforeΨ andΣ, are known, the ML estimates forβ and its (co)variance matrix conditional onξ are obtained by maximizing (6). In this case, closed-form equations are given by generalized least squares estimators:

(7)

WhenΨ is not known, the joint likelihood function in (6) needs to be maximized with respect to bothβ andξ, and iterative methods are required. However, the ML estimator of the (co)variance parametersξ is usually biased downward, as it does not account for the loss of degrees of freedom from the estimation ofβ. An alternative estimator can be obtained from the maximization of the log-likelihood function on the basis of a set ofn −q linearly independent error contrasts, withq as the number of fixed effects coefficients inβ. This restricted log-likelihood (REML) functionℓ_R(ξ), not dependent onβ, may be conveniently expressed as follows [27,33]:

(8)

where is defined in (7) andΣ_i identifyingΣ_i(ξ) as above.

The ML estimates ofβ in fixed-effect meta-analysis are simply obtained by (7), given that, as discussed in Section 4.1,Σ_i equalsS_i and is therefore completely known. The ML and REML estimates in random-effect models can be instead obtained through Newton-type iterative algorithms. For computational purposes, the objective functions in (6) and (8) are both expressed with respect toξ only, and maximization ofℓ(ξ) andℓ_R(ξ) can be achieved by plugging in at each iteration the conditional estimate of in (7) using the current estimate ofξ, until convergence. Additional information on the estimation algorithms used here are provided in the Supplementary Web Appendix.

4.3. Hypothesis testing and model comparison

We can separate inferences about the parameters in models (2) and (3)–(5) into those about the fixed effectsβ, which will typically be of prime interest, and the between-study (co)variance matrixΨ. Inferential procedures, again, follow the theory of linear mixed models [[27], Chap. 6].

Regarding fixed effects, under the marginal model and assumingξ as known, follows a multivariate normal distribution with mean and (co)variance matrix given in (7). By replacingξ by its ML or REML estimator, the corresponding entries of the resulting and may be used for obtaining significance tests or confidence intervals for single coefficients. However, in this specific context, the inference is focused on the set of coefficients defining the association of a specific study-level variable with the joint distribution of the outcome parameters. Multivariate extension of Wald or likelihood ratio (LR) tests are easily derived for this purpose. As usual for REML models, general likelihood theory does hold for LR tests only when comparing models with identical fixed effects structures [33]. Suitable adjustments for the underestimation of (co)variance error matrices due to the uncertainty in the estimation ofΨ are available [35,36] although not yet been implemented in mvmeta.

For random effects, the focus is on comparing models involving different choices about the structure of the between-study (co)variance matrix. In this setting, an interesting hypothesis to test isΨ =0, namely that no heterogeneity between studies exists beyond that explained by sampling variability. Similarly, a LR test between nested models may be performed, which is appropriate in REML models as well, given the identical fixed effects structures. Note, however, that for alternative hypotheses that constrain (co)variance matrices to be positive-definite (Supplementary Web Appendix), the null value lies on a boundary of the parameter space. Under these conditions, the conventional null asymptotic distribution does not hold, and some adjustment has been proposed [37]. A test for the same null hypothesis and distribution has also been developed as the multivariate extension of the CochranQ test for (residual) heterogeneity [7,29]. The test is based on the statistic:

(9)

where are estimated by the correspondent fixed-effect model. An extension of this heterogeneity test for a subset ofβ has also been proposed [29].

In addition, in this meta-analytical setting, the quantification of the heterogeneity among studies, or the residual amount beyond that explained by specific covariates, is also of interest. Indices of heterogeneity analogous to the univariate case may be easily derived, such as theH² andI² statistics [38]. These measures are interpreted as the relative excess in heterogeneity (above that explained by sampling error) and the proportion of total variation attributable to heterogeneity, respectively. Although recently criticized for being dependent on precision of the estimates from the first-stage model [39], these statistics provide simple summaries on the extent of heterogeneity. These measures may be produced for the joint multivariate distribution of thek parameters and interpreted as the amount of heterogeneity in the exposure–response relationships. In this case, from (9), we derive the following:

(10)

The same statistics may also be defined for single-outcome parameters [40] by simply computing the related contribution to the quantityQ anddf =n −q while applying (10). However, as previously discussed, the interpretation of these measures for single parameter of a spline function is not very informative.

More broadly, non-nested models may be compared using fit statistics, in particular Akaike information criterion,, and Bayesian information criterion,, where is the maximum log-likelihood. These statistics may also be used with REML models, with the additional requirement that fixed effects structure be held constant.

4.4. Prediction

In the context of multi-parameter associations, the general tests and fit criteria described previously, although important, are usually insufficient for interpretation. Coefficients inβ refer to single-outcome parameters that are rarely interpretable on their own, and the tests only offer a statistical belief onwhether the multivariate distribution of outcome parameters depends on study-level covariates. However, these procedures fail to inform onhow the latter modifies the former.

In the current setting, prediction represents an important tool to extend the inference from multivariate meta-regression models, offering a method to link specific values of study-level variables with outcome parameters expectations. Given a set of meta-predictor valuesx₀, the model-predicted mean and (co)variance matrix are obtained by the following:

(11)

withX₀ computed fromx₀ following (4). The equations in (11) may be used to predict specific parameters, which are used to define the average exposure–response curve predicted for the set of meta-predictor valuesx₀. The same equations may be used to predict the association in a new study characterized by a specific set of study-level variables by simply increasing the uncertainty in the estimates by adding to in (11) [41].

In addition, the assumptions outlined in Section 4.1 regarding the random-effect multivariate distribution may be exploited to extend the inference regarding study-specific outcome parametersθ_i estimated in the first-stage model, computing the (asymptotic) best linear unbiased prediction (BLUP) [[27], Section 7.4]. The predicted and associated asymptotic (co)variance matrix are as follows:

(12)

for. The BLUP equations in (12) represent the sum of two components: the predicted averaged outcome parameters in (11) and study-specific deviations, predicted as random effectsu_i in (5). Exposure–response associations predicted with BLUP in (12) represent a trade-off between city-specific and average estimates, with weights inversely proportional to the two componentsΨ andS_i of the total variabilityΣ_i, respectively. It is noteworthy, in this multivariate setting, that the BLUP estimates of missing parameters from the first stage exploit the information about the other study-specific parameters obtained through the between-study (co)variance matrix and may be therefore different from predicted values from (11).

5.1. Pooled exposure–response relationship

The main result of multivariate meta-analysis, applying the model in (2), consists of an estimate of the outcome parametersθ, representing the population-average coefficients of the quadratic B-spline. These are used to compute the average exposure–response relationship depicted inFigure 2, together with confidence intervals obtained through, over the average exposure range. As expected, the plot indicates an increase in risk for non-accidental mortality for both high and low temperatures, with a point of minimum mortality at about 20.6°C. The risk increases approximately linearly in the cold tail but shows a steep non-linear increase at high temperatures. The model estimates a relative risk (RR) of 1.085 (95%CI: 1.069–1.102) and 1.150 (95%CI: 1.076–1.230) for temperatures of − 10°C and 29°C versus 20°C, respectively. The multivariate CochranQ test andI² obtained from (9)–(10) are reported inTable II. These statistics reveal a significant heterogeneity across studies, where 70.2% of the variability in the multivariate outcome, and therefore in the derived exposure–response association, is attributed to true between-study differences.

5.2. Best linear unbiased prediction

As discussed in Section 4.4, the assumptions about the between-city variability, namely the distribution of the random effects in (5), can be used to extend the inference regarding study-specific estimates. Variability around the average is due to both heterogeneity between cities and uncertainty in the first-stage model. Applying (12), we derived the BLUP estimates of the unknown city-specific parametersθ_i of the quadratic spline, representing a shrunk version of the estimated towards the population average.Figure 3 illustrates the associated exposure–response curves: the left panel shows the original estimates from the first-stage models in the 20 studies, whereas the related BLUP relationships are displayed in the right panel. The shrinkage effect is stronger in the cold extreme, possibly because of the higher uncertainty in the within-study estimates.

Given the large difference in population size shown inTable I and in the related within-study precision, the shrinkage effect is expected to vary considerably across studies.Figure 4 shows the first-stage and BLUP temperature–mortality relationships in two studies, together with the population average as depicted inFigure 2. As expected, the BLUP estimate is closer to the original curve in Chicago, the largest city characterized by a high number of daily deaths (above 5 million inhabitants, 115 deaths on average per day). By contrast, the study in the smallest city of Muskegon (170 000 inhabitants, five deaths per day) produces very imprecise estimates (confidence intervals not shown), in particular for cold, which are heavily shrunk toward the population average.

5.3. Multivariate meta-regression

Part of the large heterogeneity in temperature–mortality curves between studies, displayed inFigure 3, may be explained by different patterns of study-level factors that modify the association. Multivariate meta-regression models in (3)–(5), including meta-predictors, are applied to extend the analysis and partly characterize such differences. From a fitted model, it is possible to compute the outcome parameters predicted for a set of meta-predictor valuesx₀ using (11). As discussed in Section 4.4, these parameters are interpreted as the average coefficients of the functions for the subset of studies from the hypothetical population that are characterized byx =x₀. A test on the significance of the multivariate association between the outcome parameters and each study-level variable is carried out by applying a Wald test on the related subset of coefficientsβ. Alternatively, an LR test or information criteria such as AIC and BIC can be used to compare the models with and without specific meta-predictors.

As an illustration, univariable multivariate meta-regression models, each including a single meta-predictor, are specified for four study-level variables, namely latitude, percentage of population living in urban settings, percentage of population with high education (high school or higher) and percentage of population living in poverty. The exposure–response associations are predicted for the values of the approximate 25th and 75th percentiles of the study-level variables (reported inTable I) and shown inFigure 5. The top-left panel illustrates the two predicted curves for latitude. The effect of hot temperature, on average, is higher in studies for northern cities, with a steeper raise in risk: the RR at 29°C versus 20°C increases from 1.101 (95%CI: 1.049–1.156) for the 25th percentile of latitude to 1.297 (95%CI: 1.205–1.396) for the 75th percentile, although there is also a suggestion of a lower point of minimum mortality when increasing latitude. No substantial effect modification is suggested for cold temperatures. The importance of the overall association with the study-level variable can be assessed by the comparison of the AIC and BIC criteria, reported inTable II, with the corresponding model with no meta-predictor inFigure 2. The lower AIC, from −322.8 to −326.5, suggests a better fit of the model with the predictor, whereas the BIC, highly penalized by the number of observations in the model, shows a preference for the more parsimonious model. The Wald and LR tests, also reported inTable II, are both specified with 6df, corresponding to the six additional coefficients used to model the linear effects of latitude on each outcome parameters. The tests are both significant at the standard 5% level, although thep-values differ somewhat. Latitude seems to explain a substantial part of heterogeneity between studies, with an overallI² of 41.9% compared with the 70.2% of the model with no predictor. However, the test for the residual amount of heterogeneity is still significant.

The other three panels inFigure 5 show the association with the other study-level variables. The percentage of population living in urban settings seems to increase the effect of both cold and hot temperatures, although fit statistics and tests inTable II suggest that this effect modification is not significant. Also, the decrease in heterogeneity is negligible. There is no evidence that percentage of population with higher education or living in poverty modify the temperature–mortality association, as shown by the two identical curves predicted for their 25th and 75th percentiles in the bottom panels ofFigure 5 and related tests and fit statistics inTable II.

The extension to multivariable multivariate meta-regression is straightforward: tests and statistics are defined in exactly the same way, and predicted effects as those showed inFigure 5, controlled for the effect of other study-level variables, may be computed in a similar way. The interpretation, as with standard multivariable models, refers to the effect of a meta-predictor for constant values of the other ones. Standard model selection procedures for regression models may be applied to choose the set of meta-predictor to be included. An example is included in the code provided in the Supplementary Web Appendix.

5.4. Adopting a relative scale

As shown inTable I and inFigure 3, the exposure ranges are slightly different in each study. In this specific example, range variation was too small to cause problems, but it is likely to do so when exposure ranges differ more substantially. An alternative is to adopt a relative scale based on percentiles, as described in Section 3.4. As an illustration, we repeat the analysis placing the knots at the 5th, 35th, 65th and 95th percentiles of study-specific temperature distribution and centering on the 75th percentile. The predicted exposure–response relationship is displayed inFigure 6. To aid interpretation in this standardized scale, thex-axis is scaled so that percentiles match those of the average temperature distribution across all the cities.

Not surprisingly, given the similar temperature ranges, the shape of the pooled association shown in the left panel ofFigure 6 is very close to that predicted in absolute scale inFigure 2. However, the results from the two approaches are not easily comparable because of the different interpretation of the results. The exposure–response relationship predicted in an absolute scale informs about the average RR between specific temperatures, for example 30°C versus 20°C. By contrast, the curve predicted in relative scale measures the RR of relative temperatures, for example 95th versus 75th percentiles, based on study-specific temperature distributions. The right panel ofFigure 6 illustrates the results from the meta-regression model including latitude. Although the risk of hot temperatures is still higher in northern cities, the difference is smaller if compared with the related plot predicted in absolute scale inFigure 5, top-left panel, and no longer significant, as shown inTable II. Also, the point of minimum mortality is similar in cities at different latitudes, approximately at the 76th percentile. This may be explained by a partial adaptation of populations living in different cities to their own climate.

5.5. Additional analysis

The selection procedure regarding the functions described in Section 3.3 indicates a preference for a spline with 4 knots. As a sensitivity analysis, the results in Section 5.1 are compared with those obtained using a spline with 3 knots, placed at − 1.8°C, 11.3°C and 23.3°C, corresponding to the 10th, 50th and 90th average distribution percentiles across studies. The exposure–response curve is very similar to that showed inFigure 2 (result not shown), suggesting that the choice of number and position of knots is not critical in this case.

All the models presented previously have been estimated through ML. The choice of this estimator allows the comparison of models with different fixed effects structures through LR test and information criteria, as described in Sections 4.2–4.3. The estimation can be potentially improved through the use of REML, which accounts for some additional uncertainty in the estimation of the between-study (co)variance matrixΨ not accounted by ML estimates. However, even with this relatively limited number of studies, standard errors of coefficients for REML models only increase in the range of 1–2%, and the two predicted exposure–response curves and confidence intervals are almost identical (Supplementary Web Appendix).

As previously pointed out, estimation difficulties may arise in this complex modelling framework, particularly regarding between-study correlations close to the boundary of their parameter space of − 1 or 1 [11]. The procedure implemented in mvmeta, based on Cholesky decomposition to assure positive-definiteness, seems to perform well under several modelling choices using these data, and estimation problems have not occurred.

In meta-regression models presented in Section 5.3, 120 correlated outcomes from 20 studies are used to estimate 12 coefficients and 21 (co)variance parameters. These sophisticated models require several assumptions that are commonly hard to verify, such as the multivariate normality of the random effectsu_i in (5). The results are therefore expected to be potentially sensitive to model mis-specifications. In particular, we have noticed some discrepancies between the Wald and LR tests in meta-regression models inFigure 5, probably because of different assumptions of the two tests. In a preliminary assessment, the quadratic approximation to the likelihood surface of the Wald test seems to perform poorly in badly conditioned models, for example meta-regression in the presence of outliers. Fitted models may be checked with appropriate regression diagnostics, such as residuals and influence measures already available for univariate meta-analysis [42], which still need to be developed in the multivariate meta-analytical techniques.

• 1.Normand SLT. Meta-analysis: formulating, evaluating, combining, and reporting. Statistics in Medicine. 1999;18(3):321–359. doi: 10.1002/(sici)1097-0258(19990215)18:3<321::aid-sim28>3.0.co;2-p. [DOI] [PubMed] [Google Scholar]
• 2.Sutton AJ, Higgins JP. Recent developments in meta-analysis. Statistics in Medicine. 2008;27(5):625–650. doi: 10.1002/sim.2934. [DOI] [PubMed] [Google Scholar]
• 3.Stukel TA, Demidenko E, Dykes J, Karagas MR. Two-stage methods for the analysis of pooled data. Statistics in Medicine. 2001;20(14):2115–2130. doi: 10.1002/sim.852. [DOI] [PubMed] [Google Scholar]
• 4.Simmonds MC, Higgins J, Stewart LA, Tierney JF, Clarke MJ, Thompson SG. Meta-analysis of individual patient data from randomized trials: a review of methods used in practice. Clinical Trials. 2005;2(3):209–217. doi: 10.1191/1740774505cn087oa. [DOI] [PubMed] [Google Scholar]
• 5.Thompson S, Kaptoge S, White I, Wood A, Perry P, Danesh J Emerging Risk Factors Collaboration. Statistical methods for the time-to-event analysis of individual participant data from multiple epidemiological studies. International Journal of Epidemiology. 2010;39:1345–1359. doi: 10.1093/ije/dyq063. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 6.Dominici F, Samet JM, Zeger SL. Combining evidence on air pollution and daily mortality from the 20 largest US cities: a hierarchical modelling strategy. Journal of the Royal Statistical Society: Series A. 2000;163(3):263–302. [Google Scholar]
• 7.Berkey CS, Anderson JJ, Hoaglin DC. Multiple-outcome meta-analysis of clinical trials. Statistics in Medicine. 1996;15(5):537–557. doi: 10.1002/(SICI)1097-0258(19960315)15:5<537::AID-SIM176>3.0.CO;2-S. [DOI] [PubMed] [Google Scholar]
• 8.Berkey CS, Hoaglin DC, Antczak-Bouckoms A, Mosteller F, Colditz GA. Meta-analysis of multiple outcomes by regression with random effects. Statistics in Medicine. 1998;17(22):2537–2550. doi: 10.1002/(sici)1097-0258(19981130)17:22<2537::aid-sim953>3.0.co;2-c. [DOI] [PubMed] [Google Scholar]
• 9.Arends LR, Vokó Z, Stijnen T. Combining multiple outcome measures in a meta-analysis: an application. Statistics in Medicine. 2003;22(8):1335–1353. doi: 10.1002/sim.1370. [DOI] [PubMed] [Google Scholar]
• 10.Gasparrini A, Armstrong B, Kenward MG. Multivariate meta-analysis: a method to summarize non-linear associations. Statistics in Medicine. 2011;30(20):2504–2506. doi: 10.1002/sim.4226. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 11.Jackson D, Riley R, White IR. Multivariate meta-analysis: potential and promise. Statistics in Medicine. 2011;30(20):2481–2498. doi: 10.1002/sim.4172. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 12.R Development Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2012. URLhttp://www.R-project.org/, ISBN 3-900051-0-0. [Google Scholar]
• 13.Gasparrini A, Armstrong B. Time series analysis on the health effects of temperature: advancements and limitations. Environmental Research. 2010;110(6):633–638. doi: 10.1016/j.envres.2010.06.005. [DOI] [PubMed] [Google Scholar]
• 14.Basu R, Samet JM. Relation between elevated ambient temperature and mortality: a review of the epidemiologic evidence. Epidemiologic Reviews. 2002;24(2):190–202. doi: 10.1093/epirev/mxf007. [DOI] [PubMed] [Google Scholar]
• 15.Armstrong B. Models for the relationship between ambient temperature and daily mortality. Epidemiology. 2006;17(6):624–631. doi: 10.1097/01.ede.0000239732.50999.8f. [DOI] [PubMed] [Google Scholar]
• 16.McMichael AJ, Wilkinson P, Kovats RS, Pattenden S, Hajat S, Armstrong B, Vajanapoom N, Niciu EM, Mahomed H, Kingkeow C, Kosnik M, O'Neill MS, Romieu I, Ramirez-Aguilar M, Barreto ML, Gouveia N, Nikiforov B. International study of temperature, heat and urban mortality: the ‘ISOTHURM’ project. International Journal of Epidemiology. 2008;37(5):1121. doi: 10.1093/ije/dyn086. [DOI] [PubMed] [Google Scholar]
• 17.Baccini M, Biggeri A, Accetta G, Kosatsky T, Katsouyanni K, Analitis A, Anderson HR, Bisanti L, D'Ippoliti D, Danova J, Forsberg B, Medina S, Paldy A, Rabczenko D, Schindler C, Michelozzi P. Heat effects on mortality in 15 European cities. Epidemiology. 2008;19(5):711–719. doi: 10.1097/EDE.0b013e318176bfcd. [DOI] [PubMed] [Google Scholar]
• 18.Anderson BG, Bell ML. Weather-related mortality: how heat, cold, and heat waves affect mortality in the United States. Epidemiology. 2009;20(2):205–213. doi: 10.1097/EDE.0b013e318190ee08. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 19.Ruppert D, Wand MP, Carroll RJ. Semiparametric Regression. New York: Cambridge University Press; 2003. [Google Scholar]
• 20.Royston P, Altman DG. Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling. Journal of the Royal Statistical Society: Series C. 1994;43(3):429–467. [Google Scholar]
• 21.Muggeo VM. Estimating regression models with unknown break-points. Statistics in Medicine. 2003;22(19):3055–3071. doi: 10.1002/sim.1545. [DOI] [PubMed] [Google Scholar]
• 22.Dominici F, Sheppard L, Clyde M. Health effects of air pollution: a statistical review. International Statistical Review. 2003;71(2):243–276. [Google Scholar]
• 23.Cao J, Valois MF, Goldberg MS. An S-Plus function to calculate relative risks and adjusted means for regression models using natural splines. Computer Methods and Programs in Biomedicine. 2006;84(1):58–62. doi: 10.1016/j.cmpb.2006.08.004. [DOI] [PubMed] [Google Scholar]
• 24.Peng RD, Dominici F, Louis TA. Model choice in time series studies of air pollution and mortality. Journal of the Royal Statistical Society: Series A. 2006;169(2):179–203. [Google Scholar]
• 25.Harrell FE. Regression Modeling Strategies: With Applications to Linear Models, Logistic Regression, and Survival Analysis. New York: Springer Verlag; 2001. [Google Scholar]
• 26.Higgins J, Thompson SG, Spiegelhalter DJ. A re-evaluation of random-effects meta-analysis. Journal of the Royal Statistical Society: Series A. 2009;172(1):137–159. doi: 10.1111/j.1467-985X.2008.00552.x. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 27.Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. New York: Springer Verlag; 2000. [Google Scholar]
• 28.White IR. Multivariate random-effects meta-analysis. Stata Journal. 2009;9(1):40–56. [Google Scholar]
• 29.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;139(7):1919–1933. [Google Scholar]
• 30.Dear KBG. Iterative generalized least squares for meta-analysis of survival data at multiple times. Biometrics. 1994;50:989–1002. [PubMed] [Google Scholar]
• 31.Nam IS, Mengersen K, Garthwaite P. Multivariate meta-analysis. Statistics in Medicine. 2003;22(14):2309–2333. doi: 10.1002/sim.1410. [DOI] [PubMed] [Google Scholar]
• 32.Jackson D, White IR, Thompson SG. Extending DerSimonian and Laird's methodology to perform multivariate random effects meta-analyses. Statistics in Medicine. 2010;29(12):1282–1297. doi: 10.1002/sim.3602. [DOI] [PubMed] [Google Scholar]
• 33.Harville DA. Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association. 1977;72(358):320–338. [Google Scholar]
• 34.Pinheiro JC, Bates DM. Mixed-Effects Models in S and S-PLUS. New York: Springer Verlag; 2000. [Google Scholar]
• 35.Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997;53(3):983–997. [PubMed] [Google Scholar]
• 36.Kenward MG, Roger JH. An improved approximation to the precision of fixed effects from restricted maximum likelihood. Computational Statistics and Data Analysis. 2009;53(7):2583–2595. [Google Scholar]
• 37.Stram DO, Lee JW. Variance components testing in the longitudinal mixed effects model. Biometrics. 1994;50(4):1171–1177. [PubMed] [Google Scholar]
• 38.Higgins JPT, Thompson SG. Quantifying heterogeneity in a meta-analysis. Statistics in Medicine. 2002;21(11):1539–1558. doi: 10.1002/sim.1186. [DOI] [PubMed] [Google Scholar]
• 39.Rucker G, Schwarzer G, Carpenter JR, Schumacher M. Undue reliance on I2 in assessing heterogeneity may mislead. BMC Medical Research Methodology. 2008;8(1):79. doi: 10.1186/1471-2288-8-79. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 40.White IR. Multivariate random-effects meta-regression: updates to mvmeta. Stata Journal. 2011;11(2):255–270. [Google Scholar]
• 41.Riley R, Higgins J, Deeks JJ. The interpretation of random-effects meta-analysis. British Medical Journal. 2011;342:d549. doi: 10.1136/bmj.d549. [DOI] [PubMed] [Google Scholar]
• 42.Viechtbauer W, Cheung MWL. Outlier and influence diagnostics for meta-analysis. Research Synthesis Methods. 2010;1(2):112–125. doi: 10.1002/jrsm.11. [DOI] [PubMed] [Google Scholar]
• 43.Ades AE. A chain of evidence with mixed comparisons: models for multi-parameter synthesis and consistency of evidence. Statistics in Medicine. 2003;22(19):2995–3016. doi: 10.1002/sim.1566. [DOI] [PubMed] [Google Scholar]
• 44.Ades AE, Sutton AJ. Multiparameter evidence synthesis in epidemiology and medical decision-making: current approaches. Journal of the Royal Statistical Society: Series A. 2006;169(1):5–35. [Google Scholar]
• 45.Hedges LV, Tipton E, Johnson MC. Robust variance estimation in meta-regression with dependent effect size estimates. Research Synthesis Methods. 2010;1(1):39–65. doi: 10.1002/jrsm.5. [DOI] [PubMed] [Google Scholar]
• 46.Berlin JA, Santanna J, Schmid CH, Szczech LA, Feldman HI. Individual patient- versus group-level data meta-regressions for the investigation of treatment effect modifiers: ecological bias rears its ugly head. Statistics in Medicine. 2002;21(3):371–387. doi: 10.1002/sim.1023. [DOI] [PubMed] [Google Scholar]
• 47.Lambert PC, Sutton AJ, Abrams KR, Jones DR. A comparison of summary patient-level covariates in meta-regression with individual patient data meta-analysis. Journal of Clinical Epidemiology. 2002;55(1):86–94. doi: 10.1016/s0895-4356(01)00414-0. [DOI] [PubMed] [Google Scholar]
• 48.Olkin I, Sampson A. Comparison of meta-analysis versus analysis of variance of individual patient data. Biometrics. 1998;54(1):317–322. [PubMed] [Google Scholar]
• 49.Mathew T, Nordstrom K. Comparison of one-step and two-step meta-analysis models using individual patient data. Biometrical Journal. 2010;52(2):271–287. doi: 10.1002/bimj.200900143. [DOI] [PubMed] [Google Scholar]
• 50.Riley RD, Thompson JR, Abrams KR. An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown. Biostatistics. 2008;9(1):172. doi: 10.1093/biostatistics/kxm023. [DOI] [PubMed] [Google Scholar]
• 51.Riley RD. Multivariate meta-analysis: the effect of ignoring within-study correlation. Journal of the Royal Statistical Society: Series A. 2009;172(4):789–811. [Google Scholar]
• 52.Sauerbrei W, Royston P. A new strategy for meta-analysis of continuous covariates in observational studies. Statistics in Medicine. 2011;30(28):3341–3360. doi: 10.1002/sim.4333. [DOI] [PubMed] [Google Scholar]
• 53.Schwartz J, Zanobetti A. Using meta-smoothing to estimate dose–response trends across multiple studies, with application to air pollution and daily death. Epidemiology. 2000;11(6):666–672. doi: 10.1097/00001648-200011000-00009. [DOI] [PubMed] [Google Scholar]
• 54.Dominici F, Daniels MJ, Zeger SL, Samet JM. Air pollution and mortality: estimating regional and national dose–response relationships. Journal of the American Statistical Association. 2002;97:100–111. [Google Scholar]
• 55.Samoli E, Analitis A, Touloumi G, Schwartz J, Anderson HR, Sunyer J, Bisanti L, Zmirou D, Vonk JM, Pekkanen J, Goodman P, Paldy A, Schindler C, Katsouyanni K. Estimating the exposure–response relationships between particulate matter and mortality within the APHEA multicity project. Environmental Health Perspectives. 2005;113(1):88–97. doi: 10.1289/ehp.7387. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 56.Greenland S, Longnecker MP. Methods for trend estimation from summarized dose–response data, with applications to meta-analysis. American Journal of Epidemiology. 1992;135(11):1301–1309. doi: 10.1093/oxfordjournals.aje.a116237. [DOI] [PubMed] [Google Scholar]
• 57.Berlin JA, Longnecker MP, Greenland S. Meta-analysis of epidemiologic dose–response data. Epidemiology. 1993;4(3):218–228. doi: 10.1097/00001648-199305000-00005. [DOI] [PubMed] [Google Scholar]
• 58.Bagnardi V, Zambon A, Quatto P, Corrao G. Flexible meta-regression functions for modeling aggregate dose–response data, with an application to alcohol and mortality. American Journal of Epidemiology. 2004;159(11):1077–1086. doi: 10.1093/aje/kwh142. [DOI] [PubMed] [Google Scholar]
• 59.Hartemink N, Boshuizen HC, Nagelkerke NJ, Jacobs MA, van Houwelingen HC. Combining risk estimates from observational studies with different exposure cutpoints: a meta-analysis on body mass index and diabetes type 2. American Journal of Epidemiology. 2006;163(11):1042–1052. doi: 10.1093/aje/kwj141. [DOI] [PubMed] [Google Scholar]
• 60.Vlaanderen J, Portengen L, Rothman N, Lan Q, Kromhout H, Vermeulen R. Flexible meta-regression to assess the shape of the benzene-leukemia exposure–response curve. Environmental Health Perspectives. 2010;118(4):526–532. doi: 10.1289/ehp.0901127. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 61.Shi JQ, Copas JB. Meta-analysis for trend estimation. Statistics in Medicine. 2004;23(1):3–19. doi: 10.1002/sim.1595. [DOI] [PubMed] [Google Scholar]
• 62.Orsini N. Meta-analysis for linear and nonlinear dose–response relations: examples, an evaluation of approximations, and software. American Journal of Epidemiology. 2012;175(1):66–73. doi: 10.1093/aje/kwr265. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 63.Liu Q, Cook NR, Bergstrom A, Hsieh CC. A two-stage hierarchical regression model for meta-analysis of epidemiologic nonlinear dose–response data. Computational Statistics and Data Analysis. 2009;53(12):4157–4167. [Google Scholar]
• 64.Rota M, Bellocco R, Scotti L, Tramacere I, Jenab M, Corrao G, La Vecchia C, Boffetta P, Bagnardi V. Random-effects meta-regression models for studying nonlinear dose–response relationship, with an application to alcohol and esophageal squamous cell carcinoma. Statistics in Medicine. 2010;29:2679–2687. doi: 10.1002/sim.4041. [DOI] [PubMed] [Google Scholar]
• 65.Arends LR, Hunink MG, Stijnen T. Meta-analysis of summary survival curve data. Statistics in Medicine. 2008;27(22):4381–4396. doi: 10.1002/sim.3311. [DOI] [PubMed] [Google Scholar]
• 66.Ishak KJ, Platt RW, Joseph L, Hanley JA, Caro JJ. Meta-analysis of longitudinal studies. Clinical Trials. 2007;4(5):525–539. doi: 10.1177/1740774507083567. [DOI] [PubMed] [Google Scholar]
• 67.Arends LR, Hamza TH, Van Houwelingen JC, Heijenbrok-Kal MH, Hunink MGM, Stijnen T. Bivariate random effects meta-analysis of ROC curves. Medical Decision Making. 2008;28(5):621–638. doi: 10.1177/0272989X08319957. [DOI] [PubMed] [Google Scholar]
• 68.Gasparrini A, Armstrong B. The impact of heat waves on mortality. Epidemiology. 2011;22(1):68–73. doi: 10.1097/EDE.0b013e3181fdcd99. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 69.Zanobetti A, Schwartz J, Samoli E, Gryparis A, Touloumi G, Atkinson R, Le Tertre A, Bobros J, Celko M, Goren A, Forsberg B, Michelozzi P, Rabczenko D, Aranguez Ruiz E, Katsouyanni K. The temporal pattern of mortality responses to air pollution: a multicity assessment of mortality displacement. Epidemiology. 2002;13(1):87–93. doi: 10.1097/00001648-200201000-00014. [DOI] [PubMed] [Google Scholar]
• 70.Zanobetti A, Schwartz J. Mortality displacement in the association of ozone with mortality: an analysis of 48 cities in the United States. American Journal of Respiratory and Critical Care Medicine. 2008;177(2):184–189. doi: 10.1164/rccm.200706-823OC. [DOI] [PubMed] [Google Scholar]
• 71.Samoli E, Zanobetti A, Schwartz J, Atkinson R, Le Tertre A, Schindler C, Perez L, Cadum E, Pekkanen J, Paldy A, Touloumi G, Katsouyanni K. The temporal pattern of mortality responses to ambient ozone in the APHEA project. Journal of Epidemiology and Community Health. 2009;63:960–966. doi: 10.1136/jech.2008.084012. [DOI] [PubMed] [Google Scholar]
• 72.Analitis A, Katsouyanni K, Biggeri A, Baccini M, Forsberg B, Bisanti L, Kirchmayer U, Ballester F, Cadum E, Goodman PG, Hojs A, Sunyer J, Tiittanen P, Michelozzi P. Effects of cold weather on mortality: results from 15 European cities within the PHEWE Project. American Journal of Epidemiology. 2008;168(12):1397. doi: 10.1093/aje/kwn266. [DOI] [PubMed] [Google Scholar]
• 73.Zanobetti A, Schwartz J. Temperature and mortality in nine US cities. Epidemiology. 2008;19(4):563–570. doi: 10.1097/EDE.0b013e31816d652d. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 74.Gasparrini A. Distributed lag linear and non-linear models in R: the package dlnm. Journal of Statistical Software. 2011;43(8):1–20. [PMC free article] [PubMed] [Google Scholar]
• 75.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. doi: 10.1002/sim.1040. [DOI] [PubMed] [Google Scholar]

Supplementary Materials