Bei der Hauptkomponentenanalyse wird das Koordinatensystem so gedreht, dass die 1. Hauptkomponente (x-Achse im Plot) in Richtung der größten Varianz der mehrdimensionalen Punktwolke zeigt. Die nächsten 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 der Daten gebildet werden. In den folgenden Formeln werden, um Schreibarbeit zu sparen, nur Daten mit 3 Variablen und 5 Observationen betrachtet, also z.B. ein solcher Datensatz:

set.seed(999)
X = matrix(sample(1:5, size = 15, replace = TRUE), byrow = TRUE, nrow = 5)
colnames(X) = c("var1", "var2", "var3")
rownames(X) = c("obs1", "obs2", "obs3", "obs4", "obs5")
X
##      var1 var2 var3
## obs1    3    4    5
## obs2    1    1    2
## obs3    2    3    5
## obs4    3    3    2
## obs5    1    5    3

Der Datensatz ist dreidimensional (denn er hat 3 Variablen). Er besteht aus 5 Punkten (= Observationen) in einem drei-dimensionalen Raum:

suppressWarnings(library(rgl))
plot3d(x = X[,1], y = X[,2], z = X[,3], col = "green", type = 's', 
       radius = .15, xlab = "x", ylab = "y", zlab = "z", aspect = c(1,1,1))
rglwidget()

Wir indizieren die Eigenvektoren der Kovarianzmatrix mit den Zahlen \(1\), \(2\) und \(3\). Der erste Eigenvektor lautet z.B.:

\[ \begin{pmatrix} d_{1,1} \\ d_{2,1} \\ d_{3,1} \end{pmatrix} \]

Die Drehmatrix ist damit

\[ \boldsymbol{D} = \begin{pmatrix} d_{1,1} & d_{1,2} & d_{1,3} \\ d_{2,1} & d_{2,2} & d_{2,3} \\ d_{3,1} & d_{3,2} & d_{3,3} \end{pmatrix} \]


Wir bezeichnen die neuen (“gedrehten”) Koordinaten mit \(\boldsymbol{U}\), die Originaldaten hatten wir \(\boldsymbol{X}\) genannt. In Matrixschreibweise können wir die Drehung so schreiben:

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

Dasselbe in Koordinatenschreibweise:

\[ \begin{pmatrix} u_{1,1} & u_{1,2} & u_{1,3} \\ u_{2,1} & u_{2,2} & u_{2,3} \\ u_{3,1} & u_{3,2} & u_{3,3} \\ u_{4,1} & u_{4,2} & u_{4,3} \\ u_{5,1} & u_{5,2} & u_{5,3} \end{pmatrix} = \begin{pmatrix} x_{1,1} & x_{1,2} & x_{1,3} \\ x_{2,1} & x_{2,2} & x_{2,3} \\ x_{3,1} & x_{3,2} & x_{3,3} \\ x_{4,1} & x_{4,2} & x_{4,3} \\ x_{5,1} & x_{5,2} & x_{5,3} \end{pmatrix} \cdot \begin{pmatrix} d_{1,1} & d_{1,2} & d_{1,3} \\ d_{2,1} & d_{2,2} & d_{2,3} \\ d_{3,1} & d_{3,2} & d_{3,3} \end{pmatrix} \hspace{2cm} (1) \]

So ergibt sich z.B. für \(u_{1,1}\):

\[ u_{1,1} = x_{1,1} \cdot d_{1,1} + x_{1,2} \cdot d_{2,1} + x_{1,3} \cdot d_{3,1} \hspace{2cm} (2) \]

oder für \(u_{4,2}\):

\[ u_{4,2} = x_{4,1} \cdot d_{1,2} + x_{4,2} \cdot d_{2,2} + x_{4,3} \cdot d_{3,2} \hspace{2cm} (3) \]

Allgemein kann man schreiben:

\[ u_{i,k} = \sum_l x_{i,l} \cdot d_{l,k} \hspace{2cm} (4) \]

Im Folgenden probieren wir die Drehung in R.

Zunächst zentrieren und skalieren wir die Daten, wie das im Allgemeinen auch bei der PCA gemacht wird:

zscore <- function(vector) (vector - mean(vector))/sd(vector) # scaling to mean = 0 and sd = 1 
X = apply(X, 2, zscore)
# X
apply(X, 2, mean)  # should be 0, at least very close to it
##          var1          var2          var3 
##  0.000000e+00 -1.443507e-16  9.989839e-17
apply(X, 2, sd)    # should be 1
## var1 var2 var3 
##    1    1    1

Wir berechnen nun die Eigenvektoren der Kovarianzmatrix des Beispieldatensatzes. Die Kovarianzmatrix wird aus den empirischen Daten mit Hilfe der Funktion cov berechnet. Um die Eigenwerte zu berechnen, wenden wir anschließend die Funktion eigen an:

cov(X)  # (we see in the main diagonal that the variance of all variables is 1)
##           var1      var2      var3
## var1 1.0000000 0.1685500 0.3296902
## var2 0.1685500 1.0000000 0.4000988
## var3 0.3296902 0.4000988 1.0000000
D = eigen(cov(X))$vectors    # matrix of eigenvectors = rotation matrix 
colnames(D) = c("PC1", "PC2", "PC3")
D
##            PC1         PC2        PC3
## [1,] 0.5082675  0.78485992 -0.3544843
## [2,] 0.5670023 -0.61478916 -0.5482177
## [3,] 0.6482072 -0.07764783  0.7574947

\(\boldsymbol{D}\) ist die berechntete Drehmatrix.

Um das Koordinatensystem zu drehen, führen wir nun die Matrixmultiplikation (Gleichung 1) in R aus:

# dim(X)
# dim(D)
U = X %*% D
U
##               PC1           PC2        PC3
## obs1  1.497950033  0.3713481596  0.1489927
## obs2 -1.947649031  0.1986990256  0.4683544
## obs3  0.607409690  0.0009790082  0.8730852
## obs4 -0.166568344  0.9394373284 -0.9798307
## obs5  0.008857653 -1.5104635217 -0.5106016

In der Matrix \(\boldsymbol{U}\) sind jetzt die Koordinaten im PCA-Raum enthalten.

Wie überprüfen unser Ergebnis mit dem durch prcomp erhaltenen:

pca <- prcomp(X, center = TRUE, scale = TRUE)
pca$x
##               PC1           PC2        PC3
## obs1 -1.497950033  0.3713481596 -0.1489927
## obs2  1.947649031  0.1986990256 -0.4683544
## obs3 -0.607409690  0.0009790082 -0.8730852
## obs4  0.166568344  0.9394373284  0.9798307
## obs5 -0.008857653 -1.5104635217  0.5106016

Es ist möglich, dass hier die Vorzeichen einiger Eigenvektoren genau umgekehrt sind. Dies ist kein Problem, da die Hauptachsen trotzdem die gleiche Lage haben. Dazu die prcomp-Dokumentation: “The signs of the columns of the rotation matrix are arbitrary, and so may differ between different programs for PCA, and even between different builds of R”.

Wir überprüfen noch kurz, ob auch die Ausdrücke (2) und (3) korrekt waren:

u11 = X[1,1]*D[1,1] + X[1,2]*D[2,1] + X[1,3]*D[3,1] 
u11
##     PC1 
## 1.49795
u42 = X[4,1]*D[1,2] + X[4,2]*D[2,2] + X[4,3]*D[3,2] 
u42
##       PC2 
## 0.9394373

Im neuen Koordiantensystem ist die Kovarianzmatrix dann eine Diagonalmatrix. Die neuen Daten, d.h. die Hauptkomponenten (Principal Components) sind also unkorreliert. In der Hauptdiagonale der Kovarianzmatrix stehen die Eigenwerte der ursprünglichen Kovarianzmatrix (wir nennen die Eigenwerte aus gutem Grund \(\sigma^2_{i}\)):


\[ Cov(U) = \begin{pmatrix} \sigma^2_{1} & 0 & 0 \\ 0 & \sigma^2_{2} & 0 \\ 0 & 0 & \sigma^2_{3} \end{pmatrix} \]

Auch dies können wir in R überprüfen:

round(cov(U), 10)     # covariance matrix of U; diagonal elements = variances of U = eigenvalues of X
##         PC1       PC2       PC3
## PC1 1.60849 0.0000000 0.0000000
## PC2 0.00000 0.8353561 0.0000000
## PC3 0.00000 0.0000000 0.5561537
eigen(cov(X))$values  # eigenvalues of the covariance matrix of X
## [1] 1.6084903 0.8353561 0.5561537

Die Eigenwerte der Kovarianzmatrix der Originaldaten \(\boldsymbol{X}\) sind ebenfalls identisch mit den Varianzen der neuen Daten \(\boldsymbol{U}\):

apply(U, 2, var)
##       PC1       PC2       PC3 
## 1.6084903 0.8353561 0.5561537

Desto größer der erste Eigenwert im Vergleich zu den anderen ist, desto größeren Anteil hat die erste Hauptkomponente an der Gesamtvarianz. Komponenten 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 die sogenannte explained variance der jeweiligen Hauptkomponenten berechnet. Dazu können wir benutzen, dass die Standardabweichungen in pca$sdev gespeichert sind, siehe die prcomp-Dokumentation.

pca$sdev
## [1] 1.2682627 0.9139782 0.7457571
variance_explained = (pca$sdev^2 / sum(pca$sdev^2)) * 100    # scale to percent 
variance_explained
## [1] 53.61634 27.84520 18.53846

Entsprechend dem oben Gesagten kann dasselbe auch mit Hilfe der Eigenwerte der Kovarianzmatrix berechnet werden:

eigenvalues = eigen(cov(X))$values  # eigenvalues of the covariance matrix of X
variance_explained = eigenvalues/sum(eigenvalues) * 100    # scale to percent 
variance_explained
## [1] 53.61634 27.84520 18.53846

Diese Werte werden beim Plotten der Hauptkomponenten üblicherweise an den Achsen eingetragen, so z.B. bei der autoplot- Funktion aus dem Package ggfortify:

suppressWarnings(suppressPackageStartupMessages(library(ggfortify)))
autoplot(pca, colour = "red", size = 2.5)


Anhand der Elemente \(d_{i,j}\) der Drehmatrix 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 die Varianz eines Datensatz gewinnen.

Mehr Theorie im Wikipedia-Eintrag zur Hauptkomponentenanalyse .


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


uwe.menzel@matstat.org