# *******************************************************************
# 1. Definição da trilha de dados
# *******************************************************************
setwd("C:\\dados\\")


# *******************************************************************
# 2. Leitura dos dados
# *******************************************************************
X<-read.table("c:\\dados\\plano2.dat",header=F)
# Listar dados
X


# *******************************************************************
# 3.1 Anova em Blocos ao acaso = Yijk = Bk + Gi + Aj + GAij + Eijk
# *******************************************************************
Amb<-as.factor(X[,1])
Trat<-as.factor(X[,2])
Blocos<-as.factor(X[,3])
#Número de fatores (ambientes, tratamentos e blocos)
col_inicial<-4
#rotina ANOVA todas as variáveis
for(i in col_inicial:ncol(X))
{
cat("Variavel = ",i-col_inicial+1)
Modelo<-aov(X[,i]~Blocos+Trat+Amb+Trat*Amb)
print(summary(Modelo))
}


# *******************************************************************
# 3.2 Anova em Blocos ao acaso = Yijk =  Gi + Aj + GAij + Eijk
# *******************************************************************
Amb<-as.factor(X[,1])
Trat<-as.factor(X[,2])
Blocos<-as.factor(X[,3])
#Número de fatores (ambientes, tratamentos e blocos)
col_inicial<-4
#rotina ANOVA todas as variáveis
for(i in col_inicial:ncol(X))
{
cat("Variavel = ",i-col_inicial+1)
Modelo<-aov(X[,i]~Trat+Amb+Trat*Amb)
print(summary(Modelo))
}


# *******************************************************************
# 3.3 Anova em Blocos ao acaso = Yijk = Bk/Aj + Gi + Aj + GAij + Eijk
# *******************************************************************
Amb<-as.factor(X[,1])
Trat<-as.factor(X[,2])
Blocos<-as.factor(X[,3])
#Número de fatores (ambientes, tratamentos e blocos)
col_inicial<-4
#rotina ANOVA todas as variáveis
for(i in col_inicial:ncol(X))
{
cat("Variavel = ",i-col_inicial+1)
Modelo<-aov(X[,i]~Blocos:Amb+Trat+Amb+Trat*Amb)
print(summary(Modelo))
}


# *******************************************************************
# 4a. Teste de Normalidade  (Shapiro-Wilk normality test)
# *******************************************************************
variavel<-4
Amb<-as.factor(X[,1])
Trat<-as.factor(X[,2])
Blocos<-as.factor(X[,3])
Modelo<-aov(X[,variavel]~Blocos+Trat+Amb+Trat*Amb)
summary (Modelo)
resid<-residuals(Modelo)
shapiro.test (resid)


# *******************************************************************
# 5. Análise de Regressão Polinomial
# *******************************************************************
variavelY<-X[,4]
variavelX<-X[,1]
modelo<- variavelY ~ variavelX + I(variavelX^2) + I(variavelX^3)
anova<-aov(modelo, X)
summary (anova)
coef<-anova$coefficients
coef
plot(variavelX,variavelY)
lines(x=variavelX,y=coef[1] + coef[2]*variavelX + coef[3]*variavelX^2 + coef[4]*variavelX^3)


# *******************************************************************
# 6. Superfície de resposta
# *******************************************************************
variavelY<-X[,4]
vX1<-X[,1]
vX2<-X[,2]
modelo<- variavelY ~ vX1 + I(vX1^2)+ vX2 + I(vX2^2) + I(vX1*vX2) + I(vX1*vX2^2) +  I(vX1^2*vX2) + I(vX1^2*vX2^2)
anova<-aov(modelo, X)
summary (anova)
coef<-anova$coefficients
coef
equacao<-function(vX1,vX2){coef[1]+ coef[2]*vX1 + coef[3]*I(vX1^2)+ coef[4]*vX2 + coef[5]*I(vX2^2) + coef[6]*I(vX1*vX2) + coef[7]*I(vX1*vX2^2) +  coef[8]*I(vX1^2*vX2) + coef[9]*I(vX1^2*vX2^2)}
valX1<-seq(min(vX1), max(vX1), .5)
valX2<-seq(min(vX2), max(vX2), .5)
superficie<-outer(valX1, valX2, equacao)
persp(superficie, theta = 310, phi=30, exmpand = 0.5, col = 4, shade = 0.5, ticktype = "simple")
#superficie: valores interpolados
#theta: ângulo horizontal de exibição
#phi: ângulo vertical de exibição
#expand: expandir o gráfico
#col: cor
#shade:  sombreamento
#ticktype = simple  eixo sem escala