Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
mpaer
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Walmes Marques Zeviani
mpaer
Commits
12554eae
Commit
12554eae
authored
5 years ago
by
Walmes Marques Zeviani
Browse files
Options
Downloads
Patches
Plain Diff
Improves chapter.
parent
8b653fbe
No related branches found
No related tags found
No related merge requests found
Changes
1
Show whitespace changes
Inline
Side-by-side
Showing
1 changed file
04-dql.Rmd
+232
-52
232 additions, 52 deletions
04-dql.Rmd
with
232 additions
and
52 deletions
04-dql.Rmd
+
232
−
52
View file @
12554eae
...
...
@@ -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
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment