#Challenger Disaster Redo

#==============================================================
#=================================================================
#==================================================================

# https://www.youtube.com/watch?v=ZOzoLdfWyKw

    setwd("~/workshops/challenger")
    remove(list = ls())
    options(digits = 3, scipen = 9999, knitr.table.format = "rst", length = 120)
    pacman::p_load(Amelia, sjPlot, knitr,readr,dplyr,ggplot2, stargazer)

    #read data
    launch = read.csv("challenger_V2.csv", header = T)

              # Replacing values higher than 1 by 1
              launch$distress_ct = ifelse(launch$distress_ct<1,0,1)
              table(launch$distress_ct)
## 
##  0  1 
## 16  7
      # 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 1 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 1 75 200 6 0 76 200 6 1 58 200 ======= =========== =========== ====================

      stargazer(launch, type = "text")
## 
## ==================================================================
## Statistic            N   Mean   St. Dev. Min Pctl(25) Pctl(75) Max
## ------------------------------------------------------------------
## n.rings              23  6.000   0.000    6     6        6      6 
## distress_ct          23  0.304   0.470    0     0        1      1 
## temperature          23 69.600   7.060   53     67       75    81 
## field_check_pressure 23 152.000  68.200  50     75      200    200
## ------------------------------------------------------------------
# Fitting a logistic regression model
model <- glm(distress_ct ~temperature + field_check_pressure, family = "binomial", data = launch)
summary(model)
## 
## Call:
## glm(formula = distress_ct ~ temperature + field_check_pressure, 
##     family = "binomial", data = launch)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.199  -0.578  -0.425   0.352   2.145  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)          13.29236    7.66397    1.73    0.083 .
## temperature          -0.22867    0.10999   -2.08    0.038 *
## field_check_pressure  0.01040    0.00898    1.16    0.247  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 18.782  on 20  degrees of freedom
## AIC: 24.78
## 
## Number of Fisher Scoring iterations: 5
    #plot the likelihood of failure
    plot_model(model, type = "pred", terms = "temperature")

                  #how good is the model? ######################
                  (  mod_fit = ifelse(model$fitted.values >= 0.5, 1,0)  )
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
##  0  0  0  0  0  0  0  0  1  1  0  0  1  1  1  0  0  0  0  0  0  0  1
                        table(mod_fit, launch$distress_ct)
##        
## mod_fit  0  1
##       0 14  3
##       1  2  4
                  (  Accuracy= mean(mod_fit == launch$distress_ct)  )
## [1] 0.783
                  paste("Accuracy=", round(Accuracy,2),"%")
## [1] "Accuracy= 0.78 %"
                  ##############################################
#But what about the likelihood of failure at 30 deg C?
                  
launch_30 = launch %>% dplyr::select(-distress_ct) %>% 
  rbind(c(n.rings = 6, temperature = 30, field_check_pressure = 200))
      launch_30
##    n.rings temperature field_check_pressure
## 1        6          66                   50
## 2        6          70                   50
## 3        6          69                   50
## 4        6          68                   50
## 5        6          67                   50
## 6        6          72                   50
## 7        6          73                  100
## 8        6          70                  100
## 9        6          57                  200
## 10       6          63                  200
## 11       6          70                  200
## 12       6          78                  200
## 13       6          67                  200
## 14       6          53                  200
## 15       6          67                  200
## 16       6          75                  200
## 17       6          70                  200
## 18       6          81                  200
## 19       6          76                  200
## 20       6          79                  200
## 21       6          75                  200
## 22       6          76                  200
## 23       6          58                  200
## 24       6          30                  200
          #predict likelihod of failure
          model_predict = predict(model, newdata = launch_30, type = "response")

            launch_30 = data.frame(cbind(model_predict, launch_30))
        launch_30 = arrange(launch_30, temperature) 
        
        #Graph No. 1
      plot(smooth(launch_30$model_predict) ~smooth(launch_30$temperature),  type = "l", col = "darkred", 
        ylab = "Probability of Failure", xlab = "Temperature (F)", 
            main = "Launch Probability of Failure")

#Graph No. 2
# 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 = "darkred", 
      main = "Launch Probability of Failure")

#===========================================================