bmr: Bayesian random-effects meta-regression
In bayesmeta: Bayesian Random-Effects Meta-Analysis and Meta-Regression
bmr 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 escalc object.
sigma
vector of standard errors associated with y.
labels
(optional) a vector of labels corresponding to y and sigma.
X
(optional) the regressor matrix for the regression.
tau.prior
a function returning the prior density for the heterogeneity
parameter (\tau) or a character string
specifying one of the default ‘non-informative’
priors; possible choices for the latter case are:
-
"uniform": a uniform prior in \tau
-
"sqrt": a uniform prior in \sqrt{\tau}
-
"Jeffreys": the Jeffreys prior for \tau
-
"BergerDeely": the prior due to Berger and Deely (1988)
-
"conventional": the conventional prior
-
"DuMouchel": the DuMouchel prior
-
"shrinkage": the ‘uniform shrinkage’ prior
-
"I2": a uniform prior on the ‘relative heterogeneity’ I^2
The default is "uniform" (which should be used with
caution). The above priors are described in some more detail in the
bayesmeta() help.
beta.prior.mean, beta.prior.sd, beta.prior.cov
the mean and standard deviations, or covariance of the normal prior
distribution for the effects \beta. If unspecified, an
(improper) uniform prior is used.
interval.type
the type of (credible, prediction, shrinkage) interval to be
returned by default; either "shortest" for shortest
intervals, or "central" for central, equal-tailed intervals.
delta, epsilon
the parameters specifying the desired accuracy for approximation of
the \beta posterior(s), and with that determining the number of
\tau support points being used internally. Smaller values
imply greater accuracy and greater computational burden (Roever and
Friede, 2017).
rel.tol.integrate, abs.tol.integrate, tol.uniroot
the rel.tol, abs.tol and tol
‘accuracy’ arguments that are passed to
the integrate() or uniroot() functions
for internal numerical integration or root finding
(see also the help there).
...
other bmr arguments.
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
codeformula 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 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.
\Sexpr[results=rd]{tools:::Rd_expr_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.
\Sexpr[results=rd]{tools:::Rd_expr_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.
\Sexpr[results=rd]{tools:::Rd_expr_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.
\Sexpr[results=rd]{tools:::Rd_expr_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)
bayesmeta documentation built on July 9, 2023, 5:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| bmr | 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
codeformula 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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.