Als Beispiel dienen Daten der Challenger-Katastrophe [Am 28. Januar 1986 explodierte die
Raumfähre Challenger 73 Sekunden nach dem Start. Der Ausfall eines oder mehrerer Dichtungsringe in einer der seitlichen Feststoffrakten wurde als Grund der Katastrophe ermittelt.]
Probleme mit den Dichtungsringen sind auch vorher am Boden aufgetreten. So liegen Erhebungsdaten vor, die den Defekt eines Rings in Zusammenhang mit der Außentemperatur darstellen können: Je kälter es draußen ist, desto höher die Ausfallwahrscheinlichkeit eines Rings. Bedauerlicher Weise wurden diese Daten erst rückwirkend ausgewertet.
Diese Daten übertragen wir wie folgt in R:
Temp = c(66,67,68,70,72,75,76,79,53,58,70,75,67,67,69,70,73,76,78,81,57,63,70)
Ausfall = factor(c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1))
shuttle = data.frame(Temp, Ausfall)
colnames(shuttle) = c("Temp", "Ausfall") |
Die Kodierung des factor Ausfall lautet: 0=kein Ausfall | 1=Ausfall
Die Einheit für Temp lautet Fahrenheit.
Der Aufruf zur logistischen Regression [Ausfall erklärt durch Temperatur] lautet:
fit = glm(Ausfall~Temp,data=shuttle, family=binomial); fit |
Dies liefert folgenden Output:
Call: glm(formula = Ausfall ~ Temp, family = binomial, data = shuttle)
Coefficients:
(Intercept) Temp
15.0429 -0.2322
Degrees of Freedom: 22 Total (i.e. Null); 21 Residual
Null Deviance: 28.27
Residual Deviance: 20.32 AIC: 24.32
Die Odds Ratio (OR) errechnet sich wie folgt:
exp(-0.2322)
[1] 0.7927875 |
Dies liefert das Ergebnis OR = 0.7927875. Die Interpretation lautet: Mit jedem höhren Grad an Temparatur ändert sich die Chance auf einem Ausfall um den Faktor 0.7927. Das bedeutet, dass die Chancen eines Ausfalls niedriger werden, je höher die Außentemperatur steigt.
In der Nacht vor dem Challengerstart konnte eine Temperatur von 37°Fahrenheit gemessen werden. Auf Grundlage der logistischen Regression lässt sich nun die Ausfallwahrscheinlichkeit bei 37° Außentemperatur berechnen
# zuerst eine Hilfsfunktion zum umrechnen
y = function(x){
blup = exp(15.0429 + (x*(-0.2322)))
bla = blup/(1+blup)
return(bla)} |
# jetzt setzen wir 37 ein:
y(37)
[1] 0.9984243 |
Das Modell sagt bei einer Außentemperatur von 37° eine Ausfallwahrscheinlichkeit von 99,84% voraus.
Die Funktion hätten wir uns auch sparen können, da der Befehl predict diese Arbeit übernimmt.
predict(fit, data.frame(Temp=37),type="response") |
Jetzt kann man dies für jede mögliche Temperatur durchmachen, und schon erhält man dieses schöne Plot.
y = predict(fit, data.frame(Temp=30:90),type="response")
plot(30:90, y, type="l", ylab="Ausfallswahrscheinlichkeit", xlab="Temperatur in F")
abline(v=seq(20,100,1),h=seq(-0.2,1.2,0.05),col="gray") #Mathegitter
abline(v=37,col="red") # 37-Grad
abline(h=0.5,col="black") # Mittellinie
points(30:90, y, type="l",lwd=2,col="blue") # Schätzwerte in "blau" |
