Curso de R no IPÊ
Script do dia 4.3: Mantel
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 quarto 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 parte do guia nós vamos realizar um teste de Mantel, um teste de Mantel parcial e calcular um correlograma.
Como sempre, todo o documento foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais, acesse:
Antes de começar, carregue os dados necessários:
Mentira, desta vez vamos usar dados que já estão disponíveis no R, que tal?
library(vegan)## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.0-10data(mite)
data(mite.env)
data(mite.xy)A função data carrega pacotes que fazem parte dos exemplos do R.
O teste de Mantel
Sua execução é bem simples no R, basta usar a função mantel, então vejamos um exemplo:
dist.amb<-dist(scale(mite.env[,1:2]), method="euclid")
dist.jac<-vegdist(mite, method="jac", binary=T)
mantel(dist.amb, dist.jac)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = dist.amb, ydis = dist.jac)
##
## Mantel statistic r: 0.4263
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0611 0.0801 0.0971 0.1128
##
## Based on 999 permutationsO gráfico a seguir pode ser útil para se pensar no resultado:
plot(dist.amb, dist.jac, pch=20)
O teste de Mantel Parcial
Também não há complicações aqui, basta usar a função mantel.partial:
dist.geo<-dist(mite.xy, method="euclid")
mantel.partial(dist.amb, dist.jac, dist.geo)##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = dist.amb, ydis = dist.jac, zdis = dist.geo)
##
## Mantel statistic r: 0.3323
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0863 0.1113 0.1344 0.1589
##
## Based on 999 permutations
E o correlograma de mantel
De novo, nada de complicado, basta usar a função mantel.correlog (mas você pode querer explorar mais a função para determinar as classes de distância manualmente, por exemplo).
Segue um exemplo rápido de execução:
correlograma<-mantel.correlog(dist.amb, dist.geo)
correlograma##
## Mantel Correlogram Analysis
##
## Call:
##
## mantel.correlog(D.eco = dist.amb, D.geo = dist.geo)
##
## class.index n.dist Mantel.cor Pr(Mantel) Pr(corrected)
## D.cl.1 0.5141820 358.0000000 0.0994010 0.001 0.001 ***
## D.cl.2 1.2425460 650.0000000 0.1152941 0.001 0.002 **
## D.cl.3 1.9709099 796.0000000 0.1050915 0.001 0.003 **
## D.cl.4 2.6992739 696.0000000 0.0419353 0.026 0.026 *
## D.cl.5 3.4276379 500.0000000 0.0053958 0.422 0.422
## D.cl.6 4.1560019 468.0000000 -0.0338560 0.054 0.108
## D.cl.7 4.8843659 364.0000000 -0.0538239 0.005 0.020 *
## D.cl.8 5.6127298 326.0000000 NA NA NA
## D.cl.9 6.3410938 260.0000000 NA NA NA
## D.cl.10 7.0694578 184.0000000 NA NA NA
## D.cl.11 7.7978218 130.0000000 NA NA NA
## D.cl.12 8.5261858 66.0000000 NA NA NA
## D.cl.13 9.2545497 32.0000000 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(correlograma)
Extra! Índices de diversidade de brinde!
Como em sala acabamos falando um pouco sobre os índices e os perfis de diversidade, seguem aqui os procedimentos que realizamos.
Calculando o índice de Shannon:
diversity(colSums(mite), index="shannon")## [1] 2.65044Calculando o índice de Simpson:
diversity(colSums(mite), index="simpson")## [1] 0.885881Deixando os índices pra la, e calculando um perfil de diversidade:
plot(renyi(colSums(mite), hill=T), type="b")
E é isso aí! Conseguimos avançar quatro dias.
Agora nos vemos nos documentos do último dia!