Triglycerides Network Meta (TNM) data

Description

A systemically reviewed network meta data set on tryglyceride (TG) lowering drugs

Usage

data(TNM)

Format

A data frame with 73 rows and 15 variables

trial

trial identifier

treat

treatment indicator for placebo (PBO), simvastatin (S), atorvastatin (A), lovastatin (L), rosuvastatin (R), pravastatin (P), ezetimibe (E), simvastatin+ezetimibe (SE), atorvastatin+ezetimibe (AE), lovastatin+ezetimibe (LE), or pravastatin+ezetimibe (PE)

n

the number of participants in the study corresponding to the trial and treatment

ptg

percentage change from baseline in triglycerides (TG)

sdtg

sample standard deviation of percentage change in triglycerides (TG)

bldlc

baseline LDL-C

bhdlc

baseline HDL-C

btg

baseline triglycerides (TG)

age

age in years

white

the proportion of white participants

male

the proportion of male participants

bmi

body fat index

potencymed

the proportion of medium statin potency

potencyhigh

the proportion of high statin potency

durat

duration in weeks

Examples

data(TNM)

Fit Bayesian Network Meta-Regression Models

Description

This is a function the fits the model introduced in Bayesian Network Meta-Regression Models Using Heavy-Tailed Multivariate Random Effects with Covariate-Dependent Variances. The first seven arguments are required except ZCovariate. If not provided, ZCovariate will be assigned a vector of ones, rep(1, length(Outcome)). ZCovariate is the centerpiece of the modeling of variances and the heavy-tailed random effects distribution.

Usage

bayes_nmr(
  Outcome,
  SD,
  XCovariate,
  ZCovariate,
  Treat,
  Trial,
  Npt,
  prior = list(),
  mcmc = list(),
  control = list(),
  init = list(),
  Treat_order = NULL,
  Trial_order = NULL,
  scale_x = FALSE,
  verbose = FALSE
)

Arguments

Value

bayes_nmr returns an object of class "bayesnmr". The functions summary or print are used to obtain and print a summary of the results. The generic accessor function fitted extracts the posterior mean, posterior standard deviation, and the interval estimates of the value returned by bayes_nmr.

An object of class bayesnmr is a list containing the following components:

• Outcome - the aggregate response used in the function call.

• SD - the standard deviation used in the function call.

• Npt - the number of participants for ⁠(k,t)⁠ used in the function call.

• XCovariate - the aggregate design matrix for fixed effects used in the function call. Depending on scale_x, this may differ from the matrix provided at function call.

• ZCovariate - the aggregate design matrix for random effects. bayes_nmr will assign rep(1, length(Outcome)) if it was not provided at function call.

• Trial - the renumbered trial indicators. Depending on Trial_order, it may differ from the vector provided at function call.

• Treat - the renumbered treatment indicators. Depending on Treat_order, it may differ from the vector provided at function call.

• TrtLabels - the vector of treatment labels corresponding to the renumbered Treat. This is equivalent to Treat_order if it was given at function call.

• TrialLabels - the vector of trial labels corresponding to the renumbered Trial. This is equivalent to Trial_order if it was given at function call.

• K - the total number of trials.

• nT - the total number of treatments.

• scale_x - a Boolean indicating whether XCovariate has been scaled/standardized.

• prior - the list of hyperparameters used in the function call.

• control - the list of tuning parameters used for MCMC in the function call.

• mcmctime - the elapsed time for the MCMC algorithm in the function call. This does not include all the other preprocessing and post-processing outside of MCMC.

• mcmc - the list of MCMC specification used in the function call.

• mcmc.draws - the list containing the MCMC draws. The posterior sample will be accessible here.

Author(s)

Daeyoung Lim, daeyoung.lim@uconn.edu

References

Li, H., Chen, M. H., Ibrahim, J. G., Kim, S., Shah, A. K., Lin, J., & Tershakovec, A. M. (2019). Bayesian inference for network meta-regression using multivariate random effects with applications to cholesterol lowering drugs. Biostatistics, 20(3), 499-516.

Li, H., Lim, D., Chen, M. H., Ibrahim, J. G., Kim, S., Shah, A. K., & Lin, J. (2021). Bayesian network meta-regression hierarchical models using heavy-tailed multivariate random effects with covariate-dependent variances. Statistics in Medicine.

See Also

bmeta_analyze for using the Formula interface

Examples

library(metapack)
data(TNM)
groupInfo <- list(c("PBO"), c("R"))
nz <- length(groupInfo)
ns <- nrow(TNM)
XCovariate <- model.matrix(~ 0 + bldlc + bhdlc + btg + age +
 white + male + bmi + potencymed + potencyhigh + durat, data = TNM)
XCovariate <- scale(XCovariate, center = TRUE, scale = FALSE)
ZCovariate <- matrix(0, ns, nz)
for (j in 1:length(groupInfo)) {
    for (i in 1:ns) {
        if (TNM$treat[i] %in% groupInfo[[j]]) {
            ZCovariate[i, j] <- 1
        }
    }
}
addz <- scale(cbind(TNM$bldlc, TNM$btg), center=TRUE, scale=TRUE)
ZCovariate <- cbind(1, ZCovariate, addz)
theta_init <- c(0.05113, -1.38866, 1.09817, -0.85855, -1.12056, -1.14133,
             -0.22435, 3.63453, -2.09322, 1.07858, 0.80566, -40.76753,
             -45.07127, -28.27232, -44.14054, -28.13203, -19.19989,
             -47.21824, -51.31234, -48.46266, -47.71443)
set.seed(2797542)
fit <- bayes_nmr(TNM$ptg, TNM$sdtg, XCovariate, ZCovariate, TNM$treat,
    TNM$trial, TNM$n, prior = list(c01 = 1.0e05, c02 = 4, df = 3),
    mcmc = list(ndiscard = 1, nskip = 1, nkeep = 1),
    init = list(theta = theta_init),
    Treat_order = c("PBO", "S", "A", "L", "R", "P", "E", "SE",
         "AE", "LE", "PE"),
    scale_x = TRUE, verbose = FALSE)

Fit Bayesian Inference for Meta-Regression

Description

This is a function for running the Markov chain Monte Carlo algorithm for the Bayesian inference for multivariate meta-regression with a partially observed within-study sample covariance matrix model. The first six arguments are required. fmodel can be one of 5 numbers: 1, 2, 3, 4, and 5. The first model, fmodel = 1 denoted by M1, indicates that the \Sigma_{kt} are diagonal matrices with zero covariances. M2 indicates that \Sigma_{kt} are all equivalent but allowed to be full symmetric positive definite. M3 is where \Sigma_{kt} are allowed to differ across treatments, i.e., \Sigma_{kt}=\Sigma_t. M4 assumes thata the correlation matrix, \rho, is identical for all trials/treatments, but the variances are allowed to vary. Finally, M5 assumes a hierarchical model where (\Sigma_{kt} | \Sigma) follows an inverse-Wishart distribution with fixed degrees of freedom and scale matrix \Sigma. \Sigma then follows another inverse-Wishart distribution with fixed parameters.

Usage

bayes_parobs(
  Outcome,
  SD,
  XCovariate,
  WCovariate,
  Treat,
  Trial,
  Npt,
  fmodel = 1,
  prior = list(),
  mcmc = list(),
  control = list(),
  init = list(),
  Treat_order = NULL,
  Trial_order = NULL,
  group = NULL,
  group_order = NULL,
  scale_x = FALSE,
  verbose = FALSE
)

Arguments

Value

bayes_parobs returns an object of class "bayesparobs". The functions summary or print are used to obtain and print a summary of the results. The generic accessor function fitted extracts the posterior mean, posterior standard deviation, and the interval estimates of the value returned by bayes_parobs.

An object of class bayesparobs is a list containing the following components:

• Outcome - the aggregate response used in the function call.

• SD - the standard deviation used in the function call.

• Npt - the number of participants for ⁠(k,t)⁠ used in the function call.

• XCovariate - the aggregate design matrix for fixed effects used in the function call. Depending on scale_x, this may differ from the matrix provided at function call.

• WCovariate - the aggregate design matrix for random effects.

• Treat - the renumbered treatment indicators. Depending on Treat_order, it may differ from the vector provided at function call.

• Trial - the renumbered trial indicators. Depending on Trial_order, it may differ from the vector provided at function call.

• group - the renumbered grouping indicators in the function call. Depending on group_order, it may differ from the vector provided at function call. If group was missing at function call, bayes_parobs will assign NULL for group.

• TrtLabels - the vector of treatment labels corresponding to the renumbered Treat. This is equivalent to Treat_order if it was given at function call.

• TrialLabels - the vector of trial labels corresponding to the renumbered Trial. This is equivalent to Trial_order if it was given at function call.

• GroupLabels - the vector of group labels corresponding to the renumbered group. This is equivalent to group_order if it was given at function call. If group was missing at function call, bayes_parobs will assign NULL for GroupLabels.

• K - the total number of trials.

• T - the total number of treatments.

• fmodel - the model number as described here.

• scale_x - a Boolean indicating whether XCovariate has been scaled/standardized.

• prior - the list of hyperparameters used in the function call.

• control - the list of tuning parameters used for MCMC in the function call.

• mcmctime - the elapsed time for the MCMC algorithm in the function call. This does not include all the other preprocessing and post-processing outside of MCMC.

• mcmc - the list of MCMC specification used in the function call.

• mcmc.draws - the list containing the MCMC draws. The posterior sample will be accessible here.

Author(s)

Daeyoung Lim, daeyoung.lim@uconn.edu

References

Yao, H., Kim, S., Chen, M. H., Ibrahim, J. G., Shah, A. K., & Lin, J. (2015). Bayesian inference for multivariate meta-regression with a partially observed within-study sample covariance matrix. Journal of the American Statistical Association, 110(510), 528-544.

See Also

bmeta_analyze for using the Formula interface

Examples

library(metapack)
data("cholesterol")
Outcome <- model.matrix(~ 0 + pldlc + phdlc + ptg, data = cholesterol)
SD <- model.matrix(~ 0 + sdldl + sdhdl + sdtg, data = cholesterol)
Trial <- cholesterol$trial
Treat <- cholesterol$treat
Npt <- cholesterol$n
XCovariate <- model.matrix(~ 0 + bldlc + bhdlc + btg + age + durat +
 white + male + dm, data = cholesterol)
WCovariate <- model.matrix(~ treat, data = cholesterol)

fmodel <- 1
set.seed(2797542)
fit <- bayes_parobs(Outcome, SD, XCovariate, WCovariate, Treat, Trial,
   Npt, fmodel, mcmc = list(ndiscard = 1, nskip = 1, nkeep = 1),
   scale_x = TRUE, group = cholesterol$onstat, verbose = FALSE)

bmeta_analyze supersedes the previous two functions: bayes_parobs, bayes_nmr

Description

All other worker functions are superseded by this function, so that users can forget about the implementation details and focus on modeling. Meta-analytic data can be either aggregate or individual participant data (IPD). Aggregate data implies that the response consists of estimated effect sizes and their corresponding standard errors, whereas IPD is raw data. Data sets to be used for metapack should be formatted as follows:

The first treatment indicator is intentionally selected to be 3, a number greater than 1, to indicate that this data format works for both meta-regression and network meta-regression. Meta-regression refers to when trials included have 2 treatments (i.e., t = 0, 1 for all k), and the treatments are compared head to head. On the other hand, network meta-regression includes more than two treatments, where each trial can have a different set of treatments, allowing indirect comparison between treatments that are not compared head to head as long as consistency holds (see Higgins et al. (2012) for consistency).

bmeta_analyze() and bmeta_analyse() are synonyms.

Usage

bmeta_analyze(
  formula,
  data,
  prior = list(),
  mcmc = list(),
  control = list(),
  init = list()
)

bmeta_analyse(
  formula,
  data,
  prior = list(),
  mcmc = list(),
  control = list(),
  init = list()
)

Arguments

Details

bmeta_analyze currently subsumes two worker functions: bayes_parobs and bayes_nmr. bmeta_analyze offers a formula interface. All formulas are parsed using Formula. Formulas for bmeta_analyze are constrained to have a strict structure: one or two LHS, and two or three RHS. That is, lhs_1 ~ rhs_1 | rhs2 | rhs3 or lhs_1 | lhs_2 ~ rhs_1 | rhs2 | rhs3 (see Examples for more). The tilde (~) separates the LHS's and RHS's, each side further separated into parts by vertical bars (|). The meaning of each part is syntactically determined by its location inside the formula, like an English sentence. Therefore, all parts must come in the exact order as prescribed for bmeta_analyze to correctly configure your model.

• The first LHS, the responses, is required for all models.

• The second LHS is only required for aggregate models, corresponding to the standard deviations of the responses.

• The first RHS corresponds to fixed-effects covariates.

• The second RHS corresponds to the variables in either the random-effects matrix (w_{tk}' * \gamma_{k}') for multivariate meta-analysis or modeling the variances (\log\tau_{tk} = z_{tk}' * \phi) for univariate network meta-analysis.

• The third RHS corresponds to the treatment and trial indicators, and optionally the grouping variable if it exists. The order must be treat + trial + group, or treat + trial if no grouping exists. Variables here must be supplied in the exact order described; otherwise, model will not be correctly identified.

Internally, bmeta_analyze looks for three things: multivariate/univariate, meta-analyis/network meta-analysis, and aggregate/IPD.

• multivariate/univariate: the dimension of the response is explicit in the formula, and determines univariate versus multivariate.

• meta-analysis/network meta-analysis: the number of levels (nlevels) of treatments determines this. If treat is not already a factor variable, it is coerced to one.

• aggregate/IPD: bmeta_analyze looks for ns() in the first RHS. Aggregate models must provide the arm sample sizes using the function ns() (e.g., if n is the sample sizes, ⁠y1 + y2 | sd1 + sd2 ~ x1 + x2 + ns(n))⁠). If there is no ns(), IPD is assumed. Currently, IPD models are a work in progress and not supported yet.

Currently, only univariate/multivariate + meta-analysis and univariate + ⁠network meta-analysis⁠ are allowed. More models will be added in the future.

Value

bmeta_analyze returns a classed object of bsynthesis for Bayesian synthesis

Author(s)

Daeyoung Lim, daeyoung.lim@uconn.edu

References

Yao, H., Kim, S., Chen, M. H., Ibrahim, J. G., Shah, A. K., & Lin, J. (2015). Bayesian inference for multivariate meta-regression with a partially observed within-study sample covariance matrix. Journal of the American Statistical Association, 110(510), 528-544.

Li, H., Chen, M. H., Ibrahim, J. G., Kim, S., Shah, A. K., Lin, J., & Tershakovec, A. M. (2019). Bayesian inference for network meta-regression using multivariate random effects with applications to cholesterol lowering drugs. Biostatistics, 20(3), 499-516.

Li, H., Lim, D., Chen, M. H., Ibrahim, J. G., Kim, S., Shah, A. K., & Lin, J. (2021). Bayesian network meta-regression hierarchical models using heavy-tailed multivariate random effects with covariate-dependent variances. Statistics in Medicine.

See Also

bayes_parobs for multivariate meta-analysis, and bayes_nmr for univariate network meta-analysis.

Examples

set.seed(2797542)
data("cholesterol")
f_1 <- 'pldlc + phdlc + ptg | sdldl + sdhdl + sdtg ~ 0 + bldlc + bhdlc + btg +
  age + durat + white + male + dm + ns(n) | treat | treat + trial + onstat'
out_1 <- bmeta_analyze(as.formula(f_1), data = cholesterol,
  prior = list(model="NoRecovery"),
  mcmc = list(ndiscard = 3, nskip = 1, nkeep = 1),
  control=list(scale_x = TRUE, verbose=FALSE))

set.seed(2797542)
data("TNM")
TNM$group <- factor(match(TNM$treat, c("PBO", "R"), nomatch = 0))
f_2 <- 'ptg | sdtg ~
  0 + bldlc + bhdlc + btg + age + white + male + bmi +
  potencymed + potencyhigh + durat + ns(n) |
  scale(bldlc) + scale(btg) + group | treat  + trial'
out_2 <- bmeta_analyze(as.formula(f_2), data = TNM,
  mcmc = list(ndiscard = 1, nskip = 1, nkeep = 1),
  control=list(scale_x = TRUE, verbose=FALSE))

26 double-blind, randomized, active, or placebo-controlled clinical trials on patients with primary hypercholesterolemia sponsored by Merck & Co., Inc., Kenilworth, NJ, USA.

Description

A data set containing clinical trial on hypercholesterolemia including 26 trials and 2 treatment arms each, and other attributes of the participants

Usage

data(cholesterol)

Format

A data frame with 52 rows and 19 variables

study

study identifier

trial

trial identifier

treat

treatment indicator for Statin or Statin+Ezetimibe

n

the number of participants in the study arms corresponding to the trial and treatment

pldlc

aggregate percentage change in LDL-C

phdlc

aggregate percentage change from baseline in HDL-C

ptg

aggregate percentage change from baseline in triglycerides (TG)

sdldl

sample standard deviation of percentage change in LDL-C

sdhdl

sample standard deviation of percentage change in HDL-C

sdtg

sample standard deviation of percentage change in triglycerides (TG)

onstat

whether the participants were on Statin prior to the trial

bldlc

baseline LDL-C

bhdlc

baseline HDL-C

btg

baseline triglycerides (TG)

age

age in years

white

the proportion of white participants

male

the proportion of male participants

dm

the proportion of participants with diabetes mellitus

durat

duration in weeks

Examples

data(cholesterol)

get the posterior mean of fixed-effect coefficients

Description

get the posterior mean of fixed-effect coefficients

Usage

## S3 method for class 'bsynthesis'
coef(object, ...)

Arguments

Value

Coefficients extracted from the model object object

get fitted values

Description

get fitted values

Usage

## S3 method for class 'bayesnmr'
fitted(object, level = 0.95, HPD = TRUE, ...)

Arguments

Value

a list of fitted values

get fitted values

Description

get fitted values

Usage

## S3 method for class 'bayesparobs'
fitted(object, level = 0.95, HPD = TRUE, ...)

Arguments

Value

a list of fitted values

get the highest posterior density (HPD) interval

Description

get the highest posterior density (HPD) interval

Usage

hpd(object, parm, level = 0.95, HPD = TRUE)

Arguments

Details

A 100(1-\alpha)% HPD interval for \theta is given by

R(\pi_\alpha) = {\theta: \pi(\theta| D) \ge \pi_\alpha},

where \pi_\alpha is the largest constant that satisfies P(\theta \in R(\pi_\alpha)) \ge 1-\alpha. hpd computes the HPD interval from an MCMC sample by letting \theta_{(j)} be the jth smallest of the MCMC sample, {\theta_i} and denoting

R_j(n) = (\theta_{(j)}, \theta_{(j+[(1-\alpha)n])}),

for j=1,2,\ldots,n-[(1-\alpha)n]. Once \theta_i's are sorted, the appropriate j is chosen so that

\theta_{(j+[(1-\alpha)n])} - \theta_{(j)} = \min_{1\le j \leq n-[(1-\alpha)n]} (\theta_{(j+[(1-\alpha)n])} - \theta_{(j)}).

Value

dataframe containing HPD intervals for the parameters

References

Chen, M. H., & Shao, Q. M. (1999). Monte Carlo estimation of Bayesian credible and HPD intervals. Journal of Computational and Graphical Statistics, 8(1), 69-92.

get the highest posterior density (HPD) interval

Description

get the highest posterior density (HPD) interval

Usage

## S3 method for class 'bayesnmr'
hpd(object, parm, level = 0.95, HPD = TRUE)

Arguments

Value

dataframe containing HPD intervals for the parameters

get the highest posterior density (HPD) interval or equal-tailed credible interval

Description

get the highest posterior density (HPD) interval or equal-tailed credible interval

Usage

## S3 method for class 'bayesparobs'
hpd(object, parm, level = 0.95, HPD = TRUE)

Arguments

Value

dataframe containing HPD intervals for the parameters

metapack: a package for Bayesian meta-analysis and network meta-analysis

Description

The metapack package provides one category of functions: bayes.parobs and bayes.nmr

Multivariate Meta-Regression function

The bayes.parobs function fits the multivariate meta-regression model with partially observed sample covariance matrix to the given data.

Network Meta-Regression function

The bayes.nmr function fits the network meta-regression model with heavy-tailed random effects distribution to the given data.

compute the model comparison measures: DIC, LPML, or Pearson's residuals

Description

model_comp is a generic function that computes the model comparison measures (DIC and LPML) or the Pearson's residuals. Note that the Pearson's residuals are not available for bayes.nmr when df is either random or fixed but smaller than 2 since the variance of the random effects is not finite.

Usage

model_comp(object, type = "lpml", verbose = FALSE, ncores = NULL)

Arguments

Value

dataframe containing the compute the model comparison measures

get compute the model comparison measures

Description

get compute the model comparison measures

Usage

## S3 method for class 'bayesnmr'
model_comp(object, type = "lpml", verbose = FALSE, ncores = NULL)

Arguments

Value

dataframe containing the compute the model comparison measures

compute the model comparison measures

Description

compute the model comparison measures

Usage

## S3 method for class 'bayesparobs'
model_comp(object, type = "lpml", verbose = FALSE, ncores = NULL)

Arguments

Value

dataframe containing the compute the model comparison measures

helper function encoding trial sample sizes in formulas

Description

helper function encoding trial sample sizes in formulas

Usage

ns(x)

Arguments

get goodness of fit

Description

get goodness of fit

Usage

## S3 method for class 'bayesnmr'
plot(x, ...)

Arguments

Value

No return value

get goodness of fit

Description

get goodness of fit

Usage

## S3 method for class 'bayesparobs'
plot(x, ...)

Arguments

Value

No return value

plot the surface under the cumulative ranking curve (SUCRA)

Description

plot the surface under the cumulative ranking curve (SUCRA)

Usage

## S3 method for class 'sucra'
plot(x, legend.position = "none", ...)

Arguments

Value

No return value

Print results

Description

Print results

Usage

## S3 method for class 'bayesnmr'
print(x, level = 0.95, HPD = TRUE, ...)

Arguments

Value

No return value; print a summary of the output

Print results

Description

Print results

Usage

## S3 method for class 'bayesparobs'
print(x, level = 0.95, HPD = TRUE, ...)

Arguments

Value

No return value; print a summary of the output

get surface under the cumulative ranking curve (SUCRA)

Description

get surface under the cumulative ranking curve (SUCRA)

Usage

sucra(object)

Arguments

Value

a list containing SUCRA and the discrete rank probability matrix of size T by T

get surface under the cumulative ranking curve (SUCRA)

Description

get surface under the cumulative ranking curve (SUCRA)

Usage

## S3 method for class 'bayesnmr'
sucra(object)

Arguments

Value

a list containing SUCRA and the discrete rank probability matrix of size T by T

'summary' method for class "'bayesnmr'"

Description

'summary' method for class "'bayesnmr'"

Usage

## S3 method for class 'bayesnmr'
summary(object, level = 0.95, HPD = TRUE, ...)

Arguments

Value

does not return anything; print a summary of the output

summary method for class "bayesparobs"

Description

summary method for class "bayesparobs"

Usage

## S3 method for class 'bayesparobs'
summary(object, level = 0.95, HPD = TRUE, ...)

Arguments

Value

print summary for the model fit