Package 'bayesmeta'

Bayesian Random-Effects Meta-Analysis and Meta-Regression

Description

A collection of functions allowing to derive the posterior distribution of the model parameters in random-effects meta-analysis or meta-regression, and providing functionality to evaluate joint and marginal posterior probability distributions, predictive distributions, shrinkage effects, posterior predictive p-values, etc.

Details

Package:	bayesmeta
Type:	Package
Version:	3.4
Date:	2024-02-15
License:	GPL (>=2)

The main functionality is provided by the bayesmeta() and bmr() function. It takes the data (estimates and associated standard errors) and prior information (effect and heterogeneity priors), and returns an object containing functions that allow to derive posterior quantities like joint or marginal densities, quantiles, etc. The bmr() function extends the approach to meta-regression by allowing to specify covariables (moderators) in addition.

Author(s)

Christian Roever <[email protected]>

References

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

See Also

forestplot.bayesmeta, plot.bayesmeta, bmr.

Examples

# example data by Snedecor and Cochran:
data("SnedecorCochran")

## Not run: 
# analysis using improper uniform prior
# (may take a few seconds to compute!):
bma <- bayesmeta(y=SnedecorCochran[,"mean"],
                 sigma=sqrt(SnedecorCochran[,"var"]),
                 label=SnedecorCochran[,"no"])

# show some summary statistics:
bma

# show a bit more details:
summary(bma)

# show a forest plot:
forestplot(bma)

# show some more plots:
plot(bma)

## End(Not run)

---

Ankylosing spondylitis example data

Description

Numbers of cases (patients) and events (responders) in the placebo control groups of eight studies.

Usage

data("BaetenEtAl2013")

Format

The data frame contains the following columns:

study	character	study label
year	numeric	publication year
events	numeric	number of responders
total	numeric	total number of patients

Details

A study was conducted in order to investigate a novel treatment in ankylosing spondylitis (Baeten et al., 2013). The primary endpoint related to treatment response. In order to formulate an informative prior distribution for the response rate to be expected in the control group of the new study, a systematic review of previous studies was consulted (McLeod et al., 2007), and, after a meta-analysis of the estimated response probabilities, the predictive distribution for the new study's response probability was derived. The predictive distribution here constitutes the meta-analytic-predictive (MAP) prior distribution (Schmidli et al., 2014). The data set contains the relevant data from the eight “historical” studies' placebo groups.

Note that the original analysis (Baeten et al., 2013) involved a binomial model, and the resulting posterior predictive distribution was eventually approximated by a mixture of beta distributions.

Source

D. Baeten et al. Anti-interleukin-17A monoclonal antibody secukinumab in treatment of ankylosing spondylitis: a randomised, double-blind, placebo-controlled trial. The Lancet, 382(9906):1705-1713, 2013. doi:10.1016/S0140-6736(13)61134-4.

References

C. McLeod et al. Adalimumab, etanercept, and infliximab for the treatment of ankylosing spondylitis: a systematic review and economic evaluation. Health Technology Assessment, 11(28), 2007. doi:10.3310/hta11280.

H. Schmidli, S. Gsteiger, S. Roychoudhury, A. O'Hagan, D. Spiegelhalter, B. Neuenschwander. Robust meta-analytic-predictive priors in clinical trials with historical control information. Biometrics, 70(4):1023-1032, 2014. doi:10.1111/biom.12242.

See Also

uisd, ess.

Examples

# load data:
data("BaetenEtAl2013")

# show data:
BaetenEtAl2013

## Not run: 
# compute effect sizes (logarithmic odds) from the count data:
as <- escalc(xi=events, ni=total, slab=study,
             measure="PLO", data=BaetenEtAl2013)

# compute the unit information standard deviation (UISD):
uisd(as)

# perform meta-analysis
# (using uniform priors for effect and heterogeneity):
bm <- bayesmeta(as)

# show results (log-odds):
forestplot(bm, xlab="log-odds", zero=NA)
# show results (odds):
forestplot(bm, exponentiate=TRUE, xlab="odds", zero=NA)

# show posterior predictive distribution --
# in terms of log-odds:
bm$summary[,"theta"]
logodds <- bm$summary[c(2,5,6), "theta"]
logodds
# in terms of odds:
exp(logodds)
# in terms of probabilities:
(exp(logodds) / (exp(logodds) + 1))

# illustrate MAP prior density:
x <- seq(-3, 1, by=0.01)
plot(x, bm$dposterior(theta=x, predict=TRUE), type="l",
     xlab="log-odds (response)", ylab="posterior predictive density")
abline(h=0, v=0, col="grey")

## End(Not run)

---

Bayesian random-effects meta-analysis

Description

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

Usage

bayesmeta(y, ...)
  ## Default S3 method:
bayesmeta(y, sigma, labels = names(y),
          tau.prior = "uniform",
          mu.prior = c("mean"=NA,"sd"=NA),
          mu.prior.mean = mu.prior[1], mu.prior.sd = mu.prior[2],
          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'
bayesmeta(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.
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; see remarks below). The above priors are described in some more detail below.
mu.prior	the mean and standard deviation of the normal prior distribution for the effect $\mu$. If unspecified, an (improper) uniform prior is used.
mu.prior.mean, mu.prior.sd	alternative parameters to specify the prior distribution for the effect $\mu$ (see above).
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 $\mu$ 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 bayesmeta arguments.

Details

The common random-effects meta-analysis model may be stated as

$y_i|\mu,\sigma_i,\tau \;\sim\; \mathrm{Normal}(\mu, \, \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$, and $i=1,...,k$. The model includes two unknown parameters, namely the (mean) effect $\mu$, 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|\mu,\tau \;\sim\; \mathrm{Normal}(\mu, \, \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 and are distributed around a common mean $\mu$ with standard deviation $\tau$.

The bayesmeta() function utilizes the fact that the joint posterior distribution $p(\mu, \tau | y, \sigma)$ may be partly analytically integrated out to determine the integrals necessary for coherent Bayesian inference on one or both of the parameters.

As long as we assume that the prior probability distribution may be factored into independent marginals $p(\mu,\tau)=p(\mu)\times p(\tau)$ and either an (improper) uniform or a normal prior is used for the effect $\mu$, the joint likelihood $p(y|\mu,\tau)$ may be analytically marginalized over $\mu$, yielding the marginal likelihood function $p(y|\tau)$. Using any prior $p(\tau)$ for the heterogeneity, the 1-dimensional marginal posterior $p(\tau|y) \propto p(y|\tau)\times p(\tau)$ may then be treated numerically. As the conditional posterior $p(\mu|\tau,y)$ is a normal distribution, inference on the remaining joint ($p(\mu,\tau|y)$) or marginal ($p(\mu|y)$) posterior may be approached numerically (using the DIRECT algorithm) as well. Accuracy of the computations is determined by the delta (maximum divergence $\delta$) and epsilon (tail probability $\epsilon$) parameters (Roever and Friede, 2017).

What constitutes a sensible prior for both $\mu$ and $\tau$ depends (as usual) very much on the context. Potential candidates for informative (or weakly informative) heterogeneity ($\tau$) priors may include the half-normal, half-Student-$t$, half-Cauchy, exponential, or Lomax distributions; we recommend the use of heavy-tailed prior distributions if in case of prior/data conflict the prior is supposed to be discounted (O'Hagan and Pericchi, 2012). A sensible informative prior might also be a log-normal or a scaled inverse $\chi^2$ distribution. One might argue that an uninformative prior for $\tau$ may be uniform or monotonically decreasing in $\tau$. Another option is to use the Jeffreys prior (see the tau.prior argument above). The Jeffreys prior implemented here is the variant also described by Tibshirani (1989) that results from fixing the location parameter ($\mu$) and considering the Fisher information element corresponding to the heterogeneity $\tau$ only. This prior also constitutes the Berger/Bernardo reference prior for the present problem (Bodnar et al., 2016). The uniform shrinkage and DuMouchel priors are described in Spiegelhalter et al. (2004, Sec. 5.7.3). The procedure is able to handle improper priors (and does so by default), but of course the usual care must be taken here, as the resulting posterior may or may not be a proper probability distribution.

Note that a wide range of different types of endpoints may be treated on a continuous scale after an appropriate transformation; for example, count data may be handled by considering corresponding logarithmic odds ratios. Many such transformations are implemented in the metafor package's escalc function and may be directly used as an input to the bayesmeta() function; see also the example below. Alternatively, the compute.es package also provides a range of effect sizes to be computed from different data types.

The bayesmeta() function eventually generates some basic summary statistics, but most importantly it provides an object containing a range of functions allowing to evaluate posterior distributions; this includes joint and marginal posteriors, prior and likelihood, predictive distributions, densities, cumulative distributions and quantile functions. For more details see also the documentation and examples below. Use of the individual argument allows to access posteriors of study-specific (shrinkage-) effects ($\theta_i$). The predict argument may be used to access the predictive distribution of a future study's effect ($\theta_{k+1}$), facilitating a meta-analytic-predictive (MAP) approach (Neuenschwander et al., 2010).

Prior specification details

When specifying the tau.prior argument as a character string (and not as a prior density function), then the actual prior probability density functions corresponding to the possible choices of the tau.prior argument are given by:

• "uniform" - the (improper) uniform prior in $\tau$:

  $p(\tau) \;\propto\; 1$

• "sqrt" - the (improper) uniform prior in $\sqrt{\tau}$:

  $p(\tau) \;\propto\; \tau^{-1/2} \;=\; \frac{1}{\sqrt{\tau}}$

• "Jeffreys" - Tibshirani's noninformative prior, a variation of the (‘general’ or ‘overall’) Jeffreys prior, which here also constitutes the Berger/Bernardo reference prior for $\tau$:

  $p(\tau) \;\propto\; \sqrt{\sum_{i=1}^k\Bigl(\frac{\tau}{\sigma_i^2+\tau^2}\Bigr)^2}$

  This is also an improper prior whose density does not integrate to 1. This prior results from applying Jeffreys' general rule (Kass and Wasserman, 1996), and in particular considering that location parameters (here: the effect $\mu$) should be treated separately (Roever, 2020).

• "overallJeffreys" - the ‘general’ or ‘overall’ form of the Jeffreys prior; this is derived based on the Fisher information terms corresponding to both the effect ($\mu$) and heterogeneity ($\tau$). The resulting (improper) prior density is

  $p(\tau) \;\propto\; \sqrt{\sum_{i=1}^k\frac{1}{\sigma_i^2+\tau^2} \; \sum_{i=1}^k\Bigl(\frac{\tau}{\sigma_i^2+\tau^2}\Bigr)^2}$

  Use of this specification is generally not recommended; see, e.g., Jeffreys (1946), Jeffreys (1961, Sec. III.3.10), Berger (1985, Sec. 3.3.3) and Kass and Wasserman (1996, Sec. 2.2). Since the effect $\mu$ constitutes a location parameter here, it should be treated separately (Roever, 2020).

• "BergerDeely" - the (improper) Berger/Deely prior:

  $p(\tau) \;\propto\; \prod_{i=1}^k \Bigl(\frac{\tau}{\sigma_i^2+\tau^2}\Bigr)^{1/k}$

  This is a variation of the above Jeffreys prior, and both are equal in case all standard errors ($\sigma_i$) are the same.

• "conventional" - the (proper) conventional prior:

  $p(\tau) \;\propto\; \prod_{i=1}^k \biggl(\frac{\tau}{(\sigma_i^2+\tau^2)^{3/2}}\biggr)^{1/k}$

  This is a proper variation of the above Berger/Deely prior intended especially for testing and model selection purposes.

• "DuMouchel" - the (proper) DuMouchel prior:

  $p(\tau) \;=\; \frac{s_0}{(s_0+\tau)^2}$

  where $s_0=\sqrt{s_0^2}$ and $s_0^2$ again is the harmonic mean of the standard errors (as above).

• "shrinkage" - the (proper) uniform shrinkage prior:

  $p(\tau) \;=\; \frac{2 s_0^2 \tau}{(s_0^2+\tau^2)^2}$

  where $s_0^2=\frac{k}{\sum_{i=1}^k \sigma_i^{-2}}$ is the harmonic mean of the squared standard errors $\sigma_i^2$.

• "I2" - the (proper) uniform prior in $I^2$:

  $p(\tau) \;=\; \frac{2 \hat{\sigma}^2 \tau}{(\hat{\sigma}^2 + \tau^2)^2}$

  where $\hat{\sigma}^2 = \frac{(k-1)\; \sum_{i=1}^k\sigma_i^{-2}}{\bigl(\sum_{i=1}^k\sigma_i^{-2}\bigr)^2 - \sum_{i=1}^k\sigma_i^{-4}}$. This prior is similar to the uniform shrinkage prior, except for the use of $\hat{\sigma}^2$ instead of $s_0^2$.

For more details on the above priors, see Roever (2020). For more general information especially on (weakly) informative heterogeneity prior distributions, see Roever et al. (2021). Regarding empirically motivated informative heterogeneity priors see also the TurnerEtAlPrior() and RhodesEtAlPrior() functions, or Roever et al. (2023) and Lilienthal et al. (2023).

Credible intervals

Credible intervals (as well as prediction and shrinkage intervals) may be determined in different ways. By default, shortest intervals are returned, which for unimodal posteriors (the usual case) is equivalent to the highest posterior density region (Gelman et al., 1997, Sec. 2.3). Alternatively, central (equal-tailed) intervals may also be derived. The default behaviour may be controlled via the interval.type argument, or also by using the method argument with each individual call of the $post.interval() function (see below). A third option, although not available for prediction or shrinkage intervals, and hence not as an overall default choice, but only for the $post.interval() function, is to determine the evidentiary credible interval, which has the advantage of being parameterization invariant (Shalloway, 2014).

Bayes factors

Bayes factors (Kass and Raftery, 1995) for the two hypotheses of $\tau=0$ and $\mu=0$ are provided in the $bayesfactor element; low or high values here constitute evidence against or in favour of the hypotheses, respectively. Bayes factors are based on marginal likelihoods and can only be computed if the priors for heterogeneity and effect are proper. Bayes factors for other hypotheses can be computed using the marginal likelihood (as provided through the $marginal element) and the $likelihood() function. Bayes factors must be interpreted with very much caution, as they are susceptible to Lindley's paradox (Lindley, 1957), which especially implies that variations of the prior specifications that have only minuscule effects on the posterior distribution may have a substantial impact on Bayes factors (via the marginal likelihood). For more details on the problems and challenges related to Bayes factors see also Gelman et al. (1997, Sec. 7.4).

Besides the ‘actual’ Bayes factors, minimum Bayes factors are also provided (Spiegelhalter et al., 2004; Sec. 4.4). The minimum Bayes factor for the hypothesis of $\mu=0$ constitutes the minimum across all possible priors for $\mu$ and hence gives a measure of how much (or how little) evidence against the hypothesis is provided by the data at most. It is independent of the particular effect prior used in the analysis, but still dependent on the heterogeneity prior. Analogously, the same is true for the heterogeneity's minimum Bayes factor. A minimum Bayes factor can also be computed when only one of the priors is proper.

Numerical accuracy

Accuracy of the numerical results is determined by four parameters, namely, the accuracy of numerical integration as specified through the rel.tol.integrate and abs.tol.integrate arguments (which are internally passed on to the integrate function), and the accuracy of the grid approximation used for integrating out the heterogeneity as specified through the delta and epsilon arguments (Roever and Friede, 2017). As these may also heavily impact on the computation time, be careful when changing these from their default values. You can monitor the effect of different settings by checking the number and range of support points returned in the $support element.

Study weights

Conditional on a given $\tau$ value, the posterior expectations of the overall effect ($\mu$) as well as the shrinkage estimates ($\theta_i$) result as convex combinations of the estimates $y_i$. The weights associated with each estimate $y_i$ are commonly quoted in frequentist meta-analysis results in order to quantify (arguably somewhat heuristically) each study's contribution to the overall estimates, often in terms of percentages.

In a Bayesian meta-analysis, these numbers to not immediately arise, since the heterogeneity is marginalized over. However, due to linearity, the posterior mean effects may still be expressed in terms of linear combinations of initial estimates $y_i$, with weights now given by the posterior mean weights, marginalized over the heterogeneity $\tau$ (Roever and Friede, 2020). The posterior mean weights are returned in the $weights and $weights.theta elements, for the overall effect $\mu$ as well as for the shrinkage estimates $\theta_i$.

Value

A list of class bayesmeta containing the following elements:

y	vector of estimates (the input data).
sigma	vector of standard errors corresponding to y (input data).
labels	vector of labels corresponding to y and sigma.
k	number of data points (in y).
tau.prior	the prior probability density function for $\tau$.
mu.prior.mean	the prior mean of $\mu$.
mu.prior.sd	the prior standard deviation of $\mu$.
dprior	a function(tau=NA, mu=NA, log=FALSE) of two parameters, tau and/or mu, returning either the joint or marginal prior probability density, depending on which parameter(s) is/are provided.
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).
likelihood	a function(tau=NA, mu=NA, log=FALSE) of two parameters, tau and/or mu, returning either the joint or marginal likelihood, depending on which parameter(s) is/are provided.
dposterior	a function(tau=NA, mu=NA, theta=mu, log=FALSE, predict=FALSE, individual=FALSE) of two parameters, tau and/or mu, returning either the joint or marginal posterior probability density, depending on which parameter(s) is/are provided. Using the argument predict=TRUE yields the posterior predictive distribution for $\theta$. Using the individual argument, you can request individual effects' ($\theta_i$) posterior distributions. May be an integer number (1,...,k) giving the index, or a character string giving the label.
pposterior	a function(tau=NA, mu=NA, theta=mu, predict=FALSE, individual=FALSE) of one parameter (either tau or mu) returning the corresponding marginal posterior cumulative distribution function. Using the argument predict=TRUE yields the posterior predictive distribution for $\theta$. Using the individual argument, you can request individual effects' ($\theta_i$) posterior distributions. May be an integer number (1,...,k) giving the index, or a character string giving the label.
qposterior	a function(tau.p=NA, mu.p=NA, theta.p=mu.p, predict=FALSE, individual=FALSE) of one parameter (either tau.p or mu.p) returning the corresponding marginal posterior quantile function. Using the argument predict=TRUE yields the posterior predictive distribution for $\theta$. Using the individual argument, you can request individual effects' ($\theta_i$) posterior distributions. May be an integer number (1,...,k) giving the index, or a character string giving the label.
rposterior	a function(n=1, predict=FALSE, individual=FALSE, tau.sample=TRUE) generating n independent random draws from the (2-dimensional) posterior distribution. Using the argument predict=TRUE yields the posterior predictive distribution for $\theta$. Using the individual argument, you can request individual effects' ($\theta_i$) posterior distributions. May be an integer number (1,...,k) giving the index, or a character string giving the label. In general, this via the inversion method, so it is rather slow. However, if one is not interested in sampling the heterogeneity parameter ($\tau$), using ‘tau.sample=FALSE’ will speed up the function substantially.
post.interval	a function(tau.level=NA, mu.level=NA, theta.level=mu.level, method=c("shortest","central","evidentiary"), predict=FALSE, individual=FALSE) returning a credible interval, depending on which of the two parameters is provided (either tau.level or mu.level). The additional parameter method may be used to specify the desired type of interval: method = "shortest" returns the shortest interval, method = "central" returns a central interval, and method = "evidentiary" returns an evidentiary interval (Shalloway, 2014); the former is the default option. Using the argument predict=TRUE yields a posterior predictive interval for $\theta$. Using the individual argument, you can request individual effects' ($\theta_i$) posterior distributions. May be an integer number (1,...,k) giving the index, or a character string giving the label.
cond.moment	a function(tau, predict=FALSE, individual=FALSE, simplify=TRUE) returning conditional moments (mean and standard deviation) of $\mu$ as a function of $\tau$. Using the argument predict=TRUE yields (conditional) posterior predictive moments for $\theta$. Using the individual argument, you can request individual effects' ($\theta_i$) posterior distributions. May be an integer number (1,...,k) giving the index, or a character string giving the label.
I2	a function(tau) returning the ‘relative’ heterogeneity $I^2$ 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 $\mu$, and of the posterior predictive distribution of $\theta$.
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,\mu)$.
MAP	a matrix giving joint and marginal maximum-a-posteriori estimates of $(\tau,\mu)$.
theta	a matrix giving the ‘shrinkage estimates’, i.e, summary statistics of the trial-specific means $\theta_i$.
weights	a vector giving the posterior expected inverse-variance weights for each study (and for the effect prior mean, if the effect prior was proper).
weights.theta	a matrix whose columns give the posterior expected weights of each study (and of the effect prior mean, if the effect prior was proper) for all shrinkage estimates.
marginal.likelihood	the marginal likelihood of the data (this number is only computed if proper effect and heterogeneity priors are specified).
bayesfactor	Bayes factors and minimum Bayes factors for the two hypotheses of $\tau=0$ and $\mu=0$. These depend on the marginal likelihood and hence can only be computed if proper effect and/or heterogeneity priors are specified; see also remark above.
support	a matrix giving the $\tau$ support points used internally in the grid approximation, along with their associated weights, conditional mean and standard deviation of $\mu$, and the standard deviation of the (conditional) predictive distribution of $\theta$.
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 bayesmeta object.
init.time	the computation time (in seconds) used to generate the bayesmeta object.

Author(s)

Christian Roever [email protected]

References

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.

C. Roever, T. Friede. Bounds for the weight of external data in shrinkage estimation. Biometrical Journal, 65(5):1131-1143, 2021. doi:10.1002/bimj.202000227.

J. O. Berger. Statistical Decision Theory and Bayesian Analysis. 2nd edition. Springer-Verlag, 1985. doi:10.1007/978-1-4757-4286-2.

J.O. Berger, J. Deely. A Bayesian approach to ranking and selection of related means with alternatives to analysis-of-variance methodology. Journal of the American Statistical Association, 83(402):364-373, 1988. doi:10.1080/01621459.1988.10478606.

O. Bodnar, A. Link, C. Elster. Objective Bayesian inference for a generalized marginal random effects model. Bayesian Analysis, 11(1):25-45, 2016. doi:10.1214/14-BA933.

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

A. Gelman. Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3):515-534, 2006. doi:10.1214/06-BA117A.

J. Hartung, G. Knapp, B.K. Sinha. Statistical meta-analysis with applications. Wiley, Hoboken, 2008.

L.V. Hedges, I. Olkin. Statistical methods for meta-analysis. Academic Press, San Diego, 1985

H. Jeffreys. An invariant form for the prior probability in estimation problems. Proceedings of the Royal Society of london, Series A, 186(1007):453-462, 1946. doi:10.1098/rspa.1946.0056.

H. Jeffreys. Theory of Probability. 3rd edition. Clarendon Press, Oxford, 1961.

A.E. Kass, R.E. Raftery. Bayes factors. Journal of the American Statistical Association, 90(430):773-795, 1995. doi:10.2307/2291091.

A.E. Kass, L. Wasserman. The selection of prior distributions by formal rules. Journal of the American Statistical Association. 91(453):1243-1370, 1996. doi:10.1080/01621459.1996.10477003.

J. Lilienthal, S. Sturtz, C. Schuermann, M. Maiworm, C. Roever, T. Friede, R. Bender. Bayesian random-effects meta-analysis with empirical heterogeneity priors for application in health technology assessment with very few studies. Research Synthesis Methods, 2023. doi:10.1002/jrsm.1685.

D.V. Lindley. A statistical paradox. Biometrika, 44(1/2):187-192, 1957. doi:10.1093/biomet/44.1-2.187.

B. Neuenschwander, G. Capkun-Niggli, M. Branson, D.J. Spiegelhalter. Summarizing historical information on controls in clinical trials. Trials, 7(1):5-18, 2010. doi:10.1177/1740774509356002.

A. O'Hagan, L. Pericchi. Bayesian heavy-tailed models and conflict resolution: A review. Brazilian Journal of Probability and Statistics, 26(4):372-401, 2012. doi:10.1214/11-BJPS164.

C. Roever, S. Sturtz, J. Lilienthal, R. Bender, T. Friede. Summarizing empirical information on between-study heterogeneity for Bayesian random-effects meta-analysis. Statistics in Medicine, 42(14):2439-2454, 2023. doi:10.1002/sim.9731.

S. Shalloway. The evidentiary credible region. Bayesian Analysis, 9(4):909-922, 2014. doi:10.1214/14-BA883.

D.J. Spiegelhalter, K.R. Abrams, J.P.Myles. Bayesian approaches to clinical trials and health-care evaluation. Wiley & Sons, 2004.

R. Tibshirani. Noninformative priors for one parameter of many. Biometrika, 76(3):604-608, 1989. doi:10.1093/biomet/76.3.604.

W. Viechtbauer. Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3):1-48, 2010. doi:10.18637/jss.v036.i03.

See Also

forestplot.bayesmeta, plot.bayesmeta, escalc, bmr, compute.es.

Examples

########################################
# example data by Snedecor and Cochran:
data("SnedecorCochran")

## Not run: 
# analysis using improper uniform prior
# (may take a few seconds to compute!):
ma01 <- bayesmeta(y=SnedecorCochran[,"mean"], sigma=sqrt(SnedecorCochran[,"var"]),
                  label=SnedecorCochran[,"no"])

# analysis using an informative prior
# (normal for mu and half-Cauchy for tau (scale=10))
# (may take a few seconds to compute!):
ma02 <- bayesmeta(y=SnedecorCochran[,"mean"], sigma=sqrt(SnedecorCochran[,"var"]),
                  label=SnedecorCochran[,"no"],
                  mu.prior.mean=50, mu.prior.sd=50,
                  tau.prior=function(x){return(dhalfcauchy(x, scale=10))})

# show some summary statistics:
print(ma01)
summary(ma01)

# show some plots:
forestplot(ma01)
plot(ma01)

# compare resulting marginal densities;
# the effect parameter (mu):
mu <- seq(30, 80, le=200)
plot(mu, ma02$dposterior(mu=mu), type="l", lty="dashed",
     xlab=expression("effect "*mu),
     ylab=expression("marginal posterior density"),
     main="Snedecor/Cochran example")
lines(mu, ma01$dposterior(mu=mu), lty="solid")

# the heterogeneity parameter (tau):
tau <- seq(0, 50, le=200)
plot(tau, ma02$dposterior(tau=tau), type="l", lty="dashed",
     xlab=expression("heterogeneity "*tau),
     ylab=expression("marginal posterior density"),
     main="Snedecor/Cochran example")
lines(tau, ma01$dposterior(tau=tau), lty="solid")

# compute posterior median relative heterogeneity I-squared:
ma01$I2(tau=ma01$summary["median","tau"])

# determine 90 percent upper limits on the heterogeneity tau:
ma01$qposterior(tau=0.90)
ma02$qposterior(tau=0.90)
# determine shortest 90 percent credible interval for tau:
ma01$post.interval(tau.level=0.9, method="shortest")
## End(Not run)


#####################################
# example data by Sidik and Jonkman:
data("SidikJonkman2007")
# add log-odds-ratios and corresponding standard errors:
sj <- SidikJonkman2007
sj <- cbind(sj, "log.or"=log(sj[,"lihr.events"])-log(sj[,"lihr.cases"]-sj[,"lihr.events"])
                             -log(sj[,"oihr.events"])+log(sj[,"oihr.cases"]-sj[,"oihr.events"]),
                "log.or.se"=sqrt(1/sj[,"lihr.events"] + 1/(sj[,"lihr.cases"]-sj[,"lihr.events"])
                                 + 1/sj[,"oihr.events"] + 1/(sj[,"oihr.cases"]-sj[,"oihr.events"])))

## Not run: 
# analysis using weakly informative half-normal prior
# (may take a few seconds to compute!):
ma03a <- bayesmeta(y=sj[,"log.or"], sigma=sj[,"log.or.se"],
                   label=sj[,"id.sj"],
                   tau.prior=function(t){dhalfnormal(t,scale=1)})

# alternatively: may utilize "metafor" package's "escalc()" function
# to compute log-ORs and standard errors:
require("metafor")
es <- escalc(measure="OR",
             ai=lihr.events, n1i=lihr.cases,
             ci=oihr.events, n2i=oihr.cases,
             slab=id, data=SidikJonkman2007)
# apply "bayesmeta()" function directly to "escalc" object:
ma03b <- bayesmeta(es, tau.prior=function(t){dhalfnormal(t,scale=1)})
# "ma03a" and "ma03b" should be identical:
print(ma03a$summary)
print(ma03b$summary)
# compare to metafor's (frequentist) random-effects meta-analysis:
rma03a <- rma.uni(es)
rma03b <- rma.uni(es, method="EB", knha=TRUE)
# compare mu estimates (estimate and confidence interval):
plot(ma03b, which=3)
abline(v=c(rma03a$b, rma03a$ci.lb, rma03a$ci.ub), col="red", lty=c(1,2,2))
abline(v=c(rma03b$b, rma03b$ci.lb, rma03b$ci.ub), col="green3", lty=c(1,2,2))
# compare tau estimates (estimate and confidence interval):
plot(ma03b, which=4)
abline(v=confint(rma03a)$random["tau",], col="red", lty=c(1,2,2))       
abline(v=confint(rma03b)$random["tau",], col="green3", lty=c(1,3,3))       

# show heterogeneity's posterior density:
plot(ma03a, which=4, main="Sidik/Jonkman example")

# show some numbers (mode, median and mean):
abline(v=ma03a$summary[c("mode","median","mean"),"tau"], col="blue")

# compare with Sidik and Jonkman's estimates:
sj.estimates <- sqrt(c("MM"  = 0.429,   # method of moments estimator
                       "VC"  = 0.841,   # variance component type estimator
                       "ML"  = 0.562,   # maximum likelihood estimator
                       "REML"= 0.598,   # restricted maximum likelihood estimator
                       "EB"  = 0.703,   # empirical Bayes estimator
                       "MV"  = 0.818,   # model error variance estimator
                       "MVvc"= 0.747))  # a variation of the MV estimator
abline(v=sj.estimates, col="red", lty="dashed")
## End(Not run)


###########################
# example data by Cochran:
data("Cochran1954")

## Not run: 
# analysis using improper uniform prior
# (may take a few seconds to compute!):
ma04 <- bayesmeta(y=Cochran1954[,"mean"], sigma=sqrt(Cochran1954[,"se2"]),
                  label=Cochran1954[,"observer"])

# show joint posterior density:
plot(ma04, which=2, main="Cochran example")
# show (known) true parameter value:
abline(h=161)

# pick a point estimate for tau:
tau <- ma04$summary["median","tau"]
# highlight two point hypotheses (fixed vs. random effects):
abline(v=c(0, tau), col="orange", lty="dotted", lwd=2)

# show marginal posterior density:
plot(ma04, which=3)
abline(v=161)
# show the conditional distributions of the effect mu
# at two tau values corresponding to fixed and random effects models:
cm <- ma04$cond.moment(tau=c(0,tau))
mu <- seq(130,200, le=200)
lines(mu, dnorm(mu, mean=cm[1,"mean"], sd=cm[1,"sd"]), col="orange", lwd=2)
lines(mu, dnorm(mu, mean=cm[2,"mean"], sd=cm[2,"sd"]), col="orange", lwd=2)

# determine a range of tau values:
tau <- seq(0, ma04$qposterior(tau=0.99), length=100)
# compute conditional posterior moments:
cm.overall <- ma04$cond.moment(tau=tau)
# compute study-specific conditional posterior moments:
cm.indiv <- ma04$cond.moment(tau=tau, individual=TRUE)
# show forest plot along with conditional posterior means:
par(mfrow=c(1,2))
  plot(ma04, which=1, main="Cochran 1954 example")
  matplot(tau, cm.indiv[,"mean",], type="l", lty="solid", col=1:ma04$k,
          xlim=c(0,max(tau)*1.2), xlab=expression("heterogeneity "*tau),
          ylab=expression("(conditional) shrinkage estimate E["*
                           theta[i]*"|"*list(tau, y, sigma)*"]"))
  text(rep(max(tau)*1.01, ma04$k), cm.indiv[length(tau),"mean",],
       ma04$label, col=1:ma04$k, adj=c(0,0.5))
  lines(tau, cm.overall[,"mean"], lty="dashed", lwd=2)
  text(max(tau)*1.01, cm.overall[length(tau),"mean"],
       "overall", adj=c(0,0.5))
par(mfrow=c(1,1))

# show the individual effects' posterior distributions:
theta <- seq(120, 240, le=300)
plot(range(theta), c(0,0.1), type="n", xlab=expression(theta[i]), ylab="")
for (i in 1:ma04$k) {
  # draw estimate +/- uncertainty as a Gaussian:
  lines(theta, dnorm(theta, mean=ma04$y[i], sd=ma04$sigma[i]), col=i+1, lty="dotted")
  # draw effect's posterior distribution:
  lines(theta, ma04$dposterior(theta=theta, indiv=i), col=i+1, lty="solid")
}
abline(h=0)
legend(max(theta), 0.1, legend=ma04$label, col=(1:ma04$k)+1, pch=15, xjust=1, yjust=1)

## End(Not run)

---

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 $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),$

$\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 [email protected]

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)

---

Direct and indirect comparison example data

Description

Numbers of subjects and events in the different treatment arms of 22 studies.

Usage

data("BucherEtAl1997")

Format

The data frame contains the following columns:

study	character	publication identifier (first author and publication year)
treat.A	factor	treatment in first study arm (“TMP-SMX” or “AP”)
treat.B	factor	treatment in second study arm (“D/P” or “AP”)
events.A	numeric	number of events in first study arm
events.B	numeric	number of events in second study arm
total.A	numeric	total number of patients in first study arm
total.B	numeric	total number of patients in second study arm

Details

Bucher et al. (1997) discussed the example case of the comparison of sulphametoxazole-trimethoprim (TMP-SMX) versus dapsone/pyrimethamine (D/P) for the prophylaxis of Pneumocystis carinii pneumonia in HIV patients. Eight studies had undertaken a head-to-head comparison of both medications, but an additional 14 studies were available investigating one of the two medications with aerosolized pentamidine (AP) as a comparator. Nine studies compared TMP-SMX vs. AP, and five studies compared D/P vs. AP. Together these provide indirect evidence on the effect of TMP-SMX compared to D/P (Kiefer et al., 2015).

The example constitutes a simple case of a network meta-analysis (NMA) setup, where only two-armed studies are considered, and analysis is based on pairwise comparisons of treatments (or contrasts). In this case, the joint analysis of direct and indirect evidence may be implemented as a special case of a meta-regression (Higgins et al., 2019; Sec. 11.4.2). The original data in fact included some three-armed studies, in which case one of the arms was deliberately omitted (Bucher et al.; 1997).

Source

H.C. Bucher, G.H. Guyatt, L.E. Griffith, S.D. Walter. The results of direct and indirect treatment comparisons in meta-analysis of randomized controlled trials. Journal of Clinical Epidemiology, 50(6):683-691, 1997. doi:10.1016/S0895-4356(97)00049-8.

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.

J.P.T. Higgins, J. Thomas, J. Chandler, M. Cumpston, T. Li, M.J. Page, V.A. Welch (eds.). Cochrane handbook for systematic reviews of interventions. Wiley and Sons, 2nd edition, 2019. doi:10.1002/9781119536604. http://training.cochrane.org/handbook.

C. Kiefer, S. Sturtz, R. Bender. Indirect comparisons and network meta-analyses. Deutsches Aerzteblatt International, 112(47):803-808, 2015. doi:10.3238/arztebl.2015.0803.

Examples

# load data:
data("BucherEtAl1997")

# show data:
head(BucherEtAl1997)

## Not run: 
# compute effect sizes (log-ORs for pairwise comparisons)
# from the count data:
es <- escalc(measure="OR",
             ai=events.A, n1i=total.A,   # "exposure group"
             ci=events.B, n2i=total.B,   # "control group"
             slab=study, data=BucherEtAl1997)

# specify regressor matrix:
X <- cbind("TMP.DP" = rep(c(1, 0, 1), c(8,5,9)),
           "AP.DP"  = rep(c(0, 1,-1), c(8,5,9)))

# perform Bayesian meta-regression:
bmr01 <- bmr(es, X=X)

# show default output:
print(bmr01)

# specify contrast matrix:
contrastX <- rbind("TMP-SMX vs. D/P"=c(1,0),
                   "AP vs. D/P"     =c(0,1),
                   "TMP-SMX vs. AP" =c(1,-1))
# show summary including contrast estimates:
summary(bmr01, X.mean=contrastX)
# show forest plot including contrast estimates:
forestplot(bmr01, X.mean=contrastX, xlab="log-OR")


# perform frequentist meta-regression:
fmr01 <- rma(es, mods=X, intercept=FALSE)
print(fmr01)

# compare Bayesian and frequentist results;
# estimated log-OR for "TMP-SMX" vs. "D/P"
rbind("bayesmeta"=bmr01$summary[c("mean","sd"),"TMP.DP"],
      "rma"      =c(fmr01$beta["TMP.DP",], fmr01$se[1]))

# estimated log-OR for "AP" vs. "D/P"
rbind("bayesmeta"=bmr01$summary[c("mean","sd"),"AP.DP"],
      "rma"      =c(fmr01$beta["AP.DP",], fmr01$se[2]))

# estimated heterogeneity:
rbind("bayesmeta"=bmr01$summary["median","tau"],
      "rma"      =sqrt(fmr01$tau2))

## End(Not run)

---

Fly counts example data

Description

This data set gives average estimated counts of flies along with standard errors from 7 different observers.

Usage

data("Cochran1954")

Format

The data frame contains the following columns:

observer	character	identifier
mean	numeric	mean count
se2	numeric	squared standard error

Details

Quoting from Cochran (1954), example 3, p.119: “In studies by the U.S. Public Health Service of observers' abilities to count the number of flies which settle momentarily on a grill, each of 7 observers was shown, for a brief period, grills with known numbers of flies impaled on them and asked to estimate the numbers. For a given grill, each observer made 5 independent estimates. The data in table 9 are for a grill which actually contained 161 flies. Estimated variances are based on 4 degrees of freedom each. [...] The only point of interest in estimating the overall mean is to test whether there is any consistent bias among observers in estimating the 161 flies on the grill. Although inspection of table 9 suggests no such bias, the data will serve to illustrate the application of partial weighting.”

Source

W.G. Cochran. The combination of estimates from different experiments. Biometrics, 10(1):101-129, 1954.

Examples

data("Cochran1954")
## Not run: 
# analysis using improper uniform prior
# (may take a few seconds to compute!):
bma <- bayesmeta(y=Cochran1954[,"mean"], sigma=sqrt(Cochran1954[,"se2"]),
                 label=Cochran1954[,"observer"])

# show joint posterior density:
plot(bma, which=2, main="Cochran example")
# show (known) true parameter value:
abline(h=161)

# show forest plot:
forestplot(bma, zero=161)

## End(Not run)

---

Convolution of two probability distributions

Description

Compute the convolution of two probability distributions, specified through their densities or cumulative distribution functions (CDFs).

Usage

convolve(dens1, dens2,
           cdf1=Vectorize(function(x){integrate(dens1,-Inf,x)$value}),
           cdf2=Vectorize(function(x){integrate(dens2,-Inf,x)$value}),
           delta=0.01, epsilon=0.0001)

Arguments

dens1, dens2	the two distributions' probability density functions.
cdf1, cdf2	the two distributions' cumulative distribution functions.
delta, epsilon	the parameters specifying the desired accuracy for approximation of the convolution, and with that determining the number of support points being used internally. Smaller values imply greater accuracy and greater computational burden (Roever and Friede, 2017).

Details

The distribution of the sum of two (independent) random variables technically results as a convolution of their probability distributions. In some cases, the calculation of convolutions may be done analytically; e.g., the sum of two normally distributed random variables again turns out as normally distributed (with mean and variance resulting as the sums of the original ones). In other cases, convolutions may need to be determined numerically. One way to achieve this is via the DIRECT algorithm; the present implementation is the one discussed by Roever and Friede (2017). Accuracy of the computations is determined by the delta (maximum divergence $\delta$) and epsilon (tail probability $\epsilon$) parameters.

Convolutions here are used within the funnel() function (to generate funnel plots), but are often useful more generally. The original probability distributions may be specified via their probability density functions or their cumulative distribution functions (CDFs). The convolve() function returns the convolution's density, CDF and quantile function (inverse CDF).

Value

A list with elements

delta	the $\delta$ parameter.
epsilon	the $\epsilon$ parameter.
binwidth	the bin width.
bins	the total number of bins.
support	a matrix containing the support points used internally for the convolution approximation.
density	the probability density function.
cdf	the cumulative distribution function (CDF).
quantile	the quantile function (inverse CDF).

Author(s)

Christian Roever [email protected]

References

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.

See Also

bayesmeta, funnel.

Examples

## Not run: 
#  Skew-normal / logistic example:

dens1 <- function(x, shape=4)
# skew-normal distribution's density
# see also: http://azzalini.stat.unipd.it/SN/Intro
{
  return(2 * dnorm(x) * pnorm(shape * x))
}

dens2 <- function(x)
# logistic distribution's density
{
  return(dlogis(x, location=0, scale=1))
}

rskewnorm <- function(n, shape=4)
# skew-normal random number generation
# (according to http://azzalini.stat.unipd.it/SN/faq-r.html)
{
  delta <- shape / sqrt(shape^2+1)
  u1 <- rnorm(n); v <- rnorm(n)
  u2 <- delta * u1 + sqrt(1-delta^2) * v
  return(apply(cbind(u1,u2), 1, function(x){ifelse(x[1]>=0, x[2], -x[2])}))
}

# compute convolution:
conv <- convolve(dens1, dens2)

# illustrate convolution:
n <- 100000
x <- rskewnorm(n)
y <- rlogis(n)
z <- x + y

# determine empirical and theoretical quantiles:
p      <- c(0.001,0.01, 0.1, 0.5, 0.9, 0.99, 0.999)
equant <- quantile(z, prob=p)
tquant <- conv$quantile(p)

# show numbers:
print(cbind("p"=p, "empirical"=equant, "theoretical"=tquant))

# draw Q-Q plot:
rg <- range(c(equant, tquant))
plot(rg, rg, type="n", asp=1, main="Q-Q-plot",
     xlab="theoretical quantile", ylab="empirical quantile")
abline(0, 1, col="grey")
points(tquant, equant)

## End(Not run)

---

Pediatric liver transplant example data

Description

Numbers of cases (transplant patients) and events (acute rejections, steroid resistant rejections, PTLDs, and deaths) in experimental and control groups of six studies.

Usage

data("CrinsEtAl2014")

Format

The data frame contains the following columns:

publication	character	publication identifier (first author and publication year)
year	numeric	publication year
randomized	factor	randomization status (y/n)
control.type	factor	type of control group (‘concurrent’ or ‘historical’)
comparison	factor	type of comparison (‘IL-2RA only’, ‘delayed CNI’, or ‘no/low steroids’)
IL2RA	factor	type of interleukin-2 receptor antagonist (IL-2RA) (‘basiliximab’ or ‘daclizumab’)
CNI	factor	type of calcineurin inhibitor (CNI) (‘tacrolimus’ or ‘cyclosporine A’)
MMF	factor	use of mycofenolate mofetil (MMF) (y/n)
followup	numeric	follow-up time in months
treat.AR.events	numeric	number of AR events in experimental group
treat.SRR.events	numeric	number of SRR events in experimental group
treat.PTLD.events	numeric	number of PTLD events in experimental group
treat.deaths	numeric	number of deaths in experimental group
treat.total	numeric	number of cases in experimental group
control.AR.events	numeric	number of AR events in control group
control.SRR.events	numeric	number of SRR events in control group
control.PTLD.events	numeric	number of PTLD events in control group
control.deaths	numeric	number of deaths in control group
control.total	numeric	number of cases in control group

Details

A systematic literature review investigated the evidence on the effect of Interleukin-2 receptor antagonists (IL-2RA) and resulted in six controlled studies reporting acute rejection (AR), steroid-resistant rejection (SRR) and post-transplant lymphoproliferative disorder (PTLD) rates as well as mortality in pediatric liver transplant recipients.

Source

N.D. Crins, C. Roever, A.D. Goralczyk, T. Friede. Interleukin-2 receptor antagonists for pediatric liver transplant recipients: A systematic review and meta-analysis of controlled studies. Pediatric Transplantation, 18(8):839-850, 2014. doi:10.1111/petr.12362.

References

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

T.G. Heffron et al. Pediatric liver transplantation with daclizumab induction therapy. Transplantation, 75(12):2040-2043, 2003. doi:10.1097/01.TP.0000065740.69296.DA.

N.E.M. Gibelli et al. Basiliximab-chimeric anti-IL2-R monoclonal antibody in pediatric liver transplantation: comparative study. Transplantation Proceedings, 36(4):956-957, 2004. doi:10.1016/j.transproceed.2004.04.070.

S. Schuller et al. Daclizumab induction therapy associated with tacrolimus-MMF has better outcome compared with tacrolimus-MMF alone in pediatric living donor liver transplantation. Transplantation Proceedings, 37(2):1151-1152, 2005. doi:10.1016/j.transproceed.2005.01.023.

R. Ganschow et al. Long-term results of basiliximab induction immunosuppression in pediatric liver transplant recipients. Pediatric Transplantation, 9(6):741-745, 2005. doi:10.1111/j.1399-3046.2005.00371.x.

M. Spada et al. Randomized trial of basiliximab induction versus steroid therapy in pediatric liver allograft recipients under tacrolimus immunosuppression. American Journal of Transplantation, 6(8):1913-1921, 2006. doi:10.1111/j.1600-6143.2006.01406.x.

J.M. Gras et al. Steroid-free, tacrolimus-basiliximab immunosuppression in pediatric liver transplantation: Clinical and pharmacoeconomic study in 50 children. Liver Transplantation, 14(4):469-477, 2008. doi:10.1002/lt.21397.

See Also

GoralczykEtAl2011.

Examples

data("CrinsEtAl2014")
## Not run: 
# compute effect sizes (log odds ratios) from count data
# (using "metafor" package's "escalc()" function):
require("metafor")
crins.es <- escalc(measure="OR",
                   ai=exp.AR.events,  n1i=exp.total,
                   ci=cont.AR.events, n2i=cont.total,
                   slab=publication, data=CrinsEtAl2014)
print(crins.es)

# analyze using weakly informative half-Cauchy prior for heterogeneity:
crins.ma <- bayesmeta(crins.es, tau.prior=function(t){dhalfcauchy(t,scale=1)})

# show results:
print(crins.ma)
forestplot(crins.ma)
plot(crins.ma)

# show heterogeneity posterior along with prior:
plot(crins.ma, which=4, prior=TRUE)

# perform meta analysis using 2 randomized studies only
# but use 4 non-randomized studies to inform heterogeneity prior:
crins.nrand <- bayesmeta(crins.es[crins.es$randomized=="no",],
                         tau.prior=function(t){dhalfcauchy(t,scale=1)})
crins.rand  <- bayesmeta(crins.es[crins.es$randomized=="yes",],
                         tau.prior=function(t){crins.nrand$dposterior(tau=t)})
plot(crins.nrand, which=4, prior=TRUE,
     main="non-randomized posterior = randomized prior")
plot(crins.rand, which=4, prior=TRUE, main="randomized posterior")
plot(crins.rand, which=1)

## End(Not run)

---

Half-logistic distribution.

Description

Half-logistic density, distribution, and quantile functions, random number generation and expectation and variance.

Usage

dhalflogistic(x, scale=1, log=FALSE)
  phalflogistic(q, scale=1)
  qhalflogistic(p, scale=1)
  rhalflogistic(n, scale=1)
  ehalflogistic(scale=1)
  vhalflogistic(scale=1)

Arguments

x, q	quantile.
p	probability.
n	number of observations.
scale	scale parameter ($>0$).
log	logical; if TRUE, logarithmic density will be returned.

Details

The half-logistic distribution is simply a zero-mean logistic distribution that is restricted to take only positive values. If $X\sim\mathrm{logistic}$, then $|sX|\sim\mathrm{halflogistic}(\mathrm{scale}\!=\!s)$.

Value

‘dhalflogistic()’ gives the density function, ‘phalflogistic()’ gives the cumulative distribution function (CDF), ‘qhalflogistic()’ gives the quantile function (inverse CDF), and ‘rhalflogistic()’ generates random deviates. The ‘ehalflogistic()’ and ‘vhalflogistic()’ functions return the corresponding half-logistic distribution's expectation and variance, respectively.

Author(s)

Christian Roever [email protected]

References

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.

N.L. Johnson, S. Kotz, N. Balakrishnan. Continuous univariate distributions, volume 2, chapter 23.11. Wiley, New York, 2nd edition, 1994.

See Also

dlogis, dhalfnormal, dlomax, drayleigh, TurnerEtAlPrior, RhodesEtAlPrior, bayesmeta.

Examples

#######################
# illustrate densities:
x <- seq(0,6,le=200)
plot(x, dhalfnormal(x), type="l", col="red", ylim=c(0,1),
     xlab=expression(tau), ylab=expression("probability density "*f(tau)))
lines(x, dhalflogistic(x), col="green3")
lines(x, dhalfcauchy(x), col="blue")
lines(x, dexp(x), col="cyan")
abline(h=0, v=0, col="grey")

# show log-densities (note the differing tail behaviour):
plot(x, dhalfnormal(x), type="l", col="red", ylim=c(0.001,1), log="y",
     xlab=expression(tau), ylab=expression("probability density "*f(tau)))
lines(x, dhalflogistic(x), col="green3")
lines(x, dhalfcauchy(x), col="blue")
lines(x, dexp(x), col="cyan")
abline(v=0, col="grey")

---

Half-normal, half-Student-t and half-Cauchy distributions.

Description

Half-normal, half-Student-t and half-Cauchy density, distribution, quantile functions, random number generation, and expectation and variance.

Usage

dhalfnormal(x, scale=1, log=FALSE)
  phalfnormal(q, scale=1)
  qhalfnormal(p, scale=1)
  rhalfnormal(n, scale=1)
  ehalfnormal(scale=1)
  vhalfnormal(scale=1)

  dhalft(x, df, scale=1, log=FALSE)
  phalft(q, df, scale=1)
  qhalft(p, df, scale=1)
  rhalft(n, df, scale=1)
  ehalft(df, scale=1)
  vhalft(df, scale=1)

  dhalfcauchy(x, scale=1, log=FALSE)
  phalfcauchy(q, scale=1)
  qhalfcauchy(p, scale=1)
  rhalfcauchy(n, scale=1)
  ehalfcauchy(scale=1)
  vhalfcauchy(scale=1)

Arguments

x, q	quantile.
p	probability.
n	number of observations.
scale	scale parameter ($>0$).
df	degrees-of-freedom parameter ($>0$).
log	logical; if TRUE, logarithmic density will be returned.

Details

The half-normal distribution is simply a zero-mean normal distribution that is restricted to take only positive values. The scale parameter $\sigma$ here corresponds to the underlying normal distribution's standard deviation: if $X\sim\mathrm{Normal}(0,\sigma^2)$, then $|X|\sim\mathrm{halfNormal}(\mathrm{scale}\!=\!\sigma)$. Its mean is $\sigma \sqrt{2/\pi}$, and its variance is $\sigma^2 (1-2/\pi)$. Analogously, the half-t distribution is a truncated Student-t distribution with df degrees-of-freedom, and the half-Cauchy distribution is again a special case of the half-t distribution with df=1 degrees of freedom.

Note that (half-) Student-t and Cauchy distributions arise as continuous mixture distributions of (half-) normal distributions. If

$Y|\sigma\;\sim\;\mathrm{Normal}(0,\sigma^2)$

where the standard deviation is $\sigma = \sqrt{k/X}$ and $X$ is drawn from a $\chi^2$-distribution with $k$ degrees of freedom, then the marginal distribution of $Y$ is Student-t with $k$ degrees of freedom.

Value

‘dhalfnormal()’ gives the density function, ‘phalfnormal()’ gives the cumulative distribution function (CDF), ‘qhalfnormal()’ gives the quantile function (inverse CDF), and ‘rhalfnormal()’ generates random deviates. The ‘ehalfnormal()’ and ‘vhalfnormal()’ functions return the corresponding half-normal distribution's expectation and variance, respectively. For the ‘dhalft()’, ‘dhalfcauchy()’ and related function it works analogously.

Author(s)

Christian Roever [email protected]

References

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.

A. Gelman. Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3):515-534, 2006. doi:10.1214/06-BA117A.

F. C. Leone, L. S. Nelson, R. B. Nottingham. The folded normal distribution. Technometrics, 3(4):543-550, 1961. doi:10.2307/1266560.

N. G. Polson, J. G. Scott. On the half-Cauchy prior for a global scale parameter. Bayesian Analysis, 7(4):887-902, 2012. doi:10.1214/12-BA730.

S. Psarakis, J. Panaretos. The folded t distribution. Communications in Statistics - Theory and Methods, 19(7):2717-2734, 1990. doi:10.1080/03610929008830342.

See Also

dnorm, dt, dcauchy, dlomax, drayleigh, TurnerEtAlPrior, RhodesEtAlPrior, bayesmeta.

Examples

#######################
# illustrate densities:
x <- seq(0,6,le=200)
plot(x, dhalfnormal(x), type="l", col="red", ylim=c(0,1),
     xlab=expression(tau), ylab=expression("probability density "*f(tau)))
lines(x, dhalft(x, df=3), col="green")
lines(x, dhalfcauchy(x), col="blue")
lines(x, dexp(x), col="cyan")
abline(h=0, v=0, col="grey")

# show log-densities (note the differing tail behaviour):
plot(x, dhalfnormal(x), type="l", col="red", ylim=c(0.001,1), log="y",
     xlab=expression(tau), ylab=expression("probability density "*f(tau)))
lines(x, dhalft(x, df=3), col="green")
lines(x, dhalfcauchy(x), col="blue")
lines(x, dexp(x), col="cyan")
abline(v=0, col="grey")

---

Inverse-Chi distribution.

Description

(Scaled) inverse-Chi density, distribution, and quantile functions, random number generation and expectation and variance.

Usage

dinvchi(x, df, scale=1, log=FALSE)
  pinvchi(q, df, scale=1, lower.tail=TRUE, log.p=FALSE)
  qinvchi(p, df, scale=1, lower.tail=TRUE, log.p=FALSE)
  rinvchi(n, df, scale=1)
  einvchi(df, scale=1)
  vinvchi(df, scale=1)

Arguments

x, q	quantile.
p	probability.
n	number of observations.
df	degrees-of-freedom parameter ($>0$).
scale	scale parameter ($>0$).
log	logical; if TRUE, logarithmic density will be returned.
lower.tail	logical; if TRUE (default), probabilities are P(X <= x), otherwise, P(X > x).
log.p	logical; if TRUE, probabilities p are returned as log(p).

Details

The (scaled) inverse-Chi distribution is defined as the distribution of the (scaled) inverse of the square root of a Chi-square-distributed random variable. It is a special case of the square-root inverted-gamma distribution (with $\alpha=\nu/2$ and $\beta=1/2$) (Bernardo and Smith; 1994). Its probability density function is given by

$p(x) \;=\; \frac{2^{(1-\nu/2)}}{s \, \Gamma(\nu/2)} \Bigl(\frac{s}{x}\Bigr)^{(\nu+1)} \exp\Bigl(-\frac{s^2}{2\,x^2}\Bigr)$

where $\nu$ is the degrees-of-freedom and $s$ the scale parameter.

Value

‘dinvchi()’ gives the density function, ‘pinvchi()’ gives the cumulative distribution function (CDF), ‘qinvchi()’ gives the quantile function (inverse CDF), and ‘rinvchi()’ generates random deviates. The ‘einvchi()’ and ‘vinvchi()’ functions return the corresponding distribution's expectation and variance, respectively.

Author(s)

Christian Roever [email protected]

References

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.

J.M. Bernardo, A.F.M. Smith. Bayesian theory, Appendix A.1. Wiley, Chichester, UK, 1994.

See Also

dhalfnormal, dhalft.

Examples

#################################
# illustrate Chi^2 - connection;
# generate Chi^2-draws:
chi2 <- rchisq(1000, df=10)
# transform:
invchi <- sqrt(1 / chi2)
# show histogram:
hist(invchi, probability=TRUE, col="grey")
# show density for comparison:
x <- seq(0, 1, length=100)
lines(x, dinvchi(x, df=10, scale=1), col="red")
# compare theoretical and empirical moments:
rbind("theoretical" = c("mean" = einvchi(df=10, scale=1),
                        "var"  = vinvchi(df=10, scale=1)),
      "Monte Carlo" = c("mean" = mean(invchi),
                        "var"  = var(invchi)))

##############################################################
# illustrate the normal/Student-t - scale mixture connection;
# specify degrees-of-freedom:
df <- 5
# generate standard normal draws:
z <- rnorm(1000)
# generate random scalings:
sigma <- rinvchi(1000, df=df, scale=sqrt(df))
# multiply to yield Student-t draws:
t <- z * sigma
# check Student-t distribution via a Q-Q-plot:
qqplot(qt(ppoints(length(t)), df=df), t)
abline(0, 1, col="red")

---