Métodos de ordenamento, PCA e NMDS
Luiz Gustavo
19/09/2022
Carregando bibliotecas e dados
library(tidyverse)
library(ggplot2)
library(factoextra)
library(vegan)
#Carregar os arquivos do Numerical Ecology
load("C:/Users/luiz8/Documents/Ciências Biológicas/2022.1/Ecologia Numérica/R/NEwR-2ed_code_data/NEwR2-Data/Doubs.RData")
#Remover dfs e dis
env <- env[,-c(1,4)]PCA
#Criação da PCA
princomp(env)-> pca_env
summary(pca_env)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 267.1096679 11.533521173 6.9246972888 4.0343429098
## Proportion of Variance 0.9972202 0.001859241 0.0006702139 0.0002274875
## Cumulative Proportion 0.9972202 0.999079484 0.9997496976 0.9999771850
## Comp.5 Comp.6 Comp.7 Comp.8
## Standard deviation 1.018804e+00 7.091927e-01 2.482529e-01 1.600339e-01
## Proportion of Variance 1.450751e-05 7.029758e-06 8.613915e-07 3.579612e-07
## Cumulative Proportion 9.999917e-01 9.999987e-01 9.999996e-01 9.999999e-01
## Comp.9
## Standard deviation 6.461191e-02
## Proportion of Variance 5.834949e-08
## Cumulative Proportion 1.000000e+00Podemos notar, então, que o PC1 explica quase que a totalidade da informação.
Histogramas PCs
Utilzando PC1 (Comp.1):
pca_env$scores %>%
as_tibble() %>%
ggplot(aes(Comp.1)) + geom_histogram()
Utilizando PC2 (Comp.2):
pca_env$scores %>%
as_tibble() %>%
ggplot(aes(Comp.2)) + geom_histogram()
Vamos dividir o rio em 2 partes, médio-baixo e alto, usando o valor 10 do PC1 como divisório, já que é aproximadamente nele que há a divisão maior entre as partes do histograma. O PC2 não parece ser muito bom para as análises, mas utilizamos o valor 2 para dividir.
Gráficos PCA
#Criar as classificações
pca_env$scores %>%
as_tibble() %>%
mutate(setor1=ifelse(Comp.1<10,"medio-baixo", "alto")) %>%
mutate(setor2=ifelse(Comp.2<2,"medio-baixo", "alto")) -> env_set
#Plotar gráfico PC1
fviz_pca_ind(pca_env,
geom.ind = "point",
col.ind = env_set$setor1,
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE,
legend.title = "Grupos"
)
#Plotar gráfico PC2
fviz_pca_ind(pca_env,
geom.ind = "point",
col.ind = env_set$setor2,
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE,
legend.title = "Grupos"
)
Como esperado, a classificação de acordo com o PC1 gerou resultados bem melhores e é possível notar uma divisão entre as 2 áreas. Fato que fica ainda mais explícito quando se observa esse gráfico:
fviz_eig(pca_env)
Que nos mostra que PC1 (comp.1) realmente tem quase 100% da informação.
NMDS
Há relação entre o NMDS e as categorias criadas anteriormente?
Para saber isso, plotaremos o gráfico NMDS codificado com cores relativas a divisão criada.
#Criação do NMDS
nmds<-metaMDS(spe)## Run 0 stress 0.001878781
## Run 1 stress 9.387515e-05
## ... New best solution
## ... Procrustes: rmse 0.004582759 max resid 0.009217972
## ... Similar to previous best
## Run 2 stress 4.723814e-05
## ... New best solution
## ... Procrustes: rmse 0.0001599194 max resid 0.0003168627
## ... Similar to previous best
## Run 3 stress 9.895516e-05
## ... Procrustes: rmse 0.0001456367 max resid 0.0002722373
## ... Similar to previous best
## Run 4 stress 6.314478e-05
## ... Procrustes: rmse 3.764033e-05 max resid 7.154917e-05
## ... Similar to previous best
## Run 5 stress 9.933156e-05
## ... Procrustes: rmse 0.000162794 max resid 0.0003328481
## ... Similar to previous best
## Run 6 stress 5.62247e-05
## ... Procrustes: rmse 2.980778e-05 max resid 5.989537e-05
## ... Similar to previous best
## Run 7 stress 9.87876e-05
## ... Procrustes: rmse 0.0002489909 max resid 0.0004973338
## ... Similar to previous best
## Run 8 stress 9.413792e-05
## ... Procrustes: rmse 0.0001395671 max resid 0.0002710208
## ... Similar to previous best
## Run 9 stress 9.755819e-05
## ... Procrustes: rmse 0.0001437138 max resid 0.0002891157
## ... Similar to previous best
## Run 10 stress 6.840587e-05
## ... Procrustes: rmse 3.775781e-05 max resid 9.112487e-05
## ... Similar to previous best
## Run 11 stress 9.871132e-05
## ... Procrustes: rmse 5.156517e-05 max resid 7.847622e-05
## ... Similar to previous best
## Run 12 stress 6.547911e-05
## ... Procrustes: rmse 4.224527e-05 max resid 8.611773e-05
## ... Similar to previous best
## Run 13 stress 7.593512e-05
## ... Procrustes: rmse 3.887498e-05 max resid 7.257405e-05
## ... Similar to previous best
## Run 14 stress 8.346256e-05
## ... Procrustes: rmse 0.0001366852 max resid 0.0002651075
## ... Similar to previous best
## Run 15 stress 9.189073e-05
## ... Procrustes: rmse 0.0001521368 max resid 0.0003010742
## ... Similar to previous best
## Run 16 stress 9.329477e-05
## ... Procrustes: rmse 5.067762e-05 max resid 9.080026e-05
## ... Similar to previous best
## Run 17 stress 9.202891e-05
## ... Procrustes: rmse 0.0001497166 max resid 0.0002957239
## ... Similar to previous best
## Run 18 stress 9.718175e-05
## ... Procrustes: rmse 0.0001567333 max resid 0.0003046979
## ... Similar to previous best
## Run 19 stress 9.825972e-05
## ... Procrustes: rmse 5.43103e-05 max resid 8.001208e-05
## ... Similar to previous best
## Run 20 stress 9.866117e-05
## ... Procrustes: rmse 0.0002461024 max resid 0.0004855225
## ... Similar to previous best
## *** Solution reached#Juntando os dados para o plot
nmds_dat<-data.frame(nmds$points, env_set$setor1)
colnames(nmds_dat) <- c("MDS1","MDS2","setor")
nmds_dat<-nmds_dat[-8,] #retirando o 8, pois o local é vazio
#Plot do gráfico
nmds_dat %>%
ggplot(aes(MDS1, MDS2, color=setor))+geom_point()
Podemos notar que os resultados se mostraram bem diferentes e que não é possível encontrar padrão entre a divisão feita pelo PCA e os pontos do MDS.
K-means
Vamos utilizar agrupamentos automáticos K-means para comparar com nossos resultados.
Primeiramente, vemos quantos grupos existem.
fviz_nbclust(env, kmeans, method = "silhouette")
Temos 2 grupos, então vamos plotar o gráfico:
km_env<-kmeans(env, centers = 2)
fviz_cluster(km_env, data = env)
Vemos que realmente existem 2 grupos, assim como sugerido e observado no PCA, apesar da falta de relação com o NMDS.