###Código em R para aálise de atributos funcionais
### Fundamentos em Ecologia 2020 - Grupo 2A e 2B


#carregando bibliotecas que serão necessárias para as análises
# bibliotecas são conjuntos de funções pré desenvolvidas, mas específicas para determinadas tarefas
# estas bibliotecas precisam ser instaladas, caso não tenha sido instaladas antes
# por exemplo = install.packages("ggpubr")
library(readxl)
library(dplyr)
library(lattice)
library(ggpubr)

Dataset_Fundamentos_Grupo_A <- read_excel("D:/Desktop/Dataset_Fundamentos_Novo.xls")
head(Dataset_Fundamentos_Grupo_A)
names(Dataset_Fundamentos_Grupo_A)

# criando novas colunas para a tabela de dados
# LMA é a razão massa - área, e é obtida dividindo-se a massa da folha por sua área
# também convertemos de g/cm2 para g/ m-2, multiplicando por 10000
Dataset_Fundamentos_Grupo_A$LMA <- Dataset_Fundamentos_Grupo_A$Sample.weight.g / Dataset_Fundamentos_Grupo_A$Leaf.Area.cm2 * 10000

# N_mass é a quantidade de nitrogênio foliar expresso em unidades de massa (mg de N por grama de folha). (a unidade % é gramas de N por 100 gramas de folha)
Dataset_Fundamentos_Grupo_A$N_mass <- Dataset_Fundamentos_Grupo_A$N.percent * 10

# N_area é a quantidade de nitronênio na folha expresso em unidades de área (g cm-2), ao invés de massa
Dataset_Fundamentos_Grupo_A$N_area <- Dataset_Fundamentos_Grupo_A$N_mass * Dataset_Fundamentos_Grupo_A$LMA / 1000

# P_mass é a quantidade de fósforo foliar expresso em unidades de massa (mg de P por grama de folha)
Dataset_Fundamentos_Grupo_A$P_mass <- Dataset_Fundamentos_Grupo_A$P.percent * 10

# P_area é a quantidade de fósforo na folha expresso em unidades de área (g cm-2), ao invés de massa
Dataset_Fundamentos_Grupo_A$P_area <- Dataset_Fundamentos_Grupo_A$P_mass * Dataset_Fundamentos_Grupo_A$LMA / 1000

# Como a base de dados não é balanceada (não tem o mesmo número de amostras de cada indivíduo, ou de cada espécies)
# precisamos tirar primeiro a média de cada réplica medida em um mesmo indivíduo
# e depos, tirar as médias de para cada espécies.



Dados_Medias_Réplicas <- Dataset_Fundamentos_Grupo_A %>%
  group_by(Plot, Region, Vegetation, Species, Family, Ind.Number) %>% 
  summarise_at(vars("Leaf.Area.cm2", "Sample.weight.g", "A.Sat", "gs.at.Asat","Vcmax.Ci.area.25C", "N.percent", "C.percent", "P.percent", "LMA", "N_mass", "N_area", "P_mass", "P_area"), mean)


Dados_Medias_especies <- Dados_Medias_Réplicas %>%
  group_by(Plot, Region, Vegetation, Species, Family) %>% 
  summarise_at(vars("Leaf.Area.cm2", "Sample.weight.g", "A.Sat", "gs.at.Asat","Vcmax.Ci.area.25C", "N.percent", "C.percent", "P.percent", "LMA", "N_mass", "N_area", "P_mass", "P_area"), mean)



boxplot(Leaf.Area.cm2~Vegetation, data = Dados_Medias_especies)
t.test(Leaf.Area.cm2~Vegetation, data = Dados_Medias_especies)

boxplot(Sample.weight.g~Vegetation, data = Dados_Medias_especies)
t.test(Sample.weight.g~Vegetation, data = Dados_Medias_especies)

boxplot(LMA~Vegetation, data = Dados_Medias_especies)
t.test(LMA~Vegetation, data = Dados_Medias_especies)

boxplot(N.percent~Vegetation, data = Dados_Medias_especies)
t.test(N.percent~Vegetation, data = Dados_Medias_especies)

boxplot(C.percent~Vegetation, data = Dados_Medias_especies)
t.test(C.percent~Vegetation, data = Dados_Medias_especies)

boxplot(P.percent~Vegetation, data = Dados_Medias_especies)
t.test(P.percent~Vegetation, data = Dados_Medias_especies)

boxplot(N_mass~Vegetation, data = Dados_Medias_especies)
t.test(N_mass~Vegetation, data = Dados_Medias_especies)

boxplot(P_mass~Vegetation, data = Dados_Medias_especies)
t.test(P_mass~Vegetation, data = Dados_Medias_especies)

boxplot(N_area~Vegetation, data = Dados_Medias_especies)
t.test(N_area~Vegetation, data = Dados_Medias_especies)

boxplot(P_area~Vegetation, data = Dados_Medias_especies)
t.test(P_area~Vegetation, data = Dados_Medias_especies)

# para avaliar diferenças entre Savana e Floresta, também podemos avaliar como diferentes variáveis se relacionam entre si

model.1 <- lm(Dados_Medias_especies$Leaf.Area.cm2[Dados_Medias_especies$Vegetation == "Open"]~Dados_Medias_especies$Sample.weight.g[Dados_Medias_especies$Vegetation == "Open"])
model.2 <- lm(Dados_Medias_especies$Leaf.Area.cm2[Dados_Medias_especies$Vegetation == "Closed"]~Dados_Medias_especies$Sample.weight.g[Dados_Medias_especies$Vegetation == "Closed"])

plot(Leaf.Area.cm2~Sample.weight.g, data = Dados_Medias_especies)
points(Dados_Medias_especies$Leaf.Area.cm2[Dados_Medias_especies$Vegetation == "Open"]~Dados_Medias_especies$Sample.weight.g[Dados_Medias_especies$Vegetation == "Open"], col= 2)
points(Dados_Medias_especies$Leaf.Area.cm2[Dados_Medias_especies$Vegetation == "Closed"]~Dados_Medias_especies$Sample.weight.g[Dados_Medias_especies$Vegetation == "Closed"], col= 4)
abline(model.1, col = 2)
abline(model.2, col = 4)
# para comparar se a elevação da reta que descreve as variações em área foliar baseado na massa foliar, podemos aplicar uma análise de variância (ANOVA)
anova(model.1, model.2)

# para conteúdo de carbono na folha contra LMA:
model.3 <- lm(Dados_Medias_especies$LMA[Dados_Medias_especies$Vegetation == "Open"]~Dados_Medias_especies$C.percent[Dados_Medias_especies$Vegetation == "Open"])
model.4 <- lm(Dados_Medias_especies$LMA[Dados_Medias_especies$Vegetation == "Closed"]~Dados_Medias_especies$C.percent[Dados_Medias_especies$Vegetation == "Closed"])

plot(LMA~C.percent, data = Dados_Medias_especies)
points(Dados_Medias_especies$LMA[Dados_Medias_especies$Vegetation == "Open"]~Dados_Medias_especies$C.percent[Dados_Medias_especies$Vegetation == "Open"], col= 2)
points(Dados_Medias_especies$LMA[Dados_Medias_especies$Vegetation == "Closed"]~Dados_Medias_especies$C.percent[Dados_Medias_especies$Vegetation == "Closed"], col= 4)
abline(model.3, col = 2)
abline(model.4, col = 4)

anova(model.3, model.4)

# para conteúdo de nireogênio na folha contra fósforo:
model.5 <- lm(Dados_Medias_especies$N_mass[Dados_Medias_especies$Vegetation == "Open"]~Dados_Medias_especies$P_mass[Dados_Medias_especies$Vegetation == "Open"])
model.6 <- lm(Dados_Medias_especies$N_mass[Dados_Medias_especies$Vegetation == "Closed"]~Dados_Medias_especies$P_mass[Dados_Medias_especies$Vegetation == "Closed"])

anova(model.5, model.6)


ggscatter(Dados_Medias_especies, x = "P_mass", y = "N_mass",
          add = "reg.line",                         # Add regression line
          conf.int = TRUE,                          # Add confidence interval
          color = "Vegetation", palette = "jco",           # Color by groups "cyl"
          shape = "Vegetation"                             # Change point shape by groups "cyl"
)+
  stat_cor(aes(color = Vegetation), label.x = 3)           # Add correlation coefficient