#===============================================================================
#. Título: Fundamentos da Análise de Variância
#. Autor : José Cláudio Faria/UESC/DCET
#.         Ivan Bezerra Allaman/UESC/DCET
# Data  : 2025-10-09 12:08:18
# Objetivos:
#===============================================================================
# a) Permitir a modelagem da ANOVA com número variado de grupos e repetições
# b) Exibir resultados numéricos e gráficos
# c) Comparar os resultados por 2 métodos diferentes (Tukey e ScottKnott)
#    de discriminação dos grupos
#===============================================================================
# Pacotes necessários:
# - plotrix
# - TukeyC
# - ScottKnott
#===============================================================================
#.. ___Ini opções___
#... Médias dos grupos
M <- c(60,
#       60,
#       60,
#       73,
       60)

#... Número de repetições de cada grupo (tratamento)
r <- 6

#... Desvio padrão de cada grupo (tratamento)
s <- 6

#... Probabilidade do erro tipo I na inferência (ANOVA)
alfa <- 5/100
#.. ___Fim opções___

# Opção de cores
o_col <- rainbow(length(M))

#.. Funções
plotar_info_Ft <- function()
{  
  polygon(x=c(Ft,
              seq(Ft,
                  20,
                  length=1e3)),
          y=c(0,
              df(seq(Ft,
                     20,
                     length=1e3),
                 length(M) - 1,
                 length(M) * (r-1))),
          col='red')

  segments(x0=Ft,
           x1=Ft,
           y0=0,
           y1=1,
           lty=3,
           col='red')

  text(x=Ft,
       y=0.4,
       pos=4,
       offset=0,
       bquote(alpha == .(alfa)),
       cex=1.4,
       col='red')

  points(x=Ft,
         y=-.06,
         pch=17,
         cex=1.5,
         col='red')
}

plotar_info_Fc <- function()
{  
  ifelse(Fc == 0,
         Fc <- 0.001,
         Fc)
  polygon(x=c(Fc,
              seq(Fc,
                  20,
                  length=1e3)),
          y=c(0,
              df(seq(Fc,
                     20,
                     length=1e3),
                 length(M) - 1,
                 length(M) * (r-1))),
          col=gray(.8))

  if (Fc <= 15)
  {
    segments(x0=Fc,
             x1=Fc,
             y0=0,
             y1=.15,
             lty=3)
  }

  if (Fc <= 15)
  {
    text(x=Fc,
         y=0.12,
         pos=4,
         offset=0,
         bquote(P == .(tab[1, 5])),
         cex=1.4)
  } else {
    text(x=15,
         y=0.12,
         pos=4,
         offset=0,
         bquote(P == .(tab[1, 5])),
         cex=1.4)
  }    
}

#.. Início da modelagem
# Gerando grupos com r repetições com distribuição normal em um data.frame
dad <- data.frame(tr=gl(length(M),
                        r,
                        labels=LETTERS[1:length(M)]),
                  y=rnorm(length(M) * r,
                          rep(M, 
                              each=r),
                          sd=s))

#.. Análise de variância
av <- aov(y ~ tr,
          data=dad)

sav <- summary(av)

tab <- as.data.frame(round(sav[[1]][, 1:5],
                           2))

cv <- 100 * sqrt(tab[2, 3]) / mean(dad$y)

colnames(tab)=c('Gl',
                'SQ',
                'QM',
                'Fcal',
                'P(>Fcal)')

rownames(tab)=c('Grupos',
                'Resíduo')

#.. Calculando médias e desvios dos grupos
medias <- with(dad, 
               tapply(y, 
                      tr, 
                      mean))

desvios <- with(dad, 
                tapply(y, 
                       tr, 
                       sd))

# O ordenamento irá facilitar a padronização de cores nos plots
ord <- order(medias)

medias <- sort(medias)

desvios <- desvios[ord]

#.. Valores F
Ft <- qf(alfa,
         length(M) - 1,
         length(M) * (r-1),
         lower.tail=FALSE)

Fc <- tab[1,
          4]

# Verificar dispositivo gráfico
if (is.null(dev.list()))
  X11(w=8,
      h=8)

# Layout do dispositivo gráfico
layout(matrix(1:6,
              ncol=2,
              byrow=T),
       heights=c(3,
                 4,
                 4))

# Plot da tabela
plot(0.5, 
     type = 'n', 
     axes=F, 
     xlab='', 
     ylab='')

# Tabela da ANOVA no plot
suppressMessages(require(plotrix))
addtable2plot(.5,
              .4,
              tab,
              display.rownames=TRUE,
              bg=gray(.8),
              cex=1.8)

#mtext('ANOVA')

mtext(paste('cv = ',
            round(cv, 
                  2),
            '%',
            sep=''),
      side=1)

#.. Distribuição F e decisão da inferência
curve(df(x,
         length(M) - 1,
         length(M) * (r-1)),
      from=0,
      to=20,
      n=1e3,
      xlab='F',
      ylab='f(F)',
      ylim=c(-.09, 
             1),
      col=gray(.4))

text(x=Ft,
     y=rep(.8, 2),
     lab=c('RAHo', 
           'RRHo'),
     pos=c(2, 
           4),
     offset=.5,
     cex=1.4,
     col=c('darkgreen',
           'red'))

#.. Decisão da ANOVA
if(Fc < Ft){
  points(Fc,
         y=-.08,
         pch=17,
         cex=1.7,
         col=gray(.5))

  mtext(text = 'Não rejeita Ho',
        cex=1.5,
        col='darkgreen')
} else {
  mtext('Rejeita Ho',
        cex=1.5,
        col='red')

  if(Fc <= 15){
    points(Fc,
           y=-.08,
           pch=17,
           cex=1.7,
           col=gray(.5))
  } else {
    text(15,
         -.08,
         pos=4,
         offset=0,
         expression(Delta > 15),
         cex=1.3,
         col=gray(.5)) 
  }
}

#.. Plotando áreas sob a distribuição F
if(Fc >= Ft){
  plotar_info_Ft()
  plotar_info_Fc()
} else {
  plotar_info_Fc()
  plotar_info_Ft()
}

legend('topright',
       pch=17,
       col=c('red',
             gray(.5)),
       legend=c('Ftab', 
                'Fcal'),
       cex=1.5)

means_dec <- with(dad, 
                  reorder(tr, 
                          -y, 
                          mean))
boxplot(y ~ means_dec, 
        data=dad,
        xlab='Grupos', 
        ylab='Y',
        col=o_col)

#.. Plotando as densidades dos grupos
plot(1,
     type='n',
     ylab='Densidade',
     xlab='Médias dos grupos',
     xlim=c( .3 * min(medias),
            1.7 * max(medias)),
     ylim=c(-.005,
            max(c(dnorm(medias[1:length(M)], 
                        medias[1:length(M)], 
                        desvios[1:length(M)])))))

for(i in 1:length(M)){
  curve(dnorm(x,
              medias[i],
              desvios[i]),
        .2 * min(medias[i]),
        1.8 * max(medias[i]),
        n=1e3,
        add=TRUE,
        col=rev(o_col)[i],
        lw=1.5)
}

for(i in 1:length(M)){
  segments(medias[i],
           max(dnorm(medias[i], 
                     medias[i], 
                     desvios[i])),
           medias[i],
           0,
           col=rev(o_col)[i],
           lty=3)
}

for(i in 1:length(M)){
  text(medias[i],
       -.005,
       substitute(mu[i], 
                  list(i = LETTERS[ord][i])),
       col=rev(o_col)[i])
}

legend('topright',
       rev(LETTERS[ord]),
       lty=rep(1, 
               length(M)),
       lw=2,
       col=o_col)

#. Tukey
suppressMessages(require(TukeyC))
tk <- TukeyC(av,
             sig.level=alfa)

plot(tk,
     xlab='Grupos',
     ylab='Y',
     col=o_col,
     disp='cip',
     title=paste('Tukey',
                 ' (',
                 100*alfa,
                 '%',
                 ')',
                 sep=''),
     cex=1.2)

#.. ScottKnott
suppressMessages(require(ScottKnott))
sk <- SK(av,
         sig.level=alfa)

plot(sk,
     xlab='Grupos',
     ylab='Y',
     col=o_col,
     disp='cip',
     title=paste('Scott & Knott',
                 ' (',
                 100*alfa,
                 '%',
                 ')',
                 sep=''),
     cex=1.2)