JAGS users: runjags has been updated

Image credit: http://www.potracksgalore.com/
accessories/hanging-chains/premierpotrackchains.cfm

The package runjags by Matt Denwood has been recently updated. The runjags package has many nice facilities for running JAGS, especially a seamless interface for running parallel chains on multi-core computers and thereby greatly reducing the time needed to achieve a desired MCMC sample size from JAGS. Give it a look!

Visualization of Metropolis with rejected proposals

In response to a previous post, Tinu Schneider has embellished the Metropolis visualization of Maxwell Joseph so it includes explicit plotting of rejected proposals (shown as grey segments). Here are the animated gif and R code. Many thanks to Tinu for sending this! (The gif, below, may take a few moments to load.)

 
 
#
# Animating the Metropolis Algorithm
# -------------------------------------------
# source:
# http://mbjoseph.github.io/blog/2013/09/08/metropolis/
 
# Prepare
----------------------------------------
require(sm) # needs 'rpanel'
dir.create("metropolis_ex")
 
 
# Setup
# --------------------------------------
# population level parameters
mu    <- 6
sigma <- 3
 
# data-values
n <- 50 
 
# metropolis
iter   <- 10000
chains <-     3
sigmaJump <- 0.2
 
# plotting
seq1 <- seq(  1, 200,  by=  5)
seq2 <- seq(200, 450,  by= 10)
seq3 <- seq(450, 750,  by= 50)
seq4 <- seq(750, iter, by=150)
 
xlims <- c(4, 8)
ylims <- c(1, 5)
zlims <- c(0, 0.7)
 
widthPNG  <- 1200
heightPNG <-  750
 
 
 
# Helpers
# --------------------------------------
# collect some data (e.g. a sample of heights)
gendata <- function(n, mu, sigma) {
    rnorm(n, mu, sigma)
}
 
# log-likelihood function
ll <- function(x, muhat, sigmahat){
    sum(dnorm(x, muhat, sigmahat, log=TRUE))
}
 
# prior densities
pmu <- function(mu){
    dnorm(mu, 0, 100, log=TRUE)
}
 
psigma <- function(sigma){
    dunif(sigma, 0, 10, log=TRUE)
}
 
# posterior density function (log scale)
post <- function(x, mu, sigma){
    ll(x, mu, sigma) + pmu(mu) + psigma(sigma)
}
 
geninits <- function(){
    list(mu    = runif(1, 4, 10),
         sigma = runif(1, 2,  6))
}
 
jump <- function(x, dist = 0.1){ # must be symmetric
    x + rnorm(1, 0, dist)
}
 
 
# compute
# ----------------------------------------
 
 
# init 
proposed  <- array(dim = c(chains, 2, iter))
posterior <- array(dim = c(chains, 2, iter))
accepted  <- array(dim = c(chains, iter-1))
 
# the data
x <- gendata(n, mu, sigma)
 
for (c in 1:chains){
    props      <- array(dim=c(2, iter))
    theta.post <- array(dim=c(2, iter))
    inits <- geninits()
    theta.post[1, 1] <- inits$mu
    theta.post[2, 1] <- inits$sigma
    for (t in 2:iter){
        # theta_star = proposed next values for parameters
        jumpMu    <- jump(theta.post[1, t-1], sigmaJump)
        jumpSigma <- jump(theta.post[2, t-1], sigmaJump)
        theta_star <- c(jumpMu, jumpSigma)
        pstar <- post(x, mu = theta_star[1],      sigma = theta_star[2])
        pprev <- post(x, mu = theta.post[1, t-1], sigma = theta.post[2, t-1])
        lr <- pstar - pprev
        r  <- exp(lr)
 
        # theta_star is accepted if posterior density is higher w/ theta_star
        # if posterior density is not higher, it is accepted with probability r
        # else theta does not change from time t-1 to t
        accept <- rbinom(1, 1, prob = min(r, 1))
        accepted[c, t - 1] <- accept
        if (accept == 1){
            theta.post[, t] <- theta_star
        } else {
            theta.post[, t] <- theta.post[, t-1]
        }
        props[, t] <- theta_star
    }
    proposed[c, , ]  <- props
    posterior[c, , ] <- theta.post
}
 
 
# Plot
# ----------------------------------------------
 
oldPath <- getwd()
setwd("metropolis_ex")
 
sequence <- c(seq1, seq2, seq3, seq4)
cols <- c("blue", "purple", "red")
 
png(file = "metrop%03d.png", width=widthPNG, height=heightPNG)
    for (i in sequence){
        par(font.main=1, cex.main=1.2)
        layout(matrix(1:2, nrow=1), widths=c(1.0, 1.5))
  # prepare empty plot
        plot(0, type="n", 
             xlim=xlims, ylim=ylims,
             xlab="mu",  ylab="sigma", main="Markov chains")
        for (j in 1:chains) { # first all the gray lines
            lines(proposed[j, 1, 1:i],  proposed[j, 2, 1:i],  col="gray")
        }
        for (j in 1:chains) { # now the accepteds
            lines(posterior[j, 1, 1:i], posterior[j, 2, 1:i], col=cols[j])
        }
        text(x=6, y=1.0, paste("Iteration ", i), cex=1.2)
 
        sm.density(x=cbind(c(posterior[, 1, 1:i]), c(posterior[, 2, 1:i])),
                   xlab="mu",  ylab="sigma", zlab="", 
                   xlim=xlims, ylim=ylims,   zlim=zlims, 
                   col="white")
        title(main="Posterior density")
    }
dev.off()
 
 
# combine .png's into one .gif
# system('convert2gif.bat') # had to write a .bat-file
# system('convert -delay 20 *.png "metropolis.gif"') # doesent work
 
#file.remove(list.files(pattern=".png"))
 
setwd(oldPath)
 
#

Created by Pretty R at inside-R.org

Cool visualization of Metropolis sampling

A cool visualization of Metropolis sampling has been posted by Maxwell Joseph here. I've pasted in his animated gif below, as a teaser for you to check his post (and I figured it would be okay because his work focuses on "parasitic biodiversity" ;-). Thanks to Joseph Young for pointing me to this.

From Maxwell Joseph's post.

What would also be cool is an animation of Metropolis that shows the proposals at every step, even if rejected.

Bayesian methods article wins best paper award in Organizational Research Methods

The article, "The Time Has Come: Bayesian Methods for Data Analysis in the Organizational Sciences," won the 2012 Best Paper Award for the journal, Organizational Research Methods. The publisher has made the article freely available. You can see the announcement and find a link to the article here. The accompanying programs can be found here . If you find that article interesting, you would probably also like the article, "Bayesian estimation supersedes the t test," described here.

How much of a Bayesian posterior distribution falls inside a region of practical equivalence (ROPE)

The posterior distribution of a parameter shows explicitly the relative credibility of the parameter values, given the data. But sometimes people want to make a yes/no decision about whether a particular parameter value is credible, such as the "null" values of 0.0 difference or 0.50 chance probability. For making decisions about null values, I advocate using a region of practical equivalence (ROPE) along with a posterior highest density interval (HDI). The null value is declared to be rejected if the (95%, say) HDI falls completely outside the ROPE, and the null value is declared to be accepted (for practical purposes) if the 95% HDI falls completely inside the ROPE. This decision rule accepts the null value only when the posterior estimate is precise enough to fall within a ROPE. The decision rule rejects the null only when the posterior exceeds the buffer provided by the ROPE, which protects against hyper-inflated false alarms in sequential testing. And the decision rule is intuitive: You can graph the posterior, its HDI, its ROPE, and literally see what it all means in terms of the meaningful parameter being estimated. (For more details about the decision rule, its predecessors in the literature, and alternative approaches, see the article on this linked web page.)

But how do we set the limits of the ROPE? How big is "practically equivalent to the null value"? There is typically no uniquely correct answer to this question, because it depends on the particular domain of application and the current practice in that domain. But the rule does not need a uniquely correct ROPE, it merely needs a reasonable ROPE that is justifiable to the audience of the analysis. As a field matures, the limits of practical equivalence may change, and therefore what may be most useful to an audience is a description of how much of the posterior falls inside or outside the ROPE as a function of the width of the ROPE. The purpose of this blog post is to provide some examples and an R program for doing that.

Consider a situation in which we want to estimate an effect size parameter, which I'll call δ. Suppose we collect a lot of data so that we have a very precise estimate, and the posterior distribution shows that the effect size is very nearly zero, as plotted in the left panel below:

We see that zero is among the 95% most credible values of the effect size, but we can ask whether we should decide to "accept" the value zero (for practical purposes). If we establish a ROPE around zero, from -0.1 to +0.1 as shown above, then we see that the 95% HDI falls entirely within the  ROPE and we accept zero. What if we think a different ROPE is appropriate, either now or in the future? Simply display how much of the posterior distribution falls inside the ROPE, as a function of the ROPE size, as shown in the right panel above. The curve plots how much  of the  posterior distribution falls inside a ROPE centered at the null value, as a function of the radius (i.e., half width) of the ROPE. As a landmark, the plot shows dashed lines at the 95% HDI limit farthest from the null value. From the plot, readers can decide for themselves.

Here is another example. Suppose we are spinning a coin (or doing something more interesting that has dichotomous outcomes) and we want to estimate its underlying probability of heads, denoted θ. Suppose we collect a lot of data so that we have a precise estimate, as shown in the left panel below.

We see that the chance value of 0.50 is among the 95% most credible values of the underlying probability of heads, but we can ask whether we should decide to "accept" the value 0.50. If we establish a ROPE around zero, from 0.49 to 0.51 as shown above, we can see that the 95% HDI falls entirely within the ROPE and we would decide to accept 0.50. But we can provide more information by displaying the amount of the posterior distribution inside a ROPE centered on the null value as a function of the width of the ROPE. The right panel, above, allows the readers to decide for themselves.

The plots can be used for cases of rejecting the null, too. Suppose we spin a coin and the estimate shows a value apparently not at chance, as shown in the left panel below.

Should we reject the null value? The 95% HDI falls entirely outside the ROPE from 0.49 to 0.51, which means that the 95% most credible values are all not practically equivalent to 0.5. But what if we used a different ROPE? The posterior distribution, along with the explicit marking of the 95% HDI limits, is all we really need to see what the biggest ROPE could be and still let us reject the null. If we want more detailed information, the panel on the right, above, reveals the answer --- although we need to mentally subtract from 1.0 to get the posterior area not in the ROPE (on either side).

Here is the R script I used for creating the graphs above. Notice that it defines a function for plotting the area in the ROPE.

#------------------------------------------------------------------------------

# From http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/Programs/ get
# openGraphSaveGraph.R, plotPost.R, and HDIofMCMC.R. Put them in the working
# directory of this script, or specify the path in the next lines' source()
# commands:
source("./openGraphSaveGraph.R")
source("./plotPost.R")
source("./HDIofMCMC.R")

# Define function for computing area in ROPE and plotting it:
plotAreaInROPE = function( mcmcChain , maxROPEradius , compVal=0.0 ,
                           HDImass=0.95 , ... ) {
  ropeRadVec = seq( 0 , maxROPEradius , length=201 ) # arbitrary comb
  areaInRope = rep( NA , length(ropeRadVec) )
  for ( rIdx in 1:length(ropeRadVec) ) {
    areaInRope[rIdx] = ( sum( mcmcChain > (compVal-ropeRadVec[rIdx])
                            & mcmcChain < (compVal+ropeRadVec[rIdx]) )
                         / length(mcmcChain) )
  }
  plot( ropeRadVec , areaInRope ,
        xlab=bquote("Radius of ROPE around "*.(compVal)) ,
        ylab="Posterior in ROPE" ,
        type="l" , lwd=4 , col="darkred" , cex.lab=1.5 , ... )
  # From http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/Programs/
  # get HDIofMCMC.R. Put it in the working directory of this script, or specify
  # the path in the next line's source() command:
  # source("./HDIofMCMC.R")
  HDIlim = HDIofMCMC( mcmcChain , credMass=HDImass )
  farHDIlim = HDIlim[which.max(abs(HDIlim-compVal))]
  ropeRadHDI = abs(compVal-farHDIlim)
  areaInFarHDIlim = ( sum( mcmcChain > (compVal-ropeRadHDI)
                         & mcmcChain < (compVal+ropeRadHDI) )
                      / length(mcmcChain) )
  lines( c(ropeRadHDI,ropeRadHDI) , c(-0.5,areaInFarHDIlim) ,
         lty="dashed" , col="darkred" )
  text( ropeRadHDI , 0 ,
        bquote( atop( .(100*HDImass)*"% HDI limit" ,
                     "farthest from "*.(compVal) ) ) , adj=c(0.5,0) )
  lines( c(-0.5,ropeRadHDI) ,c(areaInFarHDIlim,areaInFarHDIlim) ,
         lty="dashed" , col="darkred" )
  text( 0 , areaInFarHDIlim , bquote(.(signif(areaInFarHDIlim,3))) ,
        adj=c(0,1.1) )
  return( list( ROPEradius=ropeRadVec , areaInROPE=areaInRope ) )
}

#------------------------------------------------------------------------------

# Generate a fictitious MCMC posterior:
m = 0.03
s = (0.1-m)/3
paramMCMCsample = rnorm(50000,m,s)

# Specify the desired comparison value, ROPE radius, and HDI mass:
compVal = 0.0
ropeRad = 0.1
HDImass = .95

# Open graphics window and specify layout:
openGraph(width=7.5,height=3.5)
par(mar=c(3.5,3.5,2,1),mgp=c(2,0.7,0))
layout( matrix( c(1,2) , nrow=1 , byrow=FALSE ) )

# Plot the posterior with HDI and ROPE:
postInfo = plotPost( paramMCMCsample , compVal=compVal ,
                     ROPE=compVal+c(-ropeRad,ropeRad) , showMode=TRUE ,
                     credMass=HDImass , xlab=expression(delta*" (parameter)") )

# Plot the area in ROPE:
ropeInfo = plotAreaInROPE( mcmcChain=paramMCMCsample , maxROPEradius=0.15 ,
                           compVal=compVal , HDImass=HDImass )

#------------------------------------------------------------------------------

Near Paris? Snap a shot of the book with Laplace. Or pose it with other Bayesians.

The book previously has been posed at Bayes' tomb, Fisher's remains, and Jacob Bernoulli's grave (in an attempt to be amusing and informative; hopefully not inadvertently being offensive). Conspicuously absent from that list is Pierre-Simon Lapace, who was the foundational developer of Bayesian methods. If you are in or near Paris, and would have fun posing the book at Laplace's grave, please do so and send me a photo! If you're among the first responders, I'll post the photo, with attribution to you, on the blog. Information about the location of Laplace's grave can be found here, with a map here.

If you pose the book with other famous Bayesians, or pre-Bayesians, or anti-Bayesians, either dead or not-quite-yet dead, please send me those photos too! (Again, the goal is to be amusing and informative, not offensive.) Thanks, and have fun!

Decisions from posterior distributions: Tail probability or highest density interval?

How should we decide whether a parameter's posterior distribution "rejects" a particular value such as zero? Should we consider the percentage of the distribution above/below the value? Should we consider the relation of the highest density interval (HDI) to the value? Here are some examples to explain why I think it makes more sense to use the HDI.

Here are the two decision rules being compared. First, the tail-probability decision rule: If there is less than 2.5% of the distribution on either side of the value, then reject the value. This is tantamount to using a 95% equal-tailed credibility interval: Values outside the 95% equal-tailed credibility interval are "rejected." Second, the HDI decision rule: Values outside the 95% HDI are "rejected." (Of course, I like to enhance the decision rule with a ROPE, to allow acceptance decisions and to provide a buffer against false alarms -- but that's a separate discusssion.)

The two histograms below represent the MCMC results for two hypothetical parameters.

The upper panel shows a posterior distribution for which the parameter value of zero falls will within the 95% HDI, but 2.0% (<2.5%) of the distribution falls below the value zero. Do we "reject" zero or not? I think it would be wrong to reject zero, because the probability density (i.e., credibility) of zero is high. Zero is clearly among the most credible values of the parameter, even though less than 2.5% of the distribution falls below zero.

The lower panel shows a posterior distribution for which the parameter value of zero falls well outside the 95% HDI (thus, even with of a modest ROPE, e.g. from -1 to +1, zero would still fall outside the 95% HDI), but a full 3.0% (>2.5%) of the distribution falls below zero. Do we "reject" zero or not? If we use a tail-probability decision rule, we do not reject zero. But clearly zero is not among the most credible values of the parameter, in that zero has low probability density.

Proponents of using equal-tailed credibility intervals will argue that percentiles of distributions are invariant under transformations of the parameter, but HDIs are not. True enough, but I think that most parameters are specifically scaled to be meaningful, and we want to know about credibility (probability density) on the meaningful scale, not on a meaningless transformed scale. But I am not saying that one decision rule is "correct" and the other is "wrong." The decision rules are merely rules with differing interpretations and characteristics; I am showing examples that convince me that a more useful, intuitive rule is using the HDI not the equal-tailed interval.

Then why do the plots in DBDA (like those above) bother to display the percentage of the distribution below/above the comparison value? Primarily merely as an additional descriptive statistic, but also to inform those who wish to think about tail probabilities for a decision rule that I eschew.

For another example, see this previous post.