Ziel

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

Simulation der 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)
Sigma
##      [,1] [,2] [,3]
## [1,]  3.0  1.5  1.0
## [2,]  1.5  2.0  0.5
## [3,]  1.0  0.5  1.0
isSymmetric(Sigma)
## [1] TRUE
library(matrixcalc) 
is.positive.definite(Sigma, tol=1e-8)  
## [1] TRUE

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")
rownames(data3) = paste0("point", 1:nrow(data3))
head(data3, 3)
##                 x        y         z
## point1  2.9281725 1.716572 3.1399360
## point2 -0.6242028 2.018428 0.3835038
## point3 -0.6926347 3.479778 1.0662602


Zentrieren und Skalieren der Daten

Die Hauptkomponentenanalyse beruht auf einer Analyse der Varianz von Punktwolken entlang verschiedener Raumrichtungen. Die Daten müssen zunächst standardisiert werden, d.h. die Daten müssen so transformiert werden, dass jede Variable den gleichen Mittelwert und die gleiche Varianz 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) # scaling to mean = 0 and sd = 1

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

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

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


Daten plotten

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

par3d(windowRect = c(819, 88, 1492, 706))  
plot3d(x = points[,1], y = points[,2], z = points[,3], col = "green", type = 's', 
       radius = .07, xlab = "x", ylab = "y", zlab = "z", aspect = c(1,1,1))
rglwidget()

In 2D hatten wir die erste Haupkomponente “experimentell” ermittelt, indem wir eine Gerade gesucht haben, so dass die Projektionen aller Datenpunkte auf diese Gerade maximale Varianz haben. Ein solches Experiment ist in 3D 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 bestimmen.


Eigenvektoren der Kovarianz-Matrix

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

covmat = cov(points)    # covariance matrix
covmat
##           x         y         z
## x 1.0000000 0.6030284 0.5890165
## y 0.6030284 1.0000000 0.3786146
## z 0.5890165 0.3786146 1.0000000
ev = eigen(covmat)  
evalues = ev$values     # eigenvalues
evectors = ev$vectors   # eigenvectors
colnames(evectors) = c("PC1", "PC2", "PC3")
rownames(evectors) = c( "x", "y", "z")
evectors
##          PC1        PC2        PC3
## x -0.6248568 -0.0176850  0.7805390
## y -0.5551150 -0.6929358 -0.4600948
## z -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(x = points[,1], y = points[,2], z = points[,3], col = "green", type = 's', 
       radius = .07, xlab = "x", ylab = "y", zlab = "z", aspect = c(1,1,1))
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 = "blue", thickness = 0.1)
arrow3d(p0 = c(x0, y0, z0), p1 = 2.5*PC2, s = 1/3, theta = pi/12, col = "red", thickness = 0.1)
arrow3d(p0 = c(x0, y0, z0), p1 = 2.5*PC3, s = 1/3, theta = pi/12, col = "darkgoldenrod", thickness = 0.1)
rglwidget()

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


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 hier und hier gesehen. Wir werden den dort beschriebenen Algorithmus später auf alle Datenpunkte anwenden, zunächst versuchen wir es aber nur mit dem ersten Punkt:

Projektion des ersten Punktes

Wir benutzen die erste und zweite Hauptkomponente als Basis für die Ebene, in die wir die Punkte projizieren wollen:

u = PC1  # the new basis vectors that define the plane (1st and 2nd principal components)
v = PC2

Wir berechnen die Flächennormale:

n = cross(u,v)  # the normal vector of this plane is the cross product
nx = n[1]
ny = n[2]
nz = n[3]       # nx, ny, and nz define the normal vector

Wir geben einen Punkt auf der Ebene an (wir wissen, dass die Ebene durch den Punkt \((0, 0, 0)\) geht, da die Daten normiert wurden).

x0 = 0  # a point on the plane (we know that the plane goes through the origin)
y0 = 0
z0 = 0
r0 = c(x0, y0, z0)

Wir wählen den zu projizierenden Punkt aus:

point.nr = 1
x1 = points[point.nr, 1] 
y1 = points[point.nr, 2]
z1 = points[point.nr, 3]
r1 = c(x1, y1, z1)      # some point in 3D to project (1st point of our dataset) 
r1
## [1]  1.1651346 -0.8979782  1.2154292

Wir berechnen die Hilfsvariablen \(d\) und \(k\), wie hier und hier gezeigt:

d = -(nx*x0 + ny*y0 + nz*z0)
k = (d + nx*x1 + ny*y1 + nz*z1) / (nx^2 + ny^2 + nz^2)

Nun können wir die Koordinaten des projizierten Punktes berechnen:

xp = x1 - k*nx  
yp = y1 - k*ny  
zp = z1 - k*nz  
rp = c(xp, yp, zp)  
rp   # coordinates of the projected point
##          y          y          y 
##  0.5342585 -0.5261034  1.5574580

Wir plotten nun die Ebene, den Originalpunkt und den projizierten Punkt:

par3d(windowRect = c(819, 88, 1492, 706)) 
plot3d(points[,1], points[,2], points[,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]
arrow3d(p0 = c(x0, y0, z0), p1 = 2.5*PC1, s = 1/3, theta = pi/12, col = "blue", thickness = 0.1, width = 1/3)
arrow3d(p0 = c(x0, y0, z0), p1 = 2.5*PC2, s = 1/3, theta = pi/12, col = "red", thickness = 0.1, width = 1/3)
planes3d(nx, ny, nz, d, col = "yellow", alpha = 0.7, lwd = 2) 
points3d(x1, y1, z1, col = "blue", size = 8, pch = 16)           # original point
points3d(xp, yp, zp, col = "red", size = 8, pch = 16)            # projected point
lines3d(c(x1, xp),  c(y1, yp), c(z1, zp), 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 in the (u,v)-system:
xq = u %*% (rp - r0)
yq = v %*% (rp - r0)
rq = c(xq, yq)
rq
## [1] -0.896832  1.477696

Projektion aller Punkte

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

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

for (point.nr in 1:nrow(points)) {   # run through all points 
  
  x1 = points[point.nr, 1]
  y1 = points[point.nr, 2]
  z1 = points[point.nr, 3]
  
  # calculate projected point:
  d = -(nx*x0 + ny*y0 + nz*z0)
  k = (d + nx*x1 + ny*y1 + nz*z1) / (nx^2 + ny^2 + nz^2) 
  xp = x1 - k*nx  
  yp = y1 - k*ny  
  zp = z1 - k*nz  
  rp = c(xp, yp, zp)   # projected point
  
  ## coordinates in the 2D plane:
  xq = u %*% (rp - r0)
  yq = v %*% (rp - r0)
  rq = c(xq, yq)
  
  new_row = data.frame(xq = xq, yq = yq)
  coord2D = rbind(coord2D, new_row)  # add next point to the frame containing the results
}

plot(coord2D$xq, coord2D$yq, col = "blue", pch = 16, xlab = "PC1", ylab = "PC2", main = "PC1 and PC2", font.main = 1, cex = 0.7)

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


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^2 & 0 & 0 \\ 0 & \sigma_2^2 & 0 \\ 0 & 0 & \sigma_3^2 \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 bzw. 100% normiert:

evalues = eigen(covmat)$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$xq, coord2D$yq, col = "blue", pch = 16, xlab = xtext, ylab = ytext, main = "PC1 and PC2", font.main = 1, cex = 0.7)

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.

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
## x -0.6248568 -0.0176850  0.7805390
## y -0.5551150 -0.6929358 -0.4600948
## z -0.5490002  0.7207823 -0.4231685

Die in pca$rotation gespeicherten Vektoren entsprechen den von uns gefundenen Hauptkomponenten.

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

Die relativen Varianzen (Proportion of Variance) stimmen mit den oben aus den Eigenwerten ermittelten überein.
Die relativen Varianzen (“variance explained”) können auch direkt aus dem Output von prcomp berechnet werden:

pca$sdev
## [1] 1.4329110 0.7883912 0.5702679
variance_explained <- (pca$sdev^2 / sum(pca$sdev^2)) * 100
variance_explained
## [1] 68.44113 20.71869 10.84018

Zum Schluss noch ein Vergleich zwischen der von uns hergestellten Projektion und der von prcomp berechneten. Die Koordinaten der prcomp-Projektion sind in pca$x gespeichert. Wir vergleichen ein paar Zeilen mit unserem Frame coord2D:

# str(pca)
head(pca$x, 3)
##               PC1        PC2        PC3
## point1 -0.8968320  1.4776958 -0.8082568
## point2  1.8152660 -0.6462040 -0.2549418
## point3  0.8682801 -0.8887834  0.5602160
head(coord2D, 3)
##           xq         yq
## 1 -0.8968320  1.4776958
## 2  1.8152660 -0.6462040
## 3  0.8682801 -0.8887834

Wir können nun beide Datensätze in einem Plot zeichnen:

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$xq, coord2D$yq, col = "blue", pch = 16, xlab = xtext, ylab = ytext, cex = 0.7)
points(pca$x[,1], pca$x[,2], col = "red", pch = 4) 
legend("topright", c("eigen","prcomp"), col = c("blue", "red"), pch = c(16,4), cex = 1.2)

Beide Plots sind identisch.

Noch genauer untersuchen wir prcomp hier.



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


uwe.menzel@matstat.org