DESCRIPTION

 Package: bayesmeta
 Type: Package
 Title: Bayesian Random-Effects Meta-Analysis and Meta-Regression
-Version: 3.3
-Date: 2023-06-27
+Version: 3.35
+Date: 2023-11-29
 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")),

R/bayesmeta.R

@@ -453,8 +453,8 @@ bayesmeta.default <- function(y, sigma, labels=names(y),
     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))
+      samp <- matrix(NA_real_, nrow=n, ncol=2, dimnames=list(NULL,c("tau","mu")))
+      if (predict | 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))})
@@ -465,7 +465,7 @@ bayesmeta.default <- function(y, sigma, labels=names(y),
       }
       samp[,2] <- apply(matrix(samp[,"tau"],ncol=1), 1, cond.sample)
     } else {           # draw marginal, univariate (mu or theta) numbers:
-      samp <- rep(NA, n)
+      samp <- rep(NA_real_, n)
       if (!predict & (is.logical(individual) && (!individual)))
         meansd <- support[,c("mean","sd")]
       else
@@ -3646,7 +3646,12 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
             is.character(meanlabel), length(meanlabel)==1,
             length(meancol)==1)
   q975 <- qnorm(0.975)
-  gridcol <- "grey85"
+  # specify colors for axes etc.:
+  colvec <- c("axis"   = "grey40",
+              "grid"   = "grey85",
+              "median" = "grey60",
+              "ci"     = "grey80",
+              "tail"   = "grey90")
   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!")
@@ -3711,9 +3716,9 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
 
     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))
+    abline(v=vertlines, col=colvec["grid"])
+    abline(h=pretty(murange), col=colvec["grid"])
+    abline(v=0, col=colvec["axis"])
     # grey CI shading:
     if (ci) {
       for (i in 1:x$k) {
@@ -3776,25 +3781,25 @@ traceplot.bayesmeta <- function(x, mulim, taulim, ci=FALSE,
     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)
+    abline(v=vertlines, col=colvec["grid"])
     # "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))
+    polygon(c(0,tau,max(tau)), c(0,dens,0), border=NA, col=colvec["tail"])
     # 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))
+            border=NA, col=colvec["ci"])
     # 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))
+    lines(rep(x$summary["median","tau"],2), c(0,x$dposterior(tau=x$summary["median","tau"])), col=colvec["median"])
     # actual density line:
     lines(tau, dens, col="black")
     # y-axis:
-    abline(v=0, col=grey(0.40))
+    abline(v=0, col=colvec["axis"])
     # x-axis:
-    lines(taurange + c(-1,1) * 0.04*diff(taurange), c(0,0), col=grey(0.40))
+    lines(taurange + c(-1,1) * 0.04*diff(taurange), c(0,0), col=colvec["axis"])
     # plot prior density (if requested):
     if (prior) {
       lines(tau, x$dprior(tau=tau), col="black", lty="dashed")

R/bmr.R

@@ -45,7 +45,9 @@ kldiv <- function(mu1, mu2, sigma1, sigma2, symmetrized=FALSE)
     KL <- KL + 0.5 * (sum(diag(sigma1inv %*% sigma2))
                       + (t(mu2-mu1) %*% sigma1inv %*% (mu2-mu1)) - k)
   } else {
-    KL <- KL + 0.5 * (log(det(sigma2)) - log(det(sigma1)))
+    #KL <- KL + 0.5 * (log(det(sigma2)) - log(det(sigma1)))
+    KL <- KL + 0.5 * (determinant(sigma2, logarithm=TRUE)$modulus
+                      - determinant(sigma1, logarithm=TRUE)$modulus)
     # (NB: this term cancels out for the symmetrized variant)
   }
   return(as.vector(KL))
@@ -369,7 +371,9 @@ bmr.default <- function(y, sigma, labels = names(y),
         betaHat     <- Vbeta %*% t(rbind(X,Xp)) %*% sigmaTauInv %*% c(y,yp)
         residual    <- c(y,yp) - rbind(X,Xp) %*% betaHat
         logdens  <- (0.5 * ((d-k) * log(2*pi)
-                            -log(det(sigmaTau)) + log(det(Vbeta))
+                            #-log(det(sigmaTau)) + log(det(Vbeta))
+                            -determinant(sigmaTau, logarithm=TRUE)$modulus
+                            +determinant(Vbeta, logarithm=TRUE)$modulus
                             -(t(residual) %*% sigmaTauInv %*% residual))
                      + log(tau.prior(t)))
       }
@@ -2237,7 +2241,12 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
             is.logical(infinity), length(infinity)==1,
             rightmargin >= 0, ((length(col)==x$k) | (length(col)==1)))
   q975 <- qnorm(0.975)
-  gridcol <- "grey85"
+  # specify colors for axes etc.:
+  colvec <- c("axis"   = "grey40",
+              "grid"   = "grey85",
+              "median" = "grey60",
+              "ci"     = "grey80",
+              "tail"   = "grey90")
   if (length(col)==1) col <- rep(col, x$k)
   if (infinity & any(x$beta.prior.proper)) {
     warning("beta prior ignored for `tau=Inf' computations!")
@@ -2348,9 +2357,9 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
   
     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))
+    abline(v=vertlines, col=colvec["grid"])
+    abline(h=pretty(murange), col=colvec["grid"])
+    abline(v=0, col=colvec["axis"])
     # grey CI shading:
     if (ci) {
       for (i in 1:x$k) {
@@ -2437,25 +2446,26 @@ traceplot.bmr <- function(x, mulim, taulim, ci=FALSE,
     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)
+    abline(v=vertlines, col=colvec["grid"])
     # "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))
+    polygon(c(0,tau,max(tau)), c(0,dens,0), border=NA, col=colvec["tail"])
     # 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))
+            border=NA, col=colvec["ci"])
     # 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))
+    lines(rep(x$summary["median","tau"],2), c(0,x$dposterior(tau=x$summary["median","tau"])),
+          col=colvec["median"])
     # actual density line:
     lines(tau, dens, col="black")
     # y-axis:
-    abline(v=0, col=grey(0.40))
+    abline(v=0, col=colvec["axis"])
     # x-axis:
-    lines(taurange + c(-1,1) * 0.04*diff(taurange), c(0,0), col=grey(0.40))
+    lines(taurange + c(-1,1) * 0.04*diff(taurange), c(0,0), col=colvec["axis"])
     # plot prior density (if requested):
     if (prior) {
       lines(tau, x$dprior(tau=tau), col="black", lty="dashed")

data/Rubin1981.R

@@ -6,6 +6,7 @@
 #
 
 Rubin1981 <- data.frame("school"=c("A","B","C","D","E","F","G","H"),
-                         "effect"=c(28.39, 7.94, -2.75, 6.82, -0.64, 0.63, 18.01, 12.16),
-                         "stderr"=c(14.9, 10.2, 16.3, 11.0, 9.4, 11.4, 10.4, 17.6),
-                         stringsAsFactors=FALSE)
+                        "n"=as.integer(c(50, 79, 39, 91, 99, 72, 94, 35)),
+                        "effect"=c(28.39, 7.94, -2.75, 6.82, -0.64, 0.63, 18.01, 12.16),
+                        "stderr"=c(14.9, 10.2, 16.3, 11.0, 9.4, 11.4, 10.4, 17.6),
+                        stringsAsFactors=FALSE)

man/Rubin1981.Rd

@@ -7,6 +7,7 @@
 \format{The data frame contains the following columns:
   \tabular{lll}{
     \bold{school} \tab \code{character} \tab school identifier \cr
+    \bold{n}      \tab \code{integer}   \tab number of students \cr
     \bold{effect} \tab \code{numeric}   \tab effect estimate \cr
     \bold{stderr} \tab \code{numeric}   \tab associated standard error 
   }

man/bayesmeta-package.Rd

@@ -15,8 +15,8 @@
   \tabular{ll}{
     Package: \tab bayesmeta\cr
     Type:    \tab Package\cr
-    Version: \tab 3.3\cr
-    Date:    \tab 2023-06-27\cr
+    Version: \tab 3.35\cr
+    Date:    \tab 2023-11-29\cr
     License: \tab GPL (>=2)
   }
   The main functionality is provided by the \code{\link{bayesmeta}()}

man/traceplot.bayesmeta.Rd

@@ -117,6 +117,11 @@
   Christian Roever \email{christian.roever@med.uni-goettingen.de}
 }
 \references{
+  C. Roever, D. Rindskopf, T. Friede.
+  How trace plots help interpret meta-analysis results.
+  \emph{(submitted for publication)}, 2023.
+  \url{https://arxiv.org/abs/2306.17043}.
+  
   C. Roever.
   Bayesian random-effects meta-analysis using the bayesmeta R package.
   \emph{Journal of Statistical Software}, \bold{93}(6):1-51, 2020.