Descrição.

Esta página tem por objetivo disponibilizar as soluções feitos dos exercícios do curso “Ecologia Numérica” ofertado pelo bacharelado em Ciências Biológicas do Centro de Biologia da UFPE, no segundo período de 2021.

Exercício 1 - Introdução

Atenção: Neste exercício trabalhamos com dados fictícios.

O objetivo desse exercício é estudar o comportamento dos índices de diversidade \(\alpha\), \(\beta\) e \(\gamma\) de uma comunidade de bromélias em função da altura em relação ao nível do mar.

Desenho Experimental

  1. Os dados representam o número de indivíduos de diferentes espécies de bromélias detectados em um transecto ao longo de uma montanha.

  2. Foram registrados, com o uso um altímetro, o número de espécies de bromélias em cinco pontos igualmente espaçados (300 metros de altura entre si) ao longo do transecto.

  3. As espécies foram registradas e o número de indivíduos encontrados de cada espécie foi tabelado para cada ponto medido do transecto.

  4. Os cincos pontos foram numerados proporcionalmente a altura de ocorrência. Assim, o ponto 1 representa o ponto mais baixo e o pontos 5 o mais alto.

  5. Em cada ponto, observou-se uma faixa de 40m x 40m .

Dados

Número de espécies observadas em cada ponto de observação.

Dados Experimentais
V1 V2 V3 V4 V5 V6
Espécie Altura 1 Altura 2 Altura 3 Altura 4 Altura 5
Sp. 1 5 3 4 6 6
Sp. 2 0 1 0 0 0
Sp. 3 12 10 15 17 19
Sp. 4 0 1 0 4 8
Sp. 5 0 0 1 0 1
Sp. 6 7 11 12 17 23
Sp. 7 0 0 0 2 3
Sp. 8 0 0 4 7 13
Sp. 9 0 2 10 12 15
Sp. 10 4 0 5 9 14
Sp. 11 0 0 0 3 7
Sp. 12 3 10 14 18 19

Resultados

** Número de espécies observados em cada altura:**

Número de Espécies
Alt. 1 Alt. 2 Alt. 3 Alt. 4 Alt. 5
5 7 8 10 11

A seguir apresentamos o gráfico da diversidade \(\alpha\)

A diversidade \(\gamma\) pode ser calculada como sendo \[ \gamma = 12 \]

Exercício 2 - Descrição de Comunidades

O objetivo desse exercício é estudar o comportamento dos índices de diversidade \(\alpha\), \(\beta\) e \(\gamma\) de uma comunidade de bromélias em função da altura em relação ao nível do mar.

Descrição dos Dados

O conjunto de dados cobre 75 espécies de Corais, coletados ao longo de 10 transectos totalizando 60 locais de coletas na Indonésia. OS dados foram coletados ao longo dos anos 1981-1988 e encontram-se publicamente disponíveis. Para maiores detalhes ver https://www.davidzeleny.net/doku.php.

Região do Experimento.

## Carregando pacotes exigidos: permute
## Carregando pacotes exigidos: lattice
## This is vegan 2.5-7
##  Y81R1  Y81R2  Y81R3  Y81R4  Y81R5  Y81R6  Y81R7  Y81R8  Y81R9 Y81R10  Y83R1 
##     30     27     26     21     16     14     14     14     11      7      6 
##  Y83R2  Y83R3  Y83R4  Y83R5  Y83R6  Y83R7  Y83R8  Y83R9 Y83R10  Y84R1  Y84R2 
##      8      6      6      3      2      1      1      2      1     15     13 
##  Y84R3  Y84R4  Y84R5  Y84R6  Y84R7  Y84R8  Y84R9 Y84R10  Y85R1  Y85R2  Y85R3 
##     12      7      8      6      3      4      3      3     19     13     15 
##  Y85R4  Y85R5  Y85R6  Y85R7  Y85R8  Y85R9 Y85R10  Y87R1  Y87R2  Y87R3  Y87R4 
##     11      9      7      6      6      5      4     12     12     12      9 
##  Y87R5  Y87R6  Y87R7  Y87R8  Y87R9 Y87R10  Y88R1  Y88R2  Y88R3  Y88R4  Y88R5 
##      5      5      4      3      2      2     14     12     13      5      6 
##  Y88R6  Y88R7  Y88R8  Y88R9 Y88R10 
##      5      5      4      5      4
##  Y81R1  Y81R2  Y81R3  Y81R4  Y81R5  Y81R6  Y81R7  Y81R8  Y81R9 Y81R10  Y83R1 
##  167.2  106.9  106.5   80.1   62.0   51.0   58.1   45.5   36.6   47.9    8.6 
##  Y83R2  Y83R3  Y83R4  Y83R5  Y83R6  Y83R7  Y83R8  Y83R9 Y83R10  Y84R1  Y84R2 
##   14.5   10.3   10.8    5.5    4.6    2.5    2.4    2.7    3.0   25.9   25.7 
##  Y84R3  Y84R4  Y84R5  Y84R6  Y84R7  Y84R8  Y84R9 Y84R10  Y85R1  Y85R2  Y85R3 
##   23.4   16.5   18.4   16.0    8.3   15.3   10.9   36.2   50.8   45.9   61.6 
##  Y85R4  Y85R5  Y85R6  Y85R7  Y85R8  Y85R9 Y85R10  Y87R1  Y87R2  Y87R3  Y87R4 
##   41.5   58.8   41.3   25.2   42.7   44.6   62.8   39.1   50.5   34.4   23.1 
##  Y87R5  Y87R6  Y87R7  Y87R8  Y87R9 Y87R10  Y88R1  Y88R2  Y88R3  Y88R4  Y88R5 
##   20.4   30.9   28.8   35.1   78.5   88.1   53.7   48.5   47.8    9.2   15.9 
##  Y88R6  Y88R7  Y88R8  Y88R9 Y88R10 
##   17.1   11.4   10.5   92.9   35.7

Análise de Dados

Curva de Acúmulo para todas as Amostras

Ajuste do Modelo log.normal para a curva de abundância (todas as amostras)

Ajuste do Modelo log.normal para a curva de abundância (10 primeiras amostras)

Ajuste do vários modelos para a primeira amostra

Ordenamento de Espécies

## Carregando pacotes exigidos: tcltk
## BiodiversityR 2.14-1: Use command BiodiversityRGUI() to launch the Graphical User Interface; 
## to see changes use BiodiversityRGUI(changeLog=TRUE, backward.compatibility.messages=TRUE)

COntribuição de Cada Espécie

Resultados

a <- colnames(tikus.spe)
b <- colSums(tikus.spe)
aa<-a[1:3]
total_species <- colSums(tikus.spe)
total_biomass <- sum(total_species)
Percentual <- total_species/total_biomass

dados_totais <- t(data.frame(total_species, Percentual, a[1:5]))

A1=sort(Percentual,decreasing = T)

plot(A1,type="o",xlim=c(1,5),ylim=c(0,0.3),lwd=2,col="red",pch=1,cex=2,xlab="Rank",ylab="Abundância Percentual")
text(A1, labels= a[1:3] , cex=0.9, font=2, pos=4)

O gráfico nos mostra que a espécie Psammocora contigua é sozinha, responsável por 22% de toda a biomassa

biomassa_cumulativo=integer(75)
Percentual=sort(Percentual,decreasing = T)
biomassa_cumulativo[1]=Percentual[1]
for(i in 1:74){ 
     biomassa_cumulativo[i+1] =biomassa_cumulativo[i] + Percentual[i+1]  }
plot(biomassa_cumulativo,type="o",xlim=c(1,75),ylim=c(0,1.0),lwd=2,col="darkblue",pch=1,cex=2,main="Abundância Percentual Acumulada",xlab="Rank",ylab="Abundância Acumulada")

hist(Percentual, 
     main="Histograma das Contribuições Percentuais por Espécie", 
     xlab="Passengers", 
     border="blue", 
     col="gray",
     xlim=c(0,0.3),
     las=1, 
     breaks=25)

Pelo histograma e pelo gráfico de abundância podemos concluir que:

  1. As 5 espécies mais abundantes respondem por 49% da biomassa.

  2. Grande parte das espécies são “raras”. De fato 50%, espécies são responsáveis por apenas 20% da biomassa total.

Exercício 3 - Índices de Diversidade

Índices de Diversidade

Curvas de Abundância

## Carregando pacotes exigidos: usethis
## 
## Attaching package: 'devtools'
## The following object is masked from 'package:permute':
## 
##     check
## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

## Warning in qt(0.975, df = n - 1): NaNs produzidos

Índices de Diversidade

# Riqueza
specnumber(composicao_especies)
##  Com_1  Com_2  Com_3  Com_4  Com_5  Com_6  Com_7  Com_8  Com_9 Com_10 
##     10     10      5      5      5      6      2      4      6      4
# Índice de Shannon
diversity(composicao_especies)
##     Com_1     Com_2     Com_3     Com_4     Com_5     Com_6     Com_7     Com_8 
## 2.3025851 0.5002880 0.9580109 1.6068659 1.4861894 1.5607038 0.6931472 1.1058899 
##     Com_9    Com_10 
## 1.7140875 1.2636544
# Índice de Simpson
diversity(composicao_especies,index="simpson")
##     Com_1     Com_2     Com_3     Com_4     Com_5     Com_6     Com_7     Com_8 
## 0.9000000 0.1710000 0.4814815 0.7989636 0.7587500 0.7674858 0.5000000 0.5850000 
##     Com_9    Com_10 
## 0.8088889 0.6942149
Índice Com. 1 Com.2 Com.3 Com. 4 Com. 5 Com. 6 Com.7 Com. 8 Com. 9 Com. 10
Riqueza 10 10 5 5 5 6 2 4 6 4
Shannon 2.30 0.50 0.99 1.60 1.49 1.59 0.69 1.11 1.71 1.26
Simpson 0.90 0.17 0.48 0.80 0.76 0.77 0.50 0.59 0.81 0.69

Variáveis Ambientais

** Número de espécies observados em cada altura:**

## Analysis of Variance Table
## 
## Response: Riqueza
##              Df Sum Sq Mean Sq F value  Pr(>F)  
## Precipitacao  1 30.622 30.6224  8.9156 0.01744 *
## Residuals     8 27.478  3.4347                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
## 
## Response: Shannon
##              Df  Sum Sq Mean Sq F value Pr(>F)
## Precipitacao  1 0.10989 0.10989  0.3627 0.5637
## Residuals     8 2.42366 0.30296
## Analysis of Variance Table
## 
## Response: Simpson
##              Df  Sum Sq  Mean Sq F value Pr(>F)
## Precipitacao  1 0.00132 0.001325  0.0252 0.8778
## Residuals     8 0.42064 0.052580
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

Discussão

Pela análise feita podemos inferir que:

  1. Os índices de Shannon e de Simpson detectam de forma satisfatória a cauda da curva de abundância. Quanto mais lento for o decrescimento da curva de abundância maiores serão os índices. Assim, números altos dos índices indicam uma distribuição de espécies mais uniforme.

  2. A riqueza não está relacionada a nenhum dos índices, nem ao formato da curva de abundância.

  3. A riqueza foi o único índice a se correlacionar de forma estatisticamente significante com a precipitação.

log ## {-}

Exercício 4 - Diversidade Verdadeira

Neste exercício trabalhamos com dados disponibilizados no artigo

Ribeiro-Neto, J.D., Arnan, X., Tabarelli, M. et al. 

Chronic anthropogenic disturbance causes homogenization of plant and ant communities in the Brazilian Caatinga.

Biodivers Conserv 25, 943–956 (2016).

https://doi.org/10.1007/s10531-016-1099-5

Dados: Descrição

Os dados referem-se a presença de espécies de plantas na Caatinga. Foram divididas de acordo com os tipos de solos em que se encontravam: Areia e Argila. O objetivo deste exercício é comparar os índices de diversidade encontrados em cada tipo de planta. Os índices são dados pela fórmula \[ {^q} D = \left(\sum_{i=1}^{s}p_i^q\right)^{\frac{1}{1-q}}\] onde \(p_i\) é o percentual da espécie \(i\), \(s\) o número de espécies (riqueza) e \(q\) o parametro do índice.

library(readxl)
plantas_caatinga <- read_excel("C:/Users/casti/Documents/Trabalho/Biologia/Sexto Periodo/Ecologia Numerica/Programas/plantas_caatinga.xlsx")
head(plantas_caatinga)
## # A tibble: 6 x 3
##   Especie                    Areia Argila
##   <chr>                      <dbl>  <dbl>
## 1 Myracrodruon urundeuva        25     21
## 2 Schinopsis brasiliensis        8      4
## 3 Spondias tuberose              1      1
## 4 Aspidosperma pyrifolium      129     54
## 5 Handroanthus impetiginosus    16      0
## 6 Handroanthus spongiosus        0      8

Índices

library (vegan)
library(ggplot2)
library(BiodiversityR)
library(hillR)

a <- seq(from = 0, to = 2, by = 0.1)
vet_areia <- integer(21)
vet_argila <- integer(21)

contador=1
for ( indice in a){  vet_areia[contador]=hill_taxa(plantas_caatinga$Areia, q  = indice)
                     vet_argila[contador]=hill_taxa(plantas_caatinga$Argila, q  = indice)
                     
                    contador=contador+1
                     }
plot(a, vet_areia, type="o", col="blue", pch="o", lty=1, ylim=c(5,50),xlab="q",ylab="Números de Hill",main="RIQUEZA X TIPO DE SOLO" )
points(a, vet_argila, col="red", pch="+")
lines(a, vet_argila, col="red",lty=2)
legend(1.5,40,legend=c("Areia","Argila"), col=c("blue","red"),pch=c("o","+"),lty=c(1,2), ncol=1)

Discussão

Lembrando que os valores dos índices

Parâmetro q Índice
0 Riqueza
1 Shannon (exponencial)
2 Simpson

Pela curva temos que:

  1. A argila apresenta maior riqueza de espécies

  2. A areia apresenta índice de Shannon e de Simpson maiores.

  3. A Argila possui valores de índices maiores para \(q\) pequeno e para $ q > 0.5$ possui valores menores.

Inferimos que

  1. A Argila possui maior quantidade total de espécies (riqueza) bem como uma quantidade maior de espécies menos comuns.

  2. A Areia apresenta um maior número de espécies mais comuns (Shannon e Simpson)

Conluimos assim que a preservação de áreas de Argila é fundamental para a preservação da diversidade de plantas da Caatinga.

Exercício 5 - Diversidade Beta

Neste exercício trabalhamos com dados disponibilizados no artigo

Plant β-diversity in fragmented rain forests: testing floristic homogenization and differentiation hypotheses

Víctor Arroyo-Rodríguez, Matthias Rös, Federico Escobar, Felipe P. L. Melo, Bráulio A. Santos, Marcelo Tabarelli, Robin Chazdon

Journal of Ecology Volume101, Issue6 November 2013 Pages 1449-1458

https://doi.org/10.1111/1365-2745.12153

Dados: Descrição

Região de estudo: Temos aqui três paisagens caracterizadas pelo nível de desmatamento: LDL, IDL e HDL - LDL = LOW deforestain level, IDL = INTERMEDIATE deforestation level e HDL = HIGH deforestation level. Região do Experimento.

library(readxl)
dados<-read.csv("https://raw.githubusercontent.com/fplmelo/ecoa/main/content/en/courses/eco_num/betadiv/com_ltx_all.csv", row.names = "X")
 dados<-as.data.frame(dados)
library(entropart)

Índices - Cálculo dos perfis.

Lembre que :

\(\alpha\)-diversidade:

Parâmetro q Índice
0 Riqueza
1 Shannon (exponencial)
2 Simpson

\(\gamma\) diversidade - Riqueza de espécies

\(\beta\)-diversidade: Definida por \[ \gamma = \bar{\alpha} \, \beta \] onde \(\bar \alpha\) é a média das \(\alpha\)-diversidade.

# Gerando o Perfil
mc<-MetaCommunity(dados)
Profile <- DivProfile(q.seq = seq(0, 2, 0.1), mc, Biased = FALSE, Correction = "None")
plot(Profile)

Análise das Comunidades

library(betapart)
# Comunidade LDL 
dadosLDL<-dados[, 1:12]
dadosLDL<-ifelse(dadosLDL=="0",0,1) # Tranformei em 0 e 1

beta.coreL <-betapart.core(dadosLDL)
beta.multiL <-beta.multi(dadosLDL)

# Comunidade IDL 
dadosIDL<-dados[, 13:24]
dadosIDL<-ifelse(dadosIDL=="0",0,1) # Tranformei em 0 e 1

beta.coreI <-betapart.core(dadosIDL)
beta.multiI <-beta.multi(dadosIDL)

# Comunidade HDL 
dadosHDL<-dados[, 25:36]
dadosHDL<-ifelse(dadosHDL=="0",0,1) # Tranformei em 0 e 1

beta.coreH <-betapart.core(dadosHDL)
beta.multiH <-beta.multi(dadosHDL)
Índice Interpretação
1 = beta.SIM Turnover (substituição)
2= beta.SNE Aninhamento
3 = beta.SOR Beta diversidade total

Discussão

Pelo gráfico temos que:

Os índices são semelhantes nter as três comunidades (?).

Baixo aninhamento para as três comunidades.

Alta diversidade beta e alto Turnover para as três comunidades.

Assim temos que: AS comunidades possuem espécies bem diferentes entre si dado o alto valor de beta e a alta substituição. Esse efeito é intensificado pelo baixo valor do aninhamento.

Exercício 6 - Diversidade Funcional

Neste exercício estudamos o conceito de Diversidade Funcional de acordo com o capítulo 8 do livro “Introdução ao R com Aplicações em Biodiversidade e Conservação”.

Carregando os Dados

# Pacotes necessários
#library(FD)
#library(ade4)
#library(ecodados)
#library(gridExtra)
#library(ggplot2)
#library(ggrepel)
#library(tidyverse)
#library(picante)

# Dados Usados :

# devtools::install_github("paternogbc/ecodados")

comun_fren_dat <- ecodados::fundiv_frenette2012a_comu
ambie_fren_dat <- ecodados::fundiv_frenette2012a_amb
trait_fren_dat <- ecodados::fundiv_frenette2012a_trait
trait_dat      <- ecodados::fundiv_barbaro2009a_trait
comun_dat      <- ecodados::fundiv_barbaro2009a_comu
ambie_dat      <- ecodados::fundiv_barbaro2009a_amb

Os dados carregados são como seguem:

comun_fren_dat : matriz de abundância 55 localidades (linhas) x 36 espécies (colunas)

comun_dat: matriz de abundância 101 localidades (linhas) x 36 espécies (colunas).

ambie_frent_dat: matriz de características ambientais 46 localidades (linhas) x 5 fatores ambientais (colunas)

ambie_dat: matriz de características ambientais 101 localidades (linhas) x 11 fatores ambientais (colunas)

trait_dat: matriz de características funcionais 36 espécies (linhas) x 12 tratos(colunas)

comun_dat: matriz de abundância 101 localidades (linhas) x 36 espécies (colunas)

Análise da Diversidade Funcional

library(ape)
## 
## Attaching package: 'ape'
## The following object is masked from 'package:ggpubr':
## 
##     rotate
trait_cont <- trait_fren_dat
trait_pad <- decostand(trait_cont, "standardize")
euclid_dis <- vegdist(trait_pad, "euclidean")

pcoa_traits_cont <- pcoa(euclid_dis, correction="cailliez") # Resultados são idênticos ao de uma PCA
eixos_cont <- as.data.frame(pcoa_traits_cont$vectors[,1:2]) # Selecionar os dois primeiros eixos

eixos_cont %>% 
  ggplot(aes(x=Axis.1, y=Axis.2)) + 
  geom_point(pch=21, size=3, color = "black", fill="#4575b4") + 
   geom_text_repel(aes(Axis.1, Axis.2, label = rownames(eixos_cont))) +
  xlab("PCO 1") + ylab("PCO 2") + 
  theme(axis.title.x = element_text(face="bold", size=14),
        axis.text.x = element_text(vjust=0.5, size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=14),
        axis.text.y = element_text(vjust=0.5, size=12)) + 
  geom_hline(yintercept = 0, linetype=2) + 
  geom_vline(xintercept = 0, linetype=2)+ 
  theme(legend.position = "top", legend.title=element_blank())+ ggtitle("Diversidade Funcional: trait_fren_dat - Componentes Principais/Euclideana") +
  theme(plot.title = element_text(hjust = 0.5)) -> plot_trait_cont
plot_trait_cont

print("Usando somente as varáveis categóricas de trait_dat para o cálculo da
distância funcional via métrica de Gowder.")
## [1] "Usando somente as varáveis categóricas de trait_dat para o cálculo da\ndistância funcional via métrica de Gowder."
library(FD)
## Carregando pacotes exigidos: ade4
## Carregando pacotes exigidos: geometry
trait_dat %>% 
  dplyr::select_if(is.character) -> trait_cat 
dist_categ <- gowdis(trait_cat)

pcoa_traits_cat <- pcoa(dist_categ, correction="cailliez")
eixos_cat <- as.data.frame(pcoa_traits_cat$vectors[,1:2]) # Selecionar os dois primeiros eixos

eixos_cat %>% 
  ggplot(aes(x=Axis.1, y=Axis.2)) + 
  geom_point(pch=21, size=3, color = "black", fill="#4575b4") + 
   geom_text_repel(aes(Axis.1, Axis.2, label = rownames(eixos_cat))) +
  xlab("PCO 1") + ylab("PCO 2") + 
  theme(axis.title.x = element_text(face="bold", size=14),
        axis.text.x = element_text(vjust=0.5, size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=14),
        axis.text.y = element_text(vjust=0.5, size=12)) + 
  geom_hline(yintercept = 0, linetype=2) + 
  geom_vline(xintercept = 0, linetype=2)+ 
  theme(legend.position = "top", legend.title=element_blank())+ ggtitle("Diversidade Funcional: trait_dat - Componentes Principais/Gowder") +
  theme(plot.title = element_text(hjust = 0.5)) -> plot_trait_cat
plot_trait_cat

eixos_cat %>% 
  ggplot(aes(x=Axis.1, y=Axis.2)) + 
  geom_point(pch=21, size=3, color = "black", fill="#4575b4") + 
   geom_text_repel(aes(Axis.1, Axis.2, label = rownames(eixos_cat))) +
  xlab("PCO 1") + ylab("PCO 2") + 
  theme(axis.title.x = element_text(face="bold", size=14),
        axis.text.x = element_text(vjust=0.5, size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=14),
        axis.text.y = element_text(vjust=0.5, size=12)) + 
  geom_hline(yintercept = 0, linetype=2) + 
  geom_vline(xintercept = 0, linetype=2)+ 
  theme(legend.position = "top", legend.title=element_blank()) -> plot_trait_cat2



print("Verificando os tipos das variáveis funcionais:")
## [1] "Verificando os tipos das variáveis funcionais:"
trait_dat %>% 
  dplyr::summarise_all(class) %>% 
  tidyr::gather(variable, class)
##    variable     class
## 1     trend character
## 2   redlist character
## 3     regio   integer
## 4      biog character
## 5     activ character
## 6      diet character
## 7    winter character
## 8     color character
## 9     breed character
## 10     body   integer
## 11     wing character
## 12   period character
print("Separando a matriz de características  funcionais em duas: uma ordinal e outra categórica.")
## [1] "Separando a matriz de características  funcionais em duas: uma ordinal e outra categórica."
trait_dat$regio <- as.ordered(trait_dat$regio)
trait_dat$body <- as.ordered(trait_dat$body)

## Combinar cada tipo de trait em um data.frame separado
# Categóricos

trait_categ <- cbind.data.frame(trend=trait_dat$trend, redlist=trait_dat$redlist, biog=trait_dat$biog, activ=trait_dat$activ,  diet=trait_dat$diet, winter=trait_dat$winter,color=trait_dat$color, breed=trait_dat$breed,wing=trait_dat$wing, period=trait_dat$period)

# Ordinais
trait_ord <- cbind.data.frame(regio=trait_dat$regio, body=trait_dat$body)
rownames(trait_categ) <- rownames(trait_dat)
rownames(trait_ord) <- rownames(trait_dat)

# Agora, combinar os dois data.frames em uma lista chamada "ktab"
ktab_list <- ktab.list.df(list(trait_categ, trait_ord))

dist_mist <- dist.ktab(ktab_list, type= c("N", "O"))


print(" Visualisando via Componentes Principais:")
## [1] " Visualisando via Componentes Principais:"
pcoa_traits_mist <- pcoa(dist_mist, correction="cailliez")
eixos_mist <- as.data.frame(pcoa_traits_mist$vectors[,1:2])  # exportar 2 eixos

eixos_mist %>% 
  ggplot(aes(x=Axis.1, y=Axis.2)) + 
  geom_point(pch=21, size=3, color = "black", fill="#d73027") + 
  geom_text_repel(aes(Axis.1, Axis.2, label = rownames(eixos_mist)))+
  xlab("PCO 1") + ylab("PCO 2") + 
  theme(axis.title.x = element_text(face="bold", size=14),
        axis.text.x = element_text(vjust=0.5, size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=14),
        axis.text.y = element_text(vjust=0.5, size=12)) + 
  geom_hline(yintercept = 0, linetype=2) + 
  geom_vline(xintercept = 0, linetype=2)+ 
  theme(legend.position = "top", legend.title=element_blank()) -> plot_trait_mist
plot_trait_mist

print("Comparando a distância funcional entre as comunidades usando-se 
somente os tratos categóricos e outra com os dados categóricos e os dados (numéricos).")
## [1] "Comparando a distância funcional entre as comunidades usando-se \nsomente os tratos categóricos e outra com os dados categóricos e os dados (numéricos)."
library(gridExtra)
library(grid)
library(lattice)
library(ggrepel)
library(ggplot2)

grid.arrange(plot_trait_cat2, plot_trait_mist, ncol=2)

print("Calculando a Riqueza Funcional : parte do convex hull funcional disponvel utilizado pela comunidade")
## [1] "Calculando a Riqueza Funcional : parte do convex hull funcional disponvel utilizado pela comunidade"
richness <- dbFD(dist_mist, comun_dat)$nbsp
## FRic: Dimensionality reduction was required. The last 17 PCoA axes (out of 19 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.3243851 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
plot(richness,xlab="Comunidade",ylab="Riqueza",main="Riqueza Funcional Alpha ")

print("Calculando a Diversidade Funcional")
## [1] "Calculando a Diversidade Funcional"
library(picante)
## Carregando pacotes exigidos: nlme
dend <- hclust(dist_mist, "average")


tree_dend <-as.phylo(dend)



FD <- pd(comun_dat, tree_dend)$PD


plot(FD,xlab="Comunidade",ylab="Diversidade",main="Diversidade Funcional")

plot(FD,richness,xlab="Diversidade Funcional",ylab="Riqueza",main="Diversidade x Riqueza")

fdiv <- dbFD(dist_mist, comun_dat)$FDiv
## FRic: Dimensionality reduction was required. The last 17 PCoA axes (out of 19 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.3243851 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
plot(FD,xlab="Comunidade",ylab="Divergência",main="Divergência Funcional")

feve <- dbFD(dist_mist, comun_dat)$FEve
## FRic: Dimensionality reduction was required. The last 17 PCoA axes (out of 19 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.3243851 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
plot(feve,xlab="Comunidade",ylab="Regularidade",main="Regularidade Funcional")

locais <- rownames(comun_dat)

fric <-dbFD(dist_mist, comun_dat)$FRic
## FRic: Dimensionality reduction was required. The last 17 PCoA axes (out of 19 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.3243851 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
fdis <- dbFD(dist_mist, comun_dat)$FDis
## FRic: Dimensionality reduction was required. The last 17 PCoA axes (out of 19 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.3243851 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
metricas <- data.frame(richness=richness,
                       FD_gp = FD,
                       fric = fric,
                       fdiv = fdiv,
                       fdis = fdis,
                       feve = feve)
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(metricas)+ ggtitle("Correlação entre as Métricas: Diagonal = Histogramas Suavisados") + theme(plot.title = element_text(hjust = 0.5))

library(betapart)
library(vegan)

### Partição da Diversidade beta (Método Baselga)

comun_fren_pa <- as.matrix(decostand(comun_fren_dat, "pa"))
trait <- as.matrix(trait_fren_dat)
 
rowSums(comun_fren_pa)>ncol(trait)
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   17 
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 
##   18   19   20   21   22   23   24   25   26   27   28   29   31   32   33   34 
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 
##   41   42   43   44   45   46   47   48   49   50   51   53   54   55 
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
cwm_fren <- functcomp(trait_pad, as.matrix(comun_fren_dat))

plot(cwm_fren, main="Composição das Características Funcionais plotadas 2 a 2")

# Como antes mas usando a métrica euclideana.

#Passo 1: calcular a distância funcional

trait_pad <- decostand(trait_fren_dat, "standardize")
euclid_dis <- vegdist(trait_pad, "euclidean")


## Passo 2: calcular a Divergência funcional (FDis) e Regularidade Funcional (FEve)

fdis <- dbFD(euclid_dis, comun_fren_dat)$FDis# Fdis=0 em locais com somente uma espécie
## FRic: No dimensionality reduction was required. All 5 PCoA axes were kept as 'traits'. 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
feve <- dbFD(euclid_dis, comun_fren_dat)$FEve
## FRic: No dimensionality reduction was required. All 5 PCoA axes were kept as 'traits'. 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
lm_dat <- data.frame(aridez = ambie_fren_dat$Aridity, fdis = fdis, feve = feve)

mod1 <- lm(fdis ~ aridez, data = lm_dat)
plot(mod1) 

anova(mod1)
## Analysis of Variance Table
## 
## Response: fdis
##           Df Sum Sq Mean Sq F value Pr(>F)
## aridez     1 0.2083 0.20834  0.9945 0.3241
## Residuals 44 9.2179 0.20950
mod2 <- lm(feve ~ aridez, data = lm_dat)
plot(mod2)

anova(mod2)
## Analysis of Variance Table
## 
## Response: feve
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## aridez     1 0.02098 0.020979  1.0447 0.3123
## Residuals 44 0.88353 0.020080
## Passo 4: gráfico para visualizar os dois resultados  

lm_dat %>% 
  ggplot(aes(x=aridez, y=fdis)) + 
  geom_point(pch=21, size=3, color = "black", fill="royalblue") + 
  xlab("Aridez") + ylab("Divergência Funcional (FDis)") + 
  theme(axis.title.x = element_text(face="bold", size=14),
        axis.text.x = element_text(vjust=0.5, size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=14),
        axis.text.y = element_text(vjust=0.5, size=12)) + 
  theme(legend.position = "top", legend.title=element_blank()) -> plot_pred1

lm_dat %>% 
  ggplot(aes(x=aridez, y=feve)) + 
  geom_point(pch=21, size=3, color = "black", fill="#d73027") + 
  xlab("Aridez") + ylab("Regularidade Funcional (FEve)") + 
  theme(axis.title.x = element_text(face="bold", size=14),
        axis.text.x = element_text(vjust=0.5, size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=14),
        axis.text.y = element_text(vjust=0.5, size=12)) + 
  theme(legend.position = "top", legend.title=element_blank()) -> plot_pred2

grid.arrange(plot_pred1, plot_pred2, ncol=2)

# Hipótese: o pastejo determina a ocorrência de espécies de plantas com diferentes atributos funcionais 
# Predição 1: a composição funcional de plantas é diferente entre áreas com e sem pastejo


## Passo 1: calcular a distância funcional

cwm_dis <- vegdist(cwm_fren, "euclidean")

## Passo 2: calcular os valores de composição funcional (CWM)

cwm_fren <- functcomp(trait_pad, as.matrix(comun_fren_dat))

## Passo 3: testar se a composição funcional varia entre as áreas com uma PERMANOVA

perman_fren <- adonis(cwm_fren~Grazing, data = ambie_fren_dat)

## Passo 4: comparar a variação dentro de cada grupo com Betadisper 

betad_fren <- betadisper(cwm_dis, ambie_fren_dat$Grazing)
permutest(betad_fren)
## 
## 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     1  0.0539 0.053858 0.1946    999  0.687
## Residuals 44 12.1763 0.276735
plot(betad_fren)

print("Como p > 0.05 a diferença entre as comunidades NÃO é estatísticamente relevante.")
## [1] "Como p > 0.05 a diferença entre as comunidades NÃO é estatísticamente relevante."