1. Vorbetrachtung

Wir werden jetzt die Hauptkomponentenanalyse (PCA) an einem 3-dimensionalen Datensatz durchführen, also diesen derart auf zwei Dimensionen projizieren, dass die größtmögliche Varianz sichtbar wird.

2. Parameter for 3D-normalverteilte Daten

Wir simulieren dreidimensionale normalverteilte Daten mit den Mittelwerten \(\mu_x = 1\) , \(\mu_y = 3\) und \(\mu_z = 2\). Die Kovarianzmatrix (\(\Sigma\)) muss symmetrisch und positiv definit sein. Letzteres wird mit Hilfe der Funktion is.positive.definite aus dem Package matrixcalc überprüft.

mu = c(1, 3, 2)
Sigma <- matrix(c(3,   1.5,    1, 
                  1.5,   2,  0.5, 
                  1,   0.5,    1), byrow = TRUE, nrow = 3)
isSymmetric(Sigma)
## [1] TRUE
library(matrixcalc) 
is.positive.definite(Sigma, tol=1e-8)  
## [1] TRUE


3. Simulation der Daten

Wir setzen einen seed um die Datenerzeugung wiederholbar zu machen. Zur Erzeugung der Daten wird mvrnorm aus dem Package MASS benutzt.

set.seed(77)
library(MASS) 
data3 = mvrnorm(n = 500, mu = mu, Sigma = Sigma)
colnames(data3) = c("x","y","z")
head(data3)
##                x        y         z
## [1,]  2.92817247 1.716572 3.1399360
## [2,] -0.62420281 2.018428 0.3835038
## [3,] -0.69263473 3.479778 1.0662602
## [4,] -0.08737466 1.367032 0.5041806
## [5,]  0.64412600 2.979015 1.7780631
## [6,] -0.52558127 1.558227 0.6274831


4. Zentrieren und Skalieren der Daten

Die Hauptkomponentenanalyse beruht im Wesentlichen auf einer Analyse der Varianz von Datensätzen entlang verschiedener Raumrichtungen. Die Daten müssen daher zuvor standardisiert werden, d.h. die Daten müssen so transformiert werden, dass jede Variable den Mittelwert 0 und die Varianz 1 hat. Dies ermöglicht, vereinfacht ausgedrückt, eine “Gleichbehandlung” aller Variablen. (Sonst könnte man ja die Varianz einer physikalischen Größe einfach vergrößern, indem man eine andere Einheit wählt, z.B. Meter statt Kilometer.)

Wir definieren zunächst eine Funktion, welche die Standardisierung durchführen kann:

zscore <- function(vector) {
  new_vector = (vector - mean(vector))/sd(vector)  
  return(new_vector)
}

Wir wenden die Transformation nun auf die Daten an und überprüfen das Ergebnis:

dfz = apply(data3, 2, zscore)
colMeans(dfz)  
##             x             y             z 
## -4.683753e-20  5.145797e-17  2.427225e-17
apply(dfz, 2, sd) 
## x y z 
## 1 1 1

Die Mittelwerte sind nahe genug bei Null, Rundungsfehler können nicht vermieden werden.


5. Daten plotten

Zum Erstellen interaktiver Grafik verwenden wir einige Funktionen aus dem Package rgl:

par3d(windowRect = c(819, 88, 1492, 706))  
plot3d(dfz[,1], dfz[,2], dfz[,3], aspect = c(1,1,1), col = "blue", size = 3, xlab="x", ylab = "y", zlab = "z")
rglwidget()

Die Abbildung kann mit Hilfe der Maus gedreht werden.
Im Zweidimensionalen hatten wir die erste Haupkomponente “experimentell” ermittelt, indem wir einen Vektor derart gesucht haben, dass die Projektion aller Datenpunkte auf diesen Vektor eine maximale Varianz hat. Als zweite Haupkomponente wurde dann ein Vektor gewählt, der ebenfalls durch den Nullpunkt geht und senkrecht auf der ersten Hauptkomponente steht.
Ein solches Experiment ist im dreidimensionalen Fall schlecht möglich. Wir wissen jedoch bereits vom zweidimensionalen Fall, dass die Eigenvektoren der Kovarianz-Matrix der standardisierten Daten identisch mit den Hauptkomponenten sind. Wir werden also im nächsten Abschnitt diese Eigenvektoren mit Hilfe von R bestimmen.


6. Eigenvektoren der Kovarianz-Matrix

Im Allgemeinen kennen wir die Kovarianz-Matrix eines Datensatzes nicht, diese muss also aus den Daten selbst ermittelt werden. (Zu beachten ist, dass die Daten - wie oben geschehen - zunächst standardisiert werden müssen !). Die Kovarianz-Matrix kann mit Hilfe der R-Funktion cov berechnet werden, während eigen zur Berechnung der Eigenwerte und Eigenvektoren dient:

covdfz = cov(dfz)       # covariance matrix
ev = eigen(covdfz)  
evalues = ev$values     # eigenvalues
evectors = ev$vectors   # eigenvectors
colnames(evectors) = c("PC1", "PC2", "PC3")
evectors
##             PC1        PC2        PC3
## [1,] -0.6248568 -0.0176850  0.7805390
## [2,] -0.5551150 -0.6929358 -0.4600948
## [3,] -0.5490002  0.7207823 -0.4231685

Wie wir bereits im 2-dimensionalen Fall gesehen haben, entsprechen die Eigenvektoren der Kovarianz-Matrix den Hauptkomponenten des Datensatzes. Wir können diese nun zum Plot der Daten hinzufügen:

par3d(windowRect = c(819, 88, 1492, 706)) 
plot3d(dfz[,1], dfz[,2], dfz[,3], aspect = c(1,1,1), col = "blue", size = 3, xlab="x", ylab = "y", zlab = "z")
x0 = 0; y0 = 0; z0 = 0    # all PC's go through zero
PC1 = evectors[,1]
PC2 = evectors[,2]
PC3 = evectors[,3]
arrow3d(p0 = c(x0, y0, z0), p1 = 2.5*PC1, s = 1/3, theta = pi/12, col = "red", thickness = 0.1)
arrow3d(p0 = c(x0, y0, z0), p1 = 2*PC2, s = 1/3, theta = pi/12, col = "green")
arrow3d(p0 = c(x0, y0, z0), p1 = 2*PC3, s = 1/3, theta = pi/12, col = "darkgoldenrod")
rglwidget()

Der erste Haupkomponenete (in rot) zeigt in die Richtung der größten Variation des Datensatzes.


7. Projektion der Daten in die Ebene der ersten beiden Hauptkomponeneten

Durch die ersten beiden Hauptkomponenten wird eine Ebene definiert. Wir werden nun alle 3-D Datenpunkte auf diese Ebene projizieren und das Ergebnis plotten. Wir sehen dann den Datensatz von der Raumrichtung, die uns die größte Ausdehnung des Datensatzes zeigt. Wie ein 3-dimensionaler Datenpunkt auf eine Ebene projiziert werden kann, hatte wir bereits in diesem Dokument gesehen. Wir werden den dort beschriebenen Algorithmus später auf alle Datenpunkte anwenden, zunächst versuchen wir es aber nur mit dem ersten Punkt:

7.1. Projektion des ersten Punktes

# the new basis vectors that define the plane: 
u = PC1 
v = PC2

# the normal to this plane:
normal = cross(u,v)
a = normal[1]
b = normal[2]
c = normal[3]

# a point in the plane (we know the plane goes through the origin):
x0 = 0
y0 = 0
z0 = 0                          
d = -(a*x0 + b*y0 + c*z0)  
d  # should be zero
## [1] 0
# some point in 3D to project (1st point of our dataset): 
point.nr = 1
x = dfz[point.nr,1]
y = dfz[point.nr,2]
z = dfz[point.nr,3]

# calculate 3D-coordinates of the projected point:
k = (d - a*x - b*y - c*z) / (a^2 + b^2 + c^2) 
x1 = x + k*a  
y1 = y + k*b  
z1 = z + k*c  
r1 = c(x1, y1,z1)   # projected point

## plot the original point and the projection:
par3d(windowRect = c(819, 88, 1492, 706)) 
plot3d(dfz[,1], dfz[,2], dfz[,3], aspect = c(1,1,1), col = "grey", size = 3, xlab="x", ylab = "y", zlab = "z")
x0 = 0; y0 = 0; z0 = 0    # all PC's go through zero
PC1 = evectors[,1]
PC2 = evectors[,2]
PC3 = evectors[,3]
arrow3d(p0 = c(x0, y0, z0), p1 = 2.5*PC1, s = 1/3, theta = pi/12, col = "red", thickness = 0.1)
arrow3d(p0 = c(x0, y0, z0), p1 = 2*PC2, s = 1/3, theta = pi/12, col = "green")
planes3d(a, b, c, -d, col = "green", alpha = 0.5, lwd = 2) 
points3d(x, y, z, col = "blue", size = 8, pch = 16)              # original point
points3d(x1, y1, z1, col = "red", size = 8, pch = 16)            # projected point
lines3d(c(x, x1),  c(y, y1), c(z, z1), col="darkgrey", lwd = 2)
rglwidget()

Dies scheint geklappt zu haben. Wir brauchen jetzt noch die Koordinaten im 2-dimensionalen Bezugssystem, welches durch PC1 und PC2 gebildet wird. Auch dies wurde schon hier gezeigt.

# get the 2D co-ordinates:
x_prime = u %*% r1
y_prime = v %*% r1
p_prime = c(x_prime, y_prime)
p_prime
## [1] -0.896832  1.477696
# plot the plane defined by PC1 and PC2 as well as the projected point:
plot(x_prime, y_prime, col = "red", pch = 18, xlab = "x", ylab = "y", main = "PC1- and PC2-plane", font.main = 1)


7.2. Projektion aller Punkte

Wir haben bereits oben die Parameter \(a\), \(b\), \(c\), \(d\), \(x_0\), \(y_0\), und \(z_0\), sowie die Hauptkomponenten \(u\) und \(v\) berechnet. In dem folgenden Plot werden nur noch die 2-dimensionalen Koordinaten im System der Hauptkomponenten \(u\) und \(v\) gezeigt:

coord2D = data.frame(x = numeric(), y = numeric())  # data frame to save the 2D-coordinates

for (point.nr in 1:nrow(dfz)) {   # run through all points 
  
  x = dfz[point.nr,1]
  y = dfz[point.nr,2]
  z = dfz[point.nr,3]
  
  # calculate projected point:
  k = (d - a*x - b*y - c*z) / (a^2 + b^2 + c^2) 
  x1 = x + k*a  
  y1 = y + k*b  
  z1 = z + k*c  
  r1 = c(x1, y1,z1)   # projected point
  
  ## coordinates in the 2D plane:
  x_prime = u %*% r1
  y_prime = v %*% r1
  p_prime = c(x_prime, y_prime)
  
  new_row = data.frame(x = x_prime, y = y_prime)
  coord2D = rbind(coord2D, new_row)  # add next point to the frame containing the results
}

plot(coord2D$x, coord2D$y, col="blue", pch=18, xlab="PC1", ylab="PC2", main="PC1 and PC2", font.main=1)

Wir wollen im nächsten Abschnitt dieses Ergebnis mit der in R eingebauten Funktion vergleichen.


8. Bedeutung der Eigenwerte der Kovarianz-Matrix

Die Eigenvektoren der Kovarianzmatrix stehen senkrecht aufeinander und zeigen, wie oben gesehen, in die Richtungen der jeweils größten Varianzen der Punktwolke. Wenn man die Koordinatenachsen in Richtung der Eigenvektoren (Hauptachsen) dreht, wird die Kovarianzmatrix diagonalisiert, d.h. dass außerhalb der Hauptdiagonale nur Nullen vorkommen. Auf der Haupdiagonale stehen dann die Eigenwerte der Kovarianzmatrix, welche mit den Varianzen entlang der durch die Eigenvektoren vorgegebenen Raumrichtungen identisch sind. Nach der Drehung des Koordinatensystems haben wir also folgenden Ausdruck für die Kovarianzmatrix:

\[ \Sigma^{\prime} = \begin{pmatrix} E_1 & 0 & 0 \\ 0 & E_2 & 0 \\ 0 & 0 & E_3 \end{pmatrix} = \begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{pmatrix} \]

Die erste Haupdiagonale hat die größte Varianz, so dass auch der erste Eigenwert der größte ist. Die Funktion eigen in R sortiert die Eigenwerte der Größe nach. Man kann nun die relativen Varianzen berechnenen, indem man die Eigenwerte auf die Summe 1 oder 100% normiert:

evalues = eigen(covdfz)$values     # eigenvalues of the covariance matrix 
evalues
## [1] 2.0532338 0.6215607 0.3252054
perc = evalues/sum(evalues)*100        # scaled eigenvalues
round(perc, 2)
## [1] 68.44 20.72 10.84

Der Plot der Projektion kann nun noch verbessert werden, indem man die relativen Varianzen einträgt:

xtext = paste("PC1 (", round(perc[1], 2), "%)")
ytext = paste("PC2 (", round(perc[2], 2), "%)")
plot(coord2D$x, coord2D$y, col="blue", pch=18, xlab=xtext, ylab=ytext, main="PC1 and PC2", font.main=1)

Wenn die Punktwolke sehr langgestreckt ist, sind die relativen Varianzen sehr unterschiedlich. Andererseits führt eine eher kugelförmige Punktwolke zu ähnlich großen Varianzen in den drei Raumrichtungen.

9. Vergleich mit der in R eingebauten Funktion

Die Funktion prcomp ist eine in R etablierte Funktion zur Hauptkomponentenanalyse. Mit den Optionen center = TRUE und scale = TRUE wird die Standardisierung veranlasst, wir können also die ursprünglichen Daten “data3” (siehe ganz oben) in die Funktion schicken:

pca <- prcomp(data3, center = TRUE, scale = TRUE) 
pca$rotation
##          PC1        PC2        PC3
## x -0.6248568  0.0176850  0.7805390
## y -0.5551150  0.6929358 -0.4600948
## z -0.5490002 -0.7207823 -0.4231685

Hier noch einmal zum Vergleich oben von uns gefundenen Hauptachsen:

df = data.frame(PC1 = PC1, PC2 = PC2, PC3 = PC3)
df
##          PC1        PC2        PC3
## 1 -0.6248568 -0.0176850  0.7805390
## 2 -0.5551150 -0.6929358 -0.4600948
## 3 -0.5490002  0.7207823 -0.4231685

Die in pca$rotation gespeicherten Vektoren entsprechen (fast) den von uns gefundenen Hauptkomponenten. Die Koordinaten des Vektors PC2 haben lediglich ein anderes Vorzeichen, der Vektor repräsentiert damit jedoch die gleiche Raumrichtung. Die gleiche Situation hatten wir schon im 2-dimensionalen Fall.

Mit dem summary-Kommando in R kann man die Standardabweichungen und die relativen Variananzen ebenfalls listen:

pca <- prcomp(data3, center = TRUE, scale = TRUE) 
summary(pca)
## Importance of components:
##                           PC1    PC2    PC3
## Standard deviation     1.4329 0.7884 0.5703
## Proportion of Variance 0.6844 0.2072 0.1084
## Cumulative Proportion  0.6844 0.8916 1.0000

Zum Schluss noch ein Vergleich der von uns hergestellten Projektion mit der im Package ggfortify angebotenen:

par(mar = c(4, 4, .2, .2))   # plot side by side 
xtext = paste("PC1 (", round(perc[1], 2), "%)")
ytext = paste("PC2 (", round(perc[2], 2), "%)")
plot(coord2D$x, coord2D$y, col="blue", pch=18, xlab=xtext, ylab=ytext)
suppressPackageStartupMessages(library(ggfortify))
autoplot(pca)

Im Vergleich zum linken Bild ist das rechte Bild, bedingt durch das andere Vorzeichen von PC2, von oben nach unten gespiegelt. Außerdem haben beide Plots unterschiedliche Skalen. Diese Unterschiede sind jedoch für unser Anliegen unwesentlich - beide Plots repräsentieren die gleichen Aussagen in Bezug auf die Varianz entlang der ersten beiden Hauptachsen.



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


uwe.menzel@matstat.org