diff --git a/04-dql.Rmd b/04-dql.Rmd
index 68ce588ced0c29360a368c74435f136fba98e876..a3a93590ac53e54eaf644d8efcf8541353c46730 100644
--- a/04-dql.Rmd
+++ b/04-dql.Rmd
@@ -353,7 +353,34 @@ desse objeto.
## Variável resposta de contagem
+Respostas do tipo contagem quantificam a ocorrência de um evento sob um
+dado contexto, como por exemplo:
+
+ * O número de gols (evento) de uma partida (contexto).
+ * O número de acidentes (evento) em uma rodovia por semana (contexto).
+ * O número de vagens (evento) de planta de soja (contexto).
+ * O número peixes capturados (capturados) em uma jogada de rede
+ (contexto).
+ * O número de ovos (evento) em ninhos de galinhas caipiras em uma
+ fazenda (contexto).
+
+Note nestes exemplos que o contexto não permite saber qual o total
+teórico máximo para o número de eventos. Se por outro lado,
+considerarmos o número de ovos quebrados (evento) em uma dúzia de ovos
+(contexto), sabemos que o máximo de ovos quebrados é 12. Respostas de
+contagem como essa são resultados de ensaios de Bernoulli e serão
+analisadas na próxima seção.
+
+É comum experimentos terem como resposta variáveis de contagem. Os dados
+dessa seção referem-se ao número de colmos atacados (evento) por uma
+praga em plantas de arroz. A documentação dos dados em @zimmermann não
+esclarece precisamente qual é o contexto, então podemos imaginar que
+refere-se a todas as parcela experimental.
+
```{r}
+# Documentação dos dados.
+# help(ZimmermannTb16.10, package = "labestData", help_type = "html")
+
# Cria objeto com os dados.
da <- labestData::ZimmermannTb16.10
str(da)
@@ -398,10 +425,17 @@ 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.
+Os gráficos usados são os mesmos considerados na seção anterior. O
+primeiro deles, o diagrama de dispersão, fornece o diagnóstico mais
+interessante que é sobre o efeito dos tratamentos para controle da
+praga.
+
+Apenas com o propósito didático, será mostradado a análise primeiro
+assumindo distribuição Gaussiana para se discutir algumas coisas. O
+modelo estatístico é o mesmo apresentado na seção anterior. É 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.
@@ -420,12 +454,34 @@ indica que usar o logaritmo da resposta dá mais adesão aos pressupostos
do modelo Gaussiano.
```{r}
+# Indica o valor de lambda para usar na transformação Box-Cox.
MASS::boxcox(m0)
abline(v = 0, col = "red")
```
+TODO: discutir a tranformação sugerida e como isso irá mitigar. Qual a
+transformação estabilizadora da variância para dados Poisson?
+
Para demonstrar o emprego de uma distribuição mais apropriada para dados
-de contagem, será usado a Poisson, ou a especificação Quasi-Poisson.
+de contagem, será usado a Poisson, ou a especificação Quasi-Poisson. A
+diferença dessas especificações é muito relevante:
+
+ * A Poisson é uma especificação completa, ou seja, baseada em uma
+ distribuição de probabilidade especificada por completo o que
+ permite escrever a função de verossimilhança para se fazer a
+ estimação. Como a distribuição tem apenas um parâmetro, a variância
+ da distribuição é função da média pela conhecida relação
+ $\text{Var}(Y) = \text{E}(Y) = \lambda$ se $Y \sim
+ \text{Poisson}(\lambda)$. Essa característica representa uma
+ suposição limitante para cenários reais que geralmente apresentam
+ uma variância amostral superior à média (superdispersão).
+ * A Quasi-Poisson é uma especificação incompleta, ou seja, não é um
+ modelo probabilístico completo. *Quasi* é um termo em latim para
+ quase, como se fosse. Ela imita a Poisson na especificação dos
+ momentos mas introduz um parâmetro de dispersão $\phi$ que
+ estimado. A estimação é feita pelo mesmo algorítmo usado para a
+ Poisson, mas os erros padrões são corrigidos pela estimativa do
+ parâmetro de dispersão.
O modelo estatístico para esse experimento é
$$
@@ -435,12 +491,27 @@ $$
\phi \propto 1
\end{aligned}
$$
-em que $\lambda_{ijk}$ é a média para a condição experimental $ijk$,
+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.
+ * $Y_{ijk}$ é a resposta observada na condição experimental $ijk$ para
+ a qual se assume a especificação Quasi-Poisson,
+ * $\lambda_{ijk}$ é a média para as observações na condição $ijk$ e
+ * $\phi$ é 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-Poisson, $\text{Var}(Y) = \phi\text{E}(Y)$.
+ * $g(.)$ é a função de ligação. A função de ligação canônica para a
+ Poisson é a logarítmo neperiano.
+ * $\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.
+
TODO: qual a diferença da Poisson e da Quasi-Poisson.
```{r, fig.height = 7}
@@ -461,68 +532,92 @@ Anova(m0, test.statistic = "F")
# Estimativas dos efeitos e medidas de ajuste.
summary(m0)
+
+# Deviance residual.
+# sum(residuals(m0, type = "deviance")^2)
+# deviance(m0)
+
+# Como estimar o parâmetro de dispersão.
+# summary(m0)$dispersion
+sum(residuals(m0, type = "pearson")^2)/df.residual(m0)
```
+Pelo quadro de deviance, tem-se uma moderada evidência a favor da
+hipótese alternativa. A hipótese nula é de que não há efeito de
+tratamentos para o controle da praga. A estimativa do parâmetro de
+dispersão `r sprintf("%0.2f", summary(m0)$dispersion)` indica uma
+variabilidade superior àquela esperada para uma variável Poisson. Isso é
+esperado visto que há sempre algum nível de perturbação devido a
+inúmeros fatores não controlados para inflar tal expectativa.
+
+Visto que tive-se indicações de efeito de tratamentos, será feito estudo
+via comparações mútiplas para detectar diferenças entre os tratamentos
+para controle da praga. Apesar de estar-se com outra classe de modelos,
+os métodos usados ainda são os mesmos.
+
```{r}
# "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.
+# Passando a inversa da função ligação.
emm %>%
as.data.frame() %>%
- mutate_at(c(2, 5, 6), m0$family$linkinv)
+ select(-(3:4)) %>%
+ mutate_at(c(2:4), m0$family$linkinv)
-# O resultado já transformado.
+# O resultado já vem transformado, inclusive o erro padrão.
emmeans(m0, specs = ~trat, type = "response")
+# O método delta é usado para obter o erro padrão na escala da resposta.
+exp(coef(m0)["(Intercept)"])
+predict(m0,
+ newdata = data.frame(linha = "1", coluna = "1", trat = "1"),
+ type = "response",
+ se.fit = TRUE) %>%
+ as.data.frame()
+car::deltaMethod(m0, "exp((Intercept))")
+
# 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).
-```
+O código abaixo faz as comparações múltiplas entre os
+tratamentos. IMPROVE!
```{r}
# Desvendando o computo das médias ajustadas.
grid <- attr(emm, "grid")
L <- attr(emm, "linfct")
rownames(L) <- grid[[1]]
-
MASS::fractions(t(L))
+# Retorna a avaliação de cada contraste.
+contrast(emm, method = "pairwise")
+
+# 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(trat, emmean), y = emmean)) +
+ 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)),
+ nudge_y = 0.25,
+ vjust = 0) +
+ labs(x = "Tratamento", y = "Número de colmos atacados")
+```
+
+```{r, eval = FALSE, include = FALSE}
# Comparações múltiplas, contrastes de Tukey.
# Método de correção de p-valor: single-step.
tk_sgl <- summary(glht(m0, linfct = mcp(trat = "Tukey")),
@@ -538,17 +633,102 @@ results <- wzRfun::apmc(X = L,
test = "single-step") %>%
mutate_at(c("fit", "lwr", "upr"),
m0$family$linkinv)
+```
-# Gráfico com estimativas, IC e texto.
-ggplot(data = results,
- mapping = aes(x = reorder(trat, fit), y = fit)) +
- geom_point() +
- geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
- width = 0) +
- geom_label(mapping = aes(label = sprintf("%0.2f %s", fit, cld)),
- nudge_y = 0.25,
- vjust = 0) +
- labs(x = "Tratamento", y = "Número de colmos atacados")
+A médias ajustadas na escala da resposta, por ter uma função de ligação
+envolvida e uma agregação para marginalizar o efeito, podem ser
+calculadas de outra forma, trocando a ordem das operações marginalizar e
+levar para a escala da resposta.
+
+ * A primeira opção consiste em 1) obter os valores preditos na escala
+ do preditor para cada um dos $6^3$ pontos experimentais, 2) fazer a
+ agregação por média para os níveis de tratamento e 3) passar as
+ estimativas para a escala da resposta aplicando $g^{-1}(.)$. Fez-se
+ exatamente isso no código acima usando a
+ `emmeans(..., type = "response")`.
+ * A segunda opção é 1) obter os valores preditos na escala do preditor
+ para cada um dos $6^3$ pontos experimentais, 2) passar as
+ estimativas para a escala da resposta aplicando $g^{-1}(.)$ e 3)
+ fazer a agregação por média para os níveis de tratamento.
+
+A diferença é inversão na ordem das etapas 2 e 3. Como `mean(exp(x))` é
+diferente de `exp(mean(x))`, os resultados serão diferentes. O código
+abaixo reproduz os resultados.
+
+```{r}
+# Duas chamadas para obter as médias ajustadas.
+emmeans(m0, specs = ~trat, type = "response")
+emmeans(m0, specs = ~trat, transform = "response")
+```
+
+As duas chamadas acima calculam as médias ajustadas das duas formas
+discutidas. O código a seguir reproduz os resultados para as estimativas
+pontuais com o qual fica claro a troca na ordem das operações de
+transformação e agregação.
+
+O cálculo dos erros padrões para `transform = "response"` aplicam o
+método delta, como mostra o código abaixo. Por outro lado, para os
+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}
+# 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 %>%
+ select(linha, coluna, trat) %>%
+ complete(linha, coluna, trat)
+grid$eta <- predict(m0, newdata = grid)
+# grid$mu <- predict(m0, newdata = grid, type = "response")
+
+# Código reproduz: emmeans(m0, specs = ~trat, type = "response").
+grid %>%
+ group_by(trat) %>%
+ summarise(eta = mean(eta)) %>%
+ mutate(mu = m0$family$linkinv(eta), eta = NULL) %>%
+ ungroup()
+
+# Código reproduz: emmeans(m0, specs = ~trat, transform = "response").
+grid %>%
+ mutate(mu = m0$family$linkinv(eta)) %>%
+ group_by(trat) %>%
+ summarise(mu = mean(mu)) %>%
+ ungroup()
+
+# Média ajustada para a condição trat = 1 e linha = 1 (média nas colunas).
+emmeans(m0,
+ specs = ~trat,
+ transform = "response",
+ at = list(linha = "1", trat = "1")) %>%
+ as.data.frame()
+
+# A expressão aplicada pela chamada da `emmeans()` acima.
+g <- sprintf("1/6 * exp(%s)",
+ c("(Intercept)",
+ paste("(Intercept)",
+ paste0("coluna", 2:6),
+ sep = " + "))) %>%
+ paste(collapse = " + ")
+cat(strwrap(g), sep = "\n")
+
+# A aplicação do método delta.
+dm <- car::deltaMethod(m0, g = g)
+rownames(dm) <- NULL
+dm
+```
+
+Conclusões: TODO!
+
+```{r, eval = FALSE}
+# Análise com modelo linear Gaussiano e relação \sigma^2 = \mu.
+m0 <- glm(colmos ~ linha + coluna + trat,
+ data = da,
+ family = quasi(variance = "mu"))
+
+# Diganóstico dos resíduos.
+par(mfrow = c(2, 2))
+plot(m0)
+layout(1)
```
## Variável resposta de proporção amostral