From a170ccb6c6997895a0725bb7ce5096b8c1e53e7b Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Sat, 4 Jan 2020 21:10:16 -0300
Subject: [PATCH] Improves last section of chapter.
---
04-dql.Rmd | 144 ++++++++++++++++++++++++++++++-----------------------
1 file changed, 81 insertions(+), 63 deletions(-)
diff --git a/04-dql.Rmd b/04-dql.Rmd
index 0fcfde4..d62001e 100644
--- a/04-dql.Rmd
+++ b/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}
--
GitLab