VeganOrdin.R

darren — Jun 5, 2014, 12:33 PM

# vegan ordination ############################################################
# baseado em http://cc.oulu.fi/~jarioksa/opetus/metodi/vegantutor.pdf
# vegan pode fazer qualquer tipo de ordenação, 3 exemplos
#### Unconstrained ordination
# 1 # Escalonamento Multidimensional Não Métrico 
# (non-metric multidimensional scaling,NMDS)
# 2 detrended correspondence analysis

#### Constrained ordination
# 3 constrained (canonical) correspondence analysis
###############################################################################
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.0-10
#### dados
# dados de uma comunidade de arbustos
data(dune)
mydune<-dune

# dados ambientais
data(dune.env)
mydune.env<-dune.env

# 1 # Escalonamento Multidimensional Não Métrico 
# (non-metric multidimensional scaling,NMDS)
# mapeamento das dessemelhanças da comunidade em uma forma não linear
# sobre o "espaço da ordenação" e pode lidar com respostas não-lineares das espécies
# de qualquer forma.

dune.mds <- metaMDS(mydune, trace = FALSE)
dune.mds

Call:
metaMDS(comm = mydune, trace = FALSE) 

global Multidimensional Scaling using monoMDS

Data:     mydune 
Distance: bray 

Dimensions: 2 
Stress:     0.1183 
Stress type 1, weak ties
Two convergent solutions found after 2 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'mydune' 
# stress, valor de 0 a 1
#1. Function used Bray-Curtis dissimilarities.
#2. Function run isoMDS with several random starts, and stopped either
#   after a certain number of tries, or after nding two similar
#   congurations with minimum stress. In any case, it returned the
#   best solution.
#3. Function rotated the solution so that the largest variance of site
#   scores will be on the first axis.
#4. Function scaled the solution so that one unit corresponds to halving
#   of community similarity from the replicate similarity.
#5. Function found species scores as weighted averages of site scores,
#   but expanded them so that species and site scores have equal variances.
#   This expansion can be undone using shrink = TRUE in display
#   commands.
plot(dune.mds, type = "t")

plot of chunk unnamed-chunk-1


# função "envfit", 
# vai encaixa/ajustar vetores ambientais ou fatores para uma ordenação
# resultado com:
# "r" = Goodness of fit (Testes de adequação (aderência))
# e "pvals" = valores de P para cada variável.
ef <- envfit(dune.mds, mydune.env, permu = 999)
ef

***VECTORS

   NMDS1 NMDS2   r2 Pr(>r)  
A1 0.965 0.263 0.36  0.022 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P values based on 999 permutations.

***FACTORS:

Centroids:
             NMDS1 NMDS2
Moisture1    -0.51 -0.04
Moisture2    -0.39  0.01
Moisture4     0.28 -0.40
Moisture5     0.66  0.15
ManagementBF -0.45 -0.01
ManagementHF -0.26 -0.13
ManagementNM  0.30  0.58
ManagementSF  0.15 -0.47
UseHayfield  -0.16  0.32
UseHaypastu  -0.04 -0.34
UsePasture    0.29  0.08
Manure0       0.30  0.58
Manure1      -0.25 -0.02
Manure2      -0.31 -0.19
Manure3       0.31 -0.25
Manure4      -0.35 -0.56

Goodness of fit:
             r2 Pr(>r)   
Moisture   0.50  0.002 **
Management 0.41  0.009 **
Use        0.19  0.110   
Manure     0.42  0.016 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P values based on 999 permutations.

# gráfico com os dados ambientais 
windows(width=8,height=8)
plot(dune.mds, display = "sites")
plot(ef, p.max = 0.05) # somente os mais significativas

plot of chunk unnamed-chunk-1

plot(ef$vectors, p.max = 0.1)
Error: 'x' is a list, but does not have components 'x' and 'y'
# mas é muito difícil entender

# vamos tentar algumas melhorias
# 3 etapas
# 1. fatores de interesse: manejo e uso
# Então tipo de manejo (Management) vai ter cores diferentes
# mas tem que definir isso antes de fazer o gráfico
#mydune.env$Management
# Levels: BF    HF    NM    SF
#       blue, green,yellow, black 
myc<-c("blue","black","green","black","green","yellow","blue","green","yellow","black","yellow","black","green")
myb<-c(1,3,1,1,2,2,2,1,1,1,3,1,1)
colf<-rep(myc,myb)
#2. espécies, pelo ordem de importância
priSpp <- diversity(mydune, index = "invsimpson", MARGIN = 2)
#3. dados ambientais, mas somente os dados contínuos
vf <- envfit(dune.mds, mydune.env[,'A1'], permu = 999)
dimnames(vf$vectors$arrows)[1]<-"A1"

# agora faço o gráfico e acrescentar as camadas
# trabalhar com camadas, acrescentando informações novos cada vez
windows(12,10)
par(mar=c(5,5,4,2)+0.1)
plot(dune.mds, type = "n",cex.lab=2.2,cex.axis=2.2,cex.main=2.6,main=paste("NMDS/Bray-Curtis - Stress =",round(dune.mds$stress,3)))
text(dune.mds, display = "sites", cex=1.2,labels=mydune.env$Use,col=colf)
ordilabel(dune.mds, display = "species", font = 2,
          priority = priSpp)
plot(vf,col = "red")

plot of chunk unnamed-chunk-1


# mas "A1" é importante?
# melhor subconjunto de variáveis ambientais
dfA1<-data.frame(A1=mydune.env$A1)
bioenv(mydune,dfA1)

Call:
bioenv(comm = mydune, env = dfA1) 

Subset of environmental variables with best correlation to community data.

Correlations:      spearman 
Dissimilarities:   bray 

Best model has 1 parameters (max. 1 allowed):
A1
with correlation  0.1967 

# significativa?
dune.dist <- vegdist(mydune) # Bray-Curtis
env.dist <- vegdist(dfA1, "euclid")
mantel(dune.dist, env.dist)

Mantel statistic based on Pearson's product-moment correlation 

Call:
mantel(xdis = dune.dist, ydis = env.dist) 

Mantel statistic r: 0.238 
      Significance: 0.048 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.183 0.233 0.278 0.314 

Based on 999 permutations
mantel(dune.dist, env.dist,method="spearman")

Mantel statistic based on Spearman's rank correlation rho 

Call:
mantel(xdis = dune.dist, ydis = env.dist, method = "spearman") 

Mantel statistic r: 0.198 
      Significance: 0.066 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.163 0.215 0.251 0.307 

Based on 999 permutations

# para avaliar a importância de variáveis categóricas (tipo de manejo etc)
# podemos adotar outras analises "classificações" ou "Analises de agrupamento"
# mas nos não temos tempo!

#### Até este ponto as informações são referentes a abundancia das espécies
#### agora vamos fazer com dados de presença e ausência
#### transformação de dados
mydune.pa<-decostand(mydune,method="pa")
# agora nmds com distancia (semelhança) calculado através da indicie "Jaccard"
# "Jaccard" é o mais comum para dados de presença/ausência
dune.mds.pa <- metaMDS(mydune.pa, distance="jaccard",trace = FALSE)
dune.mds.pa

Call:
metaMDS(comm = mydune.pa, distance = "jaccard", trace = FALSE) 

global Multidimensional Scaling using monoMDS

Data:     mydune.pa 
Distance: jaccard 

Dimensions: 2 
Stress:     0.1067 
Stress type 1, weak ties
Two convergent solutions found after 4 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'mydune.pa' 
plot(dune.mds.pa,type="t")

plot of chunk unnamed-chunk-1


#####
windows(width=8,height=5)
par(mfrow=c(1,2))
plot(dune.mds,type="t",main="nmds - abundancias")
plot(dune.mds.pa,type="t",main="nmds- presença/ausença")

plot of chunk unnamed-chunk-1


# Agora fazer os gráficos e testes
# e comparar a nmds de abundancia com nmds de presença ausência
# respondendo 2 perguntas:
# 1) nmds de ocorrência tem o mesmo padrão que abundancia?
# apresentar gráficos adequados!
# 2) A1 tem a mesma importância na padrão de ocorrência que
#     na padrão de abundancia?