# Number of rocket booster 'O' rings seen to be damaged on each previous shuttle flight prior to the Challenger disaster of 20th January 1986
# Flury, B.D. (1997) A First Course in Multivariate Statistics, New York: Springer
# 'challenger' is a data frame with 23 observations on the following 2 variables.
# Temp:  Ambient temperature (Fahrenheit) at launch time
# Damage: Number of damaged 'O' rings (out of a total of 6)

install.packages() # pick a server
install.packages('Flury')
install.packages('lme4')  # be prepared: this one takes a few minutes to install 

require('Flury')
require('lme4')

data(challenger)

# Binomial model for probability of any one ring failing;
# independence between rings and launches, adjusted for temperature
ch <- glm(cbind(Damage, 6-Damage) ~ Temp, family = binomial, data = challenger)


with(challenger, plot(Temp, Damage/6,xlim=c(20,100),ylim=c(0,1)))
lines(seq(23,100,1),predict.glm(ch,newdata=data.frame(Temp=seq(23,100,1)),type="response"))
abline(v=32, col = "red", lwd = 2) ## temp when challenger launched
lines(seq(23,100,1),1-(1-predict.glm(ch,newdata=
      data.frame(Temp=seq(23,100,1)),type="response"))^6,col='red')


# Binomial model, for probability of at least one ring failing;
#  independence between launches, adjusted for temperature
ch2 <- glm((Damage>0) ~ Temp, family = binomial, data = challenger)

x11()
with(challenger, plot(Temp, Damage/6,xlim=c(20,100),ylim=c(0,1)))
lines(seq(23,100,1),predict.glm(ch2,newdata=data.frame(Temp=seq(23,100,1)),type="response"),col='green')
abline(v=32, col = "red", lwd = 2) ## temp when challenger launched


# Equivalent Bernouli model, for probability of a ring failing;
chal <- cbind(challenger$Damage, challenger$Temp)
chalb = matrix(rep(0,6*23*3),23*6,3)
chalb[,3]=gl(23,6,23*6)
for (i in 1:23) {
for (j in 1:6) {
chalb[((i-1)*6)+j,2] = chal[i,2]
}}

for (i in 1:23) {
for (j in 1:chal[i,1]) {
chalb[((i-1)*6)+j,1] = ((chal[i,1]>0)*1)
}}

# Equivalent Bernouli model, for probability of a ring failing;
#  independence between launches, adjusted for temperature
ch25 <- glm(chalb[,1] ~ chalb[,2], family = binomial)



####### RESULTS ##############################3
#
> summary(ch)

Call:
glm(formula = cbind(Damage, 6 - Damage) ~ Temp, family = binomial, 
    data = challenger)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.95227  -0.78299  -0.54117  -0.04379   2.65152  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  5.08498    3.05247   1.666   0.0957 .
Temp        -0.11560    0.04702  -2.458   0.0140 *
---
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: 18.086  on 21  degrees of freedom
AIC: 35.647

Number of Fisher Scoring iterations: 5

> summary(ch25)

Call:
glm(formula = chalb[, 1] ~ chalb[, 2], family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7774  -0.3677  -0.3106  -0.2209   2.6879  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  5.08498    3.05248   1.666   0.0957 .
chalb[, 2]  -0.11560    0.04702  -2.458   0.0140 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 66.540  on 137  degrees of freedom
Residual deviance: 60.396  on 136  degrees of freedom
AIC: 64.396

Number of Fisher Scoring iterations: 6