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.
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
)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.
the aggregate mean of the responses for each arm of every study.
the standard deviation of the responses for each arm of every study.
the aggregate covariates for the fixed effects.
the aggregate covariates for the random effects.
the treatment identifiers. This is equivalent to the arm number of each study. The number of unique treatments must be equal across trials. The elements within will be coerced to consecutive integers.
the trial identifiers. This is equivalent to the arm labels in each study. The elements within will be coerced to consecutive integers
the number of observations/participants for a unique (k,t), or each arm of every trial.
the model number. The possible values for fmodel are 1 to 5, each indicating a different prior specification for \(\Sigma_{kt}\). It will default to M1, fmodel=1 if not specified at function call. See the following model descriptions. The objects enclosed in parentheses at the end of every bullet point are the hyperparameters associated with each model.
fmodel=1 - \(\Sigma_{kt} = diag(\sigma_{kt,11}^2,\ldots,\sigma_{kt,JJ}^2)\) where \(\sigma_{kt,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)
fmodel=2 - \(\Sigma_{kt}=\Sigma\) for every combination of \((k,t)\) 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)
fmodel=3 - \(\Sigma_{kt}=\Sigma_t\) and \(\Sigma_t^{-1}\sim Wish_{s_0}(\Sigma_0)\). This is a relaxed version of fmodel=2, allowing the correlation structure to differ across trials but forcing it to stay identical within a trial. (c0, dj0, s0, Omega0, Sigma0)
fmodel=4 - \(\Sigma_{kt}=\delta_{kt} \rho \delta_{kt}\) where \(\delta_{kt}=diag(\Sigma_{kt,11}^{1/2},\ldots,\Sigma_{kt,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)
fmodel=5 - The fifth model is hierarchical and thus may require more data than the others: \((\Sigma_{kt}^{-1}\mid \Sigma)\sim Wish_{\nu_0}((\nu_0-J-1)^{-1}\Sigma^{-1})\) and \(\Sigma \sim Wish_{d_0}(\Sigma_0)\). \(\Sigma_{kt}\) 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)
(Optional) a list of hyperparameters. Despite theta in every model, each fmodel, along with the group argument, requires a different set of hyperparameters. See fmodel for the model specifications.
(Optional) a list for MCMC specification. ndiscard is the number of burn-in iterations. nskip configures the thinning of the MCMC. For instance, if nskip=5, bayes_parobs will save the posterior sample every 5 iterations. nkeep is the size of the posterior sample. The total number of iterations will be ndiscard + nskip * nkeep.
(Optional) a list of tuning parameters for the Metropolis-Hastings algorithm. Rho, R, and delta are sampled through either localized Metropolis algorithm or delayed rejection robust adaptive Metropolis algorithm. *_stepsize with the asterisk replaced with one of the names above specifies the stepsize for determining the sample evaluation points in the localized Metropolis algorithm. sample_Rho can be set to FALSE to suppress the sampling of Rho for fmodel=4. When sample_Rho is FALSE, \(\rho\) will be fixed using the value given by the init argument, which defaults to \(0.5 I+0.511'\) where \(1\) is the vector of ones.
(Optional) a list of initial values for the parameters to be sampled: theta, gamR, Omega, and Rho. The initial value for Rho will be effective only if fmodel=4.
(Optional) a vector of unique treatments to be used for renumbering the Treat vector. The first element will be assigned treatment zero, potentially indicating placebo. If not provided, the numbering will default to an alphabetical/numerical order.
(Optional) a vector of unique trials. The first element will be assigned zero. If not provided, the numbering will default to an alphabetical/numerical order.
(Optional) a vector containing binary variables for \(u_{kt}\). If not provided, bayes_parobs will assume that there is no grouping and set \(u_{kt}=0\) for all (k,t).
(Optional) a vector of unique group labels. The first element will be assigned zero. If not provided, the numbering will default to an alphabetical/numerical order. group_order will take effect only if group is provided by the user.
(Optional) a logical variable indicating whether XCovariate should be scaled/standardized. The effect of setting this to TRUE is not limited to merely standardizing XCovariate. The following generic functions will scale the posterior sample of theta back to its original unit: plot, fitted, summary, and print.
(Optional) a logical variable indicating whether to print the progress bar during the MCMC sampling.
Daeyoung Lim, daeyoung.lim@uconn.edu
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.
bmeta_analyze for using the Formula interface
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)
Run the code above in your browser using DataLab