reg.fat=function(X,Y,grau=1){
 if(grau>3) print("Não é possível realizar")
 if(grau<=3){ 
  X1=tapply(X,X,mean)
  Y1=tapply(Y,X,mean)
  n=length(X1)
  XX=matrix(0,n,grau)
  for(i in 1:grau){
   XX[,i]=X1^i
  }
  model=lm(Y1~XX)
  Y2=predict(model)
  a=min(Y1,Y2)
  b=max(Y1,Y2)
  a=a-0.05*a
  b=b+0.05*b
  co=coef(model)
  X2=seq(min(X1),max(X1),length.out=100)
  XX1=matrix(0,100,grau)
  for(i in 1:grau){
   XX1[,i]=X2^i
  }
  aux=rep(1,100)
  XX1=cbind(aux,XX1)
  Y3=XX1%*%co
  real=cbind(X1,Y1)
  previsto=cbind(X2,Y3)
  if(grau==1) return(list(real=real,previsto=previsto))
  if(grau>1){
   if(grau==2){
    b1=co[2]
    b2=co[3]
    XX=-b1/(2*b2)
    if(2*b2<0) tipo="Máximo"
    if(2*b2>0) tipo="Minimo"
    if(2*b2==0) tipo="Inflexão"
    Tipo=tipo
    resul=cbind(round(XX,5),tipo)
   }
   if(grau==3){
    b1=co[2]
    b2=co[3]
    b3=co[4]
    delta=(2*b2)^2-(4*3*b1*b3)
    if(delta>0){
     X1=(-2*b2+sqrt(delta))/(6*b3)
     X2=(-2*b2-delta)/(6*b3)
     XX=c(X1,X2)
     Tipo=c(0,0)
     for(i in 1:2){
      y=6*b3*XX[i]+2*b2
      if(y<0) tipo="Máximo"
      if(y>0) tipo="Minimo"
      if(y==0) tipo="Inflexão"
      Tipo[i]=tipo
     }
     resul=cbind(round(XX,5),tipo)
    }
    if(delta<0){  
     resul="Não foi possivel determinar pontos criticos"
    }
   }
  }
  return(list(real=real,previsto=previsto,ponto=resul))
 }
}


grafico.regfat2=function(fatorA,fatorB,VR,grau=c(0,1,2),ylab="Y",xlab="X",nomes=NULL,local="topright"){
 fatorA=as.factor(fatorA)
 n=nlevels(fatorA)
 nivel=levels(fatorA)
 data=data.frame(VR,fatorA,fatorB)
 AA=vector("list", n)
 limites=matrix(0,n,2)
 for(i in 1:n){
  data1=subset(data,fatorA==nivel[i])
  X=data1$fatorB
  Y=data1$VR
  g=grau[i]
  if(g>0){
   AA[[i]]=reg.fat(X,Y,grau=g)
  }
  if(g==0){
   X1=tapply(X,X,mean)
   Y1=tapply(Y,X,mean)
   X2=seq(min(X1),max(X1),length.out=100)
   Y2=rep(mean(Y1),100)
   real=cbind(X1,Y1)
   previsto=cbind(X2,Y2)
   AA[[i]]=list(real=real,previsto=previsto)
  }
  previsto=AA[[i]]$previsto
  limites[i,]=c(min(Y,previsto[,2]),max(Y,previsto[,2]))
 }
 ymin=min(limites[,1])-0.05*min(limites[,1])
 ymax=max(limites[,2])+0.05*max(limites[,2])
 plot(X,Y,type="p",ylim=c(ymin,ymax),col="white",xlab=xlab,ylab=ylab,xaxt="n")
 axis(1,X)
 if(is.null(nomes)){
  nomes=nivel
 }
 legend(local,nomes,col=seq(1,n),lty = 1,bty ="n")
 for(i in 1:n){
  points(AA[[i]]$real,col=i)
  lines(AA[[i]]$previsto,col=i)
  g=grau[i]
  if(g>1){
   print(list(Nivel=nivel[i],ponto=AA[[i]]$ponto))
  }

 }
}

  

grafico.reg=function(X,Y,grau=1,xlab=NULL,ylab=NULL){
 reg=reg.fat(X,Y,grau=grau)
 real=reg$real
 previsto=reg$previsto
 a=min(real[,2],previsto[,2])-0.05*min(real[,2],previsto[,2])
 b=max(real[,2],previsto[,2])+0.05*max(real[,2],previsto[,2])
 plot(real,ylim=c(a,b),xlab=xlab,ylab=ylab,xaxt="n")
 lines(previsto)
 axis(1,real[,1])
 if(grau>1) return(reg$ponto)
  
 
}


oneillmathewsDBC=function(resp,Trat,Bloco){
 Bloco=as.factor(Bloco)
 Trat=as.factor(Trat)
 ntrat=length(levels(factor(Trat)))
 nbloc=length(levels(factor(Bloco)))
 data=data.frame(Trat,Bloco,resp)
 data=data[order(Trat),]
 zdados=y=matrix(0,ntrat,nbloc)
 for(i in 1:ntrat) {
  y[i,]=data$resp[((i-1)*nbloc + 1) : (i*nbloc)] 
 }
 trat.mean=apply(y,1,mean)
 bloc.mean=apply(y,2,mean)
 g.mean=mean(resp)
 for(i in 1:ntrat){
  for(j in 1:nbloc){
   zdados[i,j]=abs(y[i,j]- trat.mean[i] - bloc.mean[j] + g.mean) 
  } 
 }
 zdados=as.vector(zdados)
 dadosz=data.frame('z'=zdados,'blocagem'=rep(1:nbloc,each=ntrat),'tratamento'=rep(seq(1:ntrat),nbloc))
 z=dadosz$z
 tratamento=dadosz$tratamento
 blocagem=dadosz$blocagem
 Fc6=summary(aov(z ~ factor(tratamento) + factor(blocagem)))[[1]][1,4]
 rho=c(-1/(ntrat-1), -1/(nbloc-1), 1/((nbloc-1)*(ntrat-1)))
 w0=1-(2/pi)
 w1=(2/pi)*(sqrt(1-rho[1]^2)+rho[1]*asin(rho[1])-1)
 w2=(2/pi)*(sqrt(1-rho[2]^2)+rho[2]*asin(rho[2])-1)
 w3=(2/pi)*(sqrt(1-rho[3]^2)+rho[3]*asin(rho[3])-1)
 m=(w0-w1-w2+w3)/(w0-w1+(nbloc-1)*(w2-w3))
 Fc18=m*Fc6
 pvalor.hvar=(1-pf(Fc18, (ntrat-1), (nbloc-1)*(ntrat-1)))
 output=cbind(Fc18,pvalor.hvar)
 colnames(output)=c("F value","valor-p")
 rownames(output)=""
 cat("\n Teste de homogeneidade de variância ( ONeill&Mathews ).\n")
 return(output) 
}


DQL.rep=function(Quadrado,Linha,Coluna,Tratamento,Resp,quali = TRUE, mcomp = "tukey",sigT = 0.05, sigF = 0.05){
 Quadrado=as.factor(Quadrado)
 Linha=as.factor(Linha)
 Coluna=as.factor(Coluna)
 Tratamento=as.factor(Tratamento)
 model=aov(Resp~Quadrado+Coluna+Linha+Tratamento)
 tab=anova(model)
 print(tab)
 cv <- round(sqrt(tab[5, 3])/mean(Resp) * 100, 2)
 if(tab[4, 5]<sigF) {
  if(quali == TRUE){
   if(mcomp=="tukey"){
    tukey(Resp, Tratamento,tab[5,1],tab[5,2], sigT)
   }
   if(mcomp=="duncan"){
    
   }
   if(mcomp=="lsd"){
    lsd(Resp, Tratamento,tab[5,1],tab[5,2], sigT)
   }
   if(mcomp=="lsdb"){
    lsdb(Resp, Tratamento,tab[5,1],tab[5,2], sigT)
   }
   if(mcomp=="sk") {
    scottknott(Resp, Tratamento,tab[5,1],tab[5,2], sigT)
   }
   if(mcomp=="snk"){
    snk(Resp, Tratamento,tab[5,1],tab[5,2], sigT)
   }
  }
  else{
   reg.poly(Resp, Tratamento,tab[5,1],tab[5,2],tab[4,1],tab[4,2] )
  }
 }
}

   

pressuposicao=function(model,desing="DIC",BLOCO="bloco",LINHA="linha",COLUNA="coluna"){
 e=rstandard(model)
 normal=matrix(0,2,2)
 sw=shapiro.test(e) 
 normal[1,1]=sw$statistic
 normal[1,2]=sw$p.value
 m=mean(e)
 dp=sd(e)
 ks=ks.test(e,"pnorm",sd=dp)
 normal[2,1]=ks$statistic
 normal[2,2]=ks$p.value
 rownames(normal)=c("Shapiro-Wilk","Kolmogorov-Smirnov")
 colnames(normal)=c("Estatística","valor-p")
 set.seed(12345)
 dw=durbinWatsonTest(model,reps=10000)
 if(desing=="DIC"){
  homogeneidade=matrix(0,1,2)
  lev=leveneTest(model)
  homogeneidade[1,1]=lev[1,2]
  homogeneidade[1,2]=lev[1,3]
  rownames( homogeneidade)=c("Levene")
  colnames(homogeneidade)=c("Estatística","valor-p")
  resul=list(Normalidade=normal,Homogeneidade=homogeneidade,Independencia=dw)
 }
 if(desing=="DBC"){
  data=model$model
  vr=data[,1]
  data=data[,-1]
  n=ncol(data)
  nomes=names(data)
  for(i in 1:n){
   if(nomes[i]==BLOCO) j=i
  }
  bloco=data[,j]
  data=data[,-j]
  trat=interaction(data)
  homogeneidade=oneillmathews(vr,trat,bloco)
  rownames( homogeneidade)=c("Oneill Mathews")
  colnames(homogeneidade)=c("Estatística","valor-p")
  data1=data.frame(vr,bloco,trat)
  model1=aov(vr~bloco+trat + Error(bloco/trat), data1)
  aditivo=tukey.1df(model1, data1, error.term= "bloco:trat")
  aditivo1=cbind(aditivo$Tukey.F,aditivo$Tukey.p)
  colnames(aditivo1)=c("Estatística","valor-p")
  rownames( aditivo1)=c("Aditividade de Tukey")
  resul=list(Normalidade=normal,Homogeneidade=homogeneidade,Independencia=dw,aditivo=aditivo1)
 }
 if(desing=="DQL"){
  data=model$model
  vr=data[,1]
  data=data[,-1]
  n=ncol(data)
  nomes=names(data)
  for(i in 1:n){
   if(nomes[i]==LINHA) j=i
   if(nomes[i]==COLUNA) k=i
  }
  linha=data[,j]
  coluna=data[,k]
  data=data[,-c(j,k)]
  trat=interaction(data)
  homogeneidade=matrix(0,1,2)
  lev=leveneTest(vr,trat)
  homogeneidade[1,1]=lev[1,2]
  homogeneidade[1,2]=lev[1,3]
  rownames( homogeneidade)=c("Levene")
  colnames(homogeneidade)=c("Estatística","valor-p")
  data1=data.frame(vr,linha,coluna,trat)
  model1=aov(vr~linha+coluna+trat + Error(linha*coluna), data1)
  aditivo=tukey.1df(model1, data1, error.term= "linha:coluna")
  aditivo1=cbind(aditivo$Tukey.F,aditivo$Tukey.p)
  colnames(aditivo1)=c("Estatística","valor-p")
  rownames( aditivo1)=c("Aditividade de Tukey")
  resul=list(Normalidade=normal,Homogeneidade=homogeneidade,Independencia=dw,aditivo=aditivo1)
 }
 return(resul)
}


oneillmathewsDQL=function(resp,Trat,Linha,Coluna){
 Linha=as.factor(Linha)
 Coluna=as.factor(Coluna)
 Trat=as.factor(Trat)
 n=length(resp)
 LT=levels(factor(Trat))
 LL=levels(factor(Linha))
 LC=levels(factor(Coluna))
 ntrat=length(LT)
 mtrat=matrix(0,n,ntrat)
 for(i in 1:ntrat){
  j=which(Trat==LT[i])
  mtrat[j,i]=1
 }
 mlinha=matrix(0,n,ntrat)
 for(i in 1:ntrat){
  j=which(Linha==LL[i])
  mlinha[j,i]=1
 }
 mcoluna=matrix(0,n,ntrat)
 for(i in 1:ntrat){
  j=which(Coluna==LC[i])
  mcoluna[j,i]=1
 }
 mtrat1=tapply(resp,Trat,mean)
 trat.mean=mtrat1%*%t(mtrat)
 mlinha1=tapply(resp,Linha,mean)
 linha.mean=mlinha1%*%t(mlinha)
 mcoluna1=tapply(resp,Coluna,mean)
 coluna.mean=mcoluna1%*%t(mcoluna)
 g.mean=mean(resp)
 z=abs(resp-trat.mean-linha.mean-coluna.mean+g.mean)
 z=as.vector(z)
 Fc6=summary(aov(z ~ Trat + Linha+Coluna))[[1]][1,4]
 rho=c(-1/(ntrat-1), 2/((ntrat-1)*(ntrat-2)))
 w0=1-(2/pi)
 w1=(2/pi)*(sqrt(1-rho[1]^2)+rho[1]*asin(rho[1])-1)
 w2=(2/pi)*(sqrt(1-rho[2]^2)+rho[2]*asin(rho[2])-1)
  m=(w0-3*w1+2*w2)/(w0+(ntrat-3)*w1-(ntrat-2)*w2)
 Fc18=m*Fc6
 pvalor.hvar=(1-pf(Fc18, (ntrat-1), (ntrat-1)*(ntrat-2)))
 output=cbind(Fc18,pvalor.hvar)
 colnames(output)=c("F value","valor-p")
 rownames(output)=""
 cat("\n Teste de homogeneidade de variância ( ONeill&Mathews ).\n")
 return(output) 
}