MODULO 5 - Estatística Multivariada - PERMANOVA
Leonardo Reffatti
15 de maio de 2019
Aula 5.23
PERMANOVA - teste de permutação de uma matriz de distância envolvendo uma variável categórica
Objetivo: Aplicar o teste estatístico PERMANOVA para dar rigor estatístico aos os agrupamentos formados por meio de uma variável categórica e visualizados através de um NMDS/PCA/Cluster de uma matriz de distância. A lógica da PERMANOVA é a de permutar a matriz de distância, gerando os valores das variáveis categóricas ao acaso. Portanto se o resultado da matriz original for muito improvável de ser encontrado ao acaso, rejeitamos a hipótese nula, e aceitamos que de fato existe uma associação entre as categorias (fatores/classes) com a matriz de distâncias visualizados através do NMDS/Cluster/PCA.
Este material está disponível em: http://rpubs.com/leonardoreffatti.
PERMANOVA da matriz de distâncias de composição de espécies dos 2 tipos de Ambiente. Definir as variáveis componentes das espécies, formar a matriz de distâncias pelo método “Jaccard”, guardar como resultado a função NMDS desta matriz de distâncias, plotar em um gráfico para visualizar a formação de grupos conforme a variável Ambiente. Com a função adonis() realizar o teste PEMANOVA de permutação da matriz de distância para dar rigor estatístico aos agrupamentos formados.
library(permute)
library(lattice)
library(vegan)## This is vegan 2.5-2setwd("C:/R/Curso do R/MODULO_5.2")
#inserindo os nomes das UAs, argumento row.names
dados<-read.table("insetos.txt", h=T, row.names = 1)
attach(dados)
str(dados)## 'data.frame': 30 obs. of 40 variables:
## $ Ambiente : Factor w/ 2 levels "Mata_Primária",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gado : Factor w/ 2 levels "Ausente","Presente": 2 1 2 2 2 1 1 1 2 2 ...
## $ Temperatura: num 21.8 21.3 23.8 23.2 21.5 21.7 20.3 19.7 24.7 20.2 ...
## $ Cobertura : int 73 62 57 49 59 78 64 74 48 75 ...
## $ Luz : num 53.5 61 71.5 69.5 61.5 53 58 56 71 60.5 ...
## $ Flores : int 14 28 17 29 15 24 21 29 20 5 ...
## $ sp.1 : int 4 5 1 3 5 1 1 0 6 3 ...
## $ sp.2 : int 0 0 0 0 0 0 3 1 4 0 ...
## $ sp.3 : int 7 3 2 6 7 2 2 1 8 3 ...
## $ sp.4 : int 2 7 0 4 2 4 5 3 3 1 ...
## $ sp.5 : int 0 2 2 1 1 0 3 4 0 0 ...
## $ sp.6 : int 0 5 2 14 4 3 3 9 8 0 ...
## $ sp.7 : int 3 2 3 2 1 4 3 2 1 3 ...
## $ sp.8 : int 2 4 3 5 2 3 0 0 1 0 ...
## $ sp.9 : int 1 2 0 0 0 0 2 0 1 0 ...
## $ sp.10 : int 1 3 5 2 1 1 0 0 1 1 ...
## $ sp.11 : int 0 0 0 4 0 2 0 2 0 0 ...
## $ sp.12 : int 1 5 0 2 0 2 4 5 0 0 ...
## $ sp.13 : int 3 6 0 3 3 2 3 1 4 2 ...
## $ sp.14 : int 3 4 2 10 5 2 2 5 8 0 ...
## $ sp.15 : int 0 0 0 0 1 0 0 0 1 0 ...
## $ sp.16 : int 3 3 0 2 0 0 3 0 1 0 ...
## $ sp.17 : int 0 0 0 1 0 0 0 0 0 0 ...
## $ sp.18 : int 0 1 2 3 1 3 0 0 2 0 ...
## $ sp.19 : int 0 0 2 2 3 0 0 2 0 1 ...
## $ sp.20 : int 0 3 0 0 0 4 0 2 0 0 ...
## $ sp.21 : int 0 1 3 4 2 1 3 3 3 2 ...
## $ sp.22 : int 1 2 1 2 2 2 1 2 2 2 ...
## $ sp.23 : int 1 0 1 2 1 1 2 1 2 1 ...
## $ sp.24 : int 1 1 0 2 1 2 0 1 1 1 ...
## $ sp.25 : int 0 1 0 2 0 1 1 0 0 1 ...
## $ sp.26 : int 2 2 1 0 1 1 2 1 1 1 ...
## $ sp.27 : int 0 0 0 0 0 0 3 0 0 0 ...
## $ sp.28 : int 1 0 0 0 0 0 0 0 0 2 ...
## $ sp.29 : int 0 0 0 0 1 0 0 0 0 0 ...
## $ sp.30 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ sp.31 : int 0 5 0 0 0 0 1 0 0 0 ...
## $ sp.32 : int 1 1 1 0 1 1 1 0 1 1 ...
## $ sp.33 : int 0 0 0 0 1 0 0 0 1 0 ...
## $ sp.34 : int 0 0 0 0 1 0 0 0 0 0 ...especies <- dados[,7:40]
#binary = TRUE, considerar dados presença e ausência
dist.jac <- vegdist(especies, method = "jaccard", binary = TRUE)
#o resultado do NMDS apresenta valor de stress, e deve ser < 0.2
resultado <- metaMDS(dist.jac, trymax = 100)## Run 0 stress 0.1837046
## Run 1 stress 0.227572
## Run 2 stress 0.1920033
## Run 3 stress 0.2005835
## Run 4 stress 0.2130398
## Run 5 stress 0.1856681
## Run 6 stress 0.1791314
## ... New best solution
## ... Procrustes: rmse 0.0466243 max resid 0.1923374
## Run 7 stress 0.2032267
## Run 8 stress 0.2035104
## Run 9 stress 0.1837048
## Run 10 stress 0.1844605
## Run 11 stress 0.2015696
## Run 12 stress 0.233344
## Run 13 stress 0.2026423
## Run 14 stress 0.1985781
## Run 15 stress 0.1832239
## Run 16 stress 0.2398749
## Run 17 stress 0.1825411
## Run 18 stress 0.216
## Run 19 stress 0.1791306
## ... New best solution
## ... Procrustes: rmse 0.0003722673 max resid 0.001301903
## ... Similar to previous best
## Run 20 stress 0.1837046
## *** Solution reached#Através da visualização no NMDS podemos ver a separação das UA por meio da variável Ambiente. "Separação visual subjetiva".
plot(resultado, type = "t")## species scores not available#adicionando polígonos na análise de ordenação
ordihull(resultado, groups = Ambiente, show = "Mata_Primária", col="green4")
ordihull(resultado, groups = Ambiente, show = "Mata_Secundária", col="blue")
#Temos a variável Ambiente
#Será que o padrão de distâncias dentro das classes ("Mata_Primária" e "Mata_Secundária") é menor que entre as classes? (PERMANOVA = ANOVA por Permutação)
#Se os padrões de grupos visualizados no NMDS forem reais. Isto implica que a distância média entre UA do mesmo grupo deveria ser menor que a distância média de UA de um grupo para o outro.
#A lógica da PERMANOVA é a de permutar a matriz de distância, gerando os valores das classes ao acaso. Portanto se o resultado da matriz original for muito improvável de ser encontrado ao acaso, rejeitamos a hipótese nula, e aceitamos que de fato existe uma associação entre as classes ou grupos com a matriz de distâncias. Como o NMDS é a representação da matriz de distâncias, pode-se utilizar o NMDS para representar os resultados da PERMANOVA, o mesmo vale para um Cluster.
#PERMANOVA = ANOVA adaptada para uma Matriz de Distâncias
#Pode-se utilizar para Matrizes de Distâncias criadas pelos Métodos: Euclidiano, Jaccard e Bray-Curtis.
adonis(dist.jac~Ambiente, permutations = 10000)##
## Call:
## adonis(formula = dist.jac ~ Ambiente, permutations = 10000)
##
## Permutation: free
## Number of permutations: 10000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Ambiente 1 0.83305 0.83305 11.044 0.28287 9.999e-05 ***
## Residuals 28 2.11195 0.07543 0.71713
## Total 29 2.94499 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#A diferença é que este output traz o valor de R2. Este R2(correlação), se fosse 1 as classes fazem grupos perfeitos. Como o resultado foi de 0.28, significa que apesar de haver a separação dos grupos, a "Mata_Primária" é bastante heterogênia, baixando o valor de correlação.