Resumo Multivariada
Valeria
2022-10-08
Agrupamento hierárquico
head(sp_compos)## BP4 PP4 PP3 AP1 AP2 PP1 PP2 BP9 PT1 PT2 PT3 BP2 PT5
## Aper 0 3 0 0 2 0 0 0 0 0 0 181 0
## Bahe 859 14 14 0 87 312 624 641 0 0 0 14 0
## Rict 1772 1517 207 573 796 0 0 0 0 0 0 0 0
## Cleuco 0 0 0 0 0 0 0 0 0 29 369 0 84
## Dmic 0 0 6 60 4 0 0 0 2758 319 25 0 329
## Dmin 0 84 344 1045 90 0 0 0 8 0 0 0 0distBocaina <- vegdist(x = sp_compos, method = "horn")
dendro <- hclust(d = distBocaina, method = "average")
plot(dendro, main = "Dendrograma",
ylab = "Similaridade (Ãndice de Horn)",
xlab="", sub="")
Avliando qualidade do dendograma
cofresult <- cophenetic(dendro)
cor(cofresult, distBocaina)## [1] 0.9455221Coeficiente de correlação >0.7 indica uma boa representação, portanto o dendograma é adequado
plot(dendro, main = "Dendrograma",
ylab = "Similaridade (Ãndice de Horn)",
xlab="", sub="")
k <- 4
n <- ncol(sp_compos)
MidPoint <- (dendro$height[n-k] + dendro$height[n-k+1]) / 2
abline(h = MidPoint, lty=2)
É possivel observar um agrupamento de 5 grupos distintos abaixo da linha de corte e 6 quando se observa no total.
Exemplo 2
Distancia de Chord
head(t(sp_compos))## Aper Bahe Rict Cleuco Dmic Dmin Hpoly Lfur Pbar Polf Pmelano Sduar Shay
## BP4 0 859 1772 0 0 0 61 3 387 0 0 0 0
## PP4 3 14 1517 0 0 84 275 0 187 0 0 1150 6
## PP3 0 14 207 0 6 344 388 0 0 0 0 428 0
## AP1 0 0 573 0 60 1045 1054 0 0 0 0 476 92
## AP2 2 87 796 0 4 90 3002 0 0 2 0 7 0
## PP1 0 312 0 0 0 0 329 0 0 0 0 0 0
## Sobt Ssqual Bokerm
## BP4 0 0 0
## PP4 0 0 1
## PP3 0 0 0
## AP1 13 0 0
## AP2 5 0 0
## PP1 0 0 0bocaina_transf <- disttransform(t(sp_compos), "chord")
analise <- pvclust(bocaina_transf, method.hclust = "average", method.dist = "euclidean", quiet = TRUE)
plot(analise, hang=-1, main = "Dendrograma com valores de P",
ylab = "Distância Euclideana",
xlab="", sub="")
pvrect(analise)
Distância de Hellinger
bocaina_transf2 <- disttransform(t(bocaina), "hellinger")
analise2 <- pvclust(bocaina_transf2, method.hclust="average", method.dist="euclidean", quiet = TRUE)
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)
Agrupameto não-hierárquico (k-means)
head(doubs$fish)[,1:6]## Cogo Satr Phph Neba Thth Teso
## 1 0 3 0 0 0 0
## 2 0 5 4 3 0 0
## 3 0 5 5 5 0 0
## 4 0 4 5 5 0 0
## 5 0 2 3 2 0 0
## 6 0 3 4 5 0 0rowSums(doubs$fish)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 3 12 16 21 34 21 16 0 14 14 11 18 19 28 33 40 44 42 46 56 62 72 4 15 11 43
## 27 28 29 30
## 63 70 87 89spe <- doubs$fish[-8,]
spe.norm <- decostand(x = spe, method = "normalize")
spe.kmeans <- kmeans(x = spe.norm, centers = 4, nstart = 100)
spe.kmeans## K-means clustering with 4 clusters of sizes 3, 8, 12, 6
##
## Cluster means:
## Cogo Satr Phph Neba Thth Teso
## 1 0.00000000 0.000000000 0.00000000 0.00000000 0.000000000 0.000000000
## 2 0.00000000 0.006691097 0.02506109 0.06987391 0.006691097 0.006691097
## 3 0.10380209 0.542300691 0.50086515 0.43325916 0.114024105 0.075651573
## 4 0.06167791 0.122088022 0.26993915 0.35942538 0.032664966 0.135403325
## Chna Chto Lele Lece Baba Spbi Gogo
## 1 0.05205792 0.00000000 0.07647191 0.3166705 0.00000000 0.0000000 0.20500174
## 2 0.10687104 0.09377516 0.14194394 0.2011411 0.24327992 0.1326062 0.28386032
## 3 0.00000000 0.00000000 0.06983991 0.1237394 0.02385019 0.0000000 0.05670453
## 4 0.06212775 0.21568957 0.25887226 0.2722562 0.15647062 0.1574388 0.16822286
## Eslu Pefl Rham Legi Scer Cyca Titi
## 1 0.07647191 0.00000000 0.0000000 0.05205792 0.07647191 0.00000000 0.00000000
## 2 0.20630360 0.16920496 0.2214275 0.19066542 0.13171275 0.16019126 0.26230024
## 3 0.04722294 0.02949244 0.0000000 0.00000000 0.00000000 0.00000000 0.03833408
## 4 0.12276089 0.17261621 0.0793181 0.06190283 0.04516042 0.06190283 0.14539027
## Abbr Icme Acce Ruru Blbj Alal Anan
## 1 0.00000000 0.0000000 0.18058775 0.31667052 0.05205792 0.7618709 0.00000000
## 2 0.19561641 0.1331835 0.26713081 0.32103755 0.22883055 0.3326939 0.18873077
## 3 0.00000000 0.0000000 0.00000000 0.01049901 0.00000000 0.0000000 0.00000000
## 4 0.01473139 0.0000000 0.03192175 0.32201597 0.01473139 0.1095241 0.04739636
##
## 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
## 3 3 3 3 4 3 3 4 3 3 3 3 3 3 4 4 4 4 2 2 2 1 1 1 2 2
## 28 29 30
## 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 0.3560423 0.4696535 2.5101386 1.7361453
## (between_SS / total_SS = 66.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Avaliar se esse é o número ideiais de grupos
spe.KM.cascade <- cascadeKM(spe.norm, inf.gr = 2, sup.gr = 10, iter = 100, criterion = "ssi")
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.9523667 2.4840549
## ssi 0.1312111 0.1685126 0.1233266 0.1110514 0.07121006 0.1107331 0.1343565
## 9 groups 10 groups
## SSE 2.0521888 1.75992916
## ssi 0.1264591 0.07208002plot(spe.KM.cascade, sortg = TRUE)
De acordo com a análise foi possivel observar que o número ideal de grupos é 3 pois o SSI maximo é alcansado nesse número
Espécies indicadoras
head(bocaina)## BP4 PP4 PP3 AP1 AP2 PP1 PP2 BP9 PT1 PT2 PT3 BP2 PT5
## Aper 0 3 0 0 2 0 0 0 0 0 0 181 0
## Bahe 859 14 14 0 87 312 624 641 0 0 0 14 0
## Rict 1772 1517 207 573 796 0 0 0 0 0 0 0 0
## Cleuco 0 0 0 0 0 0 0 0 0 29 369 0 84
## Dmic 0 0 6 60 4 0 0 0 2758 319 25 0 329
## Dmin 0 84 344 1045 90 0 0 0 8 0 0 0 0fitofis <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4), rep(5, 4))
res_indval <- indval(t(sp_compos), fitofis)
summary(res_indval)## cluster indicator_value probability
## Rict 1 0.8364 0.016
## Sduar 1 0.7475 0.048
##
## Sum of probabilities = 8.053
##
## Sum of Indicator Values = 7.3
##
## Sum of Significant Indicator Values = 1.58
##
## Number of Significant Indicators = 2
##
## Significant Indicator Distribution
##
## 1
## 2tab_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 1.000
## Bahe 2 0.6487329 0.058
## Rict 1 0.8363823 0.016
## Cleuco 3 0.4128631 0.429
## Dmic 3 0.6645244 0.211
## Dmin 1 0.7032145 0.081
## Hpoly 2 0.6208711 0.255
## Lfur 3 0.2279412 1.000
## Pbar 1 0.2813725 0.609
## Polf 3 0.2437500 1.000
## Pmelano 2 0.2500000 1.000
## Sduar 1 0.7474527 0.048
## Shay 3 0.4930269 0.415
## Sobt 2 0.2222222 0.706
## Ssqual 3 0.2500000 1.000
## Bokerm 2 0.4583333 0.225tab_indval[tab_indval$P < 0.05, ]## maxcls ind.value P
## Rict 1 0.8363823 0.016
## Sduar 1 0.7474527 0.048Análise de ordenação
Ordenação irrestrita
Análise de Componentes Principais (PCA)
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))
aranha.cent <- as.data.frame(base::scale(aranhas, center = TRUE, scale=FALSE))
matriz_cov <- cov(aranha.cent)
eigen_aranhas <- eigen(matriz_cov)
autovalores <- eigen_aranhas$values
autovetores <- as.data.frame(eigen_aranhas$vectors)
autovalores## [1] 36.733031 27.243824 9.443805 2.962749 1.438020colnames(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.3544591matriz_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.5630683100 * (autovalores/sum(autovalores))## [1] 47.201691 35.008126 12.135225 3.807112 1.847846ggplot(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()
Exemplo 1
sum(is.na(penguins))## [1] 19penguins <- na.omit(penguins)
penguins_trait <- penguins[, 3:6]
penguins_trait %>%
dplyr::summarise(across(where(is.numeric),
~var(.x, na.rm = TRUE)))## # A tibble: 1 × 4
## comprimento_bico profundidade_bico comprimento_nadadeira massa_corporal
## <dbl> <dbl> <dbl> <dbl>
## 1 29.9 3.88 196. 648372.penguins_pad <- decostand(x = penguins_trait, method = "standardize")
penguins_pad %>%
dplyr::summarise(across(where(is.numeric),
~var(.x, na.rm = TRUE)))## comprimento_bico profundidade_bico comprimento_nadadeira massa_corporal
## 1 1 1 1 1pca.p <- PCA(X = penguins_trait, scale.unit = TRUE, graph = FALSE)
pca.p$eig ## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 2.7453557 68.633893 68.63389
## comp 2 0.7781172 19.452929 88.08682
## comp 3 0.3686425 9.216063 97.30289
## comp 4 0.1078846 2.697115 100.00000fviz_screeplot(pca.p, addlabels = TRUE, ylim = c(0, 70), main = "",
xlab = "Dimensões",
ylab = "Porcentagem de variância explicada") 
var_env <- get_pca_var(pca.p)
var_env$coord ## Dim.1 Dim.2 Dim.3 Dim.4
## comprimento_bico 0.7518288 0.52943763 -0.3900969 -0.04768208
## profundidade_bico -0.6611860 0.70230869 0.2585287 0.05252186
## comprimento_nadadeira 0.9557480 0.00510580 0.1433474 0.25684871
## massa_corporal 0.9107624 0.06744932 0.3592789 -0.19204478var_env$contrib ## Dim.1 Dim.2 Dim.3 Dim.4
## comprimento_bico 20.58919 36.023392267 41.279994 2.107420
## profundidade_bico 15.92387 63.388588337 18.130600 2.556942
## comprimento_nadadeira 33.27271 0.003350291 5.574092 61.149849
## massa_corporal 30.21423 0.584669105 35.015313 34.185789var_env$cor ## Dim.1 Dim.2 Dim.3 Dim.4
## comprimento_bico 0.7518288 0.52943763 -0.3900969 -0.04768208
## profundidade_bico -0.6611860 0.70230869 0.2585287 0.05252186
## comprimento_nadadeira 0.9557480 0.00510580 0.1433474 0.25684871
## massa_corporal 0.9107624 0.06744932 0.3592789 -0.19204478var_env$cos2## Dim.1 Dim.2 Dim.3 Dim.4
## comprimento_bico 0.5652466 2.803042e-01 0.15217561 0.002273581
## profundidade_bico 0.4371669 4.932375e-01 0.06683710 0.002758546
## comprimento_nadadeira 0.9134542 2.606919e-05 0.02054847 0.065971260
## massa_corporal 0.8294881 4.549411e-03 0.12908133 0.036881196ind_env <- get_pca_ind(pca.p)Variáveis importante eixo 1
dimdesc(pca.p)$Dim.1 ##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## comprimento_nadadeira 0.9557480 5.962756e-178
## massa_corporal 0.9107624 3.447018e-129
## comprimento_bico 0.7518288 7.830597e-62
## profundidade_bico -0.6611860 3.217695e-43Variáveis importantes eixo 2
dimdesc(pca.p)$Dim.2##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## profundidade_bico 0.7023087 8.689230e-51
## comprimento_bico 0.5294376 1.873918e-25As 4 variáves morfologicas estão associadas com o eixo 1 explicando 69%, entretanto o comprimento da nadadeira, massa corporal e comprimento do bico estão associados positivamente e largura do bico negativamente, já o eixo 2 explica 20% sendo relacionado com a largura e comprimento do bico
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()
Análises de Coordenadas Principais (PCoA)
##Exempolo 1
mite.hel <- decostand(x = mite, method = "hellinger")
sps.dis <- vegdist(x = mite.hel, method = "bray")
pcoa.sps <- pcoa(D = sps.dis, correction = "cailliez")
100 * (pcoa.sps$values[, 1]/pcoa.sps$trace)[1]## [1] 49.10564100 * (pcoa.sps$values[, 1]/pcoa.sps$trace)[2]## [1] 14.30308sum(100 * (pcoa.sps$values[, 1]/pcoa.sps$trace)[1:2])## [1] 63.40872eixos <- pcoa.sps$vectors[, 1:2]
pcoa.dat <- data.frame(topografia = mite.env$Topo, eixos)
eixos <- pcoa.sps$vectors[, 1:2]
pcoa.dat <- data.frame(topografia = mite.env$Topo, eixos)
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()
Regressão de Componentes Principais (PCR)
Exemplo 1
env_cont <- env[,-8]
env.pca <- PCA(env_cont, scale.unit = TRUE, graph = FALSE)
var_env <- get_pca_var(env.pca)
var_env$contrib ## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## mini.jan 10.93489 22.2975487 16.1607726 7.6025527 0.01782438
## maxi.jan 20.18065 3.2890767 2.1814486 4.2756350 41.05646526
## mini.jul 11.87396 21.1379132 0.3428843 0.7750666 44.70209396
## maxi.jul 18.47244 0.9159957 56.5369988 9.4368661 2.59283074
## rain.jan 9.95206 21.5387403 6.5737927 53.7375738 4.44283706
## rain.jul 16.14997 11.2368132 7.2608047 19.6972097 0.71454880
## rain.tot 12.43603 19.5839121 10.9432983 4.4750959 6.47339980var_env$cor ## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## mini.jan 0.6830371 0.6766524 -0.21924927 0.12298817 -0.004517369
## maxi.jan 0.9279073 0.2598807 -0.08055260 0.09223249 0.216804944
## mini.jul 0.7117620 0.6588220 0.03193603 -0.03926930 -0.226225907
## maxi.jul 0.8877675 0.1371462 0.41008461 -0.13702428 0.054483561
## rain.jan -0.6516187 0.6650391 -0.13983474 -0.32698110 0.071319550
## rain.jul -0.8300858 0.4803509 0.14696011 0.19796389 -0.028601865
## rain.tot -0.7284135 0.6341424 0.18041856 0.09435932 0.086088397ind_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.00000pred.env <- ind_env$coord[, 1:3]
riqueza <- specnumber(species)
dat <- data.frame(pred.env, riqueza)
mod1 <- lm(riqueza ~ Dim.1 + Dim.2 + Dim.3, data = dat)
par(mfrow = c(2, 2))
plot(mod1) 
summary(mod1)##
## 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-07ggplot(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()## `geom_smooth()` using formula 'y ~ x'
Exemplo 2
env.dist <- gowdis(mite.env)
env.mite.pco <- pcoa(env.dist, correction = "cailliez")
100 * (env.mite.pco$values[, 1]/env.mite.pco$trace)[1]## [1] 61.49635100 * (env.mite.pco$values[, 1]/env.mite.pco$trace)[2]## [1] 32.15486pred.scores.mite <- env.mite.pco$vectors[, 1:2]
mite.riqueza <- specnumber(mite)
pred.vars <- data.frame(riqueza = mite.riqueza, pred.scores.mite)
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-07g_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()
grid.arrange(g_acari_axi1, g_acari_axi2, nrow = 1)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Ordenação restrita
Análise de Redundância (RDA)
Exemplo 1
species.hel <- decostand(x = species, method = "hellinger")
env.contin <- env[, -8]
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.183000 (> 0.050000)sel.vars$variables## [1] "rain.jul" "maxi.jul"env.sel <- env[,sel.vars$variables]
env.pad <- decostand(x = env.sel, method = "standardize")
env.pad.cat <- data.frame(env.pad, altitude = env$altitude)
rda.bird <- rda(species.hel ~ rain.jul + maxi.jul + altitude, data = env.pad.cat)
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.066 .
## RDA3 1 0.007518 1.9752 0.124
## RDA4 1 0.003582 0.9410 0.467
## Residual 18 0.068510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res.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.012 *
## altitude 2 0.019071 2.5053 0.009 **
## Residual 18 0.068510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r_quadr <- RsquareAdj(rda.bird)
r_quadr## $r.squared
## [1] 0.4938685
##
## $adj.r.squared
## [1] 0.3813949bird.rda <- mso(rda.bird, xy, grain = 1, permutations = 99)
msoplot(bird.rda)
## Error variance of regression model underestimated by -2.1 percentggord(rda.bird, ptslab = TRUE, size = 1, addsize = 3, repel = TRUE) +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 0, linetype = 2) +
tema_livro()
Análise de Redundância parcial (RDAp)
head(species.hel)[, 1:6] ## Fauvette_orphee Fauvette_des_jardins Fauvette_a_tete_noire
## S01 0 0.3651484 0.3651484
## S02 0 0.3333333 0.3333333
## S03 0 0.3162278 0.3162278
## S04 0 0.4200840 0.3429972
## S05 0 0.3872983 0.3162278
## S06 0 0.3779645 0.3779645
## Fauvette_babillarde Fauvette_grisette Fauvette_pitchou
## S01 0.2581989 0.2581989 0
## S02 0.2357023 0.2357023 0
## S03 0.3162278 0.2236068 0
## S04 0.2425356 0.0000000 0
## S05 0.2236068 0.3162278 0
## S06 0.2672612 0.0000000 0head(xy)## x y
## S01 156 252
## S02 141 217
## S03 171 233
## S04 178 215
## S05 123 189
## S06 154 195head(env.pad.cat)## rain.jul maxi.jul altitude
## S01 1.333646 0.1462557 Montanhoso
## S02 1.468827 -0.6848206 Intermediário
## S03 1.505694 -0.2099199 Montanhoso
## S04 1.296778 -2.0699476 Montanhoso
## S05 1.075572 -0.3682201 Plano
## S06 1.100151 -0.6056705 Intermediáriomat_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"))
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.129000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.066000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.069000 (> 0.050000)
## Procedure stopped (alpha criteria): pvalue for variable 4 is 0.148000 (> 0.050000)spatial.pred <- as.data.frame(W_sel_mat$best$MEM.select)
rownames(spatial.pred) <- rownames(xy)
pred.vars <- data.frame(env.pad.cat, spatial.pred)
rda.p <- rda(species.hel ~
rain.jul + maxi.jul + altitude +
Condition(MEM1 + MEM2 + MEM4 + MEM5),
data = pred.vars)
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.322
## RDA2 1 0.004830 1.2189 0.799
## RDA3 1 0.003240 0.8176 0.891
## RDA4 1 0.001891 0.4773 0.932
## Residual 14 0.055477res.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.375
## maxi.jul 1 0.004446 1.1220 0.348
## altitude 2 0.009579 1.2087 0.220
## Residual 14 0.055477RsquareAdj(rda.p)## $r.squared
## [1] 0.1361661
##
## $adj.r.squared
## [1] 0.02330319pca.comp <- dudi.pca(species.hel, scale = FALSE, scannf = FALSE)
moran.comp <- moran.mc(pca.comp$li[, 1], mat_listw, 999)
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.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: greatermoran.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: greaterpv.birds <- varpart(species.hel, env.pad.cat, spatial.pred)
plot(pv.birds, cutoff = -Inf)
mtext("Diagrama de Venn: partição de variância")
Análise de Redundância baseada em distância (db-RDA)
species.hel <- decostand(species, "hellinger")
dbeta <- vegdist(species.hel, "bray")
dbrda.bird <- capscale (dbeta~rain.jul+maxi.jul+altitude, data=env.pad.cat)
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.102
## CAP3 1 0.018141 1.1205 0.656
## CAP4 1 0.009500 0.5868 0.780
## Residual 18 0.291420
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1db.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.048 *
## altitude 2 0.094966 2.9328 0.017 *
## Residual 18 0.291420
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1db_r_quadr <- RsquareAdj(dbrda.bird)
db_r_quadr## $r.squared
## [1] 0.6116758
##
## $adj.r.squared
## [1] 0.5253816PERMANOVA
Exemplo 1
species.hel <- decostand(x = species, method = "hellinger")
sps.dis <- vegdist(x = species.hel, method = "bray")
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`.
env2 <- env[, c("mini.jan", "rain.tot", "altitude")]
perm.aves <- adonis2(sps.dis ~ mini.jan + rain.tot + altitude, data = env2)
perm.aves ### Diferenças entre os tratamentos?## 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.005 **
## rain.tot 1 0.12910 0.22771 9.0118 0.001 ***
## altitude 2 0.08929 0.15749 3.1163 0.018 *
## Residual 18 0.25787 0.45483
## Total 22 0.56695 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betad.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.217
## Residuals 20 0.0290636 0.0014532as.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.0000000sol1 <- metaMDS(sps.dis)## Run 0 stress 0.1344042
## Run 1 stress 0.1344042
## ... Procrustes: rmse 7.120076e-05 max resid 0.0002396814
## ... Similar to previous best
## Run 2 stress 0.1272417
## ... New best solution
## ... Procrustes: rmse 0.07890806 max resid 0.3526067
## Run 3 stress 0.1344042
## Run 4 stress 0.1378506
## Run 5 stress 0.1272418
## ... Procrustes: rmse 4.458154e-05 max resid 0.0001702944
## ... Similar to previous best
## Run 6 stress 0.1272417
## ... New best solution
## ... Procrustes: rmse 6.480023e-06 max resid 2.014141e-05
## ... Similar to previous best
## Run 7 stress 0.1344042
## Run 8 stress 0.1344043
## Run 9 stress 0.1336482
## Run 10 stress 0.1344042
## Run 11 stress 0.1272417
## ... New best solution
## ... Procrustes: rmse 1.008143e-05 max resid 3.76644e-05
## ... Similar to previous best
## Run 12 stress 0.1336481
## Run 13 stress 0.1272418
## ... Procrustes: rmse 2.409013e-05 max resid 9.030595e-05
## ... Similar to previous best
## Run 14 stress 0.1338432
## Run 15 stress 0.1344042
## Run 16 stress 0.1338432
## Run 17 stress 0.1272418
## ... Procrustes: rmse 0.0001110389 max resid 0.0004272811
## ... Similar to previous best
## Run 18 stress 0.1344042
## Run 19 stress 0.1336481
## Run 20 stress 0.1336481
## *** Solution reachednmds.beta <- metaMDS(sps.dis, previous.best = sol1)## Starting from 2-dimensional configuration
## Run 0 stress 0.1272417
## Run 1 stress 0.1336481
## Run 2 stress 0.1272418
## ... Procrustes: rmse 2.916685e-05 max resid 7.315048e-05
## ... Similar to previous best
## Run 3 stress 0.1338432
## Run 4 stress 0.1338432
## Run 5 stress 0.1272417
## ... Procrustes: rmse 4.089858e-06 max resid 1.177283e-05
## ... Similar to previous best
## Run 6 stress 0.1344042
## Run 7 stress 0.1272417
## ... Procrustes: rmse 8.554531e-06 max resid 3.174032e-05
## ... Similar to previous best
## Run 8 stress 0.1344042
## Run 9 stress 0.1272418
## ... Procrustes: rmse 4.107754e-05 max resid 0.0001565599
## ... Similar to previous best
## Run 10 stress 0.1630389
## Run 11 stress 0.1338432
## Run 12 stress 0.1344042
## Run 13 stress 0.1272417
## ... Procrustes: rmse 4.523109e-06 max resid 1.309694e-05
## ... Similar to previous best
## Run 14 stress 0.1584787
## Run 15 stress 0.1659927
## Run 16 stress 0.1378506
## Run 17 stress 0.1272417
## ... Procrustes: rmse 1.544736e-05 max resid 5.822014e-05
## ... Similar to previous best
## Run 18 stress 0.1336481
## Run 19 stress 0.1378506
## Run 20 stress 0.1378506
## *** Solution reachednmds.beta$stress # valor ideal entre 0 e 0.2## [1] 0.1272417dat.graf <- data.frame(nmds.beta$points, altitude = env2$altitude)
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")]), ]
hull.data <- rbind(grp.mon, grp.int, grp.pla)
head(hull.data)## MDS1 MDS2 altitude
## S04 -0.10578350 -0.10682811 Montanhoso
## S01 -0.25332376 0.04198712 Montanhoso
## S11 -0.12504929 0.14477160 Montanhoso
## S15 0.09165956 0.09857158 Montanhoso
## S18 0.01968236 -0.12417321 Intermediário
## S06 -0.16053832 -0.08924391 Intermediárioggplot(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 Mantentel parcial
Exemplo 1
head(bocaina.env)## pH tur DO temp veg dep canopy conduc hydro area
## PP3 7.07 0.66 2.26 18.93 30 0.35 30.3 0.0100 per 139
## PP4 6.85 0.33 2.78 20.53 0 0.33 32.1 0.0100 per 223
## AP1 6.64 0.40 4.09 15.26 5 3.50 42.0 0.0105 per 681
## AP2 6.13 0.80 3.47 15.40 70 2.50 82.0 0.0300 per 187
## PP1 5.80 0.75 3.22 16.75 70 0.70 91.4 0.0100 per 83
## PP2 5.32 4.00 2.10 15.90 95 0.27 75.3 0.0100 per 108head(bocaina.xy)## Long Lat
## PT5 -44.6239 -22.7243
## PP3 -44.6340 -22.7228
## PP4 -44.6334 -22.7228
## BP4 -44.6212 -22.7262
## BP2 -44.6155 -22.7289
## PT1 -44.6170 -22.7350Dist.km <- as.dist(rdist.earth(bocaina.xy, miles=F))
comp.bo.pad <- decostand(t(bocaina), "hellinger")
env.bocaina <- decostand(bocaina.env[,-9], "standardize")
dissimil.bocaina <- vegdist(comp.bo.pad, "bray")
dissimil.env <- vegdist(env.bocaina, "euclidian")
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.005
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.264 0.337 0.389 0.430
## Permutation: free
## Number of permutations: 999matrix.dist.env <- data.frame(x = melt(as.matrix(Dist.km))$value,
y = melt(as.matrix(dissimil.env))$value)## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUEggplot(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() 
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.012
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.136 0.219 0.282 0.387
## Permutation: free
## Number of permutations: 999matrix.dist.bocaina <- data.frame(x = melt(as.matrix(Dist.km))$value,
y = melt(as.matrix(dissimil.bocaina))$value)## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUEggplot(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()
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.015
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.169 0.205 0.252 0.331
## Permutation: free
## Number of permutations: 999matrix.bocaina.env <- data.frame(x = melt(as.matrix(dissimil.bocaina))$value,
y = melt(as.matrix(dissimil.env))$value)## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the
## caller; using TRUEggplot(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()
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.146
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.159 0.212 0.251 0.294
## Permutation: free
## Number of permutations: 999Mantel espacial com modelo nulo restrito considerando 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.008
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 3.159665393 0.002249728 0.013828362nb.boc <- mst.nb(Dist.km)
plot(nb.boc, bocaina.xy)
lw <- nb2listw(nb.boc)
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.118
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 1.28418890 0.16441415 0.02658694compos_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.008
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 3.159665393 0.002249728 0.013828362PROCRUSTES e PROTEST
Exemolo 1
head(fish_comm)## sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10 sp11 sp12 sp13 sp14 sp15 sp16
## site1 0 0 0 1 1 0 4 0 0 0 1 0 1 1 0 5
## site2 0 0 0 1 0 1 5 0 2 1 1 2 1 0 0 1
## site3 0 1 0 4 0 0 0 1 5 0 5 0 0 5 0 1
## site4 0 5 0 1 0 1 0 0 2 0 0 0 0 2 1 0
## site5 0 2 0 0 0 0 0 0 0 0 4 0 0 1 0 4
## site6 0 2 1 0 0 0 0 1 1 0 0 2 0 2 0 0
## sp17 sp18 sp19 sp20 sp21 sp22 sp23 sp24
## site1 0 0 0 0 0 1 0 3
## site2 0 0 0 0 0 0 1 0
## site3 0 1 0 0 2 0 0 0
## site4 0 0 0 1 1 0 0 2
## site5 0 0 0 1 5 0 5 0
## site6 0 0 0 0 0 0 4 0head(macroinv)## sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10 sp11 sp12 sp13 sp14 sp15 sp16
## site1 0 1 1 0 0 0 0 0 0 1 2 0 0 0 0 1
## site2 0 2 1 0 0 0 0 1 1 0 0 2 0 1 1 0
## site3 0 0 0 1 0 1 5 0 2 1 1 2 0 4 0 1
## site4 0 1 0 4 0 0 0 1 5 0 5 0 0 5 0 1
## site5 0 2 0 0 0 0 0 0 0 0 4 0 0 2 1 0
## site6 1 0 0 1 0 0 0 0 0 0 1 0 0 2 0 0
## sp17 sp18 sp19 sp20 sp21 sp22 sp23 sp24 sp25 sp26 sp27 sp28 sp29 sp30
## site1 1 0 4 0 0 0 1 0 0 1 0 4 0 0
## site2 0 0 0 0 0 1 2 0 1 0 1 1 0 0
## site3 0 0 1 0 0 0 0 1 0 2 1 0 0 0
## site4 0 1 0 0 2 0 0 0 0 2 1 0 0 0
## site5 0 0 0 1 1 0 0 2 0 5 0 1 0 1
## site6 0 0 0 0 0 0 4 0 0 2 0 0 0 0
## sp31 sp32 sp33 sp34 sp35 sp36
## site1 0 1 5 0 5 0
## site2 1 0 0 1 0 0
## site3 0 1 1 0 0 2
## site4 0 0 0 5 1 0
## site5 0 0 2 0 0 0
## site6 0 0 0 0 4 0set.seed(1001)
d_macro <- vegdist(macroinv, "bray")
pcoa_macro <- cmdscale(d_macro)
d_fish <- vegdist(fish_comm, "bray")
pcoa_fish <- cmdscale(d_fish)
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: 999plot(concord, main = "",
xlab = "Dimensão 1",
ylab = "Dimensão 2")
Métodos multivariados baseados em modelos
head(anuros_permanova)[1:6]## Chiasmocleis_albopunctata Dendropsophus_elianae
## Araras 6 0
## Artur Nogueira 0 0
## Conchal 0 0
## Engenheiro Coelho 1 0
## Jaguari\xfana 8 0
## Limeira 0 11
## Dendropsophus_jimi Dendropsophus_nanus Dendropsophus_minutus
## Araras 0 0 1
## Artur Nogueira 0 0 6
## Conchal 0 0 0
## Engenheiro Coelho 0 0 2
## Jaguari\xfana 0 0 0
## Limeira 15 0 15
## Dendropsophus_sanborni
## Araras 0
## Artur Nogueira 0
## Conchal 0
## Engenheiro Coelho 0
## Jaguari\xfana 0
## Limeira 0anuros_abund <- as.data.frame(anuros_permanova[,-28])
grupos <- factor(anuros_permanova[,28])
abund_tr <- mvabund(anuros_abund)
meanvar.plot(abund_tr)
Exemplo 1
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).
modelo1 <- manyglm(abund_tr ~ grupos, family = "negative.binomial")
plot(modelo1)
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(modelo1, p.uni = "adjusted")## Time elapsed: 0 hr 0 min 7 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.468 0.394 0.844
## Dendropsophus_jimi Dendropsophus_nanus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 0.23 0.844 2.197 0.707
## Dendropsophus_minutus Dendropsophus_sanborni
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 4.752 0.309 2.197 0.707
## Elachistocleis_cesari Hypsiboas_albopunctatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 2.012 0.759 2.069 0.759
## Hypsiboas_faber Hypsiboas_lundii
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 0.913 0.844 0.915 0.844
## Hypsiboas_prasinus Itapotihyla_langsdorffii
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 0.913 0.844 0.911 0.844
## Leptodactylus_fuscus Leptodactylus_latrans
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 4.652 0.309 0.811 0.844
## Leptodactylus_mystaceus Leptodactylus_mystacinus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 1.869 0.786 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.786 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.468 2.64 0.639
## Scinax_fuscovarius Scinax_similis
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## grupos 1.655 0.786 2.646 0.639
## Trachycephalus_typhonius
## Dev Pr(>Dev)
## (Intercept)
## grupos 1.622 0.786
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.