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 der Anstieg 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 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 + intercept
  u = c(xs, ys)
  v = u/normvec(u)
  if(normvec(v) - 1 > 1.0e-10) stop("Vector could not be normalized")  
  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))
}

Die Länge der Projektion eines Datenpunktes \(x_i\) auf einen (auf Eins normierten) Vektor \(v\) ist gerade das Skalarprodukt \(\vec{x_i} \cdot \vec{v}\) ( siehe Anhang).

Ein kleiner Test:

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)