From 3c1be560553de31d6b5e0ba2c78d87ca11315ece Mon Sep 17 00:00:00 2001
From: Christian Roever <christian.roever@med.uni-goettingen.de>
Date: Fri, 14 Apr 2023 12:43:04 +0200
Subject: [PATCH] extended 'traceplot()' functionality

---
 DESCRIPTION                |   2 +-
 NAMESPACE                  |  10 +-
 R/bayesmeta.R              | 131 +++++++++++++++++-------
 R/bmr.R                    | 203 +++++++++++++++++++++++++++++++------
 man/bayesmeta-package.Rd   |   2 +-
 man/traceplot.bayesmeta.Rd |  59 +++++++++--
 6 files changed, 324 insertions(+), 83 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 71c1aed..dd4b7eb 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: bayesmeta
 Type: Package
 Title: Bayesian Random-Effects Meta-Analysis and Meta-Regression
 Version: 3.3
-Date: 2023-03-16
+Date: 2023-04-14
 Authors@R: c(person(given="Christian", family="Roever", role=c("aut","cre"), 
                     email="christian.roever@med.uni-goettingen.de",
                     comment=c(ORCID="0000-0002-6911-698X")),
diff --git a/NAMESPACE b/NAMESPACE
index a51e4ca..80f3794 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -50,11 +50,11 @@ importFrom("graphics",
            "lines", "matlines", "mtext", "par", "plot", "points", "polygon",
            "text")
 importFrom("stats",
-           "dcauchy", "dlnorm", "dlogis", "dnorm", "dt", "integrate", "median",
-           "optim", "optimize", "pcauchy", "plnorm", "plogis", "pnorm",
-           "pt", "qcauchy", "qchisq", "qlnorm", "qlogis", "qnorm", "qt",
-           "rcauchy", "rexp", "rlogis", "rmultinom", "rnorm", "rt", "runif",
-           "uniroot")
+           "dcauchy", "dlnorm", "dlogis", "dnorm", "dt", "integrate", "lm.fit",
+           "median", "optim", "optimize", "pcauchy", "plnorm", "plogis",
+           "pnorm", "pt", "qcauchy", "qchisq", "qlnorm", "qlogis", "qnorm",
+           "qt", "rcauchy", "rexp", "rlogis", "rmultinom", "rnorm", "rt",
+           "runif", "uniroot")
 importFrom("forestplot",
            "forestplot", "fpDrawPointCI", "fpDrawSummaryCI","fpDrawBarCI",
            "fpColors")
diff --git a/R/bayesmeta.R b/R/bayesmeta.R
index dfbf1f7..f5c8a58 100644
--- a/R/bayesmeta.R
+++ b/R/bayesmeta.R
@@ -1,6 +1,6 @@
 #
 #    bayesmeta, an R package for Bayesian random-effects meta-analysis.
-#    Copyright (C) 2021  Christian Roever
+#    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
@@ -3629,14 +3629,30 @@ traceplot <- function(x, ...)
 
 traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
                                 ylab="effect",
-                                rightmargin=8, col=rainbow(x$k), ...)
+                                prior=FALSE,
+                                infinity=FALSE,
+                                rightmargin=8, col=rainbow(x$k),
+                                meanlabel="overall mean", meancol="black",
+                                ...)
 {
   stopifnot(missing(mulim) || (length(mulim) == 2),
             missing(taulim) || (length(taulim) <= 2),
-            rightmargin >= 0, length(col) == x$k)
+            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"
-  # convert "taulim" and "mulim" arguments
+  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]))
@@ -3648,34 +3664,55 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
   } else {
     taurange <- c(0, x$qposterior(tau=0.995)*1.1)
   }
-  
-  if (!missing(mulim) && (all(is.finite(mulim)) && (mulim[1] < mulim[2]))) {
-    murange <- mulim
+  if (infinity) {
+    xlim <- taurange + c(0, 0.15) * diff(taurange)
+    infx <- xlim[2] + 0.04*diff(xlim) # the "infinity" x-coordinate
   } else {
-    cm <- x$cond.moment(tau=taurange[2], indiv=TRUE)
-    if (ci){
-      murange <- range(c(cm[,"mean",]-q975*cm[,"sd",], cm[,"mean",]+q975*cm[,"sd",]))
-    } else {
-      murange <- range(cm[,"mean",])
-    }
-    murange <- murange + c(-1,1)*diff(murange)*0.05
+    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)
-  {
-    # range of tau values:
-    tau <- seq(max(c(0,taurange[1]-0.1*diff(taurange))),
-                   taurange[2]+0.1*diff(taurange), le=200)
-    cm.overall <- x$cond.moment(tau=tau)
+  {  
+    # 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)
-    plot(taurange, murange,         
+    # 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 shading:
+    # grey CI shading:
     if (ci) {
       for (i in 1:x$k) {
         polygon(c(tau, rev(tau)),
@@ -3695,18 +3732,33 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
       matlines(tau, cm.indiv[,"mean",]+q975*cm.indiv[,"sd",], col=col, lty=3)
     }
     # overall mean:
-    lines(tau, cm.overall[,"mean"], col="black", lty=2, lwd=1.5)
+    lines(tau, cm.overall[,"mean"], col=meancol, lty=2, lwd=1.5)
     if (ci) {
-      lines(tau, cm.overall[,"mean"]-q975*cm.overall[,"sd"], col="black", lty=3, lwd=1.5)
-      lines(tau, cm.overall[,"mean"]+q975*cm.overall[,"sd"], col="black", lty=3, lwd=1.5)
+      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=cm.indiv[length(tau),"mean",i],
+      axis(side=4, at=labpos.indiv[i],
            labels=x$labels[i], tick=FALSE,
            col.axis=col[i], las=1)
-    axis(side=4, at=cm.overall[length(tau),"mean"],
-         labels="overall mean", tick=FALSE, las=1)
+    axis(side=4, at=labpos.overall,
+         labels=meanlabel, tick=FALSE,
+         col.axis=meancol, las=1)
     invisible()
   }
   
@@ -3714,13 +3766,13 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
   # NB: function is (essentially) identical to the one within "plot.bayesmeta()"
   {
     # range of tau values:
-    tau <- seq(max(c(0,taurange[1]-0.1*diff(taurange))),
-                   taurange[2]+0.1*diff(taurange), le=200)
+    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),         
+    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:
@@ -3737,16 +3789,27 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
     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, y-axis:
-    abline(h=0, v=0, col=grey(0.40))
+    # 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"))
-    axis(1)#; box()
+    if (infinity) {
+      axis(1, at=c(vertlines, infx),
+           labels=c(as.numeric(vertlines), expression(infinity)))
+    } else {
+      axis(1, at=vertlines)
+    }
     invisible()
   }
 
-  # make sure to re-set graphical parameters later:
+  # make sure to properly re-set graphical parameters later:
   prevpar <- par(no.readonly=TRUE)
   on.exit(par(prevpar))
   # generate actual plot:
diff --git a/R/bmr.R b/R/bmr.R
index 515f403..f4e9507 100644
--- a/R/bmr.R
+++ b/R/bmr.R
@@ -1,6 +1,6 @@
 #
 #    bayesmeta, an R package for Bayesian random-effects meta-analysis.
-#    Copyright (C) 2021  Christian Roever
+#    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
@@ -2222,13 +2222,55 @@ forestplot.bmr <- function(x,
 
 traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
                           ylab="effect",
-                          rightmargin=8, col=rainbow(x$k), ...)
+                          prior=FALSE,
+                          infinity=FALSE,
+                          rightmargin=8, col=rainbow(x$k),
+                          X, Xlabels, Xcols="black",
+                          ...)
 {
   stopifnot(missing(mulim) || (length(mulim) == 2),
             missing(taulim) || (length(taulim) <= 2),
-            rightmargin >= 0, length(col) == x$k)
+            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)))
   q975 <- qnorm(0.975)
   gridcol <- "grey85"
+  if (length(col)==1) col <- rep(col, x$k)
+  if (infinity & any(x$beta.prior.proper)) {
+    warning("beta prior ignored for `tau=Inf' computations!")
+  }
+  if (prior & (!x$tau.prior.proper)) {
+    warning("Note that plots of improper priors may not be sensibly scaled.")
+  }
+  # default treatment of "X" argument:
+  if (missing(X) || (is.null(X) || all(is.na(X)))) {
+      X <- NULL
+      lincoNumber <- 0
+  } else {
+    stopifnot(is.numeric(X),
+              (is.matrix(X) && (ncol(X) == x$d))
+              || (is.vector(X) && (length(X) == x$d)))
+    if (is.vector(X)) {
+      X <- matrix(X, nrow=1, ncol=x$d)
+    }
+    if (is.null(colnames(X))) colnames(X) <- x$variables      
+    lincoNumber <- nrow(X)
+  }
+  # determine "X" labels and colours:
+  if (lincoNumber > 0) {
+    if (!missing(Xlabels)) {
+      stopifnot(is.character(Xlabels), length(Xlabels) == nrow(X))
+      rownames(X) <- Xlabels
+    } else {
+      Xlabels <- rownames(X)   # take from X row names...
+      if (is.null(Xlabels)) {  # ...or make up some
+        Xlabels <- sprintf("X.%02d", 1:nrow(X))
+      }
+    }
+    stopifnot(is.vector(Xcols), ((length(Xcols)==1) | (length(Xcols)==lincoNumber)))
+    if (length(Xcols)==1) Xcols <- rep(Xcols, lincoNumber)
+  }
   # convert "taulim" and "mulim" arguments
   # to eventual "taurange" and "murange" vectors:
   if (!missing(taulim) && all(is.finite(taulim))) {
@@ -2241,38 +2283,73 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
   } else {
     taurange <- c(0, x$qposterior(tau=0.995)*1.1)
   }
-  vertlines <- pretty(taurange)
-  if (!missing(mulim) && (all(is.finite(mulim)) && (mulim[1] < mulim[2]))) {
-    murange <- mulim
+  
+  if (infinity) {
+    xlim <- taurange + c(0, 0.15) * diff(taurange)
+    infx <- xlim[2] + 0.04*diff(xlim) # the "infinity" x-coordinate
   } else {
-    cm <- matrix(NA_real_, nrow=x$k, ncol=2,
-                 dimnames=list("study"=x$labels, "moments"=c("mean", "sd")))
-    for (i in 1:x$k)
-     cm[i,] <- x$shrink.moments(tau=taurange[2], which=i)
-    if (ci){
-      murange <- range(c(cm[,"mean"]-q975*cm[,"sd"], cm[,"mean"]+q975*cm[,"sd"]))
-    } else {
-      murange <- range(cm[,"mean"])
-    }
-    murange <- murange + c(-1,1)*diff(murange)*0.05
+    xlim <- taurange
+    infx <- NA_real_
   }
   
+  vertlines <- pretty(taurange)
+  # ensure no tickmarks beyond plotted tau range:
+  if (max(vertlines) > (taurange[2] + 0.04*diff(taurange)))
+    vertlines <- vertlines[-length(vertlines)]
+  
+  if (missing(mulim)) mulim <- NULL
+  
   mutrace <- function(x)
   {
-    # range of tau values:
-    tau <- seq(max(c(0,taurange[1]-0.1*diff(taurange))),
-                   taurange[2]+0.1*diff(taurange), le=200)
-    cm.indiv   <- array(NA_real_, dim=c(length(tau), 2, x$k),
-                        dimnames=list(NULL, c("mean", "sd"), x$labels))
+    # vector of tau values:
+    tau <- seq(max(c(0,taurange[1]-0.04*diff(taurange))),
+                   taurange[2]+0.04*diff(taurange), le=200)
+    # moments for individual studies:
+    cm.indiv <- array(NA_real_, dim=c(length(tau), 2, x$k),
+                      dimnames=list(NULL, c("mean", "sd"), x$labels))
     for (i in 1:x$k) {
       cm.indiv[,,i] <- x$shrink.moment(tau=tau, which=i)
     }
-    plot(taurange, murange,         
+    # moments for linear combinations (if applicable):
+    if (lincoNumber > 0) {
+      cm.linco <- array(NA_real_, dim=c(length(tau), 2, lincoNumber),
+                        dimnames=list(NULL, c("mean", "sd"), Xlabels))
+      for (i in 1:lincoNumber) {
+        cm.linco[,,i] <- x$pred.moments(tau=tau, x=X[i,])
+      }
+    }
+
+    # 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",])))
+        if (lincoNumber > 0) {
+          murange <- range(c(murange,
+                             range(cm.linco[,"mean",]-q975*cm.linco[,"sd",]),
+                             range(cm.linco[,"mean",]+q975*cm.linco[,"sd",])))
+        }
+      } else {
+        murange <- range(cm.indiv[,"mean",])
+        if (lincoNumber > 0) {
+          murange <- range(c(murange,
+                             range(cm.linco[,"mean",])))
+        }
+      }
+      # ensure that estimates are 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 shading:
+    # grey CI shading:
     if (ci) {
       for (i in 1:x$k) {
         polygon(c(tau, rev(tau)),
@@ -2280,6 +2357,14 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
                   rev(cm.indiv[,"mean",i] + q975*cm.indiv[,"sd",i])),
                 col=grey(0.75, alpha=0.25), border=NA)
       }
+      if (lincoNumber > 0) {
+        for (i in 1:lincoNumber) {
+          polygon(c(tau, rev(tau)),
+                  c(cm.linco[,"mean",i] - q975*cm.linco[,"sd",i],
+                    rev(cm.linco[,"mean",i] + q975*cm.linco[,"sd",i])),
+                  col=grey(0.75, alpha=0.25), border=NA)
+        }
+      }
     }    
     # individual estimates:
     matlines(tau, cm.indiv[,"mean",], col=col, lty=1)
@@ -2287,11 +2372,54 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
       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)
     }
+    # linear combinations:
+    if (lincoNumber > 0) {
+      matlines(tau, cm.linco[,"mean",], col=Xcols, lty=2, lwd=1.5)
+      if (ci) {
+        matlines(tau, cm.linco[,"mean",]-q975*cm.linco[,"sd",],
+                 col=Xcols, lty=3, lwd=1.5)
+        matlines(tau, cm.linco[,"mean",]+q975*cm.linco[,"sd",],
+                 col=Xcols, lty=3, lwd=1.5)
+      }
+    }
+      
+    if (infinity) {
+      # individual studies:
+      labpos.indiv   <- 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)
+      }
+      # linear combinations:
+      if (lincoNumber > 0) {
+        # compute (unweighted) least-squares fit:
+        lscoef <- stats::lm.fit(x=x$X, y=x$y)$coefficients
+        labpos.linco <- X %*% lscoef
+        for (i in 1:lincoNumber) {
+          lines(c(max(tau), infx),
+                c(cm.linco[length(tau),"mean",i], labpos.linco[i]),
+                col=Xcols[i], lty="13", lwd=2)
+        }
+      }
+    } else {
+      labpos.indiv <- cm.indiv[length(tau),"mean",]
+      if (lincoNumber > 0) {
+        labpos.linco <- cm.linco[length(tau),"mean",]
+      }
+    }
     axis(2)
     for (i in 1:x$k)
-      axis(side=4, at=cm.indiv[length(tau),"mean",i],
+      axis(side=4, at=labpos.indiv[i],
            labels=x$labels[i], tick=FALSE,
            col.axis=col[i], las=1)
+    if (lincoNumber > 0) {
+      for (i in 1:lincoNumber) {
+        axis(side=4, at=labpos.linco[i],
+             labels=Xlabels[i], tick=FALSE,
+             col.axis=Xcols[i], las=1)
+      }
+    }
     invisible()
   }
   
@@ -2299,13 +2427,13 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
   # NB: function is (essentially) identical to the one within "plot.bayesmeta()"
   {
     # range of tau values:
-    tau <- seq(max(c(0,taurange[1]-0.1*diff(taurange))),
-                   taurange[2]+0.1*diff(taurange), le=200)
+    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),         
+    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:
@@ -2322,16 +2450,27 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
     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, y-axis:
-    abline(h=0, v=0, col=grey(0.40))
+    # 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"))
-    axis(1)#; box()
+    if (infinity) {
+      axis(1, at=c(vertlines, infx),
+           labels=c(as.numeric(vertlines), expression(infinity)))
+    } else {
+      axis(1, at=vertlines)
+    }
     invisible()
   }
   
-  # make sure to re-set graphical parameters later:
+  # make sure to properly re-set graphical parameters later:
   prevpar <- par(no.readonly=TRUE)
   on.exit(par(prevpar))
   # generate actual plot:
@@ -2342,5 +2481,5 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
   taumarginal(x)
   graphics::layout(1)
   par(mar=c(5,4,4,2)+0.1)
-  invisible()
+  invisible(X)
 }
diff --git a/man/bayesmeta-package.Rd b/man/bayesmeta-package.Rd
index 975dd70..d9d3fe1 100644
--- a/man/bayesmeta-package.Rd
+++ b/man/bayesmeta-package.Rd
@@ -16,7 +16,7 @@
     Package: \tab bayesmeta\cr
     Type:    \tab Package\cr
     Version: \tab 3.3\cr
-    Date:    \tab 2023-03-16\cr
+    Date:    \tab 2023-04-14\cr
     License: \tab GPL (>=2)
   }
   The main functionality is provided by the \code{\link{bayesmeta}()}
diff --git a/man/traceplot.bayesmeta.Rd b/man/traceplot.bayesmeta.Rd
index 0bbb6fc..41d5a9f 100644
--- a/man/traceplot.bayesmeta.Rd
+++ b/man/traceplot.bayesmeta.Rd
@@ -4,22 +4,27 @@
 \alias{traceplot.bayesmeta}
 \alias{traceplot.bmr}
 \title{
-  Illustrate conditional means of overall effect as well as study-specific
-  estimates as a function of heterogeneity.
+  Illustrate conditional means of study-specific estimates as well as
+  overall mean (or other linear combinations) as a function of
+  heterogeneity.
 }
 \description{
-  Generates a trace plot of overall mean effect and study-specific
-  (shrinkage) estimates as a function of the heterogeneity (\eqn{\tau}). The
-  heterogeneity's posterior distribution is also indicated.
+  Generates a trace plot of study-specific (shrinkage) estimates as a
+  function of the heterogeneity (\eqn{\tau}), along with conditional
+  estimates of the overall mean or other linear combinations of
+  regression parameters.
+  The heterogeneity's posterior distribution is also shown at the bottom.
 }
 \usage{
   traceplot(x, ...)
   \method{traceplot}{bayesmeta}(x, mulim, taulim, ci=FALSE,
-          ylab="effect",
-          rightmargin=8, col=rainbow(x$k), ...)
+          ylab="effect", prior=FALSE, infinity=FALSE,
+          rightmargin=8, col=rainbow(x$k),
+          meanlabel="overall mean", meancol="black", ...)
   \method{traceplot}{bmr}(x, mulim, taulim, ci=FALSE,
-          ylab="effect",
-          rightmargin=8, col=rainbow(x$k), ...)
+          ylab="effect", prior=FALSE, infinity=FALSE,
+          rightmargin=8, col=rainbow(x$k),
+          X, Xlabels, Xcols="black", ...)
 }
 \arguments{
   \item{x}{
@@ -35,7 +40,16 @@
     confidence intervals.
   }
   \item{ylab}{
-    a y-axis label.
+    a label for the effect (mu) axis.
+  }
+  \item{prior}{
+    a logical flag indicating whether to show the (heterogeneity) prior
+    density along with its posterior.
+  }
+  \item{infinity}{
+    a logical flag indicating whether add an \dQuote{infinity} tickmark
+    to the heterogeneity (tau) axis and show the corresponding limiting
+    values.
   }
   \item{rightmargin}{
     an additional margin to be added to the right side of the plot, in
@@ -45,6 +59,25 @@
   \item{col}{
     colors to be used for plotting the (\eqn{k}) estimates.
   }
+  \item{meanlabel}{
+    a label for the overall mean estimate
+    (\code{traceplot.bayesmeta()}).
+  }
+  \item{meancol}{
+    colour specification for the overall mean estimate
+    (\code{traceplot.bayesmeta()}).
+  }
+  \item{X}{
+    matrix (or vector) of coefficients defining linear combinations of
+    regression parameters to be shown (\code{traceplot.bmr()}).
+  }
+  \item{Xlabels}{
+    labels for the linear combinations (\code{traceplot.bmr()}).
+  }
+  \item{Xcols}{
+    colour specification for the linear combinations
+    (\code{traceplot.bmr()}).
+  }
   \item{...}{
     other arguments passed on to the
     \code{\link[graphics]{plot}()} function.
@@ -87,6 +120,12 @@
   Estimation in parallel randomized experiments.
   \emph{Journal of Educational Statistics}, \bold{6}(4):377-401, 1981.
   \doi{10.3102/10769986006004377}.
+  
+  DuMouchel, W. H. (1994).
+  Hierarchical Bayes linear models for meta-analysis.
+  Technical Report 27, National Institute of Statistical Sciences (NISS);
+  Research Triangle Park, NC, USA.
+  \url{https://www.niss.org/research/technical-reports/hierarchical-bayes-linear-models-meta-analysis-1994}
 } 
 \seealso{
   \code{\link{bayesmeta}}, \code{\link{bmr}}.
-- 
GitLab