oringexample.R

Jim — Nov 23, 2013, 9:42 AM

# Binomial Regression -- Outline of Approach

# Example:  O-Ring Data (The Challenger Disaster)
# Goal:  want to predict the probability of O-rank failure as a 
# function of temperature

# Prior beliefs:
# choose two temperatures, say 55 and 75, and assess prior on 
# corresponding probability of failure
# use theta1 is beta(1.6, 1), theta2 is beta(1, 1.6)

# Posterior:
# Display prior and posterior contours of (beta0, beta1)

# Predictive probabilities
# Display 90% probability intervals for range of temperatures.

library(DAAG)

Loading required package: lattice

# data setup

Failure <- ifelse(orings$Total > 0, 1, 0)
mx <- mean(orings$Temperature); sx = sd(orings$Temperature)
x <- (orings$Temperature - mx)/sx
data <- cbind(x, 1, Failure)

# graph the data

plot(data[, 1], jitter(data[, 3], factor=0.2), 
                xlab="Standardized Temperature",
                ylab="Failure",
                pch=19, cex=1.5, col="red")

# input prior information

a1 <- 1.6; b1 <- 1; n1 <- a1 + b1
a2 <- 1; b2 <- 1.6; n2 <- a2 + b2
x0 <- (c(55, 75) - mx) / sx
prior <- cbind(x0, c(n1, n2), c(a1, a2))

# the function that defines the logistic posterior for a simple
# logistic model is available on LearnBayes

require(LearnBayes)

Loading required package: LearnBayes

logisticpost

function (beta, data) 
{
    x = data[, 1]
    n = data[, 2]
    y = data[, 3]
    beta0 = beta[1]
    beta1 = beta[2]
    logf = function(x, n, y, beta0, beta1) {
        lp = beta0 + beta1 * x
        p = exp(lp)/(1 + exp(lp))
        y * log(p) + (n - y) * log(1 - p)
    }
    return(sum(logf(x, n, y, beta0, beta1)))
}
<environment: namespace:LearnBayes>

# plot prior and posterior

mycontour(logisticpost, c(-7, 7, -5, 4), prior,
          xlab="beta0", ylab="beta1")
mycontour(logisticpost, c(-7, 7, -5, 4), rbind(prior, data), 
          col="red", add=TRUE)
legend(4, 4, legend=c("Prior", "Posterior"),
       lty=1, col=c("black", "red"))

# summarizing the posterior

fit <- laplace(logisticpost, c(0,0), rbind(prior, data))
bfit <- rwmetrop(logisticpost, list(var = fit$var, scale = 2), 
                fit$mode, 10000, rbind(prior, data))
str(bfit)

List of 2
 $ par   : num [1:10000, 1:2] -0.927 -0.927 -0.927 -0.927 -0.95 ...
 $ accept: num 0.314

# want to plot predicted probs for temps between 32 and 92

pprob <- function(temp){
  x <- (temp - mx)/sx
  lp <- bfit$par[,1] + bfit$par[,2]*x
  p <- exp(lp) / (1 + exp(lp))
  quantile(p, c(.05, .5, .95))
}

pprob(52)

    5%    50%    95% 
0.5443 0.8751 0.9838 

temps <- 32:92
P <- sapply(temps, pprob)
plot(temps, P[2,], type = "l", lwd = 2, ylim = c(0, 1),
     xlab="Temperature", ylab="Prob(failure)",
     col="red",
     main="Predictive Probability of Failure")
lines(temps, P[1,], lty=2)
lines(temps, P[3,], lty=3)