Der Iris-Datensatz

Wir benutzen den Iris-Datensatz (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.Length) sowie die Breite und Länge des Blütenblatts (Petal.Width, Petal.Length) 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, 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"

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.

Korrelation der Variablen

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(iris[,-5])

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)


PCA mit prcomp

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)

Alle Ergebnisse sind nun im Objekt pca gespeichert.

Wir berechnen den Anteil der Varianz der einzelnen Komponenten an der Gesamtvarianz (explained variance).

sd = pca$sdev  # standard deviation calculated by prcomp
var_explained = sd^2 / sum(sd^2)
round(var_explained, 4)
## [1] 0.7296 0.2285 0.0367 0.0052
sum(var_explained)  # must be 1
## [1] 1

Die berechnete explained variance kann man mit dem summary-Kommando bestätigen:

# class(pca)                 # get the class of the pca object ==> "prcomp"
# methods(class = "prcomp")  # find out which functions can be applied to pca (of class "prcomp")
summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

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 colours)
# color_num

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(var_explained[1]*100, 2), "% )")
ytext = paste("PC2 (", round(var_explained[2]*100, 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 = unique(color_num), pch = 16)

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 (mit colour = "Species"):

suppressWarnings(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(var_explained[1]*100, 2), "% )")
ytext = paste("PC2 (", round(var_explained[2]*100, 2), "% )")
ztext = paste("PC3 (", round(var_explained[3]*100, 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 mit der Maus vergrößert/verkleinert oder gedreht werden.


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


uwe.menzel@matstat.org