From 48cba4975a7ed2b357ebd3e8b4378fb39592d28d Mon Sep 17 00:00:00 2001
From: Christian Roever <christian.roever@med.uni-goettingen.de>
Date: Fri, 20 Aug 2021 17:04:32 +0200
Subject: [PATCH] re-introduced 'traceplot.bmr()' function
---
NAMESPACE | 1 +
R/bmr.R | 126 +++++++++++++++++++++++++++++++++++++
man/traceplot.bayesmeta.Rd | 7 ++-
3 files changed, 132 insertions(+), 2 deletions(-)
diff --git a/NAMESPACE b/NAMESPACE
index 80b0b3a..9796f0d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -39,6 +39,7 @@ S3method(summary, bmr)
S3method(plot, bmr)
S3method(pairs, bmr)
S3method(forestplot, bmr)
+S3method(traceplot, bmr)
importFrom("grDevices",
"grey", "rainbow", "rgb")
diff --git a/R/bmr.R b/R/bmr.R
index 12f55cb..0577e7f 100644
--- a/R/bmr.R
+++ b/R/bmr.R
@@ -2097,3 +2097,129 @@ forestplot.bmr <- function(x,
"labeltext" = labeltext,
"forestplot" = fp))
}
+
+
+
+traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
+ rightmargin=8, col=rainbow(x$k), ...)
+{
+ stopifnot(missing(mulim) || (length(mulim) == 2),
+ missing(taulim) || (length(taulim) <= 2),
+ rightmargin >= 0, length(col) == x$k)
+ q975 <- qnorm(0.975)
+ gridcol <- "grey85"
+ # convert "taulim" and "mulim" 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)
+ }
+ vertlines <- pretty(taurange)
+ if (!missing(mulim) && (all(is.finite(mulim)) && (mulim[1] < mulim[2]))) {
+ murange <- mulim
+ } 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
+ }
+
+ 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))
+ for (i in 1:x$k) {
+ cm.indiv[,,i] <- x$shrink.moment(tau=tau, which=i)
+ }
+ plot(taurange, murange,
+ type="n", axes=FALSE, xlab="", ylab="effect", main="", ...)
+ abline(v=vertlines, col=gridcol)
+ abline(h=pretty(murange), col=gridcol)
+ abline(v=0, col=grey(0.40))
+ # grey 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)
+ }
+ }
+ # 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)
+ }
+ axis(2)
+ for (i in 1:x$k)
+ axis(side=4, at=cm.indiv[length(tau),"mean",i],
+ labels=x$labels[i], tick=FALSE,
+ col.axis=col[i], 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.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="")
+ 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")
+ # x-axis, y-axis:
+ abline(h=0, v=0, col=grey(0.40))
+ # 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()
+ }
+
+ # make sure to 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()
+}
diff --git a/man/traceplot.bayesmeta.Rd b/man/traceplot.bayesmeta.Rd
index 74ecb3f..99abbc5 100644
--- a/man/traceplot.bayesmeta.Rd
+++ b/man/traceplot.bayesmeta.Rd
@@ -2,6 +2,7 @@
\alias{traceplot}
\alias{traceplot.default}
\alias{traceplot.bayesmeta}
+\alias{traceplot.bmr}
\title{
Illustrate conditional means of overall effect as well as study-specific
estimates as a function of heterogeneity.
@@ -15,10 +16,12 @@
traceplot(x, ...)
\method{traceplot}{bayesmeta}(x, mulim, taulim, ci=FALSE,
rightmargin=8, col=rainbow(x$k), ...)
+ \method{traceplot}{bmr}(x, mulim, taulim, ci=FALSE,
+ rightmargin=8, col=rainbow(x$k), ...)
}
\arguments{
\item{x}{
- a \code{\link{bayesmeta}} object.
+ a \code{\link{bayesmeta}} or \code{\link{bmr}}object.
}
\item{mulim, taulim}{(optional) ranges for the effect (mu) and
heterogeneity (tau) axes. If only one value is given for
@@ -75,7 +78,7 @@
\doi{10.3102/10769986006004377}.
}
\seealso{
- \code{\link{bayesmeta}}.
+ \code{\link{bayesmeta}}, \code{\link{bmr}}.
}
\examples{
# load example data:
--
GitLab