1. Vorbetrachtung

Bei der Hauptkomponentenanalyse (engl.: PCA - Principal Component Analysis) wird ein mehrdimensionaler Datensatz vorzugsweise auf 2 oder 3 Dimensionen reduziert, um Aufschlüsse darüber zu erlangen, welche Variablen den größten Einfluß auf die Varianz einer Zielgröße haben. Liegen zweidimensionale Datensätze vor, ist daher eine Hauptkomponentenanalyse eigentlich wenig sinnvoll, da der Datensatz ja ohnehin komplett visualisiert werden kann. Trotzdem soll hier anhand eines 2D-Datensatzes das Grundprinzip der Hauptkomponentenanalyse demonstriert werden. Im Wesentlichen besteht die PCA darin, einen mehrdimensionalen Datensatz so zu rotieren, dass die größtmögliche Varianz nach Projektion auf 2 oder 3 Dimensionen sichtbar wird.

2. Parameter for 2D-normalverteilte Daten

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

mux = 1 
muy = 3
mu = c(mux, muy)    
Sigma <- matrix(c(4.0,  1.1, 1.1,  2.0), byrow = TRUE, nrow = 2)  
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) 
data2 = mvrnorm(n = 100, mu = mu, Sigma = Sigma)
data2 = data.frame(x = data2[,1], y = data2[,2])
head(data2)
##            x        y
## 1  1.7432282 4.197688
## 2 -1.9204242 3.889230
## 3 -0.1751774 2.306802
## 4 -1.3663458 2.890529
## 5 -0.1383358 4.684735
## 6 -1.3559925 2.368549


4. Zentrieren und Skalieren der Daten

Die Hauptkomponentenanalyse beruht im Wesentlichen auf einer Analyse der Varianz von Datensätzen entlang verschiedener Raumrichtungen (siehe unten). 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. 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(data2, 2, zscore)
colMeans(dfz)  
##            x            y 
## 1.932482e-17 5.141720e-17
apply(dfz, 2, sd) 
## x y 
## 1 1

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


5. Daten plotten

plot(dfz[,1], dfz[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")


6. Projektion von Punkten auf eine Gerade

Bei der Hauptkompontenanalyse wird eine Gerade gesucht, die derart ausgerichtet ist, dass die Projektionen aller Punkte auf diese Gerade eine maximale Varianz haben. Dies ist dann die erste Hauptkomponente. Zur Veranschaulichung entwerfen wir zunächst eine Funktion, die die Projektion eines Punktes \((x,y)\) auf eine beliebige Gerade - charakterisiert durch Anstieg (“slope”) und Achsenabschnitt (“intercept”) - zeichnen kann.

plotProjection <- function(slope, intercept, x, y, col = "darkgreen", lwd = 1.3) {  
  slope1 = -1/slope  
  intercept1 = y - slope1*x
  xs = (intercept1-intercept)/(slope-slope1)                  
  ys = slope*(intercept1-intercept)/(slope-slope1)+intercept    
  segments(x, y, xs, ys, col = col, lwd = lwd)
}

Die Projektion steht senkrecht auf der gewählten Geraden, daher sind die Anstiege der Geraden und der Projektion durch folgende Beziehung verknüpft:

\[ m_{Projektion} = \frac{-1}{m_{Gerade}} \] Bemerkung: Zum Verständnis des obigen Ausdrucks bedenke man, dass der Anstieg einer Geraden gleich dem Tangens des Steigungswinkels ist, und dass allgemein gilt:

\[\tan{(\phi+\pi/2)}=\frac{-1}{\tan{\phi}}\]

Die Funktion soll für den Punkt \(x_{25}\) und die Gerade \(y = 0.7 \cdot x\) probiert werden (es wird 0 als Achsenabschnitt gewählt, da die Daten standardisiert sind).

plot(dfz[,1], dfz[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
intercept = 0
slope = 0.7
abline(a = intercept, b = slope, lwd = 1.8, col = "red")  
x = dfz[25,1]
y = dfz[25,2]
points(x, y, col="green", pch = 18, cex = 1.5)
plotProjection(slope, intercept, x, y, col = "darkgreen", lwd = 1.3)

Bemerkung: Beide Achsen des Plots müssen unbedingt gleiche Skalen haben, um das Bild nicht zu verfälschen. Wir wählen daher asp = 1 im plot-Kommando.

Dies scheint geklappt zu haben. Wir können nun die Projektionen aller Punkte zeichnen:

plot(dfz[,1], dfz[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
intercept = 0
slope = 0.7
abline(a = intercept, b = slope, lwd = 1.8, col = "red") 
for (i in 1:nrow(dfz)) {
  x = dfz[i,1]
  y = dfz[i,2]
  plotProjection(slope, intercept, x, y)
}


7. Berechnung der Varianz der Projektionen

Wie bereits oben erwähnt, wird bei der Hauptkompontenanalyse zunächst eine Gerade gesucht, die derart ausgerichtet ist, dass die Projektionen aller Punkte auf diese Gerade eine maximale Varianz haben. Die empirische Varianz für \(n\) Datenpunkte \(x_i\) mit dem Mittelwert \(\mu\) wird folgendermaßen berechnet:

\[ S^2 = \frac{1}{n-1} \cdot \sum_{i=1}^n \left( x_i - \mu \right)^2 \]

Da die Daten bereits zentriert wurden (\(\mu = 0\)), vereinfacht sich die Formel zu

\[ S^2 = \frac{1}{n-1} \cdot \sum_{i=1}^n x_i^2 \]

Die Funktion proVar soll diese empirische Varianz berechnen, d.h. die Varianz der Projektionen aller Datenpunkte auf eine vorgegebene Gerade. Dazu werden noch zwei Hilfsfunktionen benötigt:
- normvec berechnet die Norm (die Länge) eines Vektors
- scalar berechnet das Skalarprodukt zweier Vektoren (die einfach mit dem c-Operator definiert werden können)

normvec <- function(vec) {  
  sum = 0
  for (i in 1:length(vec)) {
    sum = sum + vec[i]^2 
  }
  return(sqrt(sum))
}
scalar <- function(vec1, vec2) sum(vec1*vec2)

Die Gerade wird wieder durch slope und intercept definiert.

proVar<- function(points, slope, intercept, scale = FALSE) { 
  if(scale) points = apply(points, 2, zscore)   
  if(sum(colMeans(points) > 1.0e-10) > 0) {stop("Data points do not seem to be centered.\n")}  
  if(sum((apply(points, 2, sd) - 1) > 1.0e-10) > 0)  {stop("Data points do not seem to be scaled.\n")}      
  xs = 1   #  define some point on the line y = slope*x + intercept
  ys = slope*xs + intercept
  u = c(xs, ys)
  v = u/normvec(u)  # this vector of length one goes along the line
  S2 = 0
  for (i in 1:nrow(points)) {  
    x = c(points[i,1], points[i,2])  
    yi = scalar(x,v) 
    S2 = S2 + yi^2
  }
  n = nrow(points)  
  return(S2/(n-1))
}

Zur Erklärung: Wir haben hier einen auf Eins normierten Vektor \(v\) definiert, der genau in Richtung der Geraden zeigt. In dem Loop wird dann das Skalarprodukt dieses Vektors mit jedem der auf die Originalpunkte zeigenden Vektoren gebildet, also \(\vec{x_i} \cdot \vec{v}\) für alle \(i\). Dieses Skalarprokukt ist jeweils genau die Länge der Projektion des Datenpunktes \(x_i\) auf die durch slope und intercept definierte Gerade (siehe Anhang).

Ein kleiner Test mit unserem Datensatz dfz:

proVar(dfz, slope = 0.7, intercept = 0) 
## [1] 1.472078

Die ist also die Varianz der Projektionen aller in dfz gespeicherten Punkte auf die Gerade \(y = 0.7 \cdot x\).


8. Gerade mit der größten Varianz der Projektionen

Nun wollen wir die Gerade suchen, für die diese Varianz am größten wird: Wir drehen einfach eine durch den Nullpunkt gehende Gerade und registrieren für jede Neigung die empirische Varianz. Wir definieren zunächst die Funktion getVar und rufen diese dann mit verschiedenen Drehwinkeln \(\phi\) (also verschiedenen Anstiegen) auf. Die Variable speed bestimmt, wie schnell die Drehung vonstatten geht. Die Idee für eine solche Simulation stammt von dieser Webseite.

getVar <- function(phi, speed) {
  slope = tan(phi)   
  if (slope > 30) return(NA)   # phi = 90 and phi = 270 result in infinite slope since cos(phi)=0
  angle = 180*phi/pi
  S2 = proVar(dfz, slope = slope, intercept = 0)
  txt = paste("Angle = ", angle, " Slope = ", round(slope,3), " Variance of projection = ", signif(S2,4))
  varRow = data.frame(angle = angle, slope = slope, variance = S2) 
  varTable <<- rbind(varTable, varRow)  # <<- makes varTable known in the parent scope
  plot(dfz[,1], dfz[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
  abline(h=0, lty = 2, col ="darkgrey")
  abline(v=0, lty = 2, col ="darkgrey")
  abline(a = 0, b = slope, lwd = 1.8, col = "red") 
  for (i in 1:nrow(dfz)) {
    x = dfz[i,1]
    y = dfz[i,2]
    plotProjection(slope, intercept = 0, x, y)
  }
  mtext(side = 3, txt) 
  Sys.sleep(1/speed)  
} 
varTable <- data.frame(angle = numeric(0), slope = numeric(0), variance = numeric(0))
for (phi in seq(0, pi, pi/60)) getVar(phi, speed = 3)

In der Funktion getVar haben wir die Varianz der Projektionen der Punkte in Abhängigkeit vom Neigungswinkel und von Anstieg regisriert (Frame varTable), d.h. wir können jetzt ermitteln, welche Gerade die größte Varianz liefert:

# print(varTable)
plot(varTable$angle, varTable$variance, col = "red", pch = 18, xlab = "angle", ylab = "Variance of projection")  

plot(varTable$slope, varTable$variance, col = "red", pch = 18, xlab = "slope", ylab = "Variance of projection", xlim = c(-3,3)) 
opt_slope = varTable$slope[which(varTable$variance == max(varTable$variance))] 
opt_slope
## [1] 1
abline(v = opt_slope, lty = 3, col = "darkgrey")


9. Konstruktion der Hauptkomponenten

Der Anstieg 1 (Winkel 45 Grad) führt zur höchsten Varianz der Projektionen, wie aus den beiden Plots ersichtlich ist. Damit ist die gesuchte optimale Gerade einfach \(y = x\). (Das ist kein Wunder, da ganz am Anfang die Kovarianz-Matrix \(\Sigma\) dementsprechend gewählt wurde!)
Wir definieren nun einen Vektor, der in die Richtung dieser Geraden zeigt, indem wir einfach einen Punkt auf der Geraden wählen (die ja durch den Nullpunkt geht). Anschließend normieren wir den Vektor:

u = c(1,1)
u = u/normvec(u)
normvec(u)  # Norm 
## [1] 1
u
## [1] 0.7071068 0.7071068

Dieser Vektor ist die 1. Hauptkomponte. Er repräsentiert die Richtung, in der die Punktmenge am meisten variiert. Wollte man wirklich eine eindimensionale Repräsentation der Daten haben, würde man diese Richtung wählen und die Projektionen aller Punkte darauf einzeichnen. (Wie bereits erwähnt, ist das nicht sinnvoll, da die zweidimensionale Darstellung genügend aussagt - wir werden jedoch bei höherdimensionalen Daten das gleiche Prinzip anwenden.)

Eine zweidimensionale Punktmenge hat zwei Haupkomponenten. Die zweite Hauptkomponenete muss senkrecht auf der ersten stehen, d.h. deren Anstieg ist \(a = -1\) (das Produkt beider Anstiege muss \(-1\) sein, da die Geraden senkrecht aufeinander stehen, siehe dazu die Bemerkung in Abschnitt 6). Da auch die 2. Hauptkomponente durch den Ursprung geht, ist es wieder leicht, einen Vektor zu konstruieren, der in ihre Richtung zeigt und die Norm \(1\) besitzt, indem man zunächst eine beliebigen Punkt auf dieser Geraden wählt:

v = c(1,-1)
v = v/normvec(v)
normvec(v)  # Norm 
## [1] 1
v
## [1]  0.7071068 -0.7071068

Wir können nun noch die Daten zusammen mit beiden Hauptkomponenten zeichnen:

plot(dfz[,1], dfz[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 1.5)
abline(a = 0, b = -1, col = "red", lty = 2, lwd = 1.5)


10. Vergleich mit einer 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 data2 in die Funktion schicken:

pca <- prcomp(data2, center = TRUE, scale = TRUE) 
summary(pca)
## Importance of components:
##                           PC1    PC2
## Standard deviation     1.2257 0.7054
## Proportion of Variance 0.7512 0.2488
## Cumulative Proportion  0.7512 1.0000
pca$rotation
##          PC1        PC2
## x -0.7071068  0.7071068
## y -0.7071068 -0.7071068

Fällt etwas auf?
Die in pca$rotation gespeicherten Vektoren entsprechen (fast) den von uns gefundenen Hauptkomponenten. Die Koordinaten des Vektors PC1 haben lediglich beide ein negatives Vorzeichen, der Vektor repräsentiert jedoch ebenfalls die Gerade \(y = x\), genau wie der von uns gefundene.


11. Vergleich mit den Eigenvektoren der Kovarianz-Matrix

Zum Schluss sollen die Eigenvektoren der Kovarianz-Matrix der standardisierten Daten berechnet werden. Die Kovarianz 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
evectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

Die Eigenvektoren der Kovarianz-Matrix entsprechen ebenfalls den zuvor berechneten Vektoren.


12. Zusammenfassung

Wir haben:

  1. die erste Hauptkomponente eines Datensatzes experimentell gefunden, indem wir einen (normierten) Vektor ermittelt haben, der so positioniert ist, dass die Projektionen aller Datenpunkte auf diesen Vektor eine maximale Varianz haben. Die zweite Hauptkomponente wurde so gelegt, dass sie senkrecht auf der ersten steht (und ebenfalls normiert ist)
  2. die Funktion prcomp benutzt, um die Haupkomponenten zu bestimmen.
  3. die Eigenvektoren der empirischen Kovarianz-Matrix berechnet.

Alle 3 Vektorenpaare sind identisch (bis auf die Vorzeichen der mit prcomp berechneten 1. Hauptkomponente). Im Theorie-Teil werden wir sehen, wie das alles zusammenhängt.


13. Anhang


Länge der Projektion

Wir wollen einen Ausdruck für die Länge der Projektion eines Punktes \(\vec{x}\) auf einen Vektor \(\vec{v}\) finden, d.h. wir suchen die Größe \(a\) in der folgenden Grafik:

Das Skalarprodukt von \(\vec{x}\) und \(\vec{v}\) ist

\[\begin{align*} \vec{x} \cdot \vec{v} &= |\vec{x}| \cdot |\vec{v}| \cdot \cos{\theta} \\ &= |\vec{x}| \cdot |\vec{v}| \cdot \frac{a}{|\vec{x}|} \\ &= |\vec{v}| \cdot a \end{align*}\]

Wir bekommen also:

\[ a = \frac{\vec{x} \cdot \vec{v}}{|\vec{v}| } \] Wenn der Vektor \(\vec{v}\) die Norm 1 hat, vereinfacht sich dies zu

\[ a = \vec{x} \cdot \vec{v} \] Der Skalar \(a\) ist die Länge der Projektion des Vektors \(\vec{x}\) auf die Gerade \(\vec{v}\). Diese Formel wird in der Funktion proVar benutzt.


Die Norm eines Vektors

Die Norm eines n-dimensionalen Vektors ist seine Länge im n-dimensionalen Raum. Sei \(\vec{u}\) ein Vektor mit den Komponenten \(u_1, u_2, \dots u_n\):

\[ \begin{align} \vec{u} &= \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix} \end{align} \]

Dann ist seine Norm

\[ \lVert \vec{u} \rVert = \sqrt{u_1^2 + u_2^2 + \ldots + u_n^2} \]

Im Dreidimensionalen wird das

\[ \lVert \vec{u} \rVert = \sqrt{u_1^2 + u_2^2 + u_3^2} \] … und im Zweidimensionalen:

\[ \lVert \vec{u} \rVert = \sqrt{u_1^2 + u_2^2} \] Dies ist der Satz des Pythagoras.

Die Norm kann auch mit der R-internen Funktion norm berechnet werden. Allerdings muss der Input ein Matrix-Objekt sein. Der Einfachheit halber wurde daher hier eine eigene Funktion (normvec) definiert. Der Input-Vektor in normvec kann einfach mit den c-Operator generiert werden.


Die Funktion prcomp in R

Die R-Funktion prcomp führt eine Hauptkomponentenanalyse vom mehrdimensionalen Daten durch. Die offizielle Dukumentation ist hier.


Eigenwerte und Eigenvektoren

Die Eigenvektoren einer Matrix \(\boldsymbol{A}\) sind die nicht-trivialen (\(\vec{x} \ne 0\)) Lösungen der Gleichung

\[ \boldsymbol{A} \cdot \vec{x} = \lambda \cdot \vec{x} \] Diese Gleichung ist im Allgemeinen nur für bestimmte Werte von \(\lambda\) (die Eigenwerte) lösbar. Die Eigenwerte ergeben sich aus der Gleichung

\[ \det(\boldsymbol{A} - \lambda \cdot \boldsymbol{I}) = 0 \] wobei \(\boldsymbol{I}\) die Einheitsmatrix ist. Mehr dazu bei Wikipedia



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


uwe.menzel@matstat.org