PCA mit prcomp

Wir wollen hier (möglichst) genau untersuchen, was prcomp macht, indem wir (fast) alle Ergebnisse nachrechnen. Führen wir zunächst die PCA wieder am iris-Datensatz durch und schauen uns das Ergebnis an:

data(iris)
pca <- prcomp(iris[,-5], center = TRUE, scale = TRUE) 
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"

Die Beschreibung des Outputs von prcomp befindet sich hier.

Wir versuchem im Folgenden, alle Komponenten der Liste pca zu erklären. Die Komponenten sind:

names(pca)  #  "sdev"   "rotation"   "center"   "scale"   "x"
## [1] "sdev"     "rotation" "center"   "scale"    "x"

Alle 5 Komponenten der Liste werden wir uns jetzt näher anschauen.

Komponente $rotation

Die Haupkomponenten sind identisch mit den Eigenvektoren der Kovarianzmatrix der Daten. Bevor wir diese berechnen, zentrieren und skalieren wir die vier Variablen (= Spalten) der iris-Daten, wie oben im prcomp-Kommando:

zscore <- function(vector) (vector - mean(vector))/sd(vector)  # center and scale to mean 0 and sd 1
iris0 = apply(iris[,-5], 2, zscore)
apply(iris0, 2, sd)    # should be 1
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##            1            1            1            1
apply(iris0, 2, mean)  # should be zero, at least very small
##  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width 
## -4.484318e-16  2.034094e-16 -2.895326e-17 -2.989362e-17

Die Kovarianzmatrix wird aus den empirischen Daten mit Hilfe der Funktion cov berechnet. Die Eigenvektoren dieser Matrix werden dann mit der Funktion eigen ermittelt:

covmatrix = cov(iris0)
covmatrix  # (we see in the main diagonal that the variance of all variables is 1)
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
## Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
## Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
## Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000
evectors = eigen(covmatrix)$vectors    # eigenvectors
colnames(evectors) = c("PC1", "PC2", "PC3", "PC4")
rownames(evectors) = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
evectors
##                     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 vergleichen dieses Ergebniss mit pca$rotation:

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

In pca$rotation sind die Eigenvektoren der Kovarianzmatrix der Daten gespeichert. Damit ist bereits eine Komponente des Objektes pca erklärt.

Es sei hier noch angemerkt, dass prcomp nicht den gleichen Algorithmus wie hier angegeben benutzt, sondern aus Gründen der numerischen Stabilität Singulärwertzerlegung mit Hilfe der Funktion svd ausführt. Die Ergebnisse sind jedoch identisch.

Komponente $center

In $center sind einfach die Mittelwerten der Variablen (= Spalten) der Originaldaten gespeichert. Zum Beweis berechnen wir die Mittelwerte:

pca$center
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##     5.843333     3.057333     3.758000     1.199333
apply(iris[,-5], 2, mean)
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##     5.843333     3.057333     3.758000     1.199333

Komponente $scale

In $scale sind einfach die Standardabweichungen der Variablen der Originaldaten gespeichert. Zum Beweis berechnen wir diese:

pca$scale
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    0.8280661    0.4358663    1.7652982    0.7622377
apply(iris[,-5], 2, sd)
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    0.8280661    0.4358663    1.7652982    0.7622377

Komponente $x

Bei der Komponente $x handelt es sich um die Koordinaten der Punkte im PCA-Raum.

dim(pca$x)
## [1] 150   4
head(pca$x, 3)
##            PC1        PC2        PC3        PC4
## [1,] -2.257141 -0.4784238  0.1272796 0.02408751
## [2,] -2.074013  0.6718827  0.2338255 0.10266284
## [3,] -2.356335  0.3407664 -0.0440539 0.02828231

Wir haben die ersten beiden Spalten pca$x[,1:2] bereits hier benutzt, um folgenden Plot zu machen:

color_num = as.numeric(iris$Species) + 1  # specify colours (add 1 for nicer colours)
plot(pca$x[,1:2], col = color_num, xlab = "PC1", ylab = "PC2", main = "PCA results", font.main = 1, pch = 16)
legend("topright", c("setosa", "versicolor", "virginica"), col = unique(color_num), pch = 16)

Wie werden diese Koordinaten berechnet? Dazu müssen wir noch einmal auf den Abschnitt Theorie zurückgreifen. Die Originaldaten (iris0) werden mit Hilfe einer Drehmatrix in die PCA-Koordinaten überführt. Die Spalten der Drehmatrix bestehen aus den Eigenvektoren der Kovarianzmatrix. Letztere hatten wir schon oben erhalten:

D = evectors
D
##                     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

Ebenso gut hätten wir natürlich auch pca$rotation als Drehmatrix nehmen können, der Name kommt nicht von ungefähr:

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

Die neuen PCA-Koordinaten - nennen wir sie \(U\) - werden nun, wie schon hier erwähnt, durch Multiplikation der skalierten Originaldaten (\(X\)) mit der Drehmatrix (\(D\)) errechnet:

\[ \boldsymbol{U} = \boldsymbol{X} \cdot \boldsymbol{D} \]

Wir führen die Matrixmultiplikation in R aus:

U = iris0 %*% D

Wir vergleichen U mit den in pca$x gespeicherten Daten:

head(U)
##            PC1        PC2         PC3          PC4
## [1,] -2.257141 -0.4784238  0.12727962  0.024087508
## [2,] -2.074013  0.6718827  0.23382552  0.102662845
## [3,] -2.356335  0.3407664 -0.04405390  0.028282305
## [4,] -2.291707  0.5953999 -0.09098530 -0.065735340
## [5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
## [6,] -2.068701 -1.4842053 -0.02687825  0.006586116
head(pca$x)
##            PC1        PC2         PC3          PC4
## [1,] -2.257141 -0.4784238  0.12727962  0.024087508
## [2,] -2.074013  0.6718827  0.23382552  0.102662845
## [3,] -2.356335  0.3407664 -0.04405390  0.028282305
## [4,] -2.291707  0.5953999 -0.09098530 -0.065735340
## [5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
## [6,] -2.068701 -1.4842053 -0.02687825  0.006586116
all(abs(pca$x - U) < 1e-14) # we cannot expect exact agreement because of numerical errors
## [1] TRUE

Wir können also zusammenfassen, dass die neuen (PCA-) Koordinaten aus den Originaldaten durch Multiplikation mit der Drehmatrix enstehen. Die Spalten der Drehmatrix bestehen aus den Eigenvektoren der Kovarianzmatrix der Daten.

Komponente $sdev

Das Subset sdev enthält die Standardabweichungen entlang der Hauptkomponenten, also die Standardabweichungen der Spalten der neuen (PCA-) Daten \(\boldsymbol{U}\).

pca$sdev
## [1] 1.7083611 0.9560494 0.3830886 0.1439265
apply(U, 2, sd)  # U = rotated data
##       PC1       PC2       PC3       PC4 
## 1.7083611 0.9560494 0.3830886 0.1439265

Die Eigenwerte der Kovarianzmatrix der Originaldaten sind gleich den Varianzen der neuen Daten, deren Quadratwurzeln sind also ebenfalls die Standardabweichungen der neuen Daten:

evalues = eigen(cov(iris0))$values
sqrt(evalues)
## [1] 1.7083611 0.9560494 0.3830886 0.1439265

Die Diagonalelemente der Kovarinzmatrix der neuen Daten \(\boldsymbol{U}\) sind die Varianzen der neuen Daten, deren Quadratwurzeln sind also ebenfalls die Standardabweichungen der neuen Daten:

sqrt(diag(cov(U)))
##       PC1       PC2       PC3       PC4 
## 1.7083611 0.9560494 0.3830886 0.1439265

Zusammenssend stellen wir fest, dass folgende Vektoren identisch sind:

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


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


uwe.menzel@matstat.org