In Anlehnung an den Blogartikel, in welchem wir die bayessche Zeitreihenregression anhand bestehender MC FLO Formeln aufgezeigt haben, wollen wir dies hier mit der logistischen Regression fortsetzen. Vorab sei aber gesagt, dass – falls wir genügend Zuspruch bekommen sollten – die Berechnungsschritte direkt in MC FLO einzubauen gedenken.
Die Regression versucht einen numerischen Zusammenhang zwischen einer oder mehreren erklärenden Variablen und der gesuchten, erklärten Variable herzustellen, etwa der Grösse und Anzahl Badezimmer und dem Marktwert der Wohnung. Im einfachsten Fall kann dieser Zusammenhang linear wiedergegeben werden. Die gesuchte Grösse ist somit als Linearkombination der erklärenden Variablen darstellbar, wenn die Residuen (Differenz zwischen echter Grösse und berechneter Grösse anhand der Trainingsmenge) normalverteilt und die Streuung der Punkte um die Gerade in vertikaler Richtung konstant sind.
Die logistische Regression ist dann von Relevanz, wenn die gesuchte Grösse nur die Zahl 0 (trifft nicht zu) oder 1 (trifft zu) einnehmen kann, etwa beim Eintreten einer Insolvenz, welche anhand verschiedener Treiber – etwa der Eigenkapitalquote und Bestand der Lagervorräte – erklärt werden soll.
Um einen Vergleich mit R und dessen verwendeten Algorithmen herzustellen, werden wir das Wetterbeispiel aus https://www.bayesrulesbook.com/chapter-13 heranziehen. Dabei werden wir den klassischen frequentistischen als auch den bayesschen Ansatz heranziehen und beide mit MC FLO und Microsoft Excel nachstellen.
In der folgenden Excel-Datei haben wir alle Daten und Berechnungsschritte abgelegt, so dass diese nach Anstossen einer Monte-Carlo Simulation mit mindestens 1‘000 Iterationen (Beobachtungen) ausgewertet werden können.
Das Beispiel ist einfach aufgebaut: Es enthält 1‘000 Tagesdaten zu Regen (1 = es hat geregnet, 0 = es hat nicht geregnet, raintomorrow) und der Luftfeuchtigkeit in % um 9h des Vortages [humidity9am]. Die relevanten Fragen lauten: Gibt es einen Zusammenhang zwischen der Luftfeuchtigkeit um 9h und dem Aufkommen von Regen am Folgetag? Wie lautet die Regenprognose bei einer angenommenen Luftfeuchtigkeit von 50% um 9h des Vortags? Wie gut ist das Prognosemodell?
Wir werden analog der mit MC FLO mitgelieferten Beispiele vorgehen: Die Daten werden in eine Trainingsmenge und in eine Testmenge unterteilt. Die ersten 700 Datensätze (flexibel anpassbar, alle darauf verweisende Formeln werden automatisch angepasst) werden als Trainingsmenge zur Bestimmung der gesuchten Parameter herangezogen (Zelle AK22), bei den verbleibenden 300 Datensätzen werden die Parameter verprobt.
Die einfache lineare Regression hilft zur Bestimmung der gesuchten Parameter und zur Beantwortung der obigen Fragestellungen nicht weiter: Schon allein graphisch kann kein direkter linearer Zusammenhang identifiziert werden. So werden Werte von unter 0 und über 1 vorhergesagt. Formal sind aber nur Werte von 1 oder 0 möglich.
Um die Voraussetzung an eine lineare Regression zu erfüllen, ist die gesuchte Grösse zu transformieren. Die Transformation erfolgt über eine Logik Funktion, welche eine lineare Analyse erlaubt (Spalte G, Log(odds)). Die Logik-Funktion greift auf das Verhältnis der Regenwahrscheinlichkeit zu dessen Kehrwert zu.
Kernelement der logistischen Regression ist die Bestimmung von Wahrscheinlichkeiten, hier ob es am morgigen Tag auf Basis der heutigen Luftfeuchtigkeit um 9h regnen wird oder nicht. Würde das Regressionsmodell (über eine abermalige Transformation, siehe Spalten mit Bezeichung „probability“) eine Regenwahrscheinlichkeit von gleich oder über 50% ausweisen, wäre als Resultat eine 1 (= Regen) auszugeben, womit die ursprüngliche Problemstellung aufgegriffen wäre. Hingegen ist vorab vom Prüfer / Modellierer festzulegen, ab welcher Wahrscheinlichkeit eine 1 oder eine 0 zu prognostizieren ist. Ist das Tragen eines Regenschirms zumutbar, das Nasswerden hingegen nicht, könnte auch bei einer Regenwahrscheinlichkeit von ab 30% eine positive Prognose resultieren.
Die Regressionsgleichung ist bei einem Regressor (hier Luftfeuchtigkeit um 9h in %) mit y = a + b*x noch überschaubar; zu bestimmen sind nur der Achsenabschnitt a und die Steigung b. In unserem Fall sind wie bereits gesagt Wahrscheinlichkeiten in Übereinkunft mit dem binären Resultat (0 oder 1) in Einkunft zu bringen. Gesucht ist somit die Kurve, bei der die höchsten Treffer (0 oder 1) im Vergleich zum Trainingsdatensatz resultieren, womit ein Optimierungsproblem vorliegt.
Als erstes werden wir den klassischen frequentistischen Ansatz beiziehen, welcher einen „optimalen“ Wert bestimmt und dabei den Solver von Excel zur Bestimmung der Regressionsparameter a und b verwenden. Um eine konvergente Lösung, welche mindestens eine bessere Trefferquote als immer Regen = 0 vorherzusagen garantiert (hierbei wären 575 Treffer aus der Trainingsmenge maximal möglich, was einer Trefferquote von 79.5% entspricht; Zellen AK8:AK9), haben wir einige Nebenbedingungen eingebaut (Zellen AK38ff; Achtung: Wenn Sie das Beispiel mit dem Solver nachstellen, kann es sein, dass keine Lösung gefunden wird. Lassen Sie dann ggfs. mehrmals den Solver nach einer geeigneten Lösung suchen).
In unserer Lösung wurden für den Parameter a (Intercept, beta_0) ein Wert von -3.36 und für b (Slope, beta_1) ein Wert von 0.03 ermittelt, dies auf Basis der Formeln Spalten T bis W. Daraus resultieren neu 577 Treffer. Im Vergleich zur naiven, gleichlautenden Prognose „es wird kein Regen geben“, ist durch das Optimierungsverfahren eine Verbesserung von 2 zusätzlichen Treffern erwirkt worden.
Da das Optimierungsverfahren auf Basis der Trainingsdaten angewandt wurde, ist das Modell als nächster Schritt auf die Testdaten anzuwenden. Hierbei ist die Konfusionsmatrix eine Hilfestellung zur Interpretation.
237 von 239 (237+2) Daten wurden korrekt als „kein Regen“ prognostiziert. Die Spezifität liegt somit bei 99.2% (=237 / (237 + 2)). Hingegen wurden von 61 Tagen mit Regen nur 3 Tage korrekt identifiziert. Die Spezifität liegt daher bei bescheidenen 4.9% (= 3 / (58 + 3)). Die gesamthafte Treffergenauigkeit liegt bei 80% (=(237 + 3) / 300). Gegenüber dem Trainingsdatensatz wurde in den Testdaten eine Erhöhung der Genauigkeit von 0.4% gemessen.
Wird R als Referenz genommen und dabei keine zusätzlichen Nebenbedingungen aufgestellt, resultiert für Parameter a (beta_0) ein Wert von -4.28 und für b (beta_1) ein Wert von 0.04.
Im Gegensatz zur Solverlösung wird immer „kein Regen“ prognostiziert. Offensichtlich hat R keine bessere Lösung gefunden und verharrt wohl auf einem lokalen Optimum.
In der Bayesschen Statistik ist grundsätzlich ein anderer Weg einzuschlagen. Als erstes ist ein Vorwissen über die möglichen Realsierungen von a und b zum Ausdruck zu bringen. Dieses Vorwissen kann rein subjektiv aufgestellt, anhand eines Treibermodells aufbereitet oder aufgrund vergangener Daten hergeleitet sein. In unserem Fall geht die modellierende Person davon aus, dass a (beta_0) durch eine Normalverteilung mit Erwartungswert -1.4 und einer Standardabweichung von 2.1 ausgedrückt werden kann.
Bei b (beta_1) liegt der Erwartungswert bei 0.07 und einer Standardabweichung von 0.035.
Wird von einer Luftfeuchtigkeit von 0% um 9h des Vortages ausgegangen, dann beträgt im Erwartungswert die Wahrscheinlichkeit, dass das es am nächsten Tag regnen wird:
Exp(-1.4 + 0.07 * 0) / (1 + Exp(-1.4 + 0.07 * 0) = 0.2, 20%
Für eine Luftfeuchtigkeit von 88% um 9h liegt die Wahrscheinlichkeit für Regen am nächsten Tag im Erwartungswert bei:
Exp(-1.4 + 0.07 * 88) / (1 + Exp(-1.4 + 0.07 * 88) = 0.99, 99%
Da das Vorwissen als Wahrscheinlichkeitsverteilung vorliegt, können verschiedene mögliche Ausprägungen anhand der durch die Monte-Carlo Simulation generierten Daten identifiziert, ausgewertet und bessere Lösungen als die gleichlautende Prognose „es wird kein Regen geben“ identifiziert werden.
Mittels Monte-Carlo Simulation konnte im Vergleich zur Anwendung des Solvers von Microsoft Excel eine bessere Lösung in der Iteration 588 gefunden werden; neu werden aus den Trainingsdaten anhand der Parameter a (beta_0) = -6.87 und b (beta_1) = 0.07 insgesamt 580 Treffer erzielt. Auch hier kann anhand der Konfusionsmatrix eine Aussage in Bezug zu den Testdaten getroffen werden (siehe Zellen AJ25ff.):
Die Spezifität beträgt neu 98.3% (=235 / (235 + 4)). Die Gesamtgenauigkeit liegt bei 80.3%, was einer Erhöhung um 0.8% gegenüber der Strategie „immer Regen = 0 voraussagen“ entspricht.
Da mögliche Daten aus dem Vorwissen generiert wurden, können glaubhafte Intervalle zu den Parameterausprägungen getroffen werden (siehe Spalten BIff.)
Wir sind uns – bevor (neue) Daten ausgewertet wurden – zu 95% sicher, dass der Parameter a zwischen -5.35 und 2.77 und Parameter b zwischen 0.00 und 0.14 liegt. Oder auch anders: Da die Werte im 95% glaubwürdigen Intervall in Bezug auf b durchwegs positiv sind, ist von einem positiven Zusammenhang zwischen Luftfeuchtigkeit und Regenwahrscheinlichkeit auszugehen. Wir können nun mögliche Regenwahrscheinlichkeiten anhand der simulierten Daten generieren (Spalten BW-BY). Aufgrund der Unsicherheit in Bezug auf die Parameter resultiert bei einer angenommenen Luftfeuchtigkeit von 50% eine Regenwahrscheinlichkeit zwischen 5.7% und 99.9% im 95%-glaubwürdigen Intervall, der Durchschnitt liegt bei 76%, also Regen = 1. Zugegebenermassen ist die Unsicherheit hier sehr gross, aber aufgrund des Modellsettings (wir nehmen an, bisher noch keine Daten beobachtet zu haben) nicht anders zu erwarten.
In einem zweiten Schritt werden in der bayesschen Analysen Daten gesammelt und analysiert. Auch hier nehmen wir die 700 Trainingsdaten. Während wir aber anhand des Solvers die «beste» Lösung als Punktwert ermittelt haben, werden vorliegend über eine Monte-Carlo Simulation mögliche optimale Ausprägungen von a und b, gegeben die gemessenen Daten, bestimmt. Dabei werden alle Ausprägungen von a und b beigezogen, deren Trefferwerte dem Maximum ( Zahl 580) entsprechen (siehe Spalte AT). Mittels Kerndichteschätzung wird der datengenerierende Prozess simuliert, dieser entspricht der Likelihood Funktion und stellt die möglichen Ausprägungen von a und b anhand der gemessenen Daten dar (siehe Zellen AY3ff. und BC3ff.).
Im Vergleich zum Vorwissen ist aus dem datengenerierenden Prozess ein Erwartungswert von -5.11 für a (beta_0) und ein Erwartungswert von 0.06 für b (beta_1) zu erwarten. Für ein Tag mit Luftfeuchtigkeit von 50% wird im Mittel über alle möglichen Ausprägungen von a und b eine Regenwahrscheinlichkeit von 19% prognostiziert, also Regen = 0.
Der Vorteil der Bayesschen, generativen Statistik liegt in der Möglichkeit, Vorwissen mit gemessenen Daten zu kombinieren, ähnlich wie das Lernen beim Menschen funktioniert. Diese Kombination («Posterior») ist in Zellen AY50ff. und BC50ff. vorgenommen worden, das Ergebnis davon analog der Likelihood Funktion in Zellen kopiert und hiervon sind neue Eingangsvariablen erstellt worden. Nach einer erneuten Simulation sind alle Berechnungen abgeschlossen und analog der Beschreibung zum Vorwissen können die entsprechenden Auswertungen vorgenommen werden.
Wir sind uns – nachdem neue Daten mit dem Vorwissen ausgewertet wurden – zu 95% sicher, dass der Parameter a zwischen -5.92 und -2.32 und der Parameter b zwischen 0.05 und 0.09 liegt. Oder auch anders: Da die Werte im 95% glaubwürdigen Intervall in Bezug auf b durchwegs positiv sind, ist von einem positiven Zusammenhang zwischen Luftfeuchtigkeit und Regenwahrscheinlichkeit auch nach Beobachtung der Daten auszugehen. Aufgrund der Unsicherheit in Bezug auf die Parameter resultiert bei einer angenommenen Luftfeuchtigkeit von 50% eine Regenwahrscheinlichkeit zwischen 5.0% und 73.9% im 95%-glaubwürdigen Intervall, mit einem Mittelwert von 30%. Der Mittelwert liegt zwischen dem Mittelwert des Vorwissens (mit 76%) und dem Mittelwert der Likelihood-Funktion (mit 19%). Zur Bestimmung, ob Regen eintreffen wird oder nicht, wird im Regelfall aber nicht auf den Durchschnitt der ermittelten Parameterwerte abgestützt, sondern auf den Wert, welcher die Übereinstimmung mit den Trainingsdaten maximiert (Maximum a-posteriori). In unseren Fall liegt a (beta_0) bei -5.77 und b (beta_1) bei 0.06.
Wir können nun mögliche Regenwahrscheinlichkeiten anhand der simulierten Daten generieren (Spalten CA-CC). Zugegebenermassen ist die Unsicherheit weiterhin sehr gross, da das glaubhafte Intervall von 5.0% bis 73.9% sowohl die Prognose Regen = 0 und Regen = 1 zulässt. Wird hingegen in Zelle CC3 eine Luftfeuchtigkeit von 34% angegeben, dann lautet die Prognose immer Regen = 0, da die obere Grenze des Intervalls mit 49.3% unter der 50% Schwelle für Regen = 1 liegt.
Im Trainingsdatensatz sind 16 Datensätze mit Luftfeuchtigkeit (9h) <= 34% vorhanden, davon ist ein Datensatz mit Regen = 1 enthalten. 15 Datensätze enthalten Regen = 0. Gerundet sind somit ca. 94% der Datensätze mit Regen = 0 hinterlegt. Wird das Prognoseintervall des Bayesschen Modells auf 96% ausgeweitet, wird der Fall Regen = 1 ebenfalls erfasst, womit prima vista die Prognosegüte als sachgerecht eingestuft werden kann.
P.S.: Wenn Sie das Modell nachstellen, werden Sie sehr wahrscheinlich andere Zahlen bekommen, da wir das Modell mit mehreren Simulationsdurchläufen aufgebaut haben.
Schreiben Sie uns, falls Ihnen die Beschreibung gefallen hat und/oder Sie die logistische Regression auch automatisiert in MC FLO, analog der linearen Regression, vorfinden möchten.
P.S.2: Wir haben das Ergebnis mit dem Naive-Bayes Klassifikator von MC FLO nachgestellt und kommen mit dem Kreuzvalidierungsverfahren (jeweils 300 unterschiedliche Daten in der Testdatenmenge) auf eine gleich hohe Genauigkeit wie bei der bayesschen, logistischen Regression. Es stellt sich somit die Frage, ob eine Klassifikation nicht ausreichen würde. Es kommt auf den Sachverhalt an. Möchten Sie eine reine Prognose erstellen, ohne den numerischen Zusammenhang der erklärenden und der erklärten Variablen in einer Gleichung sichtbar zu machen, dann reicht ein Klassifikationsverfahren. Regressionen beschreiben hingegen explizit einen Zusammenhang (so haben wir identifizieren können, dass zwischen Regen und Luftfeuchtigkeit ein positiver Zusammenhang besteht; je höher die Luftfeuchtigkeit, desto höher steigt die Wahrscheinlichkeit für Regen an), der primär zum Verständnis eines vorher definierten Modells beiträgt und erst sekundär als Prognoseinstrument eingesetzt wird.
P.S.3: Spielen Sie mit dem Modell und fügen Sie weitere Prädiktoren (etwa die Luftfeuchtigkeit um 15h des Vortags [«humidity3pm»]) hinzu, um den Zusammenhang besser zu verstehen und um die Vorhersage zu verbessern. Was meinen Sie, welche Prädikatoren sind sinnvoll anzuwenden?
Kommentar schreiben