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")
# 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(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")
# 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")
#####
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")
# 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?