Über unsMediaKontaktImpressum
Elena Jolkver 09. Mai 2023

Datenanalyse mit R

Jeder, der schon einmal eine Tabelle erstellt und vielleicht sogar die Daten in einem Diagramm dargestellt hat, hat bereits Datenanalyse betrieben. Im Folgenden werden wir sehen, dass sich die Skriptsprache R dafür besonders gut eignet. Wir werden dies am Beispiel der Vorhersage der Ozonbelastung in der Stadt erläutern.

 

Was bedeutet "Datenanalyse"?

Datenanalyse beinhaltet praktisch jede Handlung, die man mit seinen gesammelten Daten durchführt, mit Ausnahme des Speicherns. Es umfasst zunächst die langwierige und obligatorische Datenbereinigung, die für gewöhnlich mit einer Beschreibung der Daten einhergeht. Danach taucht der Analyst ein und sucht nach Zusammenhängen zwischen verschiedenen Variablen – ein wenig wie ein Meeresforscher, der die Beziehungen zwischen den Fischen erst noch erkennen muss. Meist geschieht dies mit einem höheren Ziel, denn hat man erst die Abhängigkeiten in den Daten verstanden, kann man Modelle aufstellen und mit ihnen Vorhersagen über noch nicht beobachtete Datenpunkte treffen. Außerdem sind Modelle hilfreich, weil sie uns die komplexen Abhängigkeiten in den Daten verständlich machen und Hypothesen über unsere Welt generieren lassen.

Allerdings lauern auf dem Weg der Modellierung zahlreiche Gefahren. Das Modell soll allgemeingültig und frei von subjektiven Entscheidungen sein. Außerdem bewegt man sich unentwegt zwischen dem Overfitting bzw. der Überanpassung und dem Bias, dem systematischen Fehler. Das wichtigste Werkzeug, das den Analysten durch den Dschungel führt, ist die Statistik, insbesondere das Testen von Hypothesen. Und hierbei hat sich die Skriptsprache R als ein sehr beliebtes, da vielseitiges und einfach zu erlernendes Hilfsmittel etabliert.

Warum R?

R ist eine Sprache und eine kostenlose Open-Source-Plattform für statistische Analysen und Visualisierungen. Es gibt mehrere Paketsammlungen, die kontinuierlich aktualisiert und weiterentwickelt werden: die zentrale Sammlung CRAN, die derzeit 19.310 Pakete listet, Bioconductor, das mit 2.183 Paketen auf Bioinformatik spezialisiert ist, sowie eine unbekannte Anzahl von Paketen auf Github. Es gibt Pakete, um hoch spezialisierte Analyseabläufe aufzubauen, die dem Analysten helfen, die Daten effizient einzulesen, sie vorzuverarbeiten, zu analysieren, zu visualisieren und anschließend zu speichern. Oft veröffentlichen Statistiker neue Algorithmen zur Datenanalyse zusammen mit einem entsprechenden R-Paket. Daher findet man in R durch das reichhaltige Paketangebot so ziemlich jede Funktionalität, die man bei der Datenanalyse braucht, einschließlich der neuesten Trends und Entwicklungen.

R hat eine Reihe von positiven Eigenschaften, die sowohl in experimenteller als auch in regulierter Umgebung von Vorteil sind:

  1. Open Source erlaubt die vollständige Nachvollziehbarkeit des Programms.
  2. Es gibt eine Reihe hervorragender Bibliotheken, um Daten in ein ordentliches Format zu bringen und exploratorische Datenanalyse zu betreiben, z. B. dplyr, data.table und tidyr.
  3. Es gibt eine Vielzahl von Bibliotheken, die auf Datenmanipulation und -analyse, einschließlich linearer Regression, Hypothesentests, Korrelation, Clustering und maschinelles Lernen spezialisiert sind.
  4. Es generiert publikationstaugliche Visualisierungen mit ggplot und ggpubr und erlaubt interaktive Visualisierungen mit plotly.
  5. HTML-Reports und sogar Bücher lassen sich mit RMarkdown einfach generieren.
  6. Mittels der Shiny-Umgebung lassen sich Web-Apps ohne spezielle Kenntnisse von html/css/javascript erstellen.
  7. Es kann mit anderen Sprachen (z. B. Python, SQL, C++) kombiniert werden und lässt sich mit Docker teilen. Außerdem ist es plattformunabhängig und läuft auf jedem der gängigen Betriebssysteme.
  8. R ist sehr benutzerfreundlich und auch für Anwender ohne Programmierkenntnisse sehr leicht zu erlernen.
  9. Oft umfasst die Dokumentation nicht nur die Beschreibung der Funktionen, sondern auch eine Erklärung der statistischen Methoden und worauf es bei der Analyse ankommt, ohne dass man sich dieses Wissen anderswo aneignen muss.

Die Liste der Nachteile ist hingegen übersichtlicher:

  1. Wenn die Analysen umfangreicher oder die Datensätze größer werden, zeigt sich R’s größter Nachteil: Es ist relativ langsam. Wenn die Performance ein wichtiges Kriterium ist und es um Millisekunden geht, sollte man zu einer anderen Sprache greifen, auch wenn sich viele Prozesse in R parallelisieren lassen. 
  2. R benötigt tendenziell mehr Speicherplatz.
  3. Die Vielzahl der Bibliotheken geht mit Unübersichtlichkeit einher. Es gibt sogar Pakete, die bei der Auswahl des richtigen Pakets unterstützen. Zwar haben sich gerade für die Datenanalyse ein paar Kernpakete etabliert, dennoch muss man sich hin und wieder über neue Entwicklungen informieren und z. B. R-bloggers verfolgen [1].

Von Besen und Waschmitteln – Erst wird geputzt

Meist beginnt die Datenanalyse damit, dass man irgendwoher Daten bekommt und je nach Fragestellung bestimmte Analysen durchführt. So stellt es sich die Analystin zumindest vor. In Wahrheit verbringt sie jedoch den Großteil ihrer Arbeitszeit damit, die Daten zusammenzutragen (z. B. durch Webscraping oder Datenbankabfragen) und sie "zu putzen". Nur zu oft enthält der Datensatz Lücken oder unplausible Werte. Auch die Verteilung der Daten spielt eine große Rolle dafür, welche statistischen Verfahren anwendbar sind. Eine Datenanalyse sollte man folglich niemals blind durchführen, sondern sich permanent vergewissern, dass man nicht in die Falle "Garbage In – Garbage Out" fällt.

Zur Demonstration bedienen wir uns eines Datensatzes zur Luftqualität in New York aus dem Jahr 1973. Der Datensatz wird von R bereitgestellt und man kann direkt darauf zugreifen, wenn man R installiert hat. Der Aufruf

head(airquality)

liefert die ersten sechs Zeilen des Datensatzes, während summary(airquality) eine kurze Zusammenfassung bereitstellt.

Die Parameter sind wie folgt definiert:

  • Ozone: Ozonmittelwerte [ppb] über 1300-1500 Stunden auf Roosevelt Island.
  • Solar.R: Sonnenstrahlung [4000-7700 Ångström] von 8-12 Uhr.
  • Wind: Durchschnittliche Windgeschwindigkeit [Meilen/Stunde] von 7-10 Uhr.
  • Temp: Maximale Tagestemperatur [Fahrenheit].

Wir erkennen gleich mehrere Probleme des Datensatzes. Zum einen ist der Monat als Integer angegeben. Dies könnte später zu Problemen führen, z. B. bei der Darstellung oder bei der Modellierung. Deswegen erzeugen wir eine neue Spalte, in der wir die Monatsnamen als Strings kodieren. Für die Umbenennung von Monatsnamen gibt es einen speziellen Befehl.

airquality$Monthname <- month.abb(airquality$Month)

airquality$Monthname <- factor(airquality$Monthname, levels = c("May", "Jun", "Jul", "Aug", "Sep") ) # dieser Aufruf bringt die Monatsnamen in chronologische Reihenfolge

Zum Anderen weisen sowohl die Ozon- als auch die Solarwerte Lücken auf. Bevor wir uns jedoch damit weiter beschäftigen, betrachten wir noch mögliche Ausreißer. Ausreißer können die gesamte Analyse zunichte machen, weshalb wir die Daten daraufhin unbedingt prüfen müssen. Eine Ausreißerdefinition besagt, dass es sich dabei um Punkte handelt, die unter Q(1/4)−3/2 IQR (Q: Quantile, IQR: Interquartilabstand) oder über Q(3/4)+3/2 IQR liegen. Gibt es solche Daten, können wir die Ausreißer unterschiedlich behandeln:

  • Kappen
  • Löschen
  • Imputieren
  • Binnen + Glätten

Wir sollten nur noch bedenken, dass wir eine mögliche saisonale Komponente im Datensatz haben, weshalb es sinnvoll sein könnte, die Ausreißer pro Monat zu betrachten:

for (mo in unique(airquality$Monthname)){
  print(mo)
  lower_bound <- quantile(airquality$Ozone[airquality$Monthname==mo], na.rm=T)[2] - 3/2 * IQR(airquality$Ozone[airquality$Monthname==mo], na.rm=T)
  upper_bound <- quantile(airquality$Ozone[airquality$Monthname==mo], na.rm=T)[4] + 3/2 * IQR(airquality$Ozone[airquality$Monthname==mo], na.rm=T)
  print(airquality$Ozone[airquality$Monthname==mo][which(airquality$Ozone[airquality$Monthname==mo] < lower_bound | airquality$Ozone[airquality$Monthname==mo] > upper_bound)])
}

Wir können die Daten auch grafisch darstellen. Sehr schöne Graphen erzeugt man mit der ggplot2-Bibliothek. Diese kann man entweder einzeln laden oder als Teilpaket von "tidyverse", das noch einige andere hilfreiche Bibliotheken zur Verfügung stellt. Um die Bibliothek nutzen zu können, muss man sie zunächst mit install.packages("tidyverse") einmalig installieren und bei Bedarf aufrufen.

library(tidyverse)
airquality %>% # Auswahl des Datensatzes
  select(Ozone, Solar.R, Wind, Temp, Monthname) %>% # Auswahl Spalten
  pivot_longer(c(Ozone, Solar.R, Wind, Temp), names_to = "observation", values_to = "value") %>% # wir brauchen die Daten im long-Format
    ggplot(., aes(x = observation, y = value)) + # definiere, was geplottet werden soll
      geom_boxplot(aes(fill=Monthname)) + # definiere Plottyp
      theme_bw() + # Aufhübschen des Plots
      ggtitle("Übersicht über 'airquality'") # Plotüberschrift

Wir erkennen, dass es auch in anderen Parametern Ausreißer gibt. Diese ersetzen wir nun durch den Median des Parameters des jeweiligen Monats.

cols_to_clean <- c("Ozone", "Solar.R", "Wind", "Temp") #generiere Vektor mit Spaltennamen, über die iteriert wird
for (col in cols_to_clean){ #Schleife über Spalten
  for (mo in unique(airquality$Monthname)){ #Schleife über Monate
    lower_bound <- quantile(airquality[,col][airquality$Monthname==mo], na.rm=T)[2] - 3/2 * IQR(airquality[,col][airquality$Monthname==mo], na.rm=T)
    upper_bound <- quantile(airquality[,col][airquality$Monthname==mo], na.rm=T)[4] + 3/2 * IQR(airquality[,col][airquality$Monthname==mo], na.rm=T)
# ersetze Werte, die über/unter den Grenzwerten liegen, durch den Median des Parameters   
airquality[,col][airquality$Monthname==mo][which(airquality[,col][airquality$Monthname==mo] < lower_bound | airquality[,col][airquality$Monthname==mo] > upper_bound)] <- median(airquality[,col][airquality$Monthname==mo], na.rm=T)
  }
}

Jetzt betrachten wir die lückenhaften Werte. Es gibt verschiedene Ansätze, mit fehlenden Werten umzugehen. Zu unterscheiden sind einerseits Fälle, bei denen fehlende Werte auf Probleme in der Datenerhebung deuten und andererseits sog. "informed missingness", bei denen die Abwesenheit eines Werts eine Aussage tragen könnte (z. B. fehlende Antworten bei Umfragen). Da wir davon ausgehen können, dass das Fehlen von Ozon- bzw. Temperaturwerten das Verständnis dieses Datensatzes enorm einschränkt, haben wir folgende Möglichkeiten, damit umzugehen:

  • Löschen des Datensatzes
  • Ersetzen durch globale Konstante ("n.d")
  • Ersetzen durch Lageparameter (Mittelwert, Median)
  • Ersetzen durch wahrscheinlichen Wert, den man z. B. mittels missForest imputiert

Wir ersetzen auch die fehlenden Werte, so wie die Ausreißer, durch den Median. In einem komplexeren Datensatz, der mehr Parameter enthält, könnte man stattdessen versuchen, wahrscheinliche Werte beispielsweise mittels missForest zu schätzen.

#Ersetze NA's
airquality$Ozone[is.na(airquality$Ozone)] <- median(airquality$Ozone, na.rm=T)
airquality$Solar.R[is.na(airquality$Solar.R)] <- median(airquality$Solar.R, na.rm=T)

Jetzt haben wir einen ersten Überblick über die Daten gewonnen und den "gröbsten Dreck" entfernt. Ehe wir mit weiteren Analysen fortfahren, müssen wir uns noch einen Überblick über die Verteilung der Werte verschaffen. Einige statistische Methoden, wie z. B. der t-Test, verlangen, dass die Daten normalverteilt sind, andernfalls hat der p-Wert des Tests keine Aussagekraft. Einen schnellen Überblick verschafft uns der Dichte- sowie Korrelationsplot, der mittels ggpairs() aus der GGally-Bibliothek erstellt werden kann.

# install.packages("GGally")
library(GGally)

ggpairs(airquality,                 # Data frame
        columns = 1:4,        # Auswahl der Spalten
        aes(color = Monthname,  # Farbe nach Monat
            alpha = 0.5))     # Transparenz

Die Schiefe der Verteilung legt nahe, sie zu logarithmieren, zu quadrieren oder andere Transformationen anzuwenden, woraufhin man die Verteilung erneut prüfen sollte. Wir werden für die Modellierung die Daten weiterverwenden, wie sie sind.

Alles eine Frage der Klasse

Ozon ist ein Bestandteil des Luftqualitätsindex AQI[2]. Dieser setzt sich aus Ozon, Partikelverschmutzung, CO, SO2 und NO2 zusammen. Übersteigt der AQI bestimmte Grenzwerte, reagieren erst empfindliche Menschen auf die Belastung, bei höheren Werten bekommen alle Menschen Atemprobleme und sollten körperlich anstrengende Tätigkeiten im Freien einstellen. Wenn man die Luftqualität vorhersagen kann, kann man entsprechend warnen. Wir werden versuchen, eine grobe Klassifizierung vorzunehmen, indem wir die AQI-Klassenzuweisung zunächst auf zwei Klassen (normal, belastend) vergröbern [3].

airquality$class <- ifelse(airquality$Ozone <= 70, "normal", "belastend") # Generierung der Klasse
table(airquality$class, exclude=NULL) # zähle Einträge pro Klasse

Der Aufruf table() zählt die Einträge pro Klasse, wobei es in dem betrachteten Zeitraum an 18 Tagen zu einer Grenzwertüberschreitung kam, die man als "belastend" einstufen würde, während 135 im normalen Bereich blieben.

Es gibt sehr viele verschiedene Methoden zur Klassifizierung. Hier sind nur ein paar Beispiele für verschiedene Algorithmen, die unterschiedliche Ansätze verfolgen:

  • Nach Ähnlichkeit: kNN
  • Nach Wahrscheinlichkeit: Logistische Regression und Naive Bayes
  • Nach Entfernung: Discriminant analysis und Support Vector Machine (SVM)
  • Nach Regeln: Entscheidungsbäume
  • Kombinationen: Bagging und Boosting

Eine nützliche Auflistung der Methoden finden Sie bei caret[4].

All diese Ansätze haben gemeinsam, dass sie zum sogenannten "Supervised Machine Learning" gehören. Dabei lernt der Algorithmus zunächst anhand von Beispielen mit bekannter (deswegen supervised) Klasse die zugrundeliegenden Regeln oder Verteilungen und kann dann einen neuen Datenpunkt von unbekannter Klasse einer Klasse zuweisen.

Leider kann man von vornherein überhaupt nicht vorhersagen, welche Methode das beste Modell liefert, denn das hängt von der Datenverteilung ab. Deswegen sollte man verschiedene Ansätze ausprobieren und das eindeutig beste Modell wählen. Sollte die Auswahl schwer fallen, kann man auch ein Ensemble aus verschiedenen Modellen verwenden, indem man beispielsweise die Wahrscheinlichkeit der Klassenzugehörigkeit aus z. B. dem Mittelwert der Wahrscheinlichkeiten der Einzelmodelle ermittelt.

Wir probieren, den Datensatz mittels einiger gängiger Algorithmen zu klassifizieren. Zunächst müssen wir entscheiden, welche Parameter als Prädiktoren ins Modell eingehen werden. Natürlich nehmen wir den numerischen Ozonwert heraus, da wir ihn für die Klassifizierung benutzt haben und keine Vorhersage für einen bereits bekannten Wert benötigen. Stattdessen wollen wir die Ozongefahrenklasse basierend auf den anderen Parametern ermitteln. Wahrscheinlich würden dafür die numerischen Parameter genügen, wir beziehen jedoch im ersten Schritt die Monate mit ein, zumal sie uns den Umgang mit kategorischen Parametern üben lassen. Die meisten Algorithmen können mit einer kategorischen Spalte nicht umgehen und brauchen sog. "hot-one" kodierte Daten, bei denen der Datensatz für jede Stufe eine neue Spalte mit den Werten (0, 1) bekommt, die die An- bzw. Abwesenheit der Stufe kodieren. Diese, sowie viele weitere nützliche Funktionen für das maschinelle Lernen, sind z. B. in dem Paket "caret" zusammengefasst.

# install.packages("caret")
library(caret)
dummy <- dummyVars(" ~ Monthname", data=airquality, fullRank = TRUE)
newdata <- data.frame(predict(dummy, newdata = airquality)) #generiere Hot-one Datensatz
df <- cbind(airquality[,c("class","Solar.R","Wind","Temp" )], newdata) # füge neue Spalten an den alten Datensatz mit numerischen Werten hinzu
head(df) # zeige erste Zeilen

Jetzt müssen wir den Datensatz aufteilen: Ein Teil wird benutzt, um den Algorithmus zu trainieren, der andere, um ihn anschließend zu testen. Wenn man Hyperparametertuning betreibt, sollte man zusätzlich noch einen Validierungsdatensatz zurückhalten. Da man beim Tuning Gefahr läuft, den Algorithmus zu "overfitten", braucht man einen Satz unabhängiger Daten, der die Performance des Modells nach dem Feinschliff des Hyperparametertunings ehrlich schätzt, denn das Modell wurde ja an den Validierungssatz angepasst. Oft bedient man sich hierbei eines Verhältnisses 65-25-10, man kann aber –  abhängig von der tatsächlichen Datensatzgröße – das Verhältnis anpassen. Dabei sollte man mehrere Test- und Trainingssätze generieren und die Modellqualität als Durchschnitt der Qualitäten aller Modelle betrachten. Damit verhindert man, dass durch die zufällige Probenauswahl das Modell zu gut oder zu schlecht abschneidet. Wir werden auf einen Testdatensatz verzichten, da die Anzahl "belastender" Datenpunkte so gering ist, dass der Test nur eine sehr geringe Aussagekraft haben wird. Dafür erhöhen wir die Anzahl der Daten im Trainingsdatensatz.

Auch hierfür bietet "caret" eine Funktion an:

testnumber = 5
set.seed(1234) # da die Auswahl zufällig erfolgt, setzen wir ein Seed, um immer wieder gleiche Ergebnisse zu erzielen
flds <- createDataPartition(df$class, times=testnumber,  p = 0.75, list = TRUE) # generiert 5x einen Trainingssatz

Nun können wir den Algorithmus trainieren, wobei wir gleich eine Schleife über verschiedene Algorithmen drehen. Dabei werden wir den Datensatz mehrfach in Train/Test aufteilen, die Datensätze skalieren, ein Modell trainieren und dieses testen. Die Ergebnisse werden in einer Tabelle für die Post-Modellierungsschritte aufbewahrt.

library(pROC) # Bibliothek für ROC-Kurvenanalyse
# Funktion, um Modelleigenschaften zu extrahieren
getModelQuality <- function(confMatrix){
  data.frame(accuracy=confMatrix$overall["Accuracy"],
             kappa=confMatrix$overall["Kappa"],
             precision=confMatrix$byClass["Precision"],
             sensitivity=confMatrix$byClass["Sensitivity"],
             specificity=confMatrix$byClass["Specificity"],
            F1=confMatrix$byClass["F1"])
}

#Festlegung von Trainingsparametern
control=trainControl(method="repeatedcv", number=10, repeats=3, classProbs = TRUE, summaryFunction = twoClassSummary)
Modellqualitaet <- data.frame() #Kontainer für Ergebnisse
Modellparameter <- data.frame() # Kontainer für Parameterbedeutung
Testergebnisse <- data.frame() # Kontainer für Testergebnisse
algos <- c("knn", "rf", "svmLinear2", "glmnet")
set.seed(1234)
for (algo in algos){ # Schleife über Algorithmen
  for (i in 1:(length(flds))){ # Schleife über alle crossvalidation sets
    print(i)
    # Train/Test split
    train_unscaled <- df[ flds[[i]], ]
    test_unscaled <- df[ -flds[[i]], ]
    # Normierung
    preProcValues <- preProcess(train_unscaled, method = c("center", "scale"))
    train <- predict(preProcValues, train_unscaled)
    train$class <- as.factor(train$class )
    test <- predict(preProcValues, test_unscaled)
    test$class <- as.factor(test$class )
    # Trainiere Modell
    model <- train(class~., data=train, method=algo, trControl=control)
    # Vorhersage mit Modell
    prediction <- predict(model, newdata = test)
    test$pred <- prediction
    test$model <- algo
    test$id <- rownames(test)
    Testergebnisse <- rbind(Testergebnisse, test)
    results_predict<-getModelQuality(confusionMatrix(data = test$pred, test$class, positive="belastend"))
    results_predict$model <- algo
    results_predict$test_number <- i
    # ROC
    prediction_prob <- predict(model, newdata = test,type="prob")
    ROC_curve <- roc(predictor=prediction_prob[,"belastend"],
             response=test$class,
             levels=rev(levels(test$class)))
    ROCAUC <- ROC_curve$auc[1] # extrahiere AUC Wert
    results_predict <- cbind(results_predict, ROCAUC)
    # Confusion matrix
    ConfMatrix <- confusionMatrix(data = prediction, test$class,  positive="belastend")
    confMatrixTable_1 <- t(as.data.frame(c(ConfMatrix$table[1],ConfMatrix$table[2],ConfMatrix$table[3],ConfMatrix$table[4])))
    colnames(confMatrixTable_1) <- c("belastend.belastend_vorhergesagt","belastend.normal_vorhergesagt","normal.belastend_vorhergesagt","normal.normal_vorhergesagt")
    results_predict <- cbind(results_predict,confMatrixTable_1)  
    # Abschließen der Schleife: Zusammenfügen der Ergebnisse zu Modellqualität
    Modellqualitaet <- rbind(Modellqualitaet,results_predict)
    # Export der Modellparameter
    Modellpar <- varImp(model)$importance
    Modellpar$Parameter <- row.names(Modellpar)
    Modellpar <- as.data.frame(Modellpar[,c(1, ncol(Modellpar))])
    colnames(Modellpar) <- c("relImportance", "Parameter")
    Modellpar$model <- algo
    Modellpar$test_number <- i
    Modellparameter <- rbind(Modellparameter, Modellpar)
  }
}

# Generiere Tabelle im long-Format für die anschließende Grafik
Modellqualitaet_long <- Modellqualitaet %>%
  select(model, accuracy, precision, sensitivity, specificity,ROCAUC,F1, kappa) %>%
  gather(key = parameter, value = measurement,
       accuracy, precision, sensitivity, specificity,ROCAUC,F1, kappa)

# generiere Grafik
ggplot(Modellqualitaet_long, aes(x= parameter, y= measurement, fill= model)) + #definiere, was geplottet wird
  geom_boxplot() + # definiere Plottyp
      theme_bw() + # Aufhübschen des Plots
      ggtitle("Qualität der Klassifikatoren") # Plotüberschrift

Wir sehen, dass kNN die höchste Sensitivität mit der geringsten Präzision und Spezifität aufweist, d. h. von den als "belastend" vorhergesagten Tagen lagen eigentlich einige im Normalbereich. Wir sehen auch eine starke Streuung der Werte, je nachdem, welchen Trainings-/Testdatensatz wir verwendet haben. Das unterstreicht noch einmal, wie wichtig es sein kann, die Modellierung mit mehreren Sätzen durchzuführen, wenn man verschiedene Modelle miteinander vergleichen möchte. Jetzt liegt es an der Fragestellung, welches der Modelle das geeignetere ist. Wenn wir keinen belastenden Tag auslassen wollen, dafür aber falsch-positive in Kauf nehmen und zu oft eine Warnung aussprechen, können wir das sensitivere kNN-Modell bevorzugen. Wollen wir hingegen nicht riskieren, dass unsere Warnungen irgendwann nicht mehr ernst genommen werden, da wir zu oft die Alarmglocke läuten, sollten wir ein spezifisches Modell wählen. Aber vielleicht irren die weniger sensitiven Modelle auch, wenn die Entscheidung grenzwertig war? Dafür ist es notwendig, die Modelle besser zu verstehen. Wir werden uns zunächst ansehen, welche Parameter das Modell als wichtig erachtet hat und wie die falschen Vorhersagen zustande kamen.

Die Suche nach der Nadel im Heuhaufen

Wir wollen dem Modell auf den Grund gehen. Dafür können wir uns zunächst einen Überblick verschaffen, welche Parameter von den Modellen als wichtig erachtet wurden. Dazu haben wir in der oberen Trainingsschleife die Tabelle "Modellparameter" generiert, die die Bedeutung der Parameter für die einzelnen Algorithmen festhielt und die wir uns jetzt ansehen werden.

library(knitr)
library(kableExtra)
Modellparameter %>%
  group_by(model, Parameter) %>% # analysiere pro Modell und Parameter
  summarize(meanImportance = round(mean(relImportance, na.rm=T), 1)) %>% #bilde Mittelwert und runde auf 1 Stelle
  pivot_wider(names_from = model, values_from = meanImportance) %>% # pivotisiere Tabelle
  arrange(desc(knn)) %>% #sortiere nach Spalte 'knn'
  kbl(caption = "Relative Bedeutung einzelner Parameter") %>% #Aufhübschen der Tabelle
  kable_classic(full_width = F, html_font = "Cambria")

Wir sehen, dass alle Modelle die Temperatur als wichtigsten Einflussfaktor auf den Ozonwert favorisieren. Es folgte Wind, während sich bei den nachfolgenden Parametern die Modelle unterschieden. Dabei fällt auf, dass kNN und SVM identische Parameterbedeutung haben, sich jedoch in ihrer Sensitivität bzw. Spezifität unterscheiden. Tatsächlich sieht man eine Korrelation zwischen Temperatur und Ozonwerten bzw. Wind und Ozon, v. a. im Juli und August:

library(ggpubr)
library(gridExtra)
ozonplot <- ggplot(airquality, aes(x = Temp, y = Ozone, color = Monthname)) + # definiere, was geplottet wird
  geom_point() + # definiere Plottyp
    geom_hline(yintercept=70, linetype="dashed", color = "red") + #füge horizontale Linie hinzu, um Ozongrenzwert zu verdeutlichen
  geom_text(x=70, y= 100, label="Ozongrenzwert", color = "red") + # füge Beschriftung hinzu
  stat_cor(method = "pearson", label.x = 50, label.y = c(40, 55, 75, 90, 105), p.accuracy = 0.001, r.accuracy = 0.01) + # berechne Pearson-Korrelationskoeffizient
  geom_smooth(method='lm', formula= y~x, se=FALSE) + # füge Regressionsgerade hinzu
  theme_bw() +   # Aufhübschen des Plots
  ggtitle("Korrelaton zwischen Temperatur und Ozonwerten, nach Monat") # Überschrift

windplot <- ggplot(airquality, aes(x = Wind, y = Ozone, color = Monthname)) + # definiere, was geplottet wird
  geom_point() + # definiere Plottyp
  geom_hline(yintercept=70, linetype="dashed", color = "red") + #füge horizontale Linie hinzu, um Ozongrenzwert zu verdeutlichen
  geom_text(x=10, y= 100, label="Ozongrenzwert", color = "red") + # füge Beschriftung hinzu
  stat_cor(method = "pearson", label.x = 15, label.y = c(40, 55, 75, 90, 105), p.accuracy = 0.001, r.accuracy = 0.01) + # berechne Pearson-Korrelationskoeffizient
  geom_smooth(method='lm', formula= y~x, se=FALSE) + # füge Regressionsgerade hinzu
  theme_bw() +   # Aufhübschen des Plots
  gg title("Korrelation zwischen Wind und Ozonwerten, nach Monat")

grid.arrange(ozonplot, windplot, nrow=2)

Offenbar gibt es einen Zusammenhang zwischen Temperatur/Wind und Ozon, allerdings nur in den Monaten Juli und August. In den übrigen Monaten gab es ähnliche Wind- und Temperaturwerte, jedoch ohne dass der Ozonwert überschritten wurde. Das erklärt, warum einige Algorithmen "Monat" als hilfreichen Parameter ins Modell aufgenommen haben.

Die Kunst der Überzeugung

Nun, da wir uns ein Bild von den Daten und ihren Abhängigkeiten gemacht haben, gilt es, die Ergebnisse so präzise und überzeugend wie möglich darzustellen, sodass sich unsere Aussage durch eine geeignete Grafik von selbst aufdrängt. Dabei machen wir davon Gebrauch, dass wir es mit longitudinalen Daten zu tun haben und gewohnt sind, diese im Zeitverlauf zu betrachten. Zusätzlich werden wir uns noch die falsch vorhergesagten Werte genauer ansehen.

library(ggrepel)
airquality$id <- rownames(airquality)
Testergebnisse$prediction <- ifelse(Testergebnisse$class == Testergebnisse$pred, "correct", "wrong")
fehleranfaelligeProben <- Testergebnisse %>%
  filter(prediction == "wrong") %>% # suche alle falschen Vorhersagen heraus
  group_by(id) %>% # gruppiere nach Zeilenname als id des Datenpunkts
  tally() %>%  # zähle, wie oft die Vorhersage pro Tag falsch war
  filter(n>=5) %>% # filtere nach den wirklich schweren Fällen, die 5 mal falsch vorhergesagt wurden
  pull(id) # extrahiere die ID der problematischen Datensätze

airquality$fehleranfaellig <- ifelse(airquality$id %in% fehleranfaelligeProben, "fehlerhaft", "") # kennzeichne Tage mit fehlerhafter Vorhersage
# füge Datumspalte hinzu
airquality$date <- as.Date(paste0(airquality$Day, "-" ,airquality$Month), "%d-%m")

# Skaliere Messungen auf einheitliche Skala und generiere eine Tabelle im long Format
airquality_scaled <- airquality %>%
  mutate(Ozone = scale(Ozone), Temp = scale(Temp), Wind = scale(Wind)) %>%
  select(date, class, fehleranfaellig, Ozone, Temp, Wind) %>%
  pivot_longer(!c(date, class, fehleranfaellig), names_to = "Parameter", values_to = "Wert" )

# Bestimme Ozongrenzwert im skalierten Datensatz
yintercept <- min(airquality_scaled[airquality_scaled$Parameter == "Ozone" & airquality_scaled$class == "belastend", "Wert"])

ggplot(airquality_scaled, aes(x = date, y = Wert, color = Parameter, label = fehleranfaellig)) + # definiere, was geplottet wird
  geom_line() + # definiere Plottyp
  geom_hline(yintercept=yintercept, linetype="dashed", color = "red") + # füge horizontale Linie hinzu, um Ozongrenzwert zu verdeutlichen
  geom_vline(xintercept = unique(airquality_scaled$date[airquality_scaled$fehleranfaellig =="fehlerhaft"]), linetype="dashed", color = "black") + # füge vertikale Linie an den fehlerhaften Vorhersagen hinzu
  geom_text(x=as.Date("2023-06-01"), y= yintercept + 0.5, label="Ozongrenzwert überschritten", color = "red") + # füge Beschriftung hinzu
    geom_text_repel(data = subset(airquality_scaled, fehleranfaellig=="fehlerhaft" & Parameter=="Ozone"),
                    color = "black",
                         point.padding = 0.4,
                          nudge_x = 0.15,
                          nudge_y = 0.5,
                          segment.linetype = 2,
                          arrow = arrow(length = unit(0.025, "npc"))) +   # labeln der problematischen Vorhersagen
  theme_bw() + # Aufhübschen des Plots
  ggtitle("Ozon, Temperatur und Wind im Zeitverlauf", subtitle = "Werte skaliert") + # Titel
  ylab("Skalierter Wert") # Ändere Achsenbeschriftung

In der Grafik können wir sehr anschaulich erkennen, dass die Ozonpeaks vor allem bei Windflauten und hohen Temperaturen auftraten. Zwar wurde dieser Zusammenhang bereits durch die Korrelationsplots nahegelegt, hier wird es aber deutlich. Wir können nachvollziehen, warum die Modelle Temperatur und Wind als besonders wichtig für die Vorhersage erkannt haben. Das schafft Vertrauen, zeigt aber gleichzeitig noch die Herausforderungen, die es zu meistern gilt.

Wir sehen, dass der letzte problematische Tag tatsächlich sehr grenzwertig war und wir ihn wahrscheinlich deswegen falsch klassifiziert haben. An den anderen Tagen, mit wiederholt falsch vorausgesagten Ozonwerten, wurde der Ozongrenzwert trotz relativ starken Windes überschritten, was sonst nicht beobachtet wurde (wir erinnern uns an die negative Korrelation zwischen Wind und Ozon). Dies würde dazu einladen, sich doch noch die Solar.R-Werte anzusehen, die zumindest vom Random-Forest-Modell als wichtig erachtet wurden. Im Fall von Zeitreihen kann man auch statt herkömmlicher Klassifikationsalgorithmen solche benutzen, die die Abhängigkeit von zeitlich aufeinander folgenden Daten ausnutzen, wie z. B. ARIMA. Die Modellierung von Zeitreihen ist jedoch ein eigenes Thema, das wir hier ausklammern.

Schlussendlich

Dies war ein kleiner Ausflug in die Datenanalyse mit R. Wir haben den von R bereitgestellten Datensatz bereinigt und demonstriert, wie man relativ einfach einen Klassifikator mit maschinellen Lernansätzen erstellen kann. Die Analyse ist nicht abgeschlossen, aber wir haben einen Überblick gewonnen und einige Ideen, wie man die Vorhersage verbessern könnte. Gerade im Bereich der Modellierung muss man sich das Problem der Modelldrift vor Augen halten, d. h. ein erstelltes Modell wird mit der Zeit schlechtere Vorhersagen liefern, weil sich die Daten zunehmend von den fürs Training benutzten Daten wegbewegen. Deswegen ist man mit der Modellierung niemals fertig und tut gut daran, eine (halb-)automatisierte Pipeline aufzubauen, die regelmäßig neue Modelle generiert und das bisherige durch das am besten bewertete ersetzt. Auch dies lässt sich beispielsweise mit mlflow für R bewerkstelligen.

Ebenfalls steht der Analyst vor der Herausforderung, den Überblick sowohl über die zahlreichen Methoden der Statistik und des maschinellen Lernens, als auch über neue technische Entwicklungen und Pakete z. B. für R zu behalten. Wer seine Reise in die Welt der Datenanalyse mit R gerade beginnt, dem steht ein großes Angebot an Kursen und Büchern zur Verfügung, um die Grundlagen zu lernen oder sich in Teilaspekte zu vertiefen. Als Einstieg seien ein paar frei verfügbare Ressourcen genannt (s. u.).

Letztendlich muss die Analystin die Komfortzone der Anleitungen verlassen und beginnen, anhand von eigenen Projekten zu lernen. Erst dann begegnet man den vielen Herausforderungen und Fragen, die eine Erweiterung des Könnens ermöglichen. Dabei wünsche ich allen viel Erfolg!

Quellen
  1. R-bloggers
  2. Wikipedia: Luftqualitätsindex AQI
  3. EPA’s Ground-level Ozone Standard
  4. caret & Available Models

frei verfügbare Ressourcen:

Autorin

Elena Jolkver

Elena Jolkver arbeitet als Beraterin bei xValue und bearbeitet mit ihrem Team Forschungsfragen ihrer Kunden der Life Science Branche. Die entwickelten Modelle werden von ihnen mit RShiny-basierten Anw
>> Weiterlesen
Das könnte Sie auch interessieren
Kommentare (0)

Neuen Kommentar schreiben