Wir haben in einem vorangegangenen Abschnitt gesehen, wie man mit der PCA Prognosen vornehmen bzw. Objekte klassifizieren kann. Hier wollen wir untersuchen, was prcomp dabei genau macht.
Als Beispieldatensatz benutzen wir auch hier die Iris-Daten.

data(iris)
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

Die letzte Spalte enthält die Bezeichnung der Art:

levels(iris[,5])
## [1] "setosa"     "versicolor" "virginica"

Wie im vorangegangenen Abschnitt teilen wir die Daten in Trainings- und Testdaten ein:

set.seed(333)
ind = sample(1:150, 15)      # find 15 indices out of 150 randomly
testdata = iris[ind,-5]      # 15 observations for testing     
traindata = iris[-ind,-5]    # 135 observations for training 
dim(testdata)
## [1] 15  4
dim(traindata)
## [1] 135   4

Wir erstellen nun ein Modell mit Hilfe der Trainingsdaten:

pca_model = prcomp(traindata, center = TRUE, scale = TRUE)  # object "pca_model" contains the model

Wie wir schon wissen, rotiert prcomp das Koordiantensystem derart, dass die größte Varianz in der ersten Hauptkomponente zu liegen kommt. Die Drehmatrix befindet sich in der Komponente $rotation:

pca_model$rotation
##                     PC1        PC2        PC3        PC4
## Sepal.Length  0.5193202 0.38712778 -0.7139349  0.2659617
## Sepal.Width  -0.2734621 0.91956001  0.2536313 -0.1236894
## Petal.Length  0.5801174 0.02665727  0.1375187 -0.8023975
## Petal.Width   0.5647910 0.06189399  0.6380100  0.5197341

Die Prognose (Klassifizierung) führt prcomp mittels der Funktion predict aus. Dabei wird das Object pca_model benutzt und auf die neuen (Test-) Daten angewendet:

prediction = predict(pca_model, newdata = testdata)
prediction
##             PC1         PC2         PC3          PC4
## 142  1.88106485  0.66424411  0.11103828  0.469721435
## 55   1.07084436 -0.21259611 -0.41145016  0.112405598
## 66   0.87458297  0.48872305 -0.51041783  0.115703977
## 39  -2.37925415 -0.93275442  0.18815215 -0.012022109
## 60   0.01930294 -1.04399134  0.51469712 -0.030742794
## 2   -2.03390204 -0.69791716 -0.23453144  0.103593941
## 134  1.10721618 -0.29851189 -0.20100793 -0.175220438
## 67   0.35235897 -0.22149888  0.46931603 -0.186963996
## 44  -1.92229166  0.41516848  0.30811352  0.175412297
## 103  2.18060259  0.54720482 -0.22115318  0.070045316
## 40  -2.12178153  0.22214132 -0.17139430  0.012018883
## 63   0.56810696 -1.73258788 -0.77958523  0.051907825
## 87   1.04431753  0.50115423 -0.40521653  0.048580920
## 75   0.70436706 -0.07336571 -0.45589768  0.052738349
## 90   0.28936330 -1.32308737  0.06812567  0.009270895

Der Output von predict besteht aus den Koordinaten der Testdaten im Koordinatensystem der Trainingsdaten. Wir werden das weiter unten plotten.

Was macht predict genau?

Zunächst werden Mittelwerte und Standardabweichungen der Testdaten skaliert, d.h. die Mittelwerte der Trainingsdaten (!) werden subtrahiert und es wird durch die Standardabweichungen der Trainingsdaten (!) geteilt. Wir behandeln die Testdaten somit in gleicher Weise wie zuvor die Trainingsdaten, denn letztere wurde ja auf genau diese Weise skaliert. Wir können dazu die Funktion scale benutzen:

testdata_scaled = scale(testdata, pca_model$center, pca_model$scale) # scale testdata

Wir schauen uns die neuen Mittelwerte und Standardabweichungen der Testdaten an:

apply(testdata_scaled, 2, mean)   # the means of the scaled testdata
##  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width 
##  0.0580453220 -0.2924928057 -0.0028885609  0.0009545726
apply(testdata_scaled, 2, sd)     # the standard deviations of the scaled testdata
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    1.0156839    0.7235240    0.8492670    0.8351556

Das gleiche hätten wir auch so machen können:

testdata_scaled_2 = matrix(NA, nrow = nrow(testdata), ncol = ncol(testdata))
for (k in 1:ncol(testdata)) testdata_scaled_2[,k] = (testdata[,k] - pca_model$center[k])/ pca_model$scale[k]
apply(testdata_scaled_2, 2, mean)
## [1]  0.0580453220 -0.2924928057 -0.0028885609  0.0009545726
apply(testdata_scaled_2, 2, sd)
## [1] 1.0156839 0.7235240 0.8492670 0.8351556

Hinweis: Das bedeutet nicht, dass die neuen Mittelwerte/Standardabweichungen der Testdaten nach dem Skalieren gleich den entsprechenden Werten der Traingsdaten sind!

## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##     5.838519     3.070370     3.758519     1.199259
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    0.8294923    0.4457216    1.7950756    0.7759920

Die skalierten Testdaten werden jetzt mit der gleichen Drehmatrix rotiert wie die Trainingsdaten, d.h. wir benutzen pca_model$rotation, um die Drehung der Testdaten auszuführen. Wir hatten das im Theorieteil schon etwas näher betrachtet.

D = pca_model$rotation     # rotation matrix which has been used to establish the model with the training data
U = testdata_scaled %*% D  # rotate with the testdata using this rotation matrix    
U
##             PC1         PC2         PC3          PC4
## 142  1.88106485  0.66424411  0.11103828  0.469721435
## 55   1.07084436 -0.21259611 -0.41145016  0.112405598
## 66   0.87458297  0.48872305 -0.51041783  0.115703977
## 39  -2.37925415 -0.93275442  0.18815215 -0.012022109
## 60   0.01930294 -1.04399134  0.51469712 -0.030742794
## 2   -2.03390204 -0.69791716 -0.23453144  0.103593941
## 134  1.10721618 -0.29851189 -0.20100793 -0.175220438
## 67   0.35235897 -0.22149888  0.46931603 -0.186963996
## 44  -1.92229166  0.41516848  0.30811352  0.175412297
## 103  2.18060259  0.54720482 -0.22115318  0.070045316
## 40  -2.12178153  0.22214132 -0.17139430  0.012018883
## 63   0.56810696 -1.73258788 -0.77958523  0.051907825
## 87   1.04431753  0.50115423 -0.40521653  0.048580920
## 75   0.70436706 -0.07336571 -0.45589768  0.052738349
## 90   0.28936330 -1.32308737  0.06812567  0.009270895

Wir können dieses Ergebnis mit dem Output von predict vergleichen:

prediction
##             PC1         PC2         PC3          PC4
## 142  1.88106485  0.66424411  0.11103828  0.469721435
## 55   1.07084436 -0.21259611 -0.41145016  0.112405598
## 66   0.87458297  0.48872305 -0.51041783  0.115703977
## 39  -2.37925415 -0.93275442  0.18815215 -0.012022109
## 60   0.01930294 -1.04399134  0.51469712 -0.030742794
## 2   -2.03390204 -0.69791716 -0.23453144  0.103593941
## 134  1.10721618 -0.29851189 -0.20100793 -0.175220438
## 67   0.35235897 -0.22149888  0.46931603 -0.186963996
## 44  -1.92229166  0.41516848  0.30811352  0.175412297
## 103  2.18060259  0.54720482 -0.22115318  0.070045316
## 40  -2.12178153  0.22214132 -0.17139430  0.012018883
## 63   0.56810696 -1.73258788 -0.77958523  0.051907825
## 87   1.04431753  0.50115423 -0.40521653  0.048580920
## 75   0.70436706 -0.07336571 -0.45589768  0.052738349
## 90   0.28936330 -1.32308737  0.06812567  0.009270895

Wir haben also das gleiche Ergebnis erreicht.

Wir können die Ergebnisse noch plotten, wie schon hier getan.

Wir plotten zunächst die ersten beiden Hauptkomponenten der Trainingsdaten, undzwar als kleine Kreise (pch = 1). Die Farben werden anhand der Species- Bezeichnung gewählt. Wir fügen die ersten beiden Hauptkomponenten der Testdaten hinzu, diese als große ausgefüllte Kreise (pch = 16, cex = 1.5), um die Punkte unterscheiden zu können:

color_num = as.numeric(iris[-ind,]$Species) + 1 # specify colours for training data
plot(pca_model$x[,1:2], col = color_num, xlab = "PC1", ylab = "PC2", main = "PCA training results", font.main = 1, pch = 1)
legend(-1, 2.7, c("setosa", "versicolor", "virginica"), col = c(2,3,4), pch = 16)
test_colors = as.numeric(iris[ind,]$Species) + 1                   # colours for the test data
points(prediction[,1:2], pch = 16, cex = 1.5, col = test_colors)   # test data

Dies scheint hinreichend gut geklappt zu haben. Die Testdaten könnnen fast alle anhand ihrer Lage im Plot der korrekten Art zugeordnet werden. Das heißt: hätten wir die wahren Arten-Bezeichnungen der Testdaten nicht gekannt, hätten wir diese nun anhand der Lage der Punkte in der Abbildung erraten können. Probleme würde es eventuell bei der Unterscheidung von versicolor und virginica geben, da diese schon in den Trainingsdaten recht nahe beieinander liegen - kein Modell ist perfekt!.


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


uwe.menzel@matstat.org