#--------------------------------------------------------------------
# Com a finalidade de estudar os efeitos da administração de raízes -
# e tubérculos, como suplementação de inverno na alimentação de     -
# vacas em lactação, considerou-se um experimento em Blocos         -
# casualizados com 4 tipos de suplementos (tratamentos) e 5 raças   -
# (Blocos).                                                         -
# A variável resposta é a produção média diárias de leite (kg).     - 
#--------------------------------------------------------------------
rm(list=ls())

# Gerando o croqui do experimento no delineamento em blocos
require(agricolae)

Blocos = factor(c('Gir','Holandesa','Jersey','Nelore','Guzerá'))
Trat = factor(c('S','M','A','BD'), level=c('S','M','A','BD'))

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

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

#----------------------
# Após obter os dados -
#----------------------
Prod = c( 6.4, 12.0, 11.2, 10.9,
         11.6, 10.9,  6.2, 11.6,
         11.5, 11.4, 10.9,  6.2, 
         10.4,  7.1, 12.1, 11.1, 
         12.4,  6.6, 10.1, 11.8)

dados$resp = Prod

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

#---------------------------
# Estatísticas Descritivas -
#---------------------------
(média.g = with(dados, mean(resp)))
(variância.g = with(dados, var(resp)))

(médias.t = with(dados, tapply(resp, Trat, mean)))
(médias.b = with(dados, tapply(resp, Blocos, mean)))

(var.t = with(dados, tapply(resp, Trat, var)))
(var.b = with(dados, tapply(resp, Blocos, var)))

#-----------
# Gráficos -
#-----------
par(mai=c(1.2, 1.2, .5, .5))

with(dados, interaction.plot(Trat, Blocos, resp, las=1, xlab='Tratamentos', ylab='', col=1:5, bty='l'))
mtext('Produção média (kg)', side=2, line=3.5)

with(dados, interaction.plot(Blocos, Trat, resp, las=1, xlab='Blocos', ylab='', col=1:5, bty='l'))
mtext('Produção média (kg)', side=2, line=3.5)

with(dados, plot.default(Trat, resp, las=1, xlab='Tratamentos', ylab='', col='blue', bty='l', pch=19))
mtext('Produção média (kg)', side=2, line=3.5)
points(médias.t, pch="*", col='red', cex=2.5)
abline(h=média.g, col='red', lwd=2)

with(dados, plot.default(Blocos, resp, las=1, xlab='Blocos', pch=19, ylab='', col='blue', bty='l'))
mtext('Produção média (kg)', side=2, line=3.5)
points(médias.b, pch="*", col='red', cex=2.5)
abline(h=média.g, col='red', lwd=2)

with(dados, boxplot(resp ~ Trat, names=c("Sem Suplemento", "Mandioca", "Araruta", "Batata Doce"),
                    ylab="Produção média (kg)", xlab="Tratamentos", las=1, col='LightYellow'))
points(médias.t, pch="+", col=2, cex=1.5)

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

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

#---------------------------------------
# Pressupostos para validade do modelo -
#---------------------------------------
#
# Normalidade dos erros
plot(mod, which=c(2:2), pch=19, col='red', las=1)

require(hnp)
hnp(mod, las=1, pch=19)

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

require(nortest)
lillie.test(mod$res)     # Teste de Kolmogorov-Smirnov

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

require(car)
with(dados, leveneTest(resp, Trat, center=mean))    # Teste de Levene

# Aditividade de dos efeitos dos fatores
require(asbio)
with(dados, tukey.add.test(resp, Trat, Blocos))

# ou, usando o pacote agricolae
require(agricolae)

(gl = anova(mod)$Df[3])
(QMRes = anova(mod)$Mean[3])

nonadditivity(y=resp, factor1=Blocos, factor2=Trat, gl, QMRes)
 
#----------------------------------
# Testes de COMPARAÇÕES MÚLTIPLAS -
#----------------------------------
#
#-----------------
# Teste de Tukey -
#-----------------
(tukey = TukeyHSD(mod, 'Trat', ord=T))
plot(tukey, las=1, col='blue')

# Usando o pacote agricolae
library(agricolae)
teste.HSD = HSD.test(mod, 'Trat', main='Efeito dos tratamentos na Produção leite', console=T)
names(teste.HSD)

bar.group(teste.HSD$groups, density=10, las=1, angle=45, col='red', border='blue', ylim=c(0,15),
          main='Teste de Tukey', xlab='Tratamentos', ylab='Médias de Produção (kg)')
abline(h=0, col='black')

#------------------
# Teste de Duncan -
#------------------
teste.duncan = duncan.test(mod, 'Trat', main='Efeito dos tratamentos na Produção leite', console=T)

bar.group(teste.duncan$groups, density=10, las=1, angle=45, col='red', border='blue', ylim=c(0,15),
          main='Teste de Duncan', xlab='Tratamentos', ylab='Médias de Produção (kg)', space=.7)
abline(h=0, col='black')

#---------------------------
# Teste de Tukey e Dunnett -
#---------------------------
require(multcomp)
MTukey = glht(mod, linfct = mcp(Trat="Tukey"))
summary(MTukey)

summary((td = glht(mod, linfct=mcp(Trat="Dunnett"))))
confint(td, level = 0.95)

#---------------------------------
# Análise usando o pacote ExpDes -
#---------------------------------
require(ExpDes.pt)

with(dados, dbc(Trat, Blocos, resp, quali=TRUE, mcomp="tukey"))
with(dados, dbc(Trat, Blocos, resp, quali=TRUE, mcomp="duncan"))
with(dados, dbc(Trat, Blocos, resp, quali=TRUE, mcomp="snk"))
with(dados, dbc(Trat, Blocos, resp, quali=TRUE, mcomp="sk"))
with(dados, dbc(Trat, Blocos, resp, quali=TRUE, mcomp="lsd"))
with(dados, dbc(Trat, Blocos, resp, quali=TRUE, mcomp="lsdb"))

#----------------------------
# 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