mbnma.run: Run MBNMA dose-response models

Description Usage Arguments Details Value Dose-response parameters Dose-response function References Examples

View source: R/run.functions.R

Description

Fits a Bayesian dose-response for model-based network meta-analysis (MBNMA) that can account for multiple doses of different agents by applying a desired dose-response function. Follows the methods of \insertCitemawdsley2016;textualMBNMAdose.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
mbnma.run(
  network,
  fun = "linear",
  beta.1 = "rel",
  beta.2 = "rel",
  beta.3 = "rel",
  beta.4 = "rel",
  method = "common",
  class.effect = list(),
  UME = FALSE,
  knots = 3,
  cor = TRUE,
  var.scale = NULL,
  user.fun = NULL,
  parameters.to.save = NULL,
  pd = "pv",
  parallel = FALSE,
  likelihood = NULL,
  link = NULL,
  priors = NULL,
  model.file = NULL,
  n.iter = 10000,
  n.chains = 3,
  n.burnin = floor(n.iter/2),
  n.thin = max(1, floor((n.iter - n.burnin)/1000)),
  autojags = FALSE,
  Rhat = 1.1,
  n.update = 10,
  arg.params = NULL,
  ...
)

Arguments

network

An object of class mbnma.network.

fun

A character vector specifying a functional form to be assigned to the dose-response. Options are given in details.

beta.1

Refers to dose-parameter(s) specified within the dose-response function(s). Can take either "rel", "common", "random", or be assigned a numeric value (see details).

beta.2

Refers to dose-parameter(s) specified within the dose-response function(s). Can take either "rel", "common", "random", or be assigned a numeric value (see details).

beta.3

Refers to dose-parameter(s) specified within the dose-response function(s). Can take either "rel", "common", "random", or be assigned a numeric value (see details).

beta.4

Refers to dose-parameter(s) specified within the dose-response function(s). Can take either "rel", "common", "random", or be assigned a numeric value (see details).

method

Can take either "common" or "random" to indicate whether relative effects should be modelled with between-study heterogeneity or not (see details).

class.effect

A list of named strings that determines which dose-response parameters to model with a class effect and what that effect should be ("common" or "random"). Element names should match dose-response parameter names (which will therefore depend on whether or not a wrapper function has been used for mbnma.run()). For example: list("beta.2"="fixed", "beta.3"="random") when using mbnma.run() or list("ed50"="fixed", "hill"="random") when using mbnma.emax.hill().

UME

A boolean object to indicate whether to fit an Unrelated Mean Effects model that does not assume consistency and so can be used to test if the consistency assumption is valid.

knots

The number/location of knots if a restricted cubic spline dose-response function is fitted (fun="rcs"). If a single number is given it indicates the number of knots (they will be equally spaced across the range of doses). If a numeric vector is given it indicates the location of the knots. Minimum number of knots is 3.

cor

A boolean object that indicates whether correlation should be modelled between relative effect dose-response parameters (TRUE) or not (FALSE). This is automatically set to FALSE if class effects are modelled or if multiple dose-response functions are fitted.

var.scale

A numeric vector indicating the relative scale of variances between correlated dose-response parameters when relative effects are modelled on more than one dose-response parameter and cor=TRUE (see details). Each element of the vector refers to the relative scale of each of the dose-response parameters that is modelled using relative effects.

user.fun

A formula specifying any relationship including dose and one/several of: beta.1, beta.2, beta.3, beta.4.

parameters.to.save

A character vector containing names of parameters to monitor in JAGS

pd

Can take either:

  • pv only pV will be reported (as automatically outputted by R2jags).

  • plugin calculates pD by the plug-in method \insertCitespiegelhalter2002MBNMAdose. It is faster, but may output negative non-sensical values, due to skewed deviances that can arise with non-linear models.

  • pd.kl calculates pD by the Kullback-Leibler divergence \insertCiteplummer2008MBNMAdose. This will require running the model for additional iterations but will always produce a positive result.

  • popt calculates pD using an optimism adjustment which allows for calculation of the penalized expected deviance \insertCiteplummer2008MBNMAdose

parallel

A boolean value that indicates whether JAGS should be run in parallel (TRUE) or not (FALSE). If TRUE then the number of cores to use is automatically calculated.

likelihood

A string indicating the likelihood to use in the model. Can take either "binomial", "normal" or "poisson". If left as NULL the likelihood will be inferred from the data.

link

A string indicating the link function to use in the model. Can take any link function defined within JAGS (e.g. "logit", "log", "probit", "cloglog") or be assigned the value "identity" for and identity link function. If left as NULL the link function will be automatically assigned based on the likelihood.

priors

A named list of parameter values (without indices) and replacement prior distribution values given as strings using distributions as specified in JAGS syntax (see examples).

model.file

A JAGS model written as a character object that can be used to overwrite the JAGS model that is automatically written based on the specified options.

n.iter

number of total iterations per chain (including burn in; default: 15000)

n.chains

number of Markov chains (default: 3)

n.burnin

length of burn in, i.e. number of iterations to discard at the beginning. Default is 'n.iter/2“, that is, discarding the first half of the simulations. If n.burnin is 0, jags() will run 100 iterations for adaption.

n.thin

thinning rate. Must be a positive integer. Set n.thin > 1`` to save memory and computation time if n.iter is large. Default is max(1, floor(n.chains * (n.iter-n.burnin) / 1000))“ which will only thin if there are at least 2000 simulations.

autojags

A boolean value that indicates whether the model should be continually updated until it has converged below a specific cutoff of Rhat

Rhat

A cutoff value for the Gelman-Rubin convergence diagnostic\insertCitegelmanrubinMBNMAdose. Unless all parameters have Rhat values lower than this the model will continue to sequentially update up to a maximum of n.update. Default is 1.1

n.update

The maximum number of updates. Each update is run for 1000 iterations, after which the Rhat values of all parameters are checked against Rhat. Default maximum updates is 10 (i.e. 10,000 additional iterations in total).

arg.params

Contains a list of arguments sent to mbnma.run() by dose-response specific wrapper functions

...

Arguments to be sent to R2jags.

Details

When relative effects are modelled on more than one dose-response parameter and cor = TRUE, correlation between the dose-response parameters is automatically estimated using a vague Wishart prior. This prior can be made slightly more informative by specifying the relative scale of variances between the dose-response parameters using var.scale. cor will automatically be set to FALSE if class effects are modelled or if a model is fitted with multiple dose-response functions.

Value

An object of S3 class(c("mbnma", "rjags")) containing parameter results from the model. Can be summarized by print() and can check traceplots using R2jags::traceplot() or various functions from the package mcmcplots.

Nodes that are automatically monitored (if present in the model) have the following interpretation. These will have an additional suffix that relates to the name/number of the dose-response parameter to which they correspond (e.g. d.ed50 or d.1):

If there are errors in the JAGS model code then the object will be a list consisting of two elements - an error message from JAGS that can help with debugging and model.arg, a list of arguments provided to mbnma.run() which includes jagscode, the JAGS code for the model that can help users identify the source of the error.

Dose-response parameters

Dose-response function

Several general dose-response functions are provided, but a user-defined dose-response relationship can instead be used.

Built-in dose-response functions are:

As of version 0.2.5, separate dose-response functions can be specified for different agents in the network by passing a character vector with multiple elements to fun. Each agent in network is assigned the dose-response function in the corresponding element in fun. fun must therefore be the same length as the number of agents in network. Dose-response parameters beta.1, beta.2, beta.3 and beta.4 refer to the corresponding dose-response parameters across the multiple functions in the following order: user, linear, exponential, emax, emax.hill.

This would mean that if fun included linear, exponential and emax within it then for the corresponding agents beta.1 would refer to linear slope parameters, beta.2 to exponential rate of growth/decay parameters, beta.3 to Emax parameters, and beta.4 to ED50 parameters.

References

\insertAllCited

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
# Using the triptans data
network <- mbnma.network(HF2PPITT)


######## Dose-response functions ########

# Fit a dose-response MBNMA with a linear function and common treatment effects
result <- mbnma.run(network, fun="linear", beta.1="rel", method="common")

# Fit a dose-response MBNMA with an exponential function and random treatment effects
result <- mbnma.run(network, fun="exponential", beta.1="rel", method="random")

# Fit a user-defined function (quadratic)
fun.def <- ~ (beta.1 * dose) + (beta.2 * (dose^2))
result <- mbnma.run(network, fun="user", user.fun=fun.def,
              beta.1="rel", beta.2="rel", method="common")

# Fit an Emax function with a single random (exchangeable) parameter estimated
#for ED50 and common treatment effects on relative Emax effects
result <- mbnma.run(network, fun="emax",
              beta.1="rel", beta.2="random", method="common")

# Fit an Emax function with a Hill parameter, with a fixed value for the Hill parameter
#provided to the model and random relative effects on Emax and ED50 (which will
#therefore be modelled with a correlation between them).
result <- mbnma.run(network, fun="emax.hill",
              beta.1="rel", beta.2="rel", beta.3=5, method="random")

# Fit a model with restricted cubic splines and 3 knots
#at 10% 30% and 60% quartiles of dose ranges
depnet <- mbnma.network(ssri) # Using the sSRI depression dataset
result <- mbnma.run(depnet, fun="rcs", knots=c(0.1,0.3,0.6))

# Fit a model with different dose-response functions for each agent
multidose <- mbnma.run(network, fun=c("emax", "emax", "emax", "exponential",
                 "emax", "emax", "exponential", "emax"))


########## Class effects ##########

# Generate a dataset with one class for active treatments and one for placebo
class.df <- HF2PPITT
class.df$class <- ifelse(class.df$agent=="placebo", "placebo", "active")
netclass <- mbnma.network(class.df)

painnet <- mbnma.network(osteopain_2wkabs)

# Fit an Emax function with random relative effects on Emax and ED50 and
#a common class effect on Emax
result <- mbnma.run(painnet, fun="emax", method="random",
              beta.1="rel", beta.2="rel",
              class.effect=list(beta.1="common"))


####### Priors #######

# Obtain priors from an Emax function with random relative effects on Emax and ED50
result <- mbnma.run(network, fun="emax",
              beta.1="rel", beta.2="rel", method="random")
print(result$model.arg$priors)

# Set new more informative prior distributions
newpriors <- list(sd = "dnorm(0,0.5) T(0,)",
                 inv.R = "dwish(Omega[,],100)")

result <- mbnma.run(network, fun="emax",
              beta.1="rel", beta.2="rel", method="random",
              priors=newpriors)


########## Sampler options ##########

# Change the number of MCMC iterations, the number of chains, and the thin
result <- mbnma.run(network, fun="exponential", beta.1="rel", method="random",
              n.iter=5000, n.thin=5, n.chains=4)

# Calculate effective number of parameters via plugin method
result <- mbnma.run(network, fun="exponential", beta.1="rel", method="random",
              pd="plugin")

# Calculate effective number of parameters via Kullback-Leibler method
result <- mbnma.run(network, fun="exponential", beta.1="rel", method="random",
              pd="pd.kl")


####### Examine MCMC diagnostics (using mcmcplots package) #######

# Density plots
mcmcplots::denplot(result)

# Traceplots
mcmcplots::traplot(result)

# Caterpillar plots
mcmcplots::caterplot(result, "d.1")


####### Automatically run jags until convergence is reached #########

# Rhat of 1.08 is set as the criteria for convergence
#on all monitored parameters
conv.res <- mbnma.run(network, fun="emax",
              beta.1="rel", beta.2="rel", method="random",
              n.iter=10000, n.burnin=9000,
              autojags=TRUE, Rhat=1.08, n.update=8)


########## Output ###########

# Print R2jags output and summary
print(result)
summary(result)

# Plot forest plot of results
plot(result)

MBNMAdose documentation built on Sept. 13, 2020, 5:08 p.m.