From 154b35f2b808b8d497016b1af44e6cfd41d8ba03 Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Sat, 4 Jan 2020 17:50:57 -0300
Subject: [PATCH] Improves code of section.
---
04-dql.Rmd | 171 ++++++++++++++++++++++++++++++-----------------------
1 file changed, 98 insertions(+), 73 deletions(-)
diff --git a/04-dql.Rmd b/04-dql.Rmd
index a3a9359..0fcfde4 100644
--- a/04-dql.Rmd
+++ b/04-dql.Rmd
@@ -37,6 +37,7 @@ library(tidyverse) # Manipulação e visualização de dados.
# Carrega funções.
source("mpaer_functions.R")
+ls()
```
## Variável resposta contínua
@@ -491,11 +492,7 @@ $$
\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.
+em que:
* $Y_{ijk}$ é a resposta observada na condição experimental $ijk$ para
a qual se assume a especificação Quasi-Poisson,
@@ -672,7 +669,7 @@ contrastes -- que de fato são o que interessa -- isso é menos relevante
já que os efeitos de linha e coluna se cancelam na expressão de cada
contraste entre tratamentos.
-```{r, eval = FALSE}
+```{r}
# Grid com todas as combinações 6 x 6 x 6 = 216 e o valor predito na
# escala do preditor: \hat{\eta} = X \hat{\beta}.
grid <- da %>%
@@ -721,19 +718,47 @@ Conclusões: TODO!
```{r, eval = FALSE}
# Análise com modelo linear Gaussiano e relação \sigma^2 = \mu.
-m0 <- glm(colmos ~ linha + coluna + trat,
+m1 <- glm(colmos ~ linha + coluna + trat,
data = da,
family = quasi(variance = "mu"))
# Diganóstico dos resíduos.
par(mfrow = c(2, 2))
-plot(m0)
+plot(m2)
layout(1)
+
+# Análise com modelo linear Poisson e ligação indentidade.
+m2 <- glm(colmos ~ linha + coluna + trat,
+ data = da,
+ family = quasipoisson(link = "identity"))
+
+# Medidas de ajuste.
+c(deviance(m1), deviance(m2))
+
+# Estiativas dos parâmetros.
+cbind(coef(m1), coef(m2))
```
## Variável resposta de proporção amostral
+Assim como dados de contagem, dados de proporção são muito comuns em
+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.
+
```{r}
+# help(ZimmermannTb5.11, package = "labestData", help_type = "html")
+
+# Cria o objeto com os dados.
da <- labestData::ZimmermannTb5.11
str(da)
@@ -748,40 +773,36 @@ ggplot(data = da,
geom = "line",
fun.y = "mean",
show.legend = FALSE)
+```
-db <- da %>%
- group_by(linha, coluna, inset) %>%
- summarise_at("prop", "mean")
-
-gg1 <-
- ggplot(data = db,
- mapping = aes(x = coluna,
- y = linha,
- fill = inset)) +
- geom_tile(show.legend = FALSE) +
- geom_text(mapping = aes(label = sprintf("%s\n%0.2f",
- inset,
- prop)),
- size = 2) +
- coord_equal()
-
-gg2 <-
- ggplot(data = db,
- mapping = aes(x = coluna,
- y = linha,
- fill = prop)) +
- geom_tile(show.legend = FALSE) +
- geom_text(mapping = aes(label = sprintf("%s\n%0.2f",
- inset,
- prop)),
- size = 2) +
- scale_fill_distiller(palette = "Greens", direction = 1) +
- coord_equal()
+O modelo estatístico para esse experimento é
+$$
+ \begin{aligned}
+ Y_{ijk} &= \text{Quasi-Binomial}(\theta_{ijk}, n_{ijk}, \phi)\\
+ g(\theta_{ijk}) &= \mu + \alpha_i + \beta_j + \tau_k \\
+ \phi \propto 1
+ \end{aligned}
+$$
+em que:
-gridExtra::grid.arrange(gg1, gg2, nrow = 1)
-```
+ * $Y_{ijk}$ é a resposta observada na condição experimental $ijk$ para
+ a qual se assume a especificação Quasi-Binomial,
+ * $\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.
+ * $\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}
+# Ajusta o modelo aos dados. Resposta: proporção de sobreviventes.
m0 <- lm(prop ~ linha + coluna + inset,
data = db)
@@ -793,7 +814,7 @@ layout(1)
# ATTENTION: Box-Cox não resolve problemas de relação média-variância
# decrescente.
-# Modelo para o desfecho complementar.
+# Modelo o complementar. Resposta: proporção de não sobreviventes.
m0 <- lm(c(1 - prop) ~ linha + coluna + inset,
data = db)
@@ -802,10 +823,13 @@ 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")
-# Transformação arco seno da raíz quadrada da proporção amostral.
+# 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,
data = db)
@@ -820,13 +844,18 @@ layout(1)
# 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: comentar os resultados.
+
+```{r}
# 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)
# TIP: as inferências são invariantes ao valor do `size` para
# quasibinomial. O parâmetro de dispersão absorve o `size`.
@@ -834,8 +863,22 @@ size <- 100
m0 <- glm(cbind(size * prop, size - size * prop) ~ linha + coluna + inset,
data = db,
family = quasibinomial)
+deviance(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)
+
+# 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)
-# Var(Y|.) = \phi * p * (1 - p).
# Diganóstico dos resíduos.
par(mfrow = c(1, 2))
@@ -867,47 +910,29 @@ emm %>%
mutate_at(c("emmean", "asymp.LCL", "asymp.UCL"),
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?
-
-grid <- attr(emm, "grid")
-L <- attr(emm, "linfct")
-rownames(L) <- grid[[1]]
-
-MASS::fractions(t(L))
-
-# Comparações múltiplas, contrastes de Tukey.
-# Método de correção de p-valor: false discovery rate.
-tk_sgl <- summary(glht(m0, linfct = mcp(inset = "Tukey")),
- test = adjusted(type = "fdr"))
-tk_sgl
-
-# Resumo compacto com letras.
-cld(tk_sgl, decreasing = TRUE)
-
-results <- wzRfun::apmc(X = L,
- model = m0,
- focus = "inset",
- test = "fdr") %>%
- mutate_at(c("fit", "lwr", "upr"),
- m0$family$linkinv)
-results %>%
- arrange(desc(fit))
+# Retorna as médias ajustadas com teste de comparação múltipla aplicado.
+results <- cld(emm)
+results <- results %>%
+ as.data.frame() %>%
+ mutate_at(c(2, 5, 6), m0$family$linkinv) %>%
+ mutate(.group = as_letters(.group))
+results
# Gráfico com estimativas, IC e texto.
ggplot(data = results,
- mapping = aes(x = reorder(inset, fit), y = fit)) +
+ mapping = aes(x = reorder(inset, emmean), y = emmean)) +
geom_point() +
- geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
+ geom_errorbar(mapping = aes(ymin = asymp.LCL, ymax = asymp.UCL),
width = 0) +
- geom_text(mapping = aes(label = sprintf("%0.2f %s", fit, cld)),
- angle = 90,
- vjust = 0,
- nudge_x = -0.05) +
+ geom_label(mapping = aes(label = sprintf("%0.2f %s", emmean, .group)),
+ angle = 90,
+ vjust = -0.4) +
labs(x = "Tratamento",
y = "Proporção de hastes sobreviventes")
```
+<!------------------------------------------- -->
+
```{r, eval = FALSE, include = FALSE}
library(labestData)
--
GitLab