From 8b653fbebb8fce7f573d116e0b7d46cc0c16b143 Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Thu, 2 Jan 2020 17:29:52 -0300
Subject: [PATCH] Improves section with gaussian distribution.
---
04-dql.Rmd | 151 +++++++++++++++++++++++++++++++++++++++++++----------
1 file changed, 123 insertions(+), 28 deletions(-)
diff --git a/04-dql.Rmd b/04-dql.Rmd
index acdee53..68ce588 100644
--- a/04-dql.Rmd
+++ b/04-dql.Rmd
@@ -41,11 +41,31 @@ source("mpaer_functions.R")
## Variável resposta contínua
+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
+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.
+
```{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.
+
+```{r, eval = FALSE}
+# Consulta à documentação dos dados.
+# help(PimentelEg6.2, package = "labestData", help_type = "html")
-# Vendo o quadrado no plano.
+# Visões em fomato retangular.
da %>%
select(linha, coluna, varied) %>%
spread(key = coluna, value = varied)
@@ -55,7 +75,15 @@ da %>%
# 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.
+```{r}
# Análise exploratória.
ggplot(data = da,
mapping = aes(x = reorder(varied, prod),
@@ -88,7 +116,15 @@ gg2 <-
coord_equal()
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.
+
+```{r, eval = FALSE}
gridExtra::grid.arrange(
ggplot(data = da,
mapping = aes(x = linha, y = prod)) +
@@ -102,7 +138,14 @@ gridExtra::grid.arrange(
nrow = 1)
```
-O modelo estatístico para esse experimento é
+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
+Normal).
+
+Dito isso e considerando o delineamento utilizado, o modelo estatístico
+para esse experimento é
$$
\begin{aligned}
Y_{ijk} &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \\
@@ -110,12 +153,18 @@ $$
\sigma^2 \propto 1
\end{aligned}
$$
-em que $\alpha_i$ é o efeito da linha $i$, $\beta_j$ é o efeito da
-coluna $j$ 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 dos termos
-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.
+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,
@@ -125,7 +174,8 @@ 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.
+`lm()` estime os parâmetros -- pelo método de mínimos quadrados -- e
+compute as demais saídas relevantes.
```{r, fig.height = 7}
# Ajuste o modelo aos dados.
@@ -142,43 +192,56 @@ 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.
+TODO: descrição do gráfico de diagnóstico.
+
+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
+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.
+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.
Como as linhas e colunas são fatores experimentais de blocagem,
-prefere-se obter as médias marginais das variedades.
+prefere-se obter as médias marginais das variedades que são obtidas e a
+seguir.
```{r}
-# Comparações múltiplas.
-
# Médias marginais ajustadas.
emm <- emmeans(m0, specs = ~varied)
emm
-# Erro padrão.
+# Erro padrão das médias ajustadas.
summary(m0)$sigma/sqrt(5)
+# Extração dos pesos para obter as médias ajustadas.
grid <- attr(emm, "grid")
L <- attr(emm, "linfct")
rownames(L) <- grid[[1]]
# Para entender como se calcula as médias marginais.
t(MASS::fractions(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.
+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.
+
+```{r}
# Obtém os constrastes par a par (pairwise). P-valores são ajustados
# pelo método de Tukey.
contrast(emm, method = "pairwise")
@@ -189,8 +252,15 @@ cld(emm)
# Representação gráfica.
pwpp(emm) +
geom_vline(xintercept = 0.05, lty = 2)
+```
-# Comparações múltiplas, contrastes de Tukey.
+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.
+
+```{r}
+# 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")),
test = adjusted(type = "single-step"))
@@ -198,7 +268,31 @@ tk_sgl
# Resumo compacto com letras.
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
+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.
+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
+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
+intervalos de confiança (95%) baseados no modelo ajustado.
+
+```{r}
# Teste HSD de Tukey, baseado na amplitude total studentizada.
tk_hsd <- HSD.test(m0, trt = "varied")
@@ -231,8 +325,9 @@ 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.
+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.
```{r}
library(ScottKnott)
--
GitLab