# Run MBNMA dose-response models 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 if(!exists(".Rdpack.currefs")) .Rdpack.currefs \<-new.env();Rdpack::insert_citeOnly(keys="mawdsley2016;textual",package="MBNMAdose",cached_env=.Rdpack.currefs) . ```r mbnma.run( network, fun = dpoly(degree = 1), method = "common", regress = NULL, regress.effect = "common", class.effect = list(), UME = FALSE, sdscale = FALSE, cor = FALSE, omega = NULL, parameters.to.save = NULL, pd = "pd.kl", likelihood = NULL, link = NULL, priors = NULL, n.iter = 20000, n.chains = 3, n.burnin = floor(n.iter/2), n.thin = max(1, floor((n.iter - n.burnin)/1000)), autojags = FALSE, Rhat = 1.05, n.update = 10, model.file = NULL, jagsdata = NULL, ... ) ``` ## Arguments - `network`: An object of class `mbnma.network`. - `fun`: An object of `class("dosefun")` that specifies a functional form to be assigned to the dose-response. 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). - `regress`: A formula of effect modifiers (variables that interact with the treatment effect) to incorporate using Network Meta-Regression (E.g. `~ Population + Age`). All variables in the formula are modelled as interacting with the treatment effect (i.e. prognostic variables cannot be included in this way). Effects modifiers must be named variables in `network$data.ab` and must be identical within a study. Factor and character effect modifiers will be converted to a series of named dummy variables. - `regress.effect`: Indicates whether effect modification should be assumed to be `"common"` (assumed to be equal versus Placebo throughout the network), `"random"` (assumed to be exchangeable versus Placebo throughout the network), `"agent"` (assumed to be equal versus Placebo within each agent), or `"class"` (assumed to be equal versus Placebo within each class). - `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. Note that assuming class effects on some dose-response parameters may be unreasonable if the range of doses differ substantially across agents within a class. - `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. - `sdscale`: Logical object to indicate whether to write a model that specifies a reference SD for standardising when modelling using Standardised Mean Differences. Specifying `sdscale=TRUE` will therefore only modify the model if link function is set to SMD (`link="smd"`). - `cor`: A boolean object that indicates whether correlation should be modelled between relative effect dose-response parameters. This is automatically set to `FALSE` if class effects are modelled or if multiple dose-response functions are fitted. - `omega`: A scale matrix for the inverse-Wishart prior for the covariance matrix used to model the correlation between dose-response parameters (see Details for dose-response functions). `omega` must be a symmetric positive definite matrix with dimensions equal to the number of dose-response parameters modelled using relative effects (`"rel"`). If left as `NULL` (the default) a diagonal matrix with elements equal to 100 is used. - `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 if(!exists(".Rdpack.currefs")) .Rdpack.currefs \\<-new.env();Rdpack::insert_citeOnly(keys="spiegelhalter2002",package="MBNMAdose",cached_env=.Rdpack.currefs) . 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 if(!exists(".Rdpack.currefs")) .Rdpack.currefs \\<-new.env();Rdpack::insert_citeOnly(keys="plummer2008",package="MBNMAdose",cached_env=.Rdpack.currefs) . This will require running the model for additional iterations but is a more robust calculation for the effective number of parameters in non-linear models. * `popt` calculates pD using an optimism adjustment which allows for calculation of the penalized expected deviance if(!exists(".Rdpack.currefs")) .Rdpack.currefs \\<-new.env();Rdpack::insert_citeOnly(keys="plummer2008",package="MBNMAdose",cached_env=.Rdpack.currefs) . - `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"`), be assigned the value `"identity"` for an identity link function, or be assigned the value `"smd"` for modelling Standardised Mean Differences using an 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 if(!exists(".Rdpack.currefs")) .Rdpack.currefs \<-new.env();Rdpack::insert_citeOnly(keys="jagsmanual;textual",package="MBNMAdose",cached_env=.Rdpack.currefs) ). Note that normal distributions in JAGS are specified as $$ N(\mu, prec) $$ , where $$ prec = 1 / {\sigma^2} $$ . - `n.iter`: number of total iterations per chain (including burn in; default: 20000) - `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 diagnosticif(!exists(".Rdpack.currefs")) .Rdpack.currefs \<-new.env();Rdpack::insert_citeOnly(keys="gelmanrubin",package="MBNMAdose",cached_env=.Rdpack.currefs) . 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.05`. - `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). - `model.file`: The file path to a JAGS model (.jags file) that can be used to overwrite the JAGS model that is automatically written based on the specified options in `MBNMAdose`. Useful for adding further model flexibility. - `jagsdata`: A named list of the data objects to be used in the JAGS model. Only required if users are defining their own JAGS model using `model.file`. Format should match that of standard models fitted in `MBNMAdose` (see `mbnma$model.arg$jagsdata`) - `...`: Arguments to be sent to R2jags. ## Returns 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: ||| |:--|:--| |Parameters(s)/Parameter Prefix|Interpretation| |`` (e.g. `emax` )|The pooled effect for each dose-response parameter, as defined in dose-response functions. Will vary by agent if pooling is specified as `"rel"` in the dose-response function.| |`sd`|The between-study SD (heterogeneity) for relative effects, reported if `method="random"`| |`sd.` (e.g. `sd.emax` )|Between-study SD (heterogeneity) for absolute dose-response parameters specified as `"random"` .| |`` (e.g. `EMAX` )|The class effect within each class for a given dose-response parameter. These will be estimated by the model if specified in `class.effects` for a given dose-response parameter.| |`sd.` (e.g. `sd.EMAX` )|The within-class SD for different agents within the same class. Will be estimated by the model if any dose-response parameter in `class.effect` is set to `"random"` .| |`totresdev`|The residual deviance of the model| |`deviance`|The deviance of the model| 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. ## 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 `omega`. `cor` will automatically be set to `FALSE` if class effects are modelled. ## Dose-response parameter arguments ||| |:--|:--| |Argument|Model specification| |`"rel"`|Implies that **relative** effects should be pooled for this dose-response parameter separately for each agent in the network.| |`"common"`|Implies that all agents share the same common effect for this dose-response parameter.| |`"random"`|Implies that all agents share a similar (exchangeable) effect for this dose-response parameter. This approach allows for modelling of variability between agents.| |`numeric()`|Assigned a numeric value, indicating that this dose-response parameter should not be estimated from the data but should be assigned the numeric value determined by the user. This can be useful for fixing specific dose-response parameters (e.g. Hill parameters in Emax functions) to a single value.| ## Dose-response function Several general dose-response functions are provided, but a user-defined dose-response relationship can instead be used. As of version 0.4.0 dose-response functions are specified as an object of `class("dosefun")`. See help details for each of the functions below for the interpretation of specific dose-response parameters. Built-in dose-response functions are: * `dpoly()`: polynomial (e.g. for a linear model - `dpoly(degree=1)`) * `dloglin()`: log-linear * `dexp()`: exponential * `demax()`: (emax with/without a Hill parameter) * `dspline()`: splines (can fit B-splines (`type="bs"`), restricted cubic splines (`type="rcs"`), natural splines (`type="ns"`), or piecewise linear splines (`type="ls"`)) * `dfpoly()`: fractional polynomials * `dnonparam()`: Non-parametric monotonic function (`direction` can be either `"increasing"` or `"decreasing"`) following the method of if(!exists(".Rdpack.currefs")) .Rdpack.currefs \\<-new.env();Rdpack::insert_citeOnly(keys="owen2015;textual",package="MBNMAdose",cached_env=.Rdpack.currefs) * `duser()`: user-defined function * `dmulti()`: allows agent-specific dose-response functions to be fitted. A separate function must be provided for each agent in the network. ## Examples ```r # Using the triptans data network <- mbnma.network(triptans) ######## Dose-response functions ######## # Fit a dose-response MBNMA with a linear function # with common treatment effects result <- mbnma.run(network, fun=dpoly(degree=1), method="common") # Fit a dose-response MBNMA with a log-linear function # with random treatment effects result <- mbnma.run(network, fun=dloglin(), method="random") # Fit a dose-response MBNMA with a fractional polynomial function # with random treatment effects # with a probit link function result <- mbnma.run(network, fun=dfpoly(), method="random", link="probit") # Fit a user-defined function (quadratic) fun.def <- ~ (beta.1 * dose) + (beta.2 * (dose^2)) result <- mbnma.run(network, fun=duser(fun=fun.def), method="common") # Fit an Emax function # with a single random (exchangeable) parameter for ED50 # with common treatment effects result <- mbnma.run(network, fun=demax(emax="rel", ed50="random"), method="common") # Fit an Emax function with a Hill parameter # with a fixed value of 5 for the Hill parameter # with random relative effects result <- mbnma.run(network, fun=demax(hill=5), method="random") # Fit a model with natural cubic splines # with 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=dspline(type="ns", knots=c(0.1,0.3,0.6))) # Fit a model with different dose-response functions for each agent multifun <- dmulti(list(dloglin(), # for placebo (can be any function) demax(), # for eletriptan demax(), # for sumatriptan dloglin(), # for frovatriptan demax(), # for almotriptan demax(), # for zolmitriptan dloglin(), # for naratriptan demax())) # for rizatriptan multidose <- mbnma.run(network, fun=multifun) ########## Class effects ########## # Using the osteoarthritis dataset pain.df <- osteopain # Set a shared class (NSAID) only for Naproxcinod and Naproxen pain.df <- pain.df %>% dplyr::mutate( class = dplyr::case_when(agent %in% c("Naproxcinod", "Naproxen") ~ "NSAID", !agent %in% c("Naproxcinod", "Naproxen") ~ agent ) ) # Run an Emax MBNMA with a common class effect on emax painnet <- mbnma.network(pain.df) result <- mbnma.run(painnet, fun = demax(), class.effect = list(emax = "common")) ####### Priors ####### # Obtain priors from a fractional polynomial function result <- mbnma.run(network, fun=dfpoly(degree=1), method="random") print(result$model.arg$priors) # Change the prior distribution for the power newpriors <- list(power.1 = "dnorm(0,0.001) T(0,)") newpriors <- list(sd = "dnorm(0,0.5) T(0,)") result <- mbnma.run(network, fun=dfpoly(degree=1), method="random", priors=newpriors) ########## Sampler options ########## # Change the number of MCMC iterations, the number of chains, and the thin result <- mbnma.run(network, fun=dloglin(), method="random", n.iter=5000, n.thin=5, n.chains=4) # Calculate effective number of parameters via plugin method result <- mbnma.run(network, fun=dloglin(), method="random", pd="plugin") # Calculate effective number of parameters using penalized expected deviance result <- mbnma.run(network, fun=dloglin(), method="random", pd="popt") ####### Examine MCMC diagnostics (using mcmcplots or coda packages) ####### # Density plots mcmcplots::denplot(result) # Traceplots mcmcplots::traplot(result) # Caterpillar plots mcmcplots::caterplot(result, "rate") # Autocorrelation plots (using the coda package) coda::autocorr.plot(coda::as.mcmc(result)) ####### 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=demax(), 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) ``` ## References if(!exists(".Rdpack.currefs")) .Rdpack.currefs \<-new.env();Rdpack::insert_all_ref(.Rdpack.currefs)

mbnma.run function