1. Vorbetrachtung

Wir werden hier einen der in R enthaltenen Datensätze herunterladen und dann die PCA mit Hilfe der Funktion prcomp durchführen. Einige theoretische Grundlagen sind ja schon von den vorhergehenden Lektionen bekannt.


2. Der Iris-Datensatz

Wir laden zunächst den Iris-Datensatz herunter (Iris = Schwertlilie). Dieser enthält morphologische Daten für die drei Arten Iris setosa, Iris virginica und Iris versicolor der Gattung Iris. Von jeder Art sind bei je 50 Exemplaren jeweils die Breite und Länge des Kelchblatts (Sepal.Width, Sepal.Lenght) sowie die Breite und Länge des Blütenblatts (Petal.Width, Petal.Lenght) vermessen worden:



Anhand der morphologischen Daten können die Arten klassifiziert werden.



Werfen wir zunächst eine Blick auf den Datensatz:

data(iris)
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
dim(iris)
## [1] 150   5
head(iris)
##   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
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Die letzte Spalte enthält die Bezeichnung der Art:

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

Da die letzte Spalte vom Datentyp factor ist, können wir das levels-Kommando anwenden, um zu sehen, welche Bezeichungen in dieser Spalte vorkommen. Wenn die letze Spalte vom Typ character wäre, würden wir dafür das unique-Kommando verwenden.


3. Daten plotten

Zunächst entfernen wir die letzte Spalte aus dem Datensatz, da diese keine numerischen Daten enthält:

data4 = iris[,-5]
head(data4)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1          5.1         3.5          1.4         0.2
## 2          4.9         3.0          1.4         0.2
## 3          4.7         3.2          1.3         0.2
## 4          4.6         3.1          1.5         0.2
## 5          5.0         3.6          1.4         0.2
## 6          5.4         3.9          1.7         0.4

Einen 4-dimensionalen Datensatz können wir nicht mehr bildlich darstellen. Wir können aber z.B. die Funktion scatterplotMatrix aus dem Package car nutzen, um jede Variable gegen jede andere zu plotten:

suppressPackageStartupMessages(library("car"))
scatterplotMatrix(data4)

In der Diagonale werden die Wahrscheinlichkeitsdichtefunktion der jeweiligen Variablen angezeigt.

Es sind deutliche Korrelationen zwischen verschiedenen Variablenpaaren sichtbar, denn die Punktwolken haben eine gestreckte Form, d.h. man kann aus der Kenntnis einer Variable “ungefähr” die Größe der anderen Variable erraten. (Dies ist ja auch kein Wunder, denn eine große Blüte wird sowohl eine große Länge als auch eine große Breite haben)


4. Noch einmal etwas Theorie

Bei der Hauptkomponentenanalyse wird das Koordinatensystem so gedreht, dass die 1. Hauptkomponente in Richtung der größten Varianz der mehrdimensionalen Punktwolke zeigt. Die folgenden Hauptkomponenten stehen jeweils senkrecht zu den vorhergehenden und enthalten jeweils weniger Varianz. Die Drehung des Koordinatensystems erfolgt mathematisch mit Hilfe einer Matrix, deren Spalten aus den Eigenvektoren der Kovarianzmatrix gebildet werden. In den folgenden Formeln werden, um Schreibarbeit zu sparen, nur 3-dimensionale Daten betrachtet. Wir indizieren die Eigenvektoren mit den Buchstaben, \(A\), \(B\) und \(C\). Der erste Eigenvektor lautet also z.B.

\[ \boldsymbol{E}_1 = \begin{pmatrix} \varepsilon_{A,1} \\ \varepsilon_{A,2} \\ \varepsilon_{A,3} \end{pmatrix} \]

Die Drehmatrix ist damit

\[ \boldsymbol{E} = \begin{pmatrix} \varepsilon_{A,1} & \varepsilon_{B,1} & \varepsilon_{C,1} \\ \varepsilon_{A,2} & \varepsilon_{B,2} & \varepsilon_{C,2} \\ \varepsilon_{A,3} & \varepsilon_{B,3} & \varepsilon_{C,3} \end{pmatrix} \] In Matrixschreibweise können wir die Drehung so schreiben:

\[ \boldsymbol{X}^{\prime} = \boldsymbol{E}^T \cdot \boldsymbol{X} \] (\(\boldsymbol{E}^T\) ist die zu \(\boldsymbol{E}\) transponierte Matrix).

Dasselbe in Koordinatenschreibweise:

\[ \begin{pmatrix} X_1^{\prime} \\ X_2^{\prime} \\ X_3^{\prime} \end{pmatrix} = \begin{pmatrix} \varepsilon_{A,1} & \varepsilon_{A,2} & \varepsilon_{A,3} \\ \varepsilon_{B,1} & \varepsilon_{B,2} & \varepsilon_{B,3} \\ \varepsilon_{C,1} & \varepsilon_{C,2} & \varepsilon_{C,3} \end{pmatrix} \cdot \begin{pmatrix} X_1 \\ X_2 \\ X_3 \end{pmatrix} \hspace{2cm} (1) \]

So ergibt sich z.B. für \(X_1^{\prime}\):

\[ X_1^{\prime} = \varepsilon_{A,1} \cdot X_1 + \varepsilon_{A,2} \cdot X_2 + \varepsilon_{A,3} \cdot X_3 \]

Im neuen Koordiantensystem ist die Kovarianzmatrix dann eine Diagonalmatrix, in deren Hauptdiagonale die Eigenwerte der Kovarianzmatrix stehen:

\[ Cov(X^{\prime}) = \begin{pmatrix} \varepsilon_{1} & 0 & 0 \\ 0 & \varepsilon_{2} & 0 \\ 0 & 0 & \varepsilon_{3} \end{pmatrix} \]

Die Eigenwerte geben also die Varianzen der neuen Variablen an: \(Var(X_1^{\prime}) = \varepsilon_1\), \(Var(X_2^{\prime}) = \varepsilon_2\) und \(Var(X_3^{\prime}) = \varepsilon_3\). Desto größer der erste Eigenwert im Vergleich zu den anderen ist, desto größeren Anteil hat die erste Hauptkomponente an der Gesamtvarianz. Hauptkomponenten mit sehr kleinen zugehörigen Eigenwerten können vernachlässigt werden, ohne viel Information über den Datensatz zu verlieren. Um einzuschätzen, wie viele Hauptkomponenten in einen Datensatz mit reduzierter Dimension einbezogen werden müssen, wird oft ein Balkendiagramm mit den zugehörigen Eigenwerten (= Varianzen) geplottet (siehe unten).

Anhand der Faktoren \(\varepsilon_{i,j}\) in Gleichung (1) kann man ersehen, welche der ursprünglichen Variablen am meisten zu den einzelnen Hauptkomponenten beitragen, d.h. welche ursprünglichen Variablen am meisten zur Varianz des Datensatzes beitragen. Man kann also mit der PCA Aufschlüsse über wichtige Einflussfaktoren auf einen Datensatz gewinnen.

Mehr Theorie im Wikipedia-Eintrag zur Hauptkomponentenanalyse .


5. PCA mit prcomp

Das obligatorische Zentrieren und Skalieren der Daten wird von prcomp übernommen, wenn center = TRUE und scale = TRUE gesetzt wird:

pca <- prcomp(data4, center = TRUE, scale = TRUE)

Alle Ergebnisse sind nun in dem Objekt pca gespeichert. Wir werden weiter unten dieses Objekt mehrmals benutzen. Die Struktur eines Objektes kann mit dem str-Kommando ermittelt werden.

str(pca)
## List of 5
##  $ sdev    : num [1:4] 1.708 0.956 0.383 0.144
##  $ rotation: num [1:4, 1:4] 0.521 -0.269 0.58 0.565 -0.377 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
##   .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
##  $ center  : Named num [1:4] 5.84 3.06 3.76 1.2
##   ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
##  $ scale   : Named num [1:4] 0.828 0.436 1.765 0.762
##   ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
##  $ x       : num [1:150, 1:4] -2.26 -2.07 -2.36 -2.29 -2.38 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
##  - attr(*, "class")= chr "prcomp"

Das Subset sdev enthält die Standardabweichungen entlang der Hauptkomponenten:

pca$sdev
## [1] 1.7083611 0.9560494 0.3830886 0.1439265

… während in rotation die berechneten Eigenvektoren gespeichert sind:

pca$rotation
##                     PC1         PC2        PC3        PC4
## Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
## Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
## Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
## Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

Wir überzeugen uns noch einmal, dass die gefundenen Hauptkomponenten PC1 bis PC4 mit den Eigenwerten der Kovarianzmatrix identisch sind. Die Kovarianzmatrix wird aus den empirischen Daten mit Hilfe der Funktion cov berechnet. Wir werden vorher wieder die Daten mit zscore skalieren:

zscore <- function(vector) {return((vector - mean(vector))/sd(vector))}
dfz = apply(data4, 2, zscore)
evectors = eigen(cov(dfz))$vectors    # eigenvectors
colnames(evectors) = c("PC1", "PC2", "PC3", "PC4")
evectors
##             PC1         PC2        PC3        PC4
## [1,]  0.5210659 -0.37741762  0.7195664  0.2612863
## [2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
## [3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
## [4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

Die gefundenen Eigenvektoren stimmen genau mit dem Subset rotation überein.

Wir hatten schon in einer vorherigen Lektion sowie oben in Abschnitt 4 gesehen, dass die Eigenwerte der Kovarianzmatrix gleich den Varianzen entlang der Hauptkomponenten sind. Die Standardabweichungen entlang der Achsen entsprechen den Quadratwurzeln der Varianzen. Hier der Vergleich:

evalues = eigen(cov(dfz))$values    # eigenvalues
evalues                             # eigenvalues = variances
## [1] 2.91849782 0.91403047 0.14675688 0.02071484
sqrt(evalues)
## [1] 1.7083611 0.9560494 0.3830886 0.1439265
pca$sdev                            # standard deviations = square root of variances
## [1] 1.7083611 0.9560494 0.3830886 0.1439265

Also nochmal zum Merken:

\[ Eigenvalues = Variances = Standard\_Deviations^2 \]

Die Varianzen entlang aller Hauptkomponenten kann man mühelos plotten, indem man die Funktion plot auf das oben generierte Objekt pca anwendet:

plot(pca, col = "red")

Die relativen Anteile der einzelnen Hauptkomponenten an der Varianz kann man durch Normierung der Eigenwerte auf 1 oder auf 100% berechnen:

perc = evalues/sum(evalues)*100     # proportions (%)
round(perc, 2)
## [1] 72.96 22.85  3.67  0.52

Wir sehen, dass der weitaus größte Anteil der Gesamtvarianz auf die erste und zweite Hauptkomponente entfällt. Wir können also einen 2-dimensionalen Plot mit den ersten beiden Hauptkomponenten erstellen, ohne viel Information über den Datensatz zu verlieren.

Die transformierten Daten befinden sich im Subset x des Objektes pca:

str(pca$x)
##  num [1:150, 1:4] -2.26 -2.07 -2.36 -2.29 -2.38 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"

Wir können diesen Subset also benutzen, um die ersten 2 Hauptkomponenten zu plotten. Dabei macht es sicherlich Sinn, die einzelnen Arten der Gattung unterschiedlich einzufärben, um zu sehen, ob diese durch PCA klassifiziert werden können.

color_num = as.numeric(iris$Species) + 1  # specify colours (add 1 for nicer colouers)
color_num
##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
##  [75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4
## [112] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [149] 4 4

Das Konstrukt as.numeric(iris$Species) wandelt hier zunächst die factor-Variable in der letzten Spalte der Iris-Daten in ganze Zahlen um. (Außerdem addieren wir zu allen Zahlen eine 1, um die Farben schöner zu machen …)

Nun der Plot:

xtext = paste("PC1 (", round(perc[1], 2), "% )")
ytext = paste("PC2 (", round(perc[2], 2), "% )")
plot(pca$x[,1:2], col=color_num, xlab=xtext, ylab=ytext, main="PCA results", font.main=1, pch=16)
legend("topright", c("setosa", "versicolor", "virginica"), col = c(2,3,4), pch = 16)

Bei der Legende muss man aufpassen, dass die Farben korrekt den Arten zugeordnet werden.

Die ersten beiden Hauptkomponeneten können auch mit Hilfe der autoplot-Funktion aus dem Package ggfortify abgebildet werden. Das Zuordnen der Farben zu den Arten ist hier sehr elegant möglich:

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

Wir sehen, dass die Arten aufgrund ihrer Position im Plot der Hauptkomponenten hinreichend gut separiert werden können, zumindestens 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. Mehr dazu im nächsten Abschnitt.

Wenn es tatsächlich nötig sein sollte, die ersten drei Hauptkomponenten zu plotten, dann kann man z.B. auf plot3d aus dem Package rgl zurückgreifen:

xtext = paste("PC1 (", round(perc[1], 2), "% )")
ytext = paste("PC2 (", round(perc[2], 2), "% )")
ztext = paste("PC3 (", round(perc[3], 2), "% )")
par3d(windowRect = c(819, 88, 1492, 706))
plot3d(pca$x[,1:3], col=color_num, xlab=xtext, ylab=ytext, zlab=ztext, size = 5, main="PCA results", font.main=1, pch=16)
rglwidget()


Die Darstellung kann wieder mit der Maus vergrößert/verkleinert oder gedreht werden.


6. Prognosen mit PCA

Eine Prognose erstellen heißt in dem hier betrachteten Zusammenhang, dass bisher nicht klassifizierte Exemplare der Gattung Iris aufgrund der vorliegenden morphologischen Daten in eine der 3 Arten eingeteilt werden können. Wir wollen im Folgenden 90% der Iris-Daten verwenden, um mittels PCA ein Modell zu erstellen. Mit Modell ist bei der PCA im Prinzip die zur Projektion verwendete Drehmatrix gemeint (welche aus den Eigenvektoren der Kovarianzmatrix besteht). 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 betrachten. Es handelt sich dabei also um “supervised learning”.

Wir nehmen 10% der Exemplare als Testdaten aus den Datensatz heraus und führen die PCA zunächst mit den anderen 90% durch. Letztere sind dann die Trainingsdaten, welche zur Erstellung des Modells dienen. Mit den sample-Kommando werden zufällig 15 Zahlen aus den Zahlen 1 bis 150 ausgewählt - 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 = data4[ind,]
dim(testdata)
## [1] 15  4
head(testdata)
##     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
## 39           4.4         3.0          1.3         0.2
## 60           5.2         2.7          3.9         1.4
## 2            4.9         3.0          1.4         0.2

Die anderen 135 Exemplare dienen als Trainingsdaten:

traindata = data4[-ind,]
dim(traindata)
## [1] 135   4
head(traindata)
##   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
## 5          5.0         3.6          1.4         0.2
## 6          5.4         3.9          1.7         0.4
## 7          4.6         3.4          1.4         0.3

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_train gespeichert:

pca_train = prcomp(traindata, center = TRUE, scale = TRUE)  # new object, just for the training data

Wir plotten jetzt erste und zweite Hauptkomponente. Die Farben müssen jetzt auf die Trainingsdaten bezogen werden, also auf iris[-ind,]:

color_num = as.numeric(iris[-ind,]$Species) + 1 # specify colours for training data
color_num
##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
##  [75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [112] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

Um die Varianzen (= Eigenwerte) im Trainings-Datensatz zu ermitteln, berechnen wir einfach das Quadrat der im Objekt pca_train gespeicherten Standardabweichungen, siehe Abschnitt 5. Wir berechnen gleichzeitig die prozentualen Anteile der Hauptkomponenten an der Gesamtvarianz, wie bereits oben für den kompletten Datensatz durchgeführt:

evalues_train = pca_train$sdev^2                # eigenvalues of training set, equal to variances 
perc = evalues_train/sum(evalues_train)*100     # proportions of variances (%)

Jetzt können diese Prozentwerte auf den Achsen der Abbildung ausgewiesen werden:

xtext = paste("PC1 (", round(perc[1], 2), "% )")
ytext = paste("PC2 (", round(perc[2], 2), "% )")
plot(pca_train$x[,1:2], col=color_num, xlab=xtext, ylab=ytext, main="PCA training results", font.main=1, pch=16)
legend(-1, 2.7, c("setosa", "versicolor", "virginica"), col = c(2,3,4), pch = 16)

Mit der autoplot-Funktion von ggfortify würde man das so machen:

autoplot(pca_train, data = iris[-ind,], colour = 'Species')

Diesmal wurden nur die Trainingsdaten mittels PCA auf eine 2-dimensionale Fläche projiziert. Alle notwendigen Parameter der Projektion sind im Objekt pca_train gespeichert. Und nun der eigentliche Schritt bei der Prognose: Wir können jetzt mittels der Funktion predict genau die gleiche Projektion auf die Testdaten anwenden:

prediction = predict(pca_train, newdata = testdata)
head(prediction)
##             PC1        PC2        PC3         PC4
## 142  1.88106485  0.6642441  0.1110383  0.46972144
## 55   1.07084436 -0.2125961 -0.4114502  0.11240560
## 66   0.87458297  0.4887231 -0.5104178  0.11570398
## 39  -2.37925415 -0.9327544  0.1881522 -0.01202211
## 60   0.01930294 -1.0439913  0.5146971 -0.03074279
## 2   -2.03390204 -0.6979172 -0.2345314  0.10359394

Das Objekt prediction enthält nun 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 beide Datensätze unterscheiden zu können:

xtext = paste("PC1 (", round(perc[1], 2), "% )")
ytext = paste("PC2 (", round(perc[2], 2), "% )")
plot(pca_train$x[,1:2], col=color_num, xlab=xtext, ylab=ytext, 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 nah beieinander liegen - kein Modell ist perfekt!.


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


uwe.menzel@matstat.org