Run MBNMA time-course models — mb.run • MBNMAtime

Fits a Bayesian time-course model for model-based network meta-analysis (MBNMA) that can account for repeated measures over time within studies by applying a desired time-course function. Follows the methods of Pedder et al. (2019) .

Usage

mb.run(
  network,
  fun = tpoly(degree = 1),
  positive.scale = FALSE,
  intercept = NULL,
  link = "identity",
  sdscale = FALSE,
  parameters.to.save = NULL,
  rho = 0,
  covar = "varadj",
  omega = NULL,
  corparam = FALSE,
  class.effect = list(),
  UME = FALSE,
  pd = "pv",
  parallel = FALSE,
  priors = NULL,
  n.iter = 20000,
  n.chains = 3,
  n.burnin = floor(n.iter/2),
  n.thin = max(1, floor((n.iter - n.burnin)/1000)),
  model.file = NULL,
  jagsdata = NULL,
  ...
)

Arguments

network

An object of class "mb.network".

fun

An object of class "timefun" generated (see Details) using any of tloglin(), tpoly(), titp(), temax(), tfpoly(), tspline() or tuser()

positive.scale

A boolean object that indicates whether all continuous mean responses (y) are positive and therefore whether the baseline response should be given a prior that constrains it to be positive (e.g. for scales that cannot be <0).

intercept

A boolean object that indicates whether an intercept (written as alpha in the model) is to be included. If left as NULL (the default), an intercept will be included only for studies reporting absolute means, and will be excluded for studies reporting change from baseline (as indicated in network$cfb).

link

Can take either "identity" (the default), "log" (for modelling Ratios of Means (Friedrich et al. 2011) ) or "smd" (for modelling Standardised Mean Differences - although this also corresponds to an identity link function).

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").

parameters.to.save

A character vector containing names of parameters to monitor in JAGS

rho

The correlation coefficient when modelling within-study correlation between time points. The default is a string representing a prior distribution in JAGS, indicating that it be estimated from the data (e.g. rho="dunif(0,1)"). rho also be assigned a numeric value (e.g. rho=0.7), which fixes rho in the model to this value (e.g. for use in a deterministic sensitivity analysis). If set to rho=0 (the default) then this implies modelling no correlation between time points.

covar

A character specifying the covariance structure to use for modelling within-study correlation between time-points. This can be done by specifying one of the following:

"varadj" - a univariate likelihood with a variance adjustment to assume a constant correlation between subsequent time points (Jansen et al. 2015) . This is the default.
"CS" - a multivariate normal likelihood with a compound symmetry structure
"AR1" - a multivariate normal likelihood with an autoregressive AR1 structure

omega

DEPRECATED IN VERSION 0.2.3 ONWARDS (~uniform(-1,1) now used for correlation between parameters rather than a Wishart prior). A scale matrix for the inverse-Wishart prior for the covariance matrix used to model the correlation between time-course parameters (see Details for time-course functions). omega must be a symmetric positive definite matrix with dimensions equal to the number of time-course parameters modelled using relative effects (pool="rel"). If left as NULL (the default) a diagonal matrix with elements equal to 1 is used.

corparam

A boolean object that indicates whether correlation should be modeled between relative effect time-course parameters. Default is FALSE and this is automatically set to FALSE if class effects are modeled. Setting it to TRUE models correlation between time-course parameters. This can help identify parameters that are estimated poorly for some treatments by allowing sharing of information between parameters for different treatments in the network, but may also cause some shrinkage.

class.effect

A list of named strings that determines which time-course parameters to model with a class effect and what that effect should be ("common" or "random"). For example: list(emax="common", et50="random").

UME

Can take either TRUE or FALSE (for an unrelated mean effects model on all or no time-course parameters respectively) or can be a vector of parameter name strings to model as UME. For example: c("beta.1", "beta.2").

pd

Can take either:

pv only pV will be reported (as automatically outputted by R2jags).
plugin calculates pD by the plug-in method (Spiegelhalter et al. 2002) . It is faster, but may output negative non-sensical values, due to skewed deviances that can arise with non-linear models.
pd.kl (the default) calculates pD by the Kullback–Leibler divergence (Plummer 2008) . This will require running the model for additional iterations but will always produce a sensical result.
popt calculates pD using an optimism adjustment which allows for calculation of the penalized expected deviance (Plummer 2008)

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. Functions that involve updating the model (e.g. devplot(), fitplot()) cannot be used with models implemented in parallel.

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 Plummer (2017) ).

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.iteris large. Default ismax(1, floor(n.chains * (n.iter-n.burnin) / 1000))`` which will only thin if there are at least 2000 simulations.

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 MBNMAtime. 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 MBNMAtime (see mbnma$model.arg$jagsdata)

...

Arguments to be sent to R2jags.

Value

An object of S3 class c("mbnma", "rjags")`` containing parameter results from the model. Can be summarized by print()and can check traceplots usingR2jags::traceplot()or various functions from the packagemcmcplots`.#'

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 mb.run()

which includes jagscode, the JAGS code for the model that can help users identify the source of the error.

Time-course parameters

Nodes that are automatically monitored (if present in the model) have the same name as in the time-course function for named time-course parameters (e.g. emax). However, for named only as beta.1, beta.2, beta.3 or beta.4 parameters may have an alternative interpretation.

Details of the interpretation and model specification of different parameters can be shown by using the summary() method on an "mbnma" object generated by mb.run().

Parameters modelled using relative effects

If pooling is relative (e.g. pool.1="rel") for a given parameter then the named parameter (e.g. emax) or a numbered d parameter (e.g. d.1) corresponds to the pooled relative effect for a given treatment compared to the network reference treatment for this time-course parameter.
sd. followed by a named (e.g. emax, beta.1) is the between-study SD (heterogeneity) for relative effects, reported if pooling for a time-course parameter is relative (e.g. pool.1="rel") and the method for synthesis is random (e.g. method.1="random).
If class effects are modelled, parameters for classes are represented by the upper case name of the time-course parameter they correspond to. For example if class.effect=list(emax="random"), relative class effects will be represented by EMAX. The SD of the class effect (e.g. sd.EMAX, sd.BETA.1) is the SD of treatments within a class for the time-course parameter they correspond to.

Parameters modelled using absolute effects

If pooling is absolute (e.g. pool.1="abs") for a given parameter then the named parameter (e.g. emax) or a numbered beta parameter (e.g. beta.1) corresponds to the estimated absolute effect for this time-course parameter.
For an absolute time-course parameter if the corresponding method is common (e.g. method.1="common") the parameter corresponds to a single common parameter estimated across all studies and treatments. If the corresponding method is random (e.g. method.1="random") then parameter is a mean effect around which the study-level absolute effects vary with SD corresponding to sd. followed by the named parameter (e.g. sd.emax, sd.beta.1) .

Other model parameters

rho The correlation coefficient for correlation between time-points. Its interpretation will differ depending on the covariance structure specified in covar
totresdev The residual deviance of the model
deviance The deviance of the model

Time-course function

Several general time-course functions with up to 4 time-course parameters are provided, but a user-defined time-course relationship can instead be used. Details can be found in the respective help files for each function.

Available time-course functions are:

Log-linear: tloglin()
Polynomial: tpoly()
Integrated Two-Component Prediction (ITP): titp()
Emax: temax()
Fractional polynomial: tfpoly()
Splines (various spline types can be used): tspline()
User-defined: tuser()

Correlation between observations

When modelling correlation between observations using rho, values for rho must imply a positive semidefinite covariance matrix.

Advanced options

model.file and jagsdata can be used to run an edited JAGS model and dataset. This allows users considerably more modelling flexibility than is possible using the basic MBNMAtime syntax, though requires strong understanding of JAGS and the MBNMA modelling framework. Treatment-specific priors, meta-regression and bias-adjustment are all possible in this way, and it allows users to make use of the subsequent functions in MBNMAtime (plotting, prediction, ranking) whilst fitting these more complex models.

References

Friedrich JO, Adhikari NKJ, Beyene J (2011). “Ratio of means for analyzing continuous outcomes in meta-analysis performed as well as mean difference methods.” Journal of Clinical Epidemiology, 64(5), 556-564. doi:10.1016/j.jclinepi.2010.09.016 .

Jansen JP, Vieira MC, Cope S (2015). “Network meta-analysis of longitudinal data using fractional polynomials.” Stat Med, 34(15), 2294-311. ISSN 1097-0258 (Electronic) 0277-6715 (Linking), doi:10.1002/sim.6492 , https://pubmed.ncbi.nlm.nih.gov/25877808/.

Pedder H, Dias S, Bennetts M, Boucher M, Welton NJ (2019). “Modelling time-course relationships with multiple treatments: Model-Based Network Meta-Analysis for continuous summary outcomes.” Res Synth Methods, 10(2), 267-286.

Plummer M (2008). “Penalized loss functions for Bayesian model comparison.” Biostatistics, 9(3), 523-39. ISSN 1468-4357 (Electronic) 1465-4644 (Linking), https://pubmed.ncbi.nlm.nih.gov/18209015/.

Plummer M (2017). JAGS user manual. https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf.

Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A (2002). “Bayesian measures of model complexity and fit.” J R Statistic Soc B, 64(4), 583-639.

Examples

# \donttest{
# Create mb.network object
network <- mb.network(osteopain)
#> Reference treatment is `Pl_0`
#> Studies reporting change from baseline automatically identified from the data

# Fit a linear time-course MBNMA with:
# random relative treatment effects on the slope
mb.run(network, fun=tpoly(degree=1, pool.1="rel", method.1="random"))
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 417
#>    Unobserved stochastic nodes: 163
#>    Total graph size: 7701
#> 
#> Initializing model
#> 
#> Inference for Bugs model at "/tmp/RtmpAXd0fj/file176e64ba2561", fit using jags,
#>  3 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
#>  n.sims = 3000 iterations saved
#>            mu.vect sd.vect     2.5%      25%      50%      75%    97.5%  Rhat
#> d.1[1]       0.000   0.000    0.000    0.000    0.000    0.000    0.000 1.000
#> d.1[2]      -0.063   0.044   -0.149   -0.093   -0.063   -0.034    0.021 1.001
#> d.1[3]      -0.094   0.018   -0.129   -0.105   -0.094   -0.082   -0.061 1.001
#> d.1[4]      -0.095   0.045   -0.185   -0.124   -0.094   -0.065   -0.006 1.001
#> d.1[5]      -0.043   0.053   -0.150   -0.078   -0.044   -0.009    0.062 1.001
#> d.1[6]       0.003   0.070   -0.137   -0.042    0.003    0.051    0.138 1.001
#> d.1[7]      -0.147   0.035   -0.216   -0.170   -0.147   -0.124   -0.077 1.003
#> d.1[8]      -0.027   0.072   -0.169   -0.073   -0.028    0.022    0.111 1.001
#> d.1[9]      -0.289   0.049   -0.388   -0.321   -0.289   -0.258   -0.192 1.002
#> d.1[10]     -0.304   0.070   -0.441   -0.351   -0.304   -0.257   -0.165 1.001
#> d.1[11]     -0.076   0.037   -0.147   -0.101   -0.077   -0.051   -0.003 1.001
#> d.1[12]     -0.071   0.037   -0.146   -0.096   -0.071   -0.045    0.002 1.001
#> d.1[13]     -0.085   0.045   -0.174   -0.115   -0.085   -0.054    0.004 1.002
#> d.1[14]     -0.083   0.061   -0.199   -0.124   -0.082   -0.042    0.041 1.001
#> d.1[15]     -0.128   0.025   -0.176   -0.144   -0.128   -0.112   -0.078 1.001
#> d.1[16]     -0.133   0.039   -0.211   -0.159   -0.133   -0.107   -0.058 1.001
#> d.1[17]      0.115   0.064   -0.011    0.069    0.115    0.157    0.241 1.002
#> d.1[18]     -0.127   0.048   -0.220   -0.158   -0.126   -0.096   -0.035 1.001
#> d.1[19]     -0.122   0.079   -0.276   -0.176   -0.124   -0.069    0.032 1.002
#> d.1[20]     -0.114   0.049   -0.212   -0.147   -0.112   -0.081   -0.022 1.001
#> d.1[21]     -0.429   0.075   -0.580   -0.479   -0.428   -0.378   -0.285 1.001
#> d.1[22]     -0.170   0.042   -0.252   -0.198   -0.170   -0.142   -0.089 1.001
#> d.1[23]     -0.029   0.040   -0.107   -0.056   -0.029   -0.004    0.049 1.001
#> d.1[24]     -0.040   0.040   -0.121   -0.067   -0.039   -0.015    0.039 1.001
#> d.1[25]     -0.066   0.041   -0.148   -0.091   -0.065   -0.039    0.014 1.001
#> d.1[26]     -0.083   0.062   -0.209   -0.124   -0.083   -0.041    0.036 1.001
#> d.1[27]     -0.066   0.066   -0.195   -0.110   -0.067   -0.023    0.065 1.001
#> d.1[28]     -0.093   0.065   -0.222   -0.135   -0.093   -0.048    0.030 1.001
#> d.1[29]     -0.078   0.065   -0.206   -0.120   -0.079   -0.036    0.053 1.001
#> rho          0.000   0.000    0.000    0.000    0.000    0.000    0.000 1.000
#> sd.beta.1    0.072   0.009    0.056    0.066    0.072    0.078    0.093 1.001
#> totresdev 9299.479  16.771 9267.959 9287.866 9298.891 9310.876 9332.959 1.001
#> deviance  8419.440  16.771 8387.920 8407.827 8418.852 8430.837 8452.920 1.001
#>           n.eff
#> d.1[1]        1
#> d.1[2]     3000
#> d.1[3]     2500
#> d.1[4]     3000
#> d.1[5]     3000
#> d.1[6]     3000
#> d.1[7]     1600
#> d.1[8]     3000
#> d.1[9]     1500
#> d.1[10]    3000
#> d.1[11]    3000
#> d.1[12]    3000
#> d.1[13]    1200
#> d.1[14]    2200
#> d.1[15]    3000
#> d.1[16]    3000
#> d.1[17]    1800
#> d.1[18]    3000
#> d.1[19]    1900
#> d.1[20]    3000
#> d.1[21]    3000
#> d.1[22]    2100
#> d.1[23]    3000
#> d.1[24]    3000
#> d.1[25]    3000
#> d.1[26]    3000
#> d.1[27]    3000
#> d.1[28]    3000
#> d.1[29]    2700
#> rho           1
#> sd.beta.1  3000
#> totresdev  3000
#> deviance   3000
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 140.7 and DIC = 8559.5
#> DIC is an estimate of expected predictive error (lower deviance is better).

# Fit an emax time-course MBNMA with:
# fixed relative treatment effects on emax
# a common parameter estimated independently of treatment
# a common Hill parameter estimated independently of treatment
# a prior for the Hill parameter (normal with mean 0 and precision 0.1)
# data reported as change from baseline
result <- mb.run(network, fun=temax(pool.emax="rel", method.emax="common",
                                    pool.et50="abs", method.et50="common",
                                    pool.hill="abs", method.hill="common"),
                 priors=list(hill="dunif(0.5, 2)"),
                 intercept=TRUE)
#> 'et50' parameters must take positive values.
#>  Default half-normal prior restricts posterior to positive values.
#> 'hill' parameters must take positive values.
#>  Default half-normal prior restricts posterior to positive values.
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 417
#>    Unobserved stochastic nodes: 90
#>    Total graph size: 7740
#> 
#> Initializing model
#> 


#### commented out to prevent errors from JAGS version in github actions build ####
# Fit a log-linear MBNMA with:
# random relative treatment effects on the rate
# an autoregressive AR1 covariance structure
# modelled as standardised mean differences
# copdnet <- mb.network(copd)
# result <- mb.run(copdnet, fun=tloglin(pool.rate="rel", method.rate="random"),
#                covar="AR1", rho="dunif(0,1)", link="smd")



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

# Traceplots
# mcmcplots::traplot(result)

# Plots for assessing convergence
# mcmcplots::mcmcplot(result, c("rate", "sd.rate", "rho"))

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

# Print R2jags output and summary
print(result)
#> Inference for Bugs model at "/tmp/RtmpAXd0fj/file176e83cd95c", fit using jags,
#>  3 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
#>  n.sims = 3000 iterations saved
#>           mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
#> emax[1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
#> emax[2]    -0.642   0.107  -0.850  -0.714  -0.641  -0.572  -0.435 1.001  3000
#> emax[3]    -0.930   0.045  -1.018  -0.960  -0.930  -0.899  -0.841 1.001  3000
#> emax[4]    -1.037   0.103  -1.242  -1.107  -1.037  -0.965  -0.841 1.001  3000
#> emax[5]    -0.626   0.194  -1.018  -0.756  -0.626  -0.497  -0.249 1.001  2400
#> emax[6]    -0.120   0.142  -0.404  -0.213  -0.118  -0.026   0.153 1.001  3000
#> emax[7]    -1.228   0.073  -1.370  -1.279  -1.230  -1.177  -1.086 1.001  3000
#> emax[8]    -0.117   0.144  -0.401  -0.213  -0.118  -0.019   0.162 1.001  2900
#> emax[9]    -1.893   0.114  -2.118  -1.969  -1.892  -1.814  -1.671 1.003   850
#> emax[10]   -1.898   0.158  -2.207  -2.004  -1.895  -1.789  -1.592 1.002  1600
#> emax[11]   -0.876   0.062  -0.997  -0.918  -0.876  -0.836  -0.757 1.001  3000
#> emax[12]   -0.858   0.064  -0.982  -0.901  -0.857  -0.814  -0.734 1.001  3000
#> emax[13]   -1.061   0.077  -1.214  -1.111  -1.061  -1.009  -0.911 1.001  3000
#> emax[14]   -0.986   0.124  -1.238  -1.067  -0.983  -0.905  -0.749 1.001  3000
#> emax[15]   -1.316   0.066  -1.440  -1.362  -1.317  -1.270  -1.190 1.001  3000
#> emax[16]   -1.220   0.077  -1.374  -1.273  -1.219  -1.167  -1.071 1.001  3000
#> emax[17]    0.362   0.139   0.090   0.265   0.360   0.458   0.632 1.001  3000
#> emax[18]   -0.923   0.090  -1.097  -0.984  -0.923  -0.862  -0.753 1.001  3000
#> emax[19]   -1.289   0.240  -1.755  -1.454  -1.287  -1.130  -0.824 1.001  3000
#> emax[20]   -0.881   0.103  -1.084  -0.949  -0.882  -0.812  -0.678 1.001  2300
#> emax[21]   -2.634   0.208  -3.054  -2.773  -2.632  -2.496  -2.225 1.001  3000
#> emax[22]   -1.308   0.113  -1.531  -1.384  -1.308  -1.230  -1.086 1.002  1700
#> emax[23]   -0.214   0.069  -0.349  -0.260  -0.212  -0.170  -0.075 1.001  3000
#> emax[24]   -0.331   0.070  -0.470  -0.379  -0.329  -0.283  -0.195 1.002  2200
#> emax[25]   -0.765   0.073  -0.907  -0.812  -0.765  -0.715  -0.620 1.001  3000
#> emax[26]   -0.825   0.097  -1.020  -0.888  -0.825  -0.760  -0.632 1.001  3000
#> emax[27]   -0.785   0.120  -1.020  -0.867  -0.783  -0.703  -0.553 1.001  3000
#> emax[28]   -1.018   0.121  -1.260  -1.099  -1.016  -0.934  -0.788 1.001  3000
#> emax[29]   -0.831   0.123  -1.073  -0.916  -0.829  -0.746  -0.596 1.001  3000
#> et50        0.462   0.044   0.384   0.432   0.459   0.488   0.560 1.003   880
#> hill        0.589   0.078   0.502   0.527   0.567   0.632   0.782 1.008  3000
#> rho         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
#> totresdev 837.278  13.773 812.831 827.518 836.731 846.434 865.744 1.001  3000
#> deviance  -42.761  13.773 -67.208 -52.522 -43.308 -33.605 -14.295 1.001  3000
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 94.9 and DIC = 51.6
#> DIC is an estimate of expected predictive error (lower deviance is better).
summary(result)
#> ========================================
#> Time-course MBNMA
#> ========================================
#> 
#> Time-course function: emax
#> Data modelled with intercept
#> 
#> emax parameter
#> Pooling: relative effects
#> Method: common treatment effects
#> 
#> |Treatment |Parameter |  Median|    2.5%|   97.5%|
#> |:---------|:---------|-------:|-------:|-------:|
#> |Pl_0      |emax[1]   |  0.0000|  0.0000|  0.0000|
#> |Ce_100    |emax[2]   | -0.6409| -0.8503| -0.4353|
#> |Ce_200    |emax[3]   | -0.9304| -1.0179| -0.8407|
#> |Ce_400    |emax[4]   | -1.0367| -1.2423| -0.8411|
#> |Du_90     |emax[5]   | -0.6262| -1.0181| -0.2485|
#> |Et_10     |emax[6]   | -0.1182| -0.4044|  0.1527|
#> |Et_30     |emax[7]   | -1.2296| -1.3695| -1.0863|
#> |Et_5      |emax[8]   | -0.1185| -0.4011|  0.1623|
#> |Et_60     |emax[9]   | -1.8925| -2.1184| -1.6706|
#> |Et_90     |emax[10]  | -1.8955| -2.2072| -1.5916|
#> |Lu_100    |emax[11]  | -0.8764| -0.9972| -0.7565|
#> |Lu_200    |emax[12]  | -0.8573| -0.9816| -0.7340|
#> |Lu_400    |emax[13]  | -1.0609| -1.2138| -0.9111|
#> |Lu_NA     |emax[14]  | -0.9830| -1.2378| -0.7487|
#> |Na_1000   |emax[15]  | -1.3173| -1.4403| -1.1901|
#> |Na_1500   |emax[16]  | -1.2194| -1.3736| -1.0712|
#> |Na_250    |emax[17]  |  0.3597|  0.0904|  0.6316|
#> |Na_750    |emax[18]  | -0.9226| -1.0971| -0.7532|
#> |Ox_44     |emax[19]  | -1.2872| -1.7550| -0.8242|
#> |Ro_12     |emax[20]  | -0.8819| -1.0836| -0.6783|
#> |Ro_125    |emax[21]  | -2.6318| -3.0544| -2.2247|
#> |Ro_25     |emax[22]  | -1.3078| -1.5313| -1.0865|
#> |Tr_100    |emax[23]  | -0.2119| -0.3495| -0.0754|
#> |Tr_200    |emax[24]  | -0.3294| -0.4704| -0.1946|
#> |Tr_300    |emax[25]  | -0.7649| -0.9071| -0.6204|
#> |Tr_400    |emax[26]  | -0.8248| -1.0197| -0.6322|
#> |Va_10     |emax[27]  | -0.7834| -1.0197| -0.5529|
#> |Va_20     |emax[28]  | -1.0157| -1.2601| -0.7878|
#> |Va_5      |emax[29]  | -0.8291| -1.0729| -0.5962|
#> 
#> 
#> et50 parameter
#> Pooling: absolute effects
#> Method: common treatment effects
#> 
#> |Treatment |Parameter | Median|   2.5%|  97.5%|
#> |:---------|:---------|------:|------:|------:|
#> |Pl_0      |et50      | 0.4595| 0.3841| 0.5597|
#> |Ce_100    |et50      | 0.4595| 0.3841| 0.5597|
#> |Ce_200    |et50      | 0.4595| 0.3841| 0.5597|
#> |Ce_400    |et50      | 0.4595| 0.3841| 0.5597|
#> |Du_90     |et50      | 0.4595| 0.3841| 0.5597|
#> |Et_10     |et50      | 0.4595| 0.3841| 0.5597|
#> |Et_30     |et50      | 0.4595| 0.3841| 0.5597|
#> |Et_5      |et50      | 0.4595| 0.3841| 0.5597|
#> |Et_60     |et50      | 0.4595| 0.3841| 0.5597|
#> |Et_90     |et50      | 0.4595| 0.3841| 0.5597|
#> |Lu_100    |et50      | 0.4595| 0.3841| 0.5597|
#> |Lu_200    |et50      | 0.4595| 0.3841| 0.5597|
#> |Lu_400    |et50      | 0.4595| 0.3841| 0.5597|
#> |Lu_NA     |et50      | 0.4595| 0.3841| 0.5597|
#> |Na_1000   |et50      | 0.4595| 0.3841| 0.5597|
#> |Na_1500   |et50      | 0.4595| 0.3841| 0.5597|
#> |Na_250    |et50      | 0.4595| 0.3841| 0.5597|
#> |Na_750    |et50      | 0.4595| 0.3841| 0.5597|
#> |Ox_44     |et50      | 0.4595| 0.3841| 0.5597|
#> |Ro_12     |et50      | 0.4595| 0.3841| 0.5597|
#> |Ro_125    |et50      | 0.4595| 0.3841| 0.5597|
#> |Ro_25     |et50      | 0.4595| 0.3841| 0.5597|
#> |Tr_100    |et50      | 0.4595| 0.3841| 0.5597|
#> |Tr_200    |et50      | 0.4595| 0.3841| 0.5597|
#> |Tr_300    |et50      | 0.4595| 0.3841| 0.5597|
#> |Tr_400    |et50      | 0.4595| 0.3841| 0.5597|
#> |Va_10     |et50      | 0.4595| 0.3841| 0.5597|
#> |Va_20     |et50      | 0.4595| 0.3841| 0.5597|
#> |Va_5      |et50      | 0.4595| 0.3841| 0.5597|
#> 
#> 
#> hill parameter
#> Pooling: absolute effects
#> Method: common treatment effects
#> 
#> |Treatment |Parameter | Median|   2.5%|  97.5%|
#> |:---------|:---------|------:|------:|------:|
#> |Pl_0      |hill      | 0.5668| 0.5023| 0.7822|
#> |Ce_100    |hill      | 0.5668| 0.5023| 0.7822|
#> |Ce_200    |hill      | 0.5668| 0.5023| 0.7822|
#> |Ce_400    |hill      | 0.5668| 0.5023| 0.7822|
#> |Du_90     |hill      | 0.5668| 0.5023| 0.7822|
#> |Et_10     |hill      | 0.5668| 0.5023| 0.7822|
#> |Et_30     |hill      | 0.5668| 0.5023| 0.7822|
#> |Et_5      |hill      | 0.5668| 0.5023| 0.7822|
#> |Et_60     |hill      | 0.5668| 0.5023| 0.7822|
#> |Et_90     |hill      | 0.5668| 0.5023| 0.7822|
#> |Lu_100    |hill      | 0.5668| 0.5023| 0.7822|
#> |Lu_200    |hill      | 0.5668| 0.5023| 0.7822|
#> |Lu_400    |hill      | 0.5668| 0.5023| 0.7822|
#> |Lu_NA     |hill      | 0.5668| 0.5023| 0.7822|
#> |Na_1000   |hill      | 0.5668| 0.5023| 0.7822|
#> |Na_1500   |hill      | 0.5668| 0.5023| 0.7822|
#> |Na_250    |hill      | 0.5668| 0.5023| 0.7822|
#> |Na_750    |hill      | 0.5668| 0.5023| 0.7822|
#> |Ox_44     |hill      | 0.5668| 0.5023| 0.7822|
#> |Ro_12     |hill      | 0.5668| 0.5023| 0.7822|
#> |Ro_125    |hill      | 0.5668| 0.5023| 0.7822|
#> |Ro_25     |hill      | 0.5668| 0.5023| 0.7822|
#> |Tr_100    |hill      | 0.5668| 0.5023| 0.7822|
#> |Tr_200    |hill      | 0.5668| 0.5023| 0.7822|
#> |Tr_300    |hill      | 0.5668| 0.5023| 0.7822|
#> |Tr_400    |hill      | 0.5668| 0.5023| 0.7822|
#> |Va_10     |hill      | 0.5668| 0.5023| 0.7822|
#> |Va_20     |hill      | 0.5668| 0.5023| 0.7822|
#> |Va_5      |hill      | 0.5668| 0.5023| 0.7822|
#> 
#> 
#> 
#> Correlation between time points
#> Covariance structure: varadj 
#> Rho assigned a numeric value: 0
#> 
#> #### Model Fit Statistics ####
#> 
#> Effective number of parameters:
#> pD (pV) calculated using the rule, pD = var(deviance)/2 = 95
#> Deviance = -43
#> Residual deviance = 837
#> Deviance Information Criterion (DIC) = 52 
#> 

# Plot forest plot of results
plot(result)



###### Additional model arguments ######

# Use gout dataset
goutnet <- mb.network(goutSUA_CFBcomb)
#> Reference treatment is `Plac`
#> Studies reporting change from baseline automatically identified from the data

# Define a user-defined time-course relationship for use in mb.run
timecourse <- ~ exp(beta.1 * time) + (time^beta.2)

# Run model with:
# user-defined time-course function
# random relative effects on beta.1
# default common effects on beta.2
# default relative pooling on beta.1 and beta.2
# common class effect on beta.2
mb.run(goutnet, fun=tuser(fun=timecourse, method.1="random"),
       class.effect=list(beta.1="common"))
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 224
#>    Unobserved stochastic nodes: 121
#>    Total graph size: 4851
#> 
#> Initializing model
#> 
#> Inference for Bugs model at "/tmp/RtmpAXd0fj/file176e2e3b581b", fit using jags,
#>  3 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
#>  n.sims = 3000 iterations saved
#>               mu.vect sd.vect        2.5%         25%         50%         75%
#> D.1[1]          0.000   0.000       0.000       0.000       0.000       0.000
#> D.1[2]          8.047  21.200     -33.980      -5.781       6.886      20.366
#> D.1[3]         -6.036  31.027     -62.924     -23.709      -9.403       9.036
#> D.1[4]          2.203  27.596     -50.866     -15.753       1.816      21.495
#> D.1[5]        -19.347  26.870     -75.359     -36.599     -19.264      -1.494
#> D.1[6]        -19.204  20.516     -58.088     -32.177     -18.491      -7.257
#> D.1[7]        -16.018  30.387     -80.449     -35.759     -16.756       4.940
#> d.1[1]          0.000   0.000       0.000       0.000       0.000       0.000
#> d.1[2]          8.047  21.200     -33.980      -5.781       6.886      20.366
#> d.1[3]          8.047  21.200     -33.980      -5.781       6.886      20.366
#> d.1[4]          8.047  21.200     -33.980      -5.781       6.886      20.366
#> d.1[5]          8.047  21.200     -33.980      -5.781       6.886      20.366
#> d.1[6]         -6.036  31.027     -62.924     -23.709      -9.403       9.036
#> d.1[7]          2.203  27.596     -50.866     -15.753       1.816      21.495
#> d.1[8]          2.203  27.596     -50.866     -15.753       1.816      21.495
#> d.1[9]          2.203  27.596     -50.866     -15.753       1.816      21.495
#> d.1[10]         2.203  27.596     -50.866     -15.753       1.816      21.495
#> d.1[11]       -19.347  26.870     -75.359     -36.599     -19.264      -1.494
#> d.1[12]       -19.204  20.516     -58.088     -32.177     -18.491      -7.257
#> d.1[13]       -19.204  20.516     -58.088     -32.177     -18.491      -7.257
#> d.1[14]       -19.204  20.516     -58.088     -32.177     -18.491      -7.257
#> d.1[15]       -19.204  20.516     -58.088     -32.177     -18.491      -7.257
#> d.1[16]       -16.018  30.387     -80.449     -35.759     -16.756       4.940
#> d.1[17]       -16.018  30.387     -80.449     -35.759     -16.756       4.940
#> d.1[18]       -16.018  30.387     -80.449     -35.759     -16.756       4.940
#> d.1[19]       -16.018  30.387     -80.449     -35.759     -16.756       4.940
#> d.2[1]          0.000   0.000       0.000       0.000       0.000       0.000
#> d.2[2]         -4.124  29.585     -62.986     -23.284      -3.580      14.937
#> d.2[3]          3.214  29.965     -54.600     -17.142       3.345      23.284
#> d.2[4]        -16.218   9.861     -41.224     -21.258     -13.755      -8.741
#> d.2[5]          2.980  29.513     -52.998     -17.064       2.961      23.146
#> d.2[6]        -18.673  27.342     -72.928     -34.486     -19.041      -6.894
#> d.2[7]         -3.524  30.274     -65.512     -23.756      -3.162      16.828
#> d.2[8]         -7.311  27.497     -63.185     -25.868      -6.906      10.939
#> d.2[9]         -2.213  30.311     -64.007     -22.318      -1.467      18.877
#> d.2[10]        -4.845  28.884     -63.772     -24.042      -4.059      13.890
#> d.2[11]       -17.266  25.537     -68.699     -34.079     -16.994      -0.271
#> d.2[12]       -11.878  26.469     -65.694     -29.007     -10.998       5.989
#> d.2[13]       -17.549  26.484     -68.051     -34.409     -18.006      -2.351
#> d.2[14]       -15.545  21.625     -60.060     -29.436     -14.931      -1.564
#> d.2[15]       -16.123  24.497     -67.627     -31.838     -15.188       0.287
#> d.2[16]         3.606  30.278     -56.255     -16.396       2.859      24.123
#> d.2[17]       -24.814  21.224     -71.678     -36.976     -22.192     -11.040
#> d.2[18]       -23.021  19.986     -64.380     -35.057     -20.983     -10.501
#> d.2[19]       -26.040  21.391     -73.237     -38.837     -23.248     -11.851
#> rho             0.000   0.000       0.000       0.000       0.000       0.000
#> sd.beta.1       3.509   2.491       0.330       1.495       3.010       4.987
#> totresdev 2013066.484   6.522 2013056.406 2013062.454 2013064.973 2013068.762
#> deviance  2012372.295   6.522 2012362.218 2012368.265 2012370.784 2012374.574
#>                 97.5%  Rhat n.eff
#> D.1[1]          0.000 1.000     1
#> D.1[2]         57.647 1.197    16
#> D.1[3]         63.587 1.381     9
#> D.1[4]         53.186 1.164    19
#> D.1[5]         33.190 1.011   360
#> D.1[6]         25.992 1.183    30
#> D.1[7]         48.398 1.604     7
#> d.1[1]          0.000 1.000     1
#> d.1[2]         57.647 1.197    16
#> d.1[3]         57.647 1.197    16
#> d.1[4]         57.647 1.197    16
#> d.1[5]         57.647 1.197    16
#> d.1[6]         63.587 1.381     9
#> d.1[7]         53.186 1.164    19
#> d.1[8]         53.186 1.164    19
#> d.1[9]         53.186 1.164    19
#> d.1[10]        53.186 1.164    19
#> d.1[11]        33.190 1.011   360
#> d.1[12]        25.992 1.183    30
#> d.1[13]        25.992 1.183    30
#> d.1[14]        25.992 1.183    30
#> d.1[15]        25.992 1.183    30
#> d.1[16]        48.398 1.604     7
#> d.1[17]        48.398 1.604     7
#> d.1[18]        48.398 1.604     7
#> d.1[19]        48.398 1.604     7
#> d.2[1]          0.000 1.000     1
#> d.2[2]         54.111 1.002  1800
#> d.2[3]         62.723 1.001  3000
#> d.2[4]         -4.523 1.002  1200
#> d.2[5]         61.407 1.001  3000
#> d.2[6]         42.845 1.185    17
#> d.2[7]         55.059 1.001  2400
#> d.2[8]         46.511 1.001  3000
#> d.2[9]         56.321 1.001  2800
#> d.2[10]        52.863 1.001  3000
#> d.2[11]        30.903 1.001  3000
#> d.2[12]        39.175 1.002  1200
#> d.2[13]        38.500 1.157    18
#> d.2[14]        24.795 1.002  1000
#> d.2[15]        29.802 1.001  2600
#> d.2[16]        64.519 1.004   590
#> d.2[17]        17.332 1.096    31
#> d.2[18]        18.625 1.129    24
#> d.2[19]        15.137 1.074    47
#> rho             0.000 1.000     1
#> sd.beta.1       9.332 1.081    37
#> totresdev 2013082.871 1.000     1
#> deviance  2012388.682 1.000     1
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 18.6 and DIC = 2012389.3
#> DIC is an estimate of expected predictive error (lower deviance is better).

# Fit a log-linear MBNMA
# with variance adjustment for correlation between time-points
result <- mb.run(network, fun=tloglin(),
                 rho="dunif(0,1)", covar="varadj")
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 417
#>    Unobserved stochastic nodes: 89
#>    Total graph size: 7303
#> 
#> Initializing model
#> 
# }