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

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)

[Package bayesmeta version 3.4 Index]