Experimento em delineamento de blocos casualizados

##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.

pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
         "plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)

source("lattice_setup.R")

##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.

sessionInfo()

## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)

Importação dos dados

##-----------------------------------------------------------------------------
## Lendo arquivos de dados.

## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/pimentel_batatinha.txt"

## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")

## Para facilitar, truncar nomes para 4 caracteres.
names(da) <- substr(names(da), start=1, stop=4)

## Mostra a estrutura do objeto.
str(da)

## Tabela de frequência.
xtabs(~vari+bloc, data=da)

##-----------------------------------------------------------------------------
## Análise exploratória.

xyplot(prod~vari, groups=bloc, data=da, type=c("p","a"))

Especificação e estimação do modelo

\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+\alpha_i+\tau_j \end{aligned} \]

em que \(\alpha_i\) é o efeito do bloco \(i\) e \(\tau_j\) o efeito da variedade \(j\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito de bloco e variedade, \(\mu_{ij}\) é a média para as observações em cada combinação de níveis de bloco e variedade e \(\sigma^2\) é a variância das observações ao redor desse valor médio.

##-----------------------------------------------------------------------------
## Estimação.

m0 <- lm(prod~bloc+vari, data=da)

## Estimativas e erros padrões.
summary(m0)

## Restrição paramétrica de zerar o primeiro nível.
contrasts(da$bloc)
contrasts(da$vari)


##-----------------------------------------------------------------------------
## Existe efeito das fontes de variação? Ou seja, todos os \alpha_i
## podem ser considerados iguais à zero? Todos os \tau_j podem ser
## considerados à zero.

## Quadro de análise de variância.
anova(m0)

Diagnóstico do modelo

##-----------------------------------------------------------------------------
## Análise de diagnóstico dos pressupostos via resíduos.

## As inferências acima só passam ser valiadas após confirmação de não
## haver afastamento dos pressupostos.

par(mfrow=c(2,2))
plot(m0); layout(1)

##-----------------------------------------------------------------------------
## Transformar os dados pode diminuir a fuga dos pressupostos?

MASS::boxcox(m0)

Comparações múltiplas

##-----------------------------------------------------------------------------
## Comparações múltiplas.

## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="vari", level=0.95)

## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0

## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- summary(glht(m0, linfct=mcp(vari="Tukey")),
              test=adjusted(type="fdr"))
g0

## Resumo compacto com letras.
l0 <- cld(g0, decreasing=TRUE)
l0

## Formatando média seguida de letras.
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
p0$let

##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.

xlab <- "Produção"
ylab <- "Variedades de batata"
sub <- list(
    "Médias seguidas de mesma letra indicam diferença nula à 5%.",
    font=1, cex=0.8)

p0 <- transform(p0, vari=reorder(vari, estimate))
levels(p0$vari)
p0 <- p0[order(p0$vari),]

segplot(vari~lwr+upr, centers=estimate,
        data=p0, draw=FALSE,
        xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
        panel=function(x, y, z, centers, letras, ...){
            panel.segplot(x=x, y=y, z=z, centers=centers, ...)
            panel.text(y=as.numeric(z), x=centers,
                       label=letras, pos=3)
        })