Datensätze mit mehr Variablen als Fällen sind eine besondere Herausforderung für die Datenanalyse: p > n, p für predictors, Prädiktoren; n für die Stichprobengröße. Klassische Verfahren wie die lineare Regression sind unter diesen Bedingungen rechnerisch nicht lösbar.
Die Daten: p > n
Hier ein Beispiel: biomarker.Rda enthält 90 Fälle (Beobachtungen) von 2000 unabhängigen Variablen – hier einfach von X1 bis X2000 numeriert. Die abhängige Variable heißt result: Typ 1 kommt 24 mal vor, Typ 2 66 mal.
Die Werte der X-Variablen liegen im Bereich von 0 bis 100, sie sind in etwa normalverteilt.
Modellierung: Vor-Auswahl der Prädiktoren
2000 Prädiktoren sind ganz schön viel. Nehmen wir an, es handelt sich um aufwändige Messungen, z. B. Eigenschaften von Molekülen. Wäre es nicht sinnvoll, die Auswahl etwas einzuschränken? Viele Algorithmen des Maschinellen Lernens liefern Aussagen zur Wichtigkeit von Prädiktoren (variable importance). Wir verwenden einen einfachen Entscheidungsbaum (rpart: Recursive Partitioning and Regression Trees), der uns die 18 wichtigsten Prädiktoren heraussucht. Hier der R-Code:
Modellierung: Random Forest mit Kreuzvalidierung
Nun geht’s an die Modellbildung. Um eine realistische Einschätzung der Leistungsfähigkeit unseres Modells zu erhalten, setzen wir Kreuzvalidierung (cross validation) ein. Die Idee dahinter ist simpel: Wir beurteilen das Modell anhand von neuen Daten, d. h. Daten, die nicht zur Modellbildung verwendet wurden. Da wir keinen neuen Datensatz zur Verfügung haben, ziehen wir mehrfach Stichproben aus unseren vorhandenen Daten. Wir setzen 5-fache Kreuzvalidierung ein. Dabei werden fünf Teildatensätze gebildet, von denen vier zur Modellbildung dienen und der fünfte als Testdatensatz, um die Modellgüte zu beurteilen. Dieses Vorgehen wird fünf mal wiederholt, bis jeder der fünf Teildatensätze ein Mal Testdatensatz war.
Das R-Paket caret erleichtert uns diese Arbeit. Die train-Funktion baut eine einheitliche Brücke zu einer Vielzahl von Machine Learning-Algorithmen, die ansonsten mit uneinheitlich benannten Parametern angesprochen werden müssten. Zunächst bilden wir den Teildatenatz mit den 18 ausgewählten Prädiktoren.
biomarker2 <- biomarker[, c(names(tree$variable.importance), "result")] library(caret) ctrl <- trainControl(method = "cv", number = 5, classProbs = TRUE, summaryFunction = twoClassSummary) set.seed(2018) biomarker_rf2 <- train(result ~ ., data = biomarker2, trControl = ctrl, method = "rf", metric = "ROC")
Hier die Ergebnisse:
Das beste Modell verwendet alle 18 Prädiktoren. (Wir verzichten darauf, den mtry-Parameter feiner einzustellen – caret hat als Voreinstellung drei Modelle mit 2, 10 und 18 Variablen pro Verzweigung verglichen.) Die Fläche unter der ROC-Kurve beträgt 0,85. Ein „perfektes“ Modell liegt bei 1, ein Zufalls-Klassifizierer bei 0,5. Wir schneiden also deutlich besser ab, als wenn wir den Typ (Typ 1 oder Typ 2) raten müssten. Sensitivität und Spezifität sind nicht perfekt, aber einigermaßen brauchbar.
Modell mit allen Prädiktoren
Zum Vergleich berechnen wir auch ein Modell, das über alle 2000 Prädiktoren verfügen darf, nicht nur über die vorausgewählten 18. Die Syntax ist identisch zu oben, nur dass wir auf die Originaldaten zugreifen:
Das beste Modell verwendet tatsächlich alle 2000 Prädiktoren. Doch was ist mit der Modellgüte geschehen? Die Fläche unter der ROC-Kurve ist deutlich schwächer – mit 0,57 nur noch knapp besser als der Zufalls-Klassifizierer. Ein Modell mit 63 Prädiktoren erreicht gar nur einen Wert von 0,51. Die Sensitivität – eine glatte 0! Völlig unbrauchbar – das Modell kann gar nicht zwischen Typ 1 und Typ 2 differenzieren. Sollte das Modell mit allen Prädiktoren nicht besser sein als eines, das nur 18 vorausgewählte Prädiktoren verwenden darf? Wie konnte das passieren?
Die Daten: Reine Zufallswerte!
Zur Erklärung betrachten wir, wie die Daten entstanden sind. Ich habe sie nicht irgendwo heruntergeladen oder selbst gemessen, sondern:
set.seed(2018)
biomarker <- data.frame(replicate(2000, sample(0:100, 90, replace = TRUE)))
biomarker$result <- sample(c("Typ1", "Typ2"), size = 90, replace = TRUE, prob = c(0.3, 0.7))
biomarker$result <- factor(biomarker$result)
save(biomarker, file = "biomarker.Rda")
Die Daten sind komplett zufällig! Sowohl die 2000 Prädiktoren als auch die abhängige Variable result resultieren aus Zufalls-Verteilungen. Um die Ergebnisse reproduzierbar zu machen, habe ich per set.seed()-Funktion den Startwert des Zufallsgenerators fixiert.
Was bedeutet das für die bisherigen Ergebnisse?
Kreuzvalidierung: Gesamten Modellbildungsprozess einbeziehen!
Da es sich um Zufallsdaten handelt, ist offenbar das zweite Modell mit allen 2000 Prädiktoren das ehrlichere: Es liefert zwar eine schlechtere Modellgüte, aber das ist realistisch, denn in Wirklichkeit sollte man mit reinen Zufallsdaten kein Vorhersagemodell bauen können. Die Frage lautet also nicht mehr, warum das zweite Modell so schlecht ist, sondern vielmehr: Warum liefert uns das erste Modell trotz Kreuzvalidierung eine so unrealistische Schätzung der Modellgüte?
Die Antwort lautet: Der gesamte Modellbildungsprozess muss in die Kreuzvalidierung einbezogen werden. Das ist beim ersten Modell nicht geschehen: Wir haben 18 Prädiktoren vor-ausgewählt, doch dieser Prozess ist nicht generalisierbar. Bei so vielen Prädiktoren und einer so geringen Fallzahl lassen sich gewisse zufällige Korrelationen nicht ausschließen!
Wie man hört, wird dieser Fehler bei Daten mit p > n häufig gemacht! Darauf weisen zum Beispiel Trevor Hastie und Rob Tibshirani sowie Max Kuhn und Kjell Johnson in ihren empfehlenswerten Büchern hin.
An Introduction to Statistical Learning: with Applications in R – Hastie, Tibshirani, Witten, James
Applied Predictive Modeling: Max Kuhn, Kjell Johnson
Neuer: Tidy Modeling with R: A Framework for Modeling in the Tidyverse – Kuhn, Silge
Hier noch das Video auf Englisch: