Hands-on 2.Termin: Erweiterte Regression

Für das heutige Beispiel bearbeiten wir Daten über das Wetter am Seattle Airport im Jahr 2023, die wir aus einem separaten File einlesen (in unserem Fall einem Text-File mit der Endung .csv).

Dementsprechend brauchen wir ein paar vorbereitende Schritte. Diese führen wir heute einmal durch, damit ersparen Sie sich diese in zukünftigen Einheiten:

Im BOKUlearn-Kurs zur Lehrveranstaltung finden Sie im Abschnitt 'Unterlagen' sowohl die Hands-on-Unterlagen als auch die Datendateien.

Erstellen Sie auf Ihren PC einen Ordner, z.B. Statistik Vertiefung Hands-on. Sie können den Ordner beliebig nennen, Sie müssen sich nur merken, welchen Ordner Sie verwenden.

In diesem Ordner erstellen Sie einen Unterordner mit dem Namen Data. Für diesen Ordner verwenden Sie bitte den exakt gleichen Namen inklusive Klein- und Großschreibweise. Diesen verwenden wir nämlich später in Befehlen, und wenn Sie Ihn ebenso wie wir Data nennen, ersparen Sie sich in späteren Schritten das Abändern der von uns gezeigten Befehle.

Laden Sie die Hands-on-Unterlagen aus BOKUlearn herunter, und speichern Sie sie im von Ihnen gewählten Ordner.
Extrahieren Sie den zip-Ordner mit den Daten-Dateien im Unterordner Data.

Wählen Sie den Menüpunkt: File > New Project

Da Sie den Ordner (=directory) im vorigen Schritt schon erstellt und mit Dateien befüllt haben, wählen Sie nun: 'Existing Directory':

Klicken Sie nun auf 'Browse' und navigieren Sie zu Ihrem Ordner (also in unserem Fall Statistik Vertiefung Hands-on), und wählen Sie Ihn aus. (Achtung: Nicht den Unterordner Data auswählen!)

Klicken Sie nun auf 'Create Project', und warten Sie, bis das Project erstellt ist.

Sie sollten nun rechts unten in Rstudio im File-Browser die Hands-on-Unterlagen und einen Unterordner Data sehen:

Weiters sollte im rechten oberen Eck des Rstudio-Fensters das R-Projekt-Symbol mit dem Namen Ihres Ordners (für uns: Statistik Vertiefung Hands-on) aufscheinen:

Wenn diese beiden Dinge der Fall sind, haben Sie Ihr R-Projekt erfolgreich erstellt! Beim nächsten Mal können Sie nun einfach in Ihrem Ordner auf das File Statistik Vertiefung Hands-on.Rproj klicken, um Rstudio zu starten, und Zugriff auf die Daten im Data-Unterordner zu haben.

Nun öffnen wir ein neues R-Skript, um in die tatsächliche Analyse zu starten:

Öffnen Sie nun aus den Hands-on-Unterlagen die PDF-Datei Ex_02_GAM_.pdf zu den generalisierten Modellen, und befolgen Sie die Angaben. Die wichtigsten Punkte fassen wir hier noch einmal zusammen:

Lesen Sie die Daten ein:

rain_data<-read.csv("Data/RainSeattle2023_short.csv", 
    header=TRUE, na.strings=c(""), sep=";")

Wenn Sie eine Fehlermeldung bekommen, dass die Datei nicht gefunden wurde, achten Sie bitte auf die Pfadangabe im Befehl. Ist der Ordner Data tatsächlich ein Unterordner? Sind die Dateien darin entpackt? Falls Sie bereits in dem Ordner sind, in dem sich die Datendateien wie z.B. RainSeattle2023_short.csv befinden, müssen Sie das in der Pfadangabe berücksichtigen und das Data/ weglassen.

Wenn Sie die Datei richtig eingelesen haben, können Sie sich anschauen, was denn im File gespeichert ist, z.B. mit folgenden Befehlen:

head(rain_data) # zeigt die ersten 6 Zeilen der Datei an
str(rain_data) # zeigt die Datenstruktur an (ZeilenXSpalten der Datei+spaltenweise Infos)
summary(rain_data) # zeigt zusammenfassende Infos zu jeder Variable
View(rain_data) # öffnet die Datei in einer Tabellenkalkulationsprogramm-ähnlichen Ansicht

Die Anmerkungen, die hinter dem # stehen, sind Kommentare, die nicht ausgeführt werden, sondern Notizen für menschliche Leser*innen sind. Sie können sie auch weglassen, oder sich selbst andere relevante Notizen machen.

Beachten Sie, dass Sie mit 'Run' immer nur die Zeile ausführen, in der der Cursor gerade steht.

Falls Sie mehrere Befehle auf einmal ausführen möchten, markieren Sie alle entsprechenden Zeilen, und klicken dann auf 'Run'.

Wie wir sehen können, handelt es sich um eine Datei, die Wetterinformationen für alle Tage des Jahres 2023 am Seattle Tacoma Airport enthält.
Im Speziellen haben wir Infos zur Temperatur, der Windgeschwindigkeit und -richtung, zum Schneevorkommen, und ob es an dem Tag geregnet hat oder nicht (codiert als 0/1, wobei "1" für "es hat geregnet" steht).

Mit folgendem Befehl können wir uns die Beziehung zwischen Regen Ja/Nein und anderen Variablen anschauen, z.B. mit Windgeschwindigkeit:

plot(rain_data$WSF5, rain_data$Rain)

Wie wir sehen können, scheint es auf den ersten Blick keinen Zusammenhang zwischen Windgeschwindigkeit und Regen Ja/Nein zu geben:

Die rote Linie stellt die gedachte lineare Regressionslinie dar. Alternativ können wir uns aber auch überlegen, dass wir die Regenwahrscheinlichkeit abhängig von der Windgeschwindigkeit berechnen könnten. Dies ist durch die blaue Linie dargestellt. Relevante Unterschiede sind, dass wir einen nichtlinearen Anstieg berechnen können, und –ganz wichtig– nur Werte bekommen die zwischen 0 und 1 liegen (also Werte, die für eine Wahrscheinlichkeit auch tatsächlich möglich sind).

Diese Berechnung der Regenwahrscheinlichkeit ist mithilfe eines Generalisierten Linearen Modells (GLMs) möglich, bei welchem zuerst der lineare Prädiktor (rote Linie) berechnet wird, und dieser dann mithilfe einer weiteren Funktion (in unserem Fall der Exponentialfunktion) in die tatsächlich mögliche Form transformiert wird.

Da dieser Transformationsschritt sowohl die interne Berechnung als auch die Interpretation erschwert, testen wir die Treffsicherheit unseres Modells mithilfe einem Trainings- und Testdatensatz.

Bevor Sie folgenden Befehle ausführen, müssen Sie das Package "caTools" einmalig installieren (mithilfe install.packages('caTools'). Pakete müssen auf Ihrem PC nur einmal installiert werden. Damit Sie die Funktionen im Package nutzen können, müssen Sie es nach jedem Start von RStudio wieder laden (das geschieht mit dem Befehl library).

 library(caTools)
 set.seed(88)
 split <- sample.split(rain_data$Rain,SplitRatio = 0.8)
 train_df <- subset(rain_data,split==TRUE)
 test_df <- subset(rain_data,split==FALSE)

Mit dem obenstehenden Code-Chunk haben wir unseren Datensatz nun in 2 Teile geteilt: An 80% der Daten passen wir unser Modell an, und an den restlichen 20% vergleichen wir danach, wie häufig die Vorhersagen Regen Ja/Nein tatsächlich zugetroffen haben.
Die $365 \times 0.8 = 292$ Tage, die in unserem Trainingsdatensatz landen, werden dabei zufällig ausgewählt.

Auf dem Trainingsdatensatz passen wir nun unser Modell aus allen verfügbaren erklärenden Variablen an:

glm1 <- glm(Rain ~ Season + SNOW + SNWD + TAVG + TMAX + TMIN
            + WDF5 + WSF5, data=train_df,
            family=binomial)
summary(glm1)

Signifikanzen von erklärenden Variablen lassen sich aus der Ausgabe des Modells wie gewohnt ablesen, die Interpretation der Schätzer ist nun ein bisschen schwieriger. Wie gewohnt zeigt die Ausgabe Achsenabstand und Steigungen der 'Regressionsgeraden' an; im eindimensionalen Beispiel von oben also die rote Linie, die den linearen Prädiktor zeigt. Wir transformieren unseren linearen Prädiktor aber wie erwähnt in die Wahrscheinlichkeitskurve von Regen, also im 1D-Beispiel zur blauen Kurve. Somit müssen wir für die Interpretation unsere Schätzer ebenso transformieren. In unserem Fall einer Binomialverteilung mit Logit-Link (d.h. die Transformation von linearem Prädiktor zu Wahrscheinlickeitskurve erfolgt durch logarithmieren) können wir unsere Schätzer durch Exponenzieren¹⁾ interpretieren. Für eine detailliertere Erklärung der Interpretation der Schätzer halten Sie sich bitte an die Vorlesungseinheiten, als Kurzzusammenfassung für den Fall der Binomialverteilung (mit Logit-Link) reichen folgende Daumenregeln:

• Wenn sich eine Variable mit einem negativen Schätzer um 1 erhöht (bei Gleichbleiben aller anderen Variablen), sinkt die Regenchance²⁾.
• Wenn sich eine Variable mit einem positiven Schätzer um 1 erhöht (bei Gleichbleiben aller anderen Variablen), steigt die Regenchance³⁾.

Um zu sehen, wie stark die Regenchance genau steigt/sinkt, könnten wir uns $e^x$ für jeden Schätzer mit der Funktion exp() ausrechnen.

Zusätzlich lesen wir aus der Modell-summary noch ab, dass wir den Anteil unerklärter Abweichungen (=deviance) von 392.65 auf 266.43 reduzieren, und der AIC-Wert des Modells 288.43 beträgt.⁴⁾ Der AIC ist nicht allgemein interpretierbar, sondern nur in Vergleichen von geschachtelten Modellen, wie wir schon bei linearen Modellen mit der schrittweisen Variablenselektion (Funktion step) genutzt haben. Um somit trotzdem allgemeinere Aussagen über die Treffsicherheit unseres Modells machen zu können, nutzen wir unseren Testdatensatz:

preds <- predict(glm1, test_df, type='response')
predicted<-ifelse(preds>0.5,1,0)
table(predicted, test_df$Rain)

In der ersten Zeile berechnen wir die Regenwahrscheinlichkeiten, die unser Modell für den Testdatensatz voraussagt.
In der zweiten Zeile setzen wir einen harten Schnitt: Für alle Tage, an denen die vorausgesagte Regenwahrscheinlichkeit höher als 50% ist, sagen wir '1', also 'Regen Ja'. Je nach Anwendung könnten wir auch einen anderen Schwellenwert wählen.
In der letzten Zelle zeigen wir uns eine Vergleichstabelle an, wo unsere Voraussage Regen Ja/Nein mit dem tatsächlichen Regenvorkommen an den 'Testtagen' verglichen wird:

Wie wir sehen können, hat es an 42 Tagen nicht geregnet, und unser Modell hat auch Regen Nein vorausgesagt; und an 18 Tagen hat es tatsächlich geregnet, und unser Modell hat auch Regen Ja vorausgesagt. Falsche Voraussagen haben wir in 12+1 Fällen getroffen. Somit haben wir für $100 \times (42+18)/(42+18+1+12)= 82.2$% der 'Testtage' richtig vorausgesagt, ob es regnen wird oder nicht.

Im GLM könnten wir analog wie zur linearen Regression auch Variablen selektieren um die Modellperformance zu optimieren. Weiters könnten wir auch für zukünftige Tage, von denen wir die verwendeten erklärenden Variablen (z.B. Temperatur und Windgeschwindigkeit) wissen, die Regenwahrscheinlichkeit voraussagen. Hierfür verweise ich Sie aber auf die Hands-On-Angabe Ex_02_GAM_.pdf.

Als zweites Beispiel zu Erweiterungen schauen wir uns Generalisierte Additive Modelle (GAMs) an. Diese wurden entwickelt, um Nichtlinearitäten in den erklärenden Variablen berücksichtigen zu können (im Gegensatz dazu ist es bei den GLMs bis jetzt um komplexere Verteilungen in der erklärten Variable gegangen).⁵⁾

Als Beispieldaten für GAMs verwenden wir wieder Daten, die in einem R-package abgespeichert sind:

#install.packages(faraway) #entkommentieren, und führen Sie diese Zeile aus, wenn Sie das Paket noch nicht installiert haben
 data(ozone, package='faraway')
 summary(ozone)

Da die Daten aus einem R-Paket sind, gibt es zu dem Datensatz wieder eine Hilfeseite, die mit help(ozone) aufgerufen werden kann. Wie wir aus der Hilfeseite auslesen können, enthält der Datensatz Ozonmessungen und Wetterinformationen (plus 'Tag im Jahr') zur Großregion Los Angeles im Jahr 1976. Alle Variablen sind numerisch.

Die Beziehung zwischen den vorhandenen Variablen und der Ozonkonzentration könnten wir uns für jede Variable einzeln mit plot(ozone$*Variablenname*, ozone$O3) ausgeben⁶⁾, zur Einfachheit hier der Plot von Ozonkonzentration gegen alle anderen Variablen:

Wie wir sehen können, scheint zwischen Ozonkonzentration und einigen erklärenden Variablen ein annähernd linearer Zusammenhang zu bestehen. Bei anderen Variablen ist zwar ein Muster in den Datenpunkten zu erkennen, dieses ist aber nicht linear (z.B. bei doy, day of the year, also Tag im Jahr).

Im folgenden Plot ist in der roten Linie ein linearer Zusammenhang zwischen doy und O3 gezeichnet, und in der blauen Linie ein nichtlinear geschätzter Zusammenhang:

(Dass sich für die Ozonkonzentration ein Zusammenhang dieser Art mit Tag im Jahr ergibt, ist nicht erstaunlich, schließlich ist diese bekanntlich in den Sommermonaten am höchsten.) Würden wir also nur lineare Zusammenhänge zwischen erklärter und erklärenden Variablen erlauben, würde unser Modell also keinen Zusammenhang zwischen Ozonkonzentration und Tag im Jahr erkennen, obwohl diese im Sommer tatsächlich viel höher ist als im Winter.

Aus diesem Grund passen wir folgendes GAM an, das uns erlaubt, die erklärenden Variablen nichtlinear zu modellieren:

library(mgcv)
gam1 <- gam(O3 ~ s(vh) + s(wind) + s(humidity) + s(temp) +
              s(ibh) + s(dpg) + s(ibt) + s(vis) + s(doy),
            data=ozone)
summary(gam1)

Die Funktion gam findet sich im Paket mgcv. Sollten Sie dieses nicht auf Ihrem PC haben, müssen Sie es wie immer installieren.

Die summary des GAM sieht folgendermaßen aus:

Wir sehen also eine Aufsplittung der Ausgabeblöcke zu den Koeffizienten in Parametric coefficients (die linear geschätzt werden) und smooth terms (die geglättet geschätzt werden). Da wir in unserem Funktionsaufruf über alle erklärenden Variablen die glättende Funktion s() gelegt haben, ist der Block für die Parametric coefficients bis auf den Achsenabstand (Intercept) leer. Dieser ist laut Ausgabe signifikant von 0 verschieden.

Nun zu dem neuen Block der smooth terms:
Dieser enthält in der Spalte edf einen Schätzer für die Komplexität der angepassten Kurve, wobei sie ungefähr mit der Potenzierung der Kurve gleichgesetzt werden kann. s(vh) mit edf=1 wird also als linear geschätzt, s(doy) mit edf$\approx4$ wird also ungefähr als eine Kurve 4. Grades, d.h. stark nichtlinear geschätzt.
Die restlichen Spalten im Block enthalten wie gewohnt Infos, die zum Testen der Signifikanz nötig sind, und schließlich in der letzten Spalte das Signifikanzniveau.

Unser Modell sagt uns also, dass es für die Variablen vh und wind unnötig ist, diese als nichtlinear zu schätzen. Dementsprechend passen wir unser Modell an, indem wir um diese beiden die glättende Funktion s() entfernen:

gam2 <- gam(O3 ~ vh + wind + s(humidity) + s(temp) +
              s(ibh) + s(dpg) + s(ibt) + s(vis) + s(doy),
            data=ozone)
summary(gam2)

Nun sind diese beiden Variablen im Block zu den Parametric Coefficients gelandet, und können wie in einer normalen linearen Regression interpretiert werden. (Also: beide sind signifikant von 0 unterschiedlich; und wenn sich z.B. die Windgeschwindigkeit um eine Einheit erhöht, sinkt die Ozonkonzentration laut unserem Modell um 0.27ppm.)

Die restlichen Terme werden weiterhin als nichtlinear geschätzt, mit Kurven zwischen knapp 2. bis 5. Grades. Feuchtigkeit, Temperatur, dgp (Daggett-Luftdruckgradient), vis (Sicht) und doy (Tag im Jahr) haben signifikanten Einfluss auf die Ozonkonzentration. Der tatsächliche Einfluss der nichtlinear geschätzten Variablen lässt sich am einfachsten grafisch interpretieren:

plot(gam2, pages=1, scale=0)
#pages=1 zeigt alle Grafiken in einer Seite
#scale=0 erlaubt, dass unterschiedliche Achsenlimits für die verschiedenen Plots verwendet werden

Beispielsweise für Temperatur (temp) kann ich sehen, dass bei niedrigen Temperaturwerten die Ozonkonzentration zuerst leicht sinkt bzw. eher gleich bleibt, und dan ab einer gewissen Temperatur (knapp 60°F $\approx$ 15°C) stark ansteigt. Für Tag im Jahr (doy) ist die Ozonkonzentration um den 150. Tag im Jahr am höchsten.
Ein bisschen anders schaut das Bild bei ibt (Inversion Base Temperature) aus: die Konfidenzbänder für die geschätzte Kurve gehen weit auseinander, R ist sich also in großen Teilen des Datenbereichs sehr unsicher über die Einflüsse. Weiter kann man auch im Nullpunkt der y-Achse eine gedankliche Gerade durch die Konfidenzbänder legen: dies zeigt grafisch nochmal die Tatsache, dass der Einfluss von ibt nicht signifikant von 0 verschieden ist.
Wenn wir zurück zu summary(gam2) schauen, sehen wir, dass auch der Einfluss von ibh (Inversion Base Height) nicht signifikant von 0 verschieden ist. Auch das zeigt sich in der Grafik, da es auch hier möglich ist, eine gedankliche Linie durch den Nullpunkt zu legen, die innerhalb der Konfidenzbänder liegt.

Analog dazu können auch die Einflüsse der restlichen Variablen interpretiert werden.

Auch im Fall der GAMs ist es natürlich möglich, Faktoren einzubauen, Variablen zu selektieren, das Modell mithilfe Trainings- und Testdatensätzen zu optimieren, etc. Für das verweise ich auf Ex_02_GAM_.pdf.

Nur noch der versprochene Hinweis für HÜ2:
Analog zum Ozonbeispiel, wo wir mit Ozongehalt eine numerische Variable erklärt haben, könnten Sie auch für das erste Beispiel zur Regenwahrscheinlichkeit (wo Regen Ja/Nein binär codiert war) annehmen, dass erklärende Variablen einen nichtlinearen Einfluss haben, z.B. mit folgendem Befehl:

gam_rain <- gam(Rain ~ s(TAVG) + s(TMAX) + s(TMIN) + s(WDF5) + s(WSF5),
                data=rain_data, family=binomial)

Auch die GAM-Funktion kann also mithilfe des richtigen family-Arguments GLMs schätzen.