Perform a time-course Model-Based Network Meta-Analysis (MBNMA)"

In MBNMAtime: Run Time-Course Model-Based Network Meta-Analysis (MBNMA) Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  eval=rmarkdown::pandoc_available("1.12.3")
)

library(MBNMAtime)
library(rmarkdown)
library(knitr)
library(dplyr)
#load(system.file("extdata", "vignettedata.rda", package="MBNMAtime", mustWork = TRUE))

Analysis using mb.run()

MBNMA models are fitted using mb.run() and they model treatment effects over time (e.g. days, weeks, months). This can just as easily be performed on datasets with many different treatments (network meta-analysis) as it can on datasets comparing only two treatments (pairwise meta-analysis) - the syntax is the same.

An object or class("mb.network") must be provided as the data for mb.run(). The key arguments within mb.run() involve specifying the functional form used to model the time-course, and the time-course parameters that comprise that functional form.

Time-course functions

A number of different time-course functions can be fitted within MBNMAtime and the specific forms of the time-course parameters are defined by arguments within these functions, and this allows for a wide variety of parameterizations and time-course shapes.

For further details check the help files for each function (e.g. ?tloglin()). These functions, are then used as inputs for the fun argument in mb.run().

• tloglin() - Log-linear function
• titp() - Integrated Two-Component Prediction (ITP) function
• temax() - Emax function
• tpoly() - Polynomial function (e.g. linear, quadratic)
• tfpoly() - Fractional polynomial function, as proposed previously for time-course NMA by Jansen [-@jansen2015].
• tspline() - Spline functions (includes B-splines, restricted cubic splines, natural splines and piecewise linear splines)
• tuser() - A time-course function that can be explicitly defined by the user

Time-course parameters within time-course functions are each defined by two arguments:

pool is used to define the approach used for the pooling of a given time-course parameter and can either of:

• "rel" indicates that relative effects (or mean differences) should be pooled for this time-course parameter. This preserves randomisation within included studies, are likely to vary less between studies (only due to effect modification), and allow for testing of consistency between direct and indirect evidence. Pooling follows the general approach for Network Meta-Analysis proposed by Lu and Ades [-@lu2004].
• "abs" indicates that study arms should be pooled across the whole network for this time-course parameter independently of assigned treatment. This implies using a single absolute value across the network for this time-course parameter, and may therefore be making strong assumptions of similarity between studies.

method is used to define the model used for meta-analysis for a given time-course parameter and can take either of:

• "common" implies that all studies estimate the same true effect (sometimes called a "fixed effect" meta-analysis)
• "random" implies that all studies estimate a separate true effect, but that each of these true effects vary randomly around a true mean effect. This approach allows for modelling of between-study heterogeneity.

Specifying pooling relative effects on all time-course parameters would imply performing a contrast-based synthesis, whereas specifying pooling absolute effects on all of them would imply performing an arm-based synthesis. There has been substantial discussion in the literature regarding the strengths and limitations of both these approaches [@dias2016; @hong2016; @karahalios2017].

Additional arguments within the function may also be used to specify the degree (e.g. for polynomials) or the number of knots or knot placement for splines.

Splines and knots

For a more flexible time-course shape, various different splines can be fitted to the data by using tspline(). This model is very flexible and can allow for a variety of non-monotonic time-course relationships, though parameters can be difficult to interpret and the resulting time-course shape is often best visualised by calculating and plotting predictions Prediction.

To fit this model, the number/location of knots (the points at which the different spline pieces meet) should be specified. If a single number is given, it represents the the number of knots to be equally spaced across the range of times within the dataset. Alternatively several probabilities can be given that represent the quantiles of the time range at which knots should be located. Note that by default, a boundary knot will be placed at the maximum of the range of time-points to limit the function extrapolating to extreme values.

tspline(type="bs", knots=3)
# ...is equivalent to 
tspline(type="bs", knots=c(0.25,0.5,0.75))

Choosing a time-course function

Choice of time-course function involves a combination of data-driven, statistical information (plots, model fit statistics) and clinical/biological plausibility. timeplot() and binplot() functions can be used to observe the how the data differs at different time-points. This can help to understand, for example, whether the data monotonically increases or decreases, or whether there is a turning point or multiple turning-points within the time-course, and at which follow-up time(s) this might occur. This can also be guided by a clinical/biological knowledge/expectation of the shape of the time-course outside of the range of time-points for which data are available.

However, once several potential candidate time-course functions have been identified, selection between them (and of the specific specification of the time-course parameters in terms of absolute or relative effects) is best achieved using model fit statistics (e.g. Deviance Information Criterion (DIC), residual deviance).

Output

mb.run() returns an object of class(c("mbnma", "rjags")). summary() provides summary estimates of posterior densities for different parameters in the model, with some explanation regarding the way in which the model has been defined. Estimates are automatically reported for parameters of interest depending on the model specification (unless otherwise specified in parameters.to.save). Nodes that are automatically monitored (if present in the model) have the following interpretation:

Parameters modelled using relative effects

If pooling is relative (e.g. pool.1="rel") for a given parameter then the named parameter (e.g. emax) or a numbered d parameter (e.g. d.1) corresponds to the pooled relative effect (or mean difference) for a given treatment compared to the network reference treatment for this time-course parameter.

sd. followed by a named (e.g. emax, beta.1) is the between-study SD (heterogeneity) for relative effects, reported if pooling for a time-course parameter is relative (e.g. pool.1="rel") and the method for synthesis is random (e.g. method.1="random).

If class effects are modelled, parameters for classes are represented by the upper case name of the time-course parameter they correspond to. For example if class.effect=list(emax="random"), relative class effects will be represented by EMAX. The SD of the class effect (e.g. sd.EMAX, sd.BETA.1) is the SD of treatments within a class for the time-course parameter they correspond to.

Parameters modelled using absolute effects

If pooling is absolute (e.g. pool.1="abs") for a given parameter then the named parameter (e.g. emax) or a numbered beta parameter (e.g. beta.1) corresponds to the estimated absolute effect for this time-course parameter.

For an absolute time-course parameter if the corresponding method is common (e.g. method.1="common") the parameter corresponds to a single common parameter estimated across all studies and treatments. If the corresponding method is random (e.g. method.1="random") then parameter is a mean effect around which the study-level absolute effects vary with SD corresponding to sd. followed by the named parameter (e.g. sd.emax, sd.beta.1).

Other model parameters

rho is the correlation coefficient for correlation between time-points. Its interpretation will differ depending on the covariance structure specified in covar.

totresdev is residual deviance of the model and deviance is the deviance of the model. Model fit statistics for pD (effective number of parameters) and DIC (Deviance Information Criterion) are also reported, with an explanation as to how they have been calculated.

Examples

An example MBNMA of the alogliptin dataset using a linear time-course function and common treatment effects that pool relative effects and assumes consistency between direct and indirect evidence can be performed as follows:

# Prepare data using the alogliptin dataset
network.alog <- mb.network(alog_pcfb, reference = "placebo")

# Run a linear time-course MBNMA
mbnma <- mb.run(network.alog, fun=tpoly(degree=1, pool.1="rel", method.1="common"))

summary(mbnma)

For this model, the d.1 parameters correspond to the 1st polynomial coefficient, and therefore are the linear gradient of the response over time for each treatment versus placebo - i.e. the mean difference for the change in efficacy for each treatment versus placebo. However, note that the residual deviance of the model is very high, suggesting (as we might expect) that this linear time-course function is a poor fit.

We may want to fit a more complex time-course function with two time-course parameters, such as an Emax function, yet limitations in the data might require that we make an assumption that one of the parameters does not vary by treatment. We can specify this by setting pool to be equal to "abs" for any parameters we choose.

# Run an Emax time-course MBNMA with two parameters
mbnma <- mb.run(network.alog, fun=temax(
  pool.emax = "rel", method.emax="common",
  pool.et50 = "abs", method.et50="common"
))

summary(mbnma)

In this case, the parameters are named following the Emax function specification. emax corresponds to the maximum effect for each treatment versus placebo (interpretable as a mean difference versus placebo), whereas et50 is the log of the time at which 50% of the maximum response is achieved, across all treatments in the network. This assumes conditional constancy of absolute effects for this time-course parameter, which is typically a strong assumption. However, if there were limited data with which to inform this parameter (e.g. at earlier time-points) then such an assumption might be necessary, with the caveat that interpolation of response at time-points informed by this parameter may be more susceptible to bias. Further exploration of the degree of data required for reliable estimation of time-course parameters is given in Pedder et al. [-@pedder2020].

Additional model specification with mb.run()

Correlation between time points

Within-study correlation between time points can easily be modeled using mb.run(), though this requires some additional considerations. The simplest approach is to incorporate correlation by using a variance adjustment [@jansen2015]. This avoids the need to use a multivariate normal likelihood (which is slow to run), and it assumes a common correlation between neighbouring time-points. This is achieved by using the argument covar="varadj", which is the default in mb.run().

There are two alternative covariance structures can be modelled, though these require fitting a multivariate normal likelihood and therefore take longer to run. covar="CS" specifies fitting a Compound Symmetry covariance structure, whilst covar="AR1" specifies fitting an autoregressive AR1 covariance structure to the multivariate normal likelihood used for modelling the correlation between multiple time points within a study [@covstruct].

However, in addition to this, it's also necessary to specify a value for rho, and this can be assigned in one of two ways:

• Given as string representing a JAGS prior distribution [@jagsmanual], which indicates that the correlation should be estimated from the data. For example, to specify a prior that the correlation between time-points will be between 0 and 1 with equal probability you could set rho="dunif(0,1)".
• Given as a single numeric value, which indicates that the correlation should be fixed to that value. For example, this value could be estimated externally from another study using Individual Participant Data. This could also be used to run a deterministic sensitivity analysis using different fixed values of rho.

# Using the COPD dataset
network.copd <- mb.network(copd)

# Run an log-linear time-course MBNMA 
# that accounts for correlation between time points using variance adjustment
mbnma <- mb.run(network.copd, 
                fun=tloglin(pool.rate="rel", method.rate="random"),
                rho="dunif(0,1)", covar="varadj")

It is important to note that the covariance matrix must be positive semi-definite. This may mean that in order to satisfy this requirement for particular covariance matrix structures, the values that rho can take are limited. rho must always be bounded by -1 and 1, but even within this range some negative values for rho can result in a non positive matrix, which can lead to an error in the evaluation of the multivariate likelihood. If so, it may be necessary to further restrict the prior distribution.

Link function

Time-course MBNMA can only be used for continuous outcomes or those that can be summarised as continuous outcomes (e.g. binary data can be converted to log-odds for inclusion in the model). Typically this means that users will analyse data using an identity link function, the default given to the link argument in mb.run()m which assumes an additive treatment effect (e.g. mean difference).

However, by specifying link="log" a user can model a log link and therefore assume a multiplicative treatment effect. For continuous data this models the treatment effect as a Ratio of Means (RoM) [@friedrich2011]. This also provides an advantage as the treatment effect is scale independent (i.e. studies measuring the same outcome using different measurement scales can be analysed simultaneously). However, within-study treatment effects must all be positive (if all negative then they can all be converted to positive and the treatment effects are inverted), and change from baseline measures must be adjusted so that they are also expressed as RoMs (log(follow-up) - log(baseline)) to avoid combining additive and multiplicative assumptions within the same analysis. Note that assuming a ratio effect in time-course models implies very different effects over time to a model on the natural scale. For example, a linear time-course model on the log scale implies an exponential relationship on the natural scale.

An alternative approach for modelling a measurement scale-independent treatment effect whilst still assuming additive treatment effects is to perform the analysis using Standardised Mean Differences (SMD). Whilst not strictly a different link function, this can be specified using link="smd".

By default, MBNMAtime standardises treatment effects using the pooled standard deviation (SD) at baseline in each study ("study-specific" SDs). However, a more robust approach to minimise bias from estimation of study-specific SD is to use a "reference" SD that is the same for each scale included in the dataset. Such a reference could be estimated from all the included studies within the dataset that report outcomes on the scale (an "internal reference SD"), but the best approach is to obtain a SD from an observational study that is specific to the target population of interest.

Once this has been obtained for each scale, it can be included in the original data frame as an additional variable, standsd. For the analysis to be meaningful, stansd should be the same for all arms within a study and for all studies that report the same outcome. If this has been included in the dataset used to create an "mb.network" object, then sdscale=TRUE argument can be specified in mbnma.run() to indicate that this variable should be used to standardise SMDs.

Note that for outputting results from a model using SMDs, results should be back-transformed to a clinically relevant scale to make robust interpretation. This can be done by multiplying the effect estimates or predictions by the standard deviation of the relevant scale in the population of interest. For further details of analysis of continuous data that include discussion of both RoM and SMD see [@GMD2].

Class effects

Shared effects between treatments within the network can be modelled using class effects. This requires assuming that different treatments have some sort of shared class effect, perhaps due to different (yet clinically similar) doses of the same agent or different treatments with a similar mechanism of action. One advantage of this is that class effects can be used to connect relative effects between treatments in a network that would be disconnected at the treatment level, but can be connected via classes at the class level. However, it is important to ensure that such an effect is clinically justifiable, as making these assumptions risks introducing heterogeneity/inconsistency.

Class effects can only be applied to time-course parameters which vary by treatment (pool="rel"), and class effects are modelled separately for each time-course parameter.

In mb.run() class effects are specified as a list, in which each element is named by the time-course parameter on which it should be modelled. The class effect for each time-course parameter can be either "common", in which the effects for each treatment within the same class are constrained to a common class effect, or "random", in which the effects for each treatment within the same class are assumed to be randomly distributed around a shared class mean.

# Create network object of gout dataset
network.gout <- mb.network(goutSUA_CFBcomb)

# Run a B-spline time-course MBNMA with a knot at 0.2 times the max follow-up
# Common class effect on beta.2, the 2nd spline coefficient
mbnma <- mb.run(network.gout, 
                fun=tspline(type="bs", knots=c(0.2),
                            pool.1 = "rel", method.1="common",
                            pool.2="rel", method.2="random"),
                class.effect = list(beta.2="common"))

summary(mbnma)

Mean class effects are given in the output as D.2 parameters. These can be interpreted as the relative effect of each class versus the Plac (Placebo), for the 2nd spline coefficient (beta.2). Note the number of D.2 parameters is therefore equal to the number of classes defined in the dataset.

Additional arguments

Several additional arguments can be given to mb.run() that require further explanation.

Priors

Default vague priors for the model are as follows:

$$ \begin{aligned} &\alpha_{i} \sim N(0,10000)\ &\boldsymbol{\mu}{i} \sim N(0,10000)\ &\boldsymbol{d}{t} \sim N(0,10000)\ &beta_{\phi} \sim N(0,10000)\ &D_{\phi,c} \sim N(0,1000)\ &\tau_{\phi} \sim N(0,400) \text{ limited to } x \in [0,\infty]\ &\tau^D_{\phi} \sim N(0,400) \text{ limited to } x \in [0,\infty]\ \end{aligned} $$

• $\alpha_i$ is the response at time=0 in study $i$
• $\mu_i$ is a vector of study reference effects for each time-course parameter in study $i$. Where only a single time-course parameter is modelled using relative effects the prior is defined as $\mu_{i} \sim N(0,10000)$.
• $\boldsymbol{d}{t}$ is a vector of pooled relative effects for treatment $t$ whose length is the number of time-course parameters in the model. Where only a single time-course parameter is modelled using relative effects the prior is defined as $d{t} \sim N(0,10000)$.
• $\beta_{\phi}$ is the absolute effect for time-course parameter $\phi$ modelled independently of treatment
• $D_{\phi,c}$ is the class relative effect for time-course parameter $\phi$ in class $c$
• $\tau_{\phi}$ is the between-study SD for time-course parameter $\phi$
• $\tau^D_{\phi}$ is the within-class SD for time-course parameter $\phi$

Users may wish to change these, perhaps in order to use more/less informative priors or different prior distributions (e.g. log-normal prior rather than a truncated normal for ET50 in an Emax model). However, it may also be because the default prior distributions in some models can lead to errors when compiling/updating models if the prior includes extremely implausible values.

This can be more likely for certain types of models. For example some prior distributions may generate results that are too extreme for JAGS to compute, such as for time-course parameters that are powers (e.g. Emax functions with a Hill parameter or power parameters in fractional polynomials).

If the model fails during compilation/updating (i.e. due to a problem in JAGS), mb.run() will generate an error and return a list of arguments that mb.run() used to generate the model. Within this (as within a model that has run successfully), the priors used by the model (in JAGS syntax) are stored within "model.arg".

In this way a model can first be run with vague priors and then rerun with different priors, perhaps to allow successful computation, perhaps to provide more informative priors, or perhaps to run a sensitivity analysis with different priors.

To change priors within a model, a list of replacements can be provided to priors in mb.run(). The name of each element is the name of the parameter to change (without indices) and the value of the element is the JAGS distribution to use for the prior. See the JAGS Manual [-@jagsmanual] for syntax details regarding specifying distributions. This can include censoring or truncation if desired. Only the priors to be changed need to be specified - priors for parameters that aren't specified will take default values. Note that in JAGS, normal distributions are specified using precision (1/variance) rather than SD.

For example, we may wish to specify a tighter prior for the between-study SD:

mbnma <- mb.run(network.copd, 
                fun=tloglin(pool.rate="rel", method.rate="random"),
                priors=list(rate="dnorm(0,2) T(0,)"))

Different prior distributions can be assigned for different indices of a parameter by specifying the list element for a parameter as a character vector. This allows (for example) for the user to fit specific priors for specific treatments. The length of this vector must be equal to the number of indices of the parameter. The ordering will also be important - for example for treatment-specific priors the order of the elements within the vector must match the order of the treatments in the network.

For example we might have differnt beliefs about the long-term efficacy of a treatment for which there is no long-term data available in the dataset. In the COPD dataset we have longer term data (up to 52 weeks) on Tiotropium, but much shorter follow-up data (up to 26 weeks) on Aclidinium.

We might believe (e.g. based on clincial opinion) that the efficacy of Aclidinium returns towards baseline at longer follow-up. We could model this using a B-spline and providing informative priors only to the parameters controlling the spline for Aclidinium at later follow-up:

# Define informative priors for spline parameters
spline.priors <- list(
  d.3 = c(
  Aclidinium="dnorm(-0.5, 100)", 
  Tiotropium="dnorm(0, 0.0001)"
  ),
  d.4 = c(
  Aclidinium="dnorm(0, 100)", 
  Tiotropium="dnorm(0, 0.0001)"
  ))

# Using the COPD dataset with a B-spline MBNMA
mbnma <- mb.run(network.copd, fun=tspline(degree=2, knots=c(0.1,0.5)),
                priors=spline.priors)

# Predict and plot time-course relative effect
pred <- predict(mbnma)
plot(pred)

As can be seen from the predicted time-course, using informative priors for Aclidinium in this way allows us to predict it's efficacy at longer-term follow-up than the data alone can inform.

pD (effective number of parameters)

The default value for pd in mb.run() is pd="pv", which uses the rapid approach automatically calculated in the R2jags package as pv = var(deviance)/2. Whilst this is easy to calculate, it is only an approximation to the effective number of parameters, and may be numerically unstable [@gelman2003]. However, it has been shown to be reliable for model comparison in time-course MBNMA models in a simulation study [@pedder2020].

A more reliable method for estimating pd is pd="pd.kl", which uses the Kullback-Leibler divergence [@plummer2008]. This is more reliable than the default method used in R2jags for calculating the effective number of parameters from non-linear models. The disadvantage of this approach is that it requires running additional MCMC iterations, so can be slightly slower to calculate.

A commonly-used approach in Bayesian models for calculating pD is the plug-in method (pd="plugin") [@spiegelhalter2002]. However, this can sometimes result in negative non-sensical values due to skewed posterior distributions for deviance contributions that can arise when fitting non-linear models.

Finally, pD can also be calculated using an optimism adjustment (pd="popt") which allows for calculation of the penalized expected deviance [@plummer2008]. This adjustment allows for the fact that data used to estimate the model is the same as that used to assess its parsimony. As for pd="pd.kl", it also requires running additional MCMC iterations.

Correlation between time-course parameters

Since MBNMAtime version 0.2.2, mb.run() no longer automatically models correlation between time-course parameters modeled using relative effects.

However, this can be enabled by setting corparam=TRUE). Time-course parameters are typically correlated and this allows information on each parameter to help inform the other(s). The correlation is modeled using a multivariate normal distribution whose covariance matrix $\Sigma_t$ is modelled using a correlation defined by the parameter rhoparam. In previous versions, an inverse Wishart prior was used, though this has since been considered to overly constrain the model.

Arguments to be sent to JAGS

In addition to the arguments specific to mb.run() it is also possible to use any arguments to be sent to R2jags::jags(). Most of these relate to improving the performance of MCMC simulations in JAGS. Some of the key arguments that may be of interest are:

• n.chains The number of Markov chains to run (default is 3)
• n.iter The total number of iterations per MCMC chain
• n.burnin The number of iterations that are discarded to ensure iterations are only saved once chains have converged
• n.thin The thinning rate which ensures that results are only saved for 1 in every n.thin iterations per chain. This can be increased to reduce autocorrelation in MCMC samples

Model Selection

Detailed description of model selection based on statistical measures such as Deviance Information Criterion (DIC) and residual deviance is outside the scope of this vignette. However, the following approach for model identification and selection is recommended in Pedder et al. [-@pedder2020], which also gives details on model fit statistics used for comparison:

1. Identify candidate time-course functions based on observed data and clinical/biological reasoning
2. Compare candidate time-course functions fitted with common relative treatment effects on all time-course parameters
3. If no candidate time-course functions converge successfully, absolute effects can be modeled on parameters for which convergence is problematic
4. Compare selected common effects model to model(s) with random effects on different time-course parameters
5. Compare model fitted with a univariate likelihood to one fitted with a multivariate likelihood

Finally the validity of the consistency assumption should be explored in the selected final model (see [Consistency Testing]).

MCMC Convergence

MBNMAtime runs Bayesian models in JAGS, which uses Markov Chain Monte Carlo (MCMC) simulation [@jags]. However, the validity of results is reliant on the MCMC chains having converged successfully on the posterior densities of all parameters in the model. For highly parameterised models run on relatively limited data, as is often the case with MBNMA models, convergence can often be an challenge. Note that convergence is necessary to be able to compare models and evaluate model fit. However, successful convergence does not imply good model fit.

A full explanation of how to facilitate and to check for convergence is outside the scope of this vignette, but below are some simple steps for checking convergence. None of these steps alone ensures convergence has been successful, but interpretation of them jointly can provide strong evidence of it.

• Rhat values close to 1 (\<1.1 is considered acceptable by some). These are shown for monitored parameters in the summary statistics table when an "mbnma" or "nma" object is printed.
• Trace plots that have a "fuzzy caterpillar" look that shows good mixing of MCMC chains
• Density plots that show a smooth posterior distribution that is similar across all chains
• Autocorrelation plots that show low autocorrelation between neighbouring MCMC iterations (i.e. low correlation at higher lags)

An HTML document with all the above convergence plots can easily be generated for all parameters in the model simultaneously using mcmcplots::mcmcplot().

Two steps may improve convergence when using MBNMAtime. One step is to run models for more iterations (this can be done using the n.iter argument in mb.run()). It can take time for MCMC chains to converge, particularly for non-linear models with limited data. It is important to note that chains should only be monitored after they have converged - increasing the number of burn-in iterations ensures that this is the case (using the n.burnin argument in mb.run()). Another method to improve convergence is by providing more information to the model via informative priors.

For a detailed review of MCMC convergence assessment see Sinharay [-@sinharay].

References

MBNMAtime documentation built on Oct. 14, 2023, 5:08 p.m.