PCA e NMDS
Débora Gênesis
04/12/2021
Imagine que você precisa entender como se agrupam as comunidades do rio que estamos estudando. Essa necessidade vem do fato de que um usina hidroelétrica está para ser construída no rio e o governo encomendou estudos de impacto ambiental que precisam demonstrar várias coisas. Você como biólogo, que está preocupado com o rio, precisa convencer o governo de construir a hidroelétrica num setor do rio que tenha o menor impacto possível. Para isso, você vai precisar analisar tanto a comunidade de peixes quanto as condições fisico-químicas do rio.
Passo 1 - Faça uma PCA pra a tabela “env” (excluindo as colunas ‘dfs’ e ‘dis’). Seu intuito em vez de dividir o rio de maneira arbitrária é encontrar divisões naturais de características ambientais, que possam ser separadas pelo seu valor médio. Examine os PCs 1 e 2 e veja, através de histogramas em quantas porções o rio pode ser dividido.
Passo 2 - Faça um NMDS para os dados da comunidade biológica “spe” e use as mesmas categorias que você encontrou no passo anterior para pintar as amostras. Será que há uma coincidência?
Passo 3 - Use o métodos de K-means para gerar agrupamentos automáticos e veja se eles coincidem com os que você gerou usando o método anterior.
Carregando os dados :
library(vegan)## Carregando pacotes exigidos: permute## Carregando pacotes exigidos: lattice## This is vegan 2.5-7library(FactoMineR)## Warning: package 'FactoMineR' was built under R version 4.1.2library(factoextra)## Warning: package 'factoextra' was built under R version 4.1.2## Carregando pacotes exigidos: ggplot2## Warning: package 'ggplot2' was built under R version 4.1.2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBaload ("C:/Users/debor/OneDrive/Documentos/EcoNum/github/NEwR-2ed_code_data/NEwR2-Data/Doubs.RData")
spe # matriz de abundância de espécies de peixes## Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl
## 1 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 5 4 3 0 0 0 0 0 0 0 0 0 0 0
## 3 0 5 5 5 0 0 0 0 0 0 0 0 0 1 0
## 4 0 4 5 5 0 0 0 0 0 1 0 0 1 2 2
## 5 0 2 3 2 0 0 0 0 5 2 0 0 2 4 4
## 6 0 3 4 5 0 0 0 0 1 2 0 0 1 1 1
## 7 0 5 4 5 0 0 0 0 1 1 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 9 0 0 1 3 0 0 0 0 0 5 0 0 0 0 0
## 10 0 1 4 4 0 0 0 0 2 2 0 0 1 0 0
## 11 1 3 4 1 1 0 0 0 0 1 0 0 0 0 0
## 12 2 5 4 4 2 0 0 0 0 1 0 0 0 0 0
## 13 2 5 5 2 3 2 0 0 0 0 0 0 0 0 0
## 14 3 5 5 4 4 3 0 0 0 1 1 0 1 1 0
## 15 3 4 4 5 2 4 0 0 3 3 2 0 2 0 0
## 16 2 3 3 5 0 5 0 4 5 2 2 1 2 1 1
## 17 1 2 4 4 1 2 1 4 3 2 3 4 1 1 2
## 18 1 1 3 3 1 1 1 3 2 3 3 3 2 1 3
## 19 0 0 3 5 0 1 2 3 2 1 2 2 4 1 1
## 20 0 0 1 2 0 0 2 2 2 3 4 3 4 2 2
## 21 0 0 1 1 0 0 2 2 2 2 4 2 5 3 3
## 22 0 0 0 1 0 0 3 2 3 4 5 1 5 3 4
## 23 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 24 0 0 0 0 0 0 1 0 0 2 0 0 1 0 0
## 25 0 0 0 0 0 0 0 0 1 1 0 0 2 1 0
## 26 0 0 0 1 0 0 1 0 1 2 2 1 3 2 1
## 27 0 0 0 1 0 0 1 1 2 3 4 1 4 4 1
## 28 0 0 0 1 0 0 1 1 2 4 3 1 4 3 2
## 29 0 1 1 1 1 1 2 2 3 4 5 3 5 5 4
## 30 0 0 0 0 0 0 1 2 3 3 3 5 5 4 5
## Rham Legi Scer Cyca Titi Abbr Icme Gyce Ruru Blbj Alal Anan
## 1 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 1 0 0 0 0 0 0 0
## 5 0 0 2 0 3 0 0 0 5 0 0 0
## 6 0 0 0 0 2 0 0 0 1 0 0 0
## 7 0 0 0 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0 0 0 0
## 9 0 0 0 0 1 0 0 0 4 0 0 0
## 10 0 0 0 0 0 0 0 0 0 0 0 0
## 11 0 0 0 0 0 0 0 0 0 0 0 0
## 12 0 0 0 0 0 0 0 0 0 0 0 0
## 13 0 0 0 0 0 0 0 0 0 0 0 0
## 14 0 0 0 0 0 0 0 0 0 0 0 0
## 15 0 0 0 0 1 0 0 0 0 0 0 0
## 16 0 1 0 1 1 0 0 0 1 0 0 0
## 17 1 1 0 1 1 0 0 0 2 0 2 1
## 18 2 1 0 1 1 0 0 1 2 0 2 1
## 19 2 1 1 1 2 1 0 1 5 1 3 1
## 20 3 2 2 1 4 1 0 2 5 2 5 2
## 21 3 2 2 2 4 3 1 3 5 3 5 2
## 22 3 3 2 3 4 4 2 4 5 4 5 2
## 23 0 0 0 0 0 0 0 0 1 0 2 0
## 24 0 1 0 0 0 0 0 2 2 1 5 0
## 25 0 0 1 0 0 0 0 1 1 0 3 0
## 26 2 2 1 1 3 2 1 4 4 2 5 2
## 27 3 3 1 2 5 3 2 5 5 4 5 3
## 28 4 4 2 4 4 3 3 5 5 5 5 4
## 29 5 5 2 3 3 4 4 5 5 4 5 4
## 30 5 3 5 5 5 5 5 5 5 5 5 5spa # matriz espacial## X Y
## 1 85.678 20.000
## 2 84.955 20.100
## 3 92.301 23.796
## 4 91.280 26.431
## 5 92.005 29.163
## 6 95.954 36.315
## 7 98.201 38.799
## 8 99.455 46.406
## 9 109.782 55.865
## 10 130.641 66.576
## 11 142.748 81.258
## 12 147.270 85.839
## 13 156.817 89.516
## 14 159.435 92.791
## 15 150.820 91.084
## 16 132.662 87.956
## 17 128.298 93.918
## 18 130.560 102.446
## 19 128.459 105.428
## 20 114.862 103.129
## 21 97.163 90.245
## 22 88.200 86.373
## 23 79.596 83.508
## 24 74.753 78.734
## 25 67.146 74.683
## 26 53.770 71.598
## 27 43.637 68.673
## 28 30.514 61.166
## 29 20.495 43.848
## 30 0.000 41.562env # matriz ambiental## dfs ele slo dis pH har pho nit amm oxy bod
## 1 0.3 934 48.0 0.84 7.9 45 0.01 0.20 0.00 12.2 2.7
## 2 2.2 932 3.0 1.00 8.0 40 0.02 0.20 0.10 10.3 1.9
## 3 10.2 914 3.7 1.80 8.3 52 0.05 0.22 0.05 10.5 3.5
## 4 18.5 854 3.2 2.53 8.0 72 0.10 0.21 0.00 11.0 1.3
## 5 21.5 849 2.3 2.64 8.1 84 0.38 0.52 0.20 8.0 6.2
## 6 32.4 846 3.2 2.86 7.9 60 0.20 0.15 0.00 10.2 5.3
## 7 36.8 841 6.6 4.00 8.1 88 0.07 0.15 0.00 11.1 2.2
## 8 49.1 792 2.5 1.30 8.1 94 0.20 0.41 0.12 7.0 8.1
## 9 70.5 752 1.2 4.80 8.0 90 0.30 0.82 0.12 7.2 5.2
## 10 99.0 617 9.9 10.00 7.7 82 0.06 0.75 0.01 10.0 4.3
## 11 123.4 483 4.1 19.90 8.1 96 0.30 1.60 0.00 11.5 2.7
## 12 132.4 477 1.6 20.00 7.9 86 0.04 0.50 0.00 12.2 3.0
## 13 143.6 450 2.1 21.10 8.1 98 0.06 0.52 0.00 12.4 2.4
## 14 152.2 434 1.2 21.20 8.3 98 0.27 1.23 0.00 12.3 3.8
## 15 164.5 415 0.5 23.00 8.6 86 0.40 1.00 0.00 11.7 2.1
## 16 185.9 375 2.0 16.10 8.0 88 0.20 2.00 0.05 10.3 2.7
## 17 198.5 349 0.5 24.30 8.0 92 0.20 2.50 0.20 10.2 4.6
## 18 211.0 333 0.8 25.00 8.0 90 0.50 2.20 0.20 10.3 2.8
## 19 224.6 310 0.5 25.90 8.1 84 0.60 2.20 0.15 10.6 3.3
## 20 247.7 286 0.8 26.80 8.0 86 0.30 3.00 0.30 10.3 2.8
## 21 282.1 262 1.0 27.20 7.9 85 0.20 2.20 0.10 9.0 4.1
## 22 294.0 254 1.4 27.90 8.1 88 0.20 1.62 0.07 9.1 4.8
## 23 304.3 246 1.2 28.80 8.1 97 2.60 3.50 1.15 6.3 16.4
## 24 314.7 241 0.3 29.76 8.0 99 1.40 2.50 0.60 5.2 12.3
## 25 327.8 231 0.5 38.70 7.9 100 4.22 6.20 1.80 4.1 16.7
## 26 356.9 214 0.5 39.10 7.9 94 1.43 3.00 0.30 6.2 8.9
## 27 373.2 206 1.2 39.60 8.1 90 0.58 3.00 0.26 7.2 6.3
## 28 394.7 195 0.3 43.20 8.3 100 0.74 4.00 0.30 8.1 4.5
## 29 422.0 183 0.6 67.70 7.8 110 0.45 1.62 0.10 9.0 4.2
## 30 453.0 172 0.2 69.00 8.2 109 0.65 1.60 0.10 8.2 4.4Fazendo um PCA para as variáveis ambientais excluindo “dfs” e “dis”:
env_simp<-env[,-c(1,4)] # excluindo "dfs" e "dis"
princomp(env_simp)-> 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+00Já que o componente 1 concentra a maioria dos dados, plotamos um histograma.
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()pca_env$scores %>%
as.tibble() %>%
ggplot(aes(Comp.1))+geom_histogram()## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
O histograma mostra valores positivos e negativos, dando uma ideia de divisão de valores. Usamos essa ideia para dividir o rio em dois grupos :o Alto Rio e o Médio-Baixo Rio.
pca_env$scores %>%
as.tibble() %>%
mutate(setor=ifelse(Comp.1<0,"Médio-Baixo Rio", "Alto Rio"))-> env_set
fviz_pca_ind(pca_env,
geom.ind = "text", # ver as amostras
col.ind = env_set$setor, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Agrupamentos"
)
Testando uma outra forma de visualizar essas variáveis:
Gerando um PCA :
res.pca<- PCA(env_simp)## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
O primeiro gráfico mostra a distribuição das amostras em relação às variáveis do segundo gráfico, onde é possível observar quais as variáveis tiveram maior influência nos locais amostrados. O segundo gráfico traz informações importantes de correlação entre as variáveis como por exemplo, a correlação negativa já observada entre oxy e bod. Esse gráfico ainda traz outras informações interessantes como a aparente correlação positiva entre a amônia, o fósforo e o nitrato. Se formos combinar esses dois gráficos teremos esse tipo de gráfico:
fviz_pca(res.pca)
Aqui fica mais claro às quais variáveis os locais de amostragem estão submetidos, onde podemos visualizar uma alta heterogeineidade ambiental.
Em um outro gráfico de hierarquia das amostras em relação às variáveis podemos visualizar três clusters bem definidos :
No cluster 1 por exemplo temos amostras com altos valores para ele (elevação) e slo (declive) e baixos valores para har (calcium concentration-hardness, concentração de cálcio-dureza) e nit (nitrato). No cluster 2 temos altos valores para a variável har e baixos para ele. Já no cluster 3 as variáveis mais significativas (da maior para a menor) são bod, pho, amn e nit, com os menores valores para oxy e ele.
Apesar dessa análise trazer três clusters, não existem aparentemente diferenças gritantes entre o entre os clusters 2 e 3, havendo até partes onde eles se intersectam : Ambos apresentam baixos valores para ele e valores próximos para har, por exemplo. O que faz ter sentido a primeira análise que dividiu o rio em duas partes, Alto Rio e Médio-Baixo Rio.
Analisando agora a matriz de espécies do rio :
nmds<-metaMDS(spe[-8,])## Run 0 stress 0.07477809
## Run 1 stress 0.1128294
## Run 2 stress 0.07477801
## ... New best solution
## ... Procrustes: rmse 0.0001865905 max resid 0.0008933482
## ... Similar to previous best
## Run 3 stress 0.09288418
## Run 4 stress 0.117344
## Run 5 stress 0.1133686
## Run 6 stress 0.1203765
## Run 7 stress 0.08797362
## Run 8 stress 0.1226315
## Run 9 stress 0.0880154
## Run 10 stress 0.07477824
## ... Procrustes: rmse 0.0002743227 max resid 0.001314125
## ... Similar to previous best
## Run 11 stress 0.07477829
## ... Procrustes: rmse 0.0003001347 max resid 0.001438945
## ... Similar to previous best
## Run 12 stress 0.1125692
## Run 13 stress 0.07478419
## ... Procrustes: rmse 0.003654012 max resid 0.01453342
## Run 14 stress 0.1222462
## Run 15 stress 0.07506666
## ... Procrustes: rmse 0.01475416 max resid 0.06363312
## Run 16 stress 0.07429346
## ... New best solution
## ... Procrustes: rmse 0.0238784 max resid 0.09223832
## Run 17 stress 0.07506669
## Run 18 stress 0.07506671
## Run 19 stress 0.08797356
## Run 20 stress 0.1221737
## *** No convergence -- monoMDS stopping criteria:
## 19: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrminnmds$points## MDS1 MDS2
## 1 -1.755186345 0.85754723
## 2 -1.122441494 -0.24958002
## 3 -0.978481611 -0.27813191
## 4 -0.617315547 -0.11860201
## 5 0.072590749 0.44662406
## 6 -0.423165043 -0.17083916
## 7 -0.866329876 -0.28234888
## 9 -0.007637838 -0.87244213
## 10 -0.527345780 -0.41367157
## 11 -1.066088458 -0.04946677
## 12 -0.974735714 -0.13234358
## 13 -1.157212489 0.10455014
## 14 -0.799570811 0.10931496
## 15 -0.498053065 0.16054758
## 16 -0.186300728 0.26959115
## 17 0.084055662 0.12202834
## 18 0.295952809 0.11236282
## 19 0.449559301 0.14254262
## 20 0.769125673 0.27433703
## 21 0.855498121 0.38503357
## 22 0.944242178 0.44977139
## 23 0.768651364 -1.43281673
## 24 1.129475633 -0.62031852
## 25 0.867674694 -0.86662838
## 26 0.931321146 0.13251282
## 27 0.966574239 0.37868238
## 28 1.008287364 0.38272306
## 29 0.784551356 0.56532825
## 30 1.052304511 0.59369227
## attr(,"centre")
## [1] TRUE
## attr(,"pc")
## [1] TRUE
## attr(,"halfchange")
## [1] TRUE
## attr(,"internalscaling")
## [1] 1.020891nmds_dat<-data.frame(nmds$points, env_set$setor[-8])
colnames(nmds_dat) <- c("MDS1","MDS2","setor")
nmds_dat %>% #Plotando o gráfico
ggplot(aes(MDS1, MDS2, color=setor))+geom_point()
Esse gráfico mostra que existe uma padronização na diversidade de espécies que depende das variáveis fisico-químicas do rio. Aparentemente existe especificidade de habitat (se observar ambos os gráficos dos Agrupamentos-PCA e ggplot das espécies) o que pode dificultar a implantação de uma hidroelétrica nesse ambiente. Se a necessidade de implantação não puder ser substituída recomendaria que fosse feita na parte do médio-baixo rio, sendo melhor no Baixo Rio, já que esse é o local com um maior grupo de espécies, e possui uma maior extensão fluvial.
Verificando o mesmo dado de espécies com o k-means:
k3<-kmeans(spe[,-8], centers = 3, nstart=25) #excluindo outras variáveis que não vão ser analisadas
k3## K-means clustering with 3 clusters of sizes 8, 10, 12
##
## Cluster means:
## Cogo Satr Phph Babl Thth Teso Chna Lele Sqce Baba
## 1 0.000000 0.125 0.375 1.0 0.125 0.125000 1.625 2.25 3.125000 3.7500000
## 2 0.200000 0.800 1.400 1.7 0.200 0.400000 0.500 1.30 1.700000 0.8000000
## 3 1.083333 4.000 4.250 4.0 1.000 1.166667 0.000 1.00 1.166667 0.4166667
## Albi Gogo Eslu Pefl Rham Legi Scer Cyca
## 1 2.12500000 4.3750000 3.25 2.7500000 3.5 3.00000000 2.125 2.62500000
## 2 0.90000000 1.2000000 0.80 1.0000000 0.5 0.40000000 0.400 0.30000000
## 3 0.08333333 0.6666667 0.50 0.3333333 0.0 0.08333333 0.000 0.08333333
## Titi Abbr Icme Gyce Ruru Blbj Alal Anan
## 1 4.0000000 3.125 2.25 4.125 4.8750000 3.625 5.0 3.0
## 2 0.8000000 0.100 0.00 0.500 2.2000000 0.200 1.7 0.3
## 3 0.4166667 0.000 0.00 0.000 0.1666667 0.000 0.0 0.0
##
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 2 3 3 3 2 3 3 2 2 3 3 3 3 3 3 3 2 2 2 1 1 1 2 2 2 1
## 27 28 29 30
## 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 172.7500 305.9000 189.0833
## (between_SS / total_SS = 65.5 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"#Fazendo clusters
fviz_cluster(k3, data = spe[,-8])
#Verificando o número de clusters
fviz_nbclust(spe[,-8], kmeans, method = "silhouette") 
#Plotando novamente com o número indicado
k2<-kmeans(spe[,-8], centers =2, nstart=25)
fviz_cluster(k2, data = spe[,-8])
Já o K-means traz um tipo de agrupamento diferente. Ao que se parece o Alto-Rio se junta ao Médio, fazendo apenas um cluster. O Baixo-Rio se separa totalmente fazendo um cluster bem menor. Ao que essa análise indica, independemente das variáveis físico-químicas existe um maior intercâmbio entre as espécies do Alto(A) e Médio(M) Rio, que as do Médio e Baixo (B) Rio. Havendo a necessidade de reconsiderar as colocações anteriores quanto a posição de um usina hidroelétrica.
Será que isolar o grupo do alto rio é a melhor abordagem? Aparentemente esse seria um impacto bem maior que se considerarmos isolar o cluster menor na análise do K-means. No entanto analisando a riqueza e abundância de espécies (1 aula) temos que os setores M e B do Rio possuem uma maior representatividade que o Alto Rio, fazendo os serem mais importantes. Mas se levarmos em conta a integridade do rio e como isso pode afetar as espécies Rio-Abaixo, não é recomendável que a usina seja construída próximo à nascente. Pensando em um menor impacto ambiental é preferível que seja construída no Baixo Rio que em qualquer outro setor.
Gráficos da Aula 1 para comparação:
sit.pres <- apply(spe > 0, 1, sum)
sort(sit.pres)## 8 1 2 23 3 7 9 10 11 12 13 4 24 25 6 14 5 15 16 26 30 17 20 22 27 28
## 0 1 3 3 4 5 5 6 6 6 6 8 8 8 10 10 11 11 17 21 21 22 22 22 22 22
## 18 19 21 29
## 23 23 23 26par(mfrow = c(1, 2))
#Riqueza de espécies e locais de amostragem
plot(sit.pres,type = "s",las = 1,col = "gray", main = "Riqueza de espécies vs. \n Gradiente Rio Acima-Rio Abaixo", xlab = "Número dos locais",ylab = "Riqueza de espécies")
text(sit.pres, row.names(spe), cex = .8, col = "brown")
#Fazendo um mapa de bolhas com coordenadas geográficas
plot(spa, asp = 1, main = "Mapa de riqueza de espécies",pch = 21,col = "white",bg = "brown", cex = 5 * sit.pres / max(sit.pres), xlab = "Coordenadas x (km)",ylab = "Coordenadas Y (km)")
lines(spa, col = "light blue")