Construindo plotagem de pontuação usando componentes principais

Aug 18 2020

Estou tentando criar gráficos de pontuação dos dois primeiros componentes principais. Começo dividindo os dados em três quadros de dados baseados em class. Em seguida, transformo os dados e executo o PCA.

Meus dados são os seguintes:

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

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

Meu código é o seguinte:

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

Para cada um dos três class, quero ter um gráfico de pontuação para PC1 e PC2:

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

Isso é o melhor que consegui descobrir:

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)

Mas eu gostaria que este enredo fosse escalado, colorido, sobreposto, etc. Comecei a ler sobre ggplot, mas não tenho experiência para fazer isso. Eu gostaria de algo como o seguinte:

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

O problema com o acima é que eu dividi os dados em 3 quadros de dados separados, então não há cabeçalhos para "class1", "class2, "class3".

Respostas

3 BappaDas Aug 18 2020 at 12:51

Você pode usar factoextrae FactoMineRgostar

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

Você pode alterar Dim1 e 2 para PC1 e 2 manualmente. Para isso, você pode anotar o valor de "Dim1 (63,9%)" e "Dim2 (23,3%)" deste gráfico e usar o seguinte código para alterar Dim1 e 2 para PC1 e 2 como

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

Se você deseja registrar a transformação dos dados, pode usar

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

Para alterar Dim1 e 2 para PC1 e 2 manualmente, você pode usar o seguinte código

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

Dados

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))
2 jay.sf Aug 18 2020 at 12:53

Você pode rbind os resultados separados e adicionar uma coluna de cor que você usa em plot.

rb <- rbind(cbind(pca.data1Logged$x[,1:2], d=2),
            cbind(pca.data2Logged$x[,1:2], d=3),
            cbind(pca.data3Logged$x[,1:2], d=4))

plot(rb, col=rb[,"d"], pch=20, main="PCA Plot")
legend("bottomleft", paste("data", 1:3), col=2:4, pch=20)


Dados:

data <- read.table(header=F, text="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")

names(data) <- c("sth", "class", paste0("v", 1:6))