#==============================================================================
#  Um Fator Aleatório - Teste F - ANOVA
#==============================================================================
#                 Companhia Textil
#  Descrição do Estudo: Resistência (Y) Teares - Amostra de quatro
#                       X: 1,2,3,4
#                        Livro do Montgomery, 1997, ??
#
#===============================================================================
 rm(list=ls(all=TRUE))
#
#############################################################
# instalando alguns pacotes necessários para a análise:
#
#install.packages("multcomp", dependencies=TRUE)
#install.packages("multcompView", dependencies=TRUE)
#install.packages("agricolae", dependencies=TRUE)
#install.packages("nortest", dependencies=TRUE)
#install.packages("car", dependencies=TRUE)
#
# carrecando os pacotes instalados:
#
 require(multcomp)
 require(multcompView)
 require(agricolae)
 require(nortest)
 require(car)
#
#############################################################
#
# lendo o banco de dados:
 dados=read.table("dados3.txt", header=TRUE)
 dados
#
#verificando a classe de cada variável do nosso objeto dados:
 sapply(dados,class)
#
# obtendo as médias amostrais associadas a cada nível do fator de interesse:
 medias <- tapply(dados$resp, dados$tear, mean)
 medias
#
# explorando graficamente:
 plot(resp~tear,data=dados,ylab="resistência")
#
# Obtendo a tabela ANOVA:
 modelo=aov(resp~tear, data=dados)
 summary(modelo)
#
# valor tabelado da distribuição F:
 qf(0.05, 3, 20, lower.tail = FALSE)
#
 class(modelo)
 names(modelo)
#
 n=4
 SQTotal <- sum((dados$resp-mean(dados$resp))^2 )
 SQRes <-  sum(modelo$residuals^2)
 SQTrat <- SQTotal -  SQRes
 QMRes <- SQRes/modelo$df 
 QMTrat <- SQTrat/(modelo$rank-1)
#
 sigma2_tau=(QMTrat-QMRes)/n
 sigma2_tau
#
##########################################################################
# Análise de resíduos
#
 tear.hat=modelo$coeff
 tear.hat
#
 y.hat=modelo$fitted.values
 y.hat
#
 residuo=modelo$residuals
#
############################################
#
# Testes de homogeneidade da variância
 bartlett.test(dados$resp~dados$tear)
 leveneTest(dados$resp~dados$tear)
#
# Testes de normalidade:
 shapiro.test(residuo)
 ad.test(residuo)
#
 par(mfrow=c(2,2))
 plot(modelo)
#
 qqnorm(residuo)
 qqline(residuo, col = 2)
# Alternativamente:
 win.graph()
 par(mfrow=c(1,2))
 qqnorm(residuo)
 qqline(residuo, col = 2)
 hist(residuo)
#
############################################
#
 win.graph()
 plot(y.hat,residuo)
#lines(rep(0,length(residuo)), col="red")
#lines(rep(-1.96,length(residuo)), col="red", lty=2)
#lines(rep(1.96,length(residuo)), col="red", lty=2)
 summary(residuo)
#
