Statistik II
Maximum-Likelihood Schätzung bei der logistischen Regression
AG Forschungsmethoden und Kognitive Psychologie, Institut für Psychologie, Universität Bremen
Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.
Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.
Versionshistory:
- v1.0: erste online-gestellte Version (04.03.2024)
1 Maximum-Likelihood Schätzung
Wir haben in der Vorlesung zur logistischen Regression (Teil 9 von Statistik II) nur in Worten erläutert, was eine iterative Maximum-Likelihood (ML) Schätzung ist. In dieser Ergänzung wollen wir darauf noch einmal etwas genauer eingehen, da die ML-Schätzmethode durchaus in vielen Feldern ihre Anwendung findet. Im Allgemeinen versteht man unter ML-Schätzung, dass wir denjenigen Parameter(satz) auswählen, der am wahrscheinlichsten den empirischen Daten zugrunde liegt. Für weitere Ausführungen sei auch auf Teil 6 von Statistik II verwiesen; hier wenden wir die Grundlagen auf den Fall der logistischen Regression an.
1.1 Likelihood
Der Begriff Likelihood an sich beschreibt erstmal nur einen Wert auf einer Wahrscheinlichkeitsverteilung. Als Beispiel betrachten wir im Folgenden einen einzelnen Münzwurf, dessen Ergebnisse Kopf (kodiert mit 1) und Zahl (kodiert mit 0) sind. Wir gehen allerdings nun davon aus, dass wir nicht wissen mit welcher Wahrscheinlichkeit Kopf bzw. Zahl (also 1 bzw. 0) auftreten (uns ist also unklar, ob die Münze wirklich fair ist). Aus Teil 9 von Statistik I wissen wir schon, dass sich ein solches Zufallsexperiment durch eine bernoulli-verteilte Zufallsvariable \(Y\) mit einem unbekannten Parameter modellieren lässt; diesen unbekannten Parameter der Wahrscheinlichkeit nennen wir hier \(P\). Wir können also die Wahrscheinlichkeit für Kopf und Zahl als eine Funktion von \(P\) ausdrücken: \[\begin{equation*} \begin{aligned} \text{Kopf: }&P(Y=1|P)=P\\ \text{Zahl: }&P(Y=0|P)=1-P\\ \end{aligned} \end{equation*}\] Etwas verallgemeinert können wir somit für einen beliebigen Münzwurf \(i\) mit Ausgang \(y_i\) (also 1 oder 0) die Likelihood des Wurfs bestimmen durch: \[ P(Y = y_i | P) = P^{y_i}\cdot (1-P)^{(1-y_i)} \] Wenn \(y_i=1\) ist, dann steht dort genau \(P(Y=1|P)=P\) und wenn \(y_i=0\) ist, dann bleibt stehen \(P(Y=0|P)=1-P\).
1.2 Likelihood-Funktion
Typischerweise haben wir aber nicht nur eine Beobachtung (also einen Münzwurf), sondern \(n\)-viele Beobachtungen \(y_1,\ldots,y_n\), wobei alle Beobachtungen voneinander unabhängig seien. Die Variable \(y\) beinhalte nun alle diese einzelnen Beobachtungen. Die gemeinsame Wahrscheinlichkeit aller Würfe in Abhängigkeit des unbekannten Parameters \(P\) kann durch das Produkt der einzelnen Likelihoods jeden einzelnen Wurfs beschrieben werden. Hier bedienen wir uns des Produktzeichens \(\Pi\) (ein großes Pi), welches im Prinzip wie das Summenzeichen \(\Sigma\) zu verstehen ist, nur dass die einzelnen Terme nicht addiert, sondern multipliziert werden: \[P(Y=y | P) = \left[P^{y_1}\cdot(1-P)^{1-y_1}\right]\cdot\,\ldots\,\cdot \left[P^{y_n}\cdot(1-P)^{1-y_n}\right] =\prod_{i=1}^n \left[P^{y_i}\cdot(1-P)^{1-y_i}\right]\]
Wir schreiben diese Gleichung nun noch einmal auf und lassen dabei den Mittelteil weg, der nur deutlich machen sollte, wie das Produktzeichen zu verstehen ist: \[P(Y=y | P) = \prod_{i=1}^n \left[P^{y_i}\cdot(1-P)^{1-y_i}\right]\] Hier haben wir rechts vom Gleichheitszeichen einen Ausdruck, der sowohl \(P\) als auch \(y\) beinhaltet und der Teil links vom Gleichheitszeichen besagt, dass wir die Werte auf \(y\) unter der Bedingung \(P\) betrachten, also für bestimmte Werte für \(P\). Davon kann man sich aber auch lösen, denn eigentlich kennen wir ja \(y\), also die Ergebnisse aller Würfe in der Stichprobe. Das heißt aber auch, die einzige unbekannte Variable ist \(P\) und wir können somit unsere gemeinsam verteilte Wahrscheinlichkeitsfunktion (also den Teil rechts vom Gleichheitszeichen) auch als eine bedingte Funktion von \(P\) (gegeben die bekannten Daten \(y\)) schreiben: \[\mathcal{L}(P | y) = \prod_{i=1}^n [P^{y_i}\cdot(1-P)^{1-y_i}]\] Diese wird dann als sog. Likelihood-Funktion \(\mathcal{L}\) bezeichnet, da sie die “Likelihood” eines Parameters gegeben eine bestimmte Stichprobe darstellt.
1.3 Maximum-Likelihood
Was wir gerade getan haben heißt nun: Unter der Bedingung, dass wir die Ausgänge des Zufallsexperiments kennen, können wir bestimmen, wie wahrscheinlich es ist, dass ein bestimmter Parameter \(P\) diesen zugrunde liegt! Falls wir also tippen müssten, mit welcher Wahrscheinlichkeit wir Kopf werfen, dann würden wir denjenigen Wert \(\widehat{P}\) wählen, für den die Likelihood-Funktion \(\mathcal{L}\) ihr Maximum erreicht. Dieser Wert \(\widehat{P}\) wird dann als Schätzer für den gesuchten Parameter \(P\) herangezogen und daher kommt der Name Maximum-Likelihood Schätzung.
Diesen letzten Punkt schauen wir noch an einem kurzen Beispiel an: Nehmen wir an, wir haben 120-mal Kopf und 80-mal Zahl geworfen (d.h. die Variable \(y\) enthält 120-mal eine 1 und 80-mal eine 0). Wenn wir diese Werte in die Likelihood-Funktion \(\mathcal{L}\) einsetzen, ergibt sich: \[\mathcal{L}(P|y) = P^{120} \cdot (1-P)^{80} \] Nun könnten wir (hier noch sehr einfach) das Maximum dieser Funktion analytisch bestimmen; wir gehen allerdings den einfacheren Weg und plotten diese Funktion einfach (beachten Sie die sehr kleine Skalierung der \(y\)-Achse):
P <- seq(0,1,0.01) # Vektor für P von 0 bis 1
L <- P^120 *(1-P)^80 # Werte auf Likelihood-Funktion
plot(P,
L,
type = "l",
ylab = "L(P|y)",
xlab = "P")
index <- which(L == max(L)) # wo ist...
P_max <- P[index] # ...das Maximum?
abline(v = P_max) # einzeichnen
Wir sehen recht deutlich, dass die Likelihood-Funktion ihr Maximum bei \(P = 0.6\) annimmt. Dieser Parameter scheint also angesichts unserer Datenlage am wahrscheinlichsten zu sein und wir würden daher davon ausgehen, dieser Parameter würde unseren Daten zugrunde liegen. Im vorliegenden Fall hätten wir auch einfach \(\frac{n(\text{Kopf})}{n} = \frac{120}{200} = 0.6\) rechnen können, was tatsächlich die analytische Lösung wäre.
Allerdings ist es auch so, dass ML-Schätzungen nicht immer durch eine analytische Lösung gefunden werden können. Das heißt, es gibt häufig keine geschlossene Formel, die immer eine optimale Lösung darstellt (wie bspw. im Zuge des ALM die Formel \(\hat{\beta} = \bf (X'X)^{-1}X'y\)). In solchen Fällen ist die Form der Likelihood-Funktion von Stichprobe zu Stichprobe verschieden und wir benötigen computergestützte Lösungen, welche iterativ die gegebene Oberfläche der Likelihood-Funktion systematisch nach ihrem Maximum absuchen. Iterativ meint dabei, dass die Algorithmen mit einem bestimmten Satz an Parameterwerten beginnen und i.W. schauen, ob sie durch systematische Veränderungen dieser Parameter ein “besseres Ergebnis” erzielen können, also eine höhere Likelihood. Dies ist beispielsweise bei der logistischen Regression der Fall, wie wir jetzt kurz sehen werden.
Abschließend wollen wir noch anmerken, dass häufig von der logarithmierten Likelihood-Funktion die Rede ist, der “log-Likelihood” \(\ell\). Diese erhält man, indem der natürliche Logarithmus auf die Likelihood-Funktion angewendet wird: \[ \ell=\ln(\mathcal{L}) \] Der Grund hierfür ist, dass dies aus den \(n\)-vielen Produkten der Likelihood-Funktion \(n\)-viele Summen macht, ohne dabei das gobale Maximum zu verändern. Dies hat vor allem mathematische Vorteile. Ein wichtiger Vorteil ist, dass die Werte der Likelihood-Funktion sehr klein sind (siehe die visualisierte Likelihood-Funktion weiter oben) und durch den Logarithmus wieder größer werden und dadurch für Computer auch handhabbarer werden. Daneben wird es dadurch manchmal möglich, analytische Lösungen zu finden.
1.4 Maximum-Likelihood Schätzung bei der logistischen Regression
Wir können das obige Beispiel sehr einfach erweitern, indem wir uns in Erinnerung rufen, dass wir den Prädiktorwert \(x_i\) einer Person mit Hilfe einer Regressionsgleichung direkt auf einen für die Person passenden Parameter der Wahrscheinlichkeit \(P_i\) beziehen können. Hierzu greifen wir auf die Formulierung der Regressionsgleichung bezogen auf Wahrscheinlichkeiten zurück, \[P_i = \frac{e^{a+bx_i}}{1+e^{a+bx_i}}\]
welche abhängig ist von den Parametern \(a\) und \(b\). Vorhin hatten wir die Likelihood-Funktion \(L\) geschrieben als: \[\mathcal{L}(P | y) = \prod_{i=1}^n P^{y_i}\cdot(1-P)^{1-y_i}\] Nun betrachten wir diese Funktion für \(a\) und \(b\) gegeben \(X\) und \(Y\), also Prädiktor- und Kriteriumswerten. Wir schreiben die Likelihood-Funktion daher als \(L(a,b|X,Y)\) und setzen für \(P\) die logistische Regressionsgleichung direkt ein:
\[\mathcal{L}(a,b | X,Y) = \prod_{i=1}^n \left[P_i\right]^{y_i}\cdot \left[1-P_i\right]^{1-y_i} = \prod_{i=1}^n \left[\frac{e^{a+bx_i}}{1+e^{a+bx_i}}\right]^{y_i}\cdot\left[1-\frac{e^{a+bx_i}}{1+e^{a+bx_i}}\right]^{1-y_i}\] Auch diese Likelihood-Funktion beschreibt die Wahrscheinlichkeit von Parametern (\(a\) und \(b\)) in Abhängigkeit der gegebenen Stichprobenwerte (\(X\) und \(Y\)). Geometrisch haben wir allerdings nun keine eindimensionale Kurve mehr, sondern eine zweidimensionale Oberfläche, für die wir ein Maximum suchen. Hier sieht man bereits, dass die Sache dann deutlich kniffeliger wird, sodass wir auf Optimierungsalgorithmen zurückgreifen müssen.
Im folgenden R-Code demonstrieren wir einmal, wie eine solche Optimierung aussehen kann (vielleicht findet das jemand schön :)) Dazu benutzen wir die Daten, die wir auch in Teil 9 von Statistik II schon genutzt haben:
daten <- read.table("./Daten/log_regression_daten.dat",
header = TRUE,
sep = "\t")Zur Bestimmung der Parameter nutzen wir hier die Funktion
optim(), die allgemeine Optimierungsfunktionen in R zur
Verfügung stellt. Im Wesentlichen formulieren wir dazu zunächst eine
Funktion log_like() die händisch die logarithmierte
Likelihood für einen bestimmten Satz an Parametern \(a\) und \(b\) berechnet und zurückgibt:
# Log-Likelihood-Funktion
# Erhält die Koeffizienten a und b (innerhalb von coefs)
# sowie die empirischen Prädiktor- (X) und Kriteriumswerte (Y)
# gibt die log-Likelihood zurück
log_like <- function(coefs,Y,X) {
a <- coefs[1] # Koeffizient a
b <- coefs[2] # Koeffizient b
# Bezug der Regressionsgleichung auf Wahrscheinlichkeiten (s.o.)
P_i <- exp(a+b*X)/(1+exp(a+b*X))
# Wert der Log-Likelihood-Funktion bei gegebenem a und b (s.o.)
log_like_Wert <- log(prod(P_i^Y*(1-P_i)^(1-Y)))
# Zwischenausgabe, zum "Beobachten" der Optimierung; nicht relevant
cat(round(a,2), "\t", round(b,2), "\t", round(log_like_Wert,2),"\n")
#Rückgabe des Wertes
return(log_like_Wert)
}Diese Funktion gibt also den Wert zurück, der maximiert werden soll.
Nun können wir optim() aufrufen und sehen auch, was hier
passiert (weil wir in die Funktion log_like() eine
Zwischenausgabe eingebaut haben): Der Algorithmus startet mit den
Startwerten \(a=0\) und \(b=0\) und es wird eine log-Likelihood von
\(-27.73\) berechnet. Dann wird \(a\) verändert und eine neue Log-Likelihood
berechnet. Im Laufe der Zeit wird die Log-Likelihood immer größer, bis
am Ende die kleinen Veränderungen (teilweise auf weiteren
Nachkommastellen) nichts mehr bringen und die Log-Likelihood den Wert
nicht mehr verändert. Dieser Wert wird dann als das Maximum
angenommen:
# Suche nach dem Maximum mit einem generellen Optimierungsalgorithmus
# 'optim' (Nelder-Mead)
mle.res <- optim(par = c(0, 0), # Startwerte coefs (wird von optim variiert)
fn = log_like, # Funktion die maximiert werden soll
Y = daten$verheiratet, # Kriterium --> wird an log_like übergeben (bleibt fix)
X = daten$alter, # Prädiktor --> wird an log_kike übergeben (bleibt fix)
control = list(fnscale = -1)) # führt zur Maximierung statt Minimierung## 0 0 -27.73
## 0.1 0 -27.78
## 0 0.1 -54.5
## 0.1 -0.1 -67.83
## 0.03 0.05 -34.75
## 0.08 -0.05 -41.06
## 0.04 0.03 -28.92
## 0.06 -0.03 -31.99
## 0.04 0.01 -27.67
## -0.06 0.01 -27.31
## -0.13 0.02 -27.26
## -0.09 0.03 -28.94
## -0.02 0.01 -27.37
## -0.2 0.01 -26.93
## -0.32 0.01 -26.71
## -0.43 0.03 -26.43
## -0.64 0.03 -26.05
## -0.83 0.03 -25.57
## -1.18 0.04 -25.14
## -1.49 0.06 -24.07
## -2.08 0.08 -22.93
## -2.62 0.08 -22.76
## -3.61 0.1 -22.67
## -4.51 0.14 -19.69
## -6.17 0.2 -17.91
## -7.7 0.22 -21.22
## -6.29 0.19 -20.18
## -8.86 0.28 -16.34
## -11.49 0.37 -14.78
## -11.37 0.38 -13.97
## -13.91 0.48 -13.47
## -19.23 0.65 -12.65
## -25.76 0.87 -13.04
## -21.65 0.76 -14.5
## -19.11 0.66 -13.42
## -24.43 0.83 -13.05
## -21.8 0.74 -12.93
## -21.92 0.73 -12.56
## -23.33 0.77 -12.81
## -19.35 0.64 -12.77
## -19.96 0.66 -12.59
## -22.65 0.75 -12.74
## -20.09 0.67 -12.57
## -22.04 0.74 -12.6
## -20.48 0.68 -12.56
## -22.32 0.74 -12.59
## -20.64 0.69 -12.56
## -22.08 0.74 -12.56
## -21.68 0.72 -12.55
## -22.96 0.77 -12.57
## -21.22 0.71 -12.55
## -20.98 0.7 -12.56
## -21.69 0.72 -12.55
## -21.23 0.71 -12.55
## -21.34 0.71 -12.55
## -20.88 0.7 -12.55
## -21.48 0.72 -12.55
## -21.6 0.72 -12.55
## -21.32 0.71 -12.55
## -21.46 0.72 -12.55
## -21.43 0.72 -12.55
## -21.26 0.71 -12.55
## -21.43 0.72 -12.55
## -21.54 0.72 -12.55
## -21.37 0.71 -12.55
## -21.37 0.71 -12.55
## -21.39 0.71 -12.55
## -21.33 0.71 -12.55
## -21.36 0.71 -12.55
## -21.34 0.71 -12.55
## -21.38 0.71 -12.55
## -21.39 0.71 -12.55
## -21.38 0.71 -12.55
## -21.38 0.71 -12.55
## -21.38 0.71 -12.55
## -21.37 0.71 -12.55
## -21.38 0.71 -12.55Nun berechnen wir die logistische Regression mit glm()
und vergleichen dessen Koeffizienten für \(a\) und \(b\) mit denen aus der Optimierung:
glm.res <- glm(verheiratet ~ alter,
data = daten,
family = binomial("logit"))
coef(glm.res) # Koeffizienten von glm()## (Intercept) alter
## -21.3794357 0.7140123mle.res$par # a und b nach Optimierung mit optim()## [1] -21.3808584 0.7140625Bis auf Nachkommastellen kommen beide Verfahren zum gleichen
Ergebnis. Die minimale Verschiedenheit rührt daher, dass
glm() einen anderen Optimierungsalgorithmus verwendet.