---
title: "Análise de Experimentos em Delineamento Inteiramente Casualizado"
author: "PET Estatística UFPR"
vignette: >
  %\VignetteIndexEntry{Experimentos em Delineamento Inteiramente Casualizado}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
source("_setup.R")
```

## Análise exploratória

Nessa *vignette* vamos considerar `BanzattoQd3.2.1`, já disponível no
pacote [**labestData**].

```{r}
library(labestData)
ls("package:labestData")
```
```{r, eval=FALSE}
help(BanzattoQd3.2.1, help_type = "html")
```

Este conjunto de dados é de um experimento instalado em delineamento
inteiramente casualizado (DIC) para estudar diferentes produtos para 
controle de pulgão na cultura do pepino. Foram avaliados 4 produtos 
(`trat`), e mais uma testemunha, para controle do pulgão 
(*Aphis gossypii* Glover). A variável resposta desse experimento foi o 
número de pulgões encontrados, uma variável resposta do tipo contagem.

Em caso de não ter o pacote **labestData**, você pode ler a tabela de
dados diretamente do repositório. No entanto, recomendamos que instale o
pacote para ter acesso a todos os datasets e às documentações

Antes de iniciar a análise dos dados do experimento, é muito
recomendável fazer uma análise exploratória que serve basicamente para

  * Verificar se os dados não possuem erros de digitação;
  * Antecipar a presença de observações discrepantes (outliers);
  * Verificar as escalas das variáveis;
  * Antecipar a existência de efeito dos fatores estudados;
  * Compreender os dados.

```{r}
#-----------------------------------------------------------------------
# Ler a partir do repositório do labestData.

# url <- paste0("https://gitlab.c3sl.ufpr.br/pet-estatistica",
#               "/labestData/raw/devel/data-raw/BanzattoQd3.2.1.txt")
#
# BanzattoQd3.2.1 <-
#     read.table(file = url, sep = "\t", header = TRUE,
#                colClasses = c("factor", "integer", "integer"))

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

# Estrutura do objeto.
str(BanzattoQd3.2.1)

# Tabela de frequência para os tratamentos.
xtabs(~trat, data = BanzattoQd3.2.1)

# Dados desemplilhados.
t(unstack(x = BanzattoQd3.2.1, form = pulgoes ~ trat))

# Média e desvio-padrão das observações em cada nível.
aggregate(pulgoes ~ trat,  data = BanzattoQd3.2.1,
          FUN = function(y) {
              c(mean = mean(y), dp = sd(y))
          })

library(lattice)
packageVersion("lattice")

# Diagrama de dispersão.
xyplot(pulgoes ~ reorder(trat, pulgoes), data = BanzattoQd3.2.1,
       xlab = "Produtos para controle de pulgão",
       ylab = "Número de pulgões 36hs após pulverização",
       scales = list(x = list(rot = 90)))

```

Pela análise exploratória temos que esse experimento é balanceado pois
tem número igual de repetições de cada produto. O diagrama de dispersão
mostra os valores observados do número de pulgões em função dos
produtos. Os níveis foram ordenados de acordo com os números médios de 
pulgões. É imediato perceber a partir do gráfico que a variabilidade do 
número de pulgões aumenta com a média. Isso é importante porque adianta um 
possível problema com o pressuposto de homogeneidade de variâncias do modelo.

## Especificação e ajuste do modelo

O modelo estatístico correspondente ao delineamento é

$$
\begin{aligned}
  Y|\text{produtos} &\,\sim \text{Normal}(\mu_i, \sigma^2) \newline
  \mu_i &= \mu+\tau_i
\end{aligned}
$$

em que \(\tau_i\) é o efeito do produto \(i\) sobre a variável resposta
\(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos produtos,
\(\mu_i\) é a média para cada nível de produto e \(\sigma^2\) é a
variância das observações ao redor desse valor médio (ou variância
residual).

Para fazer a estimação desses parâmetros é necessário que seja
considerada alguma restrição paramétrica. Para isso, no R por default
considera-se o efeito do primeiro nível como zero ($\tau_1 = 0$). É
possível mudar isso nas opções de contrastes do R. Outra parametrização,
bem conhecida e usada, é a retrição do tipo soma zero dos efeitos ($\sum
\tau_i = 0$) e, para essa situação, $\mu$ representa a média geral do
experimento.

Existem duas funções no R para ajustar esse modelo: `lm()` e
`aov()`. Embora muitos acreditem que `aov()` deva ser preferida, por
causa do nome AOV (*Analysis Of Variance*), a função `lm()` é tão
apropriada quanto. Inclusive, objetos vindos da `lm()` dispõem de mais
métodos que a da `aov()` em pacotes para fazer comparações múltiplas, por
exemplo. Sendo assim, vamos usar a `lm()`.

```{r}
#-----------------------------------------------------------------------
# Ajuste do modelo.

m0 <- lm(pulgoes ~ trat, data = BanzattoQd3.2.1)

# Estimativas dos efeitos. Restrição de zerar primeiro nível.
cbind(coef(m0))

# Matriz de contrastes sob a retrição zerar primeiro nível.
K <- cbind(1, contrasts(BanzattoQd3.2.1$trat))
K

# Médias estimadas pelo modelo.
K %*% coef(m0)

# Médias amostrais.
aggregate(pulgoes ~ trat, data = BanzattoQd3.2.1, FUN = mean)
```

Conforme comentado, devido à restrição usada, os coeficientes estimados
não são as médias em cada nível. Para a restrição do primeiro nível ter
efeito nulo, o coeficiente `(Intercept)` é a média do primeiro nível, no
caso ` `r levels(BanzattoQd3.2.1$trat)[1]` `. Os demais coeficientes
são estimativas dos constrastes entre médias de cada um dos níveis
restantes com relação ao primeiro nível.

Como nesse exemplo existe uma testemunha, seria conveniente que a
testemunha fosse o primeiro nível para então termos os contrastes como
resultado do ajuste do modelo. No entanto, como a testemunha é o
último nível, pode-se mudar o tipo de contraste para `contr.SAS`, cuja
restrição é zerar o efeito do último nível. Na sequência temos os 
resultados das duas opções.

```{r}
m1 <- update(m0, contrasts = list(trat = "contr.SAS"))
coef(m1)

BanzattoQd3.2.1 <- transform(BanzattoQd3.2.1,
                             trat = relevel(trat, ref = "Testemunha"))

levels(BanzattoQd3.2.1$trat)

m0 <- update(m0,  data = BanzattoQd3.2.1)
coef(m0)
```

## Avaliação dos pressupostos

Para fazer inferência a partir do modelo ajustado aos dados, é
necessário verificar se os pressupostos do modelo foram atendidos. Caso 
não sejam, as inferências produzidas não são fiéis, pois a distribuição
amostral das estatísticas de teste não são as mesmas sob violação dos
pressupostos.

```{r}
#-----------------------------------------------------------------------
# Exibe o quarteto fantástico da avaliação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0); layout(1)
```

O conjunto de gráficos acima exibe fuga dos pressupostos. A suposição
de homogeneidade de variâncias é claramente violada. Vemos o padrão cone no
gráfico 1-1 dessa matriz (*Residuals vs Fitted*). No gráfico 2-1
(*Scale-Location*), temos ênfase desse padrão porque a linha média mostra
uma relação positiva com a dispersão, representada pela raíz dos valores
absolutos dos resíduos padronizados, versus os valores ajustados
(média). Ou seja, a relação média-variância não é nula, pois se fosse, a
linha de tendência seria horizontal. Por fim, o gráfico 1-2 (Normal
Q-Q), exibe fuga do pressuposto de normalidade dos erros, influenciada 
pela relação média-variância.

Da forma como está, o modelo não tem os pressupostos atendidos. Logo não
é recomendado fazer inferências. Precisamos remediar a fuga dos
pressupostos antes e a opção mais comum é procurar uma tranformação que
estabilize a variância. Em casos como esse, as transformações do tipo
potência são capazes de estabilizar a variância. Vamos considerar a
transformação Box-Cox para isso.

A tranformação Box-Cox baseia-se na escolha do valor $\lambda$ de tal
maneira que se aumente a compatibilidade dos dados com a suposição de
normalidade. A transformação aplicada aos dados é

$$
\begin{aligned}
  y_{\text{transformado}} &=
    \frac{y^{\lambda}-1}{\lambda}, \quad \text{ se }
    \lambda\neq 0 \newline
  y_{\text{transformado}} &= \log(y), \quad \text{ se } \lambda=0.
\end{aligned}
$$

A transformação só pode ser aplicada para respostas positivas, \(y>0\)
. Caso tenha-se zeros ou valores negativos podemos somar uma
constante aos dados. Como o que interessa da transformação é parte não
linear (a potenciação), na prática apenas se aplica a função potência 
sem subtrair 1 e dividir por $\lambda$. No entanto, essa é a 
transformação de fato aplicada para calcular o valor da log-verossimilhança.
Se houver a curiosidade de fazer o cálculo da log-verrossimilhança passo
a passo, lembre-se de considerar o jacobiano da transformação.

```{r}
MASS::boxcox(m0)
abline(v = c(0, 1/3), lty = 2, col = 2)
```

A função perfil de log-verossilhança em $\lambda$ tem intervalo de
confiança (95%) que inclui os valores 0 e 1/3, correspondentes às
transformações log e raíz cúbica, respectivamente. Não há necessidade de
usar o valor exato correspondente ao máximo da função porque isso não
vai mudar a normalidade da variável transformada e, como o valor exato
é um valor quebrado, valores de $\lambda$ que correspondem
a transformações conhecidas são convenientes, como log, raíz quadrada e
raíz cúbica.

Encontrada a transformação a ser aplicada, devemos ajustar o modelo com
a variável transformada e avaliar novamente os pressupostos. Não existe
garantia de qua a transformação irá corrigir a fuga. O que ela faz é
remediar o problema, mas o resultado pode ainda não ser satisfatório.

```{r}
# Atualiza o modelo anterior passando o log na resposta.
m0 <- update(m0, log(.) ~ .)

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

Note agora que os gráficos não exibem padrões de fuga de
pressupostos. Pelo gráfico 2-1, a relação média variância agora está
estável (ou nula). Não há fugas de normalidade de acordo com o gráfico
1-2. Portanto, foi bem sucedida a aplicação da transformação log ao
número de pulgões para que o modelo tivesse os pressupostos
atendidos. Agora é possível proceder com as inferências.

## Inferências

O objetivo de ter feito o experimento é saber se existe efeito dos
produtos para o controle do pulgão na cultura do pepino. Ou seja, se o
número médio de pulgões muda conforme o produto aplicado.

O efeito dos produtos é representado pelos coeficientes $\tau_i$ no
modelo estatístico do experimento. Se não existir efeito dos produtos, os
valores estimados para $\tau_{i}$ serão próximos a zero. A hipóteses nula 
de não haver efeito dos produtos é representada por

$$
  \text{H}_{0}: \tau_i = 0, \text{para todo}\,i.
$$

Essa hipótese é avaliada pelo estatística F do quadro de análise de
variância.

```{r}
anova(m0)
```

Pelo quadro de anova (*analysis of variance*), temos que o valor da
estatística F produziu um $p$-valor bem pequeno, menor que qualquer nível 
de significancia adotado na prática. Dessa maneira, rejeitamos a hipótese
nula de não haver efeito dos produtos no log do número de pulgões ou,
por consequência, no número de pulgões.

Uma vez que existe efeito dos produtos, faz-se necessário
aprofundar o conhecimento sobre eles. As hipóteses a serem testadas
devem ser estabalecidas antes da coleta dos dados e ajuste do modelo.
Como foram estudados 4 níveis de
produto e uma testemunha, as hipóteses vão agora considerar diferenças
entre cada produto e a testemunha. Como a testemunha é o primeiro nível,
temos que testar as hipóteses individuais:

$$
  \text{H}_{0}: \mu_i - \mu_1 = 0, \text{para cada}\,i \neq 1.
$$

Da forma como ajustamos o modelo, temos as estimativas dos contrastes
disponíveis no vetor de coeficientes. O teste individual para cada
parâmetro é exibido com o `summary()` do modelo.

```{r}
summary(m0)
```
```{r, include=FALSE}
k <- nlevels(BanzattoQd3.2.1$trat) - 1
p <- 1 - (1 - 0.05)^k
```

Na tabela dos coeficientes tem-se o valor t correspondente a cada contraste. 
O `(Intercept)` é a média do log do número de pulgões para a
testemunha. Os demais representam contrastes dos produtos com relação à
testemunha. Baseado nos $p$-valores, todas as hipóteses de contraste nulo 
com a testemunha são rejeitadas. Ou seja, todos os produtos apresentam número
médio diferente da testemunha, no caso um menor valor para média
se comparado à testemunha.

Deve-se mencionar que se o nível de significância adotado é de 5% para
cada teste. Assim, a probabilidade de rejeitar pelo menos uma hipótese é $1 -
(1 - 0,05)^{`r k`} = `r p`$, ou seja, superior a 5%. Com o aumento
do número de hipóteses testadas, é consequência aumentar o nível de
significância global (probabilidade de rejeitar pelo menos uma
hipótese). Existem diferentes formas de corrigir a elevação do nível 
global de significancia, mas esse não é um objetivo desta *vignette*.

Uma adequada representação dos resultados é um gráfico com os valores
médios e intervalo de confiança para os níveis dos produtos.

```{r}
pred <- data.frame(trat = levels(BanzattoQd3.2.1$trat))
pred <- cbind(pred,
              as.data.frame(predict(m0,
                                    newdata = pred,
                                    interval = "confidence")))
pred$trat <- reorder(pred$trat, pred$fit)
head(pred)

library(latticeExtra)

segplot(trat ~ lwr + upr, centers = fit, data = pred,
        draw = FALSE, horizontal = FALSE,
        xlab = "Produtos para controle de pulgão",
        ylab = "Log do número de pulgões",
        scales = list(x = list(rot = 90)),
        panel = function(x, y, z, centers, ...) {
            panel.segplot(x = x, y = y, z = z, centers = centers, ...)
            panel.text(x = as.numeric(z), y = centers,
                       label = sprintf("%0.2f", centers),
                       srt = 90, cex = 0.8, adj = c(0.5, -0.5))
    })
```

## Gerando experimento em DIC

É simples produzir o plano de um delineamento inteiramente
casualizado. Vamos gerar o plano para executar um experimento
semelhante ao apresentado, mas com 8 repetições e mais um produto,
o `Novo` produto.

```{r, eval=FALSE}
trat <- rep(x = c(levels(BanzattoQd3.2.1$trat), "Novo"),
            times = 8)

# Opção 1: sorteie os níveis para as unidades experimentais ordenadas.
data.frame(trat = sample(trat), ue = 1:length(trat))

# Opção 2: sorteie as unidades experimentais ordenadas para os níveis.
data.frame(trat = trat, ue = sample(1:length(trat)))
```

## Informações da sessão

```{r}
sessionInfo()
```

<!------------------------------------------- -->

[**labestData**]: https://gitlab.c3sl.ufpr.br/pet-estatistica/labestData
