diff --git a/04-dql.Rmd b/04-dql.Rmd
index afc6f6407fa604bad9d445762fed06c203616a49..06f824000cf31d25b20d3faa2ca1fe8b6d56fddf 100644
--- a/04-dql.Rmd
+++ b/04-dql.Rmd
@@ -364,15 +364,15 @@ pela média. Assim, a média marginal para uma variedade $k$ é obtida por
$$
\begin{align*}
\hat{\mu}_{..k} &=
- \frac{1}{5 \cdot 5} \sum_{i} \sum_{j} \hat{\mu}_{ijk} \\
+ \frac{1}{I \cdot J} \sum_{i} \sum_{j} \hat{\mu}_{ijk} \\
\hat{\mu}_{..k} &= \hat{\mu} +
- \frac{1}{5} \sum_{i} \hat{\alpha}_{i} +
- \frac{1}{5} \sum_{j} \hat{\theta}_{i} +
+ \frac{1}{I} \sum_{i} \hat{\alpha}_{i} +
+ \frac{1}{J} \sum_{j} \hat{\theta}_{j} +
\hat{\tau}_k,
\end{align*}
$$
-em que o 5 nos denominadores deve-se ao fato de ser um quadrado latino
-$5 \times 5$.
+em que o $I = J = 5$ nos denominadores deve-se ao fato de ser um
+quadrado latino $5 \times 5$.
A função `emmeans()` calcula tais médias com erro-padrão e intervalo de
confiança embutidos. O código abaixo extrai a matriz usada para o
@@ -640,8 +640,6 @@ sobre uma relação média-variância não nula. Não é proveitoso julgar a
normalidade diante de uma suspeita sobre a violação de outros
pressupostos, no entanto, não se vê grave afastamento dessa suposição.
-Uma solução possível é transformar
-
Para diminuir a violação da suposição de variância constante pode-se
tentar aplicar uma transformação na variável resposta. Busca-se aqui
então uma transformação que estabilize a relação média-variância. A
@@ -693,7 +691,7 @@ delineamento de quadrado latino é
$$
\begin{aligned}
Y_{ijk} &= \text{Quasi-Poisson}(\lambda_{ijk}, \phi) \\
- g(\lambda_{ijk}) &= \mu + \alpha_i + \theta_j + \tau_k \\
+ g(\lambda_{ijk}) &= \eta_{ijk} = \mu + \alpha_i + \theta_j + \tau_k \\
\phi &\propto 1 \\
\text{Var}(Y_{ijk}) &= \phi \text{E}(Y_{ijk})
\end{aligned}
@@ -730,13 +728,93 @@ plot(m0)
layout(1)
```
-TODO discutir os resíduos.
+Os gráficos para diagnóstico não apontam violação dos pressupostos para
+a especificação Quasi-Poisson.
+
+O gráfico *Residuals vs Fitted* exibe os resíduos de deviance contra os
+valores ajustados (ou preditos) na escala do preditor linear, ou seja,
+na escala $\eta = \log(\mu)$ e não $\mu$. Não há padrões que apontem a
+violão da suposição de correta especificação do modelo para a
+média. Além disso, os resíduos ocupam as metades superior e inferior de
+forma igualmente distribuída ao longo do eixo das coordenadas.
+
+O gráfico *Scale-Location* apresenta a raíz quadrada do absoluto dos
+resíduos de deviance padronizados contra os valores ajustados na escala
+do preditor linear. Se os pressupostos não forem violados, os resíduos
+de *deviance* padronizados apresentarão assintoticamente distribuição
+normal. Os resíduos de deviance (e também os de Pearson) padronizados
+não devem apresentar relação com a média, porque a relação
+média-variância esperada foi usada para padronizá-los. Conclui-se que
+não há evidências contra a suposição assumida para a média-variância.
+
+O gráfico *Normal Q-Q* usa os resíduos de deviance padronizados com os
+quais não se vê violação da suposição de normalidade.
+
+Por fim, o gráfico *Residuals vs Leverage* mostra os resíduos de Pearson
+padronizados contra os valores de alavancagem. Sem entrar muito em
+detalhes, conclui-se que não existem observações influentes. Um aspecto
+interessante é que, diferente do que foi visto para o modelo Gaussiano,
+as alavancagens não são todas iguais. Isso decorre da função de
+variância que acomoda o peso das observações conforme o valor de cada
+uma.
+
+Com o compromisso de ser didático, o bloco de código abaixo mostra o que
+é usado nos gráficos de resíduos que foram discutidos.
+
+<!-- .
+ https://www.sagepub.com/sites/default/files/upm-binaries/38503_Chapter6.pdf
+ https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/BinomTools/inst/ResidualsGLM.pdf?revision=6&root=binomtools&pathrev=6
+-->
+
+```{r, eval = FALSE}
+# Diferentes tipos de resíduos em GLM.
+R <- cbind(
+ raw = m0$model[, 1] - fitted(m0),
+ res = residuals(m0, type = "response"),
+ dev = residuals(m0, type = "deviance"),
+ dev_std = residuals(m0, type = "deviance")/
+ (sqrt(summary(m0)$dispersion) * sqrt(1 - hatvalues(m0))),
+ pea = residuals(m0, type = "pearson"),
+ pea_std = residuals(m0, type = "pearson")/
+ (sqrt(summary(m0)$dispersion) * sqrt(1 - hatvalues(m0))))
+
+plot(R[, "raw"] ~ R[, "res"])
+
+# Residuals vs Fitted.
+plot(m0, which = 1)
+points(R[, "dev"] ~ predict(m0), pch = 3, col = 2)
+
+# Scale-Location.
+plot(m0, which = 3)
+points(sqrt(abs(R[, "dev_std"])) ~ predict(m0),
+ pch = 3, col = 2)
+
+# Normal Q-Q.
+plot(m0, which = 2)
+qq <- qqnorm(R[, "dev_std"], plot.it = FALSE)
+points(qq$x, qq$y, pch = 3, col = 2)
+
+# Residuals vs Leverage
+plot(m0, which = 5)
+points(R[, "pea_std"] ~ hatvalues(m0),
+ pch = 3, col = 2)
+```
+
+Uma vez que não se constatou afastamento das suposições, dá-se
+continuidade com a avaliação da significância dos termos
+experimentais. O quadro de análise de deviance pode ser produzido com a
+função método `anova()` (com hipóteses sequenciais) e com a função
+método `drop1()` (com hipóteses marginais).
```{r}
# Quadro de análise de deviance.
# drop1(m0, test = "F") # Hipóteses marginais.
anova(m0, test = "F") # Hipóteses sequenciais.
+# O valor da estatística F para `anova()` e `drop1()`.
+(25.095/5)/summary(m0)$dispersion # Usa resíduos de Pearson.
+(25.095/5)/(deviance(m0)/df.residual(m0)) # Usa resíduos de deviance.
+
# Estimativas dos efeitos e medidas de ajuste.
summary(m0)
@@ -744,25 +822,24 @@ summary(m0)
# sum(residuals(m0, type = "deviance")^2)
# deviance(m0)
-# Como estimar o parâmetro de dispersão.
+# Como é estimado o parâmetro de dispersão.
# summary(m0)$dispersion
sum(residuals(m0, type = "pearson")^2)/df.residual(m0)
```
-TODO WALMES FIXME
-
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.
+inúmeros fatores não controlados para inflacionar 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.
+Visto que teve-se indicações de efeito de tratamentos, será feito estudo
+via comparações mútiplas para detectar diferenças entre eles para
+subsidiar decisões sobre o controle da praga. Apesar de ser outra classe
+de modelos, os métodos usados ainda são os mesmos, ou seja, `emmeans()`
+e associadas.
```{r}
# "Médias" ajustadas na escala do preditor linear.
@@ -793,15 +870,48 @@ da %>%
summarise_at(vars(colmos), mean)
```
-O código abaixo faz as comparações múltiplas entre os
-tratamentos. IMPROVE!
+Por padrão, as médias ajustadas são obtidas na escala do preditor linear
+($\eta$). Com uso da opção `type = "response"`, a `emmeans()` apresenta
+o resultado na escala da resposta aplicando a função de ligação inversa
+($g^{-1}(.)$) e usa o método delta para determinar o erro-padrão. As
+médias ajustadas são calculadas por
+$$
+\begin{align*}
+ \hat{\eta}_{..k} &=
+ \frac{1}{I \cdot J} \sum_{i} \sum_{j} \hat{\eta}_{ijk} \\
+ \hat{\eta}_{..k} &= \hat{\mu} +
+ \frac{1}{I} \sum_{i} \hat{\alpha}_{i} +
+ \frac{1}{J} \sum_{j} \hat{\theta}_{j} +
+ \hat{\tau}_k \\
+ \hat{\mu}_{..k} &= g^{-1}(\hat{\eta}_{..k}).
+\end{align*}
+$$
+
+É importante salientar que, apesar do experimento ser balanceado, as
+médias ajustadas tem erros-padrões diferentes. Quanto maior o valor da
+estimativa, maior o erro-padrão. Isso é decorrente da função de
+variância da especificação Quasi-Poisson.
+
+O código abaixo faz as comparações múltiplas par a par entre os
+tratamentos. O teste é feito com as médias na escala do preditor linear,
+ou seja, os contrastes avaliados são do tipo $\eta_{..k} - \eta_{..k'} =
+0$.
+
+A `contrast()` mostra a estimativa e significância de cada um dos
+contrastes. O p-valores foram corrigidos para um erro global de
+significância de 5%. Apesar do experimento ser balanceado, os
+erros-padrões dos contrastes não são iguais como visto no modelo
+Gaussiano. Isso decorre da existência da função de variância. Isso deixa
+claro que abordagens de comparações múltiplas baseadas em diferença
+mínima significativa (DMS) não são imediatamente aplicáveis porque os
+contrastes não tem a mesma variância.
```{r}
# Desvendando o computo das médias ajustadas.
grid <- attr(emm, "grid")
L <- attr(emm, "linfct")
rownames(L) <- grid[[1]]
-MASS::fractions(t(L))
+# MASS::fractions(t(L))
# Retorna a avaliação de cada contraste.
contrast(emm, method = "pairwise")
@@ -811,7 +921,8 @@ results <- cld(emm)
results <- results %>%
as.data.frame() %>%
mutate_at(c(2, 5, 6), m0$family$linkinv) %>%
- mutate(.group = as_letters(.group))
+ mutate(.group = as_letters(.group)) %>%
+ arrange(desc(emmean))
results
# Gráfico com estimativas, IC e texto.
@@ -823,46 +934,39 @@ ggplot(data = results,
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")),
- test = adjusted(type = "single-step"))
-tk_sgl
-
-# Resumo compacto com letras.
-cld(tk_sgl, decreasing = TRUE)
-
-results <- wzRfun::apmc(X = L,
- model = m0,
- focus = "trat",
- test = "single-step") %>%
- mutate_at(c("fit", "lwr", "upr"),
- m0$family$linkinv)
+ 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.
+envolvida e uma agregação para nivelar o efeito de termos acessórios,
+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
+ do preditor linear, 2) fazer a agregação por média 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")`.
+ exatamente isso no código acima usando a `emmeans(..., type =
+ "response")`. Matematicamente isso é
+ $$
+ \hat{\mu}_{..k} =
+ g^{-1}\left(\frac{1}{I \cdot J} \sum_{i} \sum_{j}
+ \hat{\eta}_{ijk}\right).
+ $$
* 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.
+ linear 2) passar as estimativas para a escala da resposta aplicando
+ $g^{-1}(.)$ e 3) fazer a agregação por média. Matematicamente isso
+ é
+ $$
+ \hat{\mu}_{..k} =
+ \frac{1}{I \cdot J} \sum_{i} \sum_{j}
+ g^{-1}\left(\hat{\eta}_{ijk}\right).
+ $$
+
+A diferença é inversão na ordem das etapas 2 e 3. Como $\sum_t
+g^{-1}(u_t)$ é diferente de $g^{-1}(\sum_t u_t)$, os resultados serão
+diferentes. O código abaixo mostra como obter as médias ajustadas nas
+duas opções.
```{r}
# Duas chamadas para obter as médias ajustadas.
@@ -870,20 +974,17 @@ 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
+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"`
+aplica o método delta, como será mostrado. 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}
# 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}.
+# escala do preditor: \eta.
grid <- da %>%
select(linha, coluna, trat) %>%
complete(linha, coluna, trat)
@@ -904,7 +1005,7 @@ grid %>%
summarise(mu = mean(mu)) %>%
ungroup()
-# Média ajustada para a condição trat = 1 e linha = 1 (média nas colunas).
+# Média ajustada para a condição trat = 1 e linha = 1.
emmeans(m0,
specs = ~trat,
transform = "response",
@@ -926,7 +1027,17 @@ rownames(dm) <- NULL
dm
```
-Conclusões: TODO!
+Com o que foi mostrado e discutido acima, encerra-se a análise dos dados
+de contagem. Porém, outras possibilidades de análise existem, como usar
+outra distribuição (Binomial Negativa, Gamma-Count, COM-Poisson, etc) ou
+outras especificações na `glm()`. No fragmento de código abaixo está
+indicado duas especificações equivalentes em que a função de ligação é
+indentidade e a função de variância é linear. Isso corresponde a
+
+ * Um modelo Gaussiano com função de variância do tipo $\text{Var}(Y) =
+ \phi \text{E}(Y)$, ou
+ * Um modelo Quasi-Poisson com função de ligação identidade, $g(\mu) =
+ \mu = \eta$.
```{r, eval = FALSE}
# Análise com modelo linear Gaussiano e relação \sigma^2 = \mu.
@@ -951,6 +1062,12 @@ c(deviance(m1), deviance(m2))
cbind(coef(m1), coef(m2))
```
+Como todos os fatores do experimento são qualitativos, usar a função de
+ligação identidade para um modelo Quasi-Poisson não apresenta grandes
+riscos de obter predições fora do espaço paramétrico. Sendo assim, a
+análise se torna um pouco mais fácil por não precisar se preocupar com
+escalas.
+
## Variável resposta de proporção amostral
Assim como dados de contagem, dados de proporção são muito comuns em