#----------------------------------------------------
#' Experimento simulado de alimentação de poedeiras -
#'                                                  -
#' A formação dos blocos foi:                       -
#' Bloco 1 - as melhores poedeiras                  -
#' Bloco 2 - as segundas melhores poedeiras e assim -
#' por diante.                                      -
#' Número de aves por parcela: 4 a 8                -
#' Os tratamentos são os diferentes tipos de ração  -
#' A variável resposta é o número médio de ovos.    -
#'--------------------------------------------------- 
rm(list=ls())

(croqui = expand.grid(blocos=c('I','II','III','IV'), Trat=c('A','B','C','D','E')))

ovos = c(202.5, 200.4, 180.9, 190.3, 
         220.3, 215.4, 219.6, 210.5,
         210.7, 205.6, 200.4, 190.8,
         230.4, 225.6, 215.7, 220.1,
         200.0, 194.1, 180.5, 190.0)

(dados = data.frame(croqui, resp=ovos))


#---------------------------------------------
# Outra forma, usando o croqui gerado pelo R -
#---------------------------------------------
#
# Gerando o croqui do experimento no delineamento em blocos
require(agricolae)

Trat = factor(c('A','B','C','D','E'))

(croqui = design.rcbd(Trat, r=4, seed=105, kinds="Super-Duper"))

dados = croqui$book
colnames(dados) = c('Parcelas', 'Blocos', 'Trat')
dados

#----------------------
# Após obter os dados -
#----------------------
Ovos = c(220.3, 200.0, 210.7, 230.4, 202.5, 
         205.6, 200.4, 225.6, 194.1, 215.4, 
         180.5, 180.9, 215.7, 200.4, 219.6, 
         190.0, 210.5, 220.1, 190.8,  190.3)

dados$resp = Ovos

head(dados)
tail(dados)
str(dados)

#---------------------------
# Estatísticas descritivas -
#---------------------------
#
# Geral
(média = with(dados, mean(resp)))
(variância = with(dados, var(resp)))
(desvio = with(dados, sd(resp)))
(CV = desvio / média * 100)

# Por Tratamento
(médias = with(dados, tapply(resp, Trat, mean)))
(variâncias = with(dados, tapply(resp, Trat, var)))
(desvios = with(dados, tapply(resp, Trat, sd)))
(cv = desvios / médias * 100)

# Obs.: Caso se queira, repetir os comandos anteriores para blocos.

#------------------
# Análise Gráfica -
#------------------
par(mai=c(1.2, 1.2, .5, .5))

with(dados, interaction.plot(Blocos, Trat, resp, las=1, xlab='Blocos', ylab='', col=1:5, bty='l'))
mtext('Número médio de ovos', side=2, line=3.5)

with(dados, interaction.plot(Trat, Blocos, resp, las=1, xlab='Tratamentos', ylab='', col=1:4, bty='l'))
mtext('Número médio de ovos', side=2, line=3.5)

with(dados, plot.default(Trat, resp, las=1, xlab='Tratamentos', ylab='', col='blue', bty='l'))
mtext('Número médio de ovos', side=2, line=3.5)
points(médias, pch="*", col='red', cex=1.5)

with(dados, plot.default(Blocos, resp, las=1, xlab='Blocos', ylab='', col='blue', bty='l'))
mtext('Número médio de ovos', side=2, line=3.5)
with(dados, points(tapply(resp, Blocos, mean), pch=19, col='red', cex=1.3))

require(car)
with(dados, Boxplot(resp ~ Trat, ylab="Número médio de ovos", xlab="Tratamentos", las=1, col='LightYellow'))
points(médias, pch="+", col=2, cex=1.5)

#-------------------------------------------------
# Visualização de todas as médias do experimento -
#           em relação à média geral             -
#-------------------------------------------------
plot.design(dados[ , -1], xlab='', ylab="", las=1, bty='l', col='blue')
mtext('Número médio de ovos', side=2, line=3.5)
mtext('Fatores', side=1, line=2)

#-----------------------
# Análise de variância -
#-----------------------
mod = with(dados, aov(resp ~ Trat + Blocos))
anova(mod)
summary(mod)

#--------------------------------------------------------
# Obs.: a significância de blocos indica que a formação -
# dos mesmos foi uma prática vantajosa. A ausência de   -
# blocos, neste caso, aumentaria a variância residual e,-
# em consequência, reduziria a precisão do experimento  -
#--------------------------------------------------------

#---------------
# Pressupostos -
#---------------
#
# Normalidade dos erros
qqnorm(mod$res, xlab="Quantis Teóricos", ylab="Residuos", main=NULL, las=1, bty='l', pch=19, col='red')
qqline(mod$res)

# De forma direta:
par(mfrow=c(2,2))
plot(mod)

# Teste de Shapiro-Wilk
shapiro.test(mod$res)
layout(1)

# Homogeneidade de variâncias
with(dados, bartlett.test(mod$res ~ Trat))

# Teste de aditividade de Tukey
require(asbio)
with(dados, tukey.add.test(resp, Blocos, Trat))


#-------------------------------
# Teste de comparação múltipla -
#-------------------------------
#
# Tukey
(tukey = TukeyHSD(mod, 'Trat', ord=T))
plot(tukey, las=1, col='blue')

# Usando o pacote agricolae
require(agricolae)
HSD.test(mod, 'Trat', console=T)

# Teste de Duncan
duncan.test(mod, 'Trat', console=T)

# Teste de Dunnett
# Supondo que o Trat. A seja a testemunha
require(asbio)
with(Dados, pairw.anova(y=resp, x=Trat, control="A", method="dunnett"))


#----------------------------
# Usando o pacote easyanova -
#----------------------------
#
# É necessário deixar apenas as colunas de
# Tratamentos, Blocos e Resposta, nesta ordem
(dados.t = dados[ , c(3,2,4)])

require(easyanova)

mod.1 = ea1(dados.t, design=2, plot=2)       # plot=1,2,3
names(mod.1)
mod.1