Analysis of the Challenger disaster data using logistic regression

The description of this data can be found at:

Challenger Explosion

A brief copy-paste description is as follows:

The NASA space shuttle Challenger exploded on January 28, 1986, just 73 seconds after liftoff, bringing a devastating end to the spacecraft’s 10th mission. The disaster claimed the lives of all seven astronauts aboard, including Christa McAuliffe, a teacher from New Hampshire who would have been the first civilian in space. It was later determined that two rubber O-rings, which had been designed to separate the sections of the rocket booster, had failed due to cold temperatures on the morning of the launch. The tragedy and its aftermath received extensive media coverage and prompted NASA to temporarily suspend all shuttle missions.

# Delete memory
rm(list=ls())

# Required packages
library(knitr)

# Challenger data set
distress_ct = c(0,1,0,0,0,0,0,0,1,1,1,0,0,2,0,0,0,0,0,0,2,0,1)
temperature = c(66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81,76,79,75,76,58)
field_check_pressure = c(50,50,50,50,50,50,100,100,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200)
n.rings <- rep(6,length(distress_ct))

launch <- data.frame(cbind(n.rings,distress_ct, temperature, field_check_pressure))

#launch <- launch[9:23,]

# Replacing values higher than 1 by 1 (as we are only interested in whether or not a failure happened)
#launch$distress_ct = ifelse(launch$distress_ct<1,0,1)

# Displaying the data
kable(launch)
n.rings distress_ct temperature field_check_pressure
6 0 66 50
6 1 70 50
6 0 69 50
6 0 68 50
6 0 67 50
6 0 72 50
6 0 73 100
6 0 70 100
6 1 57 200
6 1 63 200
6 1 70 200
6 0 78 200
6 0 67 200
6 2 53 200
6 0 67 200
6 0 75 200
6 0 70 200
6 0 81 200
6 0 76 200
6 0 79 200
6 2 75 200
6 0 76 200
6 1 58 200 
# Fitting a logistic regression model
model <- glm(cbind(distress_ct,n.rings-distress_ct) ~temperature + field_check_pressure, family = binomial(link='logit'), data = launch)
summary(model)
## 
## Call:
## glm(formula = cbind(distress_ct, n.rings - distress_ct) ~ temperature + 
##     field_check_pressure, family = binomial(link = "logit"), 
##     data = launch)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0361  -0.6434  -0.5308  -0.1625   2.3418  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           2.520195   3.486784   0.723   0.4698  
## temperature          -0.098297   0.044890  -2.190   0.0285 *
## field_check_pressure  0.008484   0.007677   1.105   0.2691  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 24.230  on 22  degrees of freedom
## Residual deviance: 16.546  on 20  degrees of freedom
## AIC: 36.106
## 
## Number of Fisher Scoring iterations: 5
#------------------------------------------------------------------
# Confidence intervals based on the profile likelihood
#------------------------------------------------------------------
kable(confint(model, level=0.95, type="LR"),digits = 4)
## Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) -4.3229 9.7726
temperature -0.1941 -0.0136
field_check_pressure -0.0043 0.0289 
#------------------------------------------------------------------
# Wald Confidence intervals
#------------------------------------------------------------------
kable(confint.default(model, level=0.95),digits = 4)
2.5 % 97.5 %
(Intercept) -4.3138 9.3542
temperature -0.1863 -0.0103
field_check_pressure -0.0066 0.0235 
#------------------------------------------------------------------
# Normal confidence intervals
#------------------------------------------------------------------
NCI <- cbind(coefficients(model)- 1.96*summary(model)$coefficients[,2],
             coefficients(model)+ 1.96*summary(model)$coefficients[,2] )

colnames(NCI) <- c("2.5 %", "97.5 %")
rownames(NCI) <- c("(Intercept)", "temperature", "pressure")
kable(NCI,digits = 4)
2.5 % 97.5 %
(Intercept) -4.3139 9.3543
temperature -0.1863 -0.0103
pressure -0.0066 0.0235 
##################################################################################
# Brief interpretation of the model
##################################################################################

# Predicted probability of failure for the Challenger
plogis(coefficients(model)%*%c(1,31,200))
##           [,1]
## [1,] 0.7631088
# Predicted probability of failure for minimum temperature experienced previously
plogis(coefficients(model)%*%c(1,53,200))
##           [,1]
## [1,] 0.2703734
# The following plot presents the estimated probability 
# as a function of the temperature at pressure = 200
curve(predict(model, data.frame(temperature = x, field_check_pressure = 200),
              type = "response"),  from = 30, to = 85, ylab = "Probability of Failure", xlab = "Temperature (F)",
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, col = 'red')