added 'traceplot.bmr()' function

R/bmr.R

@@ -1753,3 +1753,134 @@ pairs.bmr <- function(x, ...)
   }
   invisible()
 }
+
+
+
+traceplot <- function(x, ...)
+{
+  UseMethod("traceplot")
+}
+
+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()
+}