From 9e77c0a05ea9b6f9b1ffed66e6b73ecbcfcd8177 Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Fri, 17 Jan 2020 19:55:40 -0300
Subject: [PATCH] Improves first section of the chapter.

---
 04-dql.Rmd | 607 ++++++++++++++++++++++++++++++++++-------------------
 1 file changed, 394 insertions(+), 213 deletions(-)

diff --git a/04-dql.Rmd b/04-dql.Rmd
index d62001e..afc6f64 100644
--- a/04-dql.Rmd
+++ b/04-dql.Rmd
@@ -24,14 +24,17 @@ vôo (colunas, e.g. vento, visibilidade).
 knitr::include_graphics("./img/quadrado-latino-drone.png")
 ```
 
+Para demonstração da análise desse delineamento serão utilizados os
+pacotes ... TODO FIXME.
+
 ```{r, message = FALSE}
 # ATTENTION: Limpa espaço de trabalho.
 rm(list = objects())
 
 # Carrega pacotes.
-library(agricolae) # HSD.test().
 library(car)       # linearHypothesis(), Anova().
 library(multcomp)  # glht().
+library(agricolae) # HSD.test().
 library(emmeans)   # emmeans().
 library(tidyverse) # Manipulação e visualização de dados.
 
@@ -40,49 +43,74 @@ source("mpaer_functions.R")
 ls()
 ```
 
-## Variável resposta contínua
+## Resposta contínua e fator qualitativo
 
 Nesta seção será feita a análise de um conjunto de dados disponível em
 @pimentel sobre um ensaio realizado no delineamento de quadrado latino
-para comparar a produtividade de 5 variedades de cana-de-açucar. Antes
-de partir para análise é fortemente recomendável fazer a análise
+para comparar a produtividade de 5 variedades de cana-de-açúcar (fator é
+qualitativo).
+
+Antes de partir para análise é fortemente recomendável fazer a análise
 exploratória dos dados que serve como fase de inspeção contra anomalias
-além de já antecipar as conclusões mais salientes que serão,
-possivemente, certificadas pela análise.
+além de já antecipar os padrões mais salientes que serão, possivemente,
+certificadas pela análise.
 
+Os dados estão armazenados no objeto `PimentelEg6.2` do pacote
+[labestData]. Caso não tenha o pacote instalado, os dados podem ser
+lidos a partir da URL do repositório como mostra o código a seguir.
+
+```{r, eval = FALSE}
+# Lendo os dados do repositório.
+url <- paste0("https://raw.githubusercontent.com/pet-estatistica/",
+              "labestData/devel/data-raw/PimentelEg6.2.txt")
+da <- read_tsv(url) %>%
+    mutate_at(1:3, factor)
+attr(da, "spec") <- NULL
+str(da)
+```
+
+Se você possui o labestData instalado, então os dados já estão
+prontamente disponíveis, incluindo a documentação.
+
+```{r, eval = FALSE}
+# Acessa a documentação do objeto.
+help(PimentelEg6.2, package = "labestData", help_type = "html")
+```
 ```{r}
 # Cria objeto com a tabela de dados.
 da <- labestData::PimentelEg6.2
 str(da)
 ```
 
-No data frame tem-se apenas os fatores experimentais linha, coluna e
-variedade já declarados como fator (classe `factor`) e a variável
-resposta que é numérica. O código abaixo mostra uma visão com
-organização retangular dos dados que é reproduzida de forma gráfica
-também.
+No *data frame* tem-se os 25 registros do quadrado latino $5 \times
+5$. As variáveis são os fatores experimentais linha, coluna e variedade
+declarados como fator (classe `factor`) e a variável resposta que é
+numérica. O código abaixo mostra uma visão com organização retangular
+dos dados que é reproduzida de forma gráfica também, por isso essa saída
+foi omitida. O interessante é que essa formatação dos dados é capaz de
+revelar valores ausente, o que não acontece neste experimento.
 
 ```{r, eval = FALSE}
-# Consulta à documentação dos dados.
-# help(PimentelEg6.2, package = "labestData", help_type = "html")
-
-# Visões em fomato retangular.
+# Visões em fomato retangular dos tratamentos.
 da %>%
     select(linha, coluna, varied) %>%
     spread(key = coluna, value = varied)
+
+# Visões em fomato retangular da resposta.
 da %>%
     select(linha, coluna, prod) %>%
     spread(key = coluna, value = prod)
-
-# Tem registros perdidos?
-sum(is.na(da$prod))
 ```
 
-O diagrama de dispersão antecipa alguma diferença entre variedades de
-cana de açucar. Apesar de cores s símbolos mapearem os níveis de linha e
-coluna, não é tão imediato perceber algum efeito esperado deles. Além do
-mais, como são termos de blocagem, existe apenas o interesse de acomodar
-o seu efeito para haver redução do erro experimental.
+A análise gráfica exploratória permite extrair padrões marcantes dos
+dados. Por exemplo, o diagrama de dispersão antecipa alguma diferença
+entre variedades de cana-de-açúcar, com a variedade C sendo a mais
+produtiva.
+
+Apesar de cores e símbolos mapearem os níveis de linha e coluna, não é
+tão imediato reconhecer algum efeito sistemático deles. Além do mais,
+como são termos de blocagem, existe apenas o interesse de acomodar o seu
+efeito para haver redução do erro experimental.
 
 ```{r}
 # Análise exploratória.
@@ -96,34 +124,42 @@ ggplot(data = da,
                  geom = "line",
                  fun.y = "mean",
                  show.legend = FALSE)
+```
 
-gg1 <-
-    ggplot(data = da,
+Os gráficos de mosaico mostram uma representação bastante comum de
+experimentos agronômicos em quadrado latino. Para experimentos em campo,
+eles são como uma visão aérea do experimento. Neste caso as linhas
+poderiam blocar a declividade do terreno e as colunas a qualidade/lote
+dos caules de propagação.
+
+Apesar de terem um certo apelo visual, tais gráficos são subótimos para
+apreciar o efeito de variedades, sendo preferido o diagrama de
+dispersão. Eles foram incluídos por completitude didática ao
+exemplificar como confeccioná-los.
+
+```{r, fig.height = 4.0}
+# Gráficos de mosaico.
+gg0 <- ggplot(data = da,
            mapping = aes(x = coluna,
-                         y = linha,
-                         fill = varied)) +
-    geom_tile() +
-    geom_text(mapping = aes(label = sprintf("%s\n%d", varied, prod))) +
+                         y = linha)) +
     coord_equal()
 
-gg2 <-
-    ggplot(data = da,
-           mapping = aes(x = coluna,
-                         y = linha,
-                         fill = prod)) +
-    geom_tile() +
-    geom_text(mapping = aes(label = sprintf("%s\n%d", varied, prod))) +
-    scale_fill_distiller(palette = "Greens", direction = 1) +
-    coord_equal()
+gg1 <- gg0 +
+    geom_tile(mapping = aes(fill = varied)) +
+    geom_text(mapping = aes(label = sprintf("%s\n%d", varied, prod)))
+
+gg2 <- gg0 +
+    geom_tile(mapping = aes(fill = prod)) +
+    geom_text(mapping = aes(label = sprintf("%s\n%d", varied, prod)))
 
 gridExtra::grid.arrange(gg1, gg2, nrow = 1)
 ```
 
-Os gráficos de mosaico mostram uma representação bastante comum de
-experimentos agronômicos em quadrado latino. As linhas poderiam blocar,
-a topografia e colunas a qualidade/lote dos caules de propagação. No
-entanto, tais gráficos são subótimos para apreciar o efeito de
-variedades, sendo preferido o diagrama de dispersão.
+Vários gráficos podem ser confeccionados até mesmo para um experimento
+simples como este de apenas uma variável resposta. No fragmento de
+código a seguir tem-se uma representação gráfica dos totais de linha,
+coluna e variedade. Ainda que sejam gráficos ainda não vistos, a
+informação tem redundancia com os anteriores, então serão omitidos.
 
 ```{r, eval = FALSE}
 gridExtra::grid.arrange(
@@ -139,19 +175,24 @@ gridExtra::grid.arrange(
                nrow = 1)
 ```
 
-A variável resposta a ser análisada é a produção da cultura. É uma
-resposta contínua, contrala por vários fatores genéticos e
-ambientais. Dessa forma, como na grande maioria das respostas
-biométricas, considera como apropriada a distribuição Gaussiana (ou
+Em resumo, a análise exploratória não indicou nenhuma anomalia. Ainda
+pode-se antecipar existência de diferença entre variedades quanto a
+produção média. Feita essa vistoria inicial, agora será formulado o
+modelo para análise deste experimento.
+
+A variável resposta é a produção da cultura. É uma resposta contínua,
+influenciada por vários fatores genéticos e ambientais. Dessa forma,
+como na grande maioria das respostas biométricas, considera como
+apropriada, pelo menos inicialmente, a distribuição Gaussiana (ou
 Normal).
 
-Dito isso e considerando o delineamento utilizado, o modelo estatístico
-para esse experimento é
+Dito isso e considerando o delineamento utilizado neste experimento, o
+modelo estatístico para essa variável resposta é
 $$
   \begin{aligned}
     Y_{ijk} &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \\
-    \mu_{ijk} &= \mu + \alpha_i + \beta_j + \tau_k \\
-    \sigma^2 \propto 1
+    \mu_{ijk} &= \mu + \alpha_i + \theta_j + \tau_k \\
+    \sigma^2 &\propto 1
   \end{aligned}
 $$
 em que:
@@ -159,24 +200,36 @@ em que:
   * $Y_{ijk}$ é a resposta observada na condição experimental $ijk$ para
     a qual se assume a distribuição Normal,
   * $\mu_{ijk}$ é a média para as observações na condição $ijk$ e
-  * $\sigma^2$ é a variância das observações ao redor desse valor médio,
-  * $\alpha_i$ é o efeito da linha $i$ (fator qualitativo),
-  * $\beta_j$ é o efeito da coluna $j$ (fator qualitativo) e
-  * $\tau_j$ o efeito da variedade $j$ (fator qualitativo) sobre a
-    variável resposta $Y$ para a qual deseja-se testar a hipótese nula
-    $H_0: \tau_j = 0$ para todo $j$,
-  * $\mu$ é a média de $Y$ na ausência do efeito dos termos mencionados.
-
-Neste caso, todos os fatores são qualitativos. Por outro lado, se o
-fator experimental de interesse fosse quatitativo, como adubação,
-poderia-se acomodar o efeito por meio de um polinômio, um modelo não
-linear ou qualquer outra função de uma variável real. Geralmente,
-procura-se obter uma descrição que utilize menos parâmetros que o número
-de níveis do fator.
-
-O código abaixo declara o modelo acima equacionado para que a função
-`lm()` estime os parâmetros -- pelo método de mínimos quadrados -- e
-compute as demais saídas relevantes.
+  * $\sigma^2$ é a variância das observações ao redor desse valor médio
+    para qual considera-se ser constante,
+  * $\alpha_i$ é o parâmetro que acomoda o efeito da linha $i$ (fator
+    qualitativo),
+  * $\gamma_j$ é o parâmetro que acomoda o efeito efeito da coluna $j$
+    (fator qualitativo) e
+  * $\tau_k$ é o parâmetro que acomoda o efeito da variedade $k$ (fator
+    qualitativo) sobre a média da variável resposta $Y$ para a qual
+    deseja-se testar a hipótese nula $H_0: \tau_k = 0$ para todo $k$,
+  * $\mu$ é a média de $Y$ na ausência do efeito dos termos
+    mencionados. No caso de se considerar a restrição tipo tratamento em
+    que $\alpha_1 = \gamma_1 = \tau_1 = 0$, $\mu$ corresponde à média da
+    condição experimental de referência $111$, ou seja, $\mu = \mu_{111}$.
+
+Neste experimento todos os fatores são qualitativos. Assume-se que os
+efeitos de linha, coluna e variedade atuam de forma aditiva na
+média. Todavia, essa suposição pode não ser atendida e isso pode ser
+investigado na análise de resíduos. Ela também ser examinada pelo Teste
+da Não Aditividade de Tukey, visto no capítulo anterior, mas usando os
+fatores experimentais aos pares.
+
+Para o caso do fator experimental de interesse for quatitativo, como
+adubação, poderia-se acomodar seu efeito por meio de um polinômio, um
+modelo não linear ou qualquer outra função de uma variável
+real. Geralmente, procura-se obter uma descrição que utilize menos
+parâmetros que o número de níveis do fator.
+
+O código abaixo declara o modelo especificado acima para que a função
+`lm()` estime os parâmetros -- pelo método de mínimos quadrados
+ordinários -- e compute as demais saídas relevantes.
 
 ```{r, fig.height = 7}
 # Ajuste o modelo aos dados.
@@ -186,33 +239,145 @@ m0 <- lm(prod ~ linha + coluna + varied, data = da)
 par(mfrow = c(2, 2))
 plot(m0)
 layout(1)
+```
+
+De forma geral, os gráficos não fornecem suspeitas sobre o atendimento
+dos pressupostos.
+
+O gráfico *Residuals vs Fitted* exibe os resíduos crus $\hat{e} =
+y_{ijk} - \hat{y}_{ijk}$ contra os valores ajustados $\hat{y}_{ijk}$ que
+são as médias ajustadas pelo modelo para cada condição experimental. Não
+existe padrão não aleatório neste gráfico. Os resíduos ocupam as duas
+metades superior inferior na mesma proporção em todas as faixas de
+valores do eixo das coordenadas. Portanto, conclui-se que a suposição de
+especificação correta do modedo para a média não é violada.
+
+O gráfico *Scale-Location* coloca os a raíz quadrada do valor absolutos
+dos resíduos padronizados contra os valores ajustados. Os resíduos
+padronizados são os resíduos crus divididos pela desvio-padrão do erro
+experimental, $\hat{r} = \hat{e}/\sqrt{\hat{\sigma}^2}$. O gráfico não
+mostra relação de dependência entre a dispersão e os valores
+ajustados. Isso quer dizer que não existe evidência a favor de uma
+relação não nula entre dispersão e média. Logo, conclui-se que a
+suposição de homocedasticidade não é violada.
+
+O gráfico *Normal Q-Q* apresenta os resíduos padronizados em função dos
+quantis teóricos da distribuição normal padrão. Os pontos distribuem-se
+ao redor da linha diagonal não exibindo padrões que rejeitem a suposição
+de normalidade dos erros experimentais. Então, conclui-se que não a
+suposição de normalidade não é violada.
+
+O último gráfico mostra as alavancagens. No entanto, os valores de
+alavancagem são todos iguais, o que decorre do experimento ter apenas
+fatores qualitativos e ser balanceado. O exame desse gráfico é relevante
+para experimentos com fatores quantitativos.
+
+```{r, eval = FALSE, include = FALSE}
+# Sobre o teste da aditividade.
+
+# Usando as médias marginais. ---------------
+m0a <- lm(prod ~ linha + varied, da)
+anova(m0a)
+
+# Médias ajustadas centradas na média geral.
+a1 <- emmeans(m0a, ~linha) %>%
+    as.data.frame() %>%
+    select(1:2) %>%
+    mutate(z_1 = emmean - mean(db$prod),
+           emmean = NULL)
+a2 <- emmeans(m0a, ~varied) %>%
+    as.data.frame() %>%
+    select(1:2) %>%
+    mutate(z_2 = emmean - mean(db$prod),
+           emmean = NULL)
+a3 <- inner_join(db, a1) %>%
+    inner_join(a2) %>%
+    mutate(z = z_1 * z_2)
+
+m0a <- lm(prod ~ linha + varied, data = a3)
+m1a <- lm(prod ~ linha + varied + z, data = a3)
+anova(m0a, m1a)
+
+# Usando valores ajustados. -----------------
+m0b <- lm(prod ~ 1, a3)
+m1b <- lm(prod ~ linha, a3)
+m2b <- lm(prod ~ varied, a3)
+a3$zx <- (fitted(m1b) - fitted(m0b)) *
+    (fitted(m2b) - fitted(m0b))
+m3b <- lm(prod ~ linha + varied + zx, a3)
+anova(m3b, m1a)
+cbind(coef(m3b), coef(m1a))
+
+a3 %>%
+    select(z, zx)
+
+# Usando a função pronta para isso. ---------
+with(a3, nonadditivity(prod,
+                       linha,
+                       varied,
+                       df = df.residual(m0a),
+                       MSerror = deviance(m0a)/df.residual(m0a)))
+```
+
+Como não existem evidências sobre a violação dos pressupostos, dá-se
+continuidade na análise sem precisar de reparação no modelo. Com isso,
+examina-se o quadro de análise de variância e o quadro de resumo do
+modelo.
 
+```{r}
 # Quadro de análise de variância.
 # anova(m0)
 Anova(m0)
 
 # Estimativas dos efeitos e medidas de ajuste.
 summary(m0)
-```
 
-TODO: descrição do gráfico de diagnóstico.
+summary(m0)$sigma/summary(m0)$coefficients[1, 2]
+
+X <- model.matrix(m0)
+MASS::fractions(solve(crossprod(X))[1, 1])
+
+2 * summary(m0)$sigma/sqrt(2 * 5)
+
+```
 
 O quadro de análise de variância contém evidência a favor para rejeição
-da hipótese nula de igualdade entre as variedades de cana-de-açucar. Os
+da hipótese nula de igualdade entre as variedades de cana-de-açúcar. Os
 gráficos na fase exploratória já antecipavam isso.
 
 O quadro de resumo do ajuste reporta a estimativa do erro padrão
-residual `r sprintf("%0.2f", summary(m0)$sigma)` e que os efeitos
-estimados tem o mesmo erro-padrão que vale $2
-\hat{\sigma}/\sqrt{2r}$. Lembre-se que esses efeitos são diferenças com
-relação ao um ponto experimental de referência que é a varidade A na
-primeira linha e primeira coluna. Na realidade, esse ponto experimental
-sequer ocorre no experimento mas isso não impede que sua média seja
-estimada a partir do modelo que foi assumido.
+residual `r sprintf("%0.2f", summary(m0)$sigma)`. Os efeitos estimados
+tem o mesmo erro-padrão que vale $2 \hat{\sigma}/\sqrt{2r}$. Lembre-se
+que esses efeitos são diferenças com relação ao um ponto experimental de
+referência que é a varidade A na primeira linha e primeira coluna, ou
+seja, são erros padrões para diferenças do tipo $\mu_{11k} -
+\mu_{111}$. Ou seja, em que apenas um dos índices $ijk$ tem valor
+diferente de 1. Na realidade, o ponto experimental $111$ sequer ocorre
+no experimento mas isso não impede que sua média seja estimada a partir
+do modelo que foi assumido e seja considerada o valor de referência.
 
 Como as linhas e colunas são fatores experimentais de blocagem,
 prefere-se obter as médias marginais das variedades que são obtidas e a
-seguir.
+seguir. As médias marginais ajustadas são uma função linear das
+estimativas na qual o efeito dos fatores complementares são nivelados
+pela média. Assim, a média marginal para uma variedade $k$ é obtida por
+$$
+\begin{align*}
+  \hat{\mu}_{..k} &=
+    \frac{1}{5 \cdot 5} \sum_{i} \sum_{j} \hat{\mu}_{ijk} \\
+  \hat{\mu}_{..k} &= \hat{\mu} +
+    \frac{1}{5} \sum_{i} \hat{\alpha}_{i} +
+    \frac{1}{5} \sum_{j} \hat{\theta}_{i} +
+    \hat{\tau}_k,
+\end{align*}
+$$
+em que o 5 nos denominadores deve-se ao fato de ser um quadrado latino
+$5 \times 5$.
+
+A função `emmeans()` calcula tais médias com erro-padrão e intervalo de
+confiança embutidos. O código abaixo extrai a matriz usada para o
+cálculo das médias. Note que os pesos usados para nivelar as estimativas
+de linha e coluna é 1/5, conforme equacionado.
 
 ```{r}
 # Médias marginais ajustadas.
@@ -229,18 +394,23 @@ rownames(L) <- grid[[1]]
 
 # Para entender como se calcula as médias marginais.
 t(MASS::fractions(L))
+
+# Determina as médias ajustadas e os erros-padrões.
+(L %*% coef(m0))[, 1]
+sqrt(diag(L %*% vcov(m0) %*% t(L)))
 ```
 
 Conforme discutido no capítulo anterior, o nome médias ajustadas decorre
 do fato de que para obter as médias marginais em relação a um fator,
-faz-se a média sobre as estimativas dos efeitos dos demais fatores. Por
-isso que a matriz `L` tem entradas 1/5 (são 5 linhas e colunas) para os
-coeficientes que acomodam o efeito de linhas e colunas.
+faz-se a média para nivelar sobre as estimativas dos efeitos dos demais
+fatores. Por isso que a matriz `L` tem entradas 1/5 (são 5 linhas e
+colunas) para os coeficientes que acomodam o efeito de linhas e colunas.
 
 A partir da matriz `L` pode-se avaliar os contrastes entre médias. A
-função `constrast()` já faz isso. A opção `pairwise` indica para que
-sejam calculados todos os contrastes entre pares de médias como se faz
-no teste HSD de Tukey ou SNK.
+função `constrast()` já faz isso a partir do conteúdo do objeto
+retornado pela `emmeans()`. A opção `pairwise` indica para que sejam
+calculados todos os contrastes entre pares de médias como se faz no
+teste HSD de Tukey ou SNK.
 
 ```{r}
 # Obtém os constrastes par a par (pairwise). P-valores são ajustados
@@ -258,9 +428,10 @@ pwpp(emm) +
 O pacote `emmeans`, portanto, já retorna os contrastes e até a
 representação compacta de letras (cld) que é com números. Antes desse
 pacote, isso era feito de forma mais longa usando os recursos do pacote
-`multcomp`, conforme pode-se ver a seguir.
+`multcomp`, conforme pode-se ver a seguir. Pelo que foi exposto, não
+existe mistério na obtenção das médias ajustadas.
 
-```{r}
+```{r, eval = FALSE}
 # Comparações múltiplas via contrastes de Tukey (pairwise).
 # Método de correção de p-valor: single-step.
 tk_sgl <- summary(glht(m0, linfct = mcp(varied = "Tukey")),
@@ -271,26 +442,45 @@ tk_sgl
 cld(tk_sgl, decreasing = TRUE)
 ```
 
-Essa saída exibe os mesmos contrastes e o resumo compacto de
-letras. Estes resultados podem ser combinados em data frame para
+Essa saída (omitida) exibe os mesmos contrastes e o resumo compacto de
+letras. Estes resultados podem ser combinados em *data frame* para
 confecção de gráficos e tabelas que será mostrado na próxima seção.
 
 Deve-se tomar cuidado ao usar a opção com `mcp()` quando o modelo
 contiver termos de interação. Nos capítulos à frente isso será discutido
-em profundidade.
+em mais detalhes.
 
 O teste HSD de Tukey tem uma grande preferência entre praticantes de
 análise de experimentos. No entanto, não será observada nenhuma
 diferença prática entre as abordagem anteriores que i) avaliam a
-significância dos contrastes por um teste **t* seguido da correção dos
+significância dos contrastes por um teste *t* seguido da correção dos
 p-valores (ou correlação dos níveis de significância) para a
 multiplicidade de hipóteses e ii) da aplicação do teste HSD de Tukey que
 é baseado na distribuição da amplitude total studentizada para corrigir
 para a multiplicidade.
 
-O teste HSD de Tukey está implementado no pacote `agricolae`. O código
-abaixo aplica o teste e manuseia os resultados para apresentá-los de
-forma gráfica. A incerteza sobre as médias é representada com o
+O pacote `stats` tem a função `TukeyHSD()` para a aplicação do teste de
+HSD de Tukey. Os intervalos de confiança e p-valores são baseados na
+distribuição da amplitude total studentizada empregada no teste HSD de
+Tukey. A documentação orienta usá-la com cautela para casos de
+experimentos desbalanceados. Ela não permite aplicar o teste para o
+desdobramento de interações, então o seu uso é limitado aos modelos de
+fatores com efeitos aditivos.
+
+```{r}
+# Aplica o teste de HSD de Tukey do pacote `stats`.
+tHSD <- TukeyHSD(aov(formula(m0), data = m0$model),
+                 which = "varied",
+                 ordered = FALSE)
+tHSD
+
+# Método gráfico que exibe IC para os contrastes.
+plot(tHSD)
+```
+
+O teste HSD de Tukey também está implementado no pacote `agricolae`. O
+código abaixo aplica o teste e manuseia os resultados para apresentá-los
+de forma gráfica. A incerteza sobre as médias é representada com o
 intervalos de confiança (95%) baseados no modelo ajustado.
 
 ```{r}
@@ -301,13 +491,15 @@ tk_hsd <- HSD.test(m0, trt = "varied")
 tk_hsd$parameters
 tk_hsd$statistics
 
+# Converte resultados para data frame e junta com os IC.
 results <- tk_hsd$groups %>%
-    rownames_to_column("varied")
-
-results <- inner_join(results,
-                      as.data.frame(emm))
+    rownames_to_column("varied") %>%
+    inner_join(as.data.frame(emm))
 results
 
+cap <- paste("Pares de médias seguidas de mesma letra não diferem",
+             "entre si ao nível de significância global de 5%.")
+
 # Gráfico com estimativas, IC e texto.
 ggplot(data = results,
        mapping = aes(x = reorder(varied, prod), y = prod)) +
@@ -317,20 +509,28 @@ ggplot(data = results,
     geom_label(mapping = aes(label = sprintf("%0.1f %s", prod, groups)),
                nudge_y = 5,
                vjust = 0) +
-    labs(x = "Variedade", y = "Produção")
+    labs(x = "Variedade de cana-de-açúcar",
+         y = "Produção",
+         title = "Produção de variedades de cana-de-açúcar",
+         caption = cap)
 ```
 
-O teste de Scott-Knott é muito empregado, principalmente por
-pesquisadores de agrárias. A preferência parece ser pela "não
-abiguidade" nos resultados do teste que não apresenta resultados do tipo
-"ab", por exemplo. No entanto, muito dos usuários talvez não saibam que
-o teste de Scott-Knott **NÃO** é um teste de comparação múltipla de
-pares de médias mas sim um teste de agrupamento (divisivo) de médias. É
-por isso que ainda é muito comum as pessoas cometerem o equívoco de
-interpretarem o resultado do teste de Scott-Knott como se fosse um teste
-de pares de médias.
+Além do teste de Tukey, o teste de Scott-Knott é muito empregado,
+principalmente por pesquisadores de Agrárias. A preferência parece ser
+pela "não abiguidade" nos resultados do teste que não apresenta
+resultados do tipo "ab", por exemplo. No entanto, muito dos usuários
+talvez não saibam que o teste de Scott-Knott **NÃO** é um teste de
+comparação múltipla de pares de médias mas sim um teste de agrupamento
+(divisivo) de médias. É por isso que ainda é muito comum as pessoas
+cometerem o equívoco de interpretarem o resultado do teste de
+Scott-Knott como se fosse um teste de pares de médias.
+
+O pacote `ScottKnott` implementa o teste de Scott-Knott para modelos de
+classe `aov()`. O uso dos recursos do pacote é bastante simples, como
+pode ser visto no código abaixo.
 
 ```{r}
+# Carrega o pacote.
 library(ScottKnott)
 
 # Ajustando o modelo com a `aov()` para passar para a `SK()`.
@@ -347,12 +547,16 @@ formado pelas demais variáveis. Não podemos cometer o erro de
 interpretar que isso implica na rejeição da hipótese de igualdade entre
 as médias da variedade C e A, por exemplo.
 
-O método `summary()` aplicado ao objeto vindo da `SK()` é um data frame
-que contém as médias dos tratamentos e as letras para indicar os grupos
-formados. Dessa forma é fácil construir gráficos a partir da manipulação
-desse objeto.
+O método `summary()` aplicado ao objeto vindo da `SK()` é um *data
+frame* que contém as médias dos tratamentos e as letras para indicar os
+grupos formados. Dessa forma é fácil construir gráficos a partir da
+manipulação desse objeto.
 
-## Variável resposta de contagem
+## Resposta contínua e fator quantitativo
+
+TODO
+
+## Resposta de contagem
 
 Respostas do tipo contagem quantificam a ocorrência de um evento sob um
 dado contexto, como por exemplo:
@@ -364,6 +568,7 @@ dado contexto, como por exemplo:
     (contexto).
   * O número de ovos (evento) em ninhos de galinhas caipiras em uma
     fazenda (contexto).
+  * O número de perfilhos (evento) por touceira de capim (contexto).
 
 Note nestes exemplos que o contexto não permite saber qual o total
 teórico máximo para o número de eventos. Se por outro lado,
@@ -376,12 +581,16 @@ analisadas na próxima seção.
 dessa seção referem-se ao número de colmos atacados (evento) por uma
 praga em plantas de arroz. A documentação dos dados em @zimmermann não
 esclarece precisamente qual é o contexto, então podemos imaginar que
-refere-se a todas as parcela experimental.
+refere-se a todas as parcela experimental. Os dados do objeto
+`ZimmermannTb5.15` referem-se ao número de perfilhos em experimento com
+a cultura do arroz feito em DQL.
 
-```{r}
+```{r, eval = FALSE}
 # Documentação dos dados.
-# help(ZimmermannTb16.10, package = "labestData", help_type = "html")
+help(ZimmermannTb16.10, package = "labestData", help_type = "html")
+```
 
+```{r}
 # Cria objeto com os dados.
 da <- labestData::ZimmermannTb16.10
 str(da)
@@ -397,51 +606,27 @@ ggplot(data = da,
                  geom = "line",
                  fun.y = "mean",
                  show.legend = FALSE)
-
-gg1 <-
-    ggplot(data = da,
-           mapping = aes(x = coluna,
-                         y = linha,
-                         fill = trat)) +
-    geom_tile(show.legend = FALSE) +
-    geom_text(mapping = aes(label = sprintf("[%s]\n%d",
-                                            trat,
-                                            colmos)),
-              size = 3) +
-    coord_equal()
-
-gg2 <-
-    ggplot(data = da,
-           mapping = aes(x = coluna,
-                         y = linha,
-                         fill = colmos)) +
-    geom_tile(show.legend = FALSE) +
-    geom_text(mapping = aes(label = sprintf("[%s]\n%d",
-                                            trat,
-                                            colmos)),
-              size = 3) +
-    scale_fill_distiller(palette = "Greens", direction = 1) +
-    coord_equal()
-
-gridExtra::grid.arrange(gg1, gg2, nrow = 1)
 ```
 
-Os gráficos usados são os mesmos considerados na seção anterior. O
-primeiro deles, o diagrama de dispersão, fornece o diagnóstico mais
-interessante que é sobre o efeito dos tratamentos para controle da
-praga.
+Os gráficos a serem usados podem ser reaproveitados da seção anterior. O
+diagrama de dispersão fornece o diagnóstico mais interessante que é
+sobre o efeito dos tratamentos para controle da praga. Observa-se que a
+distribuição é mais assimétrica se comparado com os dados da sessão
+anterior. Embora os valores extremos para cima possam preocupar
+mostrando-se atípicos, pode ser que eles sejam ocasionados pelo efeito
+de linha e/ou coluna. Portanto, um julgamento mais assertivo pode ser
+feito na análise dos resíduos.
 
 Apenas com o propósito didático, será mostradado a análise primeiro
 assumindo distribuição Gaussiana para se discutir algumas coisas. O
 modelo estatístico é o mesmo apresentado na seção anterior. É comum
 dados de contagem apresentarem alguma relação positiva de
-média-variância o que significa que a suposição de homocedasticidade não
-é verificada.
+média-variância o que significa que a suposição de homocedasticidade é
+violada.
 
 ```{r, fig.height = 7}
 # Ajusta o modelo aos dados MAS assumindo uma distribuição normal.
-m0 <- lm(colmos ~ linha + coluna + trat,
-         data = da)
+m0 <- lm(colmos ~ linha + coluna + trat, data = da)
 
 # Diganóstico dos resíduos.
 par(mfrow = c(2, 2))
@@ -449,10 +634,21 @@ plot(m0)
 layout(1)
 ```
 
-O gráfico Scale-Location indica a relação média-variância comentada. A
-normalidade, no então, não mostra afastamento. A transformação Box-Cox
-indica que usar o logaritmo da resposta dá mais adesão aos pressupostos
-do modelo Gaussiano.
+O gráfico *Residuals vs Fitted* exibe um padrão de megafone, que é
+salientado pelo gráfico *Scale-Location*. Ambos levantam a suspeita
+sobre uma relação média-variância não nula. Não é proveitoso julgar a
+normalidade diante de uma suspeita sobre a violação de outros
+pressupostos, no entanto, não se vê grave afastamento dessa suposição.
+
+Uma solução possível é transformar
+
+Para diminuir a violação da suposição de variância constante pode-se
+tentar aplicar uma transformação na variável resposta. Busca-se aqui
+então uma transformação que estabilize a relação média-variância.  A
+transformação Box-Cox indica que usar o logaritmo da resposta, ou seja,
+$y_t = \log(y)$, dá mais adesão aos pressupostos do modelo
+Gaussiano. Dessa forma, a análise poderia ser feita com a variável
+resposta transformada assumindo-se distribuição Gaussiana.
 
 ```{r}
 # Indica o valor de lambda para usar na transformação Box-Cox.
@@ -460,36 +656,46 @@ MASS::boxcox(m0)
 abline(v = 0, col = "red")
 ```
 
-TODO: discutir a tranformação sugerida e como isso irá mitigar. Qual a
-transformação estabilizadora da variância para dados Poisson?
+A literatura mostra que a função estabilizadora da variância para uma
+variável que segue distribuição Poisson, cuja relação média-variância é
+$\text{Var}(Y) = \text{E}(y)$, é $Y_t = \sqrt{Y}$
+([aqui](https://www.sciencedirect.com/science/article/pii/S0167715209001473)). Para
+uma variável aleatória com relação média-variância dada por
+$\text{Var}(Y) \propto (\text{E}(Y))^2$, a tranformação estabilizadora é
+$Y_t = \log(Y)$
+([aqui](https://en.wikipedia.org/wiki/Variance-stabilizing_transformation)).
 
 Para demonstrar o emprego de uma distribuição mais apropriada para dados
 de contagem, será usado a Poisson, ou a especificação Quasi-Poisson. A
-diferença dessas especificações é muito relevante:
+diferença dessas especificações é pontual mas muito relevante:
 
   * A Poisson é uma especificação completa, ou seja, baseada em uma
     distribuição de probabilidade especificada por completo o que
     permite escrever a função de verossimilhança para se fazer a
-    estimação. Como a distribuição tem apenas um parâmetro, a variância
-    da distribuição é função da média pela conhecida relação
+    estimação. Como se sabe, a distribuição tem apenas um parâmetro e a
+    variância da distribuição é função da média pela conhecida relação
     $\text{Var}(Y) = \text{E}(Y) = \lambda$ se $Y \sim
     \text{Poisson}(\lambda)$. Essa característica representa uma
     suposição limitante para cenários reais que geralmente apresentam
     uma variância amostral superior à média (superdispersão).
   * A Quasi-Poisson é uma especificação incompleta, ou seja, não é um
     modelo probabilístico completo. *Quasi* é um termo em latim para
-    quase, como se fosse. Ela imita a Poisson na especificação dos
-    momentos mas introduz um parâmetro de dispersão $\phi$ que
-    estimado. A estimação é feita pelo mesmo algorítmo usado para a
+    quase, como se fosse. Ela imita a Poisson na especificação dos dois
+    primeiros momentos mas introduz um parâmetro de dispersão $\phi$ que
+    é estimado. A estimação é feita pelo mesmo algorítmo usado para a
     Poisson, mas os erros padrões são corrigidos pela estimativa do
-    parâmetro de dispersão.
+    parâmetro de dispersão. Dessa forma, a Quasi-Poisson é uma
+    especificação menos restritiva pois tem um parâmetro dedicado para a
+    dispersão.
 
-O modelo estatístico para esse experimento é
+O modelo estatístico para a variável resposta de contagem considerando o
+delineamento de quadrado latino é
 $$
   \begin{aligned}
-    y_{ijk} &= \text{Quasi-Poisson}(\lambda_{ijk}, \phi)\\
-    g(\lambda_{ijk}) &= \mu + \alpha_i + \beta_j + \tau_k \\
-    \phi \propto 1
+    Y_{ijk} &= \text{Quasi-Poisson}(\lambda_{ijk}, \phi) \\
+    g(\lambda_{ijk}) &= \mu + \alpha_i + \theta_j + \tau_k \\
+    \phi &\propto 1 \\
+    \text{Var}(Y_{ijk}) &= \phi \text{E}(Y_{ijk})
   \end{aligned}
 $$
 em que:
@@ -499,17 +705,18 @@ em que:
   * $\lambda_{ijk}$ é a média para as observações na condição $ijk$ e
   * $\phi$ é o parâmetro de dispersão que acomoda a variação das
     observações ao redor desse valor médio. Com base na especificação
-    Quasi-Poisson, $\text{Var}(Y) = \phi\text{E}(Y)$.
+    Quasi-Poisson, tem-se um relação média-variância dada por
+    $\text{Var}(Y) = \phi \text{E}(Y)$ que é assumida pelo modelo.
   * $g(.)$ é a função de ligação. A função de ligação canônica para a
     Poisson é a logarítmo neperiano.
   * $\alpha_i$ é o efeito da linha $i$ (fator qualitativo),
-  * $\beta_j$ é o efeito da coluna $j$ (fator qualitativo) e
-  * $\tau_j$ o efeito da variedade $j$ (fator qualitativo) sobre a
+  * $\theta_j$ é o efeito da coluna $j$ (fator qualitativo) e
+  * $\tau_k$ o efeito da variedade $k$ (fator qualitativo) sobre a
     variável resposta $Y$ para a qual deseja-se testar a hipótese nula
-    $H_0: \tau_j = 0$ para todo $j$,
-  * $\mu$ é a média de $Y$ na ausência do efeito dos termos mencionados.
-
-TODO: qual a diferença da Poisson e da Quasi-Poisson.
+    $H_0: \tau_k = 0$ para todo $k$,
+  * $\mu$ é a média de $Y$ na ausência do efeito dos termos mencionados
+    mas representa a média na condição de referência, ou seja, $\mu =
+    \mu_{111}$.
 
 ```{r, fig.height = 7}
 # Análise com modelo linear generalizado.
@@ -521,11 +728,14 @@ m0 <- glm(colmos ~ linha + coluna + trat,
 par(mfrow = c(2, 2))
 plot(m0)
 layout(1)
+```
 
-# Quadro de análise de variância.
-# anova(m0, test = "F")
-# drop1(m0, test = "F")
-Anova(m0, test.statistic = "F")
+TODO discutir os resíduos.
+
+```{r}
+# Quadro de análise de deviance.
+# drop1(m0, test = "F") # Hipóteses marginais.
+anova(m0, test = "F") # Hipóteses sequenciais.
 
 # Estimativas dos efeitos e medidas de ajuste.
 summary(m0)
@@ -539,6 +749,8 @@ summary(m0)
 sum(residuals(m0, type = "pearson")^2)/df.residual(m0)
 ```
 
+TODO WALMES FIXME
+
 Pelo quadro de deviance, tem-se uma moderada evidência a favor da
 hipótese alternativa. A hipótese nula é de que não há efeito de
 tratamentos para o controle da praga. A estimativa do parâmetro de
@@ -949,42 +1161,11 @@ TODO: são tantas espeficicações possiveis. Qual é a mais adequada?
 Talvez nem haja diferença apreciável nos resultados. Vai do analista
 decidir.
 
-<!------------------------------------------- -->
-
-```{r, eval = FALSE, include = FALSE}
-library(labestData)
-
-labestDataView()
-
-str(ZimmermannTb16.10) # 6x6 colmos x posto.
-help(ZimmermannTb16.10,  h = "html") # 6x6 colmos x posto.
-
-str(ZimmermannTb5.11) # 11x11 de proporção.
-help(ZimmermannTb5.11, h = "html")
-
-str(ZimmermannTb5.15) # perfilhos.
-str(ZimmermannTb5.2) # 8x8.
-
-str(BarbinPg104) # aka PimentelEg6.2
-str(PimentelEg6.2)
-
-str(PimentelEx6.6.3) # Regressão.
-
-str(PimentelTb6.3.1) # Fatorial?
-
-str(DiasEg9.4)
-str(DiasEx9.6.10)
-
-str(PimentelTb9.3.1)
-str(StorckEg2.3.5)
-str(ZimmermannTb12.26)
-str(ZimmermannTb12.27)
-str(ZimmermannTb14.9)
-```
-
 ## Geração do delineamento
 
 ```{r}
 latin_square(6)
 ```
 TODO!
+
+[labestData]: https://github.com/pet-estatistica/labestData
-- 
GitLab