PCA

Um zu demonstrieren, wie mit PCA Prognosen (predictions) gemacht werden können, benutzen wir wieder den Iris-Datensatz.

data(iris)
head(iris, 3)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa

Die letzte Spalte enthält die Bezeichnung der Art:

levels(iris[,5])
## [1] "setosa"     "versicolor" "virginica"

Die PCA kann mit Hilfe der Funktion prcomp durchgeführt werden. Die Daten sollten auf jeden Fall zentriert werden (center = TRUE), hier werden wir auch skalieren (scale = TRUE):

pca <- prcomp(iris[,-5], center = TRUE, scale = TRUE)

Die ersten beiden Hauptkomponenten können mit Hilfe der autoplot-Funktion aus dem Package ggfortify abgebildet werden. Das Zuordnen der Farben zu den Arten ist hier sehr elegant möglich (mit colour = "Species"):

suppressWarnings(suppressPackageStartupMessages(library(ggfortify)))
autoplot(pca, data = iris, colour = "Species")

Wir sehen, dass die Arten aufgrund ihrer Position im Plot hinreichend gut separiert werden können, auf jeden Fall trifft das auf setosa zu. Folglich lässt sich die PCA zur Prognose (“prediction”) verwenden, d.h. wir können bisher noch nicht klassifizierte Exemplare aufgrund dieser vier morphologischen Daten einer der drei Arten zuordnen.

Prognose

Eine Prognose erstellen heißt in dem hier betrachteten Zusammenhang, dass bisher nicht klassifizierte Exemplare der Gattung Iris aufgrund der vorliegenden morphologischen Daten einer der 3 Arten zugeordnet werden.
Wir wenden dabei das sogenannte supervised learning an. Beim supervised learning benutzen wir Daten mit bekannter Artzugehörigkeit als Trainingsdaten. Mit diesen Trainingsdaten erstellen wir ein Modell.
Bei einem anderen Teil der Daten - den Testdaten - ist die Artzugehörigkeit nicht bekannt, sondern es liegen nur die morphologischen Daten vor. Mit Hilfe des Modells soll für die Testdaten nun die Artzugehörigkeit bestimmt werden.

Wir wollen im Folgenden 90% der Iris-Daten verwenden, um mittels PCA ein Modell zu erstellen. Das Modell enthält Informationen darüber, wie die Lage der Punkte im Plot der Hauptkomponenten mit einer der 3 Arten verknüpft ist.

Dann werden wir kurz vergessen, dass wir für die restlichen 10% der Exemplare eigentlich die Arten-Bezeichnungen schon kennen. Diese Arten-Bezeichnungen versuchen wir dann herauszufinden, indem wir ihre Lage in dem mit dem Modell erstellten Hauptkomponenten-Plot auswerten.

Zur programmtechnischen Durchführung wählen wir mit dem sample-Kommando zufällig 15 Zahlen aus den Zahlen 1 bis 150 aus - dies sind dann die Indizes der Testdaten:

set.seed(333)
ind = sample(1:150, 15)
ind
##  [1] 142  55  66  39  60   2 134  67  44 103  40  63  87  75  90
testdata = iris[ind,-5]
dim(testdata)
## [1] 15  4
head(testdata, 3)
##     Sepal.Length Sepal.Width Petal.Length Petal.Width
## 142          6.9         3.1          5.1         2.3
## 55           6.5         2.8          4.6         1.5
## 66           6.7         3.1          4.4         1.4

Die anderen 135 Exemplare dienen als Trainingsdaten:

traindata = iris[-ind,-5]
dim(traindata)
## [1] 135   4
head(traindata, 3)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1          5.1         3.5          1.4         0.2
## 3          4.7         3.2          1.3         0.2
## 4          4.6         3.1          1.5         0.2

Um das Modell zu erstellen, führen wir nun die PCA mit den Trainingsdaten durch. Das Ergebnis der PCA wird zunächst im Objekt pca_model gespeichert:

pca_model = prcomp(traindata, center = TRUE, scale = TRUE)  # new object, contains the model

Alle notwendigen Parameter der Projektion sind im Objekt pca_model gespeichert.

Wir plotten jetzt die ersten beiden Hauptkomponenten der Trainingsdaten. Die Farben werden anhand der Species-Bezeichnung gewählt:

color_num = as.numeric(iris[-ind,]$Species) + 1 # specify colours for training data
plot(pca_model$x[,1:2], col = color_num, xlab = "PC1", ylab = "PC2", main = "PCA training results", font.main = 1, pch = 16)
legend(-1, 2.7, c("setosa", "versicolor", "virginica"), col = c(2,3,4), pch = 16)

Nun der eigentliche Schritt bei der Prognose’: Wir wenden auf die Testdaten genau die gleiche Projektion an, welche bereits für die Trainingsdaten verwendet wurde. Genauer gesagt: es wird die gleiche Drehmatrix auf die Testdaten angewendet wie zuvor auf die Trainingsdaten. Die Drehmatrix besteht aus den Eigenvektoren der Kovarianzmatrix der Trainingsdaten und ist im Objekt pca_model enthalten.
Technisch geschieht das mittels der Funktion predict:

prediction = predict(pca_model, newdata = testdata)
head(prediction, 3)
##          PC1        PC2        PC3       PC4
## 142 1.881065  0.6642441  0.1110383 0.4697214
## 55  1.070844 -0.2125961 -0.4114502 0.1124056
## 66  0.874583  0.4887231 -0.5104178 0.1157040

Das Objekt prediction enthält die Positionen der Testdaten im Plot der Hauptkomponenten (welcher mit den Trainingsdaten erzeugt wurde). Wir werden jetzt die Testdaten zum Plot hinzugefügen. Zunächst wird noch einmal der Plot mit den Trainingsdaten erstellt. Von den projizierten Testdaten werden ebenfalls nur die ersten beiden Haupkomponenten geplottet. Wir können den Testdaten die bekannten korrekten Farben zuordnen, um zu sehen, wie gut die Prognose funktioniert (im Realfall kann man dies bei der Prognose natürlich nicht, da die Arten-Bezeichnungen ja eben nicht bekannt sind, sondern prognostiziert werden sollen). Wir plotten diesmal die Trainingsdaten als kleine Kreise (pch = 1) und die Testdaten als große ausgefüllte Kreise (pch = 16, cex = 1.5), um die Punkte unterscheiden zu können:

color_num = as.numeric(iris[-ind,]$Species) + 1 # specify colours for training data
plot(pca_model$x[,1:2], col = color_num, xlab = "PC1", ylab = "PC2", main = "PCA results", font.main = 1, pch = 1)
legend(-1, 2.7, c("setosa", "versicolor", "virginica"), col = c(2,3,4), pch = 1)
test_colors = as.numeric(iris[ind,]$Species) + 1                   # colours for the test data
points(prediction[,1:2], pch = 16, cex = 1.5, col = test_colors)   # test data

Dies scheint hinreichend gut geklappt zu haben. Die Testdaten könnnen fast alle anhand ihrer Lage im Plot der korrekten Art zugeordnet werden. Das heißt: hätten wir die wahren Arten-Bezeichnungen der Testdaten nicht gekannt, hätten wir diese nun anhand der Lage der Punkte in der Abbildung erraten können. Probleme würde es eventuell bei der Unterscheidung von versicolor und virginica geben, da diese schon in den Trainingsdaten recht nahe beieinander liegen - kein Modell ist perfekt!.


Note that all code snippets are designed for comprehensibility, not for efficiency!


uwe.menzel@matstat.org