Modellregularisierung – The Bayesian Way

Einleitung

Manchmal kommt es bei statworx zu spontanen Diskussionen über statistische Methoden. In einer solchen Diskussion erklärte einer meiner Kollegen (wenn auch scherzhaft), dass Bayessche Statistik unnötig sei. Das brachte mich ins Grübeln: Warum sollte ich überhaupt bayessche Modelle im Kontext eines klassischen Regressionsproblems verwenden? Bestehende Ansätze wie Ridge Regression sind doch sicherlich genauso gut, wenn nicht besser. Allerdings hat der bayessche Ansatz den Vorteil, dass er das Modell sowohl regulieren kann, um Overfitting zu verhindern, und die Regularisierungsparameter sinnvoll interpretierbar macht. Anders als im üblichen Verständnis der Ridge Regression sind die Regularisierungsparameter nicht länger abstrakte Zahlen, sondern lassen sich im bayesschen Paradigma als Ausdruck vorheriger Überzeugungen deuten. In diesem Beitrag zeige ich dir die formale Ähnlichkeit zwischen einem generalisierten Ridge-Schätzer und dem bayesschen Äquivalent.

Ein (sehr kurzer) Einstieg in die bayessche Statistik

Um den bayesschen Regressionsschätzer zu verstehen, braucht man ein Mindestmaß an Wissen über die bayessche Statistik – hier das Wichtigste (falls du es nicht ohnehin schon kennst): In der bayesschen Statistik betrachten wir Modellparameter (also Regressionskoeffizienten) als zufällige Größen. Mit anderen Worten: Die uns vorliegenden Daten sind fixiert, die Parameter gelten als zufällig. Das widerspricht der klassischen, frequentistischen Perspektive, bei der die zugrunde liegenden Modellparameter als fix angesehen werden, während die Daten als zufällige Realisierungen eines durch die Parameter gesteuerten stochastischen Prozesses gelten. Ziel der bayesschen Analyse ist es, die sogenannte Posterior-Verteilung zu bestimmen, die du vielleicht noch aus dem Satz von Bayes kennst:

$$p(\theta|y) = \frac{p(y|\theta) p(\theta)}{p(y)}$$

Dabei ist $p(y|\theta)$ die Likelihood, $p(y)$ ein Normierungskonstante und $p(\theta)$ unsere Prior-Verteilung, die unabhängig von den Daten $y$ ist. In der klassischen Statistik setzt man $p(\theta) = 1$ (eine sogenannte improper reference prior), sodass beim Maximieren der posterioren „Wahrscheinlichkeit“ effektiv nur die Likelihood maximiert wird, da nur sie noch von $\theta$ abhängt. In der bayesschen Statistik hingegen verwendet man eine echte Wahrscheinlichkeitsverteilung für $p(\theta)$, z. B. eine Normalverteilung. Schauen wir uns nun ein Regressionsproblem an und nehmen an, dass sowohl unser Zielwert $y$ als auch unser Prior normalverteilt sind. Das führt uns zur konjugierten bayesschen Analyse, bei der wir eine saubere Gleichung für die Posterior-Verteilung formulieren können. In vielen Fällen ist dies nicht möglich – aus diesem Grund wurden Markov-Chain-Monte-Carlo-Methoden entwickelt, um aus der Posterior-Verteilung zu sampeln – ironischerweise ein frequentistischer Ansatz.

Wir treffen nun die übliche Annahme über die Daten: $y_i$ ist i.i.d. $N(\bold{x_i \beta}, \sigma^2)$ für alle Beobachtungen $i$. Das ergibt unsere Standard-Likelihood für die Normalverteilung. Jetzt können wir den Prior für den zu schätzenden Parameter $(\beta, \sigma^2)$ festlegen. Wenn wir einen Normalverteilungsprior (bedingt auf die Varianz $\sigma^2$) für den Gewichtsvektor $\beta$ annehmen, also $N(b_0, \sigma^2 B_0)$, und einen Invers-Gamma-Prior über die Varianz, dann lässt sich zeigen, dass die Posterior-Verteilung von $\beta$ normalverteilt ist mit dem Mittelwert

$$\hat\beta_{Bayesian} = (B_0^{-1} + X'X)^{-1}(B_0^{-1} b_0 + X'X \hat\beta)$$

Falls du an einem Beweis dieser Aussage interessiert bist, siehe Jackman (2009, S. 526).

Schauen wir uns das Stück für Stück an:

• $\hat\beta$ ist der Standard-OLS-Schätzer, $(X'X)^{-1}X'y$
• $b_0$ ist der Mittelwertsvektor der (multivariat normalen) Prior-Verteilung – damit lässt sich angeben, welche mittleren Werte wir für unsere Modellparameter erwarten
• $B_0$ ist die Kovarianzmatrix und enthält unsere jeweilige Unsicherheit bezüglich der Modellparameter. Der Kehrwert der Varianz wird als Präzision bezeichnet

Was wir aus der Gleichung sehen: Der Mittelwert der Posterior-Verteilung ist ein präzisionsgewichteter Mittelwert aus dem Prior-Mittelwert (also Informationen ohne Datenbasis) und dem OLS-Schätzer (nur datenbasiert). Der zweite Term in der Klammer zeigt, dass wir den unsicherheitsgewichteten Prior-Mittelwert $B_0^{-1} b_0$ zum gewichteten OLS-Schätzer $X'X\hat\beta$ addieren. Stell dir für einen Moment vor, $B_0^{-1} = 0$. Dann ergibt sich:

$$\hat\beta_{Bayesian} = (X'X)^{-1}(X'X \hat\beta) = \hat\beta$$

Das würde bedeuten, dass wir unendlich unsicher über unsere Prior-Annahmen sind, sodass der Mittelwertvektor der Prior-Verteilung verschwindet und nichts zur Posterior-Verteilung beiträgt. Umgekehrt: Wenn unsere Unsicherheit sinkt (und damit die Präzision steigt), trägt der Prior-Mittelwert $b_0$ stärker zum Posterior-Mittelwert bei.

Nach diesem kurzen Einstieg in die bayessche Statistik können wir nun den Ridge-Schätzer formal mit dem obigen bayesschen Schätzer vergleichen. Doch zuvor werfen wir einen Blick auf eine verallgemeinerte Form des Ridge-Schätzers.

Verallgemeinerung des Ridge-Schätzers

Ein Standardwerkzeug in vielen Regressionsproblemen ist der klassische Ridge-Schätzer, der durch Minimierung eines kleinsten-Quadrate-Problems aus folgender Verlustfunktion abgeleitet wird:

$$L(\beta,\lambda) = \frac{1}{2}\sum(y-X\beta)^2 + \frac{1}{2} \lambda ||\beta||^2$$

Während diese Minimierung den klassischen Ridge-Schätzer ergibt, wie man ihn aus Lehrbüchern kennt, gibt es eine leicht verallgemeinerte Form dieser Verlustfunktion:

$$L(\beta,\lambda,\mu) = \frac{1}{2}\sum(y-X\beta)^2 + \frac{1}{2} \lambda ||\beta - \mu||^2$$

Leiten wir den Schätzer her, indem wir die Verlustfunktion zunächst in Matrixschreibweise umformulieren:

$$\begin{aligned}L(\beta,\lambda,\mu) &= \frac{1}{2}(y - X \beta)^{T}(y - X \beta) + \frac{1}{2} \lambda||\beta - \mu||^2 \\&= \frac{1}{2} y^Ty - \beta^T X^T y + \frac{1}{2} \beta^T X^T X \beta + \frac{1}{2} \lambda||\beta - \mu||^2\end{aligned}$$

Durch Ableiten nach dem Parametervektor ergibt sich der Gradient:

$$\nabla_{\beta} L (\beta, \lambda, \mu) = -X^T y + X^T X \beta + \lambda (\beta - \mu)$$

Minimieren wir nach $\beta$, erhalten wir folgenden Ausdruck für den verallgemeinerten Ridge-Schätzer:

$$\hat\beta_{Ridge} = (X'X + \lambda I )^{-1}(\lambda \mu + X'y)$$

Der klassische Ridge-Schätzer ergibt sich, wenn man $\mu = 0$ setzt. Üblicherweise wird $\lambda$ als abstrakter Parameter betrachtet, der die Strafgröße reguliert, und $\mu$ als ein Vektor von Werten (einer pro Prädiktor), der die Verlustfunktion stärker steigen lässt, je weiter die Koeffizienten von diesen Werten abweichen. Ist $\mu = 0$, werden die Koeffizienten in Richtung null gezogen.

Schauen wir uns an, wie sich der Schätzer verhält, wenn die Parameter $\mu$ und $\lambda$ verändert werden. Wir definieren ein sinnvolles „Prior“ für unser Beispiel und variieren dann den Strafparameter. Als Beispiel verwenden wir den Datensatz diamonds aus dem ggplot2-Paket und modellieren den Preis als lineare Funktion der Karatzahl sowie der Merkmale depth, table, x, y und z.

Wie man aus dem Plot erkennt, ändern sich die Koeffizientenschätzungen – sowohl mit als auch ohne Prior – sehr schnell bei den ersten Erhöhungen der Strafgröße. Wir sehen auch den Shrinkage-Effekt im oberen Plot: Mit steigender Strafe nähern sich die Koeffizienten der Null an, einige schneller als andere. Der Plot rechts zeigt, wie sich die Koeffizienten verändern, wenn wir ein sinnvolles „Prior“ setzen. Die Koeffizienten verändern sich weiterhin, tendieren aber nun in Richtung des angegebenen Priors. Das liegt daran, dass $\lambda$ Abweichungen von $\mu$ bestraft – höhere Werte von $\lambda$ ziehen die Koeffizienten also stärker in Richtung $\mu$. Du fragst dich vielleicht, wie sich das im Vergleich zum bayesschen Schätzer verhält. Finden wir es heraus!

Vergleich von Ridge- und Bayesschem Schätzer

Nachdem wir nun sowohl den Ridge- als auch den bayesschen Schätzer gesehen haben, ist es Zeit für einen Vergleich. Wir haben festgestellt, dass der bayessche Schätzer den OLS-Schätzer enthält. Da wir dessen Form kennen, setzen wir sie ein und sehen, was passiert:

$$\begin{aligned}\hat\beta_{Bayesian} &= (X'X + B_0^{-1})^{-1}(B_0^{-1} b_0 + X'X \hat\beta) \\&= (X'X + B_0^{-1})^{-1}(B_0^{-1} b_0 + X'X (X'X)^{-1}X'y) \\&= (X'X + B_0^{-1})^{-1}(B_0^{-1} b_0 + X'y)\end{aligned}$$

In dieser Form wird die Analogie deutlich klarer:

• $\lambda I$ entspricht $B_0^{-1}$, der Präzisionsmatrix. Da $I$ die Einheitsmatrix ist, geht der Ridge-Schätzer von keiner Kovarianz zwischen den Regressionskoeffizienten und einer konstanten Präzision über alle Koeffizienten aus (denn $\lambda$ ist ein Skalar)
• $\lambda \mu$ entspricht $B_0^{-1} b_0$, was Sinn ergibt, da $b_0$ der Mittelwert der Prior-Verteilung ist – dieser zieht den Schätzer zu sich hin, genau wie $\mu$ die Koeffizienten zu sich „zieht“. Die Stärke dieses Effekts hängt von der Unsicherheit ab, die durch $B_0^{-1}$ bzw. $\lambda I$ im Ridge-Schätzer bestimmt wird.

Das ist soweit nachvollziehbar, aber sehen wir uns nun an, wie sich das Verhalten des bayesschen Schätzers bei Änderung der Unsicherheit im Vergleich zum Ridge-Schätzer verhält. Wir verwenden dieselben Daten und dieselbe Modellspezifikation wie oben, setzen die Kovarianzmatrix $B_0$ gleich $\lambda I$ und variieren dann $\lambda$. Denk daran: Kleinere Werte von $\lambda$ bedeuten hier einen stärkeren Einfluss des Priors (weniger Unsicherheit), größere Werte machen den Prior weniger bedeutend.

Die obigen Plots bestätigen unser bisheriges Verständnis: Mit einem Prior-Mittelwert von null werden die Koeffizienten in Richtung null gezogen – wie bei der Ridge-Regression, wenn der Prior dominiert, d. h. wenn die Präzision hoch ist. Und wenn ein nichttrivialer Prior-Mittelwert gesetzt ist, nähern sich die Koeffizienten diesem Wert an, sobald die Präzision steigt. Soweit zu den Koeffizienten – aber wie sieht es mit der Performance aus? Werfen wir einen Blick darauf!

Leistungsvergleich

Abschließend vergleichen wir die Vorhersagegenauigkeit der beiden Modelle. Zwar könnten wir die Parameter im Modell als Hyperparameter behandeln und sie entsprechend tunen, das würde jedoch dem Zweck widersprechen, vorhandenes Vorwissen zu nutzen. Stattdessen wählen wir für beide Modelle eine feste Prior-Spezifikation und vergleichen dann die Performance auf einem Hold-Out-Datensatz (30 % der Daten). Während wir für das Ridge-Modell einfach $X\hat\beta$ als Prädiktor verwenden können, liefert das bayessche Modell eine vollständige posterior-prädiktive Verteilung, aus der wir Stichproben ziehen können, um Modellvorhersagen zu erzeugen. Zur Schätzung des Modells wurde das Paket brms verwendet.

Insgesamt schneiden beide Modelle ähnlich ab, auch wenn einzelne Fehlermaße jeweils eines der Modelle leicht bevorzugen. Betrachtet man diese Fehlerwerte, könnte man sicher die Modelle noch verbessern, z. B. durch die Wahl einer geeigneteren Wahrscheinlichkeitsverteilung für die Zielvariable. Schließlich können Preise nicht negativ sein, und dennoch liefern unsere Modelle negative Vorhersagen.

Zusammenfassung

In diesem Beitrag habe ich gezeigt, wie der Ridge-Schätzer mit dem konjugierten bayesschen linearen Modell vergleichbar ist. Du kennst nun den Zusammenhang zwischen den beiden Ansätzen und wie der bayessche Zugang eine leichter interpretierbare Möglichkeit der Regularisierung bietet. Normalerweise würde $\lambda$ als Strafterm verstanden, aber jetzt lässt er sich auch als Maß für die Unsicherheit im Prior deuten. Ebenso kann der Parametervektor $\mu$ im erweiterten Ridge-Modell als Vektor von Prior-Mittelwerten für die Modellparameter verstanden werden. Beim bayesschen Ansatz besteht außerdem die Möglichkeit, Expertenwissen in Form von Prior-Verteilungen in den Schätzprozess einzubringen. Dies regularisiert das Modell und erlaubt die Einbindung externer Informationen. Wenn du dich für den Code interessierst, schau gerne auf unserer GitHub-Seite vorbei!

Quellen

• Jackman S. 2009. Bayesian Analysis for the Social Sciences. West Sussex: Wiley.

‍