##-----------------------------------------------------------------------------
## Knitr settings. Do not run.

opts_chunk$set(
    tidy=FALSE,
    cache=FALSE,
    ## fig.path="figures/map",
    fig.width=7,
    fig.height=6,
    dpi=100,
    fig.align="center",
    dev.args=list(family="Palatino")
    )
options(width=70)

##-----------------------------------------------------------------------------
## Pacotes necessários.

pkg <- c("latticeExtra", "plyr", "doBy", "multcomp")
sapply(pkg, require, character.only=TRUE)

##-----------------------------------------------------------------------------
## Definições gráficas.

## Definições de preenchimentos, linhas, etc.
mycol <- brewer.pal(6, "Set1")
ps <- list(
    box.rectangle=list(col=1, fill=c("gray70")),
    box.umbrella=list(col=1, lty=1),
    dot.symbol=list(col=1, pch=19),
    dot.line=list(col="gray50", lty=3),
    plot.symbol=list(col=1, cex=1.2),
    plot.line=list(col=1),
    plot.polygon=list(col="gray95"),
    superpose.line=list(col=mycol),
    superpose.symbol=list(col=mycol, pch=c(1:5)),
    superpose.polygon=list(col=mycol),
    strip.background=list(col=c("gray80","gray50")))
trellis.par.set(ps)
## show.settings()

## Remove objetos da mémoria para começar vazia.
rm(list=ls())

##-----------------------------------------------------------------------------

## Endereço.
url <- "http://www.leg.ufpr.br/~walmes/data/zimmermann_fusarium.txt"
if(Sys.info()["user"]=="walmes"){
    url <- "~/Dropbox/CursoR/DADOS/zimmermann_fusarium.txt"
}

## Leitura e inspeção.
da <- read.table(file=url, header=TRUE, sep="\t")
names(da) <- substr(names(da), 1, 4)
str(da)

## Primeiras linhas.
head(da)

##-----------------------------------------------------------------------------
## Frequência.

cbind(xtabs(~trat, data=da))

cbind(xtabs(~fung+apli, data=da))

##----------------------------------------------------------------------
## IMPORTANTE: não é obrigatório mas conveniente que as testemunhas
## sejam os últimos níveis.

da$apli <- factor(da$apli,
                  levels=levels(da$apli)[c(1,4,2,3,5)])

levels(da$apli)
levels(da$fung) ## Em ordem (por sorte).
levels(da$trat) ## Em ordem (por sorte).

##----------------------------------------------------------------------
## Gráficos.

xyplot(infe/40~fung, groups=apli, data=da, jitter.x=TRUE,
       ylab="Proporção amostral de sementes infectadas (total de 40)",
       xlab="Fungicida usado no controle", auto.key=list(columns=2))

##-----------------------------------------------------------------------------
## Especificação do modelo.

## Será análisado como modelo gaussiano para comparações com o GLM
## binomial.

## O experimento foi instalado em DIC, assim, motivado pelo desenho e
## termos em estudo, o modelo tem os seguintes efeitos.
## 
##   \eta_{ij} = \mu+\alpha_{i}+\tau_{j}+\gamma_{ij}, se I
##   \eta_{k} = \mu+\delta_{k}, se II
##
## Em que I representa a porção fatorial (com dois efeitos principais e
## interação) e II a porção adicional (as testemunhas).

##----------------------------------------------------------------------
## Contrastes.

K <- C(da$fung, contr=contr.helmert)
K <- C(da$fung, contr=contr.treatment)
K <- attr(K, "contrasts")

## O que o constraste representa?
MASS::fractions(solve(cbind(1, K)))

## Os coeficientes podem ser alterados, caso as hipóteses sejam outras.
## Por exemplo, pode se usar os seguintes contrastes.

U <- rbind(c(1,0,0,0,0),  ## Média no primeiro nível.
           c(-1,1,0,0,0), ## Contraste entre 1 e 2.
           c(-1,0,1,0,0), ## Contraste entre 2 e 3.
           c(1,1,1,-3,0), ## {1,2,3} vs 1ra testemunha.
           c(1,1,1,0,-3)) ## {1,2,3} vs 2da testemunha
U

K <- solve(U)[,-1]
MASS::fractions(K)

## Atribuindo esses contrastes.
contrasts(da$fung) <- K
contrasts(da$apli) <- K

##----------------------------------------------------------------------
## Ajuste do modelo.

## ctr <- list(fung=contr.helmert, apli=contr.helmert)
## ctr <- list(fung=contr.treatment, apli=contr.treatment)

m0 <- lm(infe/40~fung*apli,
         ## contrasts=ctr,
         data=da)

## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

## Com transformação justificada para estabilização da variância para
## variáveis de proporção em que se assume binomial.
m0 <- update(m0, asin(sqrt(.))~.)

## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

## Quadro de análise de deviance.
anova(m0, test="F")

## NOTE: os graus de liberdade não são os corretos. Os tratamentros
## adicionais são parâmetros dentro do vetor de parâmetros do primeiro
## fator declarado. No segundo, por confundimento completo, se tornam
## não estimáveis e por isso não alteram os graus de liberdade.

##----------------------------------------------------------------------
## Efeitos não estimáveis.

coef(m0)

## Dimensões diferentes no vetor de estimativas e matrix de
## covariância das estimativas.
c(length(coef(m0)), nrow(vcov(m0)))

##----------------------------------------------------------------------
## TRUQUE: Obter a matriz do modelo, remover as colunas que correspondem
## aos efeitos não estimáveis e ajudar o modelo passando uma matriz.

X <- model.matrix(m0)
a <- attr(X, "assign") ## Será usado para decompor a deviance.
i <- is.na(coef(m0))

levelplot(t(X[nrow(X):1,]), scales=list(x=list(rot=90)))+
    layer(panel.abline(v=which(i), lty=2))

X <- X[,!i]
a <- a[!i]

levelplot(t(X[nrow(X):1,]), scales=list(x=list(rot=90)))

##----------------------------------------------------------------------

## Modelo declarado por matriz, apenas colunas de efeitos estimáveis.
m1 <- update(m0, .~0+X)

## Iguais deviances pois são o mesmo modelo.
c(deviance(m0), deviance(m1))

c(length(coef(m1)), nrow(vcov(m1)))

a
a[4:5] <- max(a)+1:2

L <- diag(length(a))
Ls <- by(data=L, INDICES=a, FUN=as.matrix)

Ls <- Ls[-1] ## Remove intercepto.
names(Ls) <- c("fung", "apli", "fung:apli",
               "fung vs test1", "fung vs test2")
Ls <- Ls[c(1,4,5,2,3)]
names(Ls)

lapply(Ls,
       function(L){
           summary(glht(m1, linfct=L), test=Ftest())
       })

## anova(m0, test="F") ## Sequencial.
drop1(m0, test="F", scope=.~.) ## Marginal.

## O p-valor de F do Wald é igual ao p- valor de F da drop1(). São
## iguais porque o modelo é de classe "lm" e nessa classe a LL é uma
## função quadrática, então a estatística F de LL e F de Wald são
## identicas.

## MAS CUIDADO! Essas hipóteses são marginais e algumas delas só fazem
## sentido em casos muito específicos. Por exemplo, não existe sentido
## em testar o efeito de um efeito principal ajustado para interações
## que o contenham. Também não faz sentido testar a média do fatorial
## contra uma testemunha se a interação for signifcativa.

##----------------------------------------------------------------------
## Criar um modelo de mentira que seja um fatorial completo.

## Atribui indices incrementais as linhas.
rownames(da) <- NULL
dau <- unique(subset(da, select=c(fung, apli)))
Xu <- X[rownames(dau),]

Xu

## 3x3+2 = 11 celas experimentais. Se fosse um experimento em blocos,
## não haveria linhas únicas ao considerar o termo bloco, então uma rota
## diferente seria necessária para aplicar a média no efeito dos blocos.

ic <- confint(glht(m1, linfct=Xu), calpha=univariate_calpha())$confint
ic

pred <- cbind(dau, ic)
str(pred)

##----------------------------------------------------------------------
## Gráfico.
 
segplot(fung:apli~lwr+upr, data=pred, horizontal = FALSE,
        ylab="asin(sqrt(Proprorção de infecção))",
        xlab="Tratamento", scales=list(x=list(rot=90)),
        centers=Estimate, draw.bands=FALSE)

##----------------------------------------------------------------------
## Gráfico aprimorado.

require(wzRfun)

pred$desloc <- 0.05*c(0,-1,1,0,0)[match(pred$apli,
                                        levels(da$apli))]
 
xyplot(asin(sqrt(infe/40))~fung, groups=apli, data=da,
       jitter.x=TRUE, amount=0.05,
       xlab="Fungicida",
       ylab=expression(asin(sqrt("Proporção de infecção"))),
       prepanel=prepanel.cbH,
       ly=pred$lwr, uy=pred$upr,
       auto.key=list(
           columns=2, type="o", divide=1,
           title="Aplicação", cex.title=1.2,
           lines=TRUE, points=FALSE))+
  as.layer(xyplot(Estimate~fung, groups=apli,
                  data=pred, type=c("p","l"),
                  pch=15,
                  ly=pred$lwr, uy=pred$upr,
                  desloc=pred$desloc,
                  cty="bars",
                  panel.groups=panel.cbH,
                  panel=panel.superpose))

##----------------------------------------------------------------------
## Testar os tratamentos entre si.

Xu
dau

Xm <- Xu
rownames(Xm) <- with(dau, paste(fung, apli, sep=":"))

d <- combn(nrow(Xm), 2)
ncol(d)

## Matriz com as diferenças entre linhas da anterior, ou seja, os
## contrastes.

Xc <- sapply(as.data.frame(d),
              function(l){
                  K <- list(Xm[l[1],]-Xm[l[2],])
                  names(K) <- paste(rownames(Xm)[l], collapse="-")
                  return(K)
              })
Xc <- do.call(rbind, Xc)
## Xc

## summary(glht(m1, linfct=Xc), test=adjusted(type="fdr"))

##----------------------------------------------------------------------
## Para facilitar, dado que existe interação, declarar o modelo de um
## fator é equivalente.

m2 <- update(m0, .~trat)
anova(m2)

c(deviance(m2), deviance(m1))

g0 <- summary(glht(m2, linfct=mcp(trat="Tukey")),
              test=adjusted(type="fdr"))

## cld(g0)
## str(cld(g0))
cbind(cld(g0, decreasing=TRUE)$mcletters$Letters)

##-----------------------------------------------------------------------------
## Versões dos pacotes e data do documento.

sessionInfo()
Sys.time()