Search code examples
rggplot2data-visualizationpca

Building score plot using principal components

I'm trying to create score plots of the first two principal components. I begin by splitting the data into three data frames based on class. I then transform the data and perform PCA.

My data is as follows:

14      1   82.0 12.80   7.60   1070   105   400
14      1   82.0 11.00   9.00    830   145   402
14      1  223.6 17.90  10.35   2200   135   500
15      1  164.0 14.50   9.80   1946   138   500
15      1  119.0 12.90   7.90   1190   140   400
15      1   74.5  7.50   6.30    653   177   350
15      1   74.5 11.13   8.28    930   113   402
16      1  279.5 14.30   9.40   1575   230   700
16      1   82.0  7.80   6.70    676   175   525
16      1   67.0 11.00   8.30    920   106   300
16      2  112.0 11.70   8.00   1353   140   560
16      2  149.0 12.80   8.70   1550   170   550
16      2  119.0  8.50   7.40    888   175   250
16      2  119.0 13.30   9.60   1275   157   450
16      2  238.5 14.90   8.90   1537   183   700
16      2  205.0 12.00   7.90   1292   201   600
16      2   82.0  9.40   6.20    611   209   175
16      2  119.0 15.95  10.25   1350   145   450
16      2  194.0 16.74  10.77   1700   120   450
17      2  336.0 22.20  10.90   3312   135   450
17      3  558.9 23.40  12.60   4920   152   600
17      3  287.0 14.30   9.40   1510   176   800
17      3  388.0 23.72  11.86   3625   140   500
17      3  164.0 11.90   9.80    900   190   600
17      3  194.0 14.40   9.20   1665   175   600
17      3  194.0 14.40   8.90   1640   175   600
17      3  186.3  9.70   8.00   1081   205   600
17      3  119.0  8.00   6.50    625   196   400
17      3  119.0  9.40   6.95    932   165   250
17      3   89.4 14.55   9.83   1378   146   400

Column 1: type, Column 2: class, Column 3: v1, Column 4: v2, Column 5: v3, Column 6: v4, Column 7: v5, Column 8: v6

My code is as follows:

data <- read.csv("data.csv")
result <- split(data, data$class);

data1 <- result[[1]][,3:8];
data1Logged <- log10(data1)
pca.data1Logged = prcomp( ~ v1 + 
                         v2 + 
                         v3 + 
                         v4 + 
                         v5 + 
                         v6, 
                       data = data1Logged, scale. = FALSE );

data2 <- result[[2]][,3:8];
data2Logged <- log10(data2)
pca.data2Logged = prcomp( ~ v1 + 
                         v2 + 
                         v3 + 
                         v4 + 
                         v5 + 
                         v6, 
                       data = data2Logged, scale. = FALSE );

data3 <- result[[3]][,3:8];
data3Logged <- log10(data3)
pca.data3Logged = prcomp( ~ v1 + 
                         v2 + 
                         v3 + 
                         v4 + 
                         v5 + 
                         v6, 
                       data = data3Logged, scale. = FALSE );

For each of the three class, I want to have a score plot for PC1 and PC2:

pca.data1Logged$x[,1:2]
pca.data2Logged$x[,1:2]
pca.data3Logged$x[,1:2]

This is the best I could figure out:

opar <- par(mfrow = c(1,3))
plot(pca.data1Logged$x[,1:2])
plot(pca.data2Logged$x[,1:2])
plot(pca.data3Logged$x[,1:2])
par(opar)

But I would like this plot to be scaled, coloured, superimposed, etc. I have started reading about ggplot, but I don't have the experience to do this. I would like something like the following:

enter image description here

https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html

The problem with the above is that I have broken the data into 3 separate data frames, so there are no headings for "class1", "class2, "class3".

Solution

  • You can use factoextra and FactoMineR like

    library("factoextra")
library("FactoMineR")

#PCA analysis
df.pca <- PCA(df[,-c(1,2)], graph = T)
# Visualize
# Use habillage to specify groups for coloring
fviz_pca_ind(df.pca,
             label = "none", # hide individual labels
             habillage = as.factor(df$class), # color by groups
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE # Concentration ellipses, legend.title = "Class")

    enter image description here

    You can change Dim1 and 2 to PC1 and 2 manually. For that, you can note the value of "Dim1 (63.9%)" and "Dim2 (23.3%)" from this plot and use the following code to change the Dim1 and 2 to PC1 and 2 like

    fviz_pca_ind(df.pca,
             label = "none", # hide individual labels
             habillage = as.factor(df$class), # color by groups
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
             xlab = "PC1 (63.9%)", ylab = "PC2 (23.3%)", legend.title = "Class")

    enter image description here

    If you want to log transform the data, then you can use

    df[,3:8] <- log10(df[,3:8]) 

df.pca <- PCA(df, graph = T)

fviz_pca_ind(df.pca,
             label = "none", # hide individual labels
             habillage = as.factor(df$class), # color by groups
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
legend.title = "Class")

    To change Dim1 and 2 to PC1 and 2 manually, you can use the following code

    fviz_pca_ind(df.pca,
             label = "none", # hide individual labels
             habillage = as.factor(df$class), # color by groups
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
             xlab = "PC1 (64.9%)", ylab = "PC2 (22.6%)", legend.title = "Class")

    Data

    df =
structure(list(Type = c(14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 
17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L), class = c(1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), v1 = c(82, 82, 
223.6, 164, 119, 74.5, 74.5, 279.5, 82, 67, 112, 149, 119, 119, 
238.5, 205, 82, 119, 194, 336, 558.9, 287, 388, 164, 194, 194, 
186.3, 119, 119, 89.4), v2 = c(12.8, 11, 17.9, 14.5, 12.9, 7.5, 
11.13, 14.3, 7.8, 11, 11.7, 12.8, 8.5, 13.3, 14.9, 12, 9.4, 15.95, 
16.74, 22.2, 23.4, 14.3, 23.72, 11.9, 14.4, 14.4, 9.7, 8, 9.4, 
14.55), v3 = c(7.6, 9, 10.35, 9.8, 7.9, 6.3, 8.28, 9.4, 6.7, 
8.3, 8, 8.7, 7.4, 9.6, 8.9, 7.9, 6.2, 10.25, 10.77, 10.9, 12.6, 
9.4, 11.86, 9.8, 9.2, 8.9, 8, 6.5, 6.95, 9.83), v4 = c(1070L, 
830L, 2200L, 1946L, 1190L, 653L, 930L, 1575L, 676L, 920L, 1353L, 
1550L, 888L, 1275L, 1537L, 1292L, 611L, 1350L, 1700L, 3312L, 
4920L, 1510L, 3625L, 900L, 1665L, 1640L, 1081L, 625L, 932L, 1378L
), v5 = c(105L, 145L, 135L, 138L, 140L, 177L, 113L, 230L, 175L, 
106L, 140L, 170L, 175L, 157L, 183L, 201L, 209L, 145L, 120L, 135L, 
152L, 176L, 140L, 190L, 175L, 175L, 205L, 196L, 165L, 146L), 
    v6 = c(400L, 402L, 500L, 500L, 400L, 350L, 402L, 700L, 525L, 
    300L, 560L, 550L, 250L, 450L, 700L, 600L, 175L, 450L, 450L, 
    450L, 600L, 800L, 500L, 600L, 600L, 600L, 600L, 400L, 250L, 
    400L)), class = "data.frame", row.names = c(NA, -30L))