R: bmeta_analyze supersedes the previous two functions:...

bmeta_analyze {metapack}

R Documentation

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:

Outcome	SD	DesignM1	DesignM2	Trial indicator (`k`)	Treatment indicator (`t`)	n
`y_{13}`	`S_{13}`	`x_{13}`	`w_{13}`	1	3	1000
`y_{10}`	`S_{10}`	`x_{10}`	`w_{10}`	1	0	545
`y_{20}`	`S_{20}`	`x_{20}`	`w_{20}`	2	0	1200

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

`formula`	an object of class Formula: a symbolic description of the meta-analytic model to fit. For aggregate models, the vector of arm sample sizes must be provided using the function `ns()`. For example, `y1 + y2 \| sd1 + sd2 ~ x1 + x2 + ns(n)`—an incomplete formula only for illustration purposes. If no `ns()` is found, individual participant data (IPD) model is assumed.
`data`	a data frame, list, or environment (or an object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in `data`, the variables are taken from `environment(formula)`, typically the environment from which `bmeta_analyze` is called.
`prior`	an optional object that contains the hyperparameter values for the model. To see the complete list of hyperparameters for a specific model, refer to the corresponding worker function's help page, e.g., `help(bayes_parobs)` or `help(bayes_nmr)`. For meta-analysis, `model` is required in the `prior` argument, which is passed to `fmodel` as an integer. If the response is univariate, `NoRecovery` is the only valid option. `model="NoRecovery"` - `\Sigma_{tk} = diag(\sigma_{tk,11}^2,\ldots,\sigma_{tk,JJ}^2)` where `\sigma_{tk,jj}^2 \sim IG(a_0,b_0)` and `IG(a,b)` is the inverse-gamma distribution. This specification is useful if the user does not care about the correlation recovery. (`c0`, `dj0`, `a0`, `b0`, `Omega0`) `model="EquiCovariance"` - `\Sigma_{tk}=\Sigma` for every combination of `(t,k)` and `\Sigma^{-1}\sim Wish_{s_0}(\Sigma_0)`. This specification assumes that the user has prior knowledge that the correlation structure does not change across the arms included. (`c0`, `dj0`, `s0`, `Omega0`, `Sigma0`) `model="EquiWithinTreat"` - `\Sigma_{tk}=\Sigma_t` and `\Sigma_t^{-1}\sim Wish_{s_0}(\Sigma_0)`. This is a relaxed version of `model=2`, allowing the correlation structure to differ across trials but forcing it to stay identical within a trial. (`c0`, `dj0`, `s0`, `Omega0`, `Sigma0`) `model="EquiCorrelation"` - `\Sigma_{tk}=\delta_{tk} \rho \delta_{tk}` where `\delta_{tk}=diag(\Sigma_{tk,11}^{1/2},\ldots,\Sigma_{tk,JJ}^{1/2})`, and `\rho` is the correlation matrix. This specification allows the variances to vary across arms but requires that the correlations be the same. This is due to the lack of correlation information in the data, which would in turn lead to the nonidentifiability of the correlations if they were allowed to vary. However, this still is an ambitious model which permits maximal degrees of freedom in terms of variance and correlation estimation. (`c0`, `dj0`, `a0`, `b0`, `Omega0`) `model="Hierarchical"` - The fifth model is hierarchical and thus may require more data than the others: `(\Sigma_{tk}^{-1}\mid \Sigma)\sim Wish_{\nu_0}((\nu_0-J-1)^{-1}\Sigma^{-1})` and `\Sigma \sim Wish_{d_0}(\Sigma_0)`. `\Sigma_{tk}` encodes the within-treatment-arm variation while `\Sigma` captures the between-treatment-arm variation. The hierarchical structure allows the "borrowing of strength" across treatment arms. (`c0`, `dj0`, `d0`, `nu0`, `Sigma0`, `Omega0`) For network meta-analysis, `df` - the degrees of freedom of the multivariate t-distribution for the random effects. Any positive value can be assigned; if `df=Inf`, multivariate normal random effects will be assumed. `c01` - the variance of the fixed-effect coefficients' prior distribuiton, a multivariate normal distribution, i.e., `\theta \sim N(0, c_{1}I)`. `c02` - the variance of the random-effects' variance-related coefficients' prior distribution, a multivariate normal distribution, i.e., `\phi \sim N(0, c_{2}I)`. `a4`, `b4`, `a5`, `b5` - the hyperparameters related to when the degrees of freedom for the random effects are treated as unknown/random. `df` is then considered to follow `Ga(\nu_a, \nu_a/\nu_b)`, `\nu_a \sim Ga(a_4, b_4)`, and `\nu_b \sim IG(a_5, b_5)`. All gamma and inverse-gamma distributions are rate-parameterized.
`mcmc`	an optional object containing MCMC specification. `ndiscard` is the number of burn-in iterations. `nskip` configures the thinning of the MCMC. For instance, if `nskip=5`, parameters will be saved every 5 iterations. `nkeep` is the size of the posterior sample. The total number of iterations will be `ndiscard + nskip * nkeep`.
`control`	an optional object that contains the control tuning parameters for the Metropolis-Hastings algorithm. Similar to `prior`, the complete list of control parameters for a specific model is given in the corresponding worker function's help page (see `bayes_parobs` or `bayes_nmr`). These are the lists of available tuning parameters in `control` for meta-analysis and network meta-analysis. Keep in mind that `model` will render some irrelevant tuning parameters ineffective. Meta-analysis - `model` (string), `sample_Rho` (logical), `Rho_stepsize` (double), `R_stepsize` (double), `delta_stepsize` (double), `sample_Rho` (logical) Network meta-analysis - `sample_df` (logical), `sample_Rho` (logical), `lambda_stepsize` (double), `phi_stepsize` (double), `Rho_stepsize` (double)
`init`	(Optional) a list of initial values for the parameters to be sampled. The following is the list of available parameters for meta-analysis and network meta-analysis. Meta-analysis - `theta` (vector), `gamR` (matrix), `Omega` (matrix), `Rho` (matrix) Network meta-analysis - `theta` (vector), `phi` (vector), `sig2` (vector), `Rho` (matrix) The dimensions of the initial values must be conformable for matrix operations. If dimensions don't agree, `bmeta_analyze` will tell you the correct dimension.

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.

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

[Package metapack version 0.3 Index]