Improves last section of chapter.

04-dql.Rmd

@@ -746,14 +746,28 @@ experimentos planejados.
 
 O conjunto de dados que será análise registrou a proporção de hastes
 sobreviventes ao ataque de pragas em um experimento quadrado latino $11
-\times 11$ que avaliou o efeito de inseticidas para controle da
-praga. Os dados já são a proporção amostral de hastes sobreviventes,
-portanto, não se conhece o numerador e denominador que produziram tal
-razão. Em experimentos como esse é importante ter o registro do par de
-variáveis: total de hastes sobreviventes ($y$) e total de hastes sobre
-exposição ($n$). Com esse par pode-se declarar uma variável com
-distribuição que geralmente é considerado apropriado se assumidas certas
-suposições.
+\times 11$ que avaliou o efeito de inseticidas para controle da praga. A
+resposta é a proporção amostral de hastes sobreviventes, em algum
+momento calculada pelo quociente do número de hastes sobreviventes ($y$,
+numerador) pelo total de hastes avaliadas ($n$, denominador). O ideal é
+que se registre o par $y$ e $n$, mesmo que ele seja constante para todas
+as unidades experimentais pois o $n$ determina a previsão das
+estimativas de proporção. Com esse par pode-se declarar uma variável com
+distribuição Binomial que geralmente é considerada apropriada se
+assumidas certas suposições.
+
+O experimento tinha 4 unidades amostrais por cela experimental, o que
+perfaz um vetor de valores observados de tamanho $11 \times 11 \times 4
+= 484$. Para a análise, no entanto, será feita a agregação por unidade
+experimental usando a média das 4 unidades amostrais. Poderia-se fazer a
+análise sem agregração, mas isso demandaria um modelo de efeitos
+aleatórios para dissociar o erro entre unidades experimentais do erro
+entre unidades amostrais.
+
+O gráfico abaixo indica haver alguma efeito do inseticida sobre a
+proporção de hastes sobreviventes. A variabilidade é bem pronunciada, de
+forma que talvez seja detectável diferença apenas entre os tratamentos
+extremos.
 
 ```{r}
 # help(ZimmermannTb5.11, package = "labestData", help_type = "html")
@@ -762,12 +776,14 @@ suposições.
 da <- labestData::ZimmermannTb5.11
 str(da)
 
+db <- da %>%
+    group_by(linha, coluna, inset) %>%
+    summarise(prop = mean(prop))
+
 # Análise exploratória.
-ggplot(data = da,
+ggplot(data = db,
        mapping = aes(x = reorder(inset, prop),
-                     y = prop,
-                     color = linha,
-                     size = coluna)) +
+                     y = prop)) +
     geom_jitter(width = 0.15, show.legend = FALSE) +
     stat_summary(mapping = aes(group = 1),
                  geom = "line",
@@ -775,6 +791,10 @@ ggplot(data = da,
                  show.legend = FALSE)
 ```
 
+A análise que será feita assumirá a especificação Quasi-Binomial. Para
+isso, embora o $n$ seja desconhecido, será assumido que é igual para
+todas as observações.
+
 O modelo estatístico para esse experimento é
 $$
   \begin{aligned}
@@ -790,72 +810,72 @@ em que:
   * $\theta_{ijk} \in (0, 1)$ é a proporção de sobrevivência na condição
     $ijk$,
   * $n_{ijk}$ é o número total de hastes expostas na condição $ijk$ e
-  * $\phi > 0$ é 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-Binomial, $\text{Var}(Y) \propto \phi \theta (1 - \theta)$.
-  * $g(.)$ é a função de ligação. A função de ligação canônica para a
-    Quasi-Binomial é a logit.
+  * $\phi > 0$ é o parâmetro de dispersão que acomoda a variação não
+    explicada nas observações ao redor do valor central $\theta_{ijk}$.
+  * Com base na especificação Quasi-Binomial, tem-se que os dois
+    primeiros momentos são $\text{E}(Y_{ijk}) = n_{ijk} \theta_{ijk}$ e
+    $\text{Var}(Y_{ijk}) = \phi n_{ijk} \theta_{ijk} (1 -
+    \theta_{ijk})$.
+  * $g(.)$ é a função de ligação. A função de ligação canônica é a
+    logit.
   * $\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$.
 
-```{r}
+Antes de declarar o modelo Quasi-Binomial, será averiguado como o modelo
+Gaussiano tem os pressupostos não atendidos. Uma das características
+desse material é mostrar formas alternativas de análise sempre que
+possível para ampliar o conhecimento do leitor.
+
+```{r, results = "hide", fig.show = "hide"}
 # Ajusta o modelo aos dados. Resposta: proporção de sobreviventes.
-m0 <- lm(prop ~ linha + coluna + inset,
+m1 <- lm(prop ~ linha + coluna + inset,
          data = db)
 
-# Diganóstico dos resíduos.
-par(mfrow = c(1, 2))
-plot(m0, which = 2:3)
-layout(1)
-
-# ATTENTION: Box-Cox não resolve problemas de relação média-variância
-# decrescente.
-
 # Modelo o complementar. Resposta: proporção de não sobreviventes.
-m0 <- lm(c(1 - prop) ~ linha + coluna + inset,
+m2 <- lm(c(1 - prop) ~ linha + coluna + inset,
          data = db)
 
-# Diganóstico dos resíduos.
-par(mfrow = c(1, 2))
-plot(m0, which = 2:3)
-layout(1)
-
 # Indica lambda para aplicar transformação Box-Cox.
 MASS::boxcox(m0)
 abline(v = 1/3, col = "red")
 
+# Aplicando a transformação Box-Cox.
+m3 <- lm(c(1 - prop)^(1/3) ~ linha + coluna + inset,
+         data = db)
+
 # Transformação arco seno da raíz quadrada da proporção amostral. É a
 # transformação estabilizadora da variância quando a variável resposta
 # tem distribuição binomial.
 db$asinsqrt <- asin(sqrt(db$prop))
-m0 <- lm(asinsqrt ~ linha + coluna + inset,
+m4 <- lm(asinsqrt ~ linha + coluna + inset,
          data = db)
+```
 
+```{r, fig.height = 12}
 # Diganóstico dos resíduos.
-par(mfrow = c(1, 2))
-plot(m0, which = 2:3)
+par(mfrow = c(4, 2))
+plot(m1, which = 2:3)
+plot(m2, which = 2:3)
+plot(m3, which = 2:3)
+plot(m4, which = 2:3)
 layout(1)
-
-# DISCUSSION: Pressupostos melhor atendidos.
-
-# ATTENTION FIXME: não sabemos o denominador dessa binomial! E pode ser
-# que esse `n` seja variável em função da abundância de hastes em cada
-# unidade experimental.
 ```
 
+TODO: Acima foram declarados 4 diferentes especificações.
+
 TODO: comentar os resultados.
 
-```{r}
+```{r, eval = FALSE}
 # Análise com modelo linear generalizado.
 size <- 1
 m0 <- glm(cbind(size * prop, size - size * prop) ~
               linha + coluna + inset,
           data = db,
           family = quasibinomial)
-deviance(m0)
+m0
 
 # TIP: as inferências são invariantes ao valor do `size` para
 # quasibinomial. O parâmetro de dispersão absorve o `size`.
@@ -863,22 +883,27 @@ size <- 100
 m0 <- glm(cbind(size * prop, size - size * prop) ~ linha + coluna + inset,
           data = db,
           family = quasibinomial)
-deviance(m0)
+m0
 
 # Var(Y) \propto \phi * p * (1 - p).
 m0 <- glm(prop ~ linha + coluna + inset,
           data = db,
           family = quasi(link = "logit",
                          variance = "mu(1 - mu)"))
-deviance(m0)
+m0
+```
+
+Para facilitar a análise, será usada função de ligação identidade. Como
+todos os fatores são qualitativos e não será feito extrapolação
+IMPROVE. Mas atenção, isso não impede intervalos de confiança fora do
+espaço paramétrico.
 
+```{r}
 # Var(Y) \propto \phi * p * (1 - p).
 m0 <- glm(prop ~ linha + coluna + inset,
           data = db,
           family = quasi(link = "identity",
                          variance = "mu(1 - mu)"))
-deviance(m0)
-
 
 # Diganóstico dos resíduos.
 par(mfrow = c(1, 2))
@@ -891,28 +916,16 @@ layout(1)
 Anova(m0, test.statistic = "F")
 
 # Estimativas dos efeitos e medidas de ajuste.
-summary(m0)
+m0
 ```
 
 ```{r}
-# Comparações múltiplas.
-
 # Médias marginais ajustadas.
-# emmeans(m0, specs = ~inset, transform = "response")
-# emmeans(m0, specs = ~inset, transform = "unlink") # Igual.
-# emmeans(m0, specs = ~inset, transform = "mu")     # Igual.
-
 emm <- emmeans(m0, specs = ~inset)
 emm
 
-emm %>%
-    as.data.frame() %>%
-    mutate_at(c("emmean", "asymp.LCL", "asymp.UCL"),
-              m0$family$linkinv)
-
 # Retorna as médias ajustadas com teste de comparação múltipla aplicado.
-results <- cld(emm)
-results <- results %>%
+results <- cld(emm) %>%
     as.data.frame() %>%
     mutate_at(c(2, 5, 6), m0$family$linkinv) %>%
     mutate(.group = as_letters(.group))
@@ -924,13 +937,18 @@ ggplot(data = results,
     geom_point() +
     geom_errorbar(mapping = aes(ymin = asymp.LCL, ymax = asymp.UCL),
                   width = 0) +
-    geom_label(mapping = aes(label = sprintf("%0.2f %s", emmean, .group)),
+    geom_label(mapping = aes(label = sprintf("%0.3f %s", emmean, .group)),
                angle = 90,
                vjust = -0.4) +
+    # expand_limits(y = c(NA, 1)) +
     labs(x = "Tratamento",
          y = "Proporção de hastes sobreviventes")
 ```
 
+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}