| 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,
rel.tol.integrate = 2^16*.Machine$double.eps,
abs.tol.integrate = 0.0,
tol.uniroot = rel.tol.integrate, ...)
## S3 method for class 'escalc'
bmr(y, labels = NULL, ...)
Arguments
y |
vector of estimates, or an |
sigma |
vector of standard errors associated with |
labels |
(optional) a vector of labels corresponding to |
X |
(optional) the regressor matrix for the regression. |
tau.prior |
a
The default is |
beta.prior.mean, beta.prior.sd, beta.prior.cov |
the mean and standard deviations, or covariance of the normal prior
distribution for the effects |
interval.type |
the type of (credible, prediction, shrinkage) interval to be
returned by default; either |
delta, epsilon |
the parameters specifying the desired accuracy for approximation of
the |
rel.tol.integrate, abs.tol.integrate, tol.uniroot |
the |
... |
other |
Details
The random-effects meta-regression model may be stated as
y_i|x_i,\beta,\sigma_i,\tau \;\sim\; \mathrm{Normal}(\beta_1 x_{i,1}
+ \beta_2 x_{i,2} + \ldots + \beta_d x_{i,d}, \;
\sigma_i^2 + \tau^2)
where the data (y, \sigma) enter as y_i, the
i-th estimate, that is associated with standard error
\sigma_i > 0, where i=1,...,k. In addition to
estimates and standard errors for the ith observation,
a set of covariables x_{i,j} with j=1,...,d are available
for each estimate y_i.
The model includes d+1 unknown parameters,
namely, the d coefficients (\beta_1,...,\beta_d), and the heterogeneity
\tau. Alternatively, the model may also be formulated via an
intermediate step as
y_i|\theta_i,\sigma_i \;\sim\; \mathrm{Normal}(\theta_i, \, \sigma_i^2),
\theta_i|\beta,x_i,\tau \;\sim\; \mathrm{Normal}(\beta_1 x_{i,1}
+ \ldots + \beta_d x_{i,d}, \; \tau^2),
where the \theta_i denote the trial-specific means
that are then measured through the estimate y_i with an
associated measurement uncertainty \sigma_i. The
\theta_i again differ from trial to trial (even for
identical covariable vectors x_i) and are
distributed around a mean of \beta_1 x_{i,1}
+ \ldots + \beta_d x_{i,d} with
standard deviation \tau.
It if often convenient to express the model in matrix notation, i.e.,
y|\theta,\sigma \;\sim\; \mathrm{Normal}(\theta, \,
\Sigma)
\theta|X,\beta,\tau \;\sim\; \mathrm{Normal}(X \beta, \,
\tau I)
where y, \sigma, \beta and \theta now denote
k-dimensional vectors, X is the (k \times d) regressor matrix, and \Sigma is a (k \times k) diagonal covariance matrix containing the
\sigma_i^2 values, while
I is the (k \times k) identity matrix. The
regressor matrix X plays a crucial role here, as the
‘X’ argument (with rows corresponding to studies, and
columns corresponding to covariables) is required to specify the exact
regression setup.
Meta-regression allows the consideration of (study-level) covariables in a meta-analysis. Quite often, these may also be indicator variables (‘zero/one’ variables) simply identifying subgroups of studies. See also the examples shown below.
Connection to the simple random-effects model
The meta-regression model is a generalisation of the ‘simple’
random-effects model that is implemented in the
bayesmeta() function. Meta-regression reduces to the
estimation of a single “intercept” term when the regressor
matrix (X) consists of a single column of
ones (which is also the default setting in case the ‘X’
argument is left unspecified). The single regression coefficient
\beta_1 then is equivalent to the \mu parameter
from the simple random effects model (see also the ‘Examples’
section below).
Specification of the regressor matrix
The actual
regression model is specified through the regressor matrix X,
which is supplied via the ‘X’ argument, and which
often may be specified in different ways. There usually is no unique
solution, and what serves the present purpose best then depends on
the context; see also the examples below. Sensible column names
should be specified for X, as these will subsequently
determine the labels for the associated parameters later on. Model
specification via the regressor matrix has the advantage of being
very explicit and transparent; if one prefers a
formula interface instead, a regressor matrix may
be generated via the ‘model.matrix()’
function.
Prior specification
Priors for \beta and \tau are assumed to factor into
into independent marginals p(\beta,\tau)=p(\beta)\times
p(\tau) and either (improper)
uniform or a normal priors may be specified for the regression coefficients
\beta.
For sensible prior choices for the heterogeneity parameter \tau,
see also Roever (2020), Roever et al. (2021) and the
‘bayesmeta()’ function's help.
Accessing posterior density functions, etc.
Within the bayesmeta() function, access to posterior
density, cumulative distribution function, quantile functtion,
random number generation and posterior inverval computation is
implemented via the $dposterior(), $dposterior(),
$pposterior(), $qposterior(), $rposterior()
and $post.interval() functions that are accessible as elements
in the returned list object. Prediction and shrinkage
estimation are available by setting additional arguments in the
above functions.
In the meta-regression context things get slightly more complicated,
as the \beta parameter may be of higher dimension. Hence, in the
bmr() function, the three different types of distributions
related to posterior distribution, prediction and
shrinkage are split up into three groups of
functions. For example, the posterior density is accessible via the
$dposterior() function, the predictive distribution via the
$dpredict() function, and the shrinkage estimates via the
$dshrink() function. Analogous functions are returned for
cumulative distribution, quantile function, etc.; see also the
‘Value’ section below.
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 |
X |
the regressor matrix. |
k |
number of data points (length of |
d |
number of coefficients (columns of |
labels |
vector of labels corresponding to |
variables |
variable names for the |
tau.prior |
the prior probability density function for |
tau.prior.proper |
a |
beta.prior |
a |
beta.prior.proper |
a |
dprior |
a |
likelihood |
a |
dposterior, pposterior, qposterior, rposterior, post.interval |
functions of |
dpredict, ppredict, qpredict, rpredict, pred.interval |
functions
of |
dshrink, pshrink, qshrink, rshrink, shrink.interval |
functions
of |
post.moments |
a |
pred.moments |
a |
shrink.moments |
a |
summary |
a |
interval.type |
the |
ML |
a |
MAP |
a |
theta |
a |
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 |
support |
a |
delta, epsilon |
the ‘ |
rel.tol.integrate, abs.tol.integrate, tol.uniroot |
the input
parameters determining the numerical accuracy of the internally used
|
call |
an object of class |
init.time |
the computation time (in seconds) used to generate
the |
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.
See Also
bayesmeta, escalc,
model.matrix, CrinsEtAl2014,
RobergeEtAl2017.
Examples
## Not run:
######################################################################
# (1) A simple example with two groups of studies
# load data:
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:
# load data:
data("RobergeEtAl2017")
str(RobergeEtAl2017)
head(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)