Analyse du risque de d?faillance des joints toriques de la navette Challenger - CORRECTION
Arnaud Legrand & cor.FMG
14/12/2021
CORRECTION PERSONNELLE
On va essayer de refaire l’analyse avec les éléments qui nous semblent “légers” :
- la disparition des données de malfonctions jugées inutiles, qui vont réduire mécaniquement l’incertitude et biaisé la prise de décision
- la sortie arbitraire du facteur pression dans la regression logistique
- ou alternativement (1), tenter de réduire la difficulté d’analyse des résultats de la régression par une approche plus conservatrice consistant à détecter toute malfonction comme la survenue d’un évènement à éviter, à partir du recodage de la variable en une variable binaire/booléenne 0/1 pure (=>1 pour la moindre malfonction de joint) => permettrait d’avoir exp(coefficient regression log) = odds-ratio qui associe la variable explicative à la variable expliquée)
Chargement et visualisation du dataset
On commence de la même façon par le chargement des données /création du dataset de travail, mais nous ne retirons pas les données jugées valables arbitrairement. Nous conservons toutes les données de disfonctionnement des joints. Cela augmente la taille de l’échantillon et devrait réduire l’incertitude générale.
On les affiches ensuite
data_f = read.csv("shuttle.csv",header=T)
data_f## Date Count Temperature Pressure Malfunction
## 1 4/12/81 6 66 50 0
## 2 11/12/81 6 70 50 1
## 3 3/22/82 6 69 50 0
## 4 11/11/82 6 68 50 0
## 5 4/04/83 6 67 50 0
## 6 6/18/82 6 72 50 0
## 7 8/30/83 6 73 100 0
## 8 11/28/83 6 70 100 0
## 9 2/03/84 6 57 200 1
## 10 4/06/84 6 63 200 1
## 11 8/30/84 6 70 200 1
## 12 10/05/84 6 78 200 0
## 13 11/08/84 6 67 200 0
## 14 1/24/85 6 53 200 2
## 15 4/12/85 6 67 200 0
## 16 4/29/85 6 75 200 0
## 17 6/17/85 6 70 200 0
## 18 7/29/85 6 81 200 0
## 19 8/27/85 6 76 200 0
## 20 10/03/85 6 79 200 0
## 21 10/30/85 6 75 200 2
## 22 11/26/85 6 76 200 0
## 23 1/12/86 6 58 200 1On poursuit pour l’instant la même procédure d’analyse. On va donc représenter les variations de dysfonctionnement des joints avec la température, sur les données non-tromquées
plot(data=data_f, Malfunction/Count ~ Temperature, ylim=c(0,1))
C’est pas super évident, mais il semble malgré tout que la température “réduit” les dysfonctionnement.
Ce qui est visible en revanche, si on essaie d’agrandir les échelles, c’est que nous n’avons qu’une plage très ressérées d’observations, représentant des évènements de températures assez hohmogènes
plot(data=data_f, Malfunction/Count~Temperature, xlim=c(30, 90), ylim=c(0,1))
Il était aussi possible de réaliser une visualisation plus précise, qui intégrait directement la regression logistique
library(ggplot2)## Registered S3 methods overwritten by 'tibble':
## method from
## format.tbl pillar
## print.tbl pillarggplot(data_f,aes(x=Temperature, y=Malfunction/Count)) + geom_point(alpha=.3, size=3) +
theme_bw() + geom_smooth(method = "glm", method.args = list(family = "binomial")) + xlim(30,90)## `geom_smooth()` using formula 'y ~ x'## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
On voit ici que la tendancce que nous donne la regression n’est déjà plus aussi sécurisante que les simple plots…
Estimation de l’impact de la température
Avec le modèle complet
logistic_reg1 = glm(data=data_f, Malfunction/Count ~ Temperature, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg1)##
## Call:
## glm(formula = Malfunction/Count ~ Temperature, family = binomial(link = "logit"),
## data = data_f, weights = Count)
##
## 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 .
## Temperature -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: 5La température est devenue significative ! L’analyse des coefficients restent cependant toujours compliqués (pas une véritable variable binaire qui permettrait d’avoir des odds-ratios par exponentielle des coefficients), bien que le signe de trompe pas : relation négative. Donc, plus la température monte, plus la proportion de défaillance des joints toriques baisse.
Essai twoby et Epi
Juste pour le fun
library(Epi)
malf_c <- (data_f$Malfunction/data_f$Count)
twoby2(1-malf_c, 1-data_f$Temperature)## 2 by 2 table analysis:
## ------------------------------------------------------
## Outcome : -80
## Comparing : 0.666666666666667 vs. 0.833333333333333
##
## -80 -78 P(-80) 95% conf. interval
## 0.666666666666667 0 0 NaN NaN NaN
## 0.833333333333333 0 0 NaN NaN NaN
##
## 95% conf. interval
## Relative Risk: NaN NaN NaN
## Sample Odds Ratio: NaN NaN NaN
## Conditional MLE Odds Ratio: 0 0 Inf
## Probability difference: NaN NaN NaN
##
## Exact P-value: 1.0000
## Asymptotic P-value: NaN
## ------------------------------------------------------Ben ça marche pas, pas des variables binaires… :-)
exp(-0.11560)## [1] 0.8908315Envrion 0.9… Cela caractérise bien la perte de proba à mesure que la température augmente (et inversement, 1,15 pour si la température chute)
Prédiction/Estimation de la probabilité de dysfonctionnant des joints toriques
On reprend le code proposé en modifiant les variables de calcul. Ici, on crée une séquence par “0,5” de 30 à 90 (abscisses axe), et on crée un jeu de données prédites par le modèle de la regression logistique (ordonnées axe). Ensuite, on le représente graphiqement
tempv = seq(from=30, to=90, by = .5)
rmv <- predict(logistic_reg1,list(Temperature=tempv),type="response")
plot(tempv,rmv,type="l",ylim=c(0,1))
points(data=data_f, Malfunction/Count ~ Temperature)
Là, cela devient évident! Avec des température de 30 degré farenheit, le risque de défaillance est très grand ! Nous pourrions nous en tenir là, et estimer les defauts par l’absence d’un dataset complet…
Intégration de la variable pression
logistic_reg2 = glm(data=data_f, Malfunction/Count ~ Temperature + Pressure + Temperature:Pressure, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg2)##
## Call:
## glm(formula = Malfunction/Count ~ Temperature + Pressure + Temperature:Pressure,
## family = binomial(link = "logit"), data = data_f, weights = Count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0279 -0.6600 -0.5141 -0.1744 2.3724
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.765547 46.285607 -0.470 0.638
## Temperature 0.252463 0.662989 0.381 0.703
## Pressure 0.130857 0.232643 0.562 0.574
## Temperature:Pressure -0.001769 0.003336 -0.530 0.596
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 24.230 on 22 degrees of freedom
## Residual deviance: 16.257 on 19 degrees of freedom
## AIC: 37.817
##
## Number of Fisher Scoring iterations: 6Pas d’interaction, ni d’effet simple : le retrait de la pression n’était donc pas dommageable selon les données disponibles
Modèle avec variables binaires : approche plus conservatrice
L’idée est de recoder la variable “incident” (ici la défaillance d’un joint torique) en une variablel binaire/bouléenne. Ainsi, on comptabilise le risque de défaillance d’au moins 1 joint torique (avec l’idée qu’un seul suffit à mettre en danger la navette)
Création de la variable binaire
data_f$Malfunction_b <- ifelse(data_f$Malfunction>=1, 1,0)On refait la regression
logistic_reg3 = glm(data=data_f, Malfunction_b ~ Temperature, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg3)##
## Call:
## glm(formula = Malfunction_b ~ Temperature, family = binomial(link = "logit"),
## data = data_f, weights = Count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5992 -1.8647 -0.9266 1.1081 5.4318
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.04290 3.01231 4.994 5.92e-07 ***
## Temperature -0.23216 0.04419 -5.254 1.49e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 169.60 on 22 degrees of freedom
## Residual deviance: 121.89 on 21 degrees of freedom
## AIC: 125.89
##
## Number of Fisher Scoring iterations: 5Analyse des odds-ratios
La librairie twoby ne fonctonnera pas. Mais ici, l’exp(coef) = OR!
exp(0.23216)## [1] 1.261322On trouve qu’il y a 1,3 fois plus de chance de défaillance avec la baisse de la temperature. Cette approche, comme on s’y attendait, montre plus facilemnt encore le risque qu’encourt la navette avec une plus forte probabilité signifiative qu’au moins 1 joint torique ne lâche avec la baisse de la température lors du lancement.
tempv = seq(from=30, to=90, by = .5)
rmv <- predict(logistic_reg3,list(Temperature=tempv),type="response")
plot(tempv,rmv,type="l",ylim=c(0,1))
points(data=data_f, Malfunction_b ~ Temperature)
En fait, si on en croit les prédictions du modèle, l’accident (défaillance d’au moins 1 joint torique) était inévitable en dessous de 50°F…