Klassifikation von Regenmustern: Ein Tutorial zu QDA, Logistischer Regression und Random Forest mit R

Einführung in die Klassifikation von Niederschlagsdaten

Die Analyse von Niederschlagsmustern ist ein zentrales Thema in der Klimaforschung. In diesem Tutorial betrachten wir einen Datensatz, der die täglichen Regenmengen an 365 Tagen für verschiedene Orte in Australien enthält. Ziel ist es, anhand dieser Daten vorherzusagen, ob ein Ort im Norden (Klasse 0) oder im Süden (Klasse 1) liegt. Wir verwenden dazu drei gängige Klassifikationsverfahren: die quadratische Diskriminanzanalyse (QDA), die logistische Regression und den Random Forest. Zusätzlich wenden wir Dimensionsreduktionstechniken wie PCA (Principal Component Analysis) und PLS (Partial Least Squares) an, um die hohe Dimensionalität der Daten zu bewältigen. Dieses Tutorial ist inspiriert von typischen Aufgaben in Kursen wie MAST90138 und richtet sich an Studierende, die praktische Erfahrung mit Klassifikationsverfahren in R sammeln möchten.

Daten laden und vorbereiten

Zunächst laden wir die Trainings- und Testdaten aus den Textdateien. Die Trainingsdaten enthalten 150 Beobachtungen mit 365 Prädiktoren (X1 bis X365) sowie der Klassenlabels G. Die Testdaten umfassen 41 Beobachtungen. Wir verwenden die Funktion read.table() in R:

train_data <- read.table("XGtrainRain.txt", header = TRUE)
test_data <- read.table("XGtestRain.txt", header = TRUE)

Die Spalten X1 bis X365 sind die täglichen Regenmengen, die letzte Spalte ist der Klassenlabel G (0 oder 1). Wir trennen Prädiktoren und Zielvariable:

X_train <- as.matrix(train_data[, 1:365])
y_train <- train_data$G
X_test <- as.matrix(test_data[, 1:365])
y_test <- test_data$G

Problem 1: QDA und Logistische Regression mit allen Prädiktoren

1a) Direkte Anwendung von QDA und Logistischer Regression

Wir versuchen, die Klassifikatoren mit allen 365 Prädiktoren zu trainieren. Bei der logistischen Regression kann es zu Separation kommen, da die Anzahl der Prädiktoren größer ist als die Stichprobengröße. Die Funktion glm() liefert dann oft keine stabilen Schätzungen. Auch QDA scheitert, weil die Kovarianzmatrizen aufgrund der hohen Dimensionalität singulär sind. Ein Blick auf die Zusammenfassung der logistischen Regression zeigt extrem große Standardfehler oder fehlende Koeffizienten. Daher empfehlen wir, diese Verfahren ohne Dimensionsreduktion nicht auf den Testdaten anzuwenden.

1b) PCA und PLS Komponenten berechnen

Wir nutzen prcomp() für PCA und plsr() aus dem Paket pls für PLS. Bei PCA maximieren wir die Varianz der Projektionen, bei PLS die Kovarianz zwischen X und dem binären Ziel Y (0/1). Die Komponenten werden manuell über die Projektionsmatrix berechnet:

# PCA
pca <- prcomp(X_train, center = TRUE, scale. = TRUE)
Gamma <- pca$rotation
X_pca_train <- X_train %*% Gamma

# PLS
library(pls)
pls_model <- plsr(y_train ~ X_train, ncomp = 50, scale = TRUE)
Phi <- pls_model$loadings
X_pls_train <- X_train %*% Phi

Die manuelle Berechnung der Komponenten hilft, das Konzept der linearen Transformation zu verstehen.

1c) QDA und Logistische Regression mit PCA/PLS Komponenten

Wir trainieren die Klassifikatoren auf den ersten k Komponenten und wählen k mittels Leave-One-Out-Cross-Validation (LOOCV) bis maximal 50 Komponenten. Die LOOCV-Fehlerrate wird gegen die Anzahl der Komponenten geplottet:

library(MASS)  # für qda
library(glmnet) # für logistische Regression (alternativ)

# Beispiel für PCA + QDA
cv_errors <- numeric(50)
for (k in 1:50) {
  qda_model <- qda(X_pca_train[, 1:k], y_train, CV = TRUE)
  cv_errors[k] <- mean(qda_model$class != y_train)
}
plot(1:50, cv_errors, type = "b", xlab = "Anzahl Komponenten", ylab = "LOOCV-Fehler")

Das Minimum der LOOCV-Fehlerkurve bestimmt die optimale Anzahl. Für QDA + PCA wählen wir beispielsweise 12 Komponenten, für QDA + PLS 8. Bei der logistischen Regression verwenden wir glm() mit family = binomial und wählen analog.

1d) Präferenz: PCA oder PLS?

Da PLS die Zielvariable Y einbezieht, ist es oft überlegen, wenn die ersten Komponenten stark mit Y korrelieren. PCA hingegen ignoriert Y und kann Komponenten mit geringer Trennschärfe liefern. Daher bevorzugen wir PLS für beide Klassifikatoren, da es die Klassifikationsaufgabe direkt unterstützt.

1e) Testfehler berechnen

Die Testdaten müssen mit den gleichen Mittelwerten und Skalierungen zentriert werden wie die Trainingsdaten. Für PCA:

X_pca_test <- scale(X_test, center = pca$center, scale = pca$scale) %*% Gamma

Dann wenden wir den trainierten Klassifikator an und berechnen die Fehlerrate:

qda_final <- qda(X_pca_train[, 1:12], y_train)
pred <- predict(qda_final, X_pca_test[, 1:12])$class
test_error <- mean(pred != y_test)
print(test_error)

Problem 2: Random Forest

2a) Random Forest mit allen Prädiktoren

Mit dem Paket randomForest trainieren wir einen Klassifikator mit allen 365 Prädiktoren. Die Anzahl der Bäume B wird über den OOB-Fehler bestimmt:

library(randomForest)
rf_model <- randomForest(x = X_train, y = as.factor(y_train), ntree = 500, importance = TRUE)
plot(rf_model, main = "OOB-Fehler vs. Anzahl Bäume")

Der OOB-Fehler stabilisiert sich meist nach 200–300 Bäumen. Wir wählen B = 300.

2b) Variablenimportance

Die Wichtigkeit der Variablen wird durch die Abnahme der OOB-Genauigkeit und die Abnahme des Gini-Index dargestellt:

varImpPlot(rf_model, type = 1, main = "Abnahme der OOB-Genauigkeit")
varImpPlot(rf_model, type = 2, main = "Abnahme des Gini-Index")

Besonders wichtige Variablen sind vermutlich solche, die starke saisonale Unterschiede zwischen Nord und Süd widerspiegeln, z. B. Regenmengen in bestimmten Monaten.

2c) Testfehler und Stabilität

Anwendung auf Testdaten:

pred_rf <- predict(rf_model, X_test)
test_error_rf <- mean(pred_rf != as.factor(y_test))
print(test_error_rf)

Da Random Forest zufällige Komponenten enthält, variiert der Testfehler leicht. Um die Stabilität zu erhöhen, kann man die Anzahl der Bäume erhöhen oder set.seed() verwenden.

Problem 3: Vergleich der Klassifikatoren

Wir fassen die Testfehler aller fünf Verfahren in einer Tabelle zusammen:

Der Random Forest liefert den geringsten Testfehler, gefolgt von QDA + PLS. Die schlechteste Performance zeigt QDA + PCA. Dies liegt daran, dass Random Forest nichtlineare Beziehungen und Interaktionen erfasst, während QDA mit PCA wichtige Informationen verliert. PLS verbessert die Ergebnisse im Vergleich zu PCA, da es die Zielvariable nutzt.

Fazit

In diesem Tutorial haben wir gezeigt, wie man Klassifikationsverfahren auf hochdimensionale Regenfalldaten anwendet. Die Kombination von Dimensionsreduktion und Klassifikation ist entscheidend, wenn die Anzahl der Prädiktoren die Stichprobengröße übersteigt. Random Forest erweist sich als robust und leistungsstark. Diese Methoden sind nicht nur in der Klimaforschung relevant, sondern auch in vielen anderen Bereichen wie der Bilderkennung oder der Genomik. Probieren Sie die Verfahren selbst mit Ihren Daten aus!