Einleitung

In Biplots werden die Beobachtungen (Observations) und die Variablen in einem gemeinsamen Hauptkomponenten-Plot (PCA-Plot) dargestellt. Es ist üblich, die Beobachtungen durch Punkte und die Variablen durch Vektoren (Pfeile) zu visualisieren. Eng beieinander liegende Punkte bedeuten dann, dass die entsprechenden Beobachtungen ähnlich sind. Andererseits kann man aus der Lage der Vektoren die Ähnlichkeit von Variablen und den Einfluss der Variablen auf die Hauptkomponenten ableiten.

PCA für die Iris-Daten

Die Iris-Daten hatten wir bereits hier beschrieben. Wir laden diese zunächst:

data(iris)
# str(iris)
dim(iris)
## [1] 150   5
head(iris, 3)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa

Wir führen die PCA in der gleichen Weise wie hier durch:

pca_results <- prcomp(iris[,-5], center = TRUE, scale = TRUE)
# str(pca_results)

Plotten der Variablen in der PCA-Ebene

Wir werden zunächst nur den Plot der Variablen machen. Das kann man ganz einfach mit dem plot-Kommnado in R erledigen, wie zum Beispiel schon hier oder hier gemacht.
An dieser Stelle benutzen wir jedoch das factoextra package:

suppressMessages(library(factoextra))

Der Variablen-Plot kann mit der Funktion fviz_pca_var gemacht werden. Wir speichern den Plot zunächst in der Variablen varplot und geben ihn dann mit den print-Kommando aus:

varplot <- fviz_pca_var(pca_results, axes = c(1, 2), geom = c("arrow", "text"), repel = TRUE, 
             col.var = "black", fill.var = "white", alpha.var = 1, title = "Die Variablen in der PCA-Ebene")
print(varplot)  # variable plot (loading plot, variable correlation plot, correlation circle)

Jede Variable wird durch einen Pfeil im Variablen-Plot charakterisiert. Aus dem Plot kann abgeleitet werden, wie stark eine Variable die Hauptkomponenten beeinflusst. Ist die Projektion eines Variablen-Vektors auf eine Achse groß, dann ist die Korrelation dieser Variable mit der Achse groß. In der Darstellung ist zu sehen, dass die Variablen Petal.Length und Petal.Width einen starken Einfluss auf die erste Hauptkomponente haben, während Sepal.Width mehr zur Konstruktion der zweiten Hauptkomponente beigetragen hat. Je näher ein Variablen-Pfeil an einer PC-Achse liegt, desto größer ist der Korrelationskoeffizient zwischen der Variable und dieser Hauptkomponente. Wir werden die Korrelationskoeffizienten weiter unten ausrechnen.

Zusätzlich kann man aus dem Plot ersehen, wie die Variablen untereinander korrelieren. Wenn zwei Variablen-Vektoren dicht beieinander liegen, also der Winkel zwischen ihnen klein ist, sind die Variablen stark positiv korreliert. Aus dem obigen Plot können wir ableiten, dass Petal.Length und Petal.Width stark positiv korreliert sind. Beträgt der Winkel zwischen zwei Pfeilen etwa 90 Grad, dann sind die entsprechenden Variablen schwach korreliert, d.h. der Korrelationskoeffizient ist nahe Null. Ein Beispiel dafür sind die Variablen Sepal.Length und Sepal.Width.
Schließlich, wenn zwei Variablen in unterschiedliche Richtungen zeigen, also der Winkel zwischen ihnen ungefähr 180 Grad beträgt, dann sind diese Variablen negativ korreliert. Ein Beispiel für eine starke negative Korrelation haben wir in der Abbildung nicht, jedoch dürften Petal.Length und Sepal.Width jedenfalls negativ korreliert sein, also einen Korrelationskoeffizienten haben, der unter Null liegt.

Im folgenden Abschnitt werden wir versuchern, all diese Aussagen quantitativ zu unterlegen, d.h. wir werden die genannten Korrelationskoeffizienten berechnen.

Korrelation zwischen Variablen und PC-Achsen

Die Koordinaten der Projektionen in die Ebene der Hauptkomponenten sind nach der PCA mit prcomp in pca_results$x gespeichert. (pca_results war der Output von prcomp, siehe oben):

dim(pca_results$x)
## [1] 150   4
head(pca_results$x, 3)
##            PC1        PC2        PC3        PC4
## [1,] -2.257141 -0.4784238  0.1272796 0.02408751
## [2,] -2.074013  0.6718827  0.2338255 0.10266284
## [3,] -2.356335  0.3407664 -0.0440539 0.02828231

Diese Koordinaten haben dieselbe Dimension wie die Originaldaten, sie stehen spaltenweise in pca_results$x. Um Schreibarbeit zu sparen, führen wir neue Variablen ein:

PC1 = pca_results$x[,1]  # first PC
PC2 = pca_results$x[,2]  # second PC
# colnames(iris) # "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"
Sepal.Length = iris[,1]  # Variables
Sepal.Width = iris[,2]
Petal.Length = iris[,3]
Petal.Width = iris[,4]

In R kann der Korrelationskoeffizient mit der Funktion cor berechnet werden.

Wie hatten bei der Diskussion der Grafik behauptet, dass die Variablen Petal.Length und Petal.Width einen starken Einfluss auf die erste Hauptkomponente haben, während Sepal.Width mehr zur Konstruktion der zweiten Hauptkomponente beigetragen hat. Das sollte sich in den Korrelationskoeffizienten widerspiegeln:

cor(Petal.Length, PC1)
## [1] 0.9915552
cor(Petal.Width, PC1)
## [1] 0.964979
cor(Sepal.Width, PC2)
## [1] -0.8827163

Unsere Aussagen haben sich bewahrheitet.
Wir wollen noch überprüfen, ob die Korrelationskoeffizienten tatsächlich identisch mit den Projektionen auf die Achsen sind. Wir benutzen dazu die Korrelationskoeffizienten als Koordinaten von Punkten, die wir dann in den Variablen-Plot eintragen. Diese Punkte sollten dann mit den Spitzen der Variablen-Vektoren übereinstimmen. Zum Beispiel sollte cor(Sepal.Width, PC1) die Projektion des Sepal.Width-Pfeils auf die PC1-Achse sein, und cor(Sepal.Width, PC2) sollte die Projektion des Sepal.Width- Pfeils auf die PC2-Achse sein. Das probieren wir sofort (mit fviz_add können Punkte zu bestehenden Grafiken hinzugefügt werden, wir hatten ja oben varplot gespeichert):

# print(varplot)  # the plot defined above
x = cor(Sepal.Width, PC1)  # use correlation coefficient as x-coordinate
y = cor(Sepal.Width, PC2)  # use correlation coefficient as y-coordinate
corr_df = data.frame(x = x, y = y)  # must be a data.frame
fviz_add(varplot, corr_df, color ="red", geom = "point", 
         pointsize = 3.0, addlabel = FALSE)  # add a point to a previously defined plot

Das scheint geklappt zu haben. Wir machen das noch für alle Variablen. Zunächst definieren wir die Koordinaten aller vier zu zeichnenden Punkte mit Hilfe der Korrelationen zu den Achsen und speichern diese in einem data.frame:

# colnames(iris) #  "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"
correlation_df = data.frame(  
  x = c(cor(PC1, Sepal.Length), 
        cor(PC1, Sepal.Width), 
        cor(PC1, Petal.Length), 
        cor(PC1, Petal.Width)),
  y = c(cor(PC2, Sepal.Length), 
        cor(PC2, Sepal.Width), 
        cor(PC2, Petal.Length), 
        cor(PC2, Petal.Width)))
colnames(correlation_df) = c("PC1", "PC2")
rownames(correlation_df) = colnames(iris[,-5])
correlation_df   # correlations of PC1 and PC2 with all 4 variables
##                     PC1         PC2
## Sepal.Length  0.8901688 -0.36082989
## Sepal.Width  -0.4601427 -0.88271627
## Petal.Length  0.9915552 -0.02341519
## Petal.Width   0.9649790 -0.06399985

Anschließend fügen wir die Punkte zu der Grafik hinzu:

fviz_add(varplot, correlation_df, color ="red", geom = "point", 
         pointsize = 3.0, addlabel = FALSE)  # add points stored in a data.frame to a previously defined plot

So weit, so gut.
Hier noch eine Grafik als Merkhilfe:

graph = "varplot_plus.png"
knitr::include_graphics(graph)

Korrelation zwischen Variablen

Bei der Diskussion des Variablen-Plots hatten wir Aussagen über den Zusammmenhang zwischen der gegenseitigen Lage der Vektoren und ihrer Korrelation gemacht. Diese wollen wir im Folgenden quantitativ untermauern.

Da der Winkel zwischen den Vektoren für Petal.Length und Petal.Width klein ist, sollten diese stark positiv korreliert sein:

cor(Petal.Length, Petal.Width) 
## [1] 0.9628654

Der Winkel zwischen Sepal.Length und Sepal.Width beträgt etwa 90 Grad, diese sollten schwach korreliert sein:

cor(Sepal.Length, Sepal.Width)
## [1] -0.1175698

Ein Beispiel für eine starke negative Korrelation haben wir in der Abbildung nicht, jedoch dürften Petal.Length und Sepal.Width negativ korreliert sein:

cor(Petal.Length, Sepal.Width)
## [1] -0.4284401

Bestimmtheitsmaß

Das Bestimmtheitsmaß gibt an, inwieweit eine abhängige Variable von anderen Variablen bestimmt wird, im Gegensatz zu zufälligen Einflüssen. Es gibt an, wie gut die abhängige Variable mit Hilfe der unabhängigen Variablen vorhergesagt werden kann. Das Bestimmtheitsmaß nimmt Werte zwischen 0 und 1 an. Ist eine abhängige Variable stark verrauscht, wird das Bestimmtheitsmaß klein sein (bei Null), ist der Zusammenhang zwischen den Variablen stark deterministisch, wird das Bestimmtheitsmaß groß sein (nahe Eins).

Geläufiger ist der englische Ausdruck “Coefficient of determination” oder “Explained variance”, oder einfach “R-squared”. Im Kontext der PCA spricht man teilweise von “Quality of representation”. Eigentlich wird der Begriff Bestimmtheitsmaß vorzugsweise im Zusammenhang mit linearer Regression gebraucht. Wie wir gleich sehen werden, besteht tatsächlich ein Zusammenhang mit linearer Regression.

Nicht ohne Grund verwendet man für das Bestimmtheitsmaß das Formelzeichen \(R^2\), denn es handelt sich tatsächlich um das Quadrat des Korrelationskoeffizienten (welcher oft mit \(r\) bezeichnet wird). Dies hatte wir schon einmal bei der linearen Regression mit einer unabhängigen Variablen gesehen. Wir können also das Bestimmtheitsmaß aus unserem oben konstruierten data.frame correlation_df berechnen:

Rsquared = correlation_df^2
Rsquared
##                    PC1         PC2
## Sepal.Length 0.7924004 0.130198208
## Sepal.Width  0.2117313 0.779188012
## Petal.Length 0.9831817 0.000548271
## Petal.Width  0.9311844 0.004095980

Im Kontext der Regression ist \(R^2\) der Anteil der Variation der abhängigen Variable, der durch die Variation der unabhängigen Variable erklärt werden kann. Auf unsere PCA übertragen können wir beispielsweise die erste Hauptkomponente als abhängige Variable und Sepal.Length als unabhängige Variable betrachten: \(R^2\) ist der Anteil der Variation von PC1, der durch die Variation von Sepal.Length erklärt werden kann. Aus der obigen Tabelle ist ersichtlich, dass diese Zahl \(0.792\) sein müsste, dass also rund \(79 \%\) der Variation von PC1 durch die Variation von Sepal.Length erklärt wird. Wir können die Regression von Sepal.Length auf PC1 einfach mit der Funktion lm durchführen und uns die erklärte Varianz ausgeben lassen:

res = lm(PC1 ~ Sepal.Length) # Linear regression 
# summary(res)               # Results
summary(res)$r.squared       # R-squared 
## [1] 0.7924004

Das ist genau die Zahl in der obigen Rsquared-Tabelle.

Noch ein Beispiel:

summary(lm(PC1 ~ Sepal.Width))$r.squared
## [1] 0.2117313

Beitrag der Variablen zur Varianz auf den Achsen

Wir haben im vorherigen Abschnitt hergeleitet, wie die einzelnen Variablen die Varianz auf den PC-Achsen beeinflussen. Daraus wird oft der prozentuale Anteil der einzelnen Variablen an der totalen Varianz berechnet und visualisiert. Wir berechnen den Anteil (“Contribution”) folgendermaßen:

Anteil einer Variable = die durch diese Variable erklärte Varianz / totale Varianz aud der ganzen PC-Achse

In R könnten wir das mit einem doppelten Loop realisieren:

contrib_mat = matrix(NA, nrow = nrow(Rsquared), ncol = ncol(Rsquared))  # init matrix with same dimensions as rsquared 
for (PC_index in 1:2) {           # loop through axis
  for (Variable_index in 1:4) {   # loop through variables 
    contrib_mat[Variable_index, PC_index] = Rsquared[Variable_index, PC_index]/colSums(Rsquared)[PC_index]*100
  }   # we multiplied with 100 in order to get percentages 
}
contrib_df = as.data.frame(contrib_mat)
colnames(contrib_df) = c("PC1", "PC2")
rownames(contrib_df) = colnames(iris)[1:4]
contrib_df  # contributions of the individual variables to the total varinciance on the PC-axis 
##                    PC1         PC2
## Sepal.Length 27.150969 14.24440565
## Sepal.Width   7.254804 85.24748749
## Petal.Length 33.687936  0.05998389
## Petal.Width  31.906291  0.44812296

Wir können Vertrauen in die Richtigkeit herstellen, indem wir z.B. das erste Element “zu Fuß” nachrechnen (wir listen noch einmal Rsquared aus):

Rsquared  # calcualted above
##                    PC1         PC2
## Sepal.Length 0.7924004 0.130198208
## Sepal.Width  0.2117313 0.779188012
## Petal.Length 0.9831817 0.000548271
## Petal.Width  0.9311844 0.004095980
0.7924004/(0.7924004 + 0.2117313 + 0.9831817 + 0.9311844) * 100  # compare with contrib_df[1,1] !
## [1] 27.15097

Die Funktion fviz_contrib aus dem factoextra package stellt dies grafisch dar, wobei die Variablen nach der Größe ihres Einflusses geordnet werden:

fviz_contrib(pca_results, choice = "var", color = "black", fill = "yellow", title = "Contributions of the variables to PC1", axis = 1)

Für mich funktioniert das leider nicht für die zweite Achse (axis = 2). Dies können wir jedoch auch selbst (wobei wir das ohne Sortieren machen):

par(mar = c(6.5, 4.1, 4.1, 2.1))   # increase lower margin for the long variable names 
maintxt = "Contributions of the variables to PC2"
barplot(contrib_df$PC2, col = "yellow", names.arg = rownames(contrib_df), las = 2, 
        main = maintxt, font.main = 1) 

par(mar = c(5.1, 4.1, 4.1, 2.1))   # reset margin to default

Plot der Observationen

Die Observationen in der PCA-Ebene haben wir schon oft grafisch dargestellt (z. B. hier, hier, oder hier). Wir wiederholen das hier mit der Funktion fviz_pca_ind aus dem factoextra package:

fviz_pca_ind(pca_results, label = "none", habillage = as.factor(iris[,5]), 
             addEllipses = FALSE, pch = 17, mean.point = FALSE, palette = c("blue", "red", "green"),
             title = "PCA of flowers")

Besonders einfach ist hier die Zuordnung der Farben zu den einzelnen Arten mit dem habillage-Parameter. Es es auch möglich. Konfidenz-Ellipsen hinzuzufügen:

fviz_pca_ind(pca_results, label = "none", habillage = as.factor(iris[,5]), 
             addEllipses = TRUE, pch = 17, palette = c("blue", "red", "green"))

Hier werden auch die Schwerpunkte der 3 Gruppen gezeichnet (mit größerem Font), was wir oben mit mean.point = FALSE unterdrückt haben.

Biplot

Der Biplot vereint, wie der Name schon andeutet, den Plot der Variablen und den Plot der Observationen in der PCA-Ebene. Wie oben schon erwähnt, deuten nah beieinander liegende Punkte auf ähnliche Observationen hin, während die Lage der Vektoren Aussagen über die Korrelation der Variablen untereinander und mit den Achsen macht.
Wir benutzen fviz_pca_biplot aus dem factoextra package:

fviz_pca_biplot(pca_results,
  axes = c(1, 2),
  # geom = c("point", "text"),
  geom = c("point"),  
  geom.var = c("arrow", "text"),
  col.ind = "black",
  fill.ind = "white",
  col.var = "black",
  fill.var = "white",
  gradient.cols = NULL,
  label = "all",
  invisible = "none",
  repel = TRUE,
  habillage = iris$Species,
  palette = c("blue", "red", "green"),
  addEllipses = FALSE, pointsize = 2.0,
  pch = 19, title = "PCA - Biplot of Iris Data", mean.point = FALSE
)

Bei großen Datemmengen, insbesondere wenn die Namen (bzw. Zeilennummern) der Observationen hinzugefügt werden, wird der Biplot oft unübersichtlich.

fviz_pca_biplot(pca_results,
  axes = c(1, 2),
  geom = c("point", "text"),
  # geom = c("point"),  
  geom.var = c("arrow", "text"),
  col.ind = "black",
  fill.ind = "white",
  col.var = "black",
  fill.var = "white",
  gradient.cols = NULL,
  label = "all",
  invisible = "none",
  repel = TRUE,
  habillage = iris$Species,
  palette = c("blue", "red", "green"),
  addEllipses = FALSE, pointsize = 2.0,
  pch = 19, title = "PCA - Biplot of Iris Data", mean.point = FALSE
)

Es scheint ratsam, die Variablen und die Beobachtungen separat darzustellen.


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

uwe.menzel@matstat.org