Curso de R no IPÊ

Script do dia 5.3: NMDS, adonis e envfit

Marcos V. C. Vital e participantes da turma de 2014

Para mais material, acesse http://marcosvital.wordpress.com/
 

Este documento é a 3a e última parte do quinto (e último!) de uma série que nasceu durante as atividades do “Curso de introdução ao uso do software R e suas aplicações estatísticas nas ciências biológicas”, realizado no IPÊ em dezembro de 2014.

Nesta última parte do último guia (nãããão, esta acabando…) nós vamos executar um escalonamento multidimensional não métrico (mais conhecido como NMDS), e explorar um pouco os dados relacionados a ele para brincar com as funções adonis e envfit.

Como sempre, todo o documento foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais, acesse:

http://www.r-project.org/

http://www.rstudio.com

http://rmarkdown.rstudio.com/
 

Antes de começar, carregue os dados e os pacotes necessários:

setwd("C:/R/Curso IPÊ")
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 ...

 


O NMDS

Não tem mistério, basta criar uma matriz de distância com a métrica desejada e utilizá-la dentro da função metaMDS:

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.0-10
matriz.bray<-vegdist(dados[,7:40], method="bray")
nmds.bray<-metaMDS(matriz.bray)
## Run 0 stress 0.1436911 
## Run 1 stress 0.163002 
## Run 2 stress 0.2239961 
## Run 3 stress 0.1909071 
## Run 4 stress 0.1465913 
## Run 5 stress 0.1436911 
## ... procrustes: rmse 1.935704e-05  max resid 6.765763e-05 
## *** Solution reached

Podemos plotar de maneira bem direta, mas sem nenhum elemento para interpretar:

plot(nmds.bray)
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available

Mas não tem muita graça. Fica um pouco mais informativo se usarmos os nomes das unidades amostrais:

plot(nmds.bray, type="t")
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available

  • O argumento type="t" foi usado para se chamar os nomes.

Também podemos plotar com cores relacionadas à uma variável categórica:

plot(nmds.bray, type="n")
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available
points(nmds.bray, pch=16, col=Ambiente)

  • O argumento type="n"foi usado para chamarmos uma janela vazia, e depois adicionar os pontos.

E podemos plotar os pontos e depois criar os polígonos relacionados à uma variável categórica:

plot(nmds.bray, type="n")
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available
points(nmds.bray, pch=16)
ordihull(nmds.bray, group=Ambiente, show="Mata_Primária", col="green")
ordihull(nmds.bray, group=Ambiente, show="Mata_Secundária", col="red")

 

A função adonis

A função adonis implementa uma ANOVA realizada com permutações, apropriada para analisarmos uma matriz de distância no lugar de uma variável resposta quantitativa. Podemos aplicá-la aos nossos dados com a mesma matriz de distância que criamos acima, e usando a mesma variável “Ambiente” de sempre:

resultado<-adonis(matriz.bray~Ambiente)
resultado
## 
## Call:
## adonis(formula = matriz.bray ~ Ambiente) 
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Ambiente   1    0.8660 0.86599  10.454 0.27186  0.001 ***
## Residuals 28    2.3194 0.08284         0.72814           
## Total     29    3.1854                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 

A função envfit

A função envfit busca por relações lineares entre variáveis explanatórias da nossa escolha com os dados da ordenação gerada pelo NMDS. Pode ser bastante útil para se compreender como os padrões observados na ordenação se estruturam.

Sua aplicação é simples:

nmds.env<-envfit(nmds.bray, dados[,3:6])
nmds.env
## 
## ***VECTORS
## 
##                NMDS1    NMDS2     r2 Pr(>r)    
## Temperatura  0.74326  0.66901 0.6019  0.001 ***
## Cobertura   -0.82404 -0.56653 0.6149  0.001 ***
## Luz          0.80491  0.59340 0.5450  0.001 ***
## Flores       0.82216 -0.56926 0.8014  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## P values based on 999 permutations.

E podemos visualizar no gráfico:

plot(nmds.bray, type="n")
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available
points(nmds.bray, pch=16, col=Ambiente)
plot(nmds.env)

Se você se interessar por uma função mais complexa, que busca por relações não lineares, dê uma olhada na função ordisurf.

 

Pessoal, chegamos ao fim da nossa jornada. :D

Assim como ministrar o curso foi um imenso prazer, gostei bastante de gerar este material, e espero que seja útil para todos vocês.

Um grande abraços para todos!

Marcos