R/bayesmeta.R

Defines functions traceplot.bayesmeta traceplot weightsplot.bayesmeta weightsplot ess.bayesmeta ess uisd.escalc uisd.default uisd pppvalue normalmixture funnel.bayesmeta convolve forestplot.escalc forestplot.bayesmeta forest.bayesmeta RhodesEtAlPrior TurnerEtAlPrior vinvchi einvchi rinvchi qinvchi pinvchi dinvchi vrayleigh erayleigh rrayleigh qrayleigh prayleigh drayleigh vhalflogistic ehalflogistic rhalflogistic qhalflogistic phalflogistic dhalflogistic vhalfcauchy ehalfcauchy rhalfcauchy qhalfcauchy phalfcauchy dhalfcauchy vhalft ehalft rhalft qhalft phalft dhalft vhalfnormal ehalfnormal rhalfnormal qhalfnormal phalfnormal dhalfnormal vlomax elomax rlomax qlomax plomax dlomax plot.bayesmeta summary.bayesmeta print.bayesmeta bayesmeta.escalc bayesmeta.default bayesmeta

Documented in bayesmeta bayesmeta.default bayesmeta.escalc convolve dhalfcauchy dhalflogistic dhalfnormal dhalft dinvchi dlomax drayleigh ehalfcauchy ehalflogistic ehalfnormal ehalft einvchi elomax erayleigh ess ess.bayesmeta forest.bayesmeta forestplot.bayesmeta forestplot.escalc funnel.bayesmeta normalmixture phalfcauchy phalflogistic phalfnormal phalft pinvchi plomax plot.bayesmeta pppvalue prayleigh print.bayesmeta qhalfcauchy qhalflogistic qhalfnormal qhalft qinvchi qlomax qrayleigh rhalfcauchy rhalflogistic rhalfnormal rhalft RhodesEtAlPrior rinvchi rlomax rrayleigh summary.bayesmeta traceplot traceplot.bayesmeta TurnerEtAlPrior uisd uisd.default uisd.escalc vhalfcauchy vhalflogistic vhalfnormal vhalft vinvchi vlomax vrayleigh weightsplot weightsplot.bayesmeta

#
#    bayesmeta, an R package for Bayesian random-effects meta-analysis.
#    Copyright (C) 2023  Christian Roever
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#


bayesmeta <- function(y,...)
{
  UseMethod("bayesmeta")
}


bayesmeta.default <- function(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,...)
{
  ptm <- proc.time()
  y      <- as.vector(y)
  sigma  <- as.vector(sigma)
  labels <- as.vector(labels)
  # some preliminary sanity checks:
  stopifnot(is.vector(y), is.vector(sigma),
            all(is.finite(y)), all(is.finite(sigma)),
            length(sigma)==length(y),
            all(sigma>=0), sum(sigma==0)<=1,
            length(mu.prior)==2,
            length(mu.prior.mean)==1, length(mu.prior.sd)==1,
            is.na(mu.prior.mean) || is.finite(mu.prior.mean),
            is.na(mu.prior.sd) || (is.finite(mu.prior.mean) && (mu.prior.sd>0)),
            ((is.na(mu.prior.mean) & is.na(mu.prior.sd))
             || (is.finite(mu.prior.mean) & is.finite(mu.prior.sd))),
            (is.function(tau.prior) | (is.character(tau.prior) && (length(tau.prior)==1))))
  interval.type <- match.arg(interval.type)
  stopifnot(length(interval.type)==1, is.element(interval.type, c("shortest","central")))

  zerosigma <- (sigma == 0.0)
  #if (any(zerosigma)) warning("one of the supplied 'sigma' elements is zero!")
  
  k <- length(y)
  if (is.null(labels))
    labels <- sprintf("%02d", 1:k)
  sigma2hat <- (k-1)*sum(1/sigma^2) / (sum(1/sigma^2)^2 - sum(1/sigma^4)) # Higgins/Thompson (2002), eqn. (9)

  maxratio <- max(sigma) / min(sigma)
  if ((!any(zerosigma)) && (maxratio > 1000))
    warning(paste0("Ratio of largest over smallest standard error (sigma) is ", sprintf("%.0f",maxratio), ". Extreme values may lead to computational problems."))

  tau.prior.proper <- NA
  if (is.character(tau.prior)) {
    tau.prior <- match.arg(tolower(tau.prior),
                           c("uniform","jeffreys","shrinkage","dumouchel","bergerdeely","conventional","i2","sqrt"))
    tau.prior <- c("uniform"="uniform", "jeffreys"="Jeffreys",
                   "shrinkage"="shrinkage", "dumouchel"="DuMouchel",
                   "bergerdeely"="BergerDeely", "conventional"="conventional", "i2"="I2", "sqrt"="sqrt")[tau.prior]
    stopifnot(is.element(tau.prior, c("uniform", "Jeffreys", "shrinkage", "DuMouchel", "BergerDeely", "conventional", "I2", "sqrt")))
    if (tau.prior=="uniform") {             # uniform prior on tau:
      pdens <- function(t){d<-rep(1,length(t)); d[t<0]<-0; return(d)}
      attr(pdens, "bayesmeta.label") <- "uniform(min=0, max=Inf)"
    } else if (tau.prior=="Jeffreys") {     # Jeffreys prior:
      pdens <- function(t){return(apply(matrix(t,ncol=1),1,function(x){return(sqrt(sum((x/(sigma^2+x^2))^2)))}))}
      attr(pdens, "bayesmeta.label") <- "Jeffreys prior"
    } else if (tau.prior=="BergerDeely") {  # Berger/Deely prior:
      pdens <- function(t){return(apply(matrix(t,ncol=1),1,
                                        function(x){return(exp(log(x)-sum(log(sigma^2+x^2))/k))}))}
      attr(pdens, "bayesmeta.label") <- "Berger/Deely prior"
    } else if (tau.prior=="conventional") { # conventional prior:
      pdens <- function(t){return(apply(matrix(t,ncol=1),1,
                                        function(x){return(exp(log(x)-(3/(2*k))*sum(log(sigma^2+x^2))))}))}
      attr(pdens, "bayesmeta.label") <- "conventional prior"
    } else if (tau.prior=="I2") {           # uniform on I^2:
      pdens <- function(t){return(apply(matrix(t,ncol=1),1,
                                        function(x){return(2 * exp(log(sigma2hat)+log(x) - 2*log(sigma2hat+x^2)))}))}
      attr(pdens, "bayesmeta.label") <- "uniform prior on I-squared"
    } else if (tau.prior=="sqrt") {         # sqrt-prior:
      pdens <- function(t){d <- rep(0,length(t)); d[t>=0] <- t[t>=0]^(-0.5); return(d)}
      attr(pdens, "bayesmeta.label") <- "uniform prior on sqrt(tau)"
    } else {
      # harmonic mean of squared standard errors:
      s02 <- k/sum(1/sigma^2)
      if (tau.prior=="shrinkage") {         # "uniform shrinkage" prior:
        pdens <- function(t){return(apply(matrix(t,ncol=1),1,function(x){return(2*x*s02/(s02+x^2)^2)}))}
        attr(pdens, "bayesmeta.label") <- "uniform shrinkage prior"
      } else if (tau.prior=="DuMouchel") {  # DuMouchel prior:
        s0 <- sqrt(s02)
        pdens <- function(t){return(apply(matrix(t,ncol=1),1,function(x){return(s0/(s0+x)^2)}))}
        attr(pdens, "bayesmeta.label") <- "DuMouchel prior"
      } else warning("could not make sense of 'tau.prior' argument")
    }
    if (is.element(tau.prior, c("uniform", "Jeffreys", "BergerDeely", "sqrt")))
      tau.prior.proper <- FALSE
    tau.prior <- pdens
    rm("pdens")
  }

  tau.prior.integral <- 1.0  # heterogeneity prior's normalizing constant

  dprior <- function(tau=NA, mu=NA, log=FALSE)
  # prior density (marginal or joint)
  {
    if (all(is.na(tau))) { # marginal density for mu:
      if (is.finite(mu.prior.mean))
        result <- dnorm(mu, mean=mu.prior.mean, sd=mu.prior.sd, log=log)
      else
        result <- rep(ifelse(log, 0, 1), length(mu))
    }
    else if (all(is.na(mu))) { # marginal density for tau:
      if (log) {
        if (is.element("log", names(formals(tau.prior))))
          result <- apply(matrix(tau,ncol=1), 1, tau.prior, log=TRUE) - log(tau.prior.integral)
        else
          result <- log(apply(matrix(tau,ncol=1), 1, tau.prior)) - log(tau.prior.integral)
      }
      else result <- apply(matrix(tau,ncol=1), 1, tau.prior) / tau.prior.integral
    }
    else if (is.finite(mu.prior.mean) & !all(is.na(mu))) {  # joint density:
      if (is.element("log", names(formals(tau.prior))))
        result <- (apply(matrix(tau,ncol=1), 1, tau.prior, log=TRUE) - log(tau.prior.integral)
                   + dnorm(mu, mean=mu.prior.mean, sd=mu.prior.sd, log=TRUE))
      else
        result <- (log(apply(matrix(tau,ncol=1), 1, tau.prior)) - log(tau.prior.integral)
                   + dnorm(mu, mean=mu.prior.mean, sd=mu.prior.sd, log=TRUE))
      if (!log) result <- exp(result)
    }
    else {  # joint density (uniform on mu):
      if (log) {
        if (is.element("log", names(formals(tau.prior))))
          result <- apply(matrix(tau,ncol=1), 1, tau.prior, log=TRUE) - log(tau.prior.integral)
        else
          result <- log(apply(matrix(tau,ncol=1), 1, tau.prior)) - log(tau.prior.integral)
      }
      else result <- apply(matrix(tau,ncol=1), 1, tau.prior) / tau.prior.integral
    }
    return(result)
  }

  # check whether tau prior is proper
  # (unless obvious from above specification)
  # by trying to integrate numerically:
  if (is.na(tau.prior.proper)) {
    prior.int <- integrate(function(t){return(dprior(tau=t))},
                           lower=0, upper=Inf,
                           rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate,
                           stop.on.error=FALSE)
    tau.prior.proper <- ((prior.int$message == "OK") && (prior.int$value > 0))
    # if integral is substantially different from 1.0, apply correction:
    if (tau.prior.proper && (abs(prior.int$value - 1.0) > prior.int$abs.error))
      tau.prior.integral <- prior.int$value 
    rm("prior.int")
  }
  
  likelihood <- function(tau=NA, mu=NA, log=FALSE)
  # likelihood function (marginal or joint)
  {
    if (all(is.na(mu)) & all(is.na(tau))) {
      warning("need to supply at least either 'mu' or 'tau'")
      return()
    }
    else if (all(is.na(tau))) { # return marginal likelihood (numerical):
      loglikeli <- function(taumu)
      {
        t2s2 <- taumu[1]^2 + sigma^2
        return(-(k/2)*log(2*pi) - 0.5*sum(log(t2s2)) - 0.5*sum((y-taumu[2])^2 / t2s2))
      }
      marglikeli <- function(mu)
      {
        integrand <- function(t){return(exp(loglikeli(c(t,mu)) + dprior(tau=t,log=TRUE)))}
        int <- integrate(Vectorize(integrand),
                         rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate,
                         lower=0, upper=Inf, stop.on.error=FALSE)
        return(ifelse(int$message == "OK", int$value, NA))
      }
      result <- apply(matrix(mu,ncol=1), 1, marglikeli)
      if (log) result <- log(result)
    } else if (all(is.na(mu))) {  # return marginal likelihood (analytical):
      logmarglikeli <- function(t)
      {
        if (is.na(mu.prior.mean)) { # uniform prior
          yTerm   <- y
          varTerm <- t^2 + sigma^2
        }
        else {                      # conjugate normal prior
          yTerm   <- c(mu.prior.mean, y)
          varTerm <- c(mu.prior.sd^2, t^2 + sigma^2)
        }
        if ((t==0) & any(zerosigma)) {
          conditionalmean <- y[zerosigma]
        } else {
          conditionalmean <- sum(yTerm/varTerm) / sum(1/varTerm)
        }
        return(-0.5*((length(yTerm)-1) * log(2*pi)
                     + sum(log(varTerm))
                     + sum((yTerm-conditionalmean)^2 / varTerm)
                     + log(sum(1/varTerm))))
      }
      result <- apply(matrix(tau, ncol=1), 1, logmarglikeli)
      result[tau<0] <- -Inf
      if (!log) result <- exp(result)
    }
    else {  # return joint likelihood:
      loglikeli <- function(taumu)
      {
        t2s2 <- taumu[1]^2 + sigma^2
        return(-(k/2)*log(2*pi) - 0.5*sum(log(t2s2)) - 0.5*sum((y-taumu[2])^2 / t2s2))
      }
      result <- apply(cbind(tau,mu), 1, loglikeli)
      result[tau<0] <- -Inf
      if (!log) result <- exp(result)
    }
    return(result)
  }
  
  dposterior <- function(tau=NA, mu=NA, theta=mu, log=FALSE, predict=FALSE, individual=FALSE)
  # posterior density (marginal or joint)
  {
    if (all(is.na(mu))) mu <- theta
    stopifnot(length(individual)==1)
    if (! (is.logical(individual) && (!individual))) {
      indiv.logi <- TRUE
      if (is.numeric(individual))   indiv.which <- which(is.element(1:k, individual))
      if (is.character(individual)) indiv.which <- which(is.element(labels, match.arg(individual,labels)))
      if (length(indiv.which)==0) warning("cannot make sense of 'individual' argument: empty subset.")
    }
    else indiv.logi <- FALSE
    if (predict & indiv.logi) warning("need to specify either 'predict' or 'individual' argument, but not both.")
    if (all(is.na(mu)) & all(is.na(tau))) {
      warning("need to supply at least either 'mu' or 'tau'")
      return()
    }
    else if (all(is.na(tau))) { # return marginal posterior (effect mu, numerical):
      if (all(is.na(support))) {
        warning("'support' not initialized.")
        result <- NA
      }
      else {
        result <- numeric(length(mu))
        if (predict)           # posterior predictive distribution
          for (i in 1:length(mu))
            result[i] <- sum(support[,"weight"]
                             * dnorm(mu[i], mean=support[,"mean"], sd=support[,"sd.pred"]))
        else if (indiv.logi) { # individual-effect distribution
          musigma <-   conditionalmoment(tau=support[,"tau"], individual=indiv.which)
          for (i in 1:length(mu))
            result[i] <- sum(support[,"weight"]
                             * dnorm(mu[i], mean=musigma[,"mean"], sd=musigma[,"sd"]))          
        }
        else                   # (marginal) posterior distribution
          for (i in 1:length(mu))
            result[i] <- sum(support[,"weight"]
                             * dnorm(mu[i], mean=support[,"mean"], sd=support[,"sd"]))
        if (log) result <- log(result)
      }
    }
    else {
      if (predict) warning("'predict' argument ignored!")
      if (indiv.logi) warning("'individual' argument ignored!")
      if (all(is.na(mu))) {  # return marginal posterior (heterogeneity tau, analytical):
        logmargpost <- function(t)
        {
          return(ifelse(t<0, -Inf, log(tau.prior(t)) + likelihood(mu=NA, tau=t, log=TRUE)))
        }
        result <- apply(matrix(tau, ncol=1), 1, logmargpost) - log(integral)
        result[is.nan(result)] <- -Inf
        if (!log) result <- exp(result)
      }
      else {  # return joint posterior:
        loglikeli <- function(taumu)
        {
          t2s2 <- taumu[1]^2 + sigma^2
          return(-(k/2)*log(2*pi) - 0.5*sum(log(t2s2)) - 0.5*sum((y-taumu[2])^2 / t2s2))
        }
        logpost <- function(taumu)
        {
          return(dprior(taumu[1], taumu[2], log=TRUE) + loglikeli(taumu) - log(integral))
        }
        result <- apply(cbind(tau,mu), 1, logpost)
        result[tau<0] <- -Inf
        if (!log) result <- exp(result)
      }
    }
    return(result)
  }

  # compute marginal posterior density's normalizing constant:
  integral <- 1.0
  integral <- integrate(function(x){dposterior(tau=x)}, lower=0, upper=Inf,
                        rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)$value
  if ((!is.finite(integral)) || (integral <= 0))
    warning("failed integrating marginal posterior (tau)")

  pposterior <- function(tau=NA, mu=NA, theta=mu, predict=FALSE, individual=FALSE)
  # posterior cumulative distribution function (CDF) of tau
  {
    if (all(is.na(mu))) mu <- theta
    stopifnot(length(individual)==1)
    if (! (is.logical(individual) && (!individual))) {
      indiv.logi <- TRUE
      if (is.numeric(individual))   indiv.which <- which(is.element(1:k, individual))
      if (is.character(individual)) indiv.which <- which(is.element(labels, match.arg(individual,labels)))
      if (length(indiv.which)==0) warning("cannot make sense of 'individual' argument: empty subset.")
    }
    else indiv.logi <- FALSE
    if (predict & indiv.logi) warning("need to specify either 'predict' or 'individual' argument, but not both.")
    if (all(is.na(mu))) { # numerical integration for tau
      if (predict) warning("'predict' argument ignored!")
      if (indiv.logi) warning("'individual' argument ignored!")
      cdf <- function(x)
      # marginal CDF of tau
      {
        if (x<=0) p <- 0
        else if (x==Inf) p <- 1
        else p <- integrate(dposterior, lower=0, upper=x,
                            rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)$value
        return(p)
      }
      result <- apply(matrix(tau, ncol=1), 1, cdf)
    }
    else if (all(is.na(tau))) { # grid approximation for mu
      if (all(is.na(support))) {
        warning("'support' not initialized.")
        result <- NA
      }
      else if (indiv.logi) { # individual-effect distribution
        musigma <-   conditionalmoment(tau=support[,"tau"], individual=indiv.which)
        cdf <- function(x)
        {
          p <- sum(pnorm(x, mean=musigma[,"mean"], sd=musigma[,"sd"]) * support[,"weight"])
          return(p)
        }
        result <- apply(matrix(mu, ncol=1), 1, cdf)
      }
      else { # posterior or posterior predictive CDF:
        cdf <- function(x)
        {
          p <- sum(pnorm(x, mean=support[,"mean"], sd=support[,ifelse(predict, "sd.pred", "sd")])
                   * support[,"weight"])
          return(p)
        }
        result <- apply(matrix(mu, ncol=1), 1, cdf)
      }
    }
    else {
      warning("need to supply EITHER tau OR mu, not both.")
      result <- NA
    }
    return(result)
  }

  qposterior <- function(tau.p=NA, mu.p=NA, theta.p=mu.p, predict=FALSE, individual=FALSE)
  # posterior quantile function of tau or mu
  {
    if (all(is.na(mu.p))) mu.p <- theta.p
    stopifnot(length(individual)==1)
    if (! (is.logical(individual) && (!individual))) {
      indiv.logi <- TRUE
      if (is.numeric(individual))   indiv.which <- which(is.element(1:k, individual))
      if (is.character(individual)) indiv.which <- which(is.element(labels, match.arg(individual,labels)))
      if (length(indiv.which)==0) warning("cannot make sense of 'individual' argument: empty subset.")
    }
    else indiv.logi <- FALSE
    if (all(is.na(mu.p))) { # compute tau quantile
      stopifnot(all(tau.p>=0), all(tau.p<=1))
      if (predict) warning("'predict' argument ignored!")
      if (indiv.logi) warning("'individual' argument ignored!")
      upper <- 1
      if (any((tau.p<1) & (tau.p>0))) {
        maxp <- max(tau.p[tau.p<1])
        while (pposterior(tau=upper) < maxp) upper <- upper * 2
      }
      qfun <- function(p)
      {
        stopifnot(p>=0, p<=1)
        if (p==0) quant <- 0
        else {
          if (p==1) quant <- Inf
          else
            quant <- uniroot(function(xx){return(pposterior(tau=xx)-p)},
                             interval=c(0,upper), tol=tol.uniroot)$root
        }
        return(quant)
      }
      result <- apply(matrix(tau.p,ncol=1),1,qfun)
    }
    else if (all(is.na(tau.p))) { # compute mu quantile
      if (indiv.logi && any(zerosigma) && (indiv.which == which(zerosigma))){
        result <- rep(NA, length(mu.p))
        result[is.finite(mu.p)] <- y[zerosigma]
      }
      stopifnot(all(mu.p>=0), all(mu.p<=1))
      if (any((mu.p<1) & (mu.p>0))) {
        minp <- min(c(mu.p[mu.p>0], 1-mu.p[mu.p<1]))
        # derive minimum/maximum based on variance and Chebychev inequality:
        if (predict) {
          lower <- sumstats["mean","theta"] - sqrt(1/minp)*sumstats["sd","theta"]*1.1
          upper <- sumstats["mean","theta"] + sqrt(1/minp)*sumstats["sd","theta"]*1.1
        }
        else if (indiv.logi) {
          lower <- y[indiv.which] - sqrt(1/minp)*sigma[indiv.which]*1.1
          upper <- y[indiv.which] + sqrt(1/minp)*sigma[indiv.which]*1.1
        }
        else {
          lower <- sumstats["mean","mu"] - sqrt(1/minp)*sumstats["sd","mu"]*1.1
          upper <- sumstats["mean","mu"] + sqrt(1/minp)*sumstats["sd","mu"]*1.1
        }
        while (pposterior(mu=lower,predict=predict,individual=individual) > minp)
            lower <- lower-sumstats["sd","mu"]
        while (pposterior(mu=upper,predict=predict,individual=individual) < (1-minp))
            upper <- upper+sumstats["sd","mu"]
      }
      qfun <- function(p)
      {
        stopifnot(p>=0, p<=1)
        if (p==0) quant <- -Inf
        else {
          if (p==1) quant <- Inf
          else
            quant <- uniroot(function(yy){return(pposterior(mu=yy,predict=predict,individual=individual)-p)},
                             interval=c(lower,upper), tol=tol.uniroot)$root
        }
        return(quant)
      }
      result <- apply(matrix(mu.p,ncol=1),1,qfun)
    }
    else {
      warning("need to supply EITHER tau.p OR mu.p, not both.")
      result <- NA
    }
    return(result)
  }

  rposterior <- function(n=1, predict=FALSE, individual=FALSE, tau.sample=TRUE)
  # posterior random number generation for tau and mu
  {
    stopifnot(n>0, n==round(n), length(individual)==1,
              !is.logical(individual) || !individual)
    if (tau.sample) {  # draw joint, bivariate (tau,mu) pairs:
      samp <- matrix(NA, nrow=n, ncol=2, dimnames=list(NULL,c("tau","mu")))
      if (is.numeric(individual) | is.character(individual))
          colnames(samp)[2] <- "theta"
      u <- runif(n=n)
      samp[,"tau"] <- apply(matrix(u,ncol=1), 1, function(x){return(qposterior(tau.p=x))})
      cond.sample <- function(t)
      {
        cm <- conditionalmoment(t, predict=predict, individual=individual)
        return(rnorm(n=1, mean=cm[1], sd=cm[2]))
      }
      samp[,2] <- apply(matrix(samp[,"tau"],ncol=1), 1, cond.sample)
    } else {           # draw marginal, univariate (mu or theta) numbers:
      samp <- rep(NA, n)
      if (!predict & (is.logical(individual) && (!individual)))
        meansd <- support[,c("mean","sd")]
      else
        meansd <- conditionalmoment(support[,"tau"], predict=predict, individual=individual)
      for (i in 1:n) {
        j <- sample(1:nrow(support), 1, prob=support[,"weight"])
        samp[i] <- rnorm(1, mean=meansd[j,"mean"], sd=meansd[j,"sd"])
      }          
    }
    return(samp)
  }

  post.interval <- function(tau.level=NA, mu.level=NA, theta.level=mu.level,
                            method=interval.type,
                            predict=FALSE, individual=FALSE)
  # determine credibility/prediction/shrinkage interval (tau, mu or theta)
  {
    if (is.na(mu.level)) mu.level <- theta.level
    stopifnot((is.finite(tau.level) & ((tau.level>0) & (tau.level<1)))
              | (is.finite(mu.level) & ((mu.level>0) & (mu.level<1))))
    stopifnot(length(individual)==1)
    method <- match.arg(method, c("shortest","central","evidentiary"))
    if (all(is.na(mu.level))) {      # CI for heterogeneity tau:
      if (predict) warning("'predict' argument ignored!")
      if (! (is.logical(individual) && (!individual))) warning("'individual' argument ignored!")
      if (method=="central")         # central interval
        result <- qposterior(tau.p=c((1-tau.level)/2, 1-(1-tau.level)/2))
      else if (method=="shortest") { # shortest interval
        intwidth <- function(left)
        {
          pleft <- pposterior(tau=left)
          right <- qposterior(tau=tau.level+pleft)
          return(right-left)
        }
        opti <- optimize(intwidth, lower=0, upper=qposterior(tau=1-tau.level))$minimum
        # catch marginal minimum:
        if (intwidth(0) < intwidth(opti))
          result <- c(0, qposterior(tau=tau.level))
        else
          result <- c(opti, qposterior(tau=tau.level+pposterior(tau=opti)))
      }
      else {                         # evidentiary interval (tau)
        expectedLogLikeli <- function(left)
        {
          pleft <- pposterior(tau=left)
          right <- qposterior(tau=tau.level+pleft)
          expect <- integrate(function(x){return(likelihood(tau=x,log=TRUE)*dposterior(tau=x))},
                              lower=left, upper=right,
                              rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)
          return(expect$value)
        }
        opti <- optimize(expectedLogLikeli, maximum=TRUE,
                         lower=0, upper=qposterior(tau=1-tau.level))$maximum
        # catch marginal minimum:
        if (expectedLogLikeli(0) > expectedLogLikeli(opti))
          result <- c(0, qposterior(tau=tau.level))
        else
          result <- c(opti, qposterior(tau=tau.level+pposterior(tau=opti)))          
      }
    }
    else if (all(is.na(tau.level))) {  # CI for effect mu:
      if (method=="central")           # central interval
        result <- qposterior(mu.p=c((1-mu.level)/2, 1-(1-mu.level)/2), predict=predict, individual=individual)
      else if (method=="shortest") {   # shortest interval
        intwidth <- function(left)
        {
          pleft <- pposterior(mu=left, predict=predict, individual=individual)
          right <- qposterior(mu=mu.level+pleft, predict=predict, individual=individual)
          return(right-left)
        }
        opti <- optimize(intwidth,
                         lower=qposterior(mu=(1-mu.level)/50, predict=predict, individual=individual),
                         upper=qposterior(mu=1-mu.level, predict=predict, individual=individual))$minimum
        result <- c(opti, qposterior(mu=mu.level+pposterior(mu=opti, predict=predict, individual=individual), predict=predict, individual=individual))
      }
      else {                           # evidentiary interval (mu)
        if (predict) {
          warning("evidentiary prediction intervals are not implemented!")
          result <- c(NA,NA)
        }
        else if ((length(individual)>1) || (!is.logical(individual)) || individual) {
          warning("evidentiary shrinkage intervals are not implemented!")
          result <- c(NA,NA)
        } else {
          expectedLogLikeli <- function(left)
          {
            pleft <- pposterior(mu=left)
            right <- qposterior(mu=mu.level+pleft)
            expect <- integrate(function(x){return(likelihood(mu=x)*dposterior(mu=x))},
                                lower=left, upper=right,
                                rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)
            return(expect$value)
          }
          opti <- optimize(expectedLogLikeli, maximum=TRUE,
                           lower=qposterior(mu=(1-mu.level)/100), upper=qposterior(mu=1-mu.level))$maximum
          result <- c(opti, qposterior(mu=mu.level+pposterior(mu=opti)))
        }
      }
    }
    else {
      warning("need to supply EITHER tau.level OR mu.level, not both.")
      result <- c(NA,NA)
    }
    attr(result, "interval.type") <- method
    return(result)
  }
  
  conditionalmoment <- function(tau, predict=FALSE, individual=FALSE, simplify=TRUE)
  # compute conditional posterior moments (mean, sd) of mu for given value(s) of tau
  # if (predict==TRUE), moments of the /posterior predictive distribution/ are returned.
  # if (individual==TRUE), moments of the conditional posterior /of the study-specific effects/
  # (theta[i]) are returned.
  {
    # interpret the "individual" argument, derive a "logical" and "numeric" component:
    if (all(is.logical(individual))) { #  (logical "individual" argument)
      indiv.logi  <- individual[1]
      if (indiv.logi)
        indiv.which <- 1:k
      else
        indiv.which <- NA          
    }
    else { #  (numerical/character "individual" argument)
      indiv.logi  <- TRUE
      if (all(is.numeric(individual))) # number(s) provided
        indiv.which <- which(is.element(1:k, individual))
      else if (all(is.character(individual))) # label(s) provided
        indiv.which <- which(is.element(labels, match.arg(individual, labels, several.ok=TRUE)))
      else warning("cannot make sense of 'individual' argument: funny format.")
      if (length(indiv.which)==0)
        warning("cannot make sense of 'individual' argument: empty subset.")
    }
    if (predict & indiv.logi) warning("need to specify either 'predict' or 'individual' argument, but not both.")
    cm <- function(t)
    {
      if ((t==0) && any(zerosigma)) {
        conditionalmean <- y[zerosigma]
        conditionalvar  <- 0.0
      } else {
        if (is.na(mu.prior.mean)) { # uniform prior
          yTerm   <- y
          varTerm <- t^2 + sigma^2
        }
        else {                      # conjugate normal prior
          yTerm   <- c(mu.prior.mean, y)
          varTerm <- c(mu.prior.sd^2, t^2 + sigma^2)
        }
        conditionalvar  <- 1 / sum(1/varTerm)
        conditionalmean <- sum(yTerm/varTerm) * conditionalvar
      }
      return(c("mean"=conditionalmean,"sd"=sqrt(conditionalvar)))
    }
    # compute conditional moments of mu
    # (or predictive, if requested):
    if (!indiv.logi) {
      result <- t(apply(matrix(tau,ncol=1),1,cm))
      if (predict) result[,"sd"] <- sqrt(result[,"sd"]^2 + tau^2)
    }
    # compute conditional moments of individual, study-specific means:
    else {
      result <- array(NA, dim=c(length(tau), 2, length(indiv.which)),
                      dimnames=list(NULL, c("mean","sd"), labels[indiv.which]))
      musigma <- t(apply(matrix(tau,ncol=1),1,cm))
      # loop over estimates:
      zero <- (tau==0)
      for (i in 1:length(indiv.which)) {
        if (any(!zero)) { # tau > 0
          if (sigma[indiv.which[i]] > 0.0) {
            result[!zero,"mean",i] <- (y[indiv.which[i]]/sigma[indiv.which[i]]^2 + musigma[!zero,"mean"]/tau[!zero]^2) / (1/sigma[indiv.which[i]]^2+1/tau[!zero]^2)
            result[!zero,"sd",i]   <- sqrt(1 / (1/sigma[indiv.which[i]]^2+1/tau[!zero]^2)
                                           + (musigma[!zero,"sd"] * (1/tau[!zero]^2) / (1/sigma[indiv.which[i]]^2+1/tau[!zero]^2))^2)
          } else {
            result[!zero,"mean",i] <- y[indiv.which[i]]
            result[!zero,"sd",i]   <- 0.0
          }
        }
        if (any(zero)) {  # tau = 0
          result[zero,"mean",i] <- musigma[zero,"mean"]
          result[zero,"sd",i]   <- musigma[zero,"sd"]
        }
      }
      # simplify array to matrix (if possible and desired):
      if ((dim(result)[3] == 1) && (simplify)) result <- result[,,1]
    }
    return(result)
  }

  ISquared <- function(tau)
  # compute heterogeneity measure "I-squared" as a function of tau
  {
    I2 <- rep(NA, length(tau))
    I2[tau>=0] <- tau[tau>=0]^2 / (tau[tau>=0]^2 + sigma2hat)
    return(I2)
  }

  discretize <- function(delta=0.01, epsilon=0.0001, alldivs=TRUE)
  # discretize parameter space in order to approximate marginal (mixture) distribution
  # via a finite number of mixture components; see also:
  #   C. Roever, T. Friede.
  #   Discrete approximation of a mixture distribution via restricted divergence.
  #   Journal of Computational and Graphical Statistics, 26(1): 217-222, 2017.
  #   http://doi.org/10.1080/10618600.2016.1276840
  {
    symKL <- function(mean1, sd1, mean2, sd2)
    # (general) symmetrized KL-divergence of two normal distributions
    {
      stopifnot(sd1>0, sd2>0)
      return((mean1-mean2)^2 * 0.5 * (1/sd1^2 + 1/sd2^2) + (sd1^2-sd2^2)^2 / (2*sd1^2*sd2^2))
    }
    divergence <- function(tau1, tau2)
    # (symmetrized) divergence b/w two conditionals specified through corresponding tau values
    {
      if (!alldivs) { # evaluate divergence based on (conditional) mu posterior only:
        cm <- conditionalmoment(tau=c(tau1, tau2))
        div <- symKL(cm[1,"mean"], cm[1,"sd"], cm[2,"mean"], cm[2,"sd"])
      }
      else { # evaluate divergence based also on (conditional) individual-study
             # and predictive mu posterior, and determine the maximum divergence:
        # determine array of ALL conditional moments:
        cm <- array(c(conditionalmoment(tau=c(tau1, tau2)),
                      conditionalmoment(tau=c(tau1, tau2), predict=TRUE),
                      #conditionalmoment(tau=c(tau1, tau2), individual=TRUE)),
                      conditionalmoment(tau=c(tau1, tau2), individual=(1:k)[!zerosigma])),
                    #dim=c(2, 2, k+2))
                    dim=c(2, 2, sum(!zerosigma)+2))
        # determine individual divergences and their maximum:
        #div <- rep(NA, k+2)
        div <- rep(NA, sum(!zerosigma)+2)
        #for (i in 1:(k+2))
        for (i in 1:(sum(!zerosigma)+2))
          div[i] <- symKL(cm[1,1,i], cm[1,2,i], cm[2,1,i], cm[2,2,i])
        div <- max(div)
      }
      return(div)
    }
    # determine range of interest / upper bound on tau:
    if (any(zerosigma)) maxtau <- qposterior(tau.p=1-min(c(epsilon/2, 1-epsilon/2)))
    else                maxtau <- qposterior(tau.p=1-min(c(epsilon, 1-epsilon)))
    # special procedure for 1st bin
    # (want zero to be first reference point UNLESS one of the sigma[i] is zero):
    if (any(zerosigma)) tau <- qposterior(tau.p=min(c(epsilon/2, 1-epsilon/2)))
    else                tau <- 0.0
    # search for upper bin margin:
    upper <- 1
    diverg <- divergence(tau, upper)
    while (diverg <= delta) {
      upper <- upper*2
      diverg <- divergence(tau, upper)
    }
    ur <- uniroot(function(t){return(divergence(tau, t)-delta)},
                  lower=tau, upper=upper,
                  f.lower=-delta, f.upper=diverg-delta,
                  tol=tol.uniroot)
    prob1 <- 0.0
    prob2 <- pposterior(tau=ur$root)
    # store result for 1st bin:
    result <- matrix(c(tau, prob2-prob1),
                     nrow=1, ncol=2,
                     dimnames=list(NULL,c("tau","weight")))
    tau <- ur$root
    # determine following bins (2,...):
    bin <- 2
    while ((tau <= maxtau) | (bin<=2)) {  # (at least 2 support points)
      result <- rbind(result, rep(0,2))
      # determine bin's reference point:
      diverg <- divergence(tau, upper)
      while (diverg <= delta) {
        upper <- upper*2
        diverg <- divergence(tau, upper)
      }
      ur <- uniroot(function(t){return(divergence(tau, t)-delta)},
                    lower=tau, upper=upper,
                    f.lower=-delta, f.upper=diverg-delta,
                    tol=tol.uniroot)
      tau <- ur$root
      result[bin,"tau"] <- tau
      # determine bin's upper bound:
      diverg <- divergence(tau, upper)
      while (diverg <= delta) {
        upper <- upper*2
        diverg <- divergence(tau, upper)
      }
      ur <- uniroot(function(t){return(divergence(tau, t)-delta)},
                    lower=tau, upper=upper,
                    f.lower=-delta, f.upper=diverg-delta,
                    tol=tol.uniroot)
      tau <- ur$root
      # determine bin's weight:
      prob1 <- prob2
      prob2 <- pposterior(tau=tau)
      result[bin,"weight"] <- max(c(0.0, prob2-prob1))
      # sanity check (to catch possible uniroot() failure):
      if (result[bin,"tau"] <= result[bin-1,"tau"]) {
        warning("DIRECT grid setup seems to have failed.")
        tau <- maxtau + 1.0
      }
      bin <- bin+1
    }
    # re-normalize weights (if necessary):
    sw <- sum(result[,"weight"])
    if (sw <= 0) {
      result[,"weight"] <- rep(1/nrow(result), nrow(result))
      sw <- 1.0
    }
    if (sw != 1.0)
      result[,"weight"] <- result[,"weight"] / sw
    return(result)
  }

  accelerate <- FALSE  # flag to speed up computations at the cost of a few estimates and a little accuracy

  # compute set of tau support points & corresponding weights:
  support <- discretize(delta=delta, epsilon=epsilon, alldivs=ifelse(accelerate, FALSE, TRUE))
  rm(list=c("discretize"))
  #if (nrow(support) <= 10)
  #  warning(paste0("Discretization of heterogeneity parameter space yielded very few (only ",nrow(support),") support points.\n",
  #                 "This *may* indicate numerical problems, e.g. due to a very narrow heterogeneity prior. Please double-check."))
  support <- cbind(support,
                   conditionalmoment(support[,"tau"]),
                   "sd.pred"=conditionalmoment(support[,"tau"],predict=TRUE)[,"sd"])
  
  # compute tau posterior's summary statistics:
  sumstats <- matrix(NA, nrow=6, ncol=3,
                     dimnames=list(c("mode", "median", "mean","sd", "95% lower", "95% upper"),
                                   c("tau","mu","theta")))
  expectation <- try(integrate(function(x)return(dposterior(x)*x), lower=0, upper=Inf,
                               rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)$value, silent=TRUE)
  if (inherits(expectation, "try-error")) {
    expectation <- NA
    variance <- NA
  }
  else {
    variance <- try(integrate(function(x)return(dposterior(x)*(x-expectation)^2), lower=0, upper=Inf,
                              rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)$value,
                    silent=TRUE)
    if (inherits(variance, "try-error"))
      variance <- NA
  }
  sumstats[c("mean","sd"),"tau"] <- c(expectation, sqrt(variance))
  sumstats["median","tau"] <- qposterior(0.5)
  sumstats[c("95% lower", "95% upper"),"tau"] <- post.interval(tau=0.95, method=ifelse(accelerate, "central", interval.type))
  if (!accelerate) {
    if (! is.finite(dposterior(tau=0))) {
      sumstats["mode","tau"] <- 0
    } else {
      maxi <- optimize(function(x){return(dposterior(tau=x))},
                       lower=0, upper=qposterior(tau=0.9), maximum=TRUE)
      if ((maxi$maximum <= 2 * .Machine$double.eps^0.25)
          && (dposterior(0) > maxi$objective)) maxi <- list("maximum"=0)
      sumstats["mode","tau"] <- maxi$maximum
      rm(list=c("maxi"))
    }
  }
  rm(list=c("expectation", "variance"))
  
  # compute mu posterior's summary statistics:
  if (all(is.finite(support))) {
    sumstats["mean","mu"] <- sum(support[,"mean"] * support[,"weight"])
    sumstats["sd","mu"] <- sqrt(sum(((support[,"mean"]-sumstats["mean","mu"])^2 + support[,"sd"]^2) * support[,"weight"]))
    sumstats["median","mu"] <- qposterior(mu=0.5)
    sumstats[c("95% lower", "95% upper"),"mu"] <- post.interval(mu=0.95, method=ifelse(accelerate, "central", interval.type))
    if (!accelerate) {
      maxi <- optimize(function(x){return(dposterior(mu=x))},
                       lower=qposterior(mu=0.1), upper=qposterior(mu=0.9), maximum=TRUE)
      sumstats["mode","mu"] <- maxi$maximum
      rm(list=c("maxi"))
    }
  }

  # compute mu posterior-predictive distribution's summary statistics:
  if (all(is.finite(support))) {
    sumstats["mean","theta"] <- sumstats["mean","mu"]
    sumstats["sd","theta"] <- sqrt(sum(((support[,"mean"]-sumstats["mean","theta"])^2 + support[,"sd.pred"]^2) * support[,"weight"]))
    if (!accelerate) {
      sumstats["median","theta"] <- qposterior(mu=0.5, predict=TRUE)
      sumstats[c("95% lower", "95% upper"),"theta"] <- post.interval(mu=0.95, predict=TRUE)
      maxi <- optimize(function(x){return(dposterior(mu=x, predict=TRUE))},
                       lower=qposterior(mu=0.1, predict=TRUE), upper=qposterior(mu=0.9, predict=TRUE), maximum=TRUE)
      sumstats["mode","theta"] <- maxi$maximum
      rm(list=c("maxi"))
    }
  }

  # compute joint & marginal maximum-likelihood (ML) estimates:
  ml.estimate <- rbind("joint"=c("tau"=NA, "mu"=NA), "marginal"=c("tau"=NA, "mu"=NA))
  map.estimate <- rbind("joint"=c("tau"=NA, "mu"=NA),
                        "marginal"=c("tau"=sumstats["mode","tau"],
                                     "mu"=sumstats["mode","mu"]))
  if (!accelerate) {
    # search for joint ML:
    opti <- optim(sumstats["median",c("tau","mu")],
                  function(x){return(-likelihood(tau=x[1], mu=x[2], log=TRUE))})
    # possibly catch optimum at (tau=0) parameter space margin:
    if (any(zerosigma)) {
      ml.estimate["joint",] <- opti$par
    } else {
      upper <- opti$par[2] + 1
      while (likelihood(tau=0, mu=upper+1, log=TRUE) > likelihood(tau=0, mu=upper, log=TRUE))
        upper <- upper + 1
      lower <- opti$par[2] - 1
      while (likelihood(tau=0, mu=lower-1, log=TRUE) > likelihood(tau=0, mu=lower, log=TRUE))
        lower <- lower - 1
      maxi <- optimize(function(x){return(likelihood(tau=0, mu=x, log=TRUE))},
                       lower=lower, upper=upper, maximum=TRUE)
      if (likelihood(tau=opti$par[1], mu=opti$par[2], log=TRUE) > maxi$objective)
        ml.estimate["joint",] <- opti$par
      else
        ml.estimate["joint",] <- c(0, maxi$maximum)
    }
    # marginal ML (mu):
    upper <- ml.estimate["joint","mu"] + 1
    while (likelihood(mu=upper+1) > likelihood(mu=upper))
      upper <- upper + 1
    lower <- ml.estimate["joint","mu"] - 1
    while (likelihood(mu=lower-1) > likelihood(mu=lower))
      lower <- lower - 1
    maxi <- optimize(function(x){return(likelihood(mu=x))},
                     lower=lower, upper=upper, maximum=TRUE)
    ml.estimate["marginal","mu"] <- maxi$maximum
    
    # marginal ML (tau):
    upper <- ifelse(ml.estimate["joint","tau"] > 0, 2*ml.estimate["joint","tau"], 1.0)
    while (likelihood(tau=2*upper) > likelihood(tau=upper))
      upper <- upper * 2
    maxi <- optimize(function(x){return(likelihood(tau=x))},
                     lower=0, upper=upper, maximum=TRUE)
    if (all(!zerosigma) && (likelihood(tau=0) > maxi$objective))
      ml.estimate["marginal","tau"] <- 0.0
    else
      ml.estimate["marginal","tau"] <- maxi$maximum
    
    # compute joint maximum-a-posteriori (MAP) estimate:
    opti <- optim(sumstats["median",c("tau","mu")],
                  function(x){return(-dposterior(tau=x[1], mu=x[2], log=TRUE))})
    map.value <- dposterior(tau=opti$par[1], mu=opti$par[2])
    map.estimate["joint",] <- opti$par
      
    # possibly catch optimum at (tau=0) parameter space margin:
    if (all(!zerosigma)) {
      if (is.finite(tau.prior(0))) {
        upper <- sumstats["95% upper", "mu"]
        while (dposterior(mu=upper+1, tau=0) > dposterior(mu=upper, tau=0))
          upper <- upper + 1
        lower <- sumstats["95% lower", "mu"]
        while (dposterior(mu=lower-1, tau=0) > dposterior(mu=lower, tau=0))
          lower <- lower - 1
        maxi <- optimize(function(x){return(dposterior(mu=x, tau=0))},
                         lower=lower, upper=upper, maximum=TRUE)
        if (maxi$objective > map.value) {
          map.estimate["joint",] <- c(0, maxi$maximum)
          map.value <- maxi$objective
        }
      } else {
        map.value <- Inf
        map.estimate["joint","tau"] <- 0.0
      }
    }
      
    # possibly catch diverging posterior density
    if (! is.finite(map.value)) {
      map.estimate["joint","mu"] <- conditionalmoment(tau=map.estimate["joint","tau"])[1,"mean"]
    }

    # clean up:
    rm(list=c("opti","maxi","lower","upper","map.value"))
  }
  # compute "shrinkage" estimates of theta[i]:
  shrink <- matrix(NA, nrow=8, ncol=k,
                   dimnames=list(c("y","sigma","mode", "median", "mean","sd", "95% lower", "95% upper"),
                                 labels))
  shrink["y",] <- y
  shrink["sigma",] <- sigma
  if (all(is.finite(support)) & (!accelerate)) {
    for (i in 1:k) {
      if (zerosigma[i]) {
        shrink[c("mean", "median", "mode", "95% lower", "95% upper"), i] <- y[i]
        shrink["sd", i] <- 0.0
      } else {
        musigma <- conditionalmoment(support[,"tau"], individual=i)
        shrink["mean",i] <- sum(musigma[,"mean"] * support[,"weight"])
        shrink["sd",i]   <- sqrt(sum(((musigma[,"mean"]-shrink["mean",i])^2 + musigma[,"sd"]^2) * support[,"weight"]))
        shrink["median",i] <- qposterior(mu=0.5, individual=i)
        shrink[c("95% lower", "95% upper"),i] <- post.interval(mu=0.95, individual=i)
        maxi <- optimize(function(x){return(dposterior(mu=x, individual=i))},
                         lower=qposterior(mu=0.1), upper=qposterior(mu=0.9), maximum=TRUE)
        shrink["mode",i] <- maxi$maximum
      }
    }
  }
  
  # compute marginal likelihood & Bayes factors:
  marglik <- NA_real_
  bayesfactor <- matrix(NA_real_, nrow=2, ncol=2, dimnames=list(c("actual", "minimum"), c("tau=0","mu=0")))
  if (tau.prior.proper) {
    bayesfactor["minimum","mu=0"]  <- likelihood(mu=0) / likelihood(mu=ml.estimate["marginal","mu"])
  }
  # check for proper effect prior:
  if (is.finite(mu.prior.mean) 
      && is.finite(mu.prior.sd)
      && (mu.prior.sd > 0)) {
    bayesfactor["minimum","tau=0"] <- likelihood(tau=0) / likelihood(tau=ml.estimate["marginal","tau"])
    # check for proper heterogeneity prior:
    if (tau.prior.proper) {
      marglik.int <- integrate(function(t){return(likelihood(tau=t)*dprior(tau=t))},
                               rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate,
                               lower=0, upper=Inf)
      if (marglik.int$message == "OK") {
        marglik <- marglik.int$value
        bayesfactor["actual", "tau=0"] <- likelihood(tau=0) / marglik
        bayesfactor["actual", "mu=0"]  <- likelihood(mu=0) / marglik
      } else {
        attr(marglik, "NA.reason") <- "failed computing marginal likelihood"
      }
      rm("marglik.int")
    } else {
      attr(marglik, "NA.reason") <- "improper heterogeneity prior"
    }
  } else {
    attr(marglik, "NA.reason") <- "improper effect prior"
  }

  # compute posterior mean inverse-variance weights:
  ivweights <- function(tau, idx)
  # inverse-variance weights as a function of tau
  {
    stopifnot(length(tau)==1, tau>=0)
    if (is.finite(mu.prior.mean) && is.finite(mu.prior.sd)) {
      weights <- 1 / c(tau^2 + sigma^2, mu.prior.sd^2)
      names(weights) <- c(labels, "prior mean")
    } else {
      weights <- 1 / (tau^2 + sigma^2)
      names(weights) <- labels
    }
    weights <- weights / sum(weights)
    return(weights[idx])
  }
  ivweights <- Vectorize(ivweights, vectorize.args="tau")
  
  # compute posterior means:
  if (is.finite(mu.prior.mean) && is.finite(mu.prior.sd)) {
    meanweights <- rep(NA_real_, k+1)
    names(meanweights) <- c(labels, "prior mean")
  } else {
    meanweights <- rep(NA_real_, k)
    names(meanweights) <- labels
  }
  for (i in 1:length(meanweights)) {
    meanweights[i] <- integrate(function(t){ivweights(tau=t, idx=i) * dposterior(tau=t)},
                                lower=0, upper=Inf,
                                rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)$value
  }

  # compute "shrinkage weights":
  shrinkweights <- function(tau, study.idx, shrink.idx)
  # "study.idx"  :  which study's weight
  # "shrink.idx" :  which shrinkage estimate
  {
    stopifnot(length(tau)==1, tau>=0)
    # compute inverse variance weights:
    if (is.finite(mu.prior.mean) && is.finite(mu.prior.sd)) {
      ivw <- 1 / c(tau^2 + sigma^2, mu.prior.sd^2)
      names(ivw) <- c(labels, "prior mean")
    } else {
      ivw <- 1 / (tau^2 + sigma^2)
      names(ivw) <- labels
    }
    ivw <- ivw / sum(ivw)
    if (any(sigma==0)) { # catch zero sigma case:
      s0 <- which(sigma==0)[1]
      for (i in 1:length(ivw)) ivw[i] <- 0.0
      ivw[s0] <- 1.0
    }
    # compute shrinkage weights:
    if (tau > 0)
      swt <- (sigma[shrink.idx]^(-2)) / (sigma[shrink.idx]^(-2) + tau^(-2))
    else
      swt <- 0.0
    indic <- rep(0.0, length(ivw))
    indic[shrink.idx] <- 1.0
    swt <- swt*indic + (1-swt)*ivw
    swt <- swt / sum(swt)
    if (any(sigma==0)) { # catch zero sigma case:
      s0 <- which(sigma==0)[1]
      swt <- rep(0.0, length(ivw))
      swt[s0] <- 1.0
    }
    return(swt[study.idx])
  }
  shrinkweights <- Vectorize(shrinkweights, vectorize.args="tau")
  shrinkageweights <- matrix(NA_real_, nrow=length(meanweights), ncol=k,
                             dimnames=list("study"=names(meanweights),
                                           "shrinkage estimate"=labels))
  for (i in 1:nrow(shrinkageweights)) {
    for (j in 1:ncol(shrinkageweights)) {
      shrinkageweights[i,j] <- integrate(function(t){shrinkweights(tau=t,i,j)*dposterior(tau=t)},
                                         lower=0.0, upper=Inf,
                                         rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)$value
    }
  }
  
  ptm <- proc.time()[1] - ptm[1]
  
  muprior <- c(mu.prior.mean, mu.prior.sd)
  names(muprior) <- c("mean","sd")
  
  # assemble eventual result to be returned:
  result <- list("y"                   = y,
                 "sigma"               = sigma,
                 "labels"              = labels,
                 "k"                   = k,
                 "tau.prior"           = tau.prior,
                 "mu.prior"            = muprior,
                 "dprior"              = dprior,
                 "tau.prior.proper"    = tau.prior.proper,
                 "likelihood"          = likelihood,
                 "dposterior"          = dposterior,
                 "pposterior"          = pposterior,
                 "qposterior"          = qposterior,
                 "rposterior"          = rposterior,
                 "post.interval"       = post.interval,
                 "cond.moment"         = conditionalmoment,
                 "I2"                  = ISquared,
                 "summary"             = sumstats,
                 "interval.type"       = interval.type,
                 "ML"                  = ml.estimate,
                 "MAP"                 = map.estimate,
                 "theta"               = shrink,
                 "weights"             = meanweights,
                 "weights.theta"       = shrinkageweights,
                 "marginal.likelihood" = marglik,
                 "bayesfactor"         = bayesfactor,
                 "support"             = support,
                 "delta"               = delta,
                 "epsilon"             = epsilon,
                 "rel.tol.integrate"   = rel.tol.integrate,
                 "abs.tol.integrate"   = abs.tol.integrate,
                 "tol.uniroot"         = tol.uniroot,
                 "call"                = match.call(expand.dots=FALSE),
                 "init.time"           = c("seconds"=unname(ptm)))
  class(result) <- "bayesmeta"
  return(result)
}


bayesmeta.escalc <- function(y, labels=NULL, ...)
# apply "bayesmeta()" to an "escalc" object
{
  attri <- attributes(y)
  if (all(is.element(c("yi.names", "vi.names"), names(attri)))) { # (for recent "metafor" versions)
    var.names <- c(attri$yi.names, attri$vi.names)
  } else if (is.element("var.names", names(attri))) {             # (for older "metafor" versions)
    var.names <- attri$var.names
  } else {
    stop(paste("Cannont extract \"yi.names\" and \"vi.names\" (or \"var.names\") attribute(s) from escalc object ",
               "(check use of the \"var.names\" option).", sep=""))
  }
  stopifnot(length(var.names)==2, all(is.character(var.names)))
  if (!all(is.element(var.names, names(y)))) {
    stop(paste("Cannont find columns \"",
               var.names[1],"\" and/or \"",
               var.names[2],"\" in escalc object ",
               "(check use of the \"var.names\" option).", sep=""))
  }
  if (is.null(labels)) {
    if (is.element("slab", names(attributes(y[,var.names[1]]))))
      labels <- as.character(attr(y[,var.names[1]], "slab"))
  }
  result <- bayesmeta.default(y=as.vector(y[,var.names[1]]),
                              sigma=sqrt(as.vector(y[,var.names[2]])),
                              labels=labels, ...)
  result$call <- match.call(expand.dots=FALSE)
  return(result)
}


print.bayesmeta <- function(x,...)
# print a short summary
{
  cat(" 'bayesmeta' object.\n")
  cat(paste("\n",x$k," estimates:\n", sep=""))
  if (length(x$labels)>10)
    cat(paste(c(x$labels[1:10], "..."), collapse=", "))
  else
    cat(paste(x$labels[1:length(x$labels)], collapse=", "))
  cat("\n")
  cat(paste0("\ntau prior (", ifelse(x$tau.prior.proper, "proper", "improper"), "):\n"))
  if (is.element("bayesmeta.label", names(attributes(x$tau.prior))))
    cat(paste(attributes(x$tau.prior)[["bayesmeta.label"]], "\n"))
  else
    print(x$tau.prior)
  mu.prior.proper <- (is.finite(x$mu.prior["mean"]) && is.finite(x$mu.prior["sd"]) && (x$mu.prior["sd"] > 0))
  cat(paste0("\nmu prior (", ifelse(mu.prior.proper, "proper", "improper"), "):\n"))
  if (mu.prior.proper)
    cat(paste("normal(mean=",x$mu.prior["mean"],", sd=",x$mu.prior["sd"],")\n",sep=""))
  else
    cat("uniform(min=-Inf, max=Inf)\n")
  cat("\nML and MAP estimates:\n")
  print(rbind("ML joint"=x$ML["joint",],
              "ML marginal"=x$ML["marginal",],
              "MAP joint"=x$MAP["joint",],
              "MAP marginal"=x$MAP["marginal",]))
  cat("\nmarginal posterior summary:\n")
  print(x$summary[,c("tau","mu")])
  if (x$interval.type=="shortest")
    cat("\n(quoted intervals are shortest credible intervals.)\n")
  else
    cat("\n(quoted intervals are central, equal-tailed credible intervals.)\n")      
  invisible(x)
}


summary.bayesmeta <- function(object,...)
# print a longer summary
{
  cat(" 'bayesmeta' object.\n")
  data <- matrix(c(object$y, object$sigma), ncol=2, dimnames=list(object$labels, c("y","sigma")))
  cat(paste("data (", object$k, " estimates):\n",sep=""))
  if (nrow(data)<=20)
    print(data)
  else {
    print(data[1:20,])
    cat(paste(" [...]  (truncated;", object$k, "estimates total)\n"))
  }
  cat(paste0("\ntau prior (", ifelse(object$tau.prior.proper, "proper", "improper"), "):\n"))
  if (is.element("bayesmeta.label", names(attributes(object$tau.prior))))
    cat(paste(attributes(object$tau.prior)[["bayesmeta.label"]], "\n"))
  else
    print(object$tau.prior)
  mu.prior.proper <- (is.finite(object$mu.prior["mean"]) && is.finite(object$mu.prior["sd"]) && (object$mu.prior["sd"] > 0))
  cat(paste0("\nmu prior (", ifelse(mu.prior.proper, "proper", "improper"), "):\n"))
  if (mu.prior.proper)
    cat(paste("normal(mean=",object$mu.prior["mean"],", sd=",object$mu.prior["sd"],")\n",sep=""))
  else
    cat("uniform(min=-Inf, max=Inf)\n")
  cat("\nML and MAP estimates:\n")
  print(rbind("ML joint"=object$ML["joint",],
              "ML marginal"=object$ML["marginal",],
              "MAP joint"=object$MAP["joint",],
              "MAP marginal"=object$MAP["marginal",]))
  cat("\nmarginal posterior summary:\n")
  print(object$summary)
  if (object$interval.type=="shortest")
    cat("\n(quoted intervals are shortest credible intervals.)\n")
  else
    cat("\n(quoted intervals are central, equal-tailed credible intervals.)\n")      
  if (any(is.finite(object$bayesfactor))) {
    cat("\nBayes factors:\n")
    print(object$bayesfactor)
  }
  cat("\nrelative heterogeneity I^2 (posterior median):", object$I2(tau=object$summary["median","tau"]), "\n")
  invisible(object)
}


plot.bayesmeta <- function(x, main=deparse(substitute(x)),
                           which=1:4, prior=FALSE,
                           forest.margin=8,
                           mulim=c(NA,NA), taulim=c(NA,NA),
                           violin=FALSE,...)
# generate forest plot and joint and marginal density plots.
{
  q975 <- qnorm(0.975)
  stopifnot(length(mulim)==2, length(taulim)<=2)

  if (all(is.finite(taulim))) {
    if ((length(taulim)==2) && (taulim[1]>=0) && (taulim[2]>taulim[1]))
      taurange <- taulim
    else if ((length(taulim)==1) && (taulim>0))
      taurange <- c(0, taulim)
    else
      taurange <- c(0, x$qposterior(tau=0.995)*1.1)
  } else {
    taurange <- c(0, x$qposterior(tau=0.995)*1.1)
  }

  if (all(is.finite(mulim)) && (mulim[1] < mulim[2])) {
    murange <- mulim
  } else {
    murange <- x$qposterior(mu=c(0.005,0.995))
    murange <- murange + c(-1,1)*diff(murange)*0.05
  }
  
  forestPlot <- function(x, main="", violin=violin, leftmargin=8)
  {
    original.margins <- par("mar")
    plotmargins <- original.margins
    plotmargins[2] <- leftmargin
    # determine x-axis range:
    xrange <- range(c(x$y-q975*x$sigma, x$y+q975*x$sigma,
                      x$summary[c("95% lower", "95% upper"),c("mu","theta")]))
    # empty plot:
    par("mar"=plotmargins)
    plot(xrange, c(-2, x$k)+c(-1,1)*0.5,
         type="n", axes=FALSE,
         xlab="", ylab="", main=main)
    mtext(side=1, line=par("mgp")[1], expression("effect"))
    # add horizontal line dividing data and estimates:
    abline(h=0, col="lightgrey")
    # add vertical "posterior median effect" line:
    abline(v=x$summary["median","mu"], col="black", lty="12")
    if (violin) {
      ticklength <- 0.45 # (here: maximum height of gaussian)
      maxdens <- c(dnorm(0, sd=x$sigma),
                   x$dposterior(mu=x$summary["mode","mu"]),
                   x$dposterior(theta=x$summary["mode","theta"], predict=TRUE))
      relmaxdens <- maxdens / max(maxdens)
      # density for central 95%
      narg1 <- seq(-q975, q975, le=51)
      ndens1 <- dnorm(narg1) / dnorm(0)
      # density for tails out to +/- 6 sigma
      narg2 <- seq(q975, 6.0, le=40)
      ndens2 <- dnorm(narg2) / dnorm(0)
      
      relsigma <- x$sigma / min(x$sigma)
      for (i in 1:x$k) { # loop over estimates
        # right tail:
        polygon(x$y[i] + c(narg2,rev(narg2))*x$sigma[i],
                x$k-(i-1) + c(ndens2,-rev(ndens2)) * relmaxdens[i] * ticklength,
                col="grey", border=NA)
        # left tail:
        polygon(x$y[i] - c(narg2,rev(narg2))*x$sigma[i],
                x$k-(i-1) + c(ndens2,-rev(ndens2)) * relmaxdens[i] * ticklength,
                col="grey", border=NA)
        # central chunk:
        polygon(x$y[i] + c(narg1,rev(narg1))*x$sigma[i],
                x$k-(i-1) + c(ndens1,-ndens1) * relmaxdens[i] * ticklength,
                col="black", border=NA)
        # central vertical line:
        lines(c(x$y[i], x$y[i]), x$k-(i-1)+ c(-1,1) * relmaxdens[i] * ticklength, col="grey")
      }
    }
    else {
      ticklength <- 0.3
      # draw horizontal lines for individual-study confidence intervals:
      matlines(rbind(x$y-q975*x$sigma, x$y+q975*x$sigma),
               rbind(x$k:1, x$k:1), col=grey(0.3), lty="solid")
      # draw vertical lines for individual-study effect estimates:
      matlines(rbind(x$y, x$y),
               rbind((x$k:1)+ticklength, (x$k:1)-ticklength), col="black", lty="solid")
    }
    if (violin) {
      # draw blob for mean effect estimate:
      ticklength <- 0.45 # (here: maximum height of gaussian)
      quant <- x$summary[c("95% lower","95% upper"),"mu"]
      quant <- c(quant[1]-2*(x$summary["median","mu"]-quant[1]),
                 quant,
                 quant[2]+2*(quant[2]-x$summary["median","mu"]))
      arg1 <- seq(quant[1], quant[2], le=50)  # left tail
      arg2 <- seq(quant[2], quant[3], le=51)  # central chunk
      arg3 <- seq(quant[3], quant[4], le=50)  # right tail
      dmode <- x$dposterior(mu=x$summary["mode","mu"])
      dens1 <- x$dposterior(mu=arg1) / dmode * relmaxdens[x$k+1]
      dens2 <- x$dposterior(mu=arg2) / dmode * relmaxdens[x$k+1]
      dens3 <- x$dposterior(mu=arg3) / dmode * relmaxdens[x$k+1]
      polygon(c(arg1, rev(arg1)), -1+c(dens1, -rev(dens1))*ticklength, col="grey", border=NA)
      polygon(c(arg3, rev(arg3)), -1+c(dens3, -rev(dens3))*ticklength, col="grey", border=NA)
      polygon(c(arg2, rev(arg2)), -1+c(dens2, -rev(dens2))*ticklength, col="black", border=NA)
      dm <- x$dposterior(mu=x$summary["median","mu"]) / dmode * relmaxdens[x$k+1]
      lines(rep(x$summary["median","mu"],2), -1+c(-1,1)*dm*ticklength, col="grey")
      
      # draw blob for prediction interval:
      quant <- x$summary[c("95% lower","95% upper"),"theta"]
      quant <- c(quant[1]-2*(x$summary["median","theta"]-quant[1]),
                 quant,
                 quant[2]+2*(quant[2]-x$summary["median","theta"]))
      arg1 <- seq(quant[1], quant[2], le=50)
      arg2 <- seq(quant[2], quant[3], le=51)
      arg3 <- seq(quant[3], quant[4], le=50)
      dmode <- x$dposterior(theta=x$summary["mode","theta"], predict=TRUE)
      dens1 <- x$dposterior(theta=arg1, predict=TRUE) / dmode * relmaxdens[x$k+2]
      dens2 <- x$dposterior(theta=arg2, predict=TRUE) / dmode * relmaxdens[x$k+2]
      dens3 <- x$dposterior(theta=arg3, predict=TRUE) / dmode * relmaxdens[x$k+2]
      polygon(c(arg1, rev(arg1)), -2+c(dens1, -rev(dens1))*ticklength, col="grey", border=NA)
      polygon(c(arg3, rev(arg3)), -2+c(dens3, -rev(dens3))*ticklength, col="grey", border=NA)
      polygon(c(arg2, rev(arg2)), -2+c(dens2, -rev(dens2))*ticklength, col="black", border=NA)
      dm <- x$dposterior(theta=x$summary["median","theta"], predict=TRUE) / dmode * relmaxdens[x$k+2]
      lines(rep(x$summary["median","theta"],2), -2+c(-1,1)*dm*ticklength, col="grey")
    }
    else {
      # draw diamond for mean effect estimate:
      ticklength <- 0.4
      polygon(x$summary[c("95% lower", "median", "95% upper", "median"),"mu"],
              rep(-1,4)+c(0,1,0,-1)*ticklength,
              border=NA, col=grey(0.4))
      # draw rectangle for effect prediction interval:
      ticklength <- 0.2
      polygon(x$summary[c("95% lower", "95% lower", "95% upper", "95% upper"),"theta"],
              rep(-2,4)+c(-1,1,1,-1)*ticklength,
              border=NA, col=grey(0.4))
    }
    # add axes & bounding box:
    axis(2, at=x$k:1, labels=x$labels, las=1)
    axis(2, at= c(-1,-2), labels=c(expression("effect "*mu), expression("prediction "*theta[italic(k)+1])), las=1)
    axis(1); box()
    par("mar"=original.margins)
    invisible()
  }
  
  jointdensity <- function(x, main="")
  {
    # range of tau values:
    tau <- seq(taurange[1], taurange[2],le=50)
    # range of mu values:
    mu <- seq(murange[1], murange[2], le=50)
    # grid of tau/mu value combinations:
    taumu <- expand.grid(tau,mu)
    # evaluate posterior density at grid points:
    post <- matrix(x$dposterior(tau=taumu[,1], mu=taumu[,2], log=TRUE),
                   nrow=length(tau), ncol=length(mu))
    # determine MAP value:
    map.value <- x$dposterior(x$MAP["joint",1], x$MAP["joint",2], log=TRUE)
    map.Inf <- !is.finite(map.value)
    if (map.Inf) {
      map.value <- max(post[is.finite(post)])
      warning("Non-finite posterior density, no contour lines drawn.", call.=FALSE)
    }
    # draw greyscale image:
    image(tau, mu, exp(post), axes=FALSE,
          col=grey((seq(1,0,le=128))^2),
          breaks=seq(0,exp(map.value),length=129),
          xlab="", ylab="", main=main, sub="(joint posterior density)", cex.sub=0.8)
    # add the blue lines for conditional mean & 95% confidence bounds:
    tau2 <- seq(taurange[1], taurange[2],le=200)
    cm <- x$cond.moment(tau=tau2)
    lines(tau2, cm[,"mean"], col="blue")  
    lines(tau2, cm[,"mean"]-q975*cm[,"sd"], col="blue", lty="dashed")  
    lines(tau2, cm[,"mean"]+q975*cm[,"sd"], col="blue", lty="dashed")
    # add green lines for marginal means & 95% confidence bounds:
    abline(v=x$summary[c("95% lower", "median", "95% upper"),"tau"],
           h=x$summary[c("95% lower", "median", "95% upper"),"mu"],
           col="green2", lty=c("16", "44", "16"))
    # draw ML estimate:
    points(x$ML["joint",1], x$ML["joint",2], col="magenta", pch=4)
    # draw MAP estimate:
    points(x$MAP["joint",1], x$MAP["joint",2], col="red", pch=3)
    if (!map.Inf) {
      # add contour lines:
      contour(tau, mu, post-map.value, add=TRUE, col="red",
              levels=-0.5*qchisq(p=c(0.5, 0.9, 0.95, 0.99), df=2),
              labels=paste(c(50, 90, 95, 99),"%",sep=""))
    }
    # add axes, bounding box, labels, ...
    axis(1); axis(2); box()
    mtext(side=1, line=par("mgp")[1], expression("heterogeneity "*tau))
    mtext(side=2, line=par("mgp")[1], expression("effect "*mu))
    invisible()
  }

  mumarginal <- function(x, main="", priorline=prior)
  {
    # range of mu values:
    mu <- seq(murange[1]-diff(murange)*0.05, murange[2]+diff(murange)*0.05, le=200)
    # corresponding posterior density:
    dens <- x$dposterior(mu=mu)
    # empty plot:
    plot(murange, c(0,max(dens,na.rm=TRUE)), type="n", axes=FALSE,
         xlab="", ylab="", main=main)
    # light grey shaded contour for density across whole range:
    polygon(c(min(mu), mu, max(mu)), c(0,dens,0), border=NA, col=grey(0.90))
    # dark grey shaded contour for density within 95% bounds:
    indi <- ((mu>=x$summary["95% lower","mu"]) & (mu<=x$summary["95% upper","mu"]))
    polygon(c(rep(x$summary["95% lower","mu"],2), mu[indi], rep(x$summary["95% upper","mu"],2)),
            c(0, x$dposterior(mu=x$summary["95% lower","mu"]),
              dens[indi], x$dposterior(mu=x$summary["95% upper","mu"]), 0),
            border=NA, col=grey(0.80))
    # vertical line for posterior median:
    lines(rep(x$summary["median","mu"],2),
          c(0,x$dposterior(mu=x$summary["median","mu"])), col=grey(0.6))
    # actual density line:
    lines(mu, dens, col="black")
    # x-axis:
    abline(h=0, col=grey(0.40))
    # prior density (if requested):
    if (priorline) lines(mu, x$dprior(mu=mu), col="black", lty="dashed")
    # add axes, labels, bounding box, ...
    mtext(side=1, line=par("mgp")[1], expression("effect "*mu))
    mtext(side=2, line=par("mgp")[2], expression("marginal posterior density"))
    axis(1); box()
    invisible()
  }

  taumarginal <- function(x, main="", priorline=prior)
  {
    # range of tau values:
    tau <- seq(max(c(0,taurange[1]-0.1*diff(taurange))),
                   taurange[2]+0.1*diff(taurange), le=200)
    # corresponding posterior density:
    dens <- x$dposterior(tau=tau)
    # empty plot:
    maxdens <- max(dens[is.finite(dens)],na.rm=TRUE)
    plot(c(taurange[1],taurange[2]), c(0,maxdens),         
         type="n", axes=FALSE, xlab="", ylab="", main=main)
    # "fix" diverging density:
    dens[!is.finite(dens)] <- 10*maxdens
    # light grey shaded contour for density across whole range:
    polygon(c(0,tau,max(tau)), c(0,dens,0), border=NA, col=grey(0.90))
    # dark grey shaded contour for density within 95% bounds:
    indi <- ((tau>=x$summary["95% lower","tau"]) & (tau<=x$summary["95% upper","tau"]))
    polygon(c(rep(x$summary["95% lower","tau"],2), tau[indi], rep(x$summary["95% upper","tau"],2)),
            c(0, min(c(x$dposterior(tau=x$summary["95% lower","tau"]), 10*maxdens)),
              dens[indi], x$dposterior(tau=x$summary["95% upper","tau"]), 0),
            border=NA, col=grey(0.80))
    # vertical line at posterior median:
    lines(rep(x$summary["median","tau"],2), c(0,x$dposterior(tau=x$summary["median","tau"])), col=grey(0.6))
    # actual density line:
    lines(tau, dens, col="black")
    # x-axis:
    abline(h=0, v=0, col=grey(0.40))
    # prior density (if requested):
    if (priorline) lines(tau, x$dprior(tau=tau), col="black", lty="dashed")
    # add axes, labels, bounding box, ...
    mtext(side=1, line=par("mgp")[1], expression("heterogeneity "*tau))
    mtext(side=2, line=par("mgp")[2], expression("marginal posterior density"))
    axis(1); box()
    invisible()
  }

  # main function:
  stopifnot(all(is.element(which, 1:4)),
            length(which)==length(unique(which)))  
  par.ask <- par("ask")
  for (i in 1:length(which)) {
    if (which[i]==1) forestPlot(x, main, violin=violin, leftmargin=forest.margin)
    else if (which[i]==2) jointdensity(x, main)
    else if (which[i]==3) mumarginal(x, main, priorline=prior)
    else if (which[i]==4) taumarginal(x, main, priorline=prior)
    if (i==1) {
      on.exit(par(ask=par.ask))
      par(ask=TRUE)
    }
  }
  par(ask=par.ask)
  invisible(x)
}


dlomax <- function(x, shape=1, scale=1, log=FALSE)
# probability density function for Lomax distribution
{
  stopifnot(length(shape)==1, length(scale)==1,
            all(shape>0), all(scale>0))
  result <- rep(-Inf, length(x))
  result[x>=0] <- log(shape)-log(scale)-(shape+1)*log(1+x[x>=0]/scale)
  if (!log) result <- exp(result)
  return(result)
}

plomax <- function(q, shape=1, scale=1)
# cumulative probability distribution function (CDF) of a Lomax distribution
{
  stopifnot(length(shape)==1, length(scale)==1,
            all(shape>0), all(scale>0))
  result <- rep(NA, length(q))
  result[q<=0] <- 0
  result[q>0] <- 1-(1+q[q>0]/scale)^(-shape)
  return(result)
}

qlomax <- function(p, shape=1, scale=1)
# quantile function (inverse CDF) of a Lomax distribution
{
  stopifnot(length(shape)==1, length(scale)==1,
            all(shape>0), all(scale>0))
  result <- rep(NA, length(p))
  result[p==1] <- Inf
  result[p==0] <- 0
  proper <- (is.finite(p) & (p>0) & (p<1))
  result[proper] <- ((1-p[proper])^(-1/shape) - 1) * scale
  return(result)
}

rlomax <- function(n, shape=1, scale=1)
# random number generation for a Lomax distribution
# (based on inversion method)
{
  stopifnot(length(shape)==1, length(scale)==1,
            all(shape>0), all(scale>0))
  u <- runif(n)
  result <- qlomax(u, scale=scale, shape=shape)
  return(result)
}

elomax <- function(shape=1, scale=1)
# expectation of a Lomax distribution
{
  stopifnot(length(shape)==1, length(scale)==1,
            all(shape>0), all(scale>0))
  if (shape > 1)
    expectation <- scale / (shape-1)
  else
    expectation <- Inf
  return(expectation)
}

vlomax <- function(shape=1, scale=1)
# variance of a Lomax distribution
{
  stopifnot(length(shape)==1, length(scale)==1,
            all(shape>0), all(scale>0))
  if (shape > 2)
    variance <- elomax(shape,scale)^2 * (shape/(shape-2))
  else
    variance <- Inf
  return(variance)
}


dhalfnormal <- function(x, scale=1, log=FALSE)
# probability density function of a half-normal distribution
{
  stopifnot(scale>0)
  result <- rep(-Inf, length(x))
  result[x>=0] <- log(2) + dnorm(x[x>=0], mean=0, sd=scale, log=TRUE)
  if (!log) result <- exp(result)
  return(result)  
}

phalfnormal <- function(q, scale=1)
# cumulative distribution function (CDF) of a half-normal distribution
{
  stopifnot(scale>0)
  result <- rep(0, length(q))
  result[q>0] <- 2.0 * (pnorm(q[q>0], mean=0, sd=scale) - 0.5)
  return(result)  
}

qhalfnormal <- function(p, scale=1)
# quantile function (inverse CDF) of a half-normal distribution
{
  stopifnot(scale>0)
  result <- rep(NA, length(p))
  proper <- (is.finite(p) & (p>=0) & (p<=1))
  result[proper] <- qnorm(p[proper]/2.0 + 0.5, mean=0, sd=scale)
  return(result)  
}

rhalfnormal <- function(n, scale=1)
# random number generation for half-normal distribution
{
  stopifnot(scale>0)
  return(abs(rnorm(n=n, mean=0, sd=scale)))  
}

ehalfnormal <- function(scale=1)
# expectation of a half-normal distribution
{
  stopifnot(scale>0)
  return(scale*sqrt(2/pi))  
}

vhalfnormal <- function(scale=1)
# variance of a half-normal distribution
{
  stopifnot(scale>0)
  return(scale^2*(1-2/pi))  
}


dhalft <- function(x, df, scale=1, log=FALSE)
# probability density function for half-Student-t distribution
{
  stopifnot(scale>0, df>0)
  result <- rep(-Inf, length(x))
  result[x>=0] <- log(2) - log(scale) + dt(x[x>=0]/scale, df=df, log=TRUE)
  if (!log) result <- exp(result)
  return(result)  
}

phalft <- function(q, df, scale=1)
# cumulative distribution function (CDF) of a half-Student-t distribution
{
  stopifnot(scale>0, df>0)
  result <- rep(0, length(q))
  result[q>0] <- 2.0 * (pt(q[q>0]/scale, df=df) - 0.5)
  return(result)  
}

qhalft <- function(p, df, scale=1)
# quantile function (inverse CDF) of a half-Student-t distribution
{
  stopifnot(scale>0, df>0)
  result <- rep(NA, length(p))
  proper <- (is.finite(p) & (p>=0) & (p<=1))
  result[proper] <- qt(p[proper]/2.0 + 0.5, df=df) * scale
  return(result)  
}

rhalft <- function(n, df, scale=1)
# random number generation for half-Student-t distribution
{
  stopifnot(scale>0, df>0)
  return(abs(rt(n=n, df=df))*scale)  
}

ehalft <- function(df, scale=1)
# expectation of a half-Student-t distribution
{
  stopifnot(scale>0, df>0)
  if (df>1)
    expectation <- exp(log(2)+log(scale)+0.5*log(df)-0.5*log(pi)+lgamma((df+1)/2)-lgamma(df/2)-log(df-1))
  else
    expectation <- Inf
  return(expectation)
}

vhalft <- function(df, scale=1)
# variance of a half-Student-t distribution
{
  stopifnot(scale>0, df>0)
  if (df>2)
    variance <- scale^2 * ((df / (df-2)) - ehalft(df=df, scale=1.0)^2)
  else
    variance <- Inf
  return(variance)
}


dhalfcauchy <- function(x, scale=1, log=FALSE)
# probability density function of a half-Cauchy distribution
{
  stopifnot(scale>0)
  result <- rep(-Inf, length(x))
  result[x>=0] <- log(2) + dcauchy(x[x>=0], location=0, scale=scale, log=TRUE)
  if (!log) result <- exp(result)
  return(result)  
}

phalfcauchy <- function(q, scale=1)
# cumulative distribution function (CDF) of a half-Cauchy distribution
{
  stopifnot(scale>0)
  result <- rep(0, length(q))
  result[q>0] <- 2.0 * (pcauchy(q[q>0], location=0, scale=scale) - 0.5)
  return(result)  
}

qhalfcauchy <- function(p, scale=1)
# quantile function (inverse CDF) of a half-Cauchy distribution
{
  stopifnot(scale>0)
  result <- rep(NA, length(p))
  proper <- (is.finite(p) & (p>=0) & (p<=1))
  result[proper] <- qcauchy(p[proper]/2.0 + 0.5, location=0, scale=scale)
  return(result)  
}

rhalfcauchy <- function(n, scale=1)
# random number generation for half-Cauchy distribution
{
  stopifnot(scale>0)
  return(abs(rcauchy(n=n, location=0, scale=scale)))
}

ehalfcauchy <- function(scale=1)
# expectation of a half-Cauchy distribution
{
  stopifnot(scale>0)
  return(Inf)
}

vhalfcauchy <- function(scale=1)
# variance of a half-Cauchy distribution
{
  stopifnot(scale>0)
  return(Inf)
}


dhalflogistic <- function(x, scale=1, log=FALSE)
# probability density function of a half-logistic distribution
{
  stopifnot(scale>0)
  result <- rep(-Inf, length(x))
  result[x>=0] <- log(2) + dlogis(x[x>=0], location=0, scale=scale, log=TRUE)
  if (!log) result <- exp(result)
  return(result)  
}

phalflogistic <- function(q, scale=1)
# cumulative distribution function (CDF) of a half-logistic distribution
{
  stopifnot(scale>0)
  result <- rep(0, length(q))
  result[q>0] <- 2.0 * (plogis(q[q>0], location=0, scale=scale) - 0.5)
  return(result)  
}

qhalflogistic <- function(p, scale=1)
# quantile function (inverse CDF) of a half-logistic distribution
{
  stopifnot(scale>0)
  result <- rep(NA, length(p))
  proper <- (is.finite(p) & (p>=0) & (p<=1))
  result[proper] <- qlogis(p[proper]/2.0 + 0.5, location=0, scale=scale)
  return(result)  
}

rhalflogistic <- function(n, scale=1)
# random number generation for half-logistic distribution
{
  stopifnot(scale>0)
  return(abs(rlogis(n=n, location=0, scale=scale)))
}

ehalflogistic <- function(scale=1)
# expectation of a half-logistic distribution
{
  stopifnot(scale>0)
  return(scale * log(4))
}

vhalflogistic <- function(scale=1)
# variance of a half-logistic distribution
{
  stopifnot(scale>0)
  return(scale^2 * (((pi^2)/3) - log(4)^2))
}


drayleigh <- function(x, scale=1, log=FALSE)
# probability density function of the Rayleigh distribution
{
  stopifnot(scale>0)
  result <- rep(-Inf, length(x))
  result[x>0] <- log(x[x>0])-2*log(scale)-0.5*(x[x>0]/scale)^2
  if (!log) result <- exp(result)
  return(result)
}

prayleigh <- function(q, scale=1)
# cumulative distribution function (CDF) of the Rayleigh distribution
{
  stopifnot(scale>0)
  result <- rep(0,length(q))
  result[q>0] <- 1-exp(-0.5*(q[q>0]/scale)^2)
  return(result)
}

qrayleigh <- function(p, scale=1)
# quantile function (inverse CDF) of the Rayleigh distribution
{
  stopifnot(scale>0)
  proper <- (is.finite(p) & (p>=0) & (p<=1))
  result <- rep(NA, length(p))
  result[proper] <- scale * sqrt(-2*log(1-p[proper]))
  return(result)
}

rrayleigh <- function(n, scale=1)
# random number generation for the Rayleigh distribution
{
  stopifnot(scale>0)
  return(sqrt(rexp(n, rate=1/(2*scale^2))))
}

erayleigh <- function(scale=1)
# expectation of a Rayleigh distribution
{
  stopifnot(scale>0)
  return(scale * sqrt(pi/2))
}

vrayleigh <- function(scale=1)
# variance of a Rayleigh distribution
{
  stopifnot(scale>0)
  return(scale^2 * (4-pi)/2)
}


dinvchi <- function(x, df, scale=1, log=FALSE)
# probability density function of a scaled inverse chi distribution
{
  stopifnot(length(df)==1, is.finite(df), df>0,
            length(scale)==1, is.finite(scale), scale>0)
  idx <- (is.finite(x) & (x>0))
  ldens <- rep(-Inf, length(x))
  ldens[idx] <- -(df/2-1)*log(2)-lgamma(df/2)-log(scale) -(df+1)*(log(x[idx])-log(scale)) - 0.5*(scale/x[idx])^2
  if (log) result <- ldens
  else     result <- exp(ldens)
  return(result)
}

pinvchi <- function(q, df, scale=1.0, lower.tail=TRUE, log.p=FALSE)
# cumulative distribution function (CDF) of a scaled inverse chi distribution
{
  stopifnot(length(df)==1, is.finite(df), df>0,
            length(scale)==1, is.finite(scale), scale>0)
  return(stats::pchisq(q=(q/scale)^(-2), df=df, lower.tail=(!lower.tail), log.p=log.p))
}

qinvchi <- function(p, df, scale=1.0, lower.tail=TRUE, log.p=FALSE)
# quantile function (inverse CDF) of a scaled inverse chi distribution
{
  stopifnot(length(df)==1, is.finite(df), df>0,
            length(scale)==1, is.finite(scale), scale>0)
  return(scale*stats::qchisq(p=p, df=df, lower.tail=(!lower.tail), log.p=log.p)^(-0.5))
}

rinvchi <- function(n=1, df, scale=1.0)
# random number generation for a scaled inverse chi distribution
{
  stopifnot(length(df)==1, is.finite(df), df>0,
            length(scale)==1, is.finite(scale), scale>0)
  return(scale * stats::rchisq(n=n, df=df)^(-0.5))
}

einvchi <- function(df, scale=1)
# expectation of a scaled inverse chi distribution
{
  stopifnot(length(scale)==1, is.finite(scale), scale>0,
            length(df)==1, is.finite(df), df>0)
  if (df<=1) result <- Inf
  else       result <- exp(log(scale)+lgamma((df-1)/2)-lgamma(df/2)-0.5*log(2))
  return(result)
}

vinvchi <- function(df, scale=1)
# variance of a scaled inverse chi distribution
{
  stopifnot(length(scale)==1, is.finite(scale), scale>0,
            length(df)==1, is.finite(df), df>0)
  if (df<=2) result <- Inf
  else       result <- scale^2/(df-2) - einvchi(df=df, scale=scale)^2
  return(result)
}


# Turner & al. prior data:
# ========================
# list of 5 possible intervention comparison types:
ctypes <- c("pharmacological vs. placebo / control",
            "pharmacological vs. pharmacological",
            "non-pharmacological vs. placebo / control",
            "non-pharmacological vs. pharmacological",
            "non-pharmacological vs. non-pharmacological")
# list of 16 possible outcome types:
otypes <- c("all-cause mortality",
            "obstetric outcomes",
            "cause-specific mortality / major morbidity event / composite (mortality or morbidity)",
            "resource use / hospital stay / process",
            "surgical / device related success / failure",
            "withdrawals / drop-outs",
            "internal / structure-related outcomes",
            "general physical health indicators",
            "adverse events",
            "infection / onset of new disease",
            "signs / symptoms reflecting continuation / end of condition",
            "pain",
            "quality of life / functioning (dichotomized)",
            "mental health indicators",
            "biological markers (dichotomized)",
            "subjective outcomes (various)")  
# matrix of mu values (see Tab.IV):
meanmat <- matrix(-c(3.95, 3.52, 3.71, 2.34, 2.14, 2.99, 2.71, 2.29, 1.87, 2.49, 2.06, 1.83, 2.54, 2.12, 1.77, 2.70,
                     4.18, 3.75, 3.95, 2.58, 2.37, 3.23, 2.94, 2.53, 2.10, 2.73, 2.29, 2.06, 2.78, 2.35, 2.00, 2.93,
                     4.17, 3.74, 3.93, 2.56, 2.36, 3.21, 2.93, 2.51, 2.10, 2.71, 2.28, 2.05, 2.77, 2.34, 1.99, 2.92,
                     2.92, 2.49, 2.68, 1.31, 1.11, 1.96, 1.67, 1.26, 0.84, 1.46, 1.03, 0.80, 1.51, 1.09, 0.74, 1.67,
                     3.50, 3.08, 3.27, 1.90, 1.69, 2.55, 2.26, 1.85, 1.43, 2.05, 1.61, 1.38, 2.10, 1.67, 1.33, 2.26),
                  nrow=16, ncol=5,
                 dimnames=list(otypes, ctypes))
# matrix of sigma values (see Tab.IV):
sdmat <- matrix(c(1.34, 1.74, 1.74, 1.74, 1.74, 1.74, 1.74, 1.53, 1.52, 1.52, 1.51, 1.52, 1.54, 1.53, 1.52, 1.52,
                  1.41, 1.79, 1.79, 1.79, 1.79, 1.79, 1.79, 1.58, 1.58, 1.58, 1.58, 1.58, 1.60, 1.60, 1.58, 1.58,
                  1.55, 1.91, 1.91, 1.91, 1.91, 1.91, 1.92, 1.72, 1.71, 1.71, 1.71, 1.71, 1.73, 1.72, 1.71, 1.71,
                  1.02, 1.50, 1.51, 1.50, 1.50, 1.51, 1.51, 1.25, 1.24, 1.24, 1.24, 1.25, 1.27, 1.27, 1.24, 1.25,
                  1.26, 1.68, 1.68, 1.68, 1.68, 1.68, 1.68, 1.46, 1.45, 1.45, 1.45, 1.45, 1.47, 1.47, 1.45, 1.45),
                nrow=16, ncol=5,
                dimnames=list(otypes, ctypes))
# array of mu/sigma values:
TurnerEtAlParameters <- array(c(as.vector(meanmat), as.vector(sdmat)),
                              dim=c(16,5,2),
                              dimnames=list("outcome"=otypes, "comparison"=ctypes, "parameter"=c("mu","sigma")))
# remove obsolete objects:
rm(list=c("ctypes", "otypes", "meanmat", "sdmat"))


TurnerEtAlPrior <- function(outcome=c(NA,
                                      "all-cause mortality",
                                      "obstetric outcomes",
                                      "cause-specific mortality / major morbidity event / composite (mortality or morbidity)",
                                      "resource use / hospital stay / process",
                                      "surgical / device related success / failure",
                                      "withdrawals / drop-outs",
                                      "internal / structure-related outcomes",
                                      "general physical health indicators",
                                      "adverse events",
                                      "infection / onset of new disease",
                                      "signs / symptoms reflecting continuation / end of condition",
                                      "pain",
                                      "quality of life / functioning (dichotomized)",
                                      "mental health indicators",
                                      "biological markers (dichotomized)",
                                      "subjective outcomes (various)"),
                            comparator1=c("pharmacological", "non-pharmacological", "placebo / control"),
                            comparator2=c("pharmacological", "non-pharmacological", "placebo / control"))
#
# Function to return prior parameters, densities, etc. as proposed in
#
#   Turner et al.  Predictive distributions for between-study heterogeneity
#   and simple methods for their application in Bayesian meta-analysis.
#   Statisics in Medicine 4(6):984-998, 2015.
#
# (see Table IV).
#
{
  param <- matrix(NA, nrow=2, ncol=2,
                  dimnames=list(c("tau^2","tau"),
                                "log-normal"=c("mu","sigma")))
  if (all(!is.na(outcome))) {
    # match the provided arguments:
    outcome     <- match.arg(outcome)
    comparator1 <- match.arg(comparator1)
    comparator2 <- match.arg(comparator2)
    # list of 3 possible comparators:
    clist <- c("pharmacological", "non-pharmacological", "placebo / control")
    # index matrix to match pairs of comparators to one of five possible scenarios:
    cmatrix <- matrix(c(2,4,1,4,5,3,1,3,NA), nrow=3, ncol=3,
                      dimnames=list(clist, clist))
    # list of 5 possible intervention comparison types:
    ctypes <- dimnames(TurnerEtAlParameters)[["comparison"]]
    # figure out current comparison scenario:
    comparisontype <- ctypes[cmatrix[comparator1, comparator2]]
    # assemble function output:
    param["tau^2","mu"]    <- TurnerEtAlParameters[outcome, comparisontype, "mu"]
    param["tau^2","sigma"] <- TurnerEtAlParameters[outcome, comparisontype, "sigma"]
  } else {  # the "marginal" setting, see p.993
    outcome <- comparisontype <- "any"
    param <- matrix(NA, nrow=2, ncol=2,
                    dimnames=list(c("tau^2","tau"),
                                  "log-normal"=c("mu","sigma")))
    param["tau^2","mu"]    <- -2.56
    param["tau^2","sigma"] <- 1.74
  }
  param["tau","mu"]      <- param["tau^2","mu"] / 2
  param["tau","sigma"]   <- param["tau^2","sigma"] / 2
  result <- list("parameters"      = param,
                 "outcome.type"    = outcome,
                 "comparison.type" = comparisontype,
                 "dprior"          = function(tau, log=FALSE) {return(dlnorm(tau, meanlog=param["tau","mu"],
                                                                                  sdlog=param["tau","sigma"], log=log))},
                 "pprior"          = function(tau) {return(plnorm(tau, meanlog=param["tau","mu"], sdlog=param["tau","sigma"]))},
                 "qprior"          = function(p)   {return(qlnorm(p,   meanlog=param["tau","mu"], sdlog=param["tau","sigma"]))})
  attr(result$dprior, "bayesmeta.label") <- paste("log-normal(mu=",sprintf("%1.3f",param["tau","mu"]),
                                                  ", sigma=",sprintf("%1.3f",param["tau","sigma"]),")",sep="")
  return(result)
}



# Rhodes & al. prior data:
# ========================
# list of 3 possible intervention comparison types:
ctypes <- c("pharmacological vs. placebo / control",
            "pharmacological vs. pharmacological",
            "non-pharmacological (any)")

# list of 5 possible outcome types:
otypes <- c("obstetric outcome",
            "resource use and hospital stay / process",
            "internal and external structure-related outcome",            
            "general physical health and adverse event and pain and quality of life / functioning",
            "signs / symptoms reflecting continuation / end of condition and infection / onset of new acute / chronic disease",
            "mental health outcome",
            "biological marker",
            "various subjectively measured outcomes")

# matrix of location parameter values (see Tab.3):
locationmat <- matrix(-c(4.13, 2.55, 2.43, 3.16, 3.00, 2.99, 3.41, 2.76,
                         4.40, 2.83, 2.70, 3.44, 3.27, 3.27, 3.68, 3.03,
                         3.99, 2.41, 2.29, 3.02, 2.86, 3.85, 3.27, 2.62),
                      nrow=8, ncol=3,
                      dimnames=list(otypes, ctypes))

# matrix of location parameter values for RESPIRATORY DISEASES (see Tab.A.3.1):
locationmat.r <- matrix(-c(6.03, 4.46, 4.33, 5.07, 4.90, 4.90, 5.31, 4.66,
                           6.31, 4.73, 4.61, 5.34, 5.18, 5.17, 5.59, 4.94,
                           5.89, 4.32, 4.19, 4.93, 4.76, 4.76, 5.17, 4.52),
                        nrow=8, ncol=3,
                        dimnames=list(otypes, ctypes))

# matrix of location parameter values for CANCER (see Tab.A.3.2):
locationmat.c <- matrix(-c(1.57,-0.01, 0.13, 0.60, 0.44, 0.43, 0.85, 0.20,
                           1.85, 0.27, 0.14, 0.88, 0.71, 0.71, 1.13, 0.48,
                           1.43,-0.15,-0.27, 0.46, 0.30, 0.29, 0.71, 0.06),
                        nrow=8, ncol=3,
                        dimnames=list(otypes, ctypes))

# matrix of scale parameter values (see Tab.3):
scalemat <- matrix(c(2.34, 2.73, 2.50, 2.50, 2.50, 2.16, 2.83, 2.58,
                     2.31, 2.70, 2.46, 2.44, 2.47, 2.14, 2.78, 2.59,
                     2.11, 2.57, 2.32, 2.27, 2.33, 1.93, 2.66, 2.41),
                   nrow=8, ncol=3,
                   dimnames=list(otypes, ctypes))

# matrix of scale parameter values for RESPIRATORY DISEASES (see Tab.A.3.1):
scalemat.r <- matrix(c(2.36, 2.74, 2.51, 2.51, 2.50, 2.17, 2.83, 2.59,
                       2.31, 2.70, 2.46, 2.45, 2.47, 2.14, 2.78, 2.59,
                       2.21, 2.57, 2.33, 2.28, 2.33, 1.94, 2.66, 2.41),
                     nrow=8, ncol=3,
                     dimnames=list(otypes, ctypes))

# matrix of scale parameter values for CANCER (see Tab.A.3.2):
scalemat.c <- matrix(c(2.45, 2.83, 2.61, 2.61, 2.60, 2.28, 2.93, 2.68,
                       2.41, 2.79, 2.56, 2.55, 2.57, 2.25, 2.87, 2.68,
                       2.24, 2.68, 2.45, 2.40, 2.46, 2.08, 2.78, 2.53),
                     nrow=8, ncol=3,
                     dimnames=list(otypes, ctypes))

# array of mu/sigma values:
RhodesEtAlParameters <- array(c(as.vector(locationmat), as.vector(scalemat),
                                as.vector(locationmat.r), as.vector(scalemat.r),
                                as.vector(locationmat.c), as.vector(scalemat.c)),
                              dim=c(8,3,2,3),
                              dimnames=list("outcome"=otypes, "comparison"=ctypes,
                                            "parameter"=c("location","scale"),
                                            "medical area"=c("other","respiratory","cancer")))
# remove obsolete objects:
rm(list=c("ctypes", "otypes", "locationmat", "scalemat", "locationmat.r", "scalemat.r", "locationmat.c", "scalemat.c"))

RhodesEtAlPrior <- function(outcome=c(NA,
                                      "obstetric outcome",
                                      "resource use and hospital stay / process",
                                      "internal and external structure-related outcome",            
                                      "general physical health and adverse event and pain and quality of life / functioning",
                                      paste("signs / symptoms reflecting continuation / end of condition and infection",
                                            "/ onset of new acute / chronic disease"),
                                      "mental health outcome",
                                      "biological marker",
                                      "various subjectively measured outcomes"),
                            comparator1=c("pharmacological", "non-pharmacological", "placebo / control"),
                            comparator2=c("pharmacological", "non-pharmacological", "placebo / control"),
                            area=c("other","respiratory","cancer"))
#
# Function to return prior parameters as proposed in
#
#   Rhodes et al.
#   Predictive distributions were developed for the extent of heterogeneity in meta-analyses of continuous outcome data.
#   Journal of Clinical Epidemiology 68(1):52-60, 2015.
#
# (see Table 3).
#
# Technically, this function mostly retrieves the parameters from a pre-defined array,
# the 3-dimensional "RhodesEtAlParameters" array. (This is more convenient than having
# to deal with the array itself, mostly due to partial argument matching, etc.)
#
{
  # match the provided arguments:
  if (all(!is.na(outcome))) { # settings from Tables 3, A.3.1 or A.3.2
    outcome     <- match.arg(outcome)
    comparator1 <- match.arg(comparator1)
    comparator2 <- match.arg(comparator2)
    area        <- match.arg(area)
    # list of 3 possible comparators:
    clist <- c("pharmacological", "non-pharmacological", "placebo / control")
    # index matrix to match pairs of comparators to one of 3 possible scenarios:
    cmatrix <- matrix(c(2,3,1,3,3,3,1,3,NA), nrow=3, ncol=3,
                      dimnames=list(clist, clist))
    # list of 3 possible intervention comparison types:
    ctypes <- dimnames(RhodesEtAlParameters)[["comparison"]]
    # figure out current comparison scenario:
    comparisontype <- ctypes[cmatrix[comparator1, comparator2]]
    # assemble function output:
    location <- RhodesEtAlParameters[outcome, comparisontype, "location", area]
    scale    <- RhodesEtAlParameters[outcome, comparisontype, "scale", area]
  } else { # the general (marginal) setting; see beginning of Sec. 3.3
    outcome <- comparisontype <- area <- "any"
    location <- -3.44
    scale <- 2.59
  }
  param <- matrix(NA, nrow=2, ncol=2,
                  dimnames=list(c("tau^2","tau"),
                                "log-t(5)"=c("location","scale")))
  param["tau^2","location"] <- location
  param["tau^2","scale"]    <- scale
  param["tau","location"]   <- location / 2
  param["tau","scale"]      <- scale / 2
  rm(list=c("location","scale"))
  dlt5 <- function(x, location=0, scale=1, log=FALSE)
  # density of log-t distribution with 5 d.f.
  {
    stopifnot(is.finite(location), is.finite(scale), scale>0)
    logdensity <- rep(-Inf, length(x))
    proper <- (is.finite(x) & (x>0))
    logdensity[proper] <- dt((log(x[proper])-location)/scale, df=5, log=TRUE) - log(scale) - log(x[proper])
    if (log) return(logdensity)
    else     return(exp(logdensity))
  }
  plt5 <- function(x, location=0, scale=1)
  # cumulative distribution function (CDF) of log-t distribution with 5 d.f.
  {
    stopifnot(is.finite(location), is.finite(scale), scale>0)
    return(pt((log(x)-location)/scale, df=5))
  }
  qlt5 <- function(p, location=0, scale=1)
  # quantile function (inverse CDF) of log-t distribution with 5 d.f.
  {
    stopifnot(is.finite(location), is.finite(scale), scale>0)
    return(exp((qt(p, df=5)*scale)+location))
  }
  result <- list("parameters"      = param,
                 "outcome.type"    = outcome,
                 "comparison.type" = comparisontype,
                 "medical.area"    = area,
                 "dprior"          = function(tau, log=FALSE) {return(dlt5(tau, location=param["tau","location"],
                                                                                scale=param["tau","scale"], log=log))},
                 "pprior"          = function(tau) {return(plt5(tau, location=param["tau","location"], scale=param["tau","scale"]))},
                 "qprior"          = function(p)   {return(qlt5(p, location=param["tau","location"], scale=param["tau","scale"]))})
  attr(result$dprior, "bayesmeta.label") <- paste("log-Student-t(location=",sprintf("%1.3f",param["tau","location"]),
                                                  ", scale=",sprintf("%1.3f",param["tau","scale"]),", d.f.=5)",sep="")
  return(result)
}


forest.bayesmeta <- function(x, xlab="effect size", refline=0, cex=1,...)
# forest plot for a "bayesmeta" object
# based on the "metafor" package's plotting functions
{
  if (!requireNamespace("metafor", quietly=TRUE))
    stop("required 'metafor' package not available!")
  metafor::forest.default(x=x$y, sei=x$sigma,
                          showweight=FALSE,  # (IV-weights don't make sense here)
                          ylim=c(-x$k-4, 1),
                          level=95,          # (95% level is intentionally hard-coded)
                          refline=refline,
                          xlab=xlab,
                          slab=x$labels,
                          rows=seq(-2, -x$k - 1, by = -1),
                          cex=cex, ...)
  metafor::addpoly(x$summary["median","mu"], ci.lb=x$summary["95% lower","mu"], ci.ub=x$summary["95% upper","mu"],
                   rows = -x$k-2.5, mlab=expression("mean effect ("*mu*")"), level=95, cex=cex, ...)
  metafor::addpoly(x$summary["median","theta"], ci.lb=x$summary["95% lower","theta"], ci.ub=x$summary["95% upper","theta"],
                   rows = -x$k-3.5, mlab=expression("prediction ("*vartheta[k+1]*")"), level=95, cex=cex, ...)
  
  plotdata <- cbind("95% lower"=x$y-qnorm(0.975)*x$sigma, "estimate"=x$y, "95% upper"=x$y+qnorm(0.975)*x$sigma)
  plotdata <- rbind(plotdata, t(x$summary[c("95% lower","median","95% upper"),c("mu","theta")]))
  rownames(plotdata) <- c(x$labels, c("mean effect(mu)", "prediction (theta)"))
  invisible(plotdata)
}


forestplot.bayesmeta <- function(x, labeltext,
                                 exponentiate  = FALSE,
                                 prediction    = TRUE,
                                 shrinkage     = TRUE,
                                 heterogeneity = TRUE,
                                 digits        = 2,
                                 plot          = TRUE,
                                 fn.ci_norm, fn.ci_sum, col, legend=NULL, boxsize, ...)
#
# ARGUMENTS:
#   x            :  a "bayesmeta" object.
#   labeltext    :  you may provide an alternative "labeltext" argument here
#                   (see also the "forestplot()" help).
#   exponentiate :  flag indicating whether to exponentiate numbers (figure and table).
#   prediction   :  flag indicating whether to show prediction interval.
#   shrinkage    :  flag indicating whether to show shrinkage estimates.
#   digits       :  number of significant digits to be shown (based on standard deviations).
#   plot         :  flag you can use to suppress actual plotting.
#   ...          :  further arguments passed to the "forestplot" function.
#   
# VALUE:
#   a list with components
#     $data       :  the meta-analyzed data, and mean and prediction estimates
#     $shrinkage  :  the shrinkage estimates for each study
#     $labeltext  :  the "forestplot()" function's "labeltext" argument used internally
#     $forestplot :  the "forestplot()" function's returned value
#
{
  if (!requireNamespace("forestplot", quietly=TRUE))
    stop("required 'forestplot' package not available!")
  if (utils::packageVersion("forestplot") < "1.5.2")
    warning("you may need to update 'forestplot' to a more recent version (>=1.5.2).")
  # auxiliary function:
  decplaces <- function(x, signifdigits=3)
  # number of decimal places (after decimal point)
  # to be displayed if you want at least "signifdigits" digits
  {
    return(max(c(0, -(floor(log10(x))-(signifdigits-1)))))
  }
  # some sanity checks for the provided arguments:
  stopifnot(is.element("bayesmeta", class(x)),
            length(digits)==1, digits==round(digits), digits>=0,
            length(exponentiate)==1, is.logical(exponentiate),
            length(prediction)==1, is.logical(prediction),
            length(shrinkage)==1, is.logical(shrinkage),
            length(plot)==1, is.logical(plot))
  # plotting data (1) -- the quoted estimates:
  q95 <- qnorm(0.975)
  ma.dat <- rbind(NA,
                  cbind(x$y, x$y - q95*x$sigma, x$y + q95*x$sigma),
                  x$summary[c("median", "95% lower", "95% upper"),"mu"],
                  x$summary[c("median", "95% lower", "95% upper"),"theta"])
  colnames(ma.dat) <- c("estimate", "lower", "upper")
  rownames(ma.dat) <- c("", x$label, "mean", "prediction")
  if (! prediction) ma.dat <- ma.dat[-(x$k+3),]
  # plotting data (2) -- the shrinkage estimates:
  ma.shrink <- rbind(NA,
                     t(x$theta)[,c("median","95% lower","95% upper")],
                     NA,NA)
  colnames(ma.shrink) <- c("estimate", "lower", "upper")
  rownames(ma.shrink) <- c("", x$label, "mean", "prediction")
  if (! prediction) ma.shrink <- ma.shrink[-(x$k+3),]
  if (exponentiate) {
    ma.dat    <- exp(ma.dat)
    ma.shrink <- exp(ma.shrink)
  }
  # generate "labeltext" data table for plot (unless already provided):
  if (missing(labeltext)) {
    # determine numbers of digits based on standard deviations:
    if (exponentiate) {
      stdevs <- c(exp(x$y)*x$sigma,
                  exp(x$summary["median","mu"])*x$summary["sd","mu"])
      if (prediction) stdevs <- c(stdevs, exp(x$summary["median","theta"])*x$summary["sd","theta"])
    } else {
      stdevs <- c(x$sigma, x$summary["sd","mu"])
      if (prediction) stdevs <- c(stdevs, x$summary["sd","theta"])
    }
    stdevs <- abs(stdevs[is.finite(stdevs) & (stdevs != 0)])
    formatstring <- paste0("%.", decplaces(stdevs, digits), "f")
    # fill data table:
    labeltext <- matrix(NA_character_, nrow=nrow(ma.dat), ncol=3)
    labeltext[1,] <- c("study","estimate", "95% CI")
    labeltext[,1] <- c("study", x$labels, "mean", "prediction")[1:nrow(ma.dat)]
    for (i in 2:(nrow(ma.dat))) {
      labeltext[i,2] <- sprintf(formatstring, ma.dat[i,"estimate"])
      labeltext[i,3] <- paste0("[", sprintf(formatstring, ma.dat[i,"lower"]),
                                 ", ", sprintf(formatstring, ma.dat[i,"upper"]), "]")
    }
  }
  # add horizontal lines to plot:
  horizl <- list(grid::gpar(col="grey"), grid::gpar(col="grey"))
  names(horizl) <- as.character(c(2,x$k+2))
  # specify function(s) for drawing estimates / shrinkage estimates:
  if (missing(fn.ci_norm)) {
    if (shrinkage) {
      fn.ci_norm <- list(function(...) {forestplot::fpDrawPointCI(pch=15,...)},
                         function(...) {forestplot::fpDrawPointCI(pch=18,...)})
    } else {
      fn.ci_norm <- function(...) {forestplot::fpDrawPointCI(pch=15,...)}
    }
  }
  # specify function(s) for drawing summaries (diamond / bar):
  if (missing(fn.ci_sum)) {
    fn.ci_sum <- list(NULL)
    for (i in 1:(x$k+2))
      fn.ci_sum[[i]] <- function(y.offset,...) {forestplot::fpDrawSummaryCI(y.offset=0.5,...)}
    if (prediction)
      fn.ci_sum[[x$k+3]] <- function(y.offset,...) {forestplot::fpDrawBarCI(y.offset=0.5,...)}
  }
  # specify colors:
  if (missing(col)) {
    if (shrinkage) {
      col <- forestplot::fpColors(box=c("black", "grey45"),
                                  lines=c("black","grey45"),
                                  summary="grey30")
    } else {
      col <- forestplot::fpColors(box="black", lines="black", summary="grey30")
    }
  }
  # specify plotting sizes:
  if (missing(boxsize)) {
    boxsize <- c(rep(0.25,x$k+1), 0.4)
    if (prediction) boxsize <- c(boxsize, 0.2)
  }
  if (shrinkage && is.null(legend))
    legend  <- c("quoted estimate", "shrinkage estimate")
  # specify data for plotting:
  if (shrinkage) { # (show shrinkage intervals)
    mean.arg  <- cbind(ma.dat[,1], ma.shrink[,1])
    lower.arg <- cbind(ma.dat[,2], ma.shrink[,2])
    upper.arg <- cbind(ma.dat[,3], ma.shrink[,3])
  } else {         # (no shrinkage intervals)
    mean.arg  <- ma.dat[,1]
    lower.arg <- ma.dat[,2]
    upper.arg <- ma.dat[,3]
  }
  fp <- NULL
  if (plot) {
    fp <- forestplot::forestplot(labeltext  = labeltext,
                                 mean       = mean.arg,
                                 lower      = lower.arg,
                                 upper      = upper.arg,
                                 is.summary = c(TRUE, rep(FALSE, x$k), TRUE, TRUE),
                                 hrzl_lines = horizl,
                                 fn.ci_norm = fn.ci_norm,
                                 fn.ci_sum  = fn.ci_sum,
                                 col        = col,
                                 boxsize    = boxsize,
                                 legend     = legend, ...)
    plot(fp)
    # add heterogeneity phrase at bottom left:
    if (heterogeneity) {
      tauFigures <- x$summary[c("median","95% lower", "95% upper"), "tau"]
      tauFigures <- tauFigures[tauFigures > 0]
      formatstring <- paste0("%.", decplaces(tauFigures, digits), "f")
      tauphrase <- sprintf(paste0("Heterogeneity (tau): ",formatstring,
                                  " [",formatstring,", ",formatstring,"]"),
                           x$summary["median","tau"],
                           x$summary["95% lower","tau"], x$summary["95% upper","tau"])
      grid::seekViewport("axis_margin")
      tvp <- grid::viewport(x=grid::unit(0.0, "npc"), y=grid::unit(0.0, "npc"),
                            width=grid::stringWidth(tauphrase),
                            height=grid::unit(2, "lines"),
                            just=c("left","bottom"), name="heterogeneityEstimate")
      grid::pushViewport(tvp)
      grid::grid.text(tauphrase, x=grid::unit(0.0, "npc"), y=grid::unit(0.5,"npc"),
                      just=c("left", "centre"), gp=grid::gpar(fontface="oblique"))
    }
  }
  invisible(list("data"       = ma.dat[-1,],
                 "shrinkage"  = ma.shrink[2:(x$k+1),],
                 "labeltext"  = labeltext,
                 "forestplot" = fp))
}



forestplot.escalc <- function(x, labeltext,
                              exponentiate  = FALSE,
                              digits        = 2,
                              plot          = TRUE,
                              fn.ci_norm, fn.ci_sum, col, boxsize, ...)
#
# ARGUMENTS:
#   x            :  a "bayesmeta" object.
#   labeltext    :  you may provide an alternative "labeltext" argument here
#                   (see also the "forestplot()" help).
#   exponentiate :  flag indicating whether to exponentiate numbers (figure and table).
#   prediction   :  flag indicating whether to show prediction interval.
#   shrinkage    :  flag indicating whether to show shrinkage estimates.
#   digits       :  number of significant digits to be shown (based on standard deviations).
#   plot         :  flag you can use to suppress actual plotting.
#   ...          :  further arguments passed to the "forestplot" function.
#   
# VALUE:
#   a list with components
#     $data       :  the meta-analyzed data, and mean and prediction estimates
#     $labeltext  :  the "forestplot()" function's "labeltext" argument used internally
#     $forestplot :  the "forestplot()" function's returned value
#
{
  if (!requireNamespace("forestplot", quietly=TRUE))
    stop("required 'forestplot' package not available!")
  if (utils::packageVersion("forestplot") < "1.5.2")
    warning("you may need to update 'forestplot' to a more recent version (>=1.5.2).")
  # auxiliary function:
  decplaces <- function(x, signifdigits=3)
  # number of decimal places (after decimal point)
  # to be displayed if you want at least "signifdigits" digits
  {
    return(max(c(0, -(floor(log10(x))-(signifdigits-1)))))
  }
  # some sanity checks for the provided arguments:
  stopifnot(is.element("escalc", class(x)),
            length(digits)==1, digits==round(digits), digits>=0,
            length(exponentiate)==1, is.logical(exponentiate),
            length(plot)==1, is.logical(plot))
  
  # extract relevant column names:
  attri <- attributes(x)
  if (all(is.element(c("yi.names", "vi.names"), names(attri)))) { # (for recent "metafor" versions)
    var.names <- c(attri$yi.names, attri$vi.names)
  } else if (is.element("var.names", names(attri))) {             # (for older "metafor" versions)
    var.names <- attri$var.names
  } else {
    stop(paste("Cannont extract \"yi.names\" and \"vi.names\" (or \"var.names\") attribute(s) from escalc object ",
               "(check use of the \"var.names\" option).", sep=""))
  }
  stopifnot(length(var.names)==2, all(is.character(var.names)))
  if (!all(is.element(var.names, names(x)))) {
    stop(paste("Cannont find columns \"",
               var.names[1],"\" and/or \"",
               var.names[2],"\" in escalc object ",
               "(check use of the \"var.names\" option).", sep=""))
  }
  if (is.element("slab", names(attributes(x[,var.names[1]])))) {
    labels <- as.character(attr(x[,var.names[1]], "slab"))
  } else {
    labels <- sprintf("%02d", 1:nrow(x))
  }
  
  # extract data to be plotted:
  y      <- as.vector(x[,var.names[1]])
  sigma  <- sqrt(as.vector(x[,var.names[2]]))
  k <- length(y)
  q95 <- qnorm(0.975)
  ma.dat <- rbind(NA,
                  cbind(y, y - q95*sigma, y + q95*sigma))
  colnames(ma.dat) <- c("estimate", "lower", "upper")
  rownames(ma.dat) <- c("", labels)
  if (exponentiate) {
    ma.dat    <- exp(ma.dat)
  }
  
  # generate "labeltext" data table for plot (unless already provided):
  if (missing(labeltext)) {
    # determine numbers of digits based on standard deviations:
    if (exponentiate) {
      stdevs <- exp(y)*sigma
    } else {
      stdevs <- sigma
    }
    stdevs <- abs(stdevs[is.finite(stdevs) & (stdevs != 0)])
    formatstring <- paste0("%.", decplaces(stdevs, digits), "f")
    # fill data table:
    labeltext <- matrix(NA_character_, nrow=nrow(ma.dat), ncol=3)
    labeltext[1,] <- c("study","estimate", "95% CI")
    labeltext[,1] <- c("study", labels)[1:nrow(ma.dat)]
    for (i in 2:(nrow(ma.dat))) {
      labeltext[i,2] <- sprintf(formatstring, ma.dat[i,"estimate"])
      labeltext[i,3] <- paste0("[", sprintf(formatstring, ma.dat[i,"lower"]),
                                 ", ", sprintf(formatstring, ma.dat[i,"upper"]), "]")
    }
  }
  # add horizontal lines to plot:
  horizl <- list(grid::gpar(col="grey"))
  names(horizl) <- as.character(c(2))
  # specify function(s) for drawing estimates / shrinkage estimates:
  if (missing(fn.ci_norm)) {
    fn.ci_norm <- function(...) {forestplot::fpDrawPointCI(pch=15,...)}
  }
  # specify function(s) for drawing summaries (diamond / bar):
  if (missing(fn.ci_sum)) {
    fn.ci_sum <- list(NULL)
    for (i in 1:(k+1))
      fn.ci_sum[[i]] <- function(y.offset,...) {forestplot::fpDrawSummaryCI(y.offset=0.5,...)}
  }
  # specify colors:
  if (missing(col)) {
    col <- forestplot::fpColors(box="black", lines="black", summary="grey30")
  }
  # specify plotting sizes:
  if (missing(boxsize)) {
    boxsize <- rep(0.25,k+1)
  }
  # specify data for plotting:
  mean.arg  <- ma.dat[,1]
  lower.arg <- ma.dat[,2]
  upper.arg <- ma.dat[,3]
  fp <- NULL
  if (plot) {
    fp <- forestplot::forestplot(labeltext  = labeltext,
                                 mean       = mean.arg,
                                 lower      = lower.arg,
                                 upper      = upper.arg,
                                 is.summary = c(TRUE, rep(FALSE, k)),
                                 hrzl_lines = horizl,
                                 fn.ci_norm = fn.ci_norm,
                                 fn.ci_sum  = fn.ci_sum,
                                 col        = col,
                                 boxsize    = boxsize, ...)
    plot(fp)
  }
  invisible(list("data"       = ma.dat[-1,],
                 "labeltext"  = labeltext,
                 "forestplot" = fp))
}



convolve <- function(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)
# Convolution function from:
#   C. Roever, T. Friede.
#   Discrete approximation of a mixture distribution via restricted divergence.
#   Journal of Computational and Graphical Statistics, 26(1):217-222, 2017.
#   http://doi.org/10.1080/10618600.2016.1276840
#
# (a grid is constructed in the first density's domain ("dens1")
#  so that the convolution is a sum of "dens2" conditionals,
#  weighted by "dens1".)
{
  # some basic sanity checks:
  stopifnot(is.function(dens1), is.function(dens2),
            is.function(cdf1), is.function(cdf2),
            #dens1(-1)>0, dens2(-1)>0,
            delta>0, epsilon>0)
  
  symKL <- function(d)
  # divergence w.r.t. shifting of the 2nd distribution ("dens2()")
  {
    stopifnot(d>=0)
    func1 <- function(x)
    # integrand for one directed divergence
    {
      d2md <- dens2(x-d)
      # utilize density's "log" argument, if possible:
      if (is.element("log", names(formals(dens2)))) {
        logd2   <- dens2(x, log=TRUE)
        logd2md <- dens2(x-d, log=TRUE)
      } else {
        logd2   <- log(dens2(x))
        logd2md <- log(d2md)
      }
      return(ifelse((logd2>-Inf) & (logd2md>-Inf), (logd2md-logd2)*d2md, 0.0))
    }
    
    func2 <- function(x)
    # integrand for other directed divergence
    {
      d2 <- dens2(x)
      # utilize density's "log" argument, if possible:
      if (is.element("log", names(formals(dens2)))) {
        logd2   <- dens2(x, log=TRUE)
        logd2md <- dens2(x-d, log=TRUE)
      } else {
        logd2   <- log(d2)
        logd2md <- log(dens2(x-d))
      }
      return(ifelse((logd2>-Inf) & (logd2md>-Inf), (logd2-logd2md)*d2, 0.0))
    }
    int1 <- integrate(func1, -Inf, Inf)
    if (int1$message != "OK")
      warning(paste0("Problem computing KL-divergence (1): \"", int1$message,"\""))
    int2 <- integrate(func2, -Inf, Inf)
    if (int2$message != "OK")
      warning(paste0("Problem computing KL-divergence (2): \"", int2$message,"\""))
    return(int1$value + int2$value)
  }

  # determine bin half-width:
  step <- sqrt(delta)
  while (symKL(step) < delta) step <- 2*step
  ur <- uniroot(function(X){return(symKL(X)-delta)}, lower=0, upper=step)
  step <- ur$root
  # determine grid range via  epsilon/2  and  1-epsilon/2  quantiles:
  mini <- -1
  while (cdf1(mini) > epsilon/2) mini <- 2*mini
  maxi <- 1
  while (cdf1(maxi) < 1-(epsilon/2)) maxi <- 2*maxi
  ur <- uniroot(function(X){return(cdf1(X)-epsilon/2)},
                lower=mini, upper=maxi)
  mini <- ur$root
  ur <- uniroot(function(X){return(cdf1(X)-(1-epsilon/2))},
                lower=mini, upper=maxi)
  maxi <- ur$root
  # determine number of reference points:
  k <- ceiling((maxi-mini)/(2*step))+1
  # determine reference points:
  support <- mini - ((k*2*step)-(maxi-mini))/2 + (0:(k-1))*2*step
  # determine bin margins:
  margins <- support[-1]-step
  # determine bin weights:
  weight <- rep(NA, length(support))
  for (i in 1:(k-1))
    weight[i] <- cdf1(margins[i])
  weight[k] <- 1
  for (i in k:2)
    weight[i] <- weight[i]-weight[i-1]
  # the eventual grid:
  grid <- cbind("lower"=c(-Inf, margins),
                "upper"=c(margins, Inf),
                "reference"=support,
                "prob"=weight)

  # probability density of convolution:
  density <- function(x)
  {
    return(apply(matrix(x,ncol=1), 1,
                 function(x){sum(grid[,"prob"]*dens2(x-grid[,"reference"]))}))
  }

  # cumulative distribution function (CDF) of convolution:
  cdf <- function(x)
  {
    return(apply(matrix(x,ncol=1), 1,
                 function(x){sum(grid[,"prob"]*cdf2(x-grid[,"reference"]))}))
  }

  # quantile function (inverse CDF) of convolution:
  quantile <- function(p)
  {
    quant <- function(pp)
    {
      mini <- -1
      while (cdf(mini) > pp) mini <- 2*mini
      maxi <- 1
      while (cdf(maxi) < pp) maxi <- 2*maxi
      ur <- uniroot(function(x){return(cdf(x)-pp)}, lower=mini, upper=maxi)
      return(ur$root)      
    }
    proper <- ((p>0) & (p<1))
    result <- rep(NA,length(p))
    if (any(proper)) result[proper] <- apply(matrix(p[proper],ncol=1), 1, quant)
    return(result)
  }
  
  return(list("delta"    = delta,     # tuning parameter (maximum divergence)
              "epsilon"  = epsilon,   # tuning parameter (neglected tail probability)
              "binwidth" = 2*step,    # bin width (constant across bins)
              "bins"     = k,         # number of bins
              "support"  = grid,      # bin margins, reference points and weights
              "density"  = density,   # resulting probability density function
              "cdf"      = cdf,       # resulting cumulative distribution function (CDF)
              "quantile" = quantile)) # resulting quantile function (inverse CDF)
}



funnel.bayesmeta <- function(x,
                             main=deparse(substitute(x)),
                             xlab=expression("effect "*y[i]),
                             ylab=expression("standard error "*sigma[i]),
                             zero=0.0, FE=FALSE, legend=FE, ...)
# generate funnel plot from a "bayesmeta" object.
{
  # sanity check:
  stopifnot(is.element("bayesmeta", class(x)))
  # range of standard error axis:
  yrange <- c(0.0, max(x$sigma))
  # standard error levels for prediction intervals:
  sevec <- seq(from=0, yrange[2]*1.04, le=27)
  # compute (RE) prediction intervals
  # (for observed "y" values, conditional on standard error):
  intRE <- matrix(NA_real_, nrow=length(sevec), ncol=2,
                dimnames=list(NULL, c("lower","upper")))
  intRE[1,] <- x$qposterior(theta.p=c(0.025, 0.975), predict=TRUE)
  for (i in 2:length(sevec)){
    conv <- try(convolve(dens1=function(a, log=FALSE){return(x$dposterior(theta=a, predict=TRUE, log=log))},
                         dens2=function(b, log=FALSE){return(dnorm(x=b, mean=0, sd=sevec[i], log=log))},
                         cdf1 =function(a){return(x$pposterior(theta=a, predict=TRUE))},
                         cdf2 =function(b){return(pnorm(q=b, mean=0, sd=sevec[i]))}))
    if (all(class(conv)!="try-error")) {
      intRE[i,] <- conv$quantile(p=c(0.025, 0.975))
    }
  }
  # compute _FE_ prediction intervals
  # (for observed "y" values, conditional on standard error):
  intFE <- matrix(NA_real_, nrow=length(sevec), ncol=2,
                  dimnames=list(NULL, c("lower","upper")))
  cm <- x$cond.moment(tau=0)
  for (i in 1:length(sevec)){
    intFE[i,] <- qnorm(c(0.025, 0.975), mean=cm[1,"mean"], sd=sqrt(cm[1,"sd"]^2+sevec[i]^2))
  }
  FEcol="red3"
  REcol="blue3"
  # empty plot:
  plot(range(intRE), -yrange, type="n",
       ylab=ylab, xlab=xlab, main=main, axes=FALSE)
  # funnels (grey area):
  polygon(c(intRE[,1], rev(intRE[,2])), c(-sevec, rev(-sevec)), col="grey90", border=NA)
  if (FE) polygon(c(intFE[,1], rev(intFE[,2])), c(-sevec, rev(-sevec)), col="grey80", border=NA)
  # grid lines:
  lines(c(intRE[1,1], intRE[1,1], NA, intRE[1,2], intRE[1,2]),
        c(0,-max(sevec), NA, 0, -max(sevec)), col="grey75", lty="dashed")
  abline(h=0, col="darkgrey")
  yticks <- pretty(yrange)
  abline(h=-yticks[yticks>0], col="grey75", lty="15")
  # funnels (outline):
  matlines(intRE, cbind(-sevec, -sevec), col=REcol, lty="dashed")
  if (FE) matlines(intFE, cbind(-sevec, -sevec), col=FEcol, lty="dotted")
  lines(rep(x$summary["median","theta"], 2), range(-sevec), col=REcol, lty="dashed")
  if (FE) lines(rep(cm[1,"mean"], 2), range(-sevec), col=FEcol, lty="dotted")
  # zero line:
  if (is.finite(zero))
    lines(c(zero, zero), c(-1,1)*max(sevec), col="darkgrey", lty="solid")
  # actual points:
  points(x$y, -x$sigma, pch=21, col="black", bg=grey(0.5, alpha=0.5), cex=1)
  if (FE && legend)
    legend("topleft", c("RE model", "FE model"),
           col=c(REcol, FEcol), lty=c("dashed", "dotted"), bg="white")
  axis(1); axis(2, at=-yticks, labels=yticks); box()
  invisible()
}



normalmixture <- function(density,
                          cdf = Vectorize(function(x){integrate(density,0,x)$value}),
                          mu = 0,
                          delta = 0.01, epsilon = 0.0001,
                          rel.tol.integrate=2^16*.Machine$double.eps,
                          abs.tol.integrate=rel.tol.integrate,
                          tol.uniroot=rel.tol.integrate)
# derive normal mixture where mean (mu) is given (and fixed)
# and the standard deviation (sigma) follows a distribution
# given through "density" or "cdf" argument.
# Mixture approximation is done using the 'DIRECT' algorithm
# described in http://arxiv.org/abs/1602.04060
{
  # some basic sanity checks:
  stopifnot((!missing(cdf) || !missing(density)),
            is.function(cdf), delta>0, epsilon>0)
  if (missing(cdf) && !missing(density)) {
    stopifnot(is.function(density))
    # check for properness of mixing distribution:
    density.integral <- integrate(density, lower=0, upper=Inf,
                                  rel.tol=rel.tol.integrate, abs.tol=abs.tol.integrate)
    stopifnot(density.integral$message == "OK",
              (abs(density.integral$value - 1.0) <= density.integral$abs.error))
  }

  nextsigma <- function(sigma=1, delta=0.01)
  # analytical solution for next larger sigma value
  # for which KL-divergence w.r.t. current sigma is = delta
  {
     return(sigma * sqrt(1 + delta^2 + sqrt(delta^2 + 2*delta)))
  }

  # determine (minimally required) grid range:
  s <- 1
  while (cdf(s) <= epsilon/2) s <- s*2
  ur <- uniroot(function(x){cdf(x)-epsilon/2}, lower=0, upper=s, tol=tol.uniroot)
  lower <- ur$root
  while (cdf(s) <= 1-epsilon/2) s <- s*2
  ur <- uniroot(function(x){cdf(x)-(1-epsilon/2)}, lower=lower, upper=s, tol=tol.uniroot)
  upper <- ur$root

  # now fill grid:
  grid <- matrix(c(0, NA, lower, NA),
                 nrow=1,
                 dimnames=list(NULL, c("lower","upper","reference","prob")))
  s <- lower
  i <- 1
  while (s < upper) {
    i <- i+1
    grid <- rbind(grid, NA)
    s <- nextsigma(s, delta=delta)
    grid[i-1,"upper"] <- grid[i,"lower"] <- s
    s <- nextsigma(s, delta=delta)
    grid[i,"reference"] <- s
  }
  grid[i, "upper"] <- Inf

  grid[,"prob"] <- diff(c(0, cdf(grid[,"upper"])))
  grid[,"prob"] <- grid[,"prob"] / sum(grid[,"prob"])

  # probability density of mixture:
  dens <- function(x)
  {
    return(apply(matrix(x,ncol=1), 1,
                 function(y){sum(grid[,"prob"]*dnorm(y, mean=mu, sd=grid[,"reference"]))}))
  }

  # cumulative distribution function (CDF) of mixture:
  cumul <- function(x)
  {
    return(apply(matrix(x,ncol=1), 1,
                 function(y){sum(grid[,"prob"]*pnorm(y, mean=mu, sd=grid[,"reference"]))}))
  }

  # quantile function (inverse CDF) of mixture:
  quantile <- function(p)
  {
    quant <- function(pp)
    {
      mini <- mu-1
      while (cumul(mini) > pp) mini <- mu - 2*(mu-mini)
      maxi <- mu+1
      while (cumul(maxi) < pp) maxi <- mu + 2*(maxi-mu)
      ur <- uniroot(function(x){return(cumul(x)-pp)}, lower=mini, upper=maxi, tol=tol.uniroot)
      return(ur$root)      
    }
    proper <- ((p>0) & (p<1))
    result <- rep(NA,length(p))
    if (any(proper)) result[proper] <- apply(matrix(p[proper],ncol=1), 1, quant)
    return(result)
  }  
  
  result <- list("delta"          = delta,
                 "epsilon"        = epsilon,
                 "mu"             = mu,
                 "bins"           = i,
                 "support"        = grid,
                 "density"        = dens,      # resulting probability density function
                 "cdf"            = cumul,     # resulting cumulative distribution function (CDF)
                 "quantile"       = quantile,  # resulting quantile function (inverse CDF)
                 "mixing.density" = NULL,
                 "mixing.cdf"     = cdf)  
  if (!missing(density)) result$mixing.density <- density
  return(result)
}


pppvalue<- function(x,
                    parameter = "mu",
                    value = 0.0,
                    alternative = c("two.sided", "less", "greater"),
                    statistic = "median",  # a.k.a. "discrepancy variable"
                    rejection.region,
                    n = 10,
                    prior = FALSE,
                    quietly = FALSE,
                    parallel, seed, ...)
# Posterior predictive p-values
#    * Gelman & al., Bayesian Data Analysis, 3rd edition, Chapter 6,
#      Chapman & Hall / CRC, Boca Raton, 2014.
#    * Meng, X.-L. Posterior predictive p-values.
#      The Annals of Statistics, 22(3):1142-1160, 1994.
#      http://doi.org/10.1214/aos/1176325622
#
# Parameters:
#   x                :  a "bayesmeta" object
#   parameter        :  the parameter to be tested
#   value            :  the null-hypothesized value
#   alternative      :  the alternative to be tested against
#   statistic        :  the figure to be used as "test statistic"
#   rejection.region :  the test statistic's rejection region. May be
#                       one of "upper.tail", "lower.tail" or "two.tailed".
#                       If unspecified, it is set based on the "alternative" parameter
#   n                :  the number of Monte Carlo samples
#   prior            :  flag to request _PRIOR_ predictive p-values
#   quietly          :  flag to indicate command line text output
#   parallel         :  number of parallel processes to use
#   ...              :  further arguments passed to "statistic",
#                       if "statistic" argument is a function
#
{
  # some preliminary sanity checks;
  # "x" argument:
  if (!is.element("bayesmeta",class(x)))
    warning("function applicable to objects of class 'bayesmeta' only.")
  stopifnot(is.element("bayesmeta",class(x)))
  # "parameter" argument:
  stopifnot(length(parameter)==1)
  if (is.numeric(parameter)){
    stopifnot(is.element(parameter, 1:x$k))
    parameter <- x$labels[parameter]
  }
  parameter <- match.arg(parameter, c("mu", "tau", x$labels))
  thetapar <- FALSE
  if (!is.element(parameter, c("mu","tau"))) {
    thetapar <- TRUE  # indicates that parameter concerns one of the "theta" (shrinkage) parameters
    indiv.which <- which(is.element(x$labels, parameter))  # index of concerned "theta" parameter
  }
  # "value" argument:
  stopifnot(length(value)==1, is.finite(value),
            (parameter != "tau") || (value >= 0.0))
  # "alternative" argument:
  alternative <- match.arg(alternative)
  # "statistic" argument:
  statfun <- is.function(statistic)
  statNA <- (!statfun && is.na(statistic))
  stopifnot(statfun | is.character(statistic) | statNA)
  if (statNA) {
    statname <- "statistic"
  } else if (statfun) {
    statname <- deparse(substitute(statistic))
  } else {
    statistic <- match.arg(tolower(statistic), c("t", "q", "cdf", rownames(x$summary)))
    if ((parameter=="tau") && (value==0.0)) {
      if (alternative != "greater")
        warning(paste0("a combination of 'parameter=\"tau\"' ",
                       "and 'alternative=\"", alternative, "\"' may not make sense!"))
      if (statistic == "cdf")
        warning(paste0("a combination of 'parameter=\"tau\"', ",
                       "'value=0.0' and 'statistic=\"cdf\"' does not make sense!"))
    }
    # try to speed up 'bayesmeta()' computations by omitting unnecessary CI optimizations:
    inttype <- ifelse(is.element(statistic, c("95% lower", "95% upper")),
                      x$interval.type, "central")
    statname <- ifelse(statistic=="q", "Q", statistic)
  }
  # "rejection.region" argument:
  if (!missing(rejection.region)) {
    rejection.region <- match.arg(rejection.region, c("upper.tail", "lower.tail", "two.tailed"))
  } else { # decide on rejection region based on hypothesis:
    if (statNA | (!statfun && is.element(statistic, c("q"))))
      rejection.region <- "upper.tail"
    else if (!statfun && is.element(statistic, c("cdf")))
      rejection.region <- switch(alternative,
                                 two.sided = "two.tailed",
                                 less      = "upper.tail",
                                 greater   = "lower.tail")
    else
      rejection.region <- switch(alternative,
                                 two.sided = "two.tailed",
                                 less      = "lower.tail",
                                 greater   = "upper.tail")
  }
  stopifnot(is.element(rejection.region, c("upper.tail", "lower.tail", "two.tailed")))
  # "prior" argument:
  if (prior && (!x$tau.prior.proper || !all(is.finite(x$mu.prior))))
    warning("prior predictive p-values require proper priors for effect and heterogeneity!")
  if (prior && !is.element(parameter, c("mu", "tau")))
    warning("prior predictive p-values are only available for effect (mu) and heterogeneity (tau) parameters!")
  stopifnot(!prior || (x$tau.prior.proper && all(is.finite(x$mu.prior))))
  # determine a sensible number of parallel processes:
  if (missing(parallel)) {
    # by default, use "parallel" package if available:
    if (requireNamespace("parallel")) {
      # by default, use all but one core:
      parallel <- max(c(1, parallel::detectCores() - 1))
      # but don't use more cores than MCMC samples:
      parallel <- min(c(parallel, n))
    } else {
      parallel <- 1
    }
  } else { # otherwise check number of cores & processes:
    stopifnot(parallel >= 1)
    if (requireNamespace("parallel")) {
      if (parallel > parallel::detectCores())
        warning("number of requested parallel processes (", parallel,
                ") is larger than number of cores (", parallel::detectCores(), ").")
      if (parallel > n)
        warning("number of requested parallel processes (", parallel,
                ") is larger than number of MC samples (", n, ").")
    } else {
      warning("failed to load \"parallel\" package.")
    }
  }
  seed.missing <- missing(seed)
  stopifnot(seed.missing || (seed==round(seed)))
  ptm <- proc.time()
  ##################
  if (prior) {         # ...conditional prior p(tau | mu)
    ptau <- function(tau)
    # cumulative distribution function (CDF) of tau prior
    {
      stopifnot(length(tau)==1, is.finite(tau), tau>=0)
      return(integrate(function(t){x$dprior(tau=t)}, lower=0, upper=tau,
                       rel.tol=x$rel.tol.integrate, abs.tol=x$abs.tol.integrate)$value)
    }
    ptau <- Vectorize(ptau)
    qtau <- function(p.tau)
    # quantile function (inverse CDF) of tau prior
    {
      stopifnot(length(p.tau)==1, is.finite(p.tau), p.tau>=0, p.tau<=1)
      if (p.tau==0) result <- 0.0
      else if (p.tau==1) result <- Inf
      else {
        upper <- 1.0
        while (ptau(upper) < p.tau) upper <- 2*upper
        result <- uniroot(function(t){ptau(tau=t)-p.tau},
                          lower=0, upper=upper, tol=x$tol.uniroot)$root
      }
      return(result)
    }
  }
  if (parameter=="mu") { # define functions to generate draws from conditional (tau|mu):
    # lookup table for tau conditionals' normalizing constants:
    normconst <- matrix(nrow=0, ncol=2, dimnames=list(NULL, c("mu","const")))
    dctau <- function(tau, mu)
    # conditional density of (tau | mu)
    {
      stopifnot(length(tau)==1, is.finite(tau), tau>=0)
      if ((nrow(normconst)>0) && is.element(mu, normconst[,"mu"])) {
        const <- normconst[which(normconst[,"mu"]==mu),"const"]
      } else {
        const <- integrate(function(t){x$dposterior(tau=t, mu=mu)},
                           lower=0, upper=Inf,
                           rel.tol=x$rel.tol.integrate, abs.tol=x$abs.tol.integrate)$value
        normconst <<- rbind(normconst, c(mu,const)) # (change matrix GLOBALLY)
      }
      return(x$dposterior(tau=tau, mu=mu) / const)
    }
    dctau <- Vectorize(dctau)
    pctau <- function(tau, mu)
    # conditional cumulative distribution function (CDF) of (tau | mu)
    {
      stopifnot(length(tau)==1, is.finite(tau), tau>=0)
      return(integrate(function(t){dctau(tau=t, mu=mu)}, lower=0, upper=tau,
                       rel.tol=x$rel.tol.integrate, abs.tol=x$abs.tol.integrate)$value)
    }
    pctau <- Vectorize(pctau)
    qctau <- function(p.tau, mu)
    # conditional quantile function (inverse CDF) of (tau | mu)
    {
      stopifnot(length(p.tau)==1, is.finite(p.tau), p.tau>=0, p.tau<=1)
      if (p.tau==0) result <- 0.0
      else if (p.tau==1) result <- Inf
      else {
        upper <- 1.0
        while (pctau(upper,mu) < p.tau) upper <- 2*upper
        result <- uniroot(function(t){pctau(tau=t,mu=mu)-p.tau},
                          lower=0, upper=upper, tol=x$tol.uniroot)$root
      }
      return(result)
    }
    rctau <- function(mu)
    # random number generation for (tau | mu)
    {
      return(qctau(p.tau=runif(n=1), mu=mu))
    }
  } else if (thetapar) {
    # compute conditional posterior, conditional on theta[i]==value:
    condy     <- x$y
    condsigma <- x$sigma
    condy[indiv.which]     <- value
    condsigma[indiv.which] <- 0.0
    if (alternative == "two.sided") {
      condbm <- bayesmeta(y=condy, sigma=condsigma, labels=x$labels,
                          mu.prior=x$mu.prior,
                          tau.prior=x$dprior,
                          interval.type=inttype,
                          rel.tol.integrate=x$rel.tol.integrate,
                          abs.tol.integrate=x$abs.tol.integrate,
                          tol.uniroot=x$tol.uniroot)
    }
  }
  # determine realized "test statistic" value in actual data:
  if (statNA) {
    stat.actual <- NA_real_
  } else if (statfun) {
    stat.actual <- statistic(x$y, ...)
  } else {
    if (statistic == "q") {
      cochranQ <- function(yi, si)
      {
        wi <- 1/si^2
        yhat <- sum(yi * wi) / sum(wi)
        Q <- sum(((yi-yhat)/si)^2)
        return(Q)
      }
      stat.actual <- cochranQ(x$y, x$sigma)
    } else if (statistic=="cdf") {
      if (parameter=="tau")
        stat.actual <- x$pposterior(tau=value)
      else if (parameter=="mu")
        stat.actual <- x$pposterior(mu=value)
      else
        stat.actual <- x$pposterior(theta=value, individual=indiv.which)
    } else {
      if (thetapar) {
        if (statistic == "t")
          stat.actual <- (x$theta["mean", indiv.which] - value) / x$theta["sd", indiv.which]
        else
          stat.actual <- x$theta[statistic, indiv.which]
      } else {
        if (statistic == "t")
          stat.actual <- (x$summary["mean", parameter] - value) / x$summary["sd", parameter]
        else
          stat.actual <- x$summary[statistic, parameter]
      }
    }
  }
  names(stat.actual) <- statname
  # determine at which (prior/posterior) quantile hypothesized value is situated
  # (necessary for sampling for single-tailed hypotheses):
  if (alternative != "two.sided") {
    if (parameter == "mu") {
      if (prior) p.mu <- pnorm(q=value, mean=x$mu.prior[1], sd=x$mu.prior[2])  # prior quantile
      else       p.mu <- x$pposterior(mu=value)                                # posterior quantile
    } else if (parameter == "tau") {
      if (prior) p.tau <- ptau(value)              # prior quantile
      else       p.tau <- x$pposterior(tau=value)  # posterior quantile
    } else {
      p.theta <- x$pposterior(theta=value, individual=indiv.which)  # posterior quantile
    }
  }
  
  # the main function:
  pppfun <- function(i)
  {
    if (! seed.missing) set.seed(seed + i)
    # generate data:
    sigma <- x$sigma
    if (parameter=="mu") {          # Null hypothesis concerns effect mu
      if (alternative=="two.sided") {   # fix effect mu at hypothesized value:
        rmu <- value
      } else if (alternative=="less") { # draw effect mu from H0's domain's conditional:
        if (prior) rmu <- qnorm(runif(1, p.mu, 1.0), mean=x$mu.prior[1], sd=x$mu.prior[2])  # prior
        else       rmu <- x$qposterior(mu.p=runif(1, p.mu, 1.0))                            # posterior
      } else {
        if (prior) rmu <- qnorm(runif(1, 0.0, p.mu), mean=x$mu.prior[1], sd=x$mu.prior[2])  # prior
        else       rmu <- x$qposterior(mu.p=runif(1, 0.0, p.mu))                            # posterior
      }
      # draw tau from conditional (tau|mu):
      if (prior) rtau <- qtau(runif(1, 0.0, 1.0))  # prior
      else       rtau <- rctau(mu=rmu)             # posterior
      # draw study-specific effects (theta) and effects (y):
      rtheta <- rnorm(n=x$k, mean=rmu,    sd=rtau)
      y      <- rnorm(n=x$k, mean=rtheta, sd=sigma)
    } else if (parameter=="tau") {  # Null hypothesis concerns heterogeneity tau
      if (alternative=="two.sided") {   # fix heterogeneity tau at hypothesized value:
        rtau <- value
      } else if (alternative=="less") { # draw effect mu from H0's domain's conditional:
        if (prior) rtau <- qtau(runif(1, p.tau, 1.0))
        else       rtau <- x$qposterior(tau.p=runif(1, p.tau, 1.0))
      } else {
        if (prior) rtau <- qtau(runif(1, 0.0, p.tau))
        else       rtau <- x$qposterior(tau.p=runif(1, 0.0, p.tau))
      }
      # draw mu from conditional (mu|tau):
      if (prior) {
        rmu <- qnorm(runif(1, 0.0, 1.0), mean=x$mu.prior[1], sd=x$mu.prior[2])
      } else {
        cm <- x$cond.moment(tau=rtau)
        rmu <- rnorm(n=1, mean=cm[1,"mean"], sd=cm[1,"sd"])
      }
      # draw study-specific effects (theta) and effects (y):
      rtheta <- rnorm(n=x$k, mean=rmu,    sd=rtau)
      y      <- rnorm(n=x$k, mean=rtheta, sd=sigma)
    } else {                        # Null hypothesis concerns ith "shrinkage" parameter theta
      if (alternative=="two.sided"){
        rtheta.i <- value
        taumu <- condbm$rposterior(n=1)[1,]
      } else {
        if (alternative=="less") {
          rtheta.i <- x$qposterior(theta.p=runif(1, p.theta, 1.0), indiv=indiv.which)
        } else {
          rtheta.i <- x$qposterior(theta.p=runif(1, 0.0, p.theta), indiv=indiv.which)
        }
        condy[indiv.which] <- rtheta.i
        condbm <- bayesmeta(y=condy, sigma=condsigma, labels=x$labels,
                            mu.prior=x$mu.prior,
                            tau.prior=x$dprior,
                            interval.type=inttype,
                            rel.tol.integrate=x$rel.tol.integrate,
                            abs.tol.integrate=x$abs.tol.integrate,
                            tol.uniroot=x$tol.uniroot)
        taumu <- condbm$rposterior(n=1)[1,]
      }
      rtau <- taumu["tau"]
      rmu <- taumu["mu"]
      rtheta <- rnorm(n=x$k, mean=rmu, sd=rtau)
      rtheta[indiv.which] <- rtheta.i
      #y      <- rnorm(n=x$k, mean=rmu, sd=sigma)
      y      <- rnorm(n=x$k, mean=rtheta, sd=sigma)
    }
    # data (y) generated. Now compute statistic:
    if (statNA) {
      statval <- NA_real_
    } else if (statfun) {
      statval <- statistic(y, ...)
    } else {
      if (statistic == "q") {
        statval <- cochranQ(y, x$sigma)
      } else {
        # perform meta-analysis:
        bm <- bayesmeta(y=y, sigma=x$sigma, labels=x$labels,
                        mu.prior=x$mu.prior,
                        tau.prior=x$dprior,
                        interval.type=inttype,
                        rel.tol.integrate=x$rel.tol.integrate,
                        abs.tol.integrate=x$abs.tol.integrate,
                        tol.uniroot=x$tol.uniroot)
        if (thetapar) {
          if (statistic == "t")
            statval <- (bm$theta["mean", indiv.which] - value) / bm$theta["sd", indiv.which]
          else if (statistic=="cdf")
            statval <- bm$pposterior(theta=value, individual=indiv.which)
          else
            statval <- bm$theta[statistic, indiv.which]
        } else {
          if (statistic == "t")
            statval <- (bm$summary["mean", parameter] - value) / bm$summary["sd", parameter]
          else if (statistic=="cdf") {
            if (parameter=="tau")
              statval <- bm$pposterior(tau=value)
            else if (parameter=="mu")
              statval <- bm$pposterior(mu=value)
          } else
            statval <- bm$summary[statistic, parameter]
        }
      }
    }
    # statistic computed. Return results:
    result <- unname(c(rtau, rmu, statval, rtheta, y))
    names(result) <- c("tau", "mu", "statistic",
                       sprintf("theta[%d]", 1:x$k),
                       sprintf("y[%d]", 1:x$k))
    return(result)
  } # END pppfun()
    
  # break calculations down into smaller bits
  # to allow for progress bar updates (unless suppressed)
  if (quietly) # (do computations in one go)
    idxlist <- list(1:n)
  else {
    # determine intermediate breakpoints (to update progress bar)
    # at multiples of 'parallel':
    stopfrac <- c(0.01, 0.02, 0.05, (1:9)/10)  # (1%, 2%, 5%, 10%, ...)
    stopint <- unique(ceiling(stopfrac * (n/parallel)))
    stopint <- c(stopint * parallel, n)
    stopint <- unique(stopint[stopint <= n])
    # assemble list of indices to process at each stage:
    idxlist <- list(1:stopint[1])
    if (length(stopint) > 1)
      for (i in 2:length(stopint))
        idxlist[[i]] <- (max(idxlist[[i-1]])+1):stopint[i]
  }
  
  #  the "hot loop":
  stat.repli <- NULL
  if (parallel > 1) { # initialize cluster:
    clust <- parallel::makeCluster(parallel)
    #parallel::clusterEvalQ(clust, library("bayesmeta", lib.loc="~/temp/test"))  #  <--  /!\  HACK!
    parallel::clusterEvalQ(clust, library("bayesmeta"))
  }
  if (!quietly) {      # (initialize progress bar etc.)
    cat(paste0("  Generating n=",n," Monte Carlo samples.\n"))
    if (n <= 100)
      cat(paste0("  /!\\  Caution: a sample size of  n >> 100  will usually be appropriate.\n"))
    cat(paste0("  Sampling progress",
               ifelse(parallel > 1, paste0(" (using ",parallel," parallel processes)"), ""),
               ":\n"))
    pb <- utils::txtProgressBar(style=3)
    utils::setTxtProgressBar(pb, 0.0)
  }
  for (i in 1:length(idxlist)){
    if (parallel == 1) { # simple "sapply()" call
      stat.repli <- rbind(stat.repli,
                          t(sapply(idxlist[[i]], pppfun, simplify=TRUE)))
    } else {             # parallel computation via "parSapply()":
      stat.repli <- rbind(stat.repli,
                          t(parallel::parSapply(clust, idxlist[[i]], pppfun, simplify=TRUE)))
    }
    if (!quietly) utils::setTxtProgressBar(pb, max(idxlist[[i]])/n)
  }
  if (parallel > 1) { # terminate cluster:
    parallel::stopCluster(clust)
  }
  
  # how many replicates are in the "tails" beyond the actualized statistic value:
  tail <- factor(rep("no", nrow(stat.repli)),
                 levels=c("lower", "no", "upper"), ordered=TRUE)
  lowertail <- (stat.repli[,"statistic"] <= stat.actual)
  uppertail <- (stat.repli[,"statistic"] >= stat.actual)
  if (rejection.region=="lower.tail")
    tail[lowertail] <- "lower"
  else if (rejection.region=="upper.tail")
    tail[uppertail] <- "upper"
  else {  # note: two-sided version is based on equal-tailed rejection region
    if (sum(uppertail, na.rm=TRUE) < sum(lowertail, na.rm=TRUE)) {
      tail[uppertail] <- "upper"
      tail[rank(stat.repli[,"statistic"], ties.method="min") <= sum(uppertail, na.rm=TRUE)] <- "lower"
    } else {
      tail[lowertail] <- "lower"
      tail[rank(stat.repli[,"statistic"], ties.method="max") > (n-sum(lowertail, na.rm=TRUE))] <- "upper"
    }
  }
  tail[is.na(stat.repli[,"statistic"])] <- NA
  n.tail <- sum(tail != "no", na.rm=TRUE) + sum(is.na(tail))
  # corresponding p-value:
  quant <- ifelse(statNA, NA_real_, n.tail / n)
  # null hypothesis:
  nullval <- value
  names(nullval) <- ifelse(thetapar,
                           paste0("study effect (",indiv.which,": '",x$labels[indiv.which],"')"),
                           ifelse(parameter=="mu", "effect (mu)", "heterogeneity (tau)"))
  replicates <- list("tau"       = stat.repli[,"tau"],
                     "mu"        = stat.repli[,"mu"],
                     "theta"     = stat.repli[,(3+(1:x$k))],
                     "y"         = stat.repli[,(3+x$k+(1:x$k))],
                     "statistic" = cbind.data.frame(stat.repli[,"statistic"], "tail"=tail))
  colnames(replicates$theta) <- colnames(replicates$y) <- x$labels
  colnames(replicates$statistic) <- c(statname, "tail")
  
  ptm <- proc.time() - ptm
  if (!quietly) cat(paste0("\n  (computation time: ",
                          sprintf("%.1f",ptm[3]), " seconds = ",
                          sprintf("%.1f",ptm[3]/60), " minutes.)\n"))
  
  result <- list("statistic"        = stat.actual,
                 "parameter"        = c("Monte Carlo replicates"=n),
                 "p.value"          = quant,
                 "null.value"       = nullval,
                 "alternative"      = alternative,
                 "method"           = paste0("'bayesmeta' ",
                                             ifelse(prior, "prior", "posterior"), " predictive p-value (",
                                             ifelse(alternative=="two.sided", "two-sided", "one-sided"), ")"),
                 "data.name"        = deparse(substitute(x)),
                 # nonstandard-"htest"-elements:
                 "call"             = match.call(expand.dots=FALSE),
                 "rejection.region" = rejection.region,
                 "replicates"       = replicates,
                 "computation.time" = c("seconds"=unname(ptm[3])))
  class(result) <- "htest"
  return(result)
}



uisd <- function(n, ...)
{
  UseMethod("uisd")
}


uisd.default <- function(n, sigma, sigma2=sigma^2, labels=NULL, individual=FALSE, ...)
# Compute unit information standard deviation (UISD).
# Arguments:
#   nc         :  studies' sample sizes
#   sigma      :  studies' standard errors
#   sigma2     :  studies' variances (squared standard errors)
#   individual :  if `TRUE', study-specific UISDs are returned
#                 (instead of overall)
{
  stopifnot(length(n)==length(sigma2),
            all(is.finite(n)), all(is.finite(sigma2)),
            all(n>0), all(sigma2>0),
            (!individual) |
             ((length(labels)==0) || (length(labels)==length(sigma2))))
  
  if (individual) {
    if (!is.character(labels))
      labels <- as.character(labels)
    result <- sqrt(n*sigma2)
    names(result) <- labels
  } else {
    result <- sqrt(sum(n) / sum(1/sigma2))
  }
  return(result)
}


uisd.escalc <- function(n, ...)
# Compute unit information standard deviation (UISD).
# Arguments:
#   n :  an "escalc" object
{
  stopifnot(is.element("escalc", class(n)))
  if (!is.element("ni", names(attributes(n$yi))))
    stop("An \"ni\" attribute is required for the \"escalc\" object!")
  return(uisd(n=attr(n$yi, "ni"), sigma2=n$vi,
              labels=attr(n$yi, "slab"), ...))
}


ess <- function(object, ...)
{
  UseMethod("ess")
}


ess.bayesmeta <- function(object, uisd,
                          method=c("elir", "vr", "pr", "mtm.pt"), ...)
{
  # Computation of effective sample sizes; see:
  #
  #   B. Neuenschwander, S. Weber, H. Schmidli, A. O'Hagan.
  #   Predictively consistent prior effective sample sizes.
  #   Biometrics 76(2): 578-587, 2020.
  #   https://doi.org/10.1111/biom.13252
  #
  # Note that  i_F(theta)  here equals  1 / UISD^2
  # (where the UISD may or may not vary with theta).

  # preliminary sanity checks for provided arguments:
  stopifnot(is.element("bayesmeta", class(object)))
  if (missing(uisd)) {
    stop("\"uisd\" argument need to be specified (either a numeric value or a function)!")
  }
  stopifnot(is.vector(uisd) | is.function(uisd))
  if (is.vector(uisd)) { # specify function returning a constant:
    stopifnot(length(uisd)==1, is.numeric(uisd),
              is.finite(uisd), uisd > 0)
    uisdFun <- function(theta) {return(rep(uisd, length(theta)))}
  } else {               # specify a wrapper function (with built-in sanity checks):
    uisdFun <- function(theta)
    {
      result <- uisd(theta)
      if (!is.vector(result) | !is.numeric(result)) {
        warning(paste0("Invalid output from \"uisd()\" function: uisd(",theta,")"))
      }
      if (!is.finite(result)) {
        warning(paste0("\"uisd()\" function needs to return finite values! (uisd(",theta,") <= 0)"))
      }
      if (result <= 0.0) {
        warning(paste0("\"uisd()\" function needs to return positive values! (uisd(",theta,") <= 0)"))
      }
      return(result)
    }
    uisdFun <- Vectorize(uisdFun, "theta")
  }

  ####################
  # determine method:
  method <- match.arg(tolower(method), c("elir", "vr", "pr", "mtm.pt"))
  # switch:
  if (method == "elir") { #  "expected local-information-ratio (ELIR)" method:
    # specify function "i(p(theta))"
    # (equation (3) in Neuenschwander & al. (2020)):
    ip <- function(object, theta)
    {
      stopifnot(length(theta)==1, is.finite(theta))
      # compute Hessian:
      hessi <- hessian(function(x){object$dposterior(theta=x,
                                                     predict=TRUE,
                                                     log=TRUE)},
                       theta)
      # note the slight hack here:
      margin <- 6 * diff(object$summary[c("95% lower","95% upper"),"theta"])
      # (this should -roughly- correspond to a 20-sigma margin)
      if ((!is.finite(hessi))
          && (abs(theta-object$summary["median","theta"]) > margin)) {
        hessi <- 0.0
      }
      # (Hessian is set to zero in case of numerical problems
      #  AND a ~ 20-sigma difference from median)
      return(as.vector(-hessi))
    }
    ip <- Vectorize(ip, "theta")

    # compute expectation  (equation (7)):
    integrand <- function(x)
    {
      # NB: For numerical ease, the integrand is computed in a structured way.
      #     The integrand results as a product of 3 factors.
      #     Factors are computed one-by-one, _UNLESS_ one of the factors
      #     turned out as zero already.
      factors <- matrix(0.0, nrow=length(x), ncol=3,
                        dimnames=list(NULL, c("ip", "uisd2", "density")))
      # first, compute density:
      factors[,"density"] <- object$dposterior(theta=x, predict=TRUE)
      # second, compute i(p(theta)):
      zero <- (factors[,"density"] == 0.0)
      if (any(!zero)) {
        factors[!zero, "ip"] <- ip(object, theta=x[!zero])
      }
      # third, compute uisd^2:
      zero <- (factors[,"ip"] == 0.0)
      if (any(!zero)) {
        factors[!zero, "uisd2"] <- uisdFun(x[!zero])^2
      }
      # multiply 3 factors:
      result <- apply(factors, 1, prod)
      return(result)
    }
    expect <- integrate(integrand, lower=-Inf, upper=Inf)
    if (expect$message != "OK")
      warning(paste0("Problem computing expectation (\"", expect$message,"\")."))
    ESS <- expect$value
  } else if (method == "vr") {  # "variance ratio (VR)" method:
    # compute expectation  (equation (2)):
    if (is.vector(uisd)) {
      numerator <- uisd^2
    } else {
      integrand <- function(x)
      {
        # NB: For numerical ease, the integrand is computed in a structured way.
        #     The integrand results as a product of 2 factors (uisd and density).
        #     Factors are computed one-by-one, _UNLESS_ one of the factors
        #     turned out as zero already.
        factors <- matrix(0.0, nrow=length(x), ncol=2,
                          dimnames=list(NULL, c("uisd2", "density")))
        # first, compute density:
        factors[,"density"] <- object$dposterior(theta=x, predict=TRUE)
        # second, compute i(p(theta)):
        zero <- (factors[,"density"] == 0.0)
        if (any(!zero)) {
          factors[!zero, "uisd2"] <- uisdFun(x[!zero])^2
        }
        # multiply 2 factors:
        result <- apply(factors, 1, prod)
        return(result)
      }
      expect <- integrate(integrand, lower=-Inf, upper=Inf)
      if (expect$message != "OK")
        warning(paste0("Problem computing expectation (\"", expect$message,"\")."))
      numerator <- expect$value
    }
    denominator <- object$summary["sd","theta"]^2
    ESS <- numerator / denominator
  } else if (method == "pr") {  # "precision ratio (PR)" method:
    # compute expectation  (equation (2)):
    if (is.vector(uisd)) {
      denominator <- uisd^-2
    } else {
      integrand <- function(x)
      {
        # NB: For numerical ease, the integrand is computed in a structured way.
        #     The integrand results as a product of 2 factors (uisd and density).
        #     Factors are computed one-by-one, _UNLESS_ one of the factors
        #     turned out as zero already.
        factors <- matrix(0.0, nrow=length(x), ncol=2,
                          dimnames=list(NULL, c("uisd2", "density")))
        # first, compute density:
        factors[,"density"] <- object$dposterior(theta=x, predict=TRUE)
        # second, compute i(p(theta)):
        zero <- (factors[,"density"] == 0.0)
        if (any(!zero)) {
          factors[!zero, "uisd2"] <- uisdFun(x[!zero])^-2
        }
        # multiply 2 factors:
        result <- apply(factors, 1, prod)
        return(result)
      }
      expect <- integrate(integrand, lower=-Inf, upper=Inf)
      if (expect$message != "OK")
        warning(paste0("Problem computing expectation (\"", expect$message,"\")."))
      denominator <- expect$value
    }
    numerator <- object$summary["sd","theta"]^-2
    ESS <- numerator / denominator
  } else if (method == "mtm.pt"){  # "Morita-Thall-Mueller / Pennello-Thompson (MTM.PM)" method:
    denominator <- uisdFun(theta=object$summary["mode","theta"])^-2
    hessi <- hessian(function(x){object$dposterior(theta=x,
                                                   predict=TRUE,
                                                   log=TRUE)},
                     object$summary["mode","theta"])
    numerator <- as.numeric(-hessi)
    ESS <- numerator / denominator
  }
  return(ESS)
}


weightsplot <- function(x, ...)
{
  UseMethod("weightsplot")
}

weightsplot.bayesmeta <- function(x, individual=FALSE, ordered=TRUE,
                                  extramargin=4,
                                  priorlabel="prior mean", main, xlim,...)
# Illustrate posterior mean weights (percentages)
# for overall mean or shrinkage estimates in a bar plot
# (See https://doi.org/10.1002/bimj.202000227 ).
# Numbers are taken from the bayesmeta object's
# "...$weights" or "...$weights.theta" elements.
{
  stopifnot(length(ordered)==1, is.logical(ordered),
            is.numeric(extramargin), length(extramargin)==1, all(is.finite(extramargin)),
            length(individual)==1)
  if (! (is.logical(individual) && (!individual))) { #  individual != FALSE
    indiv.logi <- TRUE       # (non-empty "individual" specification)
    if (is.numeric(individual))   indiv.which <- which(is.element(1:x$k, individual))
    if (is.character(individual)) indiv.which <- which(is.element(x$labels, match.arg(individual,x$labels)))
    if (length(indiv.which)==0) warning("cannot make sense of 'individual' argument: empty subset.")
  }
  else indiv.logi <- FALSE   # (the default (overall mean weight))
  mu.prior.proper <- all(is.finite(x$mu.prior))
  # create empty data frame to hold plotted data:
  plotdat <- cbind.data.frame("label"=rep(NA_character_, ifelse(mu.prior.proper, x$k+1, x$k)),
                              "weight"=NA_real_, "percentage"=NA_character_, stringsAsFactors=FALSE)
  if (!indiv.logi){ # fill in data for overall mean estimate
    plotdat[,"label"]  <- names(x$weights)
    plotdat[,"weight"] <- x$weights
    targetparam <- "overall mean estimate"
  } else {          # fill in data for shrinkage estimate
    plotdat[,"label"] <- rownames(x$weights.theta)
    plotdat[,"weight"] <- x$weights.theta[,indiv.which]
    targetparam <- paste0("shrinkage estimate \"", x$labels[indiv.which], "\"")
  }
  if (ordered) {    # sort rows in descending order
    plotdat[1:x$k,] <- plotdat[order(plotdat[1:x$k,"weight"], decreasing=TRUE),]
  }
  plotdat[,"percentage"] <- paste(sprintf("%.1f", plotdat[,"weight"] * 100), "%")
  maxweight <- max(plotdat[,"weight"]*100)
  if (missing(main)) {
    main <- paste0("posterior mean weights (",targetparam,")")
  }
  # set figure margins:
  if (extramargin != 0) {
    parmar <- par("mar")
    on.exit(par(mar=parmar))
    par(mar = parmar + c(0, extramargin, 0, 0))
  }
  if (missing(xlim)) {
    xlim <- c(0, maxweight * 1.15)
  }
  bp <- graphics::barplot(rev(plotdat[,"weight"])*100, horiz=TRUE,
                          names.arg=rev(plotdat[,"label"]), las=1,
                          xlim=xlim,
                          xlab="weight (%)", ylab="", main=main, ...)
  graphics::text(rev(plotdat[,"weight"])*100 + maxweight*0.02, bp[,1],
                 rev(plotdat[,"percentage"]), adj=c(0,0.5))
  invisible(plotdat)
}


traceplot <- function(x, ...)
{
  UseMethod("traceplot")
}


traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
                                ylab="effect",
                                prior=FALSE,
                                infinity=FALSE,
                                rightmargin=8,
                                col=rainbow(x$k), labcol=col,
                                meanlabel="overall mean",
                                meancol="black", meanlabcol=meancol,
                                ...)
{
  stopifnot(missing(mulim) || (length(mulim) == 2),
            missing(taulim) || (length(taulim) <= 2),
            is.character(ylab), length(ylab)==1,
            is.logical(prior), length(prior)==1,
            is.logical(infinity), length(infinity)==1,
            rightmargin >= 0, ((length(col)==x$k) | (length(col)==1)),
            is.character(meanlabel), length(meanlabel)==1,
            length(meancol)==1)
  q975 <- qnorm(0.975)
  gridcol <- "grey85"
  if (length(col)==1) col <- rep(col, x$k)
  if (infinity & any(is.finite(x$mu.prior))) {
    warning("mu prior ignored for `tau=Inf' computations!")
  }
  if (prior & (!x$tau.prior.proper)) {
    warning("Note that plots of improper priors may not be sensibly scaled.")
  }
  # convert "taulim" and "mulim" input arguments
  # to eventual "taurange" and "murange" vectors:
  if (!missing(taulim) && all(is.finite(taulim))) {
    if ((length(taulim)==2) && (taulim[1]>=0) && (taulim[2]>taulim[1]))
      taurange <- taulim
    else if ((length(taulim)==1) && (taulim>0))
      taurange <- c(0, taulim)
    else
      taurange <- c(0, x$qposterior(tau=0.995)*1.1)
  } else {
    taurange <- c(0, x$qposterior(tau=0.995)*1.1)
  }
  if (infinity) {
    xlim <- taurange + c(0, 0.15) * diff(taurange)
    infx <- xlim[2] + 0.04*diff(xlim) # the "infinity" x-coordinate
  } else {
    xlim <- taurange
    infx <- NA_real_
  }

  if (missing(mulim)) mulim <- NULL
  
  vertlines <- pretty(taurange)
  # ensure no tickmarks beyond plotted tau range:
  if (max(vertlines) > (taurange[2] + 0.04*diff(taurange)))
    vertlines <- vertlines[-length(vertlines)]

  mutrace <- function(x)
  {  
    # vector of tau values:
    tau <- seq(max(c(0,taurange[1]-0.04*diff(taurange))),
                   taurange[2]+0.04*diff(taurange), le=200)
    # conditional moments for individual studies (theta_i): 
    cm.indiv   <- x$cond.moment(tau=tau, indiv=TRUE)
    # conditional moments for overall mean (mu): 
    cm.overall <- x$cond.moment(tau=tau)
    # determine axis range for "effect" (y-) axis
    if (!is.null(mulim) && (all(is.finite(mulim)) && (mulim[1] < mulim[2]))) {
      # user-defined:
      murange <- mulim
    } else {
      # based on data:
      if (ci){
        murange <- range(c(range(cm.indiv[,"mean",]-q975*cm.indiv[,"sd",]),
                           range(cm.indiv[,"mean",]+q975*cm.indiv[,"sd",]),
                           range(cm.overall[,"mean"]-q975*cm.overall[,"sd"]),
                           range(cm.overall[,"mean"]+q975*cm.overall[,"sd"])))
      } else {
        murange <- range(c(range(cm.indiv[,"mean",]),
                           range(cm.overall[,"mean"])))
      }
      # ensure that estimates are also included:
      if (infinity) murange <- range(murange, x$y)
    }

    plot(taurange, murange, xlim=xlim,
         type="n", axes=FALSE, xlab="", ylab=ylab, main="", ...)
    abline(v=vertlines, col=gridcol)
    abline(h=pretty(murange), col=gridcol)
    abline(v=0, col=grey(0.40))
    # grey CI shading:
    if (ci) {
      for (i in 1:x$k) {
        polygon(c(tau, rev(tau)),
                c(cm.indiv[,"mean",i] - q975*cm.indiv[,"sd",i],
                  rev(cm.indiv[,"mean",i] + q975*cm.indiv[,"sd",i])),
                col=grey(0.75, alpha=0.25), border=NA)
      }
      polygon(c(tau, rev(tau)),
              c(cm.overall[,"mean"] - q975*cm.overall[,"sd"],
                rev(cm.overall[,"mean"] + q975*cm.overall[,"sd"])),
              col=grey(0.75, alpha=0.25), border=NA)
    }    
    # individual estimates:
    matlines(tau, cm.indiv[,"mean",], col=col, lty=1)
    if (ci) {
      matlines(tau, cm.indiv[,"mean",]-q975*cm.indiv[,"sd",], col=col, lty=3)
      matlines(tau, cm.indiv[,"mean",]+q975*cm.indiv[,"sd",], col=col, lty=3)
    }
    # overall mean:
    lines(tau, cm.overall[,"mean"], col=meancol, lty=2, lwd=1.5)
    if (ci) {
      lines(tau, cm.overall[,"mean"]-q975*cm.overall[,"sd"], col=meancol, lty=3, lwd=1.5)
      lines(tau, cm.overall[,"mean"]+q975*cm.overall[,"sd"], col=meancol, lty=3, lwd=1.5)
    }
    if (infinity) {
      labpos.indiv   <- x$y
      labpos.overall <- mean(x$y)
      for (i in 1:x$k)
        lines(c(max(tau), infx),
              c(cm.indiv[length(tau),"mean",i], labpos.indiv[i]),
              col=col[i], lty="13", lwd=1.5)
      lines(c(max(tau), infx),
            c(cm.overall[length(tau),"mean"], labpos.overall),
            col=meancol, lty="13", lwd=2.0)
    } else {
      labpos.indiv   <- cm.indiv[length(tau),"mean",]
      labpos.overall <- cm.overall[length(tau),"mean"]
    }
    axis(2)
    for (i in 1:x$k)
      axis(side=4, at=labpos.indiv[i],
           labels=x$labels[i], tick=FALSE,
           col.axis=labcol[i], las=1)
    axis(side=4, at=labpos.overall,
         labels=meanlabel, tick=FALSE,
         col.axis=meanlabcol, las=1)
    invisible()
  }
  
  taumarginal <- function(x)
  # NB: function is (essentially) identical to the one within "plot.bayesmeta()"
  {
    # range of tau values:
    tau <- seq(max(c(0,taurange[1]-0.04*diff(taurange))),
                   taurange[2]+0.04*diff(taurange), le=200)
    # corresponding posterior density:
    dens <- x$dposterior(tau=tau)
    # empty plot:
    maxdens <- max(dens[is.finite(dens)],na.rm=TRUE)
    plot(c(taurange[1],taurange[2]), c(0,maxdens), xlim=xlim,       
         type="n", axes=FALSE, xlab="", ylab="", main="")
    abline(v=vertlines, col=gridcol)
    # "fix" diverging density:
    dens[!is.finite(dens)] <- 10*maxdens
    # light grey shaded contour for density across whole range:
    polygon(c(0,tau,max(tau)), c(0,dens,0), border=NA, col=grey(0.90))
    # dark grey shaded contour for density within 95% bounds:
    indi <- ((tau>=x$summary["95% lower","tau"]) & (tau<=x$summary["95% upper","tau"]))
    polygon(c(rep(x$summary["95% lower","tau"],2), tau[indi], rep(x$summary["95% upper","tau"],2)),
            c(0, min(c(x$dposterior(tau=x$summary["95% lower","tau"]), 10*maxdens)),
              dens[indi], x$dposterior(tau=x$summary["95% upper","tau"]), 0),
            border=NA, col=grey(0.80))
    # vertical line at posterior median:
    lines(rep(x$summary["median","tau"],2), c(0,x$dposterior(tau=x$summary["median","tau"])), col=grey(0.6))
    # actual density line:
    lines(tau, dens, col="black")
    # y-axis:
    abline(v=0, col=grey(0.40))
    # x-axis:
    lines(taurange + c(-1,1) * 0.04*diff(taurange), c(0,0), col=grey(0.40))
    # plot prior density (if requested):
    if (prior) {
      lines(tau, x$dprior(tau=tau), col="black", lty="dashed")
    }
    # add axes, labels, bounding box, ...
    mtext(side=1, line=par("mgp")[1], expression("heterogeneity "*tau))
    #mtext(side=2, line=par("mgp")[2], expression("marginal posterior density"))
    if (infinity) {
      axis(1, at=c(vertlines, infx),
           labels=c(as.numeric(vertlines), expression(infinity)))
    } else {
      axis(1, at=vertlines)
    }
    invisible()
  }

  # make sure to properly re-set graphical parameters later:
  prevpar <- par(no.readonly=TRUE)
  on.exit(par(prevpar))
  # generate actual plot:
  graphics::layout(rbind(1,2), heights=c(2,1))
  par(mar=c(-0.1,3,0,rightmargin)+0.1, mgp=c(2.0, 0.8, 0))
  mutrace(x)
  par(mar=c(3,3,-0.1,rightmargin)+0.1)
  taumarginal(x)
  graphics::layout(1)
  par(mar=c(5,4,4,2)+0.1)
  invisible()
}

Try the bayesmeta package in your browser

Any scripts or data that you put into this service are public.

bayesmeta documentation built on July 9, 2023, 5:12 p.m.