R/bugs.plot.inferences.R

In MultiBUGS/R2MultiBUGS: Running MultiBUGS from R

bugs.plot.inferences <- function(sims, display.parallel, ...){
  if (.Device == "windows" ||
        (.Device == "null device" && options("device") == "windows")){
    cex.names <- 0.7
    cex.axis <- 0.6
    cex.tiny <- 0.4
    cex.points <- 0.7
    standard.width <- 30
    max.width <- 40
    min.width <- 0.02
  } else {
    cex.names <- 0.7
    cex.axis <- 0.6
    cex.tiny <- 0.4
    cex.points <- 0.3
    standard.width <- 30
    max.width <- 40
    min.width <- 0.01
  }
  rootnames <- sims$root.short
  n.roots <- length(rootnames)
  sims.array <- sims$sims.array
  n.chains <- sims$n.chains
  dimension.short <- sims$dimension.short
  indexes.short <- sims$indexes.short
  long.short <- sims$long.short
  height <- 0.6
  par(mar = c(0, 0, 1, 0))

  plot(c(0, 1),
       c(-n.roots - 0.5, -0.4),
       ann = FALSE,
       bty = "n",
       xaxt = "n",
       yaxt = "n",
       type = "n")

  W <- max(strwidth(rootnames, cex = cex.names))
  B <- (1 - W)/3.8
  A <- 1 - 3.5 * B
  if (display.parallel){
    text(A, -0.4, "80% interval for each chain", adj = 0, cex = cex.names)
  } else {
    text(A, -0.4, "medians and 80% intervals", adj = 0, cex = cex.names)
  }
  num.height <- strheight(1:9, cex = cex.tiny)
  for (k in 1:n.roots){
    text(0, -k, rootnames[k], adj = 0, cex = cex.names)
    J <- min(length(long.short[[k]]), max.width)
    if (k == 1){
      index <- 1:J
    } else {
      index <- sum(unlist(lapply(long.short, length))[1:(k - 1)]) + 1:J
    }
    spacing <- 3.5/max(J, standard.width)
    med <- numeric(J)
    i80 <- matrix(, J, 2)
    med.chains <- matrix(, J, sims$n.chains)
    i80.chains <- array(NA, c(J, sims$n.chains, 2))
    for (j in 1:J){
      med[j] <- median(sims.array[, , index[j]])
      i80[j, ] <- quantile(sims.array[, , index[j]], c(0.1, 0.9))
      for (m in 1:n.chains){
        med.chains[j, m] <- quantile(sims.array[, m, index[j]], 0.5)
        i80.chains[j, m, ] <- quantile(sims.array[, m, index[j]], c(0.1, 0.9))
      }
    }
    rng <- range(i80, i80.chains)
    p.rng <- pretty(rng, n = 2)
    b <- height/(max(p.rng) - min(p.rng))
    a <- -(k + height/2) - b * p.rng[1]
    lines(A + c(0, 0), -k + c(-height/2, height/2))
    # plot a line at zero (if zero is in the range of the
    # mini-plot)
    if (min(p.rng) < 0 & max(p.rng) > 0){
      lines(A + B * spacing * c(0, J + 1), rep(a, 2), lwd = 0.5, col = "gray")
    }

    for (x in p.rng){
      text(A - B * 0.2, a + b * x, x, cex = cex.axis)
      lines(A + B * c(-0.05, 0), rep(a + b * x, 2))
    }
    for (j in 1:J){
      if (display.parallel){
        for (m in 1:n.chains){
          interval <- a + b * i80.chains[j, m, ]
          if (interval[2] - interval[1] < min.width){
            interval <- mean(interval) + c(-1, 1) * min.width/2
          }
          lines(A + B * spacing *
                  rep(j + 0.6 * (m - (n.chains + 1)/2)/n.chains, 2),
                interval,
                lwd = 0.5,
                col = m + 1)
        }
      } else {
        lines(A + B * spacing * rep(j, 2), a + b * i80[j,], lwd = 0.5)
        for (m in 1:n.chains){
          # points (A+B*spacing*j, a + b*med[j], pch=20,
          # cex=cex.points)
          points(A + B * spacing * j, a + b * med.chains[j, m],
                 pch = 20,
                 cex = cex.points,
                 col = m + 1)
        }
      }
      dk <- dimension.short[k]
      if (dk > 0){
        for (m in 1:dk){
          index0 <- indexes.short[[k]][[j]][m]
          if (j == 1){
            text(A + B * spacing * j,
                 -k - height/2 - 0.05 - num.height * (m - 1),
                 index0,
                 cex = cex.tiny)
          } else if (index0 != indexes.short[[k]][[j - 1]][m] &
                       (index0%%(floor(log10(index0) + 1)) == 0)){
            text(A + B * spacing * j,
                 -k - height/2 - 0.05 - num.height * (m - 1),
                 index0,
                 cex = cex.tiny)
          }
        }
      }
    }
    if (J < length(long.short[[k]])){
      text(-0.015, -k, "*", cex = cex.names, col = "red")
    }
  }
}