# Bayesian random-effects meta-regression 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. ```r 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)y[i] | x[i], beta, sigma[i], tau ~Normal(beta[1]x[i,1] + beta[2]x[i,2] + ... + 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 $i$th 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),y[i]|theta[i],sigma[i] ~ 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),theta[i]|beta,x[i],tau ~ Normal(beta[1]x[i,1]+...+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]+...+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)y | theta, sigma ~Normal(theta, Sigma) $$ $$ \theta|X,\beta,\tau \;\sim\; \mathrm{Normal}(X \beta, \,\tau I)theta | X, beta, tau ~Normal(X beta, tau I) $$ where $y$, $\sigma$, $\beta$ and $\theta$ now denote $k$-dimensional vectors, $X$ is the (c("$k x\n$", "$ d$")) **regressor matrix**, and $\Sigma$ is a (c("$k x\n$", "$ k$")) diagonal covariance matrix containing the $sigma[i]^2$ values, while $I$ is the ($k x 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) * 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. ## Returns 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** ($t(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,...,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. ## 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. 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. 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. 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. 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. ## Author(s) Christian Roever christian.roever@med.uni-goettingen.de ## See Also `bayesmeta`, `escalc`, `model.matrix`, `CrinsEtAl2014`, `RobergeEtAl2017`. ## Examples ```r ## 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) ```

bmr function