nma: Network meta-analysis models

In multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data

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", "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(),
  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, either "consistency", "ume", or "nodesplit"
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 (see details)
link	Character string specifying a link function, if unspecified will default to the canonical link (see details)
...	Further arguments passed to sampling(), such as iter, chains, cores, etc.
nodesplit	For consistency = "nodesplit", the comparison(s) to split in the node-splitting model(s). Either a length 2 vector giving the treatments in a single comparison, or a 2 column data frame listing multiple treatment comparisons to split in turn. By default, all possible comparisons will be chosen (see get_nodesplits()).
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 (see details)
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

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:

Data type	Likelihood	Link function
Binary	bernoulli, bernoulli2	logit, probit, cloglog
Count	binomial, binomial2	logit, probit, cloglog
Rate	poisson	log
Continuous	normal	identity, log
Ordered	ordered	logit, probit, cloglog

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 \insertCite@see @methods_papermultinma.

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 \insertCite@see @TSD2multinma.

Further details on each likelihood and link function are given by \insertCiteTSD2;textualmultinma.

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.

References

\insertAllCited

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


# 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

multinma documentation built on May 31, 2023, 5:46 p.m.