bmr {bayesmeta} R Documentation

## Bayesian random-effects meta-regression

### Description

This function allows to derive the posterior distribution of the parameters in a random-effects meta-regression and provides functions to evaluate joint and marginal posterior probability distributions, etc.

### Usage

  bmr(y, ...)
## Default S3 method:
bmr(y, sigma, labels = names(y),
X = matrix(1.0, nrow=length(y), ncol=1,
dimnames=list(labels,"intercept")),
tau.prior = "uniform",
beta.prior.mean = NULL,
beta.prior.sd = NULL,
beta.prior.cov = diag(beta.prior.sd^2,
nrow=length(beta.prior.sd),
ncol=length(beta.prior.sd)),
interval.type = c("shortest", "central"),
delta = 0.01, epsilon = 0.0001,

#### Computation

The bmr() function utilizes the same computational method as the bayesmeta() function to derive the posterior distribution, namely, the DIRECT algorithm. Numerical accuracy of the computations is determined by the ‘delta’ and ‘epsilon’ arguments (Roever and Friede, 2017).

A slight difference between the bayesmeta() and bmr() implementations exists in the determination of the grid approximation within the DIRECT algorithm. While bmr() considers divergences w.r.t. the conditional posterior distributions p(\beta|\tau), bayesmeta() in addition considers divergences w.r.t. the shrinkage estimates, which in general leads to a denser binning (as one can see from the numbers of bins required; see the example below). A denser binning within the bmr() function may be achieved by reducing the ‘delta’ argument.

### Value

A list of class bmr containing the following elements:

 y vector of estimates (the input data). sigma vector of standard errors corresponding to y (input data). X the regressor matrix. k number of data points (length of y, or rows of X). d number of coefficients (columns of X). labels vector of labels corresponding to y and sigma. variables variable names for the \beta coefficients (determined by the column names of the supplied X argument). tau.prior the prior probability density function for \tau. tau.prior.proper a logical flag indicating whether the heterogeneity prior appears to be proper (which is judged based on an attempted numerical integration of the density function). beta.prior a list containing the prior mean vector and covariance matrix for the coefficients \beta. beta.prior.proper a logical vector (of length d) indicating whether the corresponding \beta coefficient's prior is proper (i.e., finite prior mean and variance were specified). dprior a function(tau, beta, which.beta, log=FALSE) of \tau and/or \beta parameters, returning either the joint or marginal prior probability density, depending on which parameter(s) is/are provided. likelihood a function(tau, beta, which.beta) \tau and/or \beta, returning either the joint or marginal likelihood, depending on which parameter(s) is/are provided. dposterior, pposterior, qposterior, rposterior, post.interval functions of \tau and/or \beta parameters, returning either the joint or marginal posterior probability density, (depending on which parameter(s) is/are provided), or cumulative distribution function, quantile function, random numbers or posterior intervals. dpredict, ppredict, qpredict, rpredict, pred.interval functions of \beta returning density, cumulative distribution function, quantiles, random numbers, or intervals for the predictive distribution. This requires specification of x values to indicate what covariable values to consider. Use of ‘mean=TRUE’ (the default) yields predictions for the mean (x'\beta values), setting it to FALSE yields predictions (\theta values). dshrink, pshrink, qshrink, rshrink, shrink.interval functions of \theta yielding density, cumulative distribution, quantiles, random numbers or posterior intervals for the shrinkage estimates of the individual \theta_i parameters corresponding to the supplied y_i data values (i=1,\ldots,k). May be identified using the ‘which’ argument via its index (i) or a character string giving the corresponding study label. post.moments a function(tau) returning conditional posterior moments (mean and covariance) of \beta as a function of \tau. pred.moments a function(tau, x, mean=TRUE) returning conditional posterior predictive moments (means and standard deviations) as a function of \tau. shrink.moments a function(tau, which) returning conditional moments (means and standard deviations of shrinkage distributions) as a function of \tau. summary a matrix listing some summary statistics, namely marginal posterior mode, median, mean, standard deviation and a (shortest) 95% credible intervals, of the marginal posterior distributions of \tau and \beta_i. interval.type the interval.type input argument specifying the type of interval to be returned by default. ML a matrix giving joint and marginal maximum-likelihood estimates of (\tau,\beta). MAP a matrix giving joint and marginal maximum-a-posteriori estimates of (\tau,\beta). theta a matrix giving the ‘shrinkage estimates’, i.e, summary statistics of the trial-specific means \theta_i. marginal.likelihood the marginal likelihood of the data (this number can only be computed if proper effect and heterogeneity priors are specified). bayesfactor Bayes factors and minimum Bayes factors for the hypotheses of \tau=0 and \beta_i=0. These depend on the marginal likelihood and hence can only be computed if proper effect and/or heterogeneity priors are specified. support a list giving the \tau support points used internally in the grid approximation, along with their associated weights, and conditional mean and covariance of \beta. delta, epsilon the ‘delta’ and ‘epsilon’ input parameter determining numerical accuracy. rel.tol.integrate, abs.tol.integrate, tol.uniroot the input parameters determining the numerical accuracy of the internally used integrate() and uniroot() functions. call an object of class call giving the function call that generated the bmr object. init.time the computation time (in seconds) used to generate the bmr object.

### Author(s)

Christian Roever christian.roever@med.uni-goettingen.de

### References

C. Roever, T. Friede. Using the bayesmeta R package for Bayesian random-effects meta-regression. Computer Methods and Programs in Biomedicine, 299:107303, 2023. doi:10.1016/j.cmpb.2022.107303.

C. Roever. Bayesian random-effects meta-analysis using the bayesmeta R package. Journal of Statistical Software, 93(6):1-51, 2020. doi:10.18637/jss.v093.i06.

C. Roever, R. Bender, S. Dias, C.H. Schmid, H. Schmidli, S. Sturtz, S. Weber, T. Friede. On weakly informative prior distributions for the heterogeneity parameter in Bayesian random-effects meta-analysis. Research Synthesis Methods, 12(4):448-474, 2021. doi:10.1002/jrsm.1475.

C. Roever, T. Friede. Discrete approximation of a mixture distribution via restricted divergence. Journal of Computational and Graphical Statistics, 26(1):217-222, 2017. doi:10.1080/10618600.2016.1276840.

A. Gelman, J.B. Carlin, H.S. Stern, D.B. Rubin. Bayesian data analysis. Chapman & Hall / CRC, Boca Raton, 1997.

bayesmeta, escalc, model.matrix, CrinsEtAl2014, RobergeEtAl2017.

### Examples

## Not run:
######################################################################
# (1)  A simple example with two groups of studies

data("CrinsEtAl2014")
# compute effect measures (log-OR):
crins.es <- escalc(measure="OR",
ai=exp.AR.events,  n1i=exp.total,
ci=cont.AR.events, n2i=cont.total,
slab=publication, data=CrinsEtAl2014)
# show data:
crins.es[,c("publication", "IL2RA", "exp.AR.events", "exp.total",
"cont.AR.events", "cont.total", "yi", "vi")]

# specify regressor matrix:
X <- cbind("bas"=as.numeric(crins.es$IL2RA=="basiliximab"), "dac"=as.numeric(crins.es$IL2RA=="daclizumab"))
print(X)
print(cbind(crins.es[,c("publication", "IL2RA")], X))
# NB: regressor matrix specifies individual indicator covariates
#     for studies with "basiliximab" and "daclizumab" treatment.

# perform regression:
bmr01 <- bmr(y=crins.es$yi, sigma=sqrt(crins.es$vi),
labels=crins.es$publication, X=X) # alternatively, one may simply supply the "escalc" object # (yields identical results): bmr01 <- bmr(crins.es, X=X) # show results: bmr01 bmr01$summary
plot(bmr01)
pairs(bmr01)

# NOTE: there are many ways to set up the regressor matrix "X"
# (also affecting the interpretation of the involved parameters).
# See the above specification and check out the following alternatives:
X <- cbind("bas"=1, "offset.dac"=c(1,0,1,0,0,0))
X <- cbind("intercept"=1, "offset"=0.5*c(1,-1,1,-1,-1,-1))
# One may also use the "model.matrix()" function
# to specify regressor matrices via the "formula" interface; e.g.:
X <- model.matrix( ~ IL2RA, data=crins.es)
X <- model.matrix( ~ IL2RA - 1, data=crins.es)

######################################################################
# (2)  A simple example reproducing a "bayesmeta" analysis:

data("CrinsEtAl2014")
crins.es <- escalc(measure="OR",
ai=exp.AR.events,  n1i=exp.total,
ci=cont.AR.events, n2i=cont.total,
slab=publication, data=CrinsEtAl2014)

# a "simple" meta-analysis:
bma02 <- bayesmeta(crins.es,
tau.prior=function(t){dhalfnormal(t, scale=0.5)},
mu.prior.mean=0, mu.prior.sd=4)

# the equivalent "intercept-only" meta-regression:
bmr02 <- bmr(crins.es,
tau.prior=function(t){dhalfnormal(t, scale=0.5)},
beta.prior.mean=0, beta.prior.sd=4)
# the corresponding (default) regressor matrix:
bmr02$X # compare computation time and numbers of bins used internally: cbind("seconds" = c("bayesmeta" = unname(bma02$init.time),
"bmr"       = unname(bmr02$init.time)), "bins" = c("bayesmeta" = nrow(bma02$support),
"bmr"       = nrow(bmr02$support$tau)))

# compare heterogeneity estimates:
rbind("bayesmeta"=bma02$summary[,1], "bmr"=bmr02$summary[,1])

# compare effect estimates:
rbind("bayesmeta"=bma02$summary[,2], "bmr"=bmr02$summary[,2])

######################################################################
# (3)  An example with binary as well as continuous covariables:

data("RobergeEtAl2017")
str(RobergeEtAl2017)
?RobergeEtAl2017

# compute effect sizes (log odds ratios) from count data:
es.pe  <- escalc(measure="OR",
ai=asp.PE.events,  n1i=asp.PE.total,
ci=cont.PE.events, n2i=cont.PE.total,
slab=study, data=RobergeEtAl2017,
subset=complete.cases(RobergeEtAl2017[,7:10]))

# show "bubble plot" (bubble sizes are
# inversely proportional to standard errors):
plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi), col=c("blue","red")[as.numeric(es.pe$onset)],
xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)

# set up regressor matrix:
# (individual intercepts and slopes for two subgroups):
X <- model.matrix(~ -1 + onset + onset:dose, data=es.pe)
colnames(X) <- c("intEarly", "intLate", "slopeEarly", "slopeLate")
# check out regressor matrix (and compare to original data):
print(X)

# perform regression:
bmr03 <- bmr(es.pe, X=X)
bmr03$summary # derive predictions from the model; # specify corresponding "regressor matrices": newx.early <- cbind(1, 0, seq(50, 150, by=5), 0) newx.late <- cbind(0, 1, 0, seq(50, 150, by=5)) # (note: columns correspond to "beta" parameters) # compute predicted medians and 95 percent intervals: pred.early <- cbind("median"=bmr03$qpred(0.5, x=newx.early),
bmr03$pred.interval(x=newx.early)) pred.late <- cbind("median"=bmr03$qpred(0.5, x=newx.late),
bmr03$pred.interval(x=newx.late)) # draw "bubble plot": plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
col=c("blue","red")[as.numeric(es.pe\$onset)],
xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
# add predictions to bubble plot:
matlines(newx.early[,3], pred.early, col="blue", lty=c(1,2,2))
matlines(newx.late[,4], pred.late, col="red", lty=c(1,2,2))

## End(Not run)


[Package bayesmeta version 3.4 Index]