multinma (version 0.1.3)

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 \(\tau\) ("sd", the default), variance \(\tau^2\) ("var"), or precision \(1/\tau^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

QR

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.

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

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