Portal de Periódicos da CAPES

Predicting personalised absolute treatment effects in individual participant data meta‐analysis: An introduction to splines

2022; Wiley; Volume: 13; Issue: 2 Linguagem: Inglês

10.1002/jrsm.1546

1759-2887

Michail Belias, Maroeska M. Rovers, Jeroen Hoogland, Johannes B. Reitsma, Thomas P. A. Debray, Joanna IntHout,

Health Systems, Economic Evaluations, Quality of Life

One of the main goals of an individual participant data meta-analysis (IPD-MA) of intervention studies is to investigate whether treatment effect differences are present, and how they are associated with patient characteristics.1 Examining treatment heterogeneity due to a continuous covariable (e.g., BMI or age) may be challenging, since there is often no prior knowledge on functional form of the conditional association between the outcome and the continuous variable. Naïve but often-used approaches are to ignore the possibly complex nature of the functional form and assume a linear effect, or to provide a step-wise approximation based on categorisation of the continuous variable. Categorisation leads to loss of information, reduced power, inflation of the type I error rates, and biased results.2-5 Linear modelling may also lead to biased results since the model may be too simplistic for the data. Therefore, there is a clear need for methods that can capture non-linear relations in IPD-MA context. Modelling treatment effect differences whilst accounting for non-linear functional shapes may provide the opportunity to accurately make inferences whether a patient should be treated or not. To account for non-linearities we may estimate the functional shape of the associations and investigate potential treatment effect differences.6 So far, a variety of methods that account for non-linear functional shapes has been proposed.7-16 In this manuscript, we focus on the use of splines since they can capture both non-linear main effects and non-linear treatment-covariable interaction effects without the need to pre-specify their functional form. The class of spline methods is still broad and includes fully parametric, semi-parametric and even non-parametric approaches, and allows for penalisation and clustering. Thus, splines are very flexible and can address a great variety of fitting problems. We focus on four widely used types of splines; restricted cubic splines as described by Harrell,17 natural B-splines,18, 19 smoothing splines20 and P-splines.21 Splines are being used in single studies, both in intervention and prediction studies, and are also available in IPD-MA context. In IPD-MA, models can be estimated in either one or two stages. A generalised additive mixed effects model (GAMM)22 provides a one-stage IPD-MA method that incorporates the flexibility provided by splines. GAMMs fit a generalised additive model using covariables with or without spline transformation, while adjusting for within-study clustering of the participants based on random effects. In two-stage IPD-MA, appropriate spline-based models are fitted per study in the first stage. Subsequently, study specific estimates are extracted and pooled in the second stage. Pooling study specific predictions and their standard errors is referred to as pointwise meta-analysis.14 Pooling of study-specific model coefficients and their variance–covariance matrix is referred to as multivariate meta-analysis.23 While the methodology has progressed, the use of splines in IPD-MA is relatively uncommon. A possible explanation is that the available guidance is limited. Perperoglou et al. provided a review of splines approaches, but limited to single studies.24 White et al.25 compared pointwise meta-analysis and multivariate meta-analysis techniques in presence of non-linear associations, but used fractional polynomials instead of splines. Gasparrini et al.23 described the use of B-splines in combination with multivariate meta-analysis. They mentioned that multivariate meta-analysis may be combined with other approaches to account for non-linearities but do not provide details. Riley et al.26 described both multivariate meta-analysis and one-stage mixed effects modelling. However, most of the examples were limited to either linear associations or a combination of restricted cubic splines based on truncated power series and multivariate meta-analysis. Our goal is to explain and illustrate how to predict a conditional absolute treatment effect, as this measure is most relevant for clinical decision-making.27 We use the four spline approaches in scenarios with multiple studies, using artificial datasets for illustration. We generated the data such that no confounder adjustment was necessary, but if needed the methods can also be implemented after individual-level confounder adjustment. We describe the various spline approaches and their application in IPD-MA using pointwise meta-analysis,14 multivariate meta-analysis,23 and GAMMs,22 and we provide the corresponding R-code. We also describe the results of the aforementioned spline and pooling methods using an empirical individual participant data-set, investigating the effect of antibiotics in children with acute otitis media (AOM).28 In order to illustrate the aforementioned spline approaches in IPD-MA we generated artificial data to mimic a previously reported non-linear association between BMI and mortality.29, 30 We consider the case where the outcome is binary, but note that splines may be used to other types of outcomes such as continuous and time-to-event outcomes. For the control group we generated a (quadratic) J-shaped association showing increased mortality for underweight and overweight participants. For the experimental group we assume a (quartic) levelled J-shaped association, where the association between BMI and mortality is much weaker, especially for those with a < 30 BMI (Figure 1). To illustrate the performance of splines in IPD-MA we generated three distinct IPD-MA scenarios, each consisting of five RCTs comparing two interventions (1:1 randomisation) with 500 participants per study. In the first scenario, which we refer to as the heterogeneous IPD-set with equal BMI ranges, the association between BMI and mortality is different across studies, see Figure 2, but the distribution and ranges of BMI are the same. In the second scenario, which we refer to as the non-heterogeneous IPD-set with different BMI ranges, the parameter values of the association for both the treated and control group are identical across all studies, but the ranges of available BMIs vary across studies (see Figure 3). In the third scenario, to which we refer as the combined IPD-set with different BMI ranges and between study differences in mortality risks, both the ranges of BMI and the association of BMI with the mortality risk vary across studies, see Figure 4. Exact equations are given in the online Appendix sections 3 and 4. "Treatment effect modification", also called "treatment effect measure modification"31, 32 is the phenomenon where the effect of a treatment varies across the levels or strata of a certain variable. We prefer the term "treatment effect measure modification" since effect modification may be present for one measure (e.g., risk difference) but not for another (e.g. odds ratio, risk ratio).31-35 The scale in which the results are presented is therefore a vital first decision. A commonly applied approach to investigate treatment effect measure modification is to model the interaction of a potential effect modifier with the treatment. In case of non-linear associations, a spline-transformed version of the continuous modifier can be used. Therefore, we model the association between the modifier and the outcome by including a spline-transformed version of the modifier, both as main effect and in interaction with the treatment. In case of a binary outcome like mortality, a logit link function can be used in the model. In order to calculate the absolute risk difference between the treatment arms, we back-transform the predicted outcome per treatment arm with the inverse logit function. To calculate the confidence interval of the difference in absolute risk, we use the approach proposed by Newcombe.36 In a setting where the association between an outcome and a continuous variable X is non-linear, one of the options is to use splines. Splines represent a continuous variable as a linear additive combination of (often)-local parts, which each have a simple mathematical form and are known as basis functions. Numerous basis functions have been developed involving various mathematical forms, such as polynomials, radials and Fourier series. In this paper, we focus on basis functions based on piecewise polynomials. As the term piecewise implies, the range of X is divided into intervals, using cut-offs called knots. Within each interval, a d-degree polynomial of X is used to model the association between the outcome Y and X. These polynomials are connected across adjacent intervals. This way, instead of estimating a global non-linear association over the full range of data, we estimate within intervals the linear association between the outcome and transformations of X. Two important choices have to be made, in addition to the degree of the basis functions: (1) the number and the position of the knots, and (2) whether a penalty should be applied. Splines calculated with the use of pre-specified knots and without penalties are often called regression splines. The most commonly used regression splines are restricted or natural splines17 and B-splines.18, 37, 38 Splines where a penalty is applied are called penalised splines. The most commonly used penalised splines are P-splines21 and smoothing splines.20 A short summary of these four types of splines is presented below. Details are presented in the online Appendix, sections 1 (Regression splines) and 2 (Penalised splines). Figure 5 shows how the aforementioned spline methods are associated with each other. Perperoglou et al. provide more details on the pros and cons of the different spline approaches.24 In order to understand the rationale for the use of splines based on piecewise polynomials such as the ones of interest here, we first shortly introduce global and piecewise polynomials. Global transformations of X (i.e., a function on the full range of X) include classical methods such as regular polynomials of X. Polynomials with successively higher degrees are able to capture successively more non-linearity, but at the cost of instability, especially near the boundaries of X. The reason is twofold. First, each non-linear part of a polynomial is more variable near the boundaries (e.g., think of a cubic function), such that small errors in the estimated coefficients have large impact. Second, each coefficient has an effect along the entire range of X, so the errors accumulate. To avoid this, polynomials fitted on different intervals of X, also called piecewise polynomials, may be preferred to global functions. Forcing the piecewise polynomials to be continuous over the knots leads to "piecewise polynomial" splines, in short "splines" (from a craftsman's tool). To accomplish continuity, several approaches have been proposed. One approach is to fit a global polynomial, whilst including terms that model the deviance from this globally defined shape within truncated parts of X, for example a truncated power series.17 Another approach is to locally generate transformations of X, for example B-splines.39 The number and location of the knots may be based on clinical knowledge or on descriptive statistics. For instance, Harrell suggests the use of quantiles.17 Depending on the available sample size and required complexity of the functional shape, we may use a different number of knots. Smoothing is improved when also the derivatives of the spline function are continuous. Typically continuity up to the second derivative offers sufficient smoothness.24 Since splines borrow information from adjacent intervals there still can be erratic behaviour near the (global) boundaries. To avoid this one may restrict the second derivative of the spline function to be zero at the boundaries. These splines are often called natural or restricted. Note that several types of regression splines may be combined with the "natural" property. For instance, Harrell describes a restricted truncated power series spline of third degree,17 Wood describes a restricted cubic spline based on cardinal basis functions,39 while Chambers and Hastie describe a natural B-spline.38 In our manuscript, we discuss restricted splines as described by Harrell and natural B-splines as described by Chambers and Hastie. A restricted spline based on truncated power series is constructed from a linear (first-degree) global polynomial over the full range of X, where the deviations from this global polynomial are modelled with truncated transformations of X. We refer to truncated power series that are constrained with the "natural" or "restricted" property, as restricted cubic splines, following the terminology of Harrell.17 In our single study example, we used restricted cubic spline transformations of X both as main effects and as interactions with the treatment. We placed five knots at values corresponding to 5%, 27.5%, 50%, 72.5%, and 95% quantiles of X. In Figure 6a, we present the predicted mortality risks per treatment arm, conditional on BMI, along with the 95% confidence intervals. Subsequently, we calculated the effect of the treatment conditional on BMI, by calculating the conditional risk for the control minus the conditional risk for the treated, see Figure 7a. To calculate the absolute risk difference's 95% confidence intervals, we followed the proposal of Newcombe, see section 3. Note that in our artificial data the boundaries and distribution of BMI-values for the treated and control group are the same. In practice, this may not be true and knots may be placed at different positions for the main effects and the interaction terms. B-splines are another commonly applied regression spline approach. As opposed to the global nature of restricted cubic splines, the basis functions of B-splines are generated locally, which improves numerical stability.19 We refer to B-splines that are constrained with the "natural" or "restricted" property as natural B-splines, following the terminology of Chambers and Hastie.38 In Figure 6b we show the results of the natural B-splines approach for the simulated single study data. In order for natural B-splines and restricted cubic splines to be comparable in terms of the degrees of freedom, we used third degree B-spline transformations of X both for the main effects and for the interactions with the treatment. We used five equidistant knots (three inner knots at BMI plus two at the boundaries 18.5 and 40). Subsequently, we calculated the effect of the treatment conditional on BMI, see Figure 7b, similar as for the restricted cubic splines. The main advantages of regression splines are their relative simplicity and the fact that they can be represented by a formula. As a consequence, the estimated regression coefficients can be reported and used in further analysis, for example, meta-analysis. Also both restricted splines and natural B-splines are readily implemented in common models such as generalised linear models (GLMs) and have low computational cost. Natural B-splines provide greater local support and numerical stability than restricted splines, since their basis functions are generated locally.24 Both restricted splines and natural B-splines with κ knots require κ + 1 degrees of freedom.24, 38 The main disadvantage of regression splines is that their model fit heavily depends on the number and position of the knots, thus careful modelling is required to avoid under- or overfitting. In some occasions, clinical knowledge on the expected curvature or descriptive statistics may be used to define the knots, but in others it is unclear how many knots should be used and where they should be placed. Commonly used criteria such as Akaike's information criterion (AIC)40 can be used for a databased choice of the number and position of the knots. Alternatively, penalised splines have been proposed to avoid these issues. The two commonly applied penalised splines that we discuss, P-splines and smoothing splines, increase the number of knots to a large set (usually, 10–40 for P-splines) or even to be equal to the number of observations (smoothing splines). This way they circumvent the problem of choosing the number and positions of the knots. Since estimating one parameter for each observation would clearly lead to a perfect fit and thus generate functional shapes with extreme variability, penalised splines introduce in their optimisation functions a penalty term ( J β ) multiplied by a non-negative λ, often called a tuning parameter. As the term, "tuning" implies, changing the value of λ changes the magnitude of the penalisation. Penalised splines circumvent the problem of knot selection, but at a cost. By using a penalty in their optimisation function, they introduce bias in their estimate in order to obtain a more stable solution. Furthermore, in both P-splines and smoothing splines the tuning parameter λ must be specified. Too high or too low values of λ may lead to over- or undersmoothing respectively. Several approaches have been proposed in order to determine the "optimal" λ, such as the AIC,40 "leave one out" generalised cross-validation (GCV)41 or mixed-effects modelling.22 These processes are automated in most of the statistical packages. Briefly, when using the AIC, a series of models fitted with different λ values are compared and the one with the lowest AIC is selected. "Leave one out" GCV is an iterative process, the algorithm goes as follows: (1) one observation is omitted, (2) a model is fitted, (3) using the model, a prediction of the omitted value is generated and (4) the distance between the observed and predicted value is calculated. This procedure is repeated for each observation and for a series of λ values. The λ that minimises the sum of the squared distances, that is, the GCV score, is selected. In Bayesian/mixed effects, modelling approaches the penalty term is estimated in a similar way as a random effects parameter.22 More details are given in the online Appendix section 2. A specific type of penalised splines, P-splines, proposed by Eilers and Marx,21 is a penalised version of B-splines, using a specific penalty term based on the sum of squared p-order differences between the coefficients of two consecutive intervals J β = ∑ Δ p β w 2 . Due to this penalty, the number of basis functions and thereby the flexibility is allowed to be large, but the penalty forces adjacent coefficients to be similar when the data do not support such flexibility. The choice of the number of knots is not as important as with regression splines as long as it is sufficiently large. Note that the degree of the underlying B-splines may be different from the order of the differences. A common combination is that of a third degree B-spline with a second order difference. Using a penalty based on a zero-order difference results in the ridge penalty.42 P-splines are usually based on equidistant knots. While it is possible to use a knot sequence that is not evenly spaced, this would require the introduction of weights,22, 24 and they are rarely used in practice. In our single study example, we used P-spline transformations of X for both its main effect and the interaction with treatment. We used 17 equidistant knots; 15 inner knots plus the boundaries, while the λ parameter was selected through a 'leave one out' GCV process as described above. In Figure 6c, we present the resulting mortality risks per treatment arm conditional on BMI, along with the 95% confidence intervals. Subsequently, the effect of the treatment conditional on BMI, calculated as the difference between the two curves in Figure 6c, is presented in Figure 7c. Smoothing splines are another member of the family of penalised spline methods. Similar to P-splines, the idea is to increase the number of knots, but this time to be equal or approximately equal to the number of observations. O′ Sullivan43 suggested a penalty based on Reinsch's integral of the second derivative of f X , where f X is a cubic spline. In our single study example, we used smoothing spline transformations of X for both the main effect and its interaction with treatment, while the λ parameter was selected through a 'leave one out' GCV process as described above. In Figure 6d, we present the resulting mortality risks per treatment arm conditional on BMI, along with the 95% confidence intervals. The effect of treatment conditional on BMI is presented in Figure 7d. Penalised splines reflect the belief that the predicted regression lines are more likely to be smooth than not. Therefore, their main advantage is that they are more likely to show smoother functional shapes as compared to unpenalised splines. Another advantage is that they circumvent the need to specify the positions and the number of knots, which in most cases are not known beforehand and may need to be estimated. Penalisation also affects the inference, due to the bias-variance trade-off. For instance, the coefficient estimates are subject to a smoothing bias, therefore their interpretation may be problematic. Note that this issue does not necessarily apply to the predicted outcomes. A related issue is that the degrees of freedom have to be modified to account for the penalisation. Wood suggests the use of effective degrees of freedom of a model. Effective degrees of freedom are calculated using the Welch-Satterthwaite approximation formula and can be used to compare models fitted with different types of splines.44 In the previous sections, we focused on estimating non-linear main effects and interactions with the treatment in a single study. As trials are typically not powered to investigate effect modifiers, exploring non-linear effects in a single study may often be problematic or yield very wide confidence intervals. Depending on the underlying curvature, splines need a high amount of data and therefore their use is more feasible in the context of an IPD-MA, where they enable the statistical modelling of complex relationships such as non-linear associations.45 They can be applied in a two-stage or one-stage meta-analysis approach. We apply the methods on three IPD-MA scenarios of five studies each. In the first scenario, the underlying curves are heterogeneous whilst the BMI ranges are the same across studies. In the second scenario, the underlying curves are homogeneous but the BMI ranges are different across studies. In the third scenario, both the underlying curves and the BMI ranges are different across studies. In pointwise meta-analysis, a separate meta-analysis is conducted per distinct value (point) of X, using the outcomes and standard errors as estimated per study. In the first stage of a two-stage pointwise meta-analysis, as proposed by Sauerbrei and Royston in method 3,14 we fit an appropriate model and estimate the predicted outcome per study. If needed, e.g. in case of observational data, the predicted outcomes should be controlled for individual-level confounders by first building a confounder model.14 Note that instead of using fractional polynomials as in Sauerbrei and Royston, we may use any of the spline approaches described in section 4. After fitting the study-specific models, we should decide, for example, by plotting the results, whether it is sensible to pool the predicted outcomes across studies. If so, the second stage consists of distinct meta-analyses, for each value of X, of the study-specific predicted values and standard errors, using either a fixed or a random effects approach.14 Given a continuous variable X the algorithm proceeds as follows: Depending on the outcome we wish to show and depending on the scale on which we wish to make inferences, we may choose to use a link function g and its inverse g − 1 : For each value within the boundaries of X we perform either a fixed or random effects meta-analysis to get the pooled outcome of choice as a function of X along with its pointwise 95% confidence interval. Note that if the available data across the studies vary over different regions of X, pooling of the predicted outcomes may produce unsmooth results, see Figures 8 and 9, especially in the second and third scenario. From our experience, this may be caused by unsmooth estimates of the heterogeneity parameter τ across the range of X. As τ plays an important role in random-effects meta-analysis, it is recommended to visualise the heterogeneity estimates over the range of X. We applied pointwise meta-analysis using all aforementioned spline approaches in all three IPD-MA scenarios. First, we estimated the mortality risk conditional on BMI per study and treatment arm (step 3.1), and we estimated per study the absolute risk difference and the confidence interval (step 3.2) over the full range of X. In the second stage, we pooled both the predicted curves per treatment arm (on logit scale) and their risk differences (on probability scale), using random effects meta-analyses with REML estimators for τ2, except when there were problems with the convergence of the Fisher scoring algorithm in the estimation of τ2. In those cases, we replaced the REML estimator with the DerSimonian-Laird estimator for τ2 which does not require an iterative procedure, and which in our example generally resulted in smoother estimates than REML. In a random effects meta-analysis the size of the estimated τ2 affects the pooling weights of the individual studies. The pooled mortality risks per treatment arm are presented in Figure 8, and the pooled treatment effects conditional on BMI in Figure 9. Detailed figures for the third scenario, including visualisations containing study data, can be found in the online Appendix, section 5.1 for the first stage, and sections 5.2–5.5 for the second stage. Since the underlying curvature (i.e., the amount by which the curve deviates from being a straight line) per study is different across the first and the two remaining scenarios, we followed different modelling strategies for the regression splines to avoid overfitting during the first stage. For the first scenario for the restricted cubic splines, we placed four knots, following Harrell's suggestion to use the 5%, 35%, 65% and 95% quantiles of BMI and for natural B-splines four equidistant knots (two inner knots plus the boundaries per study). For the second and third scenario, for the restricted cubic splines we placed three knots at 10%, 50% and 90% quantiles of BMI and for natural B-splines three equidistant knots (one inner knot plus the boundaries per study. For all scenarios for P-splines we used 17 equidistant knots (15 inner knots plus the boundaries per study). For the penalised splines (P-splines and smoothing splines) the tuning parameter λ was selected through a 'leave one out' GCV process. Instead of using pointwise meta-analysis per distinct value of X, the functional shapes can also be pooled using multivariate meta-analysis. This approach, as proposed by Gasparrini et al.,23 pools the set of regression coefficients estimated in the first stage, accounting for their within- and (if applicable) between-study correlation, using a fixed or random effects multivariate meta-analysis approach. Assumptions are a normal distribution, a constant between-studies covariance matrix, and a normal distribution for the random effect, which implies that this is symmetrical and thus does not allow for heavy or light tails.48 The use of penalised splines in multivariate meta-analysis still wants further research, e.g. with respect to possible over-penalisation in small studies and low/no data areas within studies, and is therefore beyond our current scope. Note that in order to pool the results of the first stage, each study should provide the same set of coefficients, estimated in the same domain of X. Therefore, in order to apply multivariate meta-analysis, the basis functions for the splines in the individual studies should be of the same degree, and also defined on the same intervals across studies, using the same knot positions. In case of different ranges of X across studies, the use of common positions for the knots may leave some coefficients inestimable in some studies and meta-analysing them may cause complications.23 A solution is to conduct data augmentation as a preliminary step. Data augmentation as described by White et al.49 and Riley et al.26 refers to the generation of pseudo data beyond the per study boundaries of X, with minimal weight and arbitrary outcome. Note that in multivariate meta-analysis a careful specification of the knots may be required. Use of a large number of knots may cause convergence issues during the second stage, in case of a rank-deficient variance covariance matrix. The multivariate meta-analysis algorithm proceeds as follows: We applied multivariate meta-analysis in combination with regression splines in all three scenarios. To do so we performed data augmentation as a preliminary step26, 49 in the second and third scenario. This way all studies had curves estimated over the full range of BMI. In stage 1, we fitted restricted cubic spline and B-spline transformations of BMI both as main effects and as interactions with the treatment per study. For the restricted cubic spline transformations, we used five knots, following Harrell's suggestion to use the 5%, 27.5%, 50%, 72.5% and 95% quantiles of BMI, for natural B-splines four equidistant knots (two inner knots plus the boundaries per study). Note that we positioned the knots over the full domain of BMI. Subsequently, we pooled the estimated coefficients using a random-effects meta-analysis with the REML estimation method. We calculated regression lines per treatment arm by multiplying the design (or model) matrix with the pooled coefficients. Absolute risk differences were calculated by subtracting the pooled mortality risks, conditional on the covariables, of the tr