In January 1986, the space shuttle Challenger exploded shortly after launch. An investigation was launched into the cause of the crash and attention focused on the rubber 0-ring seals in the rocket boosters. At lower temperatures, rubber becomes more brittle and is a less effective sealant. At the time of the launch, the temperature was 31F. Could the failure of the 0-rings have been predicted? In the 23 previous shuttle missions for which data exists, some evidence of damage due to blow by and erosion was recorded on some 0-rings. Each shuttle had two boosters, each with three 0-rings. For each mission, we know the number of 0-rings out of six showing some damage and the launch temperature.

Model 1: \(P(t)= \Phi{(\beta_0+\beta_1t, 1)}\)
Model 2: \(P(t)= \Phi{(\beta_0+\beta_1t+\beta_2t^2, 1)}\)

1. Theory
2. Model 1: Beta
3. Model 1: Latent Variables
4. Model 2: Beta
5. Model 2: Latent Variables
6. Plots
7. Conclusion

• Let \(Y_1,\ldots,Y_N\) be independent Binary random variables, with \(Y_i|p_i\ \sim Bernoulli(p_i)\) and \(p_i\ =\ F({\underline{x}}_i^\prime,\ \underline{\beta}), i=1,\ldots,N\), where F is a known distribution function.

• Let us consider \(F=\Phi\)=Standard normal c.d.f. That is, \(F(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{\frac{-u^2}{2}}du\) and
  \(f(x)=\frac{d}{dx}(F(x))=\frac{e^{\frac{-x^2}{2}}}{\sqrt{2\pi}}\) = standard normal density function.

• \(\underline{\beta_{p*1}}\ \epsilon\ \mathbb{R}^p\) is a vector of unknown parapters.

• Since \(Y_1,\ldots,Y_N\) are independent Bernoulli randonm variables, the joint probability function of \(Y_1,\ldots,Y_N|\underline{\beta}\) is given by \(f(\underline{y}|\underline{\beta})=\prod^{N}_{i=1} {F({x}_i^\prime \beta)^{y_{i}}}(1-F({x}_i^\prime \beta)^{1-y_{i}}\)

• Let \(\tau(\underline{\beta})\) be a prior density on \(\underline{\beta}\)
  

• Then the posterior density of \(\underline{\beta}\) is given by \(\tau(\underline{\beta}|\underline{y})=\frac{\prod^{N}_{i=1}f(\underline{y}|\underline{\beta})\tau(\beta)}{\int_{\mathbb{R}^p}{\prod^{N}_{i=1}f(\underline{y}|\underline{\upsilon })\tau(u)}du}\) where \(F=\Phi\)

• \(Z_1,\ldots,Z_N \sim N(\underline{{x}_i^\prime}\underline{\beta},1)\)

• \(Z_i\) has distribution function \(\Phi(\underline{{x}_i^\prime}\underline{\beta})\)

Fully Conditional Posterior Density: \[\tau(\underline{\beta}|\underline{Z},\underline{y})\alpha \prod^{N}_{i=1}{[I(Z_{i}>0,y_{i}=1)+I(Z_{i}\leq0,y_{i}=0)}]*e^{\frac{-(Z_{i}-\underline{{x}_i^\prime}\underline{\beta})^2}{2}}*\tau(\beta)\]

Let \(\hat{\beta_{z}}\) be equal to \((x' x)x' z\) and an unbiased estimator for \(\beta\)

\((z-x\beta)'(z-x\beta)\)=

\([(z-x\hat{\beta_{z}})+(x\hat{\beta_{z}}-x\beta)]'*[(z-x\hat{\beta_{z}})+(x\hat{\beta_{z}}-x\beta)]\)=

\((z-x\hat{\beta_{z}})'(z-x\hat{\beta_{z}})\)+ \((z-x\hat{\beta_{z}})'(x\hat{\beta_{z}}-x\beta)\)+
\((x\hat{\beta_{z}}-x\beta)'(z-x\hat{\beta_{z}})\)+ \((x\hat{\beta_{z}}-x\beta)'(x\hat{\beta_{z}}-x\beta)\)

Because \(\hat{\beta_{z}}\) is unbiased, terms that include \((x\hat{\beta_{z}}-x\beta)\) will equal zero
\((z-x\beta)'(z-x\beta)\)= \((z-x\hat{\beta_{z}})'(z-x\hat{\beta_{z}})\)+ \((\beta-\hat{\beta_{z}})x' x(\beta-\hat{\beta_{z}})\)

Because the first term doesn’t depend on \(\beta\), we can drop it.

Model 1: \(P(t)= \Phi{(\beta_0+\beta_1t, 1)}\)
Model 2: \(P(t)= \Phi{(\beta_0+\beta_1t+\beta_2t^2, 1)}\)

Model 1: \(\underline{\beta} \sim \mathcal{N}(\underline{\widehat\beta_z},\,(X'X)^{-1})\) with \(\underline{\beta}=\begin{bmatrix}\beta_0\\\beta_1\end{bmatrix}\), \(\widehat{\underline{\beta_z}}=(X'X)^{-1}X'\underline{Z}\)

\(\underline{Z_i} \sim \mathcal{N}(\underline{x_i'} \underline{\beta},\,1)\)

Model 2: \(\underline{\beta} \sim \mathcal{N}(\underline{\widehat\beta_z},\,(X'X)^{-1})\) with \(\underline{\beta}=\begin{bmatrix}\beta_0\\\beta_1\\\beta_2\end{bmatrix}\), \(\widehat{\underline{\beta_z}}=(X'X)^{-1}X'\underline{Z}\)

\(\underline{Z} \sim \mathcal{N}(\underline{x_i'} \underline{\beta},\,1)\)

y <- dat$affected # Binary list of damaged/not danaged
n_obs <- length(y)
# Model 1 X Matrix
X <- cbind(rep(1, n_obs), dat$temperature)

burn_in <- 50000
n_iter <- 90000

XtX <- t(X) %*% X
XtXi <- solve(XtX)

n_params <- 2
z0 <- rnorm(n=n_obs, mean=0, sd=1) # Random  init z
beta_hat <- solve(XtX) %*% t(X) %*% z0

# Matrices to use later
beta_out <- matrix(data = NA, 
                   nrow = n_iter, ncol = n_params)

zed_out <- matrix(data = NA, 
                  nrow = n_iter, ncol = 138)

meanz_out <- matrix(data= NA, nrow= n_iter, ncol= 138)

zed <- vector(length=138)

# Values needed for rtruncnorm function
ones <- matrix(data = 1, nrow = 1, ncol = 138)
num_ones <- sum(dat$affected)
num_zeroes <- length(dat$affected[dat$affected == 0])

for (i in (1:n_iter)) {
  beta <- mvrnorm(n = 1, mu = beta_hat, XtXi)
  zed[y == 0] <- rtruncnorm(n=num_zeroes, 
                            mean = X[y ==0,] %*% beta, 
                            sd = 1, a = -Inf, b = 0)
  
  zed[y == 1] <- rtruncnorm(n=num_ones, 
                            mean = X[y == 1,] %*% beta, 
                            sd = 1, a = 0, b = Inf)
  # Update our beta_hat based off previous iteration
  beta_hat <- solve(XtX) %*% t(X) %*% zed
  # Outputs
  beta_out[i, ] <- beta
  zed_out[i, ] <- zed
  meanz_out[i,] <- X %*% beta
}

n_samples <- 2000
mean_beta1=0
mean_beta2=0
for (h in (1:n_samples)) {
  mean_beta1 <- mean_beta1 + beta_out[burn_in+(20*h),1]
  mean_beta2 <- mean_beta2 + beta_out[burn_in+(20*h),2]
}

mean_b1 <- mean_beta1/n_samples
mean_b2 <- mean_beta2/n_samples

# Matrices to use later
beta_out <- matrix(data = NA, 
                   nrow = n_iter, ncol = n_params)

zed_out <- matrix(data = NA, 
                  nrow = n_iter, ncol = 138)

meanz_out <- matrix(data= NA, nrow= n_iter, ncol= 138)

zed <- vector(length=138)

# Values needed for rtruncnorm function
ones <- matrix(data = 1, nrow = 1, ncol = 138)
num_ones <- sum(dat$affected)
num_zeroes <- length(dat$affected[dat$affected == 0])

for (i in (1:n_iter)) {
  beta <- mvrnorm(n = 1, mu = beta_hat, XtXi)
  zed[y == 0] <- rtruncnorm(n=num_zeroes, 
                            mean = X[y ==0,] %*% beta, 
                            sd = 1, a = -Inf, b = 0)
  
  zed[y == 1] <- rtruncnorm(n=num_ones, 
                            mean = X[y == 1,] %*% beta, 
                            sd = 1, a = 0, b = Inf)
  # Update our beta_hat based off previous iteration
  beta_hat <- solve(XtX) %*% t(X) %*% zed
  # Outputs
  beta_out[i, ] <- beta
  zed_out[i, ] <- zed
  meanz_out[i,] <- X %*% beta
}

mean_z <- rep(0,138)
for(h in (1:n_samples)){
  for(i in (1:138)) {
    mean_z[i] <- mean_z[i] + 
      meanz_out[burn_in+(20*h),i]
  }
}
mean_z <- mean_z/n_samples

## Warning in kable_pipe(x = structure(c("-0.0032383", "-0.4398175",
## "-0.5489622", : The table should have a header (column names)

y <- dat$affected # Binary list of damaged/not damaged
n_obs <- length(y)

# Model 2 X Matrix
X <- cbind(rep(1, n_obs), dat$temperature, 
           (dat$temperature)^2)

XtX <- t(X) %*% X
XtXi <- solve(XtX)

n_params <- 3

z0 <- rnorm(n=n_obs, mean=0, sd=1) # Random init z
beta_hat <- solve(XtX) %*% t(X) %*% z0

# Matrices to use later
beta_out <- matrix(data = NA, nrow = n_iter, 
                   ncol = n_params)

zed_out <- matrix(data = NA, nrow = n_iter,
                  ncol = 138)

meanz_out <- matrix(data = NA, nrow = n_iter, 
                    ncol = 138)
zed <- vector(length=138)

# Values needed for rtruncnorm function
ones <- matrix(data = 1, nrow = 1, ncol = 138)
num_ones <- sum(dat$affected)
num_zeroes <- length(dat$affected[dat$affected == 0])

for (i in (1:n_iter)) {
  beta <- mvrnorm(n = 1, mu = beta_hat, XtXi)
  zed[y == 0] <- rtruncnorm(n=num_zeroes, 
                            mean = X[y ==0,] %*% beta,
                            sd = 1, a = -Inf, b = 0)
  
  zed[y == 1] <- rtruncnorm(n=num_ones, 
                            mean = X[y == 1,] %*% beta, 
                            sd = 1, a = 0, b = Inf)
  # Update out beta_hat based off previous iteration
  beta_hat <- solve(XtX) %*% t(X) %*% zed
  # Outputs
  beta_out[i, ] <- beta
  zed_out[i, ] <- zed
  meanz_out[i,] <- X %*% beta
}

n_samples <- 2000
mean_beta1 <- 0
mean_beta2 <- 0
mean_beta3 <- 0

for (h in (1:n_samples)) {
  mean_beta1 <- mean_beta1 + beta_out[burn_in+(20*h),1]
  mean_beta2 <- mean_beta2 + beta_out[burn_in+(20*h),2]
  mean_beta3 <- mean_beta3 + beta_out[burn_in+(20*h),3]
}

mean_b1 <- mean_beta1/n_samples 
mean_b2 <- mean_beta2/n_samples 
mean_b3 <- mean_beta3/n_samples

# Matrices to use later
beta_out <- matrix(data = NA, nrow = n_iter, 
                   ncol = n_params)

zed_out <- matrix(data = NA, nrow = n_iter,
                  ncol = 138)

meanz_out <- matrix(data = NA, nrow = n_iter, 
                    ncol = 138)
zed <- vector(length=138)

# Values needed for rtruncnorm function
ones <- matrix(data = 1, nrow = 1, ncol = 138)
num_ones <- sum(dat$affected)
num_zeroes <- length(dat$affected[dat$affected == 0])

for (i in (1:n_iter)) {
  beta <- mvrnorm(n = 1, mu = beta_hat, XtXi)
  zed[y == 0] <- rtruncnorm(n=num_zeroes, 
                            mean = X[y ==0,] %*% beta,
                            sd = 1, a = -Inf, b = 0)
  
  zed[y == 1] <- rtruncnorm(n=num_ones, 
                            mean = X[y == 1,] %*% beta, 
                            sd = 1, a = 0, b = Inf)
  # Update out beta_hat based off previous iteration
  beta_hat <- solve(XtX) %*% t(X) %*% zed
  # Outputs
  beta_out[i, ] <- beta
  zed_out[i, ] <- zed
  meanz_out[i,] <- X %*% beta
}

mean_z <- rep(0,138)
for(m in (1:n_samples)){
  for(j in (1:138)) {
    mean_z[j] <- mean_z[j] +
      meanz_out[burn_in+(20*m),j]
  }
}
mean_z <- mean_z/n_samples

## Warning in kable_pipe(x = structure(c("0.4177726", "-0.4177625", "-0.6014916", :
## The table should have a header (column names)

dam <- c(5,1,1,1,0,0,0,0,0,0,1,0,1,0,0,
         0,0,1,0,0,0,0,0)
exposed <- c(rep(6,23))

proportion <- dam / exposed

temp <- c(53,57,58,63,66,67,67,67,68,69,70,70,
          70,70,72,73,75,75,76,76,78,79,81)

xmat <- cbind(rep(1, length(temp)), temp)
xmat2 <- cbind(rep(1, length(temp)), temp, temp^2)

• Bayesian Analysis of Binary and Polychotomous Response Data
  Author(s): James H. Albert and Siddhartha Chib
  Source: Journal of the American Statistical Association, Vol. 88, No. 422 (Jun., 1993), pp. 669- 679
  
• Dr. E. Olusegun George’s Lectures