library(epicalc)
library(nortest)
library(multcomp)
library(multcompView)
library(Rcmdr)

EXP <- read.csv("INGESTION.STAT.csv",header=TRUE, sep=",")

#Lendo o arquivo .csv originado pela primeira função ();

#   EXP: GONY1, GONY2, GONY3 (Experimentos com Gonyaulax sp.);
#   GYM1, GYM2, GYM3 (Experimentos com Gymnodinium sp.);
#   ISO1, ISO2, ISO3 (Experimentos com Isochrysis galbana);
#   FRA1, FRA2, FRA3 (Experimentos com Alexandrium cf. fraterculus);
#   TAM1, TAM2, TAM3 (Experimentos com Alexandrium tamiyawanichii);
#   RHO1, RHO2, RHO3 (Experimentos com Rhodomonas cf. salina);
#   MIXA1, MIXA2, MIXA3 (Experimentos de misturaA: PRO/TAM);
#   MIXB1, MIXB2, MIXB3 (Experimentos de misturaB: FRA/PRO/TAM).

#############################################################################################################################################################################
#FUNÇÃO STAT#################################################################################################################################################################

EXP #Conferindo;
TAT <- EXP$C #Transforamndo a primeira coluna do arquivo .txt em um objeto numérico (tratamentos);
TAT.F <- factor(EXP$C) #Transformando a primeira coluna do arquivo .txt em um fator (tratamentos);
IR <- EXP$IR #Atribuindo um nome à segunda coluna do .txt (Ingestion Rate - IR);

TAT.F
IR

#1ºEXPLORANDO OS DADOS GRAFICAMENTE:

X11()
par(mfrow = c(1,2))
plot(TAT,IR) #Dispersão (atenção para TAT na forma numérica);
dotplot(IR, by = TAT.F, main = "", ylab = "TAT") #Posição dos pontos;

#2ºEXECUTANDO UMA REGRESSÃO ENTRE OS DADOS
REG <- aov(IR ~ TAT.F) #Executa uma regressão linear entre TAT e IR - comando similar ao
                       #comando "lm" porém adequado para a realização da ANOVA;
                       #Atribuição do resultado do "aov" ao um objeto ("reg");
REG #Conferindo;

#2ºEXPLORANDO A REGRESSÃO ENTRE OS DADOS GRAFICAMENTE:

#Análise exploratória da regressão: criação de um conjunto de dados diagnósticos para observar o comportamento dos resíduos;
RES <- rstandard(REG) #Atribuindo os resíduos padronizados de REG a um objeto (RES);
FV <- fitted.values(REG) #Atribuindo os valores ajustados de REG a um objeto (FV);
TAT.ORD <- seq(1:15) #Criação de um vetor para ordenar os tratamentos;
TAT.ORD #Conferindo;
summary(RES)
sd(RES)
X11()
par(mfrow = c(3,2))
dotplot(RES, main = "", ylab = "Frequência") #Indicação de independência: resíduos distribuídos em torno do valor "0";
boxplot(RES, main = "", ylab = "Resíduos")   #Distribuição dos resíduos;
plot(TAT.ORD, RES, main = "", xlab = "Ordem", ylab = "Resíduos") #Indicação de independência: resíduos distribuídos em torno do valor "0";
abline(h = 0)
plot(FV, RES, main = "", xlab = "Valores Ajustados", ylab = "Resíduos") #Indicação de independência: distribuição dos resíduos dentro de uma mesma faixa (largura);
abline(h = 0)
qqnorm(RES, xlab = "Quantis Teóricos", ylab = "Quantis Amostrais") #Indicação de normalidade: resíduos semelhantes aos valores ajustados;
qqline(RES, xlab = "Quantis Teóricos", ylab = "Quantis Amostrais")

#3ºVERIFICAÇÃO DOS PRESSUPOSTOS DA ANOVA:

AD <- ad.test(RES) #Teste para verificação da normalidade dos dados;
BARTLETT <- bartlett.test(RES,TAT) #Teste para verificação da homocedasticidade;

AD #Conferindo;
BARTLETT #Conferindo

#4ºSE OS PRESSUPOSTOS SÃO VERDADEIROS, EXECUÇÃO E INTERPRETAÇÃO DA ANOVA:

if(AD$p.value & BARTLETT$p.value >= 0.05){
  
  ANOVA <- anova(REG) #Excutando uma análise de variância de REG;
                      if(ANOVA$p.value >= 0.05){#Caso a ANOVA indique diferencas (Hipótese nula falsa): aplicação da correção de Tukey;
                                                #Aplicação da correção de Tukey e atribuição do resultado do teste a um obejto ("THSD");
                         THSD <- TukeyHSD(REG) 
                         X11()
                         par(mfrow = c(1,2))
                         plot(THSD) #Verificação da diferença;
                         d <- data.frame(resp = IR, TAT = TAT)
                         MULTICOMPBOXPLOT<- multcompBoxplot(resp ~ TAT, data = d, compFn = "TukeyHSD") #Verificação gráfica mais rápida da localização da diferença (Pacote:)
                         
                         tiff("MULTICOMPBOXPLOT") #Saída1: Gráfico para verificação das diferenças;
                         RESULT <- data.frame(THSD$p.value)
                         RESULT #Conferindo;
                         #Saída1: THSD para descrição das diferenças;
                      }

                      else{
                         RESULT <- data.frame(ANOVA$p.value)
                         RESULT #Conferindo;
                         #Saída2: ANOVA (sem diferença);
                      }

}

#5ºSE OS PRESSUPOSTOS NÃO SÃO VERDADEIROS, EXECUÇÃO DO TESTE DE KRUSKAL-WALLIS:

KRUSKAL <- kruskal.test(IR, TAT.F)
KRUSKAL$p.value #Conferindo o valor de p;
#Validação do Teste de Kruskal-Wallis;
LEVENE <- leveneTest(IR,TAT.F) #Verirficação da homocedasticidade (teste menos robusto por é executado
                               #com base nas mediandas).
LEVENE$Pr #Conferindo o valor de p;
RESULT <- data.frame(KRUSKAL$p.value,LEVENE$Pr)
RESULT #Conferindo;


write.table(RESULT, "RESULT.csv", sep="\t") #Saída3: KRUSKALL para descrição das diferenças.

#############################################################################################################################################################################
#############################################################################################################################################################################