Improves code in Latin Square Designs.

04-dql.Rmd

@@ -29,11 +29,11 @@ knitr::include_graphics("./img/quadrado-latino-drone.png")
 rm(list = objects())
 
 # Carrega pacotes.
-library(agricolae) # HSD.test()
-library(car)       # linearHypothesis(), Anova()
-library(multcomp)  # glht()
-library(emmeans)   # emmeans()
-library(tidyverse)
+library(agricolae) # HSD.test().
+library(car)       # linearHypothesis(), Anova().
+library(multcomp)  # glht().
+library(emmeans)   # emmeans().
+library(tidyverse) # Manipulação e visualização de dados.
 
 # Carrega funções.
 source("mpaer_functions.R")
@@ -46,8 +46,12 @@ da <- labestData::PimentelEg6.2
 str(da)
 
 # Vendo o quadrado no plano.
-reshape2::dcast(da, linha ~ coluna, value.var = "varied")
-reshape2::dcast(da, linha ~ coluna, value.var = "prod")
+da %>%
+    select(linha, coluna, varied) %>%
+    spread(key = coluna, value = varied)
+da %>%
+    select(linha, coluna, prod) %>%
+    spread(key = coluna, value = prod)
 
 # Tem registros perdidos?
 sum(is.na(da$prod))
@@ -101,8 +105,9 @@ gridExtra::grid.arrange(
 O modelo estatístico para esse experimento é
 $$
   \begin{aligned}
-    Y|i, j, k &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \\
-    \mu_{ijk} &= \mu + \alpha_i + \beta_j + \tau_k
+    Y_{ijk} &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \\
+    \mu_{ijk} &= \mu + \alpha_i + \beta_j + \tau_k \\
+    \sigma^2 \propto 1
   \end{aligned}
 $$
 em que $\alpha_i$ é o efeito da linha $i$, $\beta_j$ é o efeito da
@@ -112,7 +117,18 @@ mencionados, $\mu_{ijk}$ é a média para as observações de procedência
 $ijk$ e $\sigma^2$ é a variância das observações ao redor desse valor
 médio.
 
+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 e compute as demais saídas relevantes.
+
 ```{r, fig.height = 7}
+# Ajuste o modelo aos dados.
 m0 <- lm(prod ~ linha + coluna + varied, data = da)
 
 # Diganóstico dos resíduos.
@@ -126,8 +142,26 @@ Anova(m0)
 
 # Estimativas dos efeitos e medidas de ajuste.
 summary(m0)
+
+sprintf("%0.2f", summary(m0)$sigma)/sqrt(4)
+2 * summary(m0)$sigma/sqrt(2 * 5)
 ```
 
+O quadro de análise de variância contém evidência para rejeição da
+hipótese nula de igualdade entre as variedades de batatinha. 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 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 não ocorre no experimento mas isso
+não impede que sua média seja estimada.
+
+Como as linhas e colunas são fatores experimentais de blocagem,
+prefere-se obter as médias marginais das variedades.
+
 ```{r}
 # Comparações múltiplas.
 
@@ -135,11 +169,26 @@ summary(m0)
 emm <- emmeans(m0, specs = ~varied)
 emm
 
+# Erro padrão.
+summary(m0)$sigma/sqrt(5)
+
 grid <- attr(emm, "grid")
 L <- attr(emm, "linfct")
 rownames(L) <- grid[[1]]
 
-MASS::fractions(L)
+# Para entender como se calcula as médias marginais.
+t(MASS::fractions(L))
+
+# Obtém os constrastes par a par (pairwise). P-valores são ajustados
+# pelo método de Tukey.
+contrast(emm,  method = "pairwise")
+
+# Resumo compacto com letras.
+cld(emm)
+
+# Representação gráfica.
+pwpp(emm) +
+    geom_vline(xintercept = 0.05, lty = 2)
 
 # Comparações múltiplas, contrastes de Tukey.
 # Método de correção de p-valor: single-step.
@@ -150,7 +199,7 @@ tk_sgl
 # Resumo compacto com letras.
 cld(tk_sgl, decreasing = TRUE)
 
-# Teste HSD de Tukey.
+# Teste HSD de Tukey, baseado na amplitude total studentizada.
 tk_hsd <- HSD.test(m0, trt = "varied")
 
 # Detalhes da aplicação do teste HSD de Tukey.
@@ -176,11 +225,41 @@ ggplot(data = results,
     labs(x = "Variedade", y = "Produção")
 ```
 
-TODO! Colocar teste de Scott Knott.
+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. É
+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.
+
+```{r}
+library(ScottKnott)
+
+# Ajustando o modelo com a `aov()` para passar para a `SK()`.
+m1 <- aov(formula(m0), data = m0$model)
+
+# Faz o teste de Scott-Knott.
+sk <- SK(m1, which = "varied")
+summary(sk)
+```
+
+Conforme o que foi dito, pelo teste de Scott-Knott, rejeita-se a
+hipótese nula de que a média da variedade C seja igual a média do grupo
+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.
 
 ## Variável resposta de contagem
 
 ```{r}
+# Cria objeto com os dados.
 da <- labestData::ZimmermannTb16.10
 str(da)
 
@@ -224,7 +303,13 @@ gg2 <-
 gridExtra::grid.arrange(gg1, gg2, nrow = 1)
 ```
 
+Apenas com o propósito didático, será mostradado a análise assumindo
+distribuição Gaussiana. É 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.
+
 ```{r, fig.height = 7}
+# Ajusta o modelo aos dados MAS assumindo uma distribuição normal.
 m0 <- lm(colmos ~ linha + coluna + trat,
          data = da)
 
@@ -234,11 +319,35 @@ 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.
+
 ```{r}
 MASS::boxcox(m0)
 abline(v = 0, col = "red")
 ```
 
+Para demonstrar o emprego de uma distribuição mais apropriada para dados
+de contagem, será usado a Poisson, ou a especificação Quasi-Poisson.
+
+O modelo estatístico para esse experimento é
+$$
+  \begin{aligned}
+    y_{ijk} &= \text{Quasi-Poisson}(\lambda_{ijk}, \phi)\\
+    g(\lambda_{ijk}) &= \mu + \alpha_i + \beta_j + \tau_k \\
+    \phi \propto 1
+  \end{aligned}
+$$
+em que $\lambda_{ijk}$ é a média para a condição experimental $ijk$,
+$\alpha_i$ é o efeito da linha $i$, $\beta_j$ é o efeito da coluna $j$ e
+$\tau_j$ o efeito dos tratamentos $j$ sobre a variável resposta $Y$,
+$\phi$ é o parâmetro de dispersão e $g(.)$ é a função de ligação. A
+função de ligação canônica para a Poisson é a logarítmo neperiano.
+
+TODO: qual a diferença da Poisson e da Quasi-Poisson.
+
 ```{r, fig.height = 7}
 # Análise com modelo linear generalizado.
 m0 <- glm(colmos ~ linha + coluna + trat,
@@ -260,24 +369,59 @@ summary(m0)
 ```
 
 ```{r}
-# Comparações múltiplas.
-
-# Médias marginais ajustadas.
-emmeans(m0, specs = ~trat, transform = "response")
-# emmeans(m0, specs = ~trat, transform = "unlink") # Igual.
-# emmeans(m0, specs = ~trat, transform = "mu")     # Igual.
-
+# "Médias" ajustadas na escala do preditor linear.
 emm <- emmeans(m0, specs = ~trat)
 emm
 
+# Passando a inversa da função ligação nas 3 colunas.
 emm %>%
     as.data.frame() %>%
-    mutate_at(c("emmean", "asymp.LCL", "asymp.UCL"),
-              m0$family$linkinv)
+    mutate_at(c(2, 5, 6), m0$family$linkinv)
 
-# DANGER FIXME: o que a `transform = "response"` está fazendo não é
-# passar o inverso da função de ligação. O que está sendo feito então?
+# O resultado já transformado.
+emmeans(m0, specs = ~trat, type = "response")
+
+# Médias amostrais (apenas para verificar).
+da %>%
+    group_by(trat) %>%
+    summarise_at(vars(colmos), mean)
+```
+
+```{r, include = FALSE}
+# QUESTION: qual média ajustada é mais correto de usar? Parece que a
+# segunda está mais próxima das médias amostrais. Mas como obter o erro
+# padrão?
 
+# Na escala da resposta.
+emmeans(m0, specs = ~trat, transform = "response")
+
+# QUESTION: como obtém-se o erro padrão?
+
+# Grid com todas as combinações 6 x 6 x 6 = 216.
+grid <- da %>%
+    select(linha, coluna, trat) %>%
+    complete(linha, coluna, trat)
+grid$mu <- predict(m0, newdata = grid, type = "response")
+
+grid <-
+    bind_cols(grid,
+              predict(m0,
+                      newdata = grid,
+                      type = "response",
+                      se.fit = TRUE)[1:2] %>%
+              as_tibble())
+str(grid)
+
+grid %>%
+    group_by(trat) %>%
+    summarise_at(vars(mu, fit, se.fit), mean)
+
+# DANGER! A diferença está na ordem das operações: i) passar a inversa
+# da função de ligação e ii) fazer a agregação ou ii) seguido de i).
+```
+
+```{r}
+# Desvendando o computo das médias ajustadas.
 grid <- attr(emm, "grid")
 L <- attr(emm, "linfct")
 rownames(L) <- grid[[1]]