Package 'multinma'

Description

An R package for performing network meta-analysis and network meta-regression with aggregate data, individual patient data, or mixtures of both.

Details

Network meta-analysis (NMA) combines (aggregate) data from multiple studies on multiple treatments in order to produce consistent estimates of relative treatment effects between each pair of treatments in the network (Dias et al. 2011).

Network meta-regression (NMR) extends NMA to include covariates, allowing adjustment for differences in effect-modifying variables between studies (Dias et al. 2011). NMR is typically performed using aggregate data (AgD), which lacks power and is prone to ecological bias. NMR with individual patient data (IPD) is the gold standard, if data are available.

Multilevel network meta-regression (ML-NMR) allows IPD and AgD to be incorporated together in a network meta-regression (Phillippo et al. 2020; Phillippo 2019). As in IPD NMR, an individual-level regression model is defined. AgD studies are then fitted by integrating the individual-level model over the respective covariate distributions. This correctly links the two levels of the model (instead of "plugging in" mean covariate values), avoiding aggregation bias. Population-adjusted treatment effects (Phillippo et al. 2016) can be produced for any study population in the network, or for an external target population.

Models are estimated in a Bayesian framework using Stan (Carpenter et al. 2017). Quasi-Monte Carlo numerical integration based on Sobol' sequences is used for the integration in ML-NMR models, with a Gaussian copula to account for correlations between covariates (Phillippo et al. 2020; Phillippo 2019).

Getting Started

A good place to start is with the package vignettes which walk through example analyses, see vignette("vignette_overview") for an overview. The series of NICE Technical Support Documents on evidence synthesis gives a detailed introduction to network meta-analysis:

Dias S, Welton NJ, Sutton AJ, Caldwell DM, Lu G, Reken S, Ades AE (2011). “NICE DSU Technical Support Documents 1-7: Evidence Synthesis for Decision Making.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Multilevel network meta-regression is set out in the following methods paper:

Phillippo DM, Dias S, Ades AE, Belger M, Brnabic A, Schacht A, Saure D, Kadziola Z, Welton NJ (2020). “Multilevel Network Meta-Regression for population-adjusted treatment comparisons.” Journal of the Royal Statistical Society: Series A (Statistics in Society), 183(3), 1189–1210. doi:10.1111/rssa.12579.

Author(s)

Maintainer: David M. Phillippo [email protected] (ORCID)

References

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, Riddell A (2017). “Stan: A Probabilistic Programming Language.” Journal of Statistical Software, 76(1). doi:10.18637/jss.v076.i01.

Dias S, Sutton AJ, Welton NJ, Ades AE (2011). “NICE DSU Technical Support Document 3: Heterogeneity: subgroups, meta-regression, bias and bias-adjustment.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Phillippo DM (2019). Calibration of Treatment Effects in Network Meta-Analysis using Individual Patient Data. Ph.D. thesis, University of Bristol. Available from https://research-information.bris.ac.uk/.

Phillippo DM, Ades AE, Dias S, Palmer S, Abrams KR, Welton NJ (2016). “NICE DSU Technical Support Document 18: Methods for population-adjusted indirect comparisons in submission to NICE.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Phillippo DM, Dias S, Ades AE, Belger M, Brnabic A, Schacht A, Saure D, Kadziola Z, Welton NJ (2020). “Multilevel Network Meta-Regression for population-adjusted treatment comparisons.” Journal of the Royal Statistical Society: Series A (Statistics in Society), 183(3), 1189–1210. doi:10.1111/rssa.12579.

See Also

Useful links:

• https://dmphillippo.github.io/multinma/

• https://github.com/dmphillippo/multinma

• Report bugs at https://github.com/dmphillippo/multinma/issues

Description

The .default() function is used internally to mark certain values as default, so that the user may be notified when default values are being used. For example, choosing a default reference treatment for a network, or using default prior distributions. The function .is_default() checks whether an argument/object is set to a default value. Neither of these functions are intended to be called by the user.

Usage

.default(x = list())

.is_default(x)

Arguments

Value

For .default(), an identical object with additional attribute .default. For .is_default(), a logical value (TRUE or FALSE).

Description

The Stan control argument adapt_delta sets the target average acceptance probability for the No-U-Turn Sampler (NUTS) used by Stan.

Details

The default value of adapt_delta used by nma() is 0.8 for fixed effect models, and 0.95 for random effects models.

You should not need to change adapt_delta unless you see a warning message about divergent transitions. Increasing adapt_delta from the default to a value closer to 1 means that Stan will use a smaller step size, making sampling slower but more robust, and resulting in fewer divergent transitions.

For more details see the Stan documentation available from https://mc-stan.org/users/documentation/.

Description

The add_integration() generic creates Quasi-Monte Carlo numerical integration points using a Gaussian copula and Sobol' sequences, as described in Phillippo et al. (2020). Methods are available for networks stored in nma_data objects, and for data frames. The function unnest_integration() unnests integration points stored in a data frame, to aid plotting or other exploration.

Usage

add_integration(x, ...)

## Default S3 method:
add_integration(x, ...)

## S3 method for class 'data.frame'
add_integration(
  x,
  ...,
  cor = NULL,
  cor_adjust = NULL,
  n_int = 64L,
  int_args = list()
)

## S3 method for class 'nma_data'
add_integration(
  x,
  ...,
  cor = NULL,
  cor_adjust = NULL,
  n_int = 64L,
  int_args = list()
)

unnest_integration(data)

Arguments

Details

The arguments passed to ... specify distributions for the covariates. Argument names specify the name of the covariate, which should match a covariate name in the IPD (if IPD are present). The required marginal distribution is then specified using the function distr().

The argument cor_adjust specifies how the correlation matrix given by cor (or computed from the IPD if cor = NULL) is adjusted to obtain the correlation matrix for the Gaussian copula, using the formulae in Xiao and Zhou (2018).

• cor_adjust = "spearman" should be used when the correlations cor have been computed using Spearman's rank correlation. Correlations between continuous covariates will be reproduced exactly by the integration points. Correlations between discrete covariates will be reproduced approximately. This is the default when cor = NULL and correlations are calculated from the IPD studies.

• cor_adjust = "pearson" should be used when the correlations cor have been computed using Pearson's product-moment correlation. Correlations between Normal covariates will be reproduced exactly by the integration points, all others will be reproduced approximately. Correlations between discrete covariates will be reproduced approximately (and identically to cor_adjust = "spearman"). This is the default when cor is provided by the user, since cor() defaults to method = "pearson" and Pearson correlations are most likely reported in published data. However, we recommend providing Spearman correlations (e.g. from cor(., method = "spearman")) and using cor_adjust = "spearman" where possible.

• cor_adjust = "none" allows the user to specify the correlation matrix for the Gaussian copula directly; no adjustment is applied.

• cor_adjust = "legacy" is also available, which reproduces exactly the behaviour from version 0.3.0 and earlier. This is similar to cor_adjust = "none", but unadjusted Spearman correlations are used if cor = NULL.

When adding integration points to a network object the correlation matrix used is stored in ⁠$int_cor⁠, and the copula correlation matrix and adjustment used are stored as attributes of ⁠$int_cor⁠. If this correlation matrix is passed again to add_integration() (e.g. to reuse the correlations for an external target population) this will be detected, and the correct setting for cor_adjust will automatically be applied.

Value

For the nma_data method, an object of class nma_data. For the data.frame method, the input data frame is returned (as a tibble) with an added column for each covariate (prefixed with ".int_"), containing the numerical integration points nested as length-n_int vectors within each row. For unnest_integration(), a data frame with integration points unnested.

References

Phillippo DM, Dias S, Ades AE, Belger M, Brnabic A, Schacht A, Saure D, Kadziola Z, Welton NJ (2020). “Multilevel Network Meta-Regression for population-adjusted treatment comparisons.” Journal of the Royal Statistical Society: Series A (Statistics in Society), 183(3), 1189–1210. doi:10.1111/rssa.12579.

Xiao Q, Zhou S (2018). “Matching a correlation coefficient by a Gaussian copula.” Communications in Statistics - Theory and Methods, 48(7), 1728–1747. doi:10.1080/03610926.2018.1439962.

Examples

## Plaque psoriasis ML-NMR - network setup and adding integration points
# Set up plaque psoriasis network combining IPD and AgD
library(dplyr)
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)

pso_ipd <- pso_ipd %>%
  mutate(# Variable transformations
    bsa = bsa / 100,
    prevsys = as.numeric(prevsys),
    psa = as.numeric(psa),
    weight = weight / 10,
    durnpso = durnpso / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker"),
    # Check complete cases for covariates of interest
    complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>%
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker")
  )

# Exclude small number of individuals with missing covariates
pso_ipd <- filter(pso_ipd, complete)

pso_net <- combine_network(
  set_ipd(pso_ipd,
          study = studyc,
          trt = trtc,
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd,
              study = studyc,
              trt = trtc,
              r = pasi75_r,
              n = pasi75_n,
              trt_class = trtclass)
)

# Print network details
pso_net

# Add integration points to the network
pso_net <- add_integration(pso_net,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  n_int = 64)


## Adding integration points to a data frame, e.g. for prediction
# Define a data frame of covariate summaries
new_agd_int <- data.frame(
  bsa_mean = 0.6,
  bsa_sd = 0.3,
  prevsys = 0.1,
  psa = 0.2,
  weight_mean = 10,
  weight_sd = 1,
  durnpso_mean = 3,
  durnpso_sd = 1)

# Adding integration points, using the weighted average correlation matrix
# computed for the plaque psoriasis network
new_agd_int <- add_integration(new_agd_int,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  cor = pso_net$int_cor,
  n_int = 64)

# Here, since we reused the correlation matrix pso_net$int_cor from the
# network, the correct setting of cor_adjust = "spearman" is automatically
# applied

new_agd_int

Description

Samples (post warm-up) from a stan_nma model object can be coerced into an array, matrix, or data frame.

Usage

## S3 method for class 'stan_nma'
as.array(x, ..., pars, include = TRUE)

## S3 method for class 'stan_nma'
as.data.frame(x, ..., pars, include = TRUE)

## S3 method for class 'stan_nma'
as_tibble(x, ..., pars, include = TRUE)

## S3 method for class 'stan_nma'
as.tibble(x, ..., pars, include = TRUE)

## S3 method for class 'stan_nma'
as.matrix(x, ..., pars, include = TRUE)

Arguments

Value

The as.array() method produces a 3D array [Iteration, Chain, Parameter] containing posterior samples of each parameter (as class mcmc_array). This has the side effect of enabling bayesplot functions to seamlessly work on stan_nma objects.

The as.data.frame() method produces a data frame containing posterior samples of each parameter, combined over all chains.

The as.matrix() method produces a matrix containing posterior samples of each parameter, combined over all chains.

Description

The method as.igraph() converts nma_data objects into the form used by the igraph package. The method as_tbl_graph() converts nma_data objects into the form used by the ggraph and tidygraph packages.

Usage

## S3 method for class 'nma_data'
as.igraph(x, ..., collapse = TRUE)

## S3 method for class 'nma_data'
as_tbl_graph(x, ...)

Arguments

Value

An igraph object for as.igraph(), a tbl_graph object for as_tbl_graph().

Examples

# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net

# Convert to igraph object
igraph::as.igraph(smk_net)  # Edges combined by default
igraph::as.igraph(smk_net, collapse = FALSE)  # Without combining edges

# Convert to tbl_graph object
tidygraph::as_tbl_graph(smk_net)  # Edges combined by default
tidygraph::as_tbl_graph(smk_net, collapse = FALSE)  # Without combining edges

Description

Attempt to turn an object into a stanfit object.

Usage

as.stanfit(x, ...)

## S3 method for class 'stan_nma'
as.stanfit(x, ...)

## Default S3 method:
as.stanfit(x, ...)

Arguments

Value

A stanfit object.

Description

Data frame containing the results of 26 trials comparing 17 treatments in 4 classes for the prevention of stroke in patients with atrial fibrillation (Cooper et al. 2009). The data are the corrected versions given by van Valkenhoef and Kuiper (2016).

Usage

atrial_fibrillation

Format

A data frame with 63 rows and 11 variables:

studyc

study name

studyn

numeric study ID

trtc

treatment name

trtn

numeric treatment code

trt_class

treatment class

r

number of events

n

sample size

E

person-years at risk

stroke

proportion of individuals with prior stroke

year

year of study publication

followup

mean length of follow-up (years)

References

Cooper NJ, Sutton AJ, Morris D, Ades AE, Welton NJ (2009). “Addressing between-study heterogeneity and inconsistency in mixed treatment comparisons: Application to stroke prevention treatments in individuals with non-rheumatic atrial fibrillation.” Statistics in Medicine, 28(14), 1861–1881. doi:10.1002/sim.3594.

van Valkenhoef G, Kuiper J (2016). gemtc: Network Meta-Analysis Using Bayesian Methods. R package version 0.8-2, https://CRAN.R-project.org/package=gemtc.

Description

Data frame containing the results of 13 trials comparing BCG vaccination to no vaccination for preventing tuberculosis (TB) (Dias et al. 2011; Berkey et al. 1995). The numbers of individuals diagnosed with TB in each arm during the study follow-up period are recorded. The absolute degrees latitude at which the study was conducted are also recorded.

Usage

bcg_vaccine

Format

A data frame with 26 rows and 6 variables:

studyn

numeric study ID

trtn

numeric treatment code

trtc

treatment name

latitude

absolute degrees latitude

r

number diagnosed with TB

n

sample size

References

Berkey CS, Hoaglin DC, Mosteller F, Colditz GA (1995). “A random-effects regression model for meta-analysis.” Statistics in Medicine, 14(4), 395–411. doi:10.1002/sim.4780140406.

Dias S, Sutton AJ, Welton NJ, Ades AE (2011). “NICE DSU Technical Support Document 3: Heterogeneity: subgroups, meta-regression, bias and bias-adjustment.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Description

Data frame containing the number of deaths in 22 trials comparing beta blockers vs. control for preventing mortality after myocardial infarction (Carlin 1992; Dias et al. 2011).

Usage

blocker

Format

A data frame with 44 rows and 5 variables:

studyn

numeric study ID

trtn

numeric treatment code

trtc

treatment name

r

total number of events

n

total number of individuals

References

Carlin JB (1992). “Meta-analysis for 2 x 2 tables: A bayesian approach.” Statistics in Medicine, 11(2), 141–158. doi:10.1002/sim.4780110202.

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Description

Multiple data sources created using set_ipd(), set_agd_arm(), or set_agd_contrast() can be combined into a single network for analysis.

Usage

combine_network(..., trt_ref)

Arguments

Value

An object of class nma_data

See Also

set_ipd(), set_agd_arm(), and set_agd_contrast() for defining different data sources.

print.nma_data() for the print method displaying details of the network, and plot.nma_data() for network plots.

Examples

## Parkinson's - combining contrast- and arm-based data
studies <- parkinsons$studyn
(parkinsons_arm <- parkinsons[studies %in% 1:3, ])
(parkinsons_contr <- parkinsons[studies %in% 4:7, ])

park_arm_net <- set_agd_arm(parkinsons_arm,
                            study = studyn,
                            trt = trtn,
                            y = y,
                            se = se,
                            sample_size = n)

park_contr_net <- set_agd_contrast(parkinsons_contr,
                                   study = studyn,
                                   trt = trtn,
                                   y = diff,
                                   se = se_diff,
                                   sample_size = n)

park_net <- combine_network(park_arm_net, park_contr_net)

# Print network details
park_net

# Plot network
plot(park_net, weight_edges = TRUE, weight_nodes = TRUE)

## Plaque Psoriasis - combining IPD and AgD in a network
# Set up plaque psoriasis network combining IPD and AgD
library(dplyr)
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)

pso_ipd <- pso_ipd %>%
  mutate(# Variable transformations
    bsa = bsa / 100,
    prevsys = as.numeric(prevsys),
    psa = as.numeric(psa),
    weight = weight / 10,
    durnpso = durnpso / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker"),
    # Check complete cases for covariates of interest
    complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>%
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker")
  )

# Exclude small number of individuals with missing covariates
pso_ipd <- filter(pso_ipd, complete)

pso_net <- combine_network(
  set_ipd(pso_ipd,
          study = studyc,
          trt = trtc,
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd,
              study = studyc,
              trt = trtc,
              r = pasi75_r,
              n = pasi75_n,
              trt_class = trtclass)
)

# Print network details
pso_net


# Plot network
plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE)

Description

Density, distribution, and quantile function for the generalised t distribution with degrees of freedom df, shifted by location and scaled by scale.

Usage

dgent(x, df, location = 0, scale = 1)

pgent(q, df, location = 0, scale = 1)

qgent(p, df, location = 0, scale = 1)

Arguments

Value

dgent() gives the density, pgent() gives the distribution function, qgent() gives the quantile function.

Description

Data frame containing the number of new cases of diabetes in 22 trials of 6 antihypertensive drugs (Elliott and Meyer 2007; Dias et al. 2011). The trial duration (in years) is also recorded.

Usage

diabetes

Format

A data frame with 48 rows and 7 variables:

studyn

numeric study ID

studyc

study name

trtn

numeric treatment code

trtc

treatment name

r

total number of events

n

total number of individuals

time

trial follow-up (years)

References

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Elliott WJ, Meyer PM (2007). “Incident diabetes in clinical trials of antihypertensive drugs: a network meta-analysis.” The Lancet, 369(9557), 201–207. doi:10.1016/s0140-6736(07)60108-1.

Description

Calculate the DIC for a model fitted using the nma() function.

Usage

dic(x, penalty = c("pD", "pV"), ...)

Arguments

Value

A nma_dic object.

See Also

print.nma_dic() for printing details, plot.nma_dic() for producing plots of residual deviance contributions.

Examples

## Smoking cessation

# Run smoking FE NMA example if not already available
if (!exists("smk_fit_FE")) example("example_smk_fe", run.donttest = TRUE)


# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)


# Compare DIC of FE and RE models
(smk_dic_FE <- dic(smk_fit_FE))
(smk_dic_RE <- dic(smk_fit_RE))   # substantially better fit

# Plot residual deviance contributions under RE model
plot(smk_dic_RE)

# Check for inconsistency using UME model


# Run smoking UME NMA example if not already available
if (!exists("smk_fit_RE_UME")) example("example_smk_ume", run.donttest = TRUE)


# Compare DIC
smk_dic_RE
(smk_dic_RE_UME <- dic(smk_fit_RE_UME))  # no difference in fit

# Compare residual deviance contributions
plot(smk_dic_RE, smk_dic_RE_UME, show_uncertainty = FALSE)

Description

Data frame containing the number of deaths and person-years at risk in 10 trials comparing reduced fat diets vs. control (non-reduced fat diet) for preventing mortality (Hooper et al. 2000; Dias et al. 2011).

Usage

dietary_fat

Format

A data frame with 21 rows and 7 variables:

studyn

numeric study ID

studyc

study name

trtn

numeric treatment code

trtc

treatment name

r

number of events

n

number randomised

E

person-years at risk

References

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Hooper L, Summerbell CD, Higgins JPT, Thompson RL, Clements G, Capps N, Davey Smith G, Riemersma R, Ebrahim S (2000). “Reduced or modified dietary fat for preventing cardiovascular disease.” Cochrane Database of Systematic Reviews. ISSN 1465-1858, doi:10.1002/14651858.CD002137.

Description

distr() is used within the function add_integration() to specify marginal distributions for the covariates, via a corresponding inverse CDF. It is also used in predict.stan_nma() to specify a distribution for the baseline response (intercept) when predicting absolute outcomes.

Usage

distr(qfun, ...)

Arguments

Details

The function qfun should have a formal argument called p. This restriction serves as a crude check for inverse CDFs (e.g. an error will be given if dnorm is used instead of qnorm). If a user-written CDF is supplied, it must have an argument p which takes a vector of probabilities.

Value

An object of class distr.

See Also

add_integration() where distr() is used to specify marginal distributions for covariates to integrate over, and predict.stan_nma() where distr() is used to specify a distribution on the baseline response.

Examples

## Specifying marginal distributions for integration

df <- data.frame(x1_mean = 2, x1_sd = 0.5, x2 = 0.8)

# Distribution parameters are evaluated in the context of the data frame
add_integration(df,
                x1 = distr(qnorm, mean = x1_mean, sd = x1_sd),
                x2 = distr(qbern, prob = x2),
                cor = diag(2))

Description

Density, distribution, and quantile function for the log t distribution, whose logarithm has degrees of freedom df, mean location, and standard deviation scale.

Usage

dlogt(x, df, location = 0, scale = 1)

plogt(q, df, location = 0, scale = 1)

qlogt(p, df, location = 0, scale = 1)

Arguments

Details

If $\log(Y) \sim t_\nu(\mu, \sigma^2)$, then $Y$ has a log t distribution with location $\mu$, scale $\sigma$, and df $\nu$.

The mean and all higher moments of the log t distribution are undefined or infinite.

If df = 1 then the distribution is a log Cauchy distribution. As df tends to infinity, this approaches a log Normal distribution.

Value

dlogt() gives the density, plogt() gives the distribution function, qlogt() gives the quantile function.

Description

Density, distribution, quantile, hazard, cumulative hazard, and restricted mean survival time functions for the M-spline baseline hazards model.

Usage

dmspline(x, basis, scoef, rate, log = FALSE)

pmspline(q, basis, scoef, rate, lower.tail = TRUE, log.p = FALSE)

qmspline(p, basis, scoef, rate, lower.tail = TRUE, log.p = FALSE)

hmspline(x, basis, scoef, rate, log = FALSE)

Hmspline(x, basis, scoef, rate, log = FALSE)

rmst_mspline(t, basis, scoef, rate, start = 0)

Arguments

Details

Survival models with a flexible M-spline on the baseline hazard are described by Brilleman et al. (2020). Piecewise-exponential baseline hazards are a special case where the degree of the M-spline polynomial is 0.

The d/p/h/H functions are calculated from their definitions. qmspline() uses numerical inversion via flexsurv::qgeneric(). rmst_mspline()uses numerical integration via flexsurv::rmst_generic(), except for the special case of the piecewise-exponential hazard (i.e. degree 0 M-splines) which uses the explicit formula from Royston and Parmar (2013).

Beyond the boundary knots, the hazard is assumed to be constant. (This differs from the approach in splines2::mSpline() that extrapolates the polynomial basis functions, which is numerically unstable and highly dependent on the data just before the boundary knots.) As with all extrapolation, care should be taken when evaluating the splines at times beyond the boundary knots (either directly through the d/p/h/H/rmst functions, or indirectly by requesting quantiles with qmspline() that correspond to times beyond the boundary knots). For this reason evaluating the (unrestricted) mean survival time is not generally recommended as this requires integrating over an infinite time horizon (i.e. rmst_mspline() with t = Inf).

Value

dmspline() gives the density, pmspline() gives the distribution function (CDF), qmspline() gives the quantile function (inverse-CDF), hmspline() gives the hazard function, Hmspline() gives the cumulative hazard function, and rmst_mspline() gives restricted mean survival times.

References

Brilleman SL, Elci EM, Novik JB, Wolfe R (2020). “Bayesian Survival Analysis Using the rstanarm R Package.” arXiv. doi:10.48550/arXiv.2002.09633, 2002.09633.

Royston P, Parmar MKB (2013). “Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome.” BMC Medical Research Methodology, 13(1). doi:10.1186/1471-2288-13-152.

Description

Calling example("example_ndmm") will run a proportional hazards Weibull NMA model on the newly-diagnosed multiple myeloma data, using the code in the Examples section below. The resulting stan_nma object ndmm_fit will then be available in the global environment.

Examples

# Set up newly-diagnosed multiple myeloma network

head(ndmm_ipd)
head(ndmm_agd)

ndmm_net <- combine_network(
  set_ipd(ndmm_ipd,
          study, trt,
          Surv = Surv(eventtime / 12, status)),
  set_agd_surv(ndmm_agd,
               study, trt,
               Surv = Surv(eventtime / 12, status),
               covariates = ndmm_agd_covs))

# Fit Weibull (PH) model
ndmm_fit <- nma(ndmm_net, 
                likelihood = "weibull",
                prior_intercept = normal(scale = 100),
                prior_trt = normal(scale = 10),
                prior_aux = half_normal(scale = 10))

ndmm_fit

Description

Calling example("example_pso_mlnmr") will run a ML-NMR model with the plaque psoriasis IPD and AgD, using the code in the Examples section below. The resulting stan_nma object pso_fit will then be available in the global environment.

Details

Plaque psoriasis ML-NMR for use in examples.

Examples

# Set up plaque psoriasis network combining IPD and AgD
library(dplyr)
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)

pso_ipd <- pso_ipd %>%
  mutate(# Variable transformations
    bsa = bsa / 100,
    prevsys = as.numeric(prevsys),
    psa = as.numeric(psa),
    weight = weight / 10,
    durnpso = durnpso / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker"),
    # Check complete cases for covariates of interest
    complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>%
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker")
  )

# Exclude small number of individuals with missing covariates
pso_ipd <- filter(pso_ipd, complete)

pso_net <- combine_network(
  set_ipd(pso_ipd,
          study = studyc,
          trt = trtc,
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd,
              study = studyc,
              trt = trtc,
              r = pasi75_r,
              n = pasi75_n,
              trt_class = trtclass)
)

# Print network details
pso_net

# Add integration points to the network
pso_net <- add_integration(pso_net,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  n_int = 64)


# Fitting a ML-NMR model.
# Specify a regression model to include effect modifier interactions for five
# covariates, along with main (prognostic) effects. We use a probit link and
# specify that the two-parameter Binomial approximation for the aggregate-level
# likelihood should be used. We set treatment-covariate interactions to be equal
# within each class. We narrow the possible range for random initial values with
# init_r = 0.1, since probit models in particular are often hard to initialise.
# Using the QR decomposition greatly improves sampling efficiency here, as is
# often the case for regression models.
pso_fit <- nma(pso_net, 
               trt_effects = "fixed",
               link = "probit",
               likelihood = "bernoulli2",
               regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
               class_interactions = "common",
               prior_intercept = normal(scale = 10),
               prior_trt = normal(scale = 10),
               prior_reg = normal(scale = 10),
               init_r = 0.1,
               QR = TRUE)
pso_fit

Description

Calling example("example_smk_fe") will run a fixed effects NMA model with the smoking cessation data, using the code in the Examples section below. The resulting stan_nma object smk_fit_FE will then be available in the global environment.

Details

Smoking FE NMA for use in examples.

Examples

# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net


# Fitting a fixed effect model
smk_fit_FE <- nma(smk_net, 
                  trt_effects = "fixed",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100))

smk_fit_FE

Description

Calling example("example_smk_nodesplit") will run node-splitting models with the smoking cessation data, using the code in the Examples section below. The resulting nma_nodesplit_df object smk_fit_RE_nodesplit will then be available in the global environment.

Details

Smoking node-splitting for use in examples.

Examples

# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net


# Fitting all possible node-splitting models
smk_fit_RE_nodesplit <- nma(smk_net, 
                            consistency = "nodesplit",
                            trt_effects = "random",
                            prior_intercept = normal(scale = 100),
                            prior_trt = normal(scale = 100),
                            prior_het = normal(scale = 5))

Description

Calling example("example_smk_re") will run a random effects NMA model with the smoking cessation data, using the code in the Examples section below. The resulting stan_nma object smk_fit_RE will then be available in the global environment.

Details

Smoking RE NMA for use in examples.

Examples

# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net


# Fitting a random effects model
smk_fit_RE <- nma(smk_net, 
                  trt_effects = "random",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100),
                  prior_het = normal(scale = 5))

smk_fit_RE

Description

Calling example("example_smk_ume") will run an unrelated mean effects (inconsistency) NMA model with the smoking cessation data, using the code in the Examples section below. The resulting stan_nma object smk_fit_RE_UME will then be available in the global environment.

Details

Smoking UME NMA for use in examples.

Examples

# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net


# Fitting an unrelated mean effects (inconsistency) model
smk_fit_RE_UME <- nma(smk_net, 
                      consistency = "ume",
                      trt_effects = "random",
                      prior_intercept = normal(scale = 100),
                      prior_trt = normal(scale = 100),
                      prior_het = normal(scale = 5))

smk_fit_RE_UME

Description

This helper function constructs a ggplot2 geom to plot Kaplan-Meier curves from a network containing survival or time-to-event outcomes. This is useful for overlaying the "raw" survival data on the estimated survival functions created with plotted with plot.surv_nma_summary(), but can also be used standalone to plot Kaplan-Meier curves before fitting a model.

Usage

geom_km(
  network,
  ...,
  transform = c("identity", "cloglog", "log", "cumhaz"),
  curve_args = list(),
  cens_args = list()
)

Arguments

Value

A ggplot2 geom list that can be added to a ggplot2 plot object

Examples

# Set up newly-diagnosed multiple myeloma network

head(ndmm_ipd)
head(ndmm_agd)

ndmm_net <- combine_network(
  set_ipd(ndmm_ipd,
          study, trt,
          Surv = Surv(eventtime / 12, status)),
  set_agd_surv(ndmm_agd,
               study, trt,
               Surv = Surv(eventtime / 12, status),
               covariates = ndmm_agd_covs))
# Plot KM curves using ggplot2
library(ggplot2)

# We need to create an empty ggplot object to add the curves to
ggplot() + geom_km(ndmm_net)

# Adding plotting options, facets, axis labels, and a plot theme
ggplot() +
  geom_km(ndmm_net,
          curve_args = list(linewidth = 0.5),
          cens_args = list(size = 3, shape = 124)) +
  facet_wrap(vars(Study)) +
  labs(xlab = "Time", ylab = "Survival Probability") +
  theme_multinma()

# Using the transform argument to produce log-log plots (e.g. to assess the
# proportional hazards assumption)
ggplot() +
  geom_km(ndmm_net, transform = "cloglog") +
  facet_wrap(vars(Study)) +
  theme_multinma()

# Using the transform argument to produce cumulative hazard plots
ggplot() +
  geom_km(ndmm_net, transform = "cumhaz") +
  facet_wrap(vars(Study)) +
  theme_multinma()

# This function can also be used to add KM data to plots of estimated survival
# curves from a fitted model, in a similar manner

# Run newly-diagnosed multiple myeloma example if not already available
if (!exists("ndmm_fit")) example("example_ndmm", run.donttest = TRUE)

# Plot estimated survival curves, and overlay the KM data

plot(predict(ndmm_fit, type = "survival")) + geom_km(ndmm_net)

Description

Determine whether two treatments in a network are connected by direct and/or indirect evidence, and generate a list of comparisons with both direct and indirect evidence (i.e. potential inconsistency) for node-splitting.

Usage

get_nodesplits(network, include_consistency = FALSE)

has_direct(network, trt1, trt2)

has_indirect(network, trt1, trt2)

Arguments

Details

The list of comparisons for node-splitting is generated following the algorithm of van Valkenhoef et al. (2016). A comparison between two treatments has the potential for inconsistency, and is thus considered for node-splitting, if the comparison has both direct evidence and independent indirect evidence.

The notion of independent indirect evidence is necessary when multi-arm trials are present, since by design these trials are internally consistent. A comparison between two treatments has independent indirect evidence if, after removing all studies comparing the two treatments from the network, the two treatments are still connected by a path of evidence. This is the criterion considered by the has_indirect() function.

Value

For has_direct() and has_indirect(), a single logical value. For get_nodesplits(), a data frame with two columns giving the comparisons for node-splitting.

References

van Valkenhoef G, Dias S, Ades AE, Welton NJ (2016). “Automated generation of node-splitting models for assessment of inconsistency in network meta-analysis.” Research Synthesis Methods, 7(1), 80–93. doi:10.1002/jrsm.1167.

Examples

# Parkinsons example
park_net <- set_agd_arm(parkinsons,
                        study = studyn,
                        trt = trtn,
                        y = y,
                        se = se,
                        trt_ref = 1)

# View the network plot
plot(park_net)

# The 4 vs. 5 comparison is a spur on the network
has_direct(park_net, 4, 5)
has_indirect(park_net, 4, 5)

# 1 and 5 are not directly connected
has_direct(park_net, 1, 5)
has_indirect(park_net, 1, 5)

# The 1 vs. 2 comparison does not have independent indirect evidence, since
# the 1-2-4 loop is a multi-arm study
has_indirect(park_net, 1, 2)

# Get a list of comparisons with potential inconsistency for node-splitting
get_nodesplits(park_net)

# See van Valkenhoef (2016) for a discussion of this example

Description

Data frame containing the results of 16 trials comparing 8 treatments for moderate-to-severe plaque psoriasis from an HTA report (Woolacott et al. 2006), analysed in TSD2 (Dias et al. 2011). Outcomes are success/failure to achieve 50%, 75%, or 90% reduction in symptoms on the Psoriasis Area and Severity Index (PASI) scale. Some studies report all three ordered outcomes, others only one or two. The latter are coded as missing values (see details).

Usage

hta_psoriasis

Format

A data frame with 36 rows and 9 variables:

studyn

numeric study ID

studyc

study name

year

year of publication

trtn

numeric treatment code

trtc

treatment name

sample_size

sample size in each arm

PASI50, PASI75, PASI90

ordered multinomial outcome counts (exclusive, see details)

Details

Outcome counts are "exclusive"; that is, for a study reporting all outcomes, the counts represent the categories 50 < PASI < 75, 75 < PASI < 90, and 90 < PASI < 100, and are named by the lower end of the interval. (As opposed to "inclusive" counts, which would represent the overlapping categories PASI > 50, PASI > 70, and PASI > 90.) The count for the fourth category (the lowest), 0 < PASI < 50, is equal to sample_size - PASI50 - PASI75 - PASI90.

Missing values are used where studies only report a subset of the outcomes. For a study reporting only two outcomes, say 50 and 75, the counts represent 50 < PASI < 75 and 75 < PASI < 100. For a study reporting only one outcome, say PASI 75, the count represents 75 < PASI < 100.

References

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Woolacott N, Hawkins N, Mason A, Kainth A, Khadjesari Z, Bravo Vergel Y, Misso K, Light K, Chalmers R, Sculpher M, Riemsma R (2006). “Etanercept and efalizumab for the treatment of psoriasis: a systematic review.” Health Technology Assessment, 10(46). doi:10.3310/hta10460.

Description

Check whether a network is connected - whether there is a path of study evidence linking every pair of treatments in the network.

Usage

is_network_connected(network)

Arguments

Details

Models will still run with disconnected networks. However, estimated relative effects between treatments across disconnected parts of the network will be entirely based on the prior distribution (typically very uncertain), as there is no information to update the prior distribution. Relative effects within each connected sub-network will be estimated as if each sub-network had been analysed separately.

Value

Logical TRUE or FALSE

Examples

## Smoking cessation
# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net

is_network_connected(smk_net)  # TRUE, network is connected
## A disconnected network
disc_net <- set_agd_arm(smoking[smoking$studyn %in% c(15, 21), ],
                        study = studyn,
                        trt = trtc,
                        r = r,
                        n = n)
is_network_connected(disc_net)  # FALSE, network is disconnected
disc_net
plot(disc_net)

Description

The loo() and waic() functions from the loo package may be called directly on stan_nma and stan_mlnmr objects.

Usage

## S3 method for class 'stan_nma'
loo(x, ...)

## S3 method for class 'stan_nma'
waic(x, ...)

Arguments

Description

Several different algorithms are provided to calculate knot locations for M-spline baseline hazard models. This function is called internally within the nma() function, but may be called directly by the user for more control.

Usage

make_knots(
  network,
  n_knots = 7,
  type = c("quantile", "quantile_common", "quantile_lumped", "quantile_longest", "equal",
    "equal_common")
)

Arguments

Details

The type argument can be used to choose between different algorithms for placing the knots:

"quantile"

Creates separate knot locations for each study, internal knots are placed at evenly-spaced quantiles of the observed event times within each study.

"quantile_lumped"

Creates a common set of knots for all studies, calculated as evenly-spaced quantiles of the observed event times from all studies lumped together.

"quantile_common"

Creates a common set of knots for all studies, taking quantiles of the quantiles of the observed event times within each study. This often seems to result in a more even knot spacing than "quantile_lumped", particularly when follow-up is uneven across studies, and may handle differing behaviour in the baseline hazard across studies better than "quantile_longest".

"quantile_longest"

Creates a common set of knots for all studies, using evenly-spaced quantiles of the observed event times in the longest study.

"equal"

Creates separate knot locations for each study, at evenly-spaced times between the boundary knots in each study.

"equal_common"

Creates a common set of knots for all studies, at evenly-spaced times between the earliest entry time and last event/censoring time in the network.

Boundary knots are calculated as follows:

• For separate knot locations in each study, boundary knots are placed at the earliest entry time and last event/censoring time in each study.

• For a common set of knots across all studies, boundary knots are placed at the earliest entry time and last event/censoring time across all studies.

Models with regression on the spline coefficients (i.e. with aux_regression specified) require a common set of knots across all studies.

Provided that a sufficient number of knots are used, model fit should be largely unaffected by the knot locations. However, sampling difficulties can sometimes occur if knot placement is poor, for example if a knot is placed just before the last follow-up time in a study.

Value

A named list of vectors giving the knot locations in each study.

Examples

# Set up newly-diagnosed multiple myeloma network

head(ndmm_ipd)
head(ndmm_agd)

ndmm_net <- combine_network(
  set_ipd(ndmm_ipd,
          study, trt,
          Surv = Surv(eventtime / 12, status)),
  set_agd_surv(ndmm_agd,
               study, trt,
               Surv = Surv(eventtime / 12, status),
               covariates = ndmm_agd_covs))

# The default knot locations
make_knots(ndmm_net, type = "quantile")

# Increasing the number of knots
make_knots(ndmm_net, n_knots = 10)

# Comparing alternative knot positioning algorithms
# Visualise these with a quick function
plot_knots <- function(network, knots) {
  ggplot2::ggplot() +
    geom_km(network) +
    ggplot2::geom_vline(ggplot2::aes(xintercept = .data$knot),
                        data = tidyr::pivot_longer(as.data.frame(knots), cols = dplyr::everything(),
                                                   names_to = "Study", values_to = "knot"),
                        linetype = 2, colour = "grey60") +
    ggplot2::facet_wrap(~Study) +
    theme_multinma()
}

plot_knots(ndmm_net, make_knots(ndmm_net, type = "quantile"))
plot_knots(ndmm_net, make_knots(ndmm_net, type = "quantile_common"))
plot_knots(ndmm_net, make_knots(ndmm_net, type = "quantile_lumped"))
plot_knots(ndmm_net, make_knots(ndmm_net, type = "quantile_longest"))
plot_knots(ndmm_net, make_knots(ndmm_net, type = "equal"))
plot_knots(ndmm_net, make_knots(ndmm_net, type = "equal_common"))

Description

Generate population-average marginal treatment effects. These are formed from population-average absolute predictions, so this function is a wrapper around predict.stan_nma().

Usage

marginal_effects(
  object,
  ...,
  mtype = c("difference", "ratio", "link"),
  all_contrasts = FALSE,
  trt_ref = NULL,
  probs = c(0.025, 0.25, 0.5, 0.75, 0.975),
  predictive_distribution = FALSE,
  summary = TRUE
)

Arguments

Value

A nma_summary object if summary = TRUE, otherwise a list containing a 3D MCMC array of samples and (for regression models) a data frame of study information.

Examples

## Smoking cessation

# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)


# Marginal risk difference in each study population in the network
marginal_effects(smk_fit_RE, mtype = "difference")

# Since there are no covariates in the model, the marginal and conditional
# (log) odds ratios here coincide
marginal_effects(smk_fit_RE, mtype = "link")
relative_effects(smk_fit_RE)

# Marginal risk differences in a population with 67 observed events out of
# 566 individuals on No Intervention, corresponding to a Beta(67, 566 - 67)
# distribution on the baseline probability of response
(smk_rd_RE <- marginal_effects(smk_fit_RE,
                               baseline = distr(qbeta, 67, 566 - 67),
                               baseline_type = "response",
                               mtype = "difference"))
plot(smk_rd_RE)


## Plaque psoriasis ML-NMR

# Run plaque psoriasis ML-NMR example if not already available
if (!exists("pso_fit")) example("example_pso_mlnmr", run.donttest = TRUE)


# Population-average marginal probit differences in each study in the network
(pso_marg <- marginal_effects(pso_fit, mtype = "link"))
plot(pso_marg, ref_line = c(0, 1))

# Population-average marginal probit differences in a new target population,
# with means and SDs or proportions given by
new_agd_int <- data.frame(
  bsa_mean = 0.6,
  bsa_sd = 0.3,
  prevsys = 0.1,
  psa = 0.2,
  weight_mean = 10,
  weight_sd = 1,
  durnpso_mean = 3,
  durnpso_sd = 1
)

# We need to add integration points to this data frame of new data
# We use the weighted mean correlation matrix computed from the IPD studies
new_agd_int <- add_integration(new_agd_int,
                               durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
                               prevsys = distr(qbern, prob = prevsys),
                               bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
                               weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
                               psa = distr(qbern, prob = psa),
                               cor = pso_net$int_cor,
                               n_int = 64)

# Population-average marginal probit differences of achieving PASI 75 in this
# target population, given a Normal(-1.75, 0.08^2) distribution on the
# baseline probit-probability of response on Placebo (at the reference levels
# of the covariates), are given by
(pso_marg_new <- marginal_effects(pso_fit,
                                  mtype = "link",
                                  newdata = new_agd_int,
                                  baseline = distr(qnorm, -1.75, 0.08)))
plot(pso_marg_new)


## Progression free survival with newly-diagnosed multiple myeloma

# Run newly-diagnosed multiple myeloma example if not already available
if (!exists("ndmm_fit")) example("example_ndmm", run.donttest = TRUE)


# We can produce a range of marginal effects from models with survival
# outcomes, specified with the mtype and type arguments. For example:

# Marginal survival probability difference at 5 years, all contrasts
marginal_effects(ndmm_fit, type = "survival", mtype = "difference",
                 times = 5, all_contrasts = TRUE)

# Marginal difference in RMST up to 5 years
marginal_effects(ndmm_fit, type = "rmst", mtype = "difference", times = 5)

# Marginal median survival time ratios
marginal_effects(ndmm_fit, type = "median", mtype = "ratio")

# Marginal log hazard ratios
# With no covariates in the model, these are constant over time and study
# populations, and are equal to the log hazard ratios from relative_effects()
plot(marginal_effects(ndmm_fit, type = "hazard", mtype = "link"),
     # The hazard is infinite at t=0 in some studies, giving undefined logHRs at t=0
     na.rm = TRUE)

# The NDMM vignette demonstrates the production of time-varying marginal
# hazard ratios from a ML-NMR model that includes covariates, see
# `vignette("example_ndmm")`

# Marginal survival difference over time
plot(marginal_effects(ndmm_fit, type = "survival", mtype = "difference"))

Description

3D MCMC arrays (Iterations, Chains, Parameters) are produced by as.array() methods applied to stan_nma or nma_summary objects.

Usage

## S3 method for class 'mcmc_array'
summary(object, ..., probs = c(0.025, 0.25, 0.5, 0.75, 0.975))

## S3 method for class 'mcmc_array'
print(x, ...)

## S3 method for class 'mcmc_array'
plot(x, ...)

## S3 method for class 'mcmc_array'
names(x)

## S3 replacement method for class 'mcmc_array'
names(x) <- value

Arguments

Value

The summary() method returns a nma_summary object, the print() method returns x invisibly. The names() method returns a character vector of parameter names, and ⁠names()<-⁠ returns the object with updated parameter names. The plot() method is a shortcut for plot(summary(x), ...), passing all arguments on to plot.nma_summary().

Examples

## Smoking cessation

# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)


# Working with arrays of posterior draws (as mcmc_array objects) is
# convenient when transforming parameters

# Transforming log odds ratios to odds ratios
LOR_array <- as.array(relative_effects(smk_fit_RE))
OR_array <- exp(LOR_array)

# mcmc_array objects can be summarised to produce a nma_summary object
smk_OR_RE <- summary(OR_array)

# This can then be printed or plotted
smk_OR_RE
plot(smk_OR_RE, ref_line = 1)

# Transforming heterogeneity SD to variance
tau_array <- as.array(smk_fit_RE, pars = "tau")
tausq_array <- tau_array^2

# Correct parameter names
names(tausq_array) <- "tausq"

# Summarise
summary(tausq_array)

Description

This function aids the specification of multinomial outcome data when setting up a network with set_agd_arm() or set_ipd(). It takes a set of columns (or, more generally, numeric vectors of the same length) of outcome counts in each category, and binds these together to produce a matrix.

Usage

multi(..., inclusive = FALSE, type = c("ordered", "competing"))

Arguments

Details

When specifying ordered categorical counts, these can either be given as exclusive counts (inclusive = FALSE, the default) where individuals are only counted in the highest category they achieve, or inclusive counts (inclusive = TRUE) where individuals are counted in every category up to and including the highest category achieved. (Competing outcomes, by nature, are always specified as exclusive counts.)

NA values can be used to indicate categories/cutpoints that were not measured.

Value

A matrix of (exclusive) category counts

Examples

# These two data sets specify the same ordered categorical data for outcomes
# r0 < r1 < r2, but the first uses the "inclusive" format and the second the
# "exclusive" format.
df_inclusive <- tibble::tribble(~r0, ~r1, ~r2,
                                1, 1, 1,
                                5, 4, 1,
                                5, 2, 2,
                                10, 5, 0,
                                5, 5, 0,
                                7, NA, 6,   # Achieved r2 or not (no r1)
                                10, 4, NA)  # Achieved r1 or not (no r2)

df_exclusive <- tibble::tribble(~r0, ~r1, ~r2,
                                0, 0, 1,
                                1, 3, 1,
                                3, 0, 2,
                                5, 5, 0,
                                0, 5, 0,
                                1, NA, 6,   # Achieved r2 or not (no r1)
                                6, 4, NA)   # Achieved r1 or not (no r2)

(r_inclusive <- with(df_inclusive, multi(r0, r1, r2, inclusive = TRUE)))
(r_exclusive <- with(df_exclusive, multi(r0, r1, r2, inclusive = FALSE)))

# Counts are always stored in exclusive format
stopifnot(isTRUE(all.equal(r_inclusive, r_exclusive)))


## HTA Plaque Psoriasis
library(dplyr)

# Ordered outcomes here are given as "exclusive" counts
head(hta_psoriasis)

# Calculate lowest category count (failure to achieve PASI 50)
pso_dat <- hta_psoriasis %>%
  mutate(`PASI<50` = sample_size - rowSums(cbind(PASI50, PASI75, PASI90), na.rm = TRUE))

# Set up network
pso_net <- set_agd_arm(pso_dat,
                       study = paste(studyc, year),
                       trt = trtc,
                       r = multi(`PASI<50`, PASI50, PASI75, PASI90,
                                 inclusive = FALSE,
                                 type = "ordered"))

pso_net

Description

Three data frames, ndmm_ipd, ndmm_agd, and ndmm_agd_covs containing (simulated) individual patient data (IPD) from three studies and aggregate data (AgD) from two studies on newly diagnosed multiple myeloma. The outcome of interest is progression-free survival after autologous stem cell transplant. The IPD studies in ndmm_ipd provide event/censoring times and covariate values for each individual. The AgD studies provide reconstructed event/censoring times from digitized Kaplan-Meier curves in ndmm_agd and covariate summaries in ndmm_agd_covs, obtained from published trial reports. The data are constructed to resemble those used by Leahy and Walsh (2019).

Usage

ndmm_ipd

ndmm_agd

ndmm_agd_covs

Format

The individual patient data are contained in a data frame ndmm_ipd with 1325 rows, one per individual, and 10 variables:

study, studyf

study name

trt, trtf

treatment name

eventtime

event/censoring time

status

censoring indicator (0 = censored, 1 = event)

age

age (years)

iss_stage3

ISS stage 3 (0 = no, 1 = yes)

response_cr_vgpr

complete or very good partial response (0 = no, 1 = yes)

male

male sex (0 = no, 1 = yes)

The reconstructed Kaplan-Meier data for the aggregate studies are contained in a data frame ndmm_agd with 2819 rows and 6 variables:

study, studyf

study name

trt, trtf

treatment name

eventtime

event/censoring time

status

censoring indicator (0 = censored, 1 = event)

The covariate summaries extracted from published reportes for the aggregate studies are contained in a data frame ndmm_agd_covs with 4 rows, one per study arm, and 15 columns:

study, studyf

study name

trt, trtf

treatment name

sample_size

sample size in each arm

age_min, age_iqr_l, age_median, age_iqr_h, age_max, age_mean, age_sd

summary statistics for age (years)

iss_stage3

proportion of participants with ISS stage 3

response_cr_vgpr

proportion of participants with complete or very good partial response

male

proportion of male participants

References

Leahy J, Walsh C (2019). “Assessing the impact of a matching-adjusted indirect comparison in a Bayesian network meta-analysis.” Research Synthesis Methods, 10(4), 546–568. doi:10.1002/jrsm.1372.

Description

The nma function fits network meta-analysis and (multilevel) network meta-regression models in Stan.

Usage

nma(
  network,
  consistency = c("consistency", "ume", "nodesplit"),
  trt_effects = c("fixed", "random"),
  regression = NULL,
  class_interactions = c("common", "exchangeable", "independent"),
  likelihood = NULL,
  link = NULL,
  ...,
  nodesplit = get_nodesplits(network, include_consistency = TRUE),
  prior_intercept = .default(normal(scale = 100)),
  prior_trt = .default(normal(scale = 10)),
  prior_het = .default(half_normal(scale = 5)),
  prior_het_type = c("sd", "var", "prec"),
  prior_reg = .default(normal(scale = 10)),
  prior_aux = .default(),
  prior_aux_reg = .default(),
  aux_by = NULL,
  aux_regression = NULL,
  QR = FALSE,
  center = TRUE,
  adapt_delta = NULL,
  int_thin = 0,
  int_check = TRUE,
  mspline_degree = 3,
  n_knots = 7,
  knots = NULL,
  mspline_basis = NULL
)

Arguments

Details

When specifying a model formula in the regression argument, the usual formula syntax is available (as interpreted by model.matrix()). The only additional requirement here is that the special variable .trt should be used to refer to treatment. For example, effect modifier interactions should be specified as variable:.trt. Prognostic (main) effects and interactions can be included together compactly as variable*.trt, which expands to variable + variable:.trt (plus .trt, which is already in the NMA model).

For the advanced user, the additional specials .study and .trtclass are also available, and refer to studies and (if specified) treatment classes respectively. When node-splitting models are fitted (consistency = "nodesplit") the special .omega is available, indicating the arms to which the node-splitting inconsistency factor is added.

See ?priors for details on prior specification. Default prior distributions are available, but may not be appropriate for the particular setting and will raise a warning if used. No attempt is made to tailor these defaults to the data provided. Please consider appropriate prior distributions for the particular setting, accounting for the scales of outcomes and covariates, etc. The function plot_prior_posterior() may be useful in examining the influence of the chosen prior distributions on the posterior distributions, and the summary() method for nma_prior objects prints prior intervals.

Value

nma() returns a stan_nma object, except when consistency = "nodesplit" when a nma_nodesplit or nma_nodesplit_df object is returned. nma.fit() returns a stanfit object.

Likelihoods and link functions

Currently, the following likelihoods and link functions are supported for each data type:

The bernoulli2 and binomial2 likelihoods correspond to a two-parameter Binomial likelihood for arm-based AgD, which more closely matches the underlying Poisson Binomial distribution for the summarised aggregate outcomes in a ML-NMR model than the typical (one parameter) Binomial distribution (see Phillippo et al. 2020).

When a cloglog link is used, including an offset for log follow-up time (i.e. regression = ~offset(log(time))) results in a model on the log hazard (see Dias et al. 2011).

For survival data, all accelerated failure time models (exponential-aft, weibull-aft, lognormal, loglogistic, gamma, gengamma) are parameterised so that the treatment effects and any regression parameters are log Survival Time Ratios (i.e. a coefficient of $\log(2)$ means that the treatment or covariate is associated with a doubling of expected survival time). These can be converted to log Acceleration Factors using the relation $\log(\mathrm{AF}) = -\log(\mathrm{STR})$ (or equivalently $\mathrm{AF} = 1/\mathrm{STR}$).

Further details on each likelihood and link function are given by Dias et al. (2011).

Auxiliary parameters

Auxiliary parameters are only present in the following models.

Normal likelihood with IPD

When a Normal likelihood is fitted to IPD, the auxiliary parameters are the arm-level standard deviations $\sigma_{jk}$ on treatment $k$ in study $j$.

Ordered multinomial likelihood

When fitting a model to $M$ ordered outcomes, the auxiliary parameters are the latent cutoffs between each category, $c_0 < c_1 < \dots < c_M$. Only $c_2$ to $c_{M-1}$ are estimated; we fix $c_0 = -\infty$, $c_1 = 0$, and $c_M = \infty$. When specifying priors for these latent cutoffs, we choose to specify priors on the differences $c_{m+1} - c_m$. Stan automatically truncates any priors so that the ordering constraints are satisfied.

Survival (time-to-event) likelihoods

All survival likelihoods except the exponential and exponential-aft likelihoods have auxiliary parameters. These are typically study-specific shape parameters $\gamma_j>0$, except for the lognormal likelihood where the auxiliary parameters are study-specific standard deviations on the log scale $\sigma_j>0$.

The gengamma likelihood has two sets of auxiliary parameters, study-specific scale parameters $\sigma_j>0$ and shape parameters $k_j$, following the parameterisation of Lawless (1980), which permits a range of behaviours for the baseline hazard including increasing, decreasing, bathtub and arc-shaped hazards. This parameterisation is related to that discussed by Cox et al. (2007) and implemented in the flexsurv package with $Q = k^{-0.5}$. The parameterisation used here effectively bounds the shape parameter $k$ away from numerical instabilities as $k \rightarrow \infty$ (i.e. away from $Q \rightarrow 0$, the log-Normal distribution) via the prior distribution. Implicitly, this parameterisation is restricted to $Q > 0$ and so certain survival distributions like the inverse-Gamma and inverse-Weibull are not part of the parameter space; however, $Q > 0$ still encompasses the other survival distributions implemented in this package.

For the mspline and pexp likelihoods, the auxiliary parameters are the spline coefficients for each study. These form a unit simplex (i.e. lie between 0 and 1, and sum to 1), and are given a random walk prior distribution. prior_aux specifies the hyperprior on the random walk standard deviation $\sigma$ which controls the level of smoothing of the baseline hazard, with $\sigma = 0$ corresponding to a constant baseline hazard.

The auxiliary parameters can be stratified by additional factors through the aux_by argument. For example, to allow the shape of the baseline hazard to vary between treatment arms as well as studies, use aux_by = c(".study", ".trt"). (Technically, .study is always included in the stratification even if omitted from aux_by, but we choose here to make the stratification explicit.) This is a common way of relaxing the proportional hazards assumption. The default is equivalent to aux_by = ".study" which stratifies the auxiliary parameters by study, as described above.

A regression model may be specified on the auxiliary parameters using aux_regression. This is useful if we wish to model departures from non-proportionality, rather than allowing the baseline hazards to be completely independent using aux_by. This is necessary if absolute predictions (e.g. survival curves) are required in a population for unobserved combinations of covariates; for example, if aux_by = .trt then absolute predictions may only be produced for the observed treatment arms in each study population, whereas if aux_regression = ~.trt then absolute predictions can be produced for all treatments in any population. For mspline and pexp likelihoods, the regression coefficients are smoothed over time using a random walk prior to avoid overfitting: prior_aux_reg specifies the hyperprior for the random walk standard deviation. For other parametric likelihoods, prior_aux_reg specifies the prior for the auxiliary regression coefficients.

References

Cox C, Chu H, Schneider MF, Muñoz A (2007). “Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution.” Statistics in Medicine, 26(23), 4352–4374. doi:10.1002/sim.2836.

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Lawless JF (1980). “Inference in the Generalized Gamma and Log Gamma Distributions.” Technometrics, 22(3), 409–419. doi:10.1080/00401706.1980.10486173.

Phillippo DM, Dias S, Ades AE, Belger M, Brnabic A, Schacht A, Saure D, Kadziola Z, Welton NJ (2020). “Multilevel Network Meta-Regression for population-adjusted treatment comparisons.” Journal of the Royal Statistical Society: Series A (Statistics in Society), 183(3), 1189–1210. doi:10.1111/rssa.12579.

Examples

## Smoking cessation NMA
# Set up network of smoking cessation data
head(smoking)

smk_net <- set_agd_arm(smoking,
                       study = studyn,
                       trt = trtc,
                       r = r,
                       n = n,
                       trt_ref = "No intervention")

# Print details
smk_net


# Fitting a fixed effect model
smk_fit_FE <- nma(smk_net, 
                  trt_effects = "fixed",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100))

smk_fit_FE



# Fitting a random effects model
smk_fit_RE <- nma(smk_net, 
                  trt_effects = "random",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100),
                  prior_het = normal(scale = 5))

smk_fit_RE



# Fitting an unrelated mean effects (inconsistency) model
smk_fit_RE_UME <- nma(smk_net, 
                      consistency = "ume",
                      trt_effects = "random",
                      prior_intercept = normal(scale = 100),
                      prior_trt = normal(scale = 100),
                      prior_het = normal(scale = 5))

smk_fit_RE_UME



# Fitting all possible node-splitting models
smk_fit_RE_nodesplit <- nma(smk_net, 
                            consistency = "nodesplit",
                            trt_effects = "random",
                            prior_intercept = normal(scale = 100),
                            prior_trt = normal(scale = 100),
                            prior_het = normal(scale = 5))



# Summarise the node-splitting results
summary(smk_fit_RE_nodesplit)


## Plaque psoriasis ML-NMR
# Set up plaque psoriasis network combining IPD and AgD
library(dplyr)
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)

pso_ipd <- pso_ipd %>%
  mutate(# Variable transformations
    bsa = bsa / 100,
    prevsys = as.numeric(prevsys),
    psa = as.numeric(psa),
    weight = weight / 10,
    durnpso = durnpso / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker"),
    # Check complete cases for covariates of interest
    complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>%
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker")
  )

# Exclude small number of individuals with missing covariates
pso_ipd <- filter(pso_ipd, complete)

pso_net <- combine_network(
  set_ipd(pso_ipd,
          study = studyc,
          trt = trtc,
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd,
              study = studyc,
              trt = trtc,
              r = pasi75_r,
              n = pasi75_n,
              trt_class = trtclass)
)

# Print network details
pso_net

# Add integration points to the network
pso_net <- add_integration(pso_net,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  n_int = 64)


# Fitting a ML-NMR model.
# Specify a regression model to include effect modifier interactions for five
# covariates, along with main (prognostic) effects. We use a probit link and
# specify that the two-parameter Binomial approximation for the aggregate-level
# likelihood should be used. We set treatment-covariate interactions to be equal
# within each class. We narrow the possible range for random initial values with
# init_r = 0.1, since probit models in particular are often hard to initialise.
# Using the QR decomposition greatly improves sampling efficiency here, as is
# often the case for regression models.
pso_fit <- nma(pso_net, 
               trt_effects = "fixed",
               link = "probit",
               likelihood = "bernoulli2",
               regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
               class_interactions = "common",
               prior_intercept = normal(scale = 10),
               prior_trt = normal(scale = 10),
               prior_reg = normal(scale = 10),
               init_r = 0.1,
               QR = TRUE)
pso_fit


## Newly-diagnosed multiple myeloma NMA
# Set up newly-diagnosed multiple myeloma network

head(ndmm_ipd)
head(ndmm_agd)

ndmm_net <- combine_network(
  set_ipd(ndmm_ipd,
          study, trt,
          Surv = Surv(eventtime / 12, status)),
  set_agd_surv(ndmm_agd,
               study, trt,
               Surv = Surv(eventtime / 12, status),
               covariates = ndmm_agd_covs))

# Fit Weibull (PH) model
ndmm_fit <- nma(ndmm_net, 
                likelihood = "weibull",
                prior_intercept = normal(scale = 100),
                prior_trt = normal(scale = 10),
                prior_aux = half_normal(scale = 10))

ndmm_fit

Description

The nma_data class contains the data for a NMA in a standard format, created using the functions set_ipd(), set_agd_arm(), set_agd_contrast(), set_agd_surv(), or combine_network(). The sub-class mlnmr_data is created by the function add_integration(), and further contains numerical integration points for the aggregate data.

Details

Objects of class nma_data have the following components:

agd_arm

data from studies with aggregate data (arm format)

agd_contrast

data from studies with aggregate data (contrast format)

ipd

data from studies with individual patient data

treatments

treatment coding factor for entire network

classes

treatment class coding factor (same length as treatments for entire network)

studies

study coding factor for entire network

outcome

outcome type for each data source, named list

The agd_arm, agd_contrast, and ipd components are tibbles with the following columns:

.study

study (as factor)

.trt

treatment (as factor)

.trtclass

treatment class (as factor), if specified

.y

continuous outcome

.se

standard error (continuous)

.r

event count (discrete)

.n

event count denominator (discrete, agd_arm only)

.E

time at risk (discrete)

.Surv

survival outcome of type Surv (time-to-event), nested by study arm

.sample_size

sample size (⁠agd_*⁠ only)

...

other columns (typically covariates) from the original data frame

Objects of class mlnmr_data additionally have components:

n_int

number of numerical integration points

int_names

names of covariates with numerical integration points

int_cor

correlation matrix for covariates used to generate numerical integration points

The agd_arm and agd_contrast tibbles have additional list columns with prefix .int_, one for each covariate, which contain the numerical integration points nested as length-n_int vectors within each row.

See Also

print.nma_data() for the print method displaying details of the network, and plot.nma_data() for network plots.

Description

The nma_dic class contains details of the Deviance Information Criterion (DIC), produced using the dic() function.

Details

Objects of class nma_dic have the following components:

dic

The DIC value

pd, pv

The effective number of parameters

resdev

The total residual deviance

pointwise

A list of data frames containing the pointwise contributions for the IPD and AgD.

resdev_array

A 3D MCMC array [Iterations, Chains, Parameters] of posterior residual deviance samples.

See Also

dic(), print.nma_dic(), plot.nma_dic().

Description

The nma_nodesplit and nma_nodesplit_df classes contains the results from running a node-splitting model with the function nma().

Details

Objects of class nma_nodesplit inherit from the stan_nma class, and contain the results of fitting a single node-split model. They have one additional component, nodesplit, which gives the comparison that was node-split as a length 2 vector.

Objects of class nma_nodesplit_df are tibble data frames with one row for each node-split comparison and columns:

trt1, trt2

Treatments forming the comparison

model

A list column containing the results of each model as a nma_nodesplit object

Optionally, there will be an additional row for the consistency model if this was fitted (e.g. by get_nodesplits(., include_consistency = TRUE)) with trt1 and trt2 both NA.

Description

The nma_prior class is used to specify prior distributions.

Details

Objects of class nma_prior have the following components:

dist

Distribution name

fun

Name of constructor function, as string (e.g. "normal")

...

Parameters of the distribution

The distribution parameters, specified as named components in ..., match those in the constructor functions (see priors).

Description

The nma_summary class contains posterior summary statistics of model parameters or other quantities of interest, and the draws used to obtain these statistics.

Details

Objects of class nma_summary have the following components:

summary

A data frame containing the computed summary statistics. Column .trt indicates the corresponding treatment, or columns .trta and .trtb indicate the corresponding contrast (.trtb vs. .trta). If a regression model was fitted with effect modifier interactions with treatment, these summaries will be study-specific. In this case, the corresponding study population is indicated in the .study column. If a multinomial model was fitted, the .category column indicates the corresponding category.

sims

A 3D array [Iteration, Chain, Parameter] of MCMC simulations

studies

(Optional) A data frame containing study information, printed along with the corresponding summary statistics if summary contains a .study column. Should have a matching .study column.

The following attributes may also be set:

xlab

Label for x axis in plots, usually either "Treatment" or "Contrast".

ylab

Label for y axis in plots, usually used for the scale e.g. "log Odds Ratio".

The subclass nma_rank_probs is used by the function posterior_rank_probs(), and contains posterior rank probabilities. This subclass does not have a sims component, as the rank probabilities are themselves posterior summaries of the ranks (i.e. they do not have a posterior distribution). The posterior ranks from which the rank probabilities are calculated may be obtained from posterior_ranks().

Description

The nodesplit_summary class contains posterior summary statistics for node-splitting models, as a result of calling summary() on a nma_nodesplit or nma_nodesplit_df object.

Details

Objects of class nodesplit_summary are tibble data frames, with one row for each node-split comparison and columns:

trt1, trt2

Treatments forming the comparison

summary

A list column containing nma_summary objects with the posterior summaries and draws for each of the node-splitting parameters

p_value

Bayesian p-value for inconsistency

dic

A list column containing nma_dic objects, giving the model fit statistics

The parameters included in summary are:

d_net

Network estimate from the corresponding consistency model, if available

d_dir

Direct estimate from the node-splitting model

d_ind

Indirect estimate from the node-splitting model

omega

Inconsistency factor $\omega = d_\mathrm{dir} - d_\mathrm{ind}$

tau

Heterogeneity standard deviation from the node-splitting model, if a random effects model was fitted

tau_consistency

Heterogeneity standard deviation from the corresponding consistency model, if available and if a random effects model was fitted

Description

A pairs() method for stan_nma objects, which calls bayesplot::mcmc_pairs() on the underlying stanfit object.

Usage

## S3 method for class 'stan_nma'
pairs(x, ..., pars, include = TRUE)

Arguments

Value

A grid of ggplot objects produced by bayesplot::mcmc_pairs().

Examples

## Not run: 
## Parkinson's mean off time reduction
park_net <- set_agd_arm(parkinsons,
                        study = studyn,
                        trt = trtn,
                        y = y,
                        se = se,
                        sample_size = n)

# Fitting a RE model
park_fit_RE <- nma(park_net,
                   trt_effects = "random",
                   prior_intercept = normal(scale = 100),
                   prior_trt = normal(scale = 100),
                   prior_het = half_normal(scale = 5))

# We see a small number of divergent transition errors
# These do not go away entirely when adapt_delta is increased

# Try to diagnose with a pairs plot
pairs(park_fit_RE, pars = c("mu[4]", "d[3]", "delta[4: 3]", "tau"))

# Transforming tau onto log scale
pairs(park_fit_RE, pars = c("mu[4]", "d[3]", "delta[4: 3]", "tau"),
      transformations = list(tau = "log"))

# The divergent transitions occur in the upper tail of the heterogeneity
# standard deviation. In this case, with only a small number of studies, there
# is not very much information to estimate the heterogeneity standard deviation
# and the prior distribution may be too heavy-tailed. We could consider a more
# informative prior distribution for the heterogeneity variance to aid
# estimation.

## End(Not run)

Description

Data frame containing the mean off-time reduction in patients given dopamine agonists as adjunct therapy in Parkinson's disease, from 7 trials comparing four active drugs and placebo (Dias et al. 2011).

Usage

parkinsons

Format

A data frame with 15 rows and 7 variables:

studyn

numeric study ID

trtn

numeric treatment code (placebo = 1)

y

mean off-time reduction

se

standard error

n

sample size

diff

mean difference vs. treatment in reference arm

se_diff

standard error of mean difference, see details

Details

This dataset may be analysed using either an arm-based likelihood using y and se, or a contrast-based likelihood using diff and se_diff (or a combination of the two across different studies).

The contrast-based data is formatted as described in set_agd_contrast(). That is, for the chosen reference arm in each study, the mean difference diff is set to NA, and se_diff is set to the standard error se of the outcome on the reference arm.

References

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Description

Two data frames, plaque_psoriasis_ipd and plaque_psoriasis_agd, containing (simulated) individual patient data from four studies and aggregate data from five studies (Phillippo 2019). Outcomes are binary success/failure to achieve 75%, 90%, or 100% reduction in symptoms on the Psoriasis Area and Severity Index (PASI) scale.

Usage

plaque_psoriasis_ipd

plaque_psoriasis_agd

Format

The individual patient data are contained in a data frame plaque_psoriasis_ipd with 4118 rows, one per individual, and 16 variables:

studyc

study name

trtc_long

treatment name (long format)

trtc

treatment name

trtn

numeric treatment code

pasi75

binary PASI 75 outcome

pasi90

binary PASI 90 outcome

pasi100

binary PASI 100 outcome

age

age (years)

bmi

body mass index (BMI)

pasi_w0

PASI score at week 0

male

male sex (TRUE or FALSE)

bsa

body surface area (percent)

weight

weight (kilograms)

durnpso

duration of psoriasis (years)

prevsys

previous systemic treatment (TRUE or FALSE)

psa

psoriatic arthritis (TRUE or FALSE)

The aggregate data are contained in a data frame plaque_psoriasis_agd with 15 rows, one per study arm, and 26 variables:

studyc

study name

trtc_long

treatment name (long format)

trtc

treatment name

trtn

numeric treatment code

pasi75_r, pasi75_n

PASI 75 outcome count and denominator

pasi90_r, pasi90_n

PASI 75 outcome count and denominator

pasi100_r, pasi100_n

PASI 75 outcome count and denominator

sample_size_w0

sample size at week zero

age_mean, age_sd

mean and standard deviation of age (years)

bmi_mean, bmi_sd

mean and standard deviation of BMI

pasi_w0_mean, pasi_w0_sd

mean and standard deviation of PASI score at week 0

male

percentage of males

bsa_mean, bsa_sd

mean and standard deviation of body surface area (percent)

weight_mean, weight_sd

mean and standard deviation of weight (kilograms)

durnpso_mean, durnpso_sd

mean and standard deviation of duration of psoriasis (years)

prevsys

percentage of individuals with previous systemic treatment

psa

percentage of individuals with psoriatic arthritis

An object of class data.frame with 15 rows and 26 columns.

References

Phillippo DM (2019). Calibration of Treatment Effects in Network Meta-Analysis using Individual Patient Data. Ph.D. thesis, University of Bristol. Available from https://research-information.bris.ac.uk/.

Description

For ML-NMR models, plot the estimated numerical integration error over the entire posterior distribution, as the number of integration points increases. See (Phillippo et al. 2020; Phillippo 2019) for details.

Usage

plot_integration_error(
  x,
  ...,
  stat = "violin",
  orientation = c("vertical", "horizontal", "x", "y"),
  show_expected_rate = TRUE
)

Arguments

Details

The total number of integration points is set by the n_int argument to add_integration(), and the intervals at which integration error is estimated are set by the int_thin argument to nma(). The typical convergence rate of Quasi-Monte Carlo integration (as used here) is $1/N$, which by default is displayed on the plot output.

The integration error at each thinning interval $N_\mathrm{thin}$ is estimated for each point in the posterior distribution by subtracting the final estimate (using all n_int points) from the estimate using only the first $N_\mathrm{thin}$ points.

Value

A ggplot object.

Note for survival models

This function is not supported for survival/time-to-event models. These do not save cumulative integration points for efficiency reasons (both time and memory).

Examples

## Plaque psoriasis ML-NMR
# Set up plaque psoriasis network combining IPD and AgD
library(dplyr)
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)

pso_ipd <- pso_ipd %>%
  mutate(# Variable transformations
    bsa = bsa / 100,
    prevsys = as.numeric(prevsys),
    psa = as.numeric(psa),
    weight = weight / 10,
    durnpso = durnpso / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker"),
    # Check complete cases for covariates of interest
    complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>%
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker")
  )

# Exclude small number of individuals with missing covariates
pso_ipd <- filter(pso_ipd, complete)

pso_net <- combine_network(
  set_ipd(pso_ipd,
          study = studyc,
          trt = trtc,
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd,
              study = studyc,
              trt = trtc,
              r = pasi75_r,
              n = pasi75_n,
              trt_class = trtclass)
)

# Print network details
pso_net

# Add integration points to the network
pso_net <- add_integration(pso_net,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  n_int = 64)


# Fit the ML-NMR model
pso_fit <- nma(pso_net, 
               trt_effects = "fixed",
               link = "probit",
               likelihood = "bernoulli2",
               regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
               class_interactions = "common",
               prior_intercept = normal(scale = 10),
               prior_trt = normal(scale = 10),
               prior_reg = normal(scale = 10),
               init_r = 0.1,
               QR = TRUE,
               # Set the thinning factor for saving the cumulative results
               # (This also sets int_check = FALSE)
               int_thin = 8)
pso_fit

# Plot numerical integration error
plot_integration_error(pso_fit)

Description

Produce plots comparing the prior and posterior distributions of model parameters.

Usage

plot_prior_posterior(
  x,
  ...,
  prior = NULL,
  post_args = list(),
  prior_args = list(),
  overlay = c("prior", "posterior"),
  ref_line = NA_real_
)

Arguments

Details

Prior distributions are displayed as lines, posterior distributions are displayed as histograms.

Value

A ggplot object.

Examples

## Smoking cessation NMA

# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)


# Plot prior vs. posterior, by default all parameters are plotted
plot_prior_posterior(smk_fit_RE)

# Plot prior vs. posterior for heterogeneity SD only
plot_prior_posterior(smk_fit_RE, prior = "het")

# Customise plot
plot_prior_posterior(smk_fit_RE, prior = "het",
                     prior_args = list(colour = "darkred", size = 2),
                     post_args = list(alpha = 0.6))

Description

Create a network plot from a nma_data network object.

Usage

## S3 method for class 'nma_data'
plot(
  x,
  ...,
  layout,
  circular,
  weight_edges = TRUE,
  weight_nodes = FALSE,
  show_trt_class = FALSE,
  nudge = 0
)

Arguments

Details

The default is equivalent to layout = "linear" and circular = TRUE, which places the treatment nodes on a circle in the order defined by the treatment factor variable. An alternative layout which may give good results for simple networks is "sugiyama", which attempts to minimise the number of edge crossings.

weight_nodes = TRUE requires that sample sizes have been specified for any aggregate data in the network, using the sample_size option of ⁠set_agd_*()⁠.

Value

A ggplot object, as produced by ggraph().

Examples

## Stroke prevention in atrial fibrillation
# Setting up the network
af_net <- set_agd_arm(atrial_fibrillation,
                      study = studyc,
                      trt = abbreviate(trtc, minlength = 3),
                      r = r,
                      n = n,
                      trt_class = trt_class)
af_net

# Basic plot
plot(af_net)

# Turn off weighting edges by number of studies
plot(af_net, weight_edges = FALSE)

# Turn on weighting nodes by sample size
plot(af_net, weight_nodes = TRUE)

# Colour treatment nodes by class
plot(af_net, weight_nodes = TRUE, show_trt_class = TRUE)

# Nudge the treatment labels away from the nodes
plot(af_net, weight_nodes = TRUE, show_trt_class = TRUE, nudge = 0.1)

# Output may be customised using standard ggplot commands
# For example, to display the legends below the plot:
plot(af_net, weight_nodes = TRUE, show_trt_class = TRUE) +
  ggplot2::theme(legend.position = "bottom",
                 legend.box = "vertical",
                 legend.margin = ggplot2::margin(0, 0, 0, 0),
                 legend.spacing = ggplot2::unit(0.5, "lines"))

# Choosing a different ggraph layout, hiding some legends
plot(af_net, weight_nodes = TRUE, show_trt_class = TRUE,
     layout = "star") +
  ggplot2::guides(edge_width = "none", size = "none")

Description

The plot() method for nma_dic objects produced by dic() produces several useful diagnostic plots for checking model fit and model comparison. Further detail on these plots and their interpretation is given by Dias et al. (2011).

Usage

## S3 method for class 'nma_dic'
plot(
  x,
  y,
  ...,
  show_uncertainty = TRUE,
  stat = "pointinterval",
  orientation = c("vertical", "horizontal", "x", "y")
)

Arguments

Details

When a single nma_dic object is given, a plot of the residual deviance contribution for each data point is produced. For a good fitting model, each data point is expected to have a residual deviance of 1; larger values indicate data points that are fit poorly by the model.

When two nma_dic objects are given, a "dev-dev" plot comparing the residual deviance contributions under each model is produced. Data points with residual deviance contributions lying on the line of equality are fit equally well under either model. Data points lying below the line of equality indicate better fit under the second model (y); conversely, data points lying above the line of equality indicate better fit under the first model (x). A common use case is to compare a standard consistency model (fitted using nma() with consistency = "consistency") with an unrelated mean effects (UME) inconsistency model (fitted using nma() with consistency = "ume"), to check for potential inconsistency.

See Dias et al. (2011) for further details.

Value

A ggplot object.

References

Dias S, Welton NJ, Sutton AJ, Ades AE (2011). “NICE DSU Technical Support Document 2: A generalised linear modelling framework for pair-wise and network meta-analysis of randomised controlled trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.

Examples

## Smoking cessation

# Run smoking FE NMA example if not already available
if (!exists("smk_fit_FE")) example("example_smk_fe", run.donttest = TRUE)


# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)


# Compare DIC of FE and RE models
(smk_dic_FE <- dic(smk_fit_FE))
(smk_dic_RE <- dic(smk_fit_RE))   # substantially better fit

# Plot residual deviance contributions under RE model
plot(smk_dic_RE)

# Further customisation is possible using ggplot commands
# For example, highlighting data points with residual deviance above a certain threshold
plot(smk_dic_RE) +
  ggplot2::aes(colour = ggplot2::after_stat(ifelse(y > 1.5, "darkorange", "black"))) +
  ggplot2::scale_colour_identity()

# Or by posterior probability, for example here a central probability of 0.6
# corresponds to a lower tail probability of (1 - 0.6)/2 = 0.2
plot(smk_dic_RE, .width = c(0.6, 0.95)) +
  ggplot2::aes(colour = ggplot2::after_stat(ifelse(ymin > 1, "darkorange", "black"))) +
  ggplot2::scale_colour_identity()

# Check for inconsistency using UME model


# Run smoking UME NMA example if not already available
if (!exists("smk_fit_RE_UME")) example("example_smk_ume", run.donttest = TRUE)


# Compare DIC
smk_dic_RE
(smk_dic_RE_UME <- dic(smk_fit_RE_UME))  # no difference in fit

# Compare residual deviance contributions with a "dev-dev" plot
plot(smk_dic_RE, smk_dic_RE_UME)

# By default the dev-dev plot can be a little cluttered
# Hiding the credible intervals
plot(smk_dic_RE, smk_dic_RE_UME, show_uncertainty = FALSE)

# Changing transparency
plot(smk_dic_RE, smk_dic_RE_UME, point_alpha = 0.5, interval_alpha = 0.1)

Description

The plot method for nma_summary objects is used to produce plots of parameter estimates (when called on a stan_nma object or its summary), relative effects (when called on the output of relative_effects()), absolute predictions (when called on the output of predict.stan_nma()), posterior ranks and rank probabilities (when called on the output of posterior_ranks() or posterior_rank_probs()).

Usage

## S3 method for class 'nma_summary'
plot(
  x,
  ...,
  stat = "pointinterval",
  orientation = c("horizontal", "vertical", "y", "x"),
  ref_line = NA_real_
)

## S3 method for class 'nma_parameter_summary'
plot(
  x,
  ...,
  stat = "pointinterval",
  orientation = c("horizontal", "vertical", "y", "x"),
  ref_line = NA_real_
)

## S3 method for class 'nma_rank_probs'
plot(x, ...)

## S3 method for class 'surv_nma_summary'
plot(x, ..., stat = "lineribbon")

Arguments

Details

Plotting is handled by ggplot2 and the stats and geoms provided in the ggdist package. As a result, the output is very flexible. Any plotting stats provided by ggdist may be used, via the argument stat.

The default uses ggdist::stat_pointinterval(), to produce medians and 95% Credible Intervals with 66% inner bands. Additional arguments in ... are passed to the ggdist stat, to customise the output. For example, to produce means and Credible Intervals, specify point_interval = "mean_qi". To produce an 80% Credible Interval with no inner band, specify .width = c(0, 0.8).

Alternative stats can be specified to produce different summaries. For example, specify stat = "[half]eye" to produce (half) eye plots, or stat = "histinterval" to produce histograms with intervals.

A full list of options and examples is found in the ggdist vignette vignette("slabinterval", package = "ggdist").

For survival/hazard/cumulative hazard curves estimated from survival models, the default uses ggdist::stat_lineribbon() which produces curves of posterior medians with 50%, 80%, and 95% Credible Interval bands. Again, additional arguments in ... are passed to the ggdist stat. For example, to produce posterior means and 95% Credible bands, specify point_interval = "mean_qi" and .width = 0.95.

A ggplot object is returned which can be further modified through the usual ggplot2 functions to add further aesthetics, geoms, themes, etc.

Value

A ggplot object.

Examples

## Smoking cessation

# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)


# Produce relative effects
smk_releff_RE <- relative_effects(smk_fit_RE)
plot(smk_releff_RE, ref_line = 0)

# Customise plot options
plot(smk_releff_RE, ref_line = 0, stat = "halfeye")

# Further customisation is possible with ggplot commands
plot(smk_releff_RE, ref_line = 0, stat = "halfeye", slab_alpha = 0.6) +
  ggplot2::aes(slab_fill = ggplot2::after_stat(ifelse(x < 0, "darkred", "grey60")))

# Produce posterior ranks
smk_rank_RE <- posterior_ranks(smk_fit_RE, lower_better = FALSE)
plot(smk_rank_RE)

# Produce rank probabilities
smk_rankprob_RE <- posterior_rank_probs(smk_fit_RE, lower_better = FALSE)
plot(smk_rankprob_RE)

# Produce cumulative rank probabilities
smk_cumrankprob_RE <- posterior_rank_probs(smk_fit_RE, lower_better = FALSE,
                                           cumulative = TRUE)
plot(smk_cumrankprob_RE)

# Further customisation is possible with ggplot commands
plot(smk_cumrankprob_RE) +
  ggplot2::facet_null() +
  ggplot2::aes(colour = Treatment)

Description

Produce summary plots of node-splitting models

Usage

## S3 method for class 'nodesplit_summary'
plot(
  x,
  ...,
  pars = "d",
  stat = "dens_overlay",
  orientation = c("horizontal", "vertical", "y", "x"),
  ref_line = NA_real_
)

Arguments

Details

Plotting is handled by ggplot2 and the stats and geoms provided in the ggdist package. As a result, the output is very flexible. Any plotting stats provided by ggdist may be used, via the argument stat. The default "dens_overlay" is a special exception, which uses ggplot2::geom_density(), to plot overlaid densities. Additional arguments in ... are passed to the ggdist stat, to customise the output.

Alternative stats can be specified to produce different summaries. For example, specify stat = "[half]eye" to produce (half) eye plots, or stat = "pointinterval" to produce point estimates and credible intervals.

A full list of options and examples is found in the ggdist vignette vignette("slabinterval", package = "ggdist").

A ggplot object is returned which can be further modified through the usual ggplot2 functions to add further aesthetics, geoms, themes, etc.

Value

A ggplot object.

See Also

summary.nma_nodesplit_df()

Examples

# Run smoking node-splitting example if not already available
if (!exists("smk_fit_RE_nodesplit")) example("example_smk_nodesplit", run.donttest = TRUE)


# Summarise the node-splitting results
(smk_nodesplit_summary <- summary(smk_fit_RE_nodesplit))

# Plot the node-splitting results
plot(smk_nodesplit_summary)

# Plot the inconsistency factors instead, change the plot stat to half-eye,
# and add a reference line at 0
plot(smk_nodesplit_summary, pars = "omega", stat = "halfeye", ref_line = 0)

# Plot a comparison of the heterogeneity under the node-split models vs.
# the consistency model
plot(smk_nodesplit_summary, pars = "tau")

Description

Produce posterior treatment rankings and rank probabilities from a fitted NMA model. When a meta-regression is fitted with effect modifier interactions with treatment, these will differ by study population.

Usage

posterior_ranks(
  x,
  newdata = NULL,
  study = NULL,
  lower_better = TRUE,
  probs = c(0.025, 0.25, 0.5, 0.75, 0.975),
  sucra = FALSE,
  summary = TRUE
)

posterior_rank_probs(
  x,
  newdata = NULL,
  study = NULL,
  lower_better = TRUE,
  cumulative = FALSE,
  sucra = FALSE
)