##=============================================================================
## AF-722: Modelagem e análise de dados experimentais com o programa
##         computacional R
##                                                       www.leg.ufpr.br/af722
##
## Programa de Pós Graduação em Agronomia - 2014/1
##                                                      Prof Dr. Walmes Zeviani
##                                                      LEG - DEST - UFPR
##
## Exemplo de análise de experimento com mais de um fator. No caso é um
## experimento com 2 fatores, níveis completamente cruzados, instalado
## em DBC.
##
##=============================================================================

##-----------------------------------------------------------------------------
## Definições da sessão.

## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape")
sapply(pkg, require, character.only=TRUE)

source("http://dl.dropboxusercontent.com/u/48140237/apc.R")
source("http://dl.dropboxusercontent.com/u/48140237/equallevels.R")

##-----------------------------------------------------------------------------
## Dados de experimento com soja sob 3 níveis de umidade do solo (massa
## de água/massa de solo) e níveis de potássio fornecidos na adubação. A
## hipótese sob estudo é que presença de potássio dá suporte para a
## cultura superar a ocorrência de períodos sem chuva. Experimento
## conduzido em casa de vegetação com 5 blocos, 2 plantas por
## u.e. (vaso).

## Para acessar o artigo (pdf online).
browseURL("http://www.scielo.br/pdf/rca/v43n2/a03v43n2.pdf")

da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
                 header=TRUE, sep="\t", dec=",")
str(da)

xyplot(rengrao~potassio, groups=agua, data=da,
       type=c("p","a"))

xyplot(rengrao~potassio|agua, data=da,
       type=c("p","a"))

##-----------------------------------------------------------------------------
## Por uma questão didática, os fatores que são métricos (quatitativos) serão
## considerados como categóricos (qualitativos) para demonstrar como
## proceder a análise. Em outro script eles serão analisados como
## métricos que é mais recomendado.

## Fatores.
da$A <- factor(da$agua)
levels(da$A)
da$K <- factor(da$potassio)
levels(da$K)

## Espeficicação do modelo.
m0 <- lm(rengrao~bloco+A*K, data=da)

## Resíduos.
par(mfrow=c(2,2)); plot(m0); layout(1)

## Remove obs 74 (uma planta que definou).
m0 <- lm(rengrao~bloco+A*K, data=da[-74,])
par(mfrow=c(2,2)); plot(m0); layout(1)

## Quadro de anova.
anova(m0)

## Não inferir sobre os efeitos de menor ordem se os termos estiverem
## envolvidos em algura interação (termo de ordem maior). A interação
## mascara os efeitos principais. Não havendo interação pode-se inferir
## sobre os efeitos principais.

##-----------------------------------------------------------------------------
## Estudar o efeito dos níveis de potássio fixando o nível de água.

## Médias marginais para o efeito dos níveis de água. Marginal é porque
## cada nível de potássio contribui com 1/5.
ls <- LSmeans(m0, effect="A"); ls

## Média sob o efeito de água fixando potássio em 0.
ls <- LSmeans(m0, effect="A", at=list(K="0")); ls

## Média sob o efeito de potássio fixando água em 37.5.
ls <- LSmeans(m0, effect="K", at=list(A="37.5")); ls

## As médias para todas as combinações de efeito entre água e potássio.
ls <- LSmeans(m0, effect=c("K","A")); ls

## Informação do objeto pela sua estrutura.
typeof(ls)
names(ls)
str(ls)
ls$coef
ls$grid
ls$K

tb <- cbind(ls$grid, ls$coef); tb
tb$K <- as.integer(tb$K)

## Para dispor em um formato retangular.
cast(tb, K~A, value="estimate")

##-----------------------------------------------------------------------------
## Teste sobre as médias (contrastes par a par).

ls <- LSmeans(m0, effect=c("K"), at=list(A="37.5")); ls
## ls <- LSmeans(m0, effect=c("K"), at=list(A="50")); ls
## ls <- LSmeans(m0, effect=c("K"), at=list(A="62.5")); ls
str(ls)
X <- ls$K

## Estimativa das médias matricialmente.
X%*%coef(m0)

Xc <- apc(X)
choose(5, 2)
dim(Xc)

## Estimativa dos contrastes matricialmente.
Xc%*%coef(m0)

## Teste para hipótese nula de não haver diferenças.
g0 <- summary(glht(m0, linfct=Xc))
g0

## Para ter o resumo compacto.
g0$focus <- "K"; g0$type <- "Tukey"
cld(g0)

## Para ter as médias com IC (a LSmeans também serve).
gm <- glht(m0, linfct=X)
ci <- confint(gm, calpha=univariate_calpha()); ci

## Montar tabela com médias, ic, letras e colocar em gráfico.
grid <- ls$grid
grid <- equallevels(grid, da)
grid <- cbind(grid, ci$confint, let=cld(g0)$mcletters$Letters)
grid

## Edita a média com letra ao lado.
grid$m <- with(grid, sprintf("%0.2f %s", Estimate, let))
grid

## Gráfico de segmentos.
segplot(K~lwr+upr, centers=Estimate, data=grid,
        horizontal=FALSE, draw=FALSE)+
    layer(panel.text(x=as.integer(grid$K), y=grid$Estimate,
                     labels=grid$m, pos=2))

##-----------------------------------------------------------------------------
## Fazendo os 3 desdobramentos de uma forma automática.

a <- levels(da$A)
g <- lapply(a,
            function(a){
                ls <- LSmeans(m0, effect=c("K"), at=list(A=a))
                X <- ls$K
                Xc <- apc(X)
                g0 <- summary(glht(m0, linfct=Xc))
                g0$focus <- "K"; g0$type <- "Tukey"
                gm <- glht(m0, linfct=X)
                ci <- confint(gm, calpha=univariate_calpha())
                grid <- ls$grid
                grid <- equallevels(grid, da)
                grid <- cbind(grid, ci$confint, let=cld(g0)$mcletters$Letters)
                grid$m <- with(grid, sprintf("%0.2f %s", Estimate, let))
                return(grid)
            })
str(g)

g <- do.call(rbind, g)

panel.means <- function(x, y, z, centers, subscripts, labels, ...){
    panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...)
    panel.text(x=as.integer(z[subscripts]), y=centers[subscripts],
               labels=labels[subscripts], srt=90, adj=c(0.5,-0.4), cex=0.8)
}

sub <-
    "Médias com mesma letra não diferem entre si em cada nível de água à 5%."
fl <- sprintf("Umidade de %s", paste0(levels(da$A), "%"))

## Para fazer o gráfico em escala cinza.
dev.off()
lattice.options(
    default.theme=modifyList(
        standard.theme(color=FALSE),
        list(strip.background=list(col="gray80"))))

## Gráfico de segmentos.
segplot(K~lwr+upr|A, centers=Estimate, data=g,
        horizontal=FALSE, draw=FALSE, layout=c(NA,1),
        strip=strip.custom(factor.levels=fl),
        xlab=expression("Potássio no solo"~(mg~dm^{-3})),
        ylab=expression("Rendimento de grãos"~(g~vaso^{-1})),
        sub=list(sub, font=3, cex=0.8), labels=g$m,
        panel=panel.means)

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