Skip to content
Snippets Groups Projects
Commit ff4dffd4 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Improves second section.

parent 9e77c0a0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 181 additions and 64 deletions
04-dql.Rmd
+ 181
− 64
View file @ ff4dffd4
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment