#---------------------------------------------------------
# Exemplo obtido do caderno de mestrado da disciplina de -
# Estatística Experimental do Prof. Décio Barbin.        -
#------------------------------------------------------------------
# Um pesquisador pretende estudar ou comparar quatro variedades   -
# de pêssego, quanto ao enraizamento de estacas. O experimento    -
# será conduzido em um viveiro, sendo cada parcela constituída    -
# por 20 estacas. Optou-se por 5 repetiçoes, logo, 4 tratamentos  -
# com 5 repetições, tem-se 20 parcelas de 20 estacas cada.        -
#------------------------------------------------------------------
# Passado o período experimental, o pesquisador contou as estacas -
# enraizadas e os dados são apresentados.                         -
#------------------------------------------------------------------
rm(list=ls())

(croqui = expand.grid(rep=1:5, trat=c('A','B','C','D')))

(resp = c( 2,  2,  1,  1,  0, 
           1,  0,  0,  1,  1,
           12, 10, 14, 17, 11,
           7,  9, 15,  8, 10))

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

names(dados)
str(dados)
attach(dados)

#-------------------------
# Estatística Descritiva -
#-------------------------
#
#---------------------------------------------------------
# Estatísticas das estacas, sem levar em conta os grupos -
#---------------------------------------------------------
	(n.g = length(resp))
  (media.g = mean(resp))
	(var.g = var(resp))
	(desvio.g = sd(resp))
	(cv.g = desvio.g/media.g*100)
  (ep.g = desvio.g / sqrt(n.g))

#-------------------------------------------------------
# Estatísticas das estacas, levando em conta os grupos -
#-------------------------------------------------------
  (n =  tapply(resp, trat, length))
	(medias = tapply(resp, trat, mean))
	(vari = tapply(resp, trat, var))
	(desvios = tapply(resp, trat, sd))
  (cvs = desvios / medias * 100)
  (eps = desvios / sqrt(n))

	(quantis = tapply(resp, trat, quantile))

	(SQTrat = 6*sum((medias-mediag)^2))

	par(mai=c(1,1,.2,.2))
	plot.default(trat, resp, las=1, xlab='Variedades de Pêssego', 
      	       ylab='', col='blue', bty='l')
	mtext('Número de estacas ', side=2, line=3)
	points(medias, pch="M", col='red', cex=1.2)

#--------------------
# Gráfico de Caixas -
#--------------------
require(car)

par(mai=c(1,1,.2,.2))
Boxplot(resp ~ trat, names=c("A", "B", "C", "D"), ylab="Número de Estacas", xlab="Variedades de Pêssego",
	        las=1, col='LightYellow')

points(medias, pch="+", col=2, cex=1.5)

#------------------------------------------
# stripchart é uma alternativa ao boxplot -
# para amostras de tamanho pequeno.       -
#------------------------------------------
stripchart(resp ~ trat, method="jitter", vertical = TRUE, las=1, pch=1, 
           col='blue', group.names=c('A','B','C','D'),
           ylab="Número de Estacas", xlab="Variedades de Pêssego")

points(medias, pch = 19, cex = 1.3, col='red')

arrows(1:4, medias + eps, 1:4, medias - eps, angle = 90,
       code = 3, length = 0.1, col='red')


#-------------------------------------
# Procedemos à análise de variância: - 
#-------------------------------------
mod = aov(resp ~ trat)
anova(mod)
summary(mod)

#-------------------------------
# Verificação dos pressupostos -
#-------------------------------
#
# Normalidade dos erros:
shapiro.test(mod$res)

# Teste de Lilliefors para normdalidade dos erros:
require(nortest)
lillie.test(mod$res)

# Teste de Levene:
leveneTest(mod$res ~ trat, center=mean)

# Variâncias homogêneas:
bartlett.test(mod$res ~ trat)

# Ajuste do modelo usando o Q-Q Plot
fit.model = mod
source('C:/Livros/Modelos Lineares Generalizados/Gilberto de Paula/Arquivos R/Normal/envel_norm.R')
# source('http://www.ime.usp.br/~giapaula/envel_norm.R')


#--------------------------
# Transformação dos dados -
#--------------------------
require(MASS)

BOX = boxcox(resp + .00001 ~ trat, data=dados, lam=seq(0.3, .7, 1/10))
(lambda = BOX$x[which.max(BOX$y)])


#-------------------------------------------------
# Análise dos dados transformados: raiz quadrada -
#-------------------------------------------------
mod2 = aov(sqrt(resp) ~ trat)
anova(mod2)

# Do caderno
mod3 = aov(sqrt(resp + 0.5) ~ trat)
anova(mod3)

#-------------------------------
# Verificação dos pressupostos -
#-------------------------------
#
# Normalidade dos erros:
shapiro.test(mod2$res)

# Teste de Lilliefors para normdalidade dos erros:
require(nortest)
lillie.test(mod2$res)

# Variâncias homogêneas:
bartlett.test(sqrt(resp) ~ trat)

# Teste de Levene:
require(car)
leveneTest(sqrt(resp) ~ trat, center=mean)

#-----------------
# Teste de Tukey -
#-----------------
require(agricolae)

tukey = HSD.test(mod2, 'trat', console=T)

bar.group(tukey$groups, las=1, col='lightyellow', xlab='Tratamentos',
          ylab='Número de estacas enraizadas', ylim=c(0,4))
abline(h=0)