#===============================================================================
#. Name   : gexp_quali_psd
#. Author : J.C.Faria
# Date   : 06/10/2022 01:33:35
# Aim    : Gerar/simular e analisar experimentos qualitativos em
#          parcelas subdivididas nos delineamentos experimentais básicos:
#          - Del. Inteiramente ao Acaso - DIC
#          - Del. Blocos Casualizados - DBC
#          - Del. Quadrado Latino - DQL
# Mail   : <<< joseclaudio.faria@gmail.com >>>
#===============================================================================

#.. Instalação do pacote 'gexp'
#... do GitHub: via devtools
#install.packages('devtools', dep=TRUE)
#library(devtools)
#install_github('ivanalaman/gexp')

#... do GitHub: via pacman
#install.packages('pacman', dep=TRUE)
#p_load_gh('ivanalaman/gexp')

#... do CRAN
# install.packages('gexp')

#.. Anexando o pacote
#... Necessita gexp >= 1.0.2
library(gexp)
packageVersion('gexp') >= '1.0.2'

#.. ___PSD/DIC___
r <- 4
set.seed(1)
psd_dic <- gexp(mu=0,
                r=r,
                fl=list(P=c('p1', 'p2'),
                        SP=c('sp1', 'sp2')),
                fe=list(f1=c(0, 1),
                        f2=c(1, 0)),
                inte=c(0,  #p1 sp1*
                       0,  #p2 sp1
                       3,  #p1 sp2
                       0), #p2 sp2
#...            Descomente abaixo p/ erro na parcela = 0
#                errp=matrix(0,
#                            nrow=2*r),
#...            Descomente abaixo p/ erro experimental = 0
#                err=matrix(0,
#                           nrow=2^2*r),
                round=0,
                type='SPE',
                design='CRD'
                )

par(mfrow=c(1, 1))
plot(psd_dic,
     main='PSD sob DIC')

# Gráfico dos efeitos
par(mfrow=c(2, 1))
with(psd_dic$dfm,
     interaction.plot(P, SP, Y1,
                      xlab='P',
                      ylab='Y1',
                      col='blue3',
                      lwd=2,
                      ylim=c(min(Y1), max(Y1))))

with(psd_dic$dfm,
     interaction.plot(SP, P, Y1,
                      xlab='SP',
                      ylab='Y1',
                      col='blue3',
                      lwd=2,
                      ylim=c(min(Y1), max(Y1))))

# *: Para entender o argumento 'inte' (efeito da interação)
(expand.grid(P=c('p1', 'p2'),
             SP=c('sp1', 'sp2')))

# ANOVA
# O número de repetições é insuficiente para inferências consistentes.
# As análises tem finalidades didáticas apenas!
av_psd_dic <- aov(Y1 ~ P + SP + P*SP + Error(P:r),
                  psd_dic$dfm)  # Desconsidere o warning 'modelo Error() é singular'
summary(av_psd_dic)

# TCMM
library(TukeyC)
cv(av_psd_dic)

# p1:SP
summary(TukeyC(av_psd_dic,
               which='P:SP',
               fl1=1))

# p2:SP
summary(TukeyC(av_psd_dic,
               which='P:SP',
               fl1=2))

# sp1:P
summary(TukeyC(av_psd_dic,
            error='Within/P:r',  # Média ponderada dos erros
               which='SP:P',
               fl1=1))

# sp2:P
summary(TukeyC(av_psd_dic,
            error='Within/P:r',  # Média ponderada dos erros
               which='SP:P',
               fl1=2))

# Plot
par(mfrow=c(2, 2))
# p1:SP
plot(TukeyC(av_psd_dic,
            which='P:SP',
            fl1=1),
     main='p1:SP',
     disp='mm',
     xlab='',
     ylab='')

# p2:SP
plot(TukeyC(av_psd_dic,
            disp='mm',
            which='P:SP',
            fl1=2),
     main='p2:SP',
     disp='mm',
     xlab='',
     ylab='')

# sp1:P
plot(TukeyC(av_psd_dic,
            error='Within/P:r',  # Média ponderada dos erros
            disp='mm',
            which='SP:P',
            fl1=1),
     main='sp1:P',
     disp='mm',
     xlab='tra',
     ylab='Y1')

# sp2:P
plot(TukeyC(av_psd_dic,
            error='Within/P:r',  # Média ponderada dos erros
            disp='mm',
            which='SP:P',
            fl1=2),
     main='sp2:P',
     xlab='',
     ylab='')

#.. ___PSD/DBC___
set.seed(6)
psd_dbc <- gexp(mu=10,
                r=1,
                fl=list(P=c('p1', 'p2'),
                        SP=c('sp1', 'sp2')),
                fe=list(f1=c(0, 3),
                        f2=c(3, 0)),
                inte=c(0,  #p1 sp1 *
                       0,  #p2 sp1
                       5,  #p1 sp2
                       0), #p2 sp2
                blke=c(1, 2, 3, 4, 5),
                blkl=list(BLO=paste('b',
                                    1:5,
                                    sep='_')),
                round=1,
                type='SPE',
                design='RCBD',
               )
str(psd_dbc)
summary(psd_dbc)

par(mfrow=c(1, 1))
plot(psd_dbc,
     main='PSD sob DBC')

# Gráfico dos efeitos
par(mfrow=c(2, 1))
with(psd_dbc$dfm,
     interaction.plot(P, SP, Y1,
                      xlab='P',
                      ylab='Y1',
                      col='blue3',
                      lwd=2,
                      ylim=c(min(Y1), max(Y1))))

with(psd_dbc$dfm,
     interaction.plot(SP, P, Y1,
                      xlab='SP',
                      ylab='Y1',
                      col='blue3',
                      lwd=2,
                      ylim=c(min(Y1), max(Y1))))

# *: Para entender o argumento 'inte' (efeito da interação)
(expand.grid(P=c('p1', 'p2'),
             SP=c('sp1', 'sp2')))

# ANOVA
# O número de repetições é insuficiente para inferências consistentes.
# As análises tem finalidades didáticas apenas!
av_psd_dbc <- aov(Y1 ~ BLO + P + SP + P*SP + Error(P:BLO),
                  psd_dbc$dfm)
summary(av_psd_dbc)

library(TukeyC)
cv(av_psd_dbc)

# TCMM
# p1:SP
summary(TukeyC(av_psd_dbc,
               which='P:SP',
               fl1=1))

# p2:SP
summary(TukeyC(av_psd_dbc,
               which='P:SP',
               fl1=2))

# sp1:P
summary(TukeyC(av_psd_dbc,
               error='Within/P:BLO',  # Média ponderada dos erros
               which='SP:P',
               fl1=1))

# sp2:P
summary(TukeyC(av_psd_dbc,
               error='Within/P:BLO',  # Média ponderada dos erros
               which='SP:P',
               fl1=2))

# Plot
par(mfrow=c(2, 2))
#p1:SP
plot(TukeyC(av_psd_dbc,
            which='P:SP',
            fl1=1),
     main='p1:SP',
     disp='mm',
     xlab='',
     ylab='')

#p2:SP
plot(TukeyC(av_psd_dbc,
            disp='mm',
            which='P:SP',
            fl1=2),
     main='p2:SP',
     disp='mm',
     xlab='',
     ylab='')

#sp1:P
plot(TukeyC(av_psd_dbc,
            error='Within/P:BLO',  # Média ponderada dos erros
            disp='mm',
            which='SP:P',
            fl1=1),
     main='sp1:P',
     disp='mm',
     xlab='tra',
     ylab='Y1')

#sp2:P
plot(TukeyC(av_psd_dbc,
            error='Within/P:BLO',  # Média ponderada dos erros
            disp='mm',
            which='SP:P',
            fl1=2),
     main='sp2:P',
     xlab='',
     ylab='')

#.. ___PSD/DQL___
# O número de tratamentos, linhas e colunas desse experimento deve ser maior!
# Foi mantido pequeno apenas para finalidades didáticas.
set.seed(3)
psd_dql <- gexp(mu=10,
                r=1,
                fl=list(P=c('p1', 'p2', 'p3'),
                        SP=c('sp1', 'sp2')),
                fe=list(f1=c(0, 3, 1),
                        f2=c(3, 0)),
                inte=c(0,  #p1 sp1
                       0,  #p2 sp1
                       0,  #p3 sp1
                       5,  #p1 sp2
                       0,  #p2 sp2
                       0), #p3 sp2
                rowe=c(1, 2, 3),
                rowl=list(LIN=paste0('l',
                                      1:3)),
                cole=c(3, 2, 1),
                coll=list(COL=paste0('c',
                                     1:3)),
                round=1,
                type='SPE',
                design='LSD',
                )
str(psd_dql)
summary(psd_dql)

par(mfrow=c(1, 1))
plot(psd_dql,
     main='PSD sob DQL')

# Gráfico dos efeitos
par(mfrow=c(2, 1))
with(psd_dql$dfm,
     interaction.plot(P, SP, Y1,
                      xlab='P',
                      ylab='Y1',
                      col='blue3',
                      lwd=2,
                      ylim=c(min(Y1), max(Y1))))

with(psd_dql$dfm,
     interaction.plot(SP, P, Y1,
                      xlab='SP',
                      ylab='Y1',
                      col='blue3',
                      lwd=2,
                      ylim=c(min(Y1), max(Y1))))

# *: Para entender o argumento 'inte' (efeito da interação)
(expand.grid(P=c('p1', 'p2', 'p3'),
             SP=c('sp1', 'sp2')))

# ANOVA
# O número de repetições é insuficiente para inferências consistentes.
# As análises tem finalidades didáticas apenas!
av_psd_dql <- aov(Y1 ~ LIN + COL + P + SP + P*SP + Error(P:LIN:COL),
                  psd_dql$dfm)
summary(av_psd_dql)

library(TukeyC)
cv(av_psd_dql)

# TCMM
# p1:SP
summary(TukeyC(av_psd_dql,
               which='P:SP',
               fl1=1))

# p2:SP
summary(TukeyC(av_psd_dql,
               which='P:SP',
               fl1=2))

# p3:SP
summary(TukeyC(av_psd_dql,
               which='P:SP',
               fl1=3))

# sp1:P
summary(TukeyC(av_psd_dql,
               error='Within/P:LIN:COL',  # Média ponderada dos erros
               which='SP:P',
               fl1=1))

# sp2:P
summary(TukeyC(av_psd_dql,
               error='Within/P:LIN:COL',  # Média ponderada dos erros
               which='SP:P',
               fl1=2))

# Plot
par(mfrow=c(2, 3))
# p1:SP
plot(TukeyC(av_psd_dql,
            which='P:SP',
            fl1=1),
     main='p1:SP',
     disp='mm',
     xlab='',
     ylab='')

# p2:SP
plot(TukeyC(av_psd_dql,
            disp='mm',
            which='P:SP',
            fl1=2),
     main='p2:SP',
     disp='mm',
     xlab='',
     ylab='')

#p3:SP
plot(TukeyC(av_psd_dql,
            disp='mm',
            which='P:SP',
            fl1=3),
     main='p3:SP',
     disp='mm',
     xlab='',
     ylab='')

# sp1:P
plot(TukeyC(av_psd_dql,
            error='Within/P:LIN:COL',  # Média ponderada dos erros
            disp='mm',
            which='SP:P',
            fl1=1),
     main='sp1:P',
     disp='mm',
     xlab='tra',
     ylab='Y1')

# sp2:P
plot(TukeyC(av_psd_dql,
            error='Within/P:LIN:COL',  # Média ponderada dos erros
            disp='mm',
            which='SP:P',
            fl1=2),
     main='sp2:P',
     xlab='',
     ylab='')