r - 使用主成分构建分数图

Question

我正在尝试创建前两个主要成分的分数图。我首先根据 class 将数据分成三个数据框。然后我转换数据并执行 PCA。

我的数据如下:

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

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

我的代码如下:

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 );

对于三个 class 中的每一个，我都想要 PC1 和 PC2 的分数图:

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

这是我能想到的最好的:

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)

但我希望对这个图进行缩放、着色、叠加等操作。我已经开始阅读有关 ggplot 的内容，但我没有这样做的经验。我想要如下内容:

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

上面的问题是我将数据分成了 3 个单独的数据框，因此没有“class1”、“class2”、“class3”的标题。

Answer 1

您可以使用 factoextra 和 FactoMineR 之类的

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")

您可以手动将 Dim1 和 2 更改为 PC1 和 2。为此，您可以从此图中记下“Dim1 (63.9%)”和“Dim2 (23.3%)”的值，并使用以下代码将 Dim1 和 2 更改为 PC1 和 2，如

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")

如果你想对数据进行日志转换，那么你可以使用

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")

要手动将 Dim1 和 2 更改为 PC1 和 2，可以使用以下代码

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")

数据

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))