Lektion 6 – logistische Regression

In Lektion 5 wurde die lineare Regression behandelt. In dieser Lektion geht es nun um Modelle mit einer binären abhängigen Variable.

Logistische Regression

Im folgenden Video wird die logistische Regression erklärt. Im Anschluss folgt eine Zusammenfassung der wichtigsten Punkte.

Mit einer linearen Regression können wir kontinuierliche abhängige Variablen analysieren. Ist die abhängige Variable hingegen binär, werden die Annahmen von OLS verletzt. Die Lösung ist ein anderes Modell: die Logit Regression (oder Probit Regression, wobei diese hier nicht diskutiert wird, da sie sich nur minimal unterscheidet im Ergebnis).

Eine logistische Regression kann mit der Funktion glm() geschätzt werden, wobei family = binomial gesetzt werden muss. Der Rest der Spezifikation ist identisch zur linearen Regression.

df_selects = schlegel::selects2015
model_logit = glm(participation ~ age * gender, family = binomial,
data = df_selects)
summary(model_logit)

Dem Output können wir entnehmen, dass Alter sowie die Interaktion von Alter und Geschlecht signifikant sind. Je älter eine Person ist, desto höher die Wahlbeteiligung, wobei der Effekt für Frauen geringer ist. Über die genaue Grösse können wir mit diesem Output jedoch noch nichts aussagen, da sich die Betas auf der logistischen Skala befinden und nicht auf die linearen.

Um die Outputs richtig interpretieren zu können, haben wir zwei Möglichkeiten: Odds Ratio und vorausgesagte Wahrscheinlichkeiten. Odds Ratios sind eine Besonderheit der logistischen Regression. Die Berechnung folgt über die Formel eβ.

exp(coef(model_logit))

Die odds ratio sagen uns, wie stark sich die Wahrscheinlichkeit der abhängige Variable ändert, wenn sich die unabhängige Variable um eins ändert. Für Alter als Beispiel bedeutet das, dass ein Mann, welcher ein Jahr älter ist, eine um das 1.03-Fache höhere Wahrscheinlichkeit hat an der Wahl teilzunehmen. Noch einfacher ist die Interpretation, wenn vom Wert ein abgezogen wird. Dann können die Odds Ratio in Prozentwerten interpretiert werden. Eine um ein Jahr älterer Mann hat also eine um 3 Prozent höhere Wahrscheinlichkeit teilzunehmen.

Die zweite Möglichkeit der Interpretation ist die vorausgesagte Wahrscheinlichkeit. Dazu muss der vorausgesagte Wert auf der Logitskale berechnet werden (ŷ). Anschliessend wird auf diesen die invertierte Linkfunktion angewendet um als Ergebnis vorausgesagte Wahrscheinlichkeiten zu erhalten. Am einfachsten geht dies mit der Funktion predicts() aus glm.predict.

library(glm.predict)
predicts(model_logit, "20-60,10;F")

Zum Beispiel können wir in Zeile 8 ablesen, dass eine 40 jährige Frau eine Wahrscheinlichkeit von 63.4 Prozent hat an der Wahl teilzunehmen, wobei das 95%-Konfidenzintervall von 61.4 bis 65.4 Prozent geht. Der Output kann auch mit ggplot2 dargestellt werden. Beispiele dazu findest du im Video.

Um mehrere Modell zu vergleichen, kann das Baysianische Informationskriterium (BIC) der Modelle vergleichen werden. Voraussetzung ist, dass die Fallzahl gleich ist. Das Modell mit dem tieferen Wert ist besser. In R kann dies mit der Funktion BIC() berechnet werden.

BIC(model_logit)

Wer ein besseres Verständnis davon haben will, wir vorausgesagte Wahrscheinlichkeiten von von Hand berechnet werden können inkl. Konfidenzintervall, kann sich das folgende Video anschauen. Darin wird sowohl die Methode der Simulation als auch diejenige mit Bootstrap erklärt. Unter dem Video folgt noch ein kurzes Beispiel mit Simulation und Bootstrap.

Im ersten Schritt simulieren wir 1000 Betas, welche sich alle leicht unterscheiden. Die Simulation macht sich zu nutzte, dass die Betas bei hoher Fallzahl der multivariaten Normalverteilung folgen.

betas = coef(model_logit)
vcov = vcov(model_logit)
betas_simulated = MASS::mvrnorm(1000, betas, vcov)

Die andere Variante ist, die 1000 verschiedenen Betas mit Hilfe von Bootstrap zu berechnen. Dazu berechnen wir das Modell 1000 Mal, nehmen aber jedes Mal leicht anderen Daten, indem wir alle Fälle aus dem Datensatz ziehen, sie aber immer wieder zurücklegen. Damit kommen gewisse Fälle mehrmals vor und anderen gar nicht. In einem Modell ist der Outlayer vielleicht gar nicht enthalten, im anderen mehrmals. So ergeben sich verschiedenen Betas.

boot_function = function(x, data){
  df_sample = data[sample(seq_len(nrow(data)), replace = TRUE), ]
  model = glm(participation ~ age * gender, family = binomial,
              data = df_sample) 
  coef(model)
}

betas_boot = do.call('rbind',
                     lapply(1:1000, boot_function, 
                            data = df_selects)
                     )

Mit der Funktion do.call('rbind',...) machen wir aus der lapply-Liste eine Matrix. Nun haben wir die 1000 Betas und können mit denen weiterarbeiten. Ich nehme betas_boot, man kann es aber einfach durch betas_simulated ersetzen. Als erstes müssen wir Werte für x festlegen, wir nehmen hier eine 65-jährige Frau, Der erste Wert ist immer 1 für den Achsenabschnitt. Der zweite ist 65 für das Alter, danach 1 für Frau und 65 für Alter * Frau. Anschliessend berechnen wir ŷ auf der Logit-Skala auf welches wir anschliessend die invertierte Linkfunktion für Logit anwenden. Mit diesen 1000 vorausgesagten Wahrscheinlichkeiten können wir anschliessend Durchschnitt und 95%-Konfidenzintervall berechnen.

x = c(1, 65, 1, 65)
yhat = betas_boot %*% x
pred_boot = exp(yhat) / (1 + exp(yhat))
mean(pred_boot)
quantile(pred_boot, probs = c(0.025, 0.975))

Ich erhalte einen vorausgesagte Wahrscheinlichkeit für eine 65-jährige Frau von 75.8% und ein Konfidenzintervall von 73.5% bis 78.0%. Für diskrete Unterschiede berechnet man die 1000 vorausgesagten Wahrscheinlichkeiten 2x für verschiedenen Fälle und nimmt dann die 1000 Unterschiede. Von diesen kann wiederum Durchschnitt und Konfidenzintervall berechnet werden.