nma: Network meta-analysis models

Description

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

Usage

nma(
  network,
  consistency = c("consistency", "ume"),
  trt_effects = c("fixed", "random"),
  regression = NULL,
  class_interactions = c("common", "exchangeable", "independent"),
  likelihood = NULL,
  link = NULL,
  ...,
  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(),
  QR = FALSE,
  center = TRUE,
  adapt_delta = NULL,
  int_thin = max(network$n_int%/%10, 1)
)

Arguments

network

An nma_data object, as created by the functions set_*(), combine_network(), or add_integration()

consistency

Character string specifying the type of (in)consistency model to fit, currently either "consistency" or "ume"

trt_effects

Character string specifying either "fixed" or "random" effects

regression

A one-sided model formula, specifying the prognostic and effect-modifying terms for a regression model. Any references to treatment should use the .trt special variable, for example specifying effect modifier interactions as variable:.trt (see details).

class_interactions

Character string specifying whether effect modifier interactions are specified as "common", "exchangeable", or "independent".

likelihood

Character string specifying a likelihood, if unspecified will be inferred from the data

link

Character string specifying a link function, if unspecified will default to the canonical link

...

Further arguments passed to sampling(), such as iter, chains, cores, etc.

prior_intercept

Specification of prior distribution for the intercept

prior_trt

Specification of prior distribution for the treatment effects

prior_het

Specification of prior distribution for the heterogeneity (if trt_effects = "random")

prior_het_type

Character string specifying whether the prior distribution prior_het is placed on the heterogeneity standard deviation $τ$ ("sd", the default), variance $τ^{2}$ ("var"), or precision $1 / τ^{2}$ ("prec").

prior_reg

Specification of prior distribution for the regression coefficients (if regression formula specified)

prior_aux

Specification of prior distribution for the auxiliary parameter, if applicable (see details)

Logical scalar (default FALSE), whether to apply a QR decomposition to the model design matrix

center

Logical scalar (default TRUE), whether to center the (numeric) regression terms about the overall means

adapt_delta

See adapt_delta for details

int_thin

A single integer value, the thinning factor for returning cumulative estimates of integration error

Value

nma() returns a stan_nma object, nma.fit() returns a stanfit object.

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 $σ_{j k}$ 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.

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.

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.

Examples

Run this code

# NOT RUN {
## 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

# }
# NOT RUN {
# Fitting a fixed effect model
smk_fit_FE <- nma(smk_net, refresh = if (interactive()) 200 else 0,
                  trt_effects = "fixed",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100))

smk_fit_FE
# }
# NOT RUN {
# }
# NOT RUN {
# Fitting a random effects model
smk_fit_RE <- nma(smk_net, refresh = if (interactive()) 200 else 0,
                  trt_effects = "random",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100),
                  prior_het = normal(scale = 5))

smk_fit_RE
# }
# NOT RUN {
# }
# NOT RUN {
# Fitting an unrelated mean effects (inconsistency) model
smk_fit_RE_UME <- nma(smk_net, refresh = if (interactive()) 200 else 0,
                      consistency = "ume",
                      trt_effects = "random",
                      prior_intercept = normal(scale = 100),
                      prior_trt = normal(scale = 100),
                      prior_het = normal(scale = 5))

smk_fit_RE_UME
# }
# NOT RUN {

## 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 = 1000)

# }
# NOT RUN {
# 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, refresh = if (interactive()) 200 else 0,
               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
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning