#===============================================================================
#. Name   : tcmm_sk_psd
#. Author : J.C.Faria
# Date   : 05/08/2021 19:09:39
# Aim    : Simular experimentos qualitativos em parcelas subdivididas nos
#          delineamentos experimentais básicos (DIC, DBC e DQL) e analisar
#          pelo algoritmo de ScottKnott.
# Mail   : <<< joseclaudio.faria@gmail.com >>>
#===============================================================================

#. Carregando pacotes necessários
# gexp
#install.packages('gexp')
library(gexp)

# ScottKnott
#install.packages('ScottKnott')
library(ScottKnott)

#.. ___PSD/DIC___
r <- 4
set.seed(2019)
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
#p1:SP
summary(SK(av_psd_dic,
           error='Within',
           which='P:SP',
           fl1=1))

#p2:SP
summary(SK(av_psd_dic,
           error='Within',
           which='P:SP',
           fl1=2))

#sp1:P
summary(SK(av_psd_dic,
           error='Within', #... É preciso combinar os erros!
           which='SP:P',
           fl1=1))

#sp2:P
summary(SK(av_psd_dic,
           error='Within', #... É preciso combinar os erros!
           which='SP:P',
           fl1=2))

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

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

#sp1:P
plot(SK(av_psd_dic,
        error='Within',
        which='SP:P', #... É preciso combinar os erros!
        fl1=1),
     disp='mm',
     main='sp1:P',
     xlab='tra',
     ylab='Y1')

#sp2:P
plot(SK(av_psd_dic,
        error='Within', #... É preciso combinar os erros!
        which='SP:P',
        fl1=2),
     disp='mm',
     main='sp2:P',
     xlab='',
     ylab='')

#.. ___PSD/DBC___
set.seed(2019)
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)

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

#p2:SP
summary(SK(av_psd_dbc,
           error='Within',
           which='P:SP',
           fl1=2))

#sp1:P
summary(SK(av_psd_dbc,
           error='Within',
           which='SP:P', #... É preciso combinar os erros!
           fl1=1))

#sp2:P
summary(SK(av_psd_dbc,
           error='Within',
           which='SP:P', #... É preciso combinar os erros!
           fl1=2))

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

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

#sp1:P
plot(SK(av_psd_dbc,
        error='Within',
        which='SP:P', #... É preciso combinar os erros!
        fl1=1),
     disp='mm',
     main='sp1:P',
     xlab='tra',
     ylab='Y1')

#sp2:P
plot(SK(av_psd_dbc,
        error='Within',
        which='SP:P', #... É preciso combinar os erros!
        fl1=2),
     disp='mm',
     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(2019)
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=paste('l',
                                     1:3,
                                     sep='')),
                cole=c(3, 2, 1),
                coll=list(COL=paste('c',
                                     1:3,
                                     sep='')),
                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)

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

#p2:SP
summary(SK(av_psd_dql,
           error='Within',
           which='P:SP',
           fl1=2))

#p3:SP
summary(SK(av_psd_dql,
           error='Within',
           which='P:SP',
           fl1=3))

#sp1:P
summary(SK(av_psd_dql,
           error='Within',
           which='SP:P', #... É preciso combinar os erros!
           fl1=1))

#sp2:P
summary(SK(av_psd_dql,
           error='Within',
           which='SP:P', #... É preciso combinar os erros!
           fl1=2))

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

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

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

#sp1:P
plot(SK(av_psd_dql,
        error='Within',
        which='SP:P', #... É preciso combinar os erros!
        fl1=1),
     disp='mm',
     main='sp1:P',
     xlab='tra',
     ylab='Y1')

#sp2:P
plot(SK(av_psd_dql,
        error='Within',
        which='SP:P', #... É preciso combinar os erros!
        fl1=2),
     disp='mm',
     main='sp2:P',
     xlab='',
     ylab='')