Atividade Resumo Geral
Luiz Gustavo
07/10/2022
Pacotes e dados
# Carregando pacotes
library(ade4)
library(ecodados)
library(tidyverse)
library(vegan)
library(pvclust)
library(BiodiversityR)
library(labdsv)
library(ggplot2)
library(gridExtra)
library(ape)
library(FactoMineR)
library(factoextra)
library(FD)
library(palmerpenguins)
library(GGally)
library(fields)
library(ade4)
library(ggord)
library(udunits2)
library(adespatial)
library(spdep)
library(mvabund)
library(reshape)
#Dados utilizados
sp_compos <- ecodados::bocaina
species <- ecodados::com_birds
env <- ecodados::env_birds
xy <- ecodados::birds.xy
bocaina.env <- ecodados::bocaina.env
bocaina.xy <- ecodados::bocaina.xy
anuros_permanova <- ecodados::anuros_permanova
macroinv <- ecodados::macroinv
fish_comm <- ecodados::fish_comm
data(mite)
data(doubs)
data(mite.env)
## Traduzir nomes para português
colnames(penguins) <- c("especies", "ilha", "comprimento_bico", "profundidade_bico", "comprimento_nadadeira", "massa_corporal", "sexo", "ano")Reprodução dos exemplos do capítulo
Análises de agrupamento
Agrupamento hierárquico
Agrupamento baseado em nível de corte 50%
Serão analisados dados de girinos de anuros de diferentes espécies em 14 poças diferentes, bucando encontrar um padrão de ocorrência das espécies nesses locais. Parte-se do pressuposto que serão encontrados dois grupos: poças em locais mais cobertos pelo dossel e poças em locais menos cobertos.
## Matriz de similaridade com o coeficiente de Morisita-Horn
distBocaina <- vegdist(x = sp_compos, method = "horn")
## Agrupamento com a função hclust e o método UPGMA
dendro <- hclust(d = distBocaina, method = "average")
## Visualização os resultados
plot(dendro, main = "Dendrograma",
ylab = "Índice de Horn de similaridade)",
xlab="", sub="")
Medição da qualidade do dendograma pelo coeficiente de correlação cofenética.
## Coeficiente de correlação cofenética
cofresult <- cophenetic(dendro)
cor(cofresult, distBocaina)## [1] 0.9455221Notamos então, que há uma boa representação, que não distorce a informação.
Para o gráfico abaixo, foi utilizada uma linha de corte de 50%, ou seja, os grupos formados tem pelo menos 50% de nível de similaridade.
## Gráfico
plot(dendro, main = "Dendrograma",
ylab = "Índice de Horn de similaridade",
xlab="", sub="")
k <- 4
n <- ncol(sp_compos)
MidPoint <- (dendro$height[n-k] + dendro$height[n-k+1]) / 2
abline(h = MidPoint, lty=2)
Vemos a formação de 5 grupos abaixo da linha de corte. Logo, a hipótese sugerida de que havia 2 grupos foi refutada.
Comparação entre diferentes coeficiente de associação - Chord e Hellinger
Para ambos se utilizou uma linha de corte criada por Bootstrap.
Agrupamento baseado na distância de Chord
Será utilizada a transformação Chord e calculada a distância euclidiana; o que equivale a calcular a distância de Chord.
## Transformação dos dados com Chord
bocaina_transf <- disttransform(t(sp_compos), "chord")
## Cálculo da distância euclidiana
analise <- pvclust(bocaina_transf, method.hclust = "average", method.dist = "euclidean", quiet = TRUE)
## Plot do dendrograma
plot(analise, hang=-1, main = "Dendrograma com valores de P",
ylab = "Distância Euclideana",
xlab="", sub="")
pvrect(analise)
Foi separado apenas um grupo.
Agrupamento baseado na distância de Hillinger
Será utilizada a distância Hillinger para transformar os dados.
## Transformação dos dados com Hellinger
bocaina_transf2 <- disttransform(t(bocaina), "hellinger")
## Cálculo da distância euclidiana
analise2 <- pvclust(bocaina_transf2, method.hclust="average", method.dist="euclidean", quiet = TRUE)
## Plot do dendrograma
plot(analise2, hang=-1, main = "Dendrograma com valores de P",
ylab = "Distância Euclideana",
xlab="", sub="")
k <- 4
n <- ncol(sp_compos)
MidPoint <- (dendro$height[n-k] + dendro$height[n-k+1]) / 2
abline(h = MidPoint, lty=2)
pvrect(analise2)
Apareceu mais um grupo se comparado ao Chord. Isso porque o Hellinger dá menos peso a espécies raras que o anterior.
Independendo do agrupamento utilizado, vê-se que a hipótese inicial ainda está refutada, pois não é encontrada a divisão em dois grupos.
Agrupamento não hierárquico - K-means
Serão analisados dados de 27 espécies de peixes coletados em 30 pontos ao longo do Rio Doubs, entre a França e Suíça. Buscamos um padrão de ocorrência de espécies de peixes ao longo desse rio.
Primeiramente, os dados precisam ser padronizados. Ao se analisar os dados foi notado que no local 8 não há ocorrência de nenhuma espécie, logo ele será removido das análises.
## Retirando o local 8
spe <- doubs$fish[-8,]
## Normalizando os dados
spe.norm <- decostand(x = spe, method = "normalize") Vamos formar 4 grupos com os dados.
## Argumento centers do kmeans diz a quantidade de grupos que será formada
spe.kmeans <- kmeans(x = spe.norm, centers = 4, nstart = 100)
spe.kmeans## K-means clustering with 4 clusters of sizes 12, 3, 6, 8
##
## Cluster means:
## Cogo Satr Phph Neba Thth Teso
## 1 0.10380209 0.542300691 0.50086515 0.43325916 0.114024105 0.075651573
## 2 0.00000000 0.000000000 0.00000000 0.00000000 0.000000000 0.000000000
## 3 0.06167791 0.122088022 0.26993915 0.35942538 0.032664966 0.135403325
## 4 0.00000000 0.006691097 0.02506109 0.06987391 0.006691097 0.006691097
## Chna Chto Lele Lece Baba Spbi Gogo
## 1 0.00000000 0.00000000 0.06983991 0.1237394 0.02385019 0.0000000 0.05670453
## 2 0.05205792 0.00000000 0.07647191 0.3166705 0.00000000 0.0000000 0.20500174
## 3 0.06212775 0.21568957 0.25887226 0.2722562 0.15647062 0.1574388 0.16822286
## 4 0.10687104 0.09377516 0.14194394 0.2011411 0.24327992 0.1326062 0.28386032
## Eslu Pefl Rham Legi Scer Cyca Titi
## 1 0.04722294 0.02949244 0.0000000 0.00000000 0.00000000 0.00000000 0.03833408
## 2 0.07647191 0.00000000 0.0000000 0.05205792 0.07647191 0.00000000 0.00000000
## 3 0.12276089 0.17261621 0.0793181 0.06190283 0.04516042 0.06190283 0.14539027
## 4 0.20630360 0.16920496 0.2214275 0.19066542 0.13171275 0.16019126 0.26230024
## Abbr Icme Acce Ruru Blbj Alal Anan
## 1 0.00000000 0.0000000 0.00000000 0.01049901 0.00000000 0.0000000 0.00000000
## 2 0.00000000 0.0000000 0.18058775 0.31667052 0.05205792 0.7618709 0.00000000
## 3 0.01473139 0.0000000 0.03192175 0.32201597 0.01473139 0.1095241 0.04739636
## 4 0.19561641 0.1331835 0.26713081 0.32103755 0.22883055 0.3326939 0.18873077
##
## Clustering vector:
## 1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
## 1 1 1 1 3 1 1 3 1 1 1 1 1 1 3 3 3 3 4 4 4 2 2 2 4 4
## 28 29 30
## 4 4 4
##
## Within cluster sum of squares by cluster:
## [1] 2.5101386 0.3560423 1.7361453 0.4696535
## (between_SS / total_SS = 66.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Para descobrir o número ideal de grupos, é realizada uma repetição do k-means com uma série de valores de K pela função cascadeKM()
## Repetindo o K-Means
spe.KM.cascade <- cascadeKM(spe.norm, inf.gr = 2, sup.gr = 10, iter = 100, criterion = "ssi")
## Mostrar os resultados
spe.KM.cascade$results## 2 groups 3 groups 4 groups 5 groups 6 groups 7 groups 8 groups
## SSE 8.2149405 6.4768108 5.0719796 4.3015573 3.58561200 2.952367 2.4840549
## ssi 0.1312111 0.1685126 0.1252312 0.1286618 0.05917495 0.157523 0.1519788
## 9 groups 10 groups
## SSE 2.05218880 1.7599292
## ssi 0.09260697 0.1222963## Plot do gráfico
plot(spe.KM.cascade, sortg = TRUE)
Vemos então, que o número ideal de grupos é 3.
Espécies indicadoras
O conjunto de dados dos girinos foi utilizado novamente em busca de relações entre as espécies e condições ambientais.
## Índice para identificar as poças por fitofisionomia
fitofis <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4), rep(4,5))
## Análise de espécies indicadoras
res_indval <- indval(t(sp_compos), fitofis)
## Tabela completa com o resultado das espécies
tab_indval <- cbind.data.frame(maxcls = res_indval$maxcls,
ind.value = res_indval$indcls,
P = res_indval$pval)
tab_indval## maxcls ind.value P
## Aper 3 0.2432796 0.918
## Bahe 2 0.6487329 0.069
## Rict 1 0.8363823 0.015
## Cleuco 3 0.4571210 0.236
## Dmic 3 0.7011284 0.118
## Dmin 1 0.7032145 0.097
## Hpoly 2 0.6208711 0.296
## Lfur 3 0.2279412 0.920
## Pbar 1 0.2813725 0.608
## Polf 3 0.2437500 0.931
## Pmelano 2 0.2500000 0.933
## Sduar 1 0.7474527 0.036
## Shay 3 0.4930269 0.387
## Sobt 2 0.2222222 0.717
## Ssqual 3 0.2500000 0.930
## Bokerm 2 0.4583333 0.235# Apenas as espécies com p<0,05 que sugere o status de indicadora
tab_indval[tab_indval$P < 0.05, ]## maxcls ind.value P
## Rict 1 0.8363823 0.015
## Sduar 1 0.7474527 0.036Portanto, vemos que 2 espécies foram consideradas indicadoras, ambas em relação a fitofisionomia 1: Rhinella icterica (Rict) e Scinax duartei (Sduar)
Análises de ordenação
Ordenação irrestrita
Análise de componentes principais (PCA)
Aranhas nas cidades
Primeiramente, devemos obter a matriz centralizada
## Dados
aranhas <- data.frame(
sp1 = c(5, 7, 2, 0, 0, 0, 0, 0),
sp2 = c(0, 6, 3, 4, 0, 0, 0, 0),
sp3 = c(0, 0, 0, 9, 12, 3, 0, 0),
sp4 = c(0, 0, 0, 0, 4, 10, 8, 0),
sp5 = c(0, 0, 0, 0, 0, 6, 9, 12),
row.names = paste0("cidade", 1:8))
## Matriz centralizada
aranha.cent <- as.data.frame(base::scale(aranhas, center = TRUE, scale=FALSE))Então, calculamos a matiz de covariância e, a partir dela, os autovetores e autovalores.
## Matriz de covariância
matriz_cov <- cov(aranha.cent)
## Autovalores e autovetores
eigen_aranhas <- eigen(matriz_cov)
autovalores <- eigen_aranhas$values
autovetores <- as.data.frame(eigen_aranhas$vectors)
colnames(autovetores) <- paste("PC", 1:5, sep="")
rownames(autovetores) <- colnames(aranhas)
autovetores## PC1 PC2 PC3 PC4 PC5
## sp1 -0.2144766 0.38855265 0.29239380 -0.02330706 0.8467522
## sp2 -0.2442026 0.17463316 0.01756743 0.94587037 -0.1220204
## sp3 -0.3558368 -0.80222917 -0.27591770 0.10991178 0.3762942
## sp4 0.4159852 -0.41786654 0.78820962 0.17374202 0.0297183
## sp5 0.7711688 0.01860152 -0.46560957 0.25003826 0.3544591## Componentes principais
matriz_F <- as.data.frame(as.matrix(aranha.cent) %*% as.matrix(autovetores))
matriz_F## PC1 PC2 PC3 PC4 PC5
## cidade1 -2.979363 4.4720575 1.1533417 -3.2641923 0.5433206
## cidade2 -4.873532 6.2969618 1.8435339 2.3644158 1.5047024
## cidade3 -3.068541 3.8302991 0.3288626 -0.3566600 -2.3629973
## cidade4 -6.086322 -3.9922356 -2.7216169 1.6250305 -0.7918743
## cidade5 -4.513082 -8.7689219 -0.4668012 -1.1337476 0.9439633
## cidade6 5.812374 -3.9444494 3.9520584 0.4197281 -0.1376205
## cidade7 8.361421 -0.6462243 1.8065636 0.4926235 -0.2625625
## cidade8 7.347046 2.7525126 -5.8959421 -0.1471979 0.5630683## Porcentagem de explicação de cada eixo
100 * (autovalores/sum(autovalores))## [1] 47.201691 35.008126 12.135225 3.807112 1.847846Os dois primeiros PCs explicam 47,20% e 35,01% da variação, então serão usados para plotar um gráfico que ajudará a visualizar a similaridade entre as espécies presentes em cada cidade.
## Plot do gráfico
ggplot(matriz_F, aes(x = PC1, y = PC2, label = rownames(matriz_F))) +
geom_label() +
geom_hline(yintercept = 0, linetype=2) +
geom_vline(xintercept = 0, linetype=2) +
xlim(-10, 10) +
tema_livro()
Pinguins
Nesse exemplo, foi utilizada uma base de dados de pinguins do arquipélago Palmer (Península Antártica), que contém medidas morfométricas de três espécies. Foi proposto que a dieta irá interferir na morfologia dos animais.
Primeiramente, foi checada a existência de NAs nos dados, haviam 19, que foram deletados. Foram mantidos também apenas dados contínuos para a PCA.
## Verificar se existem NAs nos dados
sum(is.na(penguins))## [1] 19## Remover dados ausentes (NA)
penguins <- na.omit(penguins)
## Manter somentes dados contínuos
penguins_trait <- penguins[, 3:6]Atavés da função scale.unit os dados foram padronizados para terem média 0 e variância 1, algo essencial para a execussão bem sucedida da PCA.
pca.p <- PCA(X = penguins_trait, scale.unit = TRUE, graph = FALSE)
fviz_screeplot(pca.p, addlabels = TRUE, ylim = c(0, 70), main = "",
xlab = "Dimensões",
ylab = "Porcentagem de variância explicada") 
Através disso, vemos que as dimensões 1 e 2 juntas explicam 88,1% da variância dos dados. Essas PCs serão utilizadas para o gráfico abaixo.
fviz_pca_biplot(X = pca.p,
geom.ind = "point",
fill.ind = penguins$especies,
col.ind = "black",
alpha.ind = 0.7,
pointshape = 21,
pointsize = 4,
palette = c("darkorange", "darkorchid", "cyan4"),
col.var = "black",
invisible = "quali",
title = NULL) +
labs(x = "PC1 (68.63%)", y = "PC2 (19.45%)") +
xlim(c(-4, 5)) +
ylim(c(-3, 3)) +
tema_livro()
Através dele podemos interpretar que PC1 está influencia bastante positivamente na massa corporal, comprimento da nadadeira e comprimento do bico, enquanto tem influência negativa na profundidade do bico. Já PC2 quase não tem relação com o comprimento da nadadeira e a massa corporal, mas está relacionada positivamente com a profundidade do bico e o comprimento do bico.
Análises de Coordenadas Principais (PCoA)
Será utilizado um conjunto de dados da composição de ácados Oribatidae em 70 manchas de musco no intuito de encontrar uma relação entre as espécies de áraro e as diferentes topografias. Foi previsto que se encontraria pelo menos 2 grupos: que ocorrem dentro de poças dentro de florestas e que ocorrem em poças em áreas abertas.
Primeiramente, se padronizou os dados com Hellinger e foi calculado uma matriz de distância Bray-Curtis, e então calculado a PCoA.
## Hellinger para padronizar
mite.hel <- decostand(x = mite, method = "hellinger")
## Matriz de distância de Bray-Curtis
sps.dis <- vegdist(x = mite.hel, method = "bray")
## PCoA
pcoa.sps <- pcoa(D = sps.dis, correction = "cailliez")Agora, vamos checar o quanto as 2 primeiras dimensões explicam os eixos.
## Eixo 1
100 * (pcoa.sps$values[, 1]/pcoa.sps$trace)[1]## [1] 49.10564## Eixo 2
100 * (pcoa.sps$values[, 1]/pcoa.sps$trace)[2]## [1] 14.30308## Acumulado dos dois primeiros eixos
sum(100 * (pcoa.sps$values[, 1]/pcoa.sps$trace)[1:2])## [1] 63.40872É possível ver que eles explicam cerca de 63,41% da variação. Agora vamos plotar o gráfico para auxiliar na visualização.
## Selecionar os dois primeiros eixos
eixos <- pcoa.sps$vectors[, 1:2]
## Junção com o dado topográfico
pcoa.dat <- data.frame(topografia = mite.env$Topo, eixos)
## Plot do gráfico
ggplot(pcoa.dat, aes(x = Axis.1, y = Axis.2, fill = topografia,
color = topografia, shape = topografia)) +
geom_point(size = 4, alpha = 0.7) +
scale_shape_manual(values = c(21, 22)) +
scale_color_manual(values = c("black", "black")) +
scale_fill_manual(values = c("darkorange", "cyan4")) +
labs(x = "PCO 1 (49.11%)", y = "PCO 2 (14.30%)") +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 0, linetype = 2) +
tema_livro()
Nota-se que o ácaro Hummock está mais presente onde as dimensões 1 e 2 são positivas, enquanto que o Blanket está onde essas dimensões são negativas, principalmente a 1.
Regressão de Componentes Principais (PCR)
Aves dos alpes franceses
Será utilizado um conjunto de dados referentes a aves de 23 regiões dos alpes franceses, utilizando dados ambientais de variáveis climáticas e altitude, buscando encontrar variação na riqueza de aves relacionada ao gradiente climática. Espera-se que elas sejam afetadas pela umidade e temperatura.
## Coluna 8 removida por ser variável qualitativa
env_cont <- env[,-8]
## PCA
env.pca <- PCA(env_cont, scale.unit = TRUE, graph = FALSE)
fviz_screeplot(env.pca, addlabels = TRUE, ylim = c(0, 70), main = "",
xlab = "Dimensões",
ylab = "Porcentagem de variância explicada") 
ind_env <- get_pca_ind(env.pca)
env.pca$eig ## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 4.26652359 60.9503370 60.95034
## comp 2 2.05340251 29.3343216 90.28466
## comp 3 0.29745014 4.2492878 94.53395
## comp 4 0.19896067 2.8422953 97.37624
## comp 5 0.11448717 1.6355310 99.01177
## comp 6 0.04312874 0.6161248 99.62790
## comp 7 0.02604718 0.3721025 100.00000Vemos que as duas primeiras dimensões explicam 90,3% da variação dos dados. Enquanto que as 3 primeiras explicam 94,54%.
O PCR é utilizado para reduzir a dimensionalidade. Vamos agora obter os valores de escores representativos de valores convertidos para serem usados um análises.
## Dados dos primeiros 3 eixos
pred.env <- ind_env$coord[, 1:3]
## Riqueza de espécies
riqueza <- specnumber(species)
## Junção dos valores
dat <- data.frame(pred.env, riqueza) Com esses dados, será realizada uma regressão multipla.
## Regressão múltipla
mod1 <- lm(riqueza ~ Dim.1 + Dim.2 + Dim.3, data = dat)
par(mfrow = c(2, 2))
plot(mod1) # verificar pressupostos dos modelos lineares
summary(mod1) # resultados do teste##
## Call:
## lm(formula = riqueza ~ Dim.1 + Dim.2 + Dim.3, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4008 -1.1729 0.4356 1.2072 2.4571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.30435 0.37639 35.347 < 2e-16 ***
## Dim.1 0.68591 0.18222 3.764 0.00131 **
## Dim.2 -0.09961 0.26267 -0.379 0.70874
## Dim.3 -0.21708 0.69014 -0.315 0.75654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.805 on 19 degrees of freedom
## Multiple R-squared: 0.4313, Adjusted R-squared: 0.3415
## F-statistic: 4.804 on 3 and 19 DF, p-value: 0.01179dimdesc(env.pca)$Dim.1 ##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## maxi.jan 0.9279073 1.846790e-10
## maxi.jul 0.8877675 1.607390e-08
## mini.jul 0.7117620 1.396338e-04
## mini.jan 0.6830371 3.282701e-04
## rain.jan -0.6516187 7.559358e-04
## rain.tot -0.7284135 8.112903e-05
## rain.jul -0.8300858 9.588034e-07Com isso, vemos que a dimensão 1 foi a única com um impacto estatísticamente relevante na riqueza de espécies. Se nota também que há uma correlação muito forte entre as temperaturas máximas dos meses de janeiro e julho e as chuvas de julho com a riqueza. Segue um gráfico mostrando um modelo linear correlacionando esses fatores.
ggplot(dat, aes(x = Dim.1, y = riqueza)) +
geom_smooth(method = lm, fill = "#525252", color = "black") +
geom_point(size = 4, shape = 21, alpha = 0.7, color = "#1a1a1a", fill = "cyan4") +
labs(x = "Gradiente ambiental (PC1)", y = "Riqueza de aves") +
tema_livro()
Ácaros
Mais uma vez serão utilizados os dados referentes aos ácaros. Dessa vez buscando uma relação entre variáveis climáticas, vegetacionais e topográficam com a riqueza desses animais.
Primeiro, foi utilizada a distância de Gower para calcular a matriz de distância dos dados. E então a PCoA.
## Matriz de distância
env.dist <- gowdis(mite.env)
## PCoA
env.mite.pco <- pcoa(env.dist, correction = "cailliez")Medição do quanto cada eixo explica os dados.
## Eixo 1
100 * (env.mite.pco$values[, 1]/env.mite.pco$trace)[1]## [1] 61.49635## Eixo 2
100 * (env.mite.pco$values[, 1]/env.mite.pco$trace)[2]## [1] 32.15486Vemos que juntos eles explicam 93,65% dos dados. Agora exportaremos os dados para as análises.
## Seleção os dois primeiros eixos
pred.scores.mite <- env.mite.pco$vectors[, 1:2]
## Riqueza de espécies
mite.riqueza <- specnumber(mite)
## Junção dos dados
pred.vars <- data.frame(riqueza = mite.riqueza, pred.scores.mite)
### Regressão múltipla
mod.mite <- lm(riqueza ~ Axis.1 + Axis.2, data = pred.vars)
par(mfrow = c(2, 2))
plot(mod.mite)
summary(mod.mite)##
## Call:
## lm(formula = riqueza ~ Axis.1 + Axis.2, data = pred.vars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6874 -2.3960 -0.1378 2.5032 8.6873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.1143 0.4523 33.415 < 2e-16 ***
## Axis.1 -11.4303 2.0013 -5.711 2.8e-07 ***
## Axis.2 5.6832 2.7677 2.053 0.0439 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.784 on 67 degrees of freedom
## Multiple R-squared: 0.3548, Adjusted R-squared: 0.3355
## F-statistic: 18.42 on 2 and 67 DF, p-value: 4.225e-07Vemos que ambas as dimensões tem relevância na explicação dos dados, porém a dimensão 1 tem uma relação mais forte.
Agora faremos os gráficos dos modelos lineares que ajudam na visualização da relação das dimensões com a riqueza de espécies.
g_acari_axi1 <- ggplot(pred.vars, aes(x = Axis.1, y = riqueza)) +
geom_smooth(method = lm, fill = "#525252", color = "black") +
geom_point(size = 4, shape = 21, alpha = 0.7, color = "#1a1a1a", fill="cyan4") +
labs(x = "Gradiente ambiental (PC1)", y = "Riqueza de ácaros") +
tema_livro()
g_acari_axi2 <- ggplot(pred.vars, aes(x = Axis.2, y = riqueza)) +
geom_smooth(method = lm, fill = "#525252", color = "black") +
geom_point(size = 4, shape = 21, alpha = 0.7, color = "#1a1a1a", fill = "darkorange") +
labs(x = "Gradiente ambiental (PC2)", y = "Riqueza de ácaros") +
tema_livro()
## Função para combinar os dois plots em uma única janela
grid.arrange(g_acari_axi1, g_acari_axi2, nrow = 1)
Ordenação restrita
Análise de Redundância (RDA)
Será utilizada a base de dados de aves no intuito de encontrar uma relação entre o clima e altitude e a composição de espécies.
Primeiro, se transformará a matriz de espécies utilizando Hellinger e removerá a coluna 8 que é uma variável categórica.
## Transformação de hellinger
species.hel <- decostand(x = species, method = "hellinger")
## Deleção da variável categórica
env.contin <- env[, -8]Seleção das variáveis a serem utilizadas pelo método Forward Selection.
sel.vars <- forward.sel(species.hel, env.contin)## Testing variable 1
## Testing variable 2
## Testing variable 3
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.223000 (> 0.050000)sel.vars$variables## [1] "rain.jul" "maxi.jul"## As variáveis rain.jul e maxi.jul serão selecionadas
env.sel <- env[,sel.vars$variables]Padronização da matriz ambiental e matiz final com as variáveis selecionadas
## Padronização
env.pad <- decostand(x = env.sel, method = "standardize")
## Matriz final
env.pad.cat <- data.frame(env.pad, altitude = env$altitude)RDA.
## RDA
rda.bird <- rda(species.hel ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)Quais os eixos mais significativos para representar a relação entre as variáveis selecionadas e a composição de espécies?
## Análise da significância dos eixos
res.axis <- anova.cca(rda.bird, by = "axis")
res.axis## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = species.hel ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)
## Df Variance F Pr(>F)
## RDA1 1 0.045759 12.0225 0.001 ***
## RDA2 1 0.009992 2.6252 0.062 .
## RDA3 1 0.007518 1.9752 0.126
## RDA4 1 0.003582 0.9410 0.527
## Residual 18 0.068510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O eixo 1 é o mais significativo, mas também há influência do 2.
Quais variáveis mais contribuem para a variação da composição de espécies?
res.var <- anova.cca(rda.bird, by = "term")
res.var## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = species.hel ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)
## Df Variance F Pr(>F)
## rain.jul 1 0.036514 9.5936 0.001 ***
## maxi.jul 1 0.011264 2.9596 0.011 *
## altitude 2 0.019071 2.5053 0.011 *
## Residual 18 0.068510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1As três variáveis escolhidas são capazes de explicar a composição das espécies, em ordem de significância: chuvas de julho, altitude, temperatura de julho.
Vamos obter o R2 do modelo.
r_quadr <- RsquareAdj(rda.bird)
r_quadr## $r.squared
## [1] 0.4938685
##
## $adj.r.squared
## [1] 0.3813949Agoras realizaremos uma Ordenação multi-escala (MSO), que ajuda a entender os resultados de ordenação em relação à distância geográfica.
bird.rda <- mso(rda.bird, xy, grain = 1, permutations = 99)
msoplot(bird.rda)
## Error variance of regression model underestimated by -1.4 percentÉ possível notar que existem algumas autocorrelações em certas distâncias.
E um gráfico triplot que possui dados das ordenações.
ggord(rda.bird, ptslab = TRUE, size = 1, addsize = 3, repel = TRUE) +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 0, linetype = 2) +
tema_livro()
Através dele podemos ver qual a influencia de cada variável em relação a cada uma das espécies
Análise de Redundância parcial (RDAp)
Vimos que houveram autocorrelações no MSO realizado e para uma interpretação correta da RDA isso não pode ocorrer. A RDAp inclui a matriz de dados espaciais como um valor condicional dentro da RDA, podendo solucionar esse problema.
Para realizar a RDAp, primeiro é preciso gerar autovetores espaciais.
## Passo 1: Gerar um arquivo LIST W: list binária de vizinhança
mat_knn <- knearneigh(as.matrix(xy), k = 2, longlat = FALSE)
mat_nb <- knn2nb(mat_knn, sym = TRUE)
mat_listw <- nb2listw(mat_nb, style = "W")
mat_listw## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 23
## Number of nonzero links: 58
## Percentage nonzero weights: 10.96408
## Average number of links: 2.521739
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 23 529 23 18.84444 96.01111MEM_mat <- scores.listw(mat_listw, MEM.autocor = "positive")
candidates <- listw.candidates(xy, nb = c("gab", "mst", "dnear"),
weights = c("binary", "flin"))
## Passo 3: Selecionar a melhor matriz SWM e executar o MEM
W_sel_mat <- listw.select(species.hel, candidates, MEM.autocor = "positive",
p.adjust = TRUE, method = "FWD")## Procedure stopped (alpha criteria): pvalue for variable 5 is 0.110000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.069000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.061000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 4 is 0.189000 (> 0.050000)## Passo 4: Matriz dos preditores espaciais escolhidos (MEMs)
spatial.pred <- as.data.frame(W_sel_mat$best$MEM.select)
## É necessário atribuir os nomes das linhas
rownames(spatial.pred) <- rownames(xy) Agora é possível realizar o RDAp.
## Combinar variáveis ambientais e espaciais
pred.vars <- data.frame(env.pad.cat, spatial.pred)
## RDA parcial
rda.p <- rda(species.hel ~
rain.jul + maxi.jul + altitude + # Preditores ambientais
Condition(MEM1 + MEM2 + MEM4 + MEM5), # Preditores espaciais
data = pred.vars)Agora vamos realizar os mesmos processos e nos fazer as mesmas perguntas que fizemos no RDA.
Quais os eixos mais significativos para representar a relação entre as variáveis selecionadas e a composição de espécies?
res.p.axis <- anova.cca(rda.p, by = "axis")
res.p.axis## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = species.hel ~ rain.jul + maxi.jul + altitude + Condition(MEM1 + MEM2 + MEM4 + MEM5), data = pred.vars)
## Df Variance F Pr(>F)
## RDA1 1 0.008471 2.1376 0.289
## RDA2 1 0.004830 1.2189 0.787
## RDA3 1 0.003240 0.8176 0.889
## RDA4 1 0.001891 0.4773 0.899
## Residual 14 0.055477Quais variáveis mais contribuem para a variação da composição de espécies?
res.p.var <- anova.cca(rda.p, by = "term")
res.p.var## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = species.hel ~ rain.jul + maxi.jul + altitude + Condition(MEM1 + MEM2 + MEM4 + MEM5), data = pred.vars)
## Df Variance F Pr(>F)
## rain.jul 1 0.004406 1.1119 0.340
## maxi.jul 1 0.004446 1.1220 0.339
## altitude 2 0.009579 1.2087 0.237
## Residual 14 0.055477Vamos obter o R2 do modelo.
RsquareAdj(rda.p)## $r.squared
## [1] 0.1361661
##
## $adj.r.squared
## [1] 0.02330319A partir desses dados chegamos a conclusão de que não há efeito estatísticamente relevante das variáveis ambientais selecionadas com a composição de espécies.
Padrão espacial da composição de espécies e variáveis ambientais
## Composição de espécies
pca.comp <- dudi.pca(species.hel, scale = FALSE, scannf = FALSE)
moran.comp <- moran.mc(pca.comp$li[, 1], mat_listw, 999)
moran.comp##
## Monte-Carlo simulation of Moran I
##
## data: pca.comp$li[, 1]
## weights: mat_listw
## number of simulations + 1: 1000
##
## statistic = 0.62815, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater## Variáveis ambientais
env$altitude <- as.factor(env$altitude)
ca.env <- dudi.hillsmith(env, scannf = FALSE)
moran.env <- moran.mc(ca.env$li[, 1], mat_listw, 999)
moran.env##
## Monte-Carlo simulation of Moran I
##
## data: ca.env$li[, 1]
## weights: mat_listw
## number of simulations + 1: 1000
##
## statistic = 0.72714, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greaterSeguindo esses resultados, é possível que a maior influência na composição de espécies seja a distribuição espacial. Vamos explorar essa ideia com o método de partição de variância.
pv.birds <- varpart(species.hel, env.pad.cat, spatial.pred)
plot(pv.birds, cutoff = -Inf)
mtext("Diagrama de Venn: partição de variância")
Com isso, é a autocorrelação espacial das variáveis ambientais que influenciam na composição das espécies (0.36), sem que haja efeito direto delas (0.02 e -0.03). ### Análise de Redundância baseada em distância (db-RDA)
Esse tipo de análise é útil quando se utiliza outras medidas de distância que não a euclidiana. Ela será utilizada para verificar se a diversidade beta vai influenciar na dissimilaridade das espécies de aves.
## Hellinger
species.hel <- decostand(species, "hellinger")
## Diversidade beta - Bray-Curtis
dbeta <- vegdist(species.hel, "bray")
## dbRDA
dbrda.bird <- capscale (dbeta~rain.jul+maxi.jul+altitude, data=env.pad.cat)Então realizaremos as mesmas perguntas.
Quais os eixos mais significativos para representar a relação entre as variáveis selecionadas e a composição de espécies?
db.res.axis <- anova.cca(dbrda.bird, by = "axis")
db.res.axis## Permutation test for capscale under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = dbeta ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)
## Df SumOfSqs F Pr(>F)
## CAP1 1 0.274186 16.9355 0.001 ***
## CAP2 1 0.044965 2.7773 0.107
## CAP3 1 0.018141 1.1205 0.620
## CAP4 1 0.009500 0.5868 0.794
## Residual 18 0.291420
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vemos que apenas o eixo 1 é explicativo.
Quais variáveis mais contribuem para a variação da composição de espécies?
db.res.var <- anova.cca(dbrda.bird, by = "term") ## Qual variável?
db.res.var## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = dbeta ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)
## Df SumOfSqs F Pr(>F)
## rain.jul 1 0.209981 12.9697 0.001 ***
## maxi.jul 1 0.041846 2.5847 0.058 .
## altitude 2 0.094966 2.9328 0.011 *
## Residual 18 0.291420
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Todas as variáveis ambientais contribuem para a composição das espécies, porém as chuvas de julho contribuem mais seguidas pela altitudo e pela temperatura de julho.
Vamos obter o R2 do modelo.
db_r_quadr <- RsquareAdj(dbrda.bird)
db_r_quadr## $r.squared
## [1] 0.6116758
##
## $adj.r.squared
## [1] 0.5253816Esse R2 sugere uma relação forte.
PERMANOVA
Utilizada no teste de hipóteses multivariadas, comparando abundância de diferentes espécies submetidos a diferentes tratamentos ou gradientes ambientais.
Vamos usá-la para medir se diferenças climáticas e de altitude alteram a composição de espécies de aves.
## Hellinger
species.hel <- decostand(x = species, method = "hellinger")
## Bray-Curtis
sps.dis <- vegdist(x = species.hel, method = "bray")Escolha das variáveis
## Verifica correlação entre as variáveis
ggpairs(env) +
tema_livro()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
A partir disso, vamos manter somente três variáveis: mini.jan, rain.tot e altitude
env2 <- env[, c("mini.jan", "rain.tot", "altitude")]Agora vamos para a execussão da PERMANOVA em si.
perm.aves <- adonis2(sps.dis ~ mini.jan + rain.tot + altitude, data = env2)
perm.aves## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = sps.dis ~ mini.jan + rain.tot + altitude, data = env2)
## Df SumOfSqs R2 F Pr(>F)
## mini.jan 1 0.09069 0.15997 6.3307 0.002 **
## rain.tot 1 0.12910 0.22771 9.0118 0.002 **
## altitude 2 0.08929 0.15749 3.1163 0.019 *
## Residual 18 0.25787 0.45483
## Total 22 0.56695 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Será realizado uma BETADISPER, que permite saber se as diferenças entre os grupos ocorrem principalmente por diferenças na dispersão ou pela posição.
betad.aves <- betadisper(sps.dis, env2$altitude)
permutest(betad.aves)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.0042643 0.0021322 1.4672 999 0.269
## Residuals 20 0.0290636 0.0014532Como Pr = 0,218, a homogenidade de variância entre os grupos é confirmada. O que vai sugerir que a variação da composição das espécies está de acordo com os resultados da PERMANOVA e se devem pelas variáveis ambientais utilizadas.
Combinação de análises de ordenação irrestritas (PCA, PCoA, nMDS) e análises que testam hipóteses (PERMANOVA e RDA)
Implementação do nMDS na matriz de composição de espécies de ácaro.
as.matrix(sps.dis)[1:6, 1:6]## S01 S02 S03 S04 S05 S06
## S01 0.0000000 0.15133734 0.16720405 0.2559122 0.2559882 0.2588892
## S02 0.1513373 0.00000000 0.04114702 0.1190172 0.1289682 0.1391056
## S03 0.1672040 0.04114702 0.00000000 0.1420233 0.1358127 0.1410867
## S04 0.2559122 0.11901720 0.14202325 0.0000000 0.1140283 0.1199803
## S05 0.2559882 0.12896823 0.13581271 0.1140283 0.0000000 0.2129054
## S06 0.2588892 0.13910558 0.14108668 0.1199803 0.2129054 0.0000000## Cálculo da melhor solução do nMDS
sol1 <- metaMDS(sps.dis)## Run 0 stress 0.1344042
## Run 1 stress 0.1336481
## ... New best solution
## ... Procrustes: rmse 0.08548225 max resid 0.2653807
## Run 2 stress 0.1344042
## Run 3 stress 0.1272418
## ... New best solution
## ... Procrustes: rmse 0.1180675 max resid 0.3337563
## Run 4 stress 0.1344042
## Run 5 stress 0.1336481
## Run 6 stress 0.1336481
## Run 7 stress 0.1272417
## ... New best solution
## ... Procrustes: rmse 4.980776e-05 max resid 0.0001905159
## ... Similar to previous best
## Run 8 stress 0.1336481
## Run 9 stress 0.1378506
## Run 10 stress 0.1272418
## ... Procrustes: rmse 2.735292e-05 max resid 0.0001041061
## ... Similar to previous best
## Run 11 stress 0.1344042
## Run 12 stress 0.1378506
## Run 13 stress 0.1378506
## Run 14 stress 0.1344042
## Run 15 stress 0.1272417
## ... New best solution
## ... Procrustes: rmse 8.100243e-06 max resid 2.389456e-05
## ... Similar to previous best
## Run 16 stress 0.1272418
## ... Procrustes: rmse 2.184169e-05 max resid 8.074535e-05
## ... Similar to previous best
## Run 17 stress 0.1338432
## Run 18 stress 0.1378506
## Run 19 stress 0.1336481
## Run 20 stress 0.1338432
## *** Solution reachednmds.beta <- metaMDS(sps.dis, previous.best = sol1)## Starting from 2-dimensional configuration
## Run 0 stress 0.1272417
## Run 1 stress 0.1587247
## Run 2 stress 0.1272418
## ... Procrustes: rmse 3.518089e-05 max resid 0.0001322552
## ... Similar to previous best
## Run 3 stress 0.1272418
## ... Procrustes: rmse 1.111972e-05 max resid 3.193033e-05
## ... Similar to previous best
## Run 4 stress 0.1344042
## Run 5 stress 0.1272417
## ... New best solution
## ... Procrustes: rmse 9.710327e-06 max resid 3.647316e-05
## ... Similar to previous best
## Run 6 stress 0.1338432
## Run 7 stress 0.1272418
## ... Procrustes: rmse 3.461784e-05 max resid 0.0001330388
## ... Similar to previous best
## Run 8 stress 0.1272417
## ... Procrustes: rmse 1.350778e-05 max resid 5.218424e-05
## ... Similar to previous best
## Run 9 stress 0.1344042
## Run 10 stress 0.1344042
## Run 11 stress 0.1344042
## Run 12 stress 0.1272417
## ... Procrustes: rmse 6.51214e-06 max resid 2.485084e-05
## ... Similar to previous best
## Run 13 stress 0.1618991
## Run 14 stress 0.1338432
## Run 15 stress 0.1338432
## Run 16 stress 0.1344043
## Run 17 stress 0.1272418
## ... Procrustes: rmse 8.730179e-05 max resid 0.0003342361
## ... Similar to previous best
## Run 18 stress 0.1338434
## Run 19 stress 0.1272418
## ... Procrustes: rmse 6.420549e-05 max resid 0.0002461334
## ... Similar to previous best
## Run 20 stress 0.1580523
## *** Solution reachedChecar se o valor do stress sugere qualidade de ordenação.
nmds.beta$stress## [1] 0.1272417Como o número está entre 0 e 0,2, sim, a ordenação é de qualidade.
Montagem do gráfico
## Montar o dataframe
dat.graf <- data.frame(nmds.beta$points, altitude = env2$altitude)
## Definir os grupos ("HULL") para serem categorizados no gráfico
grp.mon <- dat.graf[dat.graf$altitude == "Montanhoso", ][chull(dat.graf[dat.graf$altitude == "Montanhoso", c("MDS1", "MDS2")]), ]
grp.int <- dat.graf[dat.graf$altitude == "Intermediário", ][chull(dat.graf[dat.graf$altitude == "Intermediário", c("MDS1", "MDS2")]), ]
grp.pla <- dat.graf[dat.graf$altitude == "Plano", ][chull(dat.graf[dat.graf$altitude == "Plano", c("MDS1", "MDS2")]), ]
## Combinar dados dos grupos para cada Convex Hull
hull.data <- rbind(grp.mon, grp.int, grp.pla)
## Plot do gráfico
ggplot(dat.graf, aes(x = MDS1, y = MDS2, color = altitude, shape = altitude)) +
geom_point(size = 4, alpha = 0.7) +
geom_polygon(data = hull.data, aes(fill = altitude, group = altitude), alpha = 0.3) +
scale_color_manual(values = c("darkorange", "darkorchid", "cyan4")) +
scale_fill_manual(values = c("darkorange", "darkorchid", "cyan4")) +
labs(x = "NMDS1", y = "NMDS2") +
tema_livro()
Mantel e Mantel parcial
Se aplica quando há duas ou mais matrizes de distância, testando a correlação entre elas para o efeito de uma terceira.
Vamos utilizar o conjunto de dados dos anuros para testar se há uma correlação entre a dissimilaridade na composição de espécies e dissimilaridade ambiental.
Primiero, montaremos as matrizes de dissimilaridade.
## Matriz de distância geográfica
Dist.km <- as.dist(rdist.earth(bocaina.xy, miles=F))
## Padronizações
comp.bo.pad <- decostand(t(bocaina), "hellinger")
env.bocaina <- decostand(bocaina.env[,-9], "standardize")
## Dissimilaridades
dissimil.bocaina <- vegdist(comp.bo.pad, "bray")
dissimil.env <- vegdist(env.bocaina, "euclidian")Aplicando Mantel para verificar a correlação entre dissimilaridade ambiental e distância geográfica.
mantel(Dist.km, dissimil.env) ##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = Dist.km, ydis = dissimil.env)
##
## Mantel statistic r: 0.4848
## Significance: 0.002
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.261 0.351 0.404 0.433
## Permutation: free
## Number of permutations: 999Agora plotaremos o gráfico.
## Junção dos dados
matrix.dist.env <- data.frame(x = melt(as.matrix(Dist.km))$value,
y = melt(as.matrix(dissimil.env))$value)
## Plot do gáfico
ggplot(matrix.dist.env , aes(x, y)) +
geom_point(size = 4, shape = 21, fill = "darkorange", alpha = 0.7) +
labs(x = "Distância geográfica (km)",
y = "Dissimilaridade Ambiental") +
tema_livro() 
Aplicando Mantel para verificar a correlação entre dissimilaridade na composição de espécies e distância geográfica.
mantel(Dist.km, dissimil.bocaina) ##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = Dist.km, ydis = dissimil.bocaina)
##
## Mantel statistic r: 0.3745
## Significance: 0.007
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.174 0.253 0.296 0.347
## Permutation: free
## Number of permutations: 999Agora plotaremos o gráfico.
## Junção dos dados
matrix.dist.bocaina <- data.frame(x = melt(as.matrix(Dist.km))$value,
y = melt(as.matrix(dissimil.bocaina))$value)
## Plot do gráfico
ggplot(matrix.dist.bocaina , aes(x, y)) +
geom_point(size = 4, shape = 21, fill = "darkorange", alpha = 0.7) +
labs(x = "Distância geográfica (km)",
y = "Dissimilaridade (Bray-Curtis)") +
tema_livro()
Agora, será medida a correlação entre as duas dissimilaridades: ambiental e na composição de espécies.
mantel(dissimil.env, dissimil.bocaina)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = dissimil.env, ydis = dissimil.bocaina)
##
## Mantel statistic r: 0.2898
## Significance: 0.009
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.158 0.207 0.244 0.276
## Permutation: free
## Number of permutations: 999## Junção dos dados
matrix.bocaina.env <- data.frame(x = melt(as.matrix(dissimil.bocaina))$value,
y = melt(as.matrix(dissimil.env))$value)
## Plot do gráfico
ggplot(matrix.bocaina.env, aes(x, y)) +
geom_point(size = 4, shape = 21, fill = "darkorange", alpha = 0.7) +
labs(x = "Similaridade Bocaina",
y = "Similaridade Ambiental") +
tema_livro()
Isso mostra que existe uma relação positiva fraca entre as dissimilaridades. E também foi notada uma relação entre as dissimilaridades e a localização espacial.
Para resultados mais condizentes com a realidade seria bom levar em consideração também o espaço e para isso se utiliza o teste de Mantel Parcial.
mantel.partial(dissimil.bocaina, dissimil.env, Dist.km)##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = dissimil.bocaina, ydis = dissimil.env, zdis = Dist.km)
##
## Mantel statistic r: 0.1334
## Significance: 0.13
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.154 0.196 0.239 0.285
## Permutation: free
## Number of permutations: 999Esse teste nos mostra que a relação ambiente e espécies existe, mas a força dela diminue bastante comparando com o teste anterior.
Mantel espacial com modelo nulo restrito considerando autocorrelação espacial
Utiliza-se Moran Spectral Randomization, um procedimento vindo de Mantel que lida melhor com autocorrelação espacial.
compos_espac <- mantel.randtest(sqrt(dissimil.bocaina), Dist.km)
compos_espac## Monte-Carlo test
## Call: mantel.randtest(m1 = sqrt(dissimil.bocaina), m2 = Dist.km)
##
## Observation: 0.3738076
##
## Based on 999 replicates
## Simulated p-value: 0.014
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 3.2392457797 0.0003927071 0.0132890962Construção de uma rede de vizinhança, utilizando da Minimum Spanning Tree.
nb.boc <- mst.nb(Dist.km)
plot(nb.boc, bocaina.xy)
Conversão para listw, para que possa ser utilizado na análise.
## Conversão
lw <- nb2listw(nb.boc)
## Mantel espacial
msr(compos_espac, lw, method = "pair")## Monte-Carlo test
## Call: msr.mantelrtest(x = compos_espac, listwORorthobasis = lw, method = "pair")
##
## Observation: 0.3738076
##
## Based on 999 replicates
## Simulated p-value: 0.104
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 1.3432690 0.1651373 0.0241321Esse procedimento fez com que os dados de dissimilaridade se provassem não significativos, em contraste com o do teste anterior.
PROCRUSTES e PROTEST
PROCRUSTES compara o grau de associação entre conjuntos multidimensionais de dados, que utiliza do teste PROTEST em seu processo.
Vamos utilizar conjuntos de dados de peixes e macroinvertebrados aquáticos de um riacho para verificar se existe concordância entre distribuição na composição de espécies desses animais.
## Como é um teste hipotético, fixar a amostragem
set.seed(1001)
## Matrizes de distância de PCoA
d_macro <- vegdist(macroinv, "bray")
pcoa_macro <- cmdscale(d_macro)
d_fish <- vegdist(fish_comm, "bray")
pcoa_fish <- cmdscale(d_fish)
## PROCRUSTES
concord <- procrustes(pcoa_macro, pcoa_fish)
protest(pcoa_macro, pcoa_fish)##
## Call:
## protest(X = pcoa_macro, Y = pcoa_fish)
##
## Procrustes Sum of Squares (m12 squared): 0.6185
## Correlation in a symmetric Procrustes rotation: 0.6176
## Significance: 0.019
##
## Permutation: free
## Number of permutations: 999## Gráfico
plot(concord, main = "",
xlab = "Dimensão 1",
ylab = "Dimensão 2")
Por conta do valor de m12 (0,6185) e de p (0,019), concluimos que existe concordância entre a composição de peixes e macroinvertebrados aquáticos. Ou seja, o nível de poluição cria padrões previsíveis para distribuição de ambos.
Métodos multivariados baseados em modelos
## Retirando a coluna que contém o fator e deixando apenas dados de abundância
anuros_abund <- as.data.frame(anuros_permanova[,-28])
## Incluindo apenas o fator a ser testado no modelo
grupos <- factor(anuros_permanova[,28])
## Média-variância
abund_tr <- mvabund(anuros_abund)
meanvar.plot(abund_tr)
## Gráfico
plot(abund_tr ~ grupos, cex.axis = 0.8, cex = 0.8)## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Leptodactylus_mystacinus, Leptodactylus_fuscus, Physalaemus_cuvieri, Physalaemus_nattereri, Rhinella_ornata, Scinax_fuscovarius, Hypsiboas_albopunctatus, Leptodactylus_mystaceus, Dendropsophus_minutus, Dendropsophus_jimi, Chiasmocleis_albopunctata, Itapotihyla_langsdorffii were included in the plot
## (the variables with highest total abundance).
## Modelo
modelo1 <- manyglm(abund_tr ~ grupos, family = "negative.binomial")
## Diagnose
plot(modelo1)
## Resultados
summary(modelo1)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 18.947 0.001 ***
## gruposCOLECAO 4.727 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 4.727, p-value: 0.004
## Arguments:
## Test statistics calculated assuming response assumed to be uncorrelated
## P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).## ANOVA
anova(modelo1, p.uni = "adjusted")## Time elapsed: 0 hr 0 min 5 sec## Analysis of Deviance Table
##
## Model: abund_tr ~ grupos
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 8
## grupos 7 1 94.08 0.01 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Chiasmocleis_albopunctata Dendropsophus_elianae
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 3.308 0.464 0.394 0.859
## Dendropsophus_jimi Dendropsophus_nanus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 0.23 0.859 2.197 0.721
## Dendropsophus_minutus Dendropsophus_sanborni
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 4.752 0.309 2.197 0.721
## Elachistocleis_cesari Hypsiboas_albopunctatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 2.012 0.755 2.069 0.755
## Hypsiboas_faber Hypsiboas_lundii
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 0.913 0.859 0.915 0.859
## Hypsiboas_prasinus Itapotihyla_langsdorffii
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 0.913 0.859 0.911 0.859
## Leptodactylus_fuscus Leptodactylus_latrans
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 4.652 0.309 0.811 0.859
## Leptodactylus_mystaceus Leptodactylus_mystacinus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 1.869 0.776 9.772 0.045
## Physalaemus_centralis Physalaemus_cuvieri
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 2.639 0.644 8.221 0.062
## Physalaemus_marmoratus Physalaemus_nattereri
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 1.622 0.776 17.131 0.014
## Pseudis_paradoxa Rhinella_ornata
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 2.67 0.482 11.961 0.027
## Rhinella_schneideri Scinax_fuscomarginatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 3.363 0.464 2.64 0.639
## Scinax_fuscovarius Scinax_similis
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 1.655 0.776 2.646 0.612
## Trachycephalus_typhonius
## Dev Pr(>Dev)
## (Intercept)
## grupos 1.622 0.776
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.Alguns dos exercícios propostos no capítulo
1
Dados
data(mite)
data(mite.env)
data(mite.xy)RDA
## Passo 1: transformação de hellinger da matriz de espécies
# caso tenha dados de abundância.
acaro.hel <- decostand(x = mite, method = "hellinger")
## Passo 2: selecionar variáveis importantes
# Para isso, é necessário remover a variável categórica.
acaro.conti <- mite.env[, -c(3:6)]
## Evite usar variáveis muito correlacionadas
sel.vars.acaro <- forward.sel(acaro.hel, acaro.conti)## Testing variable 1
## Testing variable 2sel.vars.acaro$variables## [1] "WatrCont" "SubsDens"env.sel.acaro <- mite.env[,sel.vars.acaro$variables]
## Passo 3: padronizar matriz ambiental (somente variáveis contínuas)
env.pad.acaro <- decostand(x = env.sel.acaro, method = "standardize")
## RDA com dados selecionados e padronizados
rda.acaro <- rda(acaro.hel ~ WatrCont + SubsDens, data = env.pad.acaro)
## Para interpretar, é necessário saber a significância dos eixos para representar a relação entre as variáveis preditoras e a composição de espécies
res.axis.acaro <- anova.cca(rda.bird, by = "axis")
res.axis.acaro## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = species.hel ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)
## Df Variance F Pr(>F)
## RDA1 1 0.045759 12.0225 0.001 ***
## RDA2 1 0.009992 2.6252 0.057 .
## RDA3 1 0.007518 1.9752 0.112
## RDA4 1 0.003582 0.9410 0.488
## Residual 18 0.068510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Relevância da dimensão 1 e pouco da 2.
## Em seguida, é possível identificar quais são as variáveis que contribuem ou que mais contribuem para a variação na composição de espécies
res.var.acaro <- anova.cca(rda.acaro, by = "term") ## Qual variável?
res.var.acaro## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = acaro.hel ~ WatrCont + SubsDens, data = env.pad.acaro)
## Df Variance F Pr(>F)
## WatrCont 1 0.107069 27.0254 0.001 ***
## SubsDens 1 0.021768 5.4945 0.001 ***
## Residual 67 0.265439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ambas as variáveis são contribuem.
r_quadr <- RsquareAdj(rda.acaro)
r_quadr## $r.squared
## [1] 0.326768
##
## $adj.r.squared
## [1] 0.3066715Indício de forte contribuição da variáveis escolhidas.
## Ordenação multi-escala (MSO) para entender os resultados da ordenação em relação à distância geográfica
acaro.rda <- mso(rda.acaro, mite.xy, grain = 1, permutations = 99)
msoplot(acaro.rda)
## Error variance of regression model underestimated by -1.1 percentVários sinais de autocorrelação.
ggord(rda.acaro, ptslab = TRUE, size = 1, addsize = 3, repel = TRUE) +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 0, linetype = 2) +
tema_livro()
RDAp (com MEM)
# Passo 1: Gerar um arquivo LIST W list binária de vizinhança
mat_knn.acaro <- knearneigh(as.matrix(mite.xy), k = 2, longlat = FALSE)
mat_nb.acaro <- knn2nb(mat_knn.acaro, sym = TRUE)
mat_listw.acaro <- nb2listw(mat_nb.acaro, style = "W")
mat_listw.acaro## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 70
## Number of nonzero links: 182
## Percentage nonzero weights: 3.714286
## Average number of links: 2.6
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 70 4900 70 55.55 289.85## Passo 2: Listar os métodos "candidatos" para obter a matriz SWM
MEM_mat.acaro <- scores.listw(mat_listw.acaro, MEM.autocor = "positive")
candidates.acaro <- listw.candidates(mite.xy, nb = c("gab", "mst", "dnear"),
weights = c("binary", "flin"))## Warning in nb2listw(nb.object, style = style, glist = lapply(nb.dist, f1, : zero
## sum general weights
## Warning in nb2listw(nb.object, style = style, glist = lapply(nb.dist, f1, : zero
## sum general weights## Passo 3: Selecionar a melhor matriz SWM e executar o MEM
W_sel_mat.acaro <- listw.select(acaro.hel, candidates.acaro, MEM.autocor = "positive", p.adjust = TRUE, method = "FWD")## Procedure stopped (alpha criteria): pvalue for variable 14 is 0.051000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 16 is 0.053000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 13 is 0.067000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 14 is 0.052000 (> 0.050000)
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.451035 with 11 variables (> 0.445310)
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.464620 with 13 variables (> 0.463888)## Passo 4: Matriz dos preditores espaciais escolhidos (MEMs)
spatial.pred.acaro <- as.data.frame(W_sel_mat.acaro$best$MEM.select)
## É necessário atribuir os nomes das linhas
rownames(spatial.pred.acaro) <- rownames(mite.xy)
## Combinar variáveis ambientais e espaciais em um único data.frame
pred.vars.acaro <- data.frame(env.pad.acaro, spatial.pred.acaro)
## RDA parcial
rda.p.acaro <- rda(acaro.hel ~ WatrCont + SubsDens + # Preditores ambientais
Condition(MEM4 + MEM2 + MEM9 + MEM1 + MEM3 + MEM6 + MEM7+ MEM13 + MEM31 + MEM17 + MEM8 + MEM30 + MEM5 + MEM26 + MEM11), # Preditores espaciais
data = pred.vars.acaro)
## Interpretação
# Para interpretar, é necessário saber a significância dos eixos para representar a relação entre as variáveis preditoras e a composição de espécies.
res.p.axis.acaro <- anova.cca(rda.p.acaro, by = "axis")
res.p.axis.acaro## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = acaro.hel ~ WatrCont + SubsDens + Condition(MEM4 + MEM2 + MEM9 + MEM1 + MEM3 + MEM6 + MEM7 + MEM13 + MEM31 + MEM17 + MEM8 + MEM30 + MEM5 + MEM26 + MEM11), data = pred.vars.acaro)
## Df Variance F Pr(>F)
## RDA1 1 0.012459 4.6299 0.001 ***
## RDA2 1 0.003495 1.2989 0.188
## Residual 52 0.139929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dimensão 1 significativa.
## Contribuição
# Em seguida, é possível identificar quais são as variáveis que contribuem ou que mais contribuem para a variação na composição de espécies.
res.p.var.acaro <- anova.cca(rda.p.acaro, by = "term") ## Qual variável?
res.p.var.acaro## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = acaro.hel ~ WatrCont + SubsDens + Condition(MEM4 + MEM2 + MEM9 + MEM1 + MEM3 + MEM6 + MEM7 + MEM13 + MEM31 + MEM17 + MEM8 + MEM30 + MEM5 + MEM26 + MEM11), data = pred.vars.acaro)
## Df Variance F Pr(>F)
## WatrCont 1 0.012028 4.4698 0.001 ***
## SubsDens 1 0.003926 1.4591 0.129
## Residual 52 0.139929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1WatrCont como variável contribuinte.
RsquareAdj(rda.p.acaro)## $r.squared
## [1] 0.04046451
##
## $adj.r.squared
## [1] 0.03426293R2 mostrando relação fraca.
dbRDA
## Passo 1: transformação de hellinger da matriz de espécies
## Passo 2: gerar matriz de distância (Diversidade beta: Bray-Curtis)
dbeta.acaro <- vegdist(acaro.hel, "bray")
## Passo 2: dbRDA
dbrda.acaro <- capscale (dbeta.acaro~WatrCont+SubsDens,
data=env.pad.acaro)
# Para interpretar, é necessário saber a significância dos eixos para representar a relação entre as variáveis preditoras e a composição de espécies
db.res.axis.acaro <- anova.cca(dbrda.acaro, by = "axis")
db.res.axis.acaro## Permutation test for capscale under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = dbeta.acaro ~ WatrCont + SubsDens, data = env.pad.acaro)
## Df SumOfSqs F Pr(>F)
## CAP1 1 3.4918 30.501 0.001 ***
## CAP2 1 0.4180 3.651 0.001 ***
## Residual 67 7.6701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dimensões 1 e 2 significantes.
# Em seguida, é possível identificar quais são as variáveis que contribuem ou que mais contribuem para a variação na composição de espécies
db.res.var.acaro <- anova.cca(dbrda.acaro, by = "term") ## Qual variável?
db.res.var.acaro## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = dbeta.acaro ~ WatrCont + SubsDens, data = env.pad.acaro)
## Df SumOfSqs F Pr(>F)
## WatrCont 1 3.2215 28.1400 0.001 ***
## SubsDens 1 0.6883 6.0122 0.001 ***
## Residual 67 7.6701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ambas variáveis contribuem para variação.
db_r_quadr <- RsquareAdj(dbrda.acaro)
db_r_quadr## $r.squared
## [1] 0.3951217
##
## $adj.r.squared
## [1] 0.3770656PERMANOVA
## Composição de espécies padronizar com método de Hellinger
## Matriz de distância com método Bray-Curtis
sps.dis.acaro <- vegdist(x = acaro.hel, method = "bray")
## Verifica correlação entre as variáveis
ggpairs(mite.env) +
tema_livro()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Após verificar a estrutura de correlação, vamos manter somente três variáveis
env2.acaro <- mite.env[, c("SubsDens", "Topo", "Shrub")]
## PERMANOVA
perm.aves <- adonis2(sps.dis.acaro ~ SubsDens + Topo + Shrub, data = env2.acaro)
perm.aves## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = sps.dis.acaro ~ SubsDens + Topo + Shrub, data = env2.acaro)
## Df SumOfSqs R2 F Pr(>F)
## SubsDens 1 0.4017 0.04059 4.2754 0.003 **
## Topo 1 1.9808 0.20018 21.0844 0.001 ***
## Shrub 2 1.4061 0.14210 7.4834 0.001 ***
## Residual 65 6.1065 0.61713
## Total 69 9.8950 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## BETADISPER
betad.acaro <- betadisper(sps.dis.acaro, env2.acaro$SubsDens)
permutest(betad.acaro)## Warning in anova.lm(lm(Distances ~ Groups, data = model.dat)): ANOVA F-tests on
## an essentially perfect fit are unreliable## Warning in summary.lm(mod): essentially perfect fit: summary may be unreliable##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 68 0.10929 0.0016072 4.1725e+30 999 0.409
## Residuals 1 0.00000 0.0000000## É homogêneo, resultado do PERMANOVA explica em parte a variação de composisão das espéciesComparando os resultados, vemos que eles são semelhantes. o RDA, RDAp, dbRDA e PERMANOVA, todos indicam que há uma relação entre as variáveis ambientais e a composição de espécies dos locais testados.
2
Método UPGA e índice de Bray-Curtis
## Matriz de similaridade com o índice de Bray-Curtis
distacaro.b <- vegdist(x = mite, method = "bray")
## Agrupamento com a função hclust e o método UPGMA
dendroacaro.b <- hclust(d = distacaro.b, method = "average")
## Visualizar os resultados
plot(dendroacaro.b, main = "Dendrograma",
ylab = "Similaridade (índice de Horn)",
xlab="", sub="")
Índice de Morisita-Horn
## Matriz de similaridade com o índice de Morisita-Horn
distacaro.h <- vegdist(x = mite, method = "horn")
## Agrupamento com a função hclust e o método UPGMA
dendroacaro.h <- hclust(d = distacaro.h, method = "average")
## Visualizar os resultados
plot(dendroacaro.h, main = "Dendrograma",
ylab = "Similaridade (índice de Horn)",
xlab="", sub="")
Os dendogramas são diferentes e poderiam acabar por atribuir similaridade a grupos mais distantes.