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.

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

Ozon ist ein Bestandteil des Luftqualitätsindex AQI[2]. Dieser setzt sich aus Ozon, Partikelverschmutzung, CO, SO₂ und NO₂ 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

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