I’m a big fan of Andrew’s Gelman blog (http://andrewgelman.com/). I think that my statistical intuition is way much better after reading it. For example, there’s a post about different types of errors in NHST, not limited to the widely known Type I and Type II errors - http://andrewgelman.com/2004/12/29/type_1_type_2_t/. You should read this before continuing because the rest of this post will be based on it, and the article which is linked in that post (http://www.stat.columbia.edu/~gelman/research/published/francis8.pdf).

The code below performs the simulation described at the beginning of the third paragraph of the “Type S error rates for classical and Bayesian single and multiple comparison procedures” article.

simulation <- function(tau, sigma, n = 2000) {
  
  # simulate the data as described at the begining of the 3rd paragraph.
  result <- t(replicate(n, {
    dif <- rnorm(n = 1, sd = 2 * tau)
    yy  <- rnorm(1, dif, 2 * sigma)
    c(dif, yy)
  })) 
  
  bThreshold <-
    1.96 * sqrt(2) * sigma * sqrt(1 + (sigma ^ 2) / (tau ^ 2)) # Equation (6)
  clasThreshold <- 1.96 * sqrt(2) * sigma # Equation (5)
  
  yy <- result[, 2]
  
  probConfClaimBayes <- mean(yy > bThreshold | yy < -bThreshold)
  probConfClaimClas  <-
    mean(yy > clasThreshold | yy < -clasThreshold)
  
  bayesConfClaim <- result[yy > bThreshold | yy < -bThreshold,]
  bayesConfClaim <- sign(bayesConfClaim)
  sErrorCondClaimBayes <-
    mean(bayesConfClaim[, 1] != bayesConfClaim[, 2])
  
  clasConfClaim <-
    result[yy > clasThreshold | yy < -clasThreshold,]
  clasConfClaim <- sign(clasConfClaim)
  sErrorCondClaimClass <-
    mean(clasConfClaim[, 1] != clasConfClaim[, 2])
  
  list(
    data = result,
    bThreshold = bThreshold,
    clasThreshold = clasThreshold,
    probConfClaimBayes = probConfClaimBayes,
    probConfClaimClas = probConfClaimClas,
    sErrorCondClaimBayes = sErrorCondClaimBayes,
    sErrorCondClaimClass = sErrorCondClaimClass
  )
  
}

format_result <- function(result) {
  
  c(
    sprintf("Prob of claim - Bayesian: %.2f%%", result$probConfClaimBayes * 100),
    sprintf("Prob of claim - Classical: %.2f%%", result$probConfClaimClas * 100),
    sprintf("Prob of S-error cond claim - Bayesian: %.2f%%", result$sErrorCondClaimBayes * 100),
    sprintf("Prob of S-error cond claim - Classical: %.2f%%", result$sErrorCondClaimClass * 100))
}

plot_tresholds <- function(result, ...) {
  plot(result$data, xlab = "theta_j - theta_k", ylab = "y_j - y_k", ...)
  abline(h = c(-result$bThreshold, result$bThreshold), lty = 2, col = "green")
  abline(h = c(-result$clasThreshold, result$clasThreshold), lty = 2, col = "red")
  abline(v = 0, lty = 2)
  legend("topleft", c("Bayesian Treshold", "Classical Treshold"), lwd = 2, lty = 1, col = c("green", "red"))
  
  legend("bottomright", format_result(result))

}

Here you can see the replicated plots from the page 7. I merged the Bayesian and Classical versions into one plot, and I also added the probabilities of claiming that there are differences and the conditional probability of S-Error given the claim.

Note that I got different values than those from figure 3, but the experimental setting might cause it (I’m not using equation described in the test to approximate the result), or I have an error in the code (but I cannot find one…). However, the main point of the article is still valid;)

set.seed(123)
plot_tresholds(simulation(0.5,1,2000), main = "tau: 0.5, sigma: 1")

plot_tresholds(simulation(1,1,2000), main = "tau: 1, sigma: 1")

plot_tresholds(simulation(2,1,2000), main = "tau: 2, sigma: 1")