diff --git a/11-fatorial-duplo-2-adicionais.Rmd b/11-fatorial-duplo-2-adicionais.Rmd
index a2bd06be4faff1e83bc0ffd520a22cf93d92be0b..d74164770b2f64beced823a623b3638c9d11078b 100644
--- a/11-fatorial-duplo-2-adicionais.Rmd
+++ b/11-fatorial-duplo-2-adicionais.Rmd
@@ -307,6 +307,49 @@ da %>% group_by(!!x) %>%
## Análise de dados de *fusarium* em Zimmermann
+Os dados que usaremos para aplicar a análise acima prototipada estão em
+@zimmermann2014. Trata-se de um experimento fatorial $3 \times 3 + 2$,
+ou seja, um fatorial duplo com mais dois pontos adicionais. O
+experimento foi realizado no delineamento inteiramente casualizado.
+
+A variável resposta é o número de sementes infectadas com *fusarium* de
+um total de 40 sementes expostas (`nsi`). Portanto, a resposta pode ser
+analisada considerando distribuição binomial se forem razoáveis as
+suposições para tal. Digo razoáveis porque nem todo experimento que
+produza uma resposta assim -- um número inteiro no numerador tendo outro
+número inteiro no denominador -- terá distribuição binomial. Existem
+suposições, como independência entre os ensaios de Bernoulli e
+probabilidade de sucesso homogênea em cada ponto
+experimental. Geralmente as sementes apresentam tamanho variável, estão
+em condições fisiológicas diferentes e são colocadas para germinar todas
+juntas em um ambiente compartilhado (gerbox, faixa de canteiro, bandeja
+de germinação). A probabilidade de sucesso pode mudar com o
+estado/tamanho/idade da semente dentro de um mesma unidade experimental,
+o que faz com que a variância observada na resposta seja maior que
+àquela esperada pela distribuição binomial. Além do mais, estarem
+germinando em um ambiente compartilhado pode não garantir independência
+das germinações individuais das sementes, ou seja, uma semente infectada
+pode produzir inóculo que aumenta a chance de infecção das demais
+sementes da mesma repetição.
+
+É bastante comum de ver em livros de Estatística Experimental o uso da
+transformação arco seno da raíz quadrada para dados de proporção como
+estes. Essa transformação é uma transformação estabilizadora da
+variância para dados de distribuição binomial[^1]. Lembre-se que uma
+variável aleatória com distribuição binomial de parâmetros $n$ e $p$ tem
+variância $p(1 - p)$. A transformação $f(x) = \arcsin(\sqrt(x))$ é
+empregada para levar os dados para uma escala com maior homogeneidade de
+variâncias e assim realizar análise assumindo-se distribuição
+gaussiana. Hoje, tem-se os modelos lineares generalizados (GLM)
+disponíveis para a análise de dados de proporção como este, o que
+dimunui a justificativa para uso de transformações estabilizadoras da
+variância. De qualquer forma, com propósito didático, eu farei a análise
+com a variável transformada e depois usando GLM.
+
+Os dados que vamos usar estão no objeto `ZimmermannTb9.22` do pacote
+`labestData`. Consulte a documentação para detalhes. A variável `y` é a
+versão transformada da proporção amostral de sementes infectadas.
+
```{r}
library(labestData)
@@ -323,41 +366,395 @@ str(tb)
# Contagem dos pontos experimentais.
xtabs(~fung + aplic, data = tb)
+# Proporção amostral de sementes infectadas em cada ponto experimental.
+with(tb, tapply(nsi/40, list(fung, aplic), FUN = mean))
+
# Análise exploratória.
ggplot(data = tb,
mapping = aes(x = aplic, y = y)) +
facet_grid(facets = . ~ fung, scale = "free", space = "free") +
- geom_jitter(alpha = 0.5, width = 0.1)
+ geom_jitter(alpha = 0.5, width = 0.1, height = 0)
```
-<!---------------------------------------------------------------------- -->
+Neste experimento o fator fungicida (`fung`) corresponde ao produto
+usado para controle de *fusarium* com 3 níveis qualitativos: Benlate,
+Captam e Derosal. O fator aplicação (`aplic`) indica a forma de
+aplicação de fungicida quando ao uso associado de polímero, também com 3
+níveis: fungicida puro, fungicida antes do polímero e fungicida
+misturado ao polímero. Os pontos adicionais são a testemunha sem
+controle de *fusarium* e a aplicação isolada de polímero.
+
+Para a análise do fatorial com tratamento adicional temos que escolher o
+conjunto de hipóteses mais adequado. Em @zimmermann2014, é feita a
+análise considerando as hipóteses da cláusula 3 que, neste caso, será:
+i) testemunha vs polímero, ii) (testemunha + polímero) vs pontos
+fatoriais e iii) entre pontos fatoriais. A hipótese (i) deseja saber se
+a aplicação de polímero altera a ocorrência de *fusarium*. A (iii)
+deseja avaliar a existência de efeito do fungicída, da forma de
+aplicação, ou da interação entre eles. A (ii) é, no caso de (i) e (iii)
+não serem rejeitadas, se o controle conseguido com qualquer aplicação
+dos fungicidas difere das testemunhas. É necessário enfatizar que a
+interpretação das hipóteses é sequêncial.
+
+Cada ponto experimental foi teve 5 repetições. A proporção média de
+sementes infectadas apresenta valor 0 para um ponto experimental
+(`derosal:misturado`). Essa característica é relevante para a análise
+com GLM binomial, portanto, discutirei adiante.
+
+O diagrama de dispersão sugere existência de difença entre os
+tratamentos para o controle de *fusarium*. Outro aspecto igualmente
+saliente é que existem muitos zeros observados. No ponto experimental
+`fung = derosal & aplic = misturado` foi observado zero em todas as
+unidades experimentais. Dois alertas decorrem disso: estimativa amostral
+de proporção de infecção é $\hat{p} = 0$, ou seja, foram do espaço
+paramétrico de $p \in (0, 1)$ da distribuição binomial e variância
+amostral igual 0. De qualquer forma, essa primeira análise assume
+distribuição para a variável transformada.
+
+As hipóteses que serão testadas já foram estabelecidas. Agora precisa-se
+criar o conjunto de contrastes que as represente. Isso é feito no bloco
+a seguir, seguido do ajuste do modelo. A 5º coluna da matriz de
+contrastes é para a hipótese (i) e a 4º coluna é para a hipótese
+(ii). Os coeficientes tem que somar sempre zero porque são contrastes. O
+conjunto de contrastes é ortogonal, ou seja, o produto interno das
+colunas é 0.
+```{r}
+# Contrates para o fator fungicida.
+ctr <- C(tb$fung, contr = contr.helmert)
+attr(ctr, "contrasts")[, 3] <- c(-2, -2, -2, 3, 3)
+attr(ctr, "contrasts")[, 4] <- c( 0, 0, 0, -1, 1)
+contrasts(tb$fung) <- attr(ctr, "contrasts")
-```{r, eval = FALSE}
-tb <- expand.grid(NR = c("Baixo", "Convencional"),
- Fitato = c("Alto", "Médio"),
- Fitase = c(0, 750, 1500))
-tb <- tb[with(tb, !(NR == "Convencional" & Fitase > 0)), ]
-ftable(xtabs(~NR + Fitato + Fitase, data = tb))
+# Contrates para o fator aplicação.
+ctr <- C(tb$aplic, contr = contr.helmert)
+attr(ctr, "contrasts")[, 3] <- c(-2, -2, -2, 3, 3)
+attr(ctr, "contrasts")[, 4] <- c( 0, 0, 0, -1, 1)
+contrasts(tb$aplic) <- attr(ctr, "contrasts")
-# Importação dados dados.
-tb <- read_tsv("./data-raw/fat2x3+2adic.txt",
- comment = "#",
- col_types = cols())
-attr(tb, "spec") <- NULL
-str(tb)
+# Mostra que o produto interno é 0.
+crossprod(contrasts(tb$fung))
-# Criando o fator NR.
-tb$nr <- "Baixo"
-tb$nr[with(tb, tratamento %in% c(1, 5))] <- "Convencional"
+# Criação da matriz de delineamento.
+X_full <- model.matrix(~fung * aplic, data = tb)
+attr(X_full, "contrasts")
-# Contagem das ocorrências dos pontos experimentais.
-ftable(xtabs(~nr + fitato + fitase, data = tb))
+# Ajuste do modelo.
+m0 <- aov(y ~ 0 + X_full, data = tb)
+cbind(coef(m0))
-# Análise exploratória.
-ggplot(data = tb,
- mapping = aes(x = fitase, y = cdaeb, color = nr)) +
- facet_wrap(facets = ~fitato) +
+# Diagnóstico não seja discutido aqui para não sair do foco.
+# par(mfrow = c(2, 2))
+# plot(m0)
+# layout(1)
+
+# Para decompor as SQ.
+ef <- list("ii) T+P vs Fat" = 4,
+ "i) T vs P" = 5,
+ "Fung" = 2:3,
+ "Aplic" = 6:7,
+ "Fung:Aplic" = 8:11)
+summary(m0, split = list(X_full = ef))
+```
+
+O quadro de anova acima está igual ao disponível em @zimmermann2014
+página 188, exceto por arredondamentos e ordem das linhas.
+
+Pela inspeção do quadro de análise de variância, rejeita-se a hipótese
+nula de não interação entre fungicida e aplicação ao nível de 5%. Dessa
+forma, não interpreta-se os efeitos principais. Pelo quadro, não se
+rejeita a hipótese (i), ou seja, em termos práticos, não houve evidência
+de que polímero isolado tem efeito sobre *fusarium*. A hipótese (ii) não
+faz sentido ser intepretada haja visto a existência de interação.
+
+Agora que usamos a `aov()`, vamos ajustar o mesmo modelo com a `lm()`
+para fazermos o estudo da interação com `multcomp::glht()`. Usaremos uma
+matriz de delineamento apenas com colunas de parâmetros estimáveis.
+
+```{r}
+X_ok <- X_full[, str_replace(names(coef(m0)), "X_full", "")]
+m1 <- lm(y ~ 0 + X_ok, data = tb)
+cbind(lm = coef(m1), aov = coef(m0))
+```
+
+Por razões didáticas, será feito o estudo da interação nos dois
+sentidos: comparar aplicações fixando o nível de fungicida e comparar
+fungicidas fixando a forma de aplicação. No entanto, apesar de sempre
+ser possível estudar a interação em todos os sentidos, na prática nem
+todos são igualmente relevantes ou de significado. Tente responder o que
+faz mais sentido da prática para o controle de **fusarium**:
+
+ 1. É primeiro feita a escolha do fungicida e depois a forma de
+ aplicação?
+ 2. É primeiro feita a escolha da forma de aplicação e depois o
+ fungicida?
+ 3. É possível escolher simultâneamente o fungicida e forma de
+ aplicação.
+
+Estou inclinado a pensar que a condição 1 é mais realista que
+a 2. Todavia, o caso 3 parece ser também apropriado. Para o caso 3
+podemos comparar todos os 9 pontos da porção fatorial sem precisar fazer
+o estudo de um fator fixando-se os demais. Começarei por 1 e 2.
+
+```{r}
+# Pontos experimentais únicos de A * B.
+i <- tb %>% select(fung, aplic) %>% duplicated()
+X_un <- X_ok[!i, ] # Estimativas.
+pts <- tb[!i, c("fung", "aplic")] # Pontos experimentais.
+cbind(pts, unname(X_un))
+
+# Estudo dos níveis de aplicação dentro de fungicida.
+L_AinF <- head(by(X_un, pts$fung, FUN = as.matrix), n = -2)
+L_AinF <- lapply(L_AinF,
+ FUN = "rownames<-",
+ head(levels(tb$aplic), n = -2))
+lapply(L_AinF,
+ FUN = function(s) {
+ K <- all_pairwise(s, collapse = " vs ")
+ summary(glht(m1, linfct = K))
+ })
+
+# Estudo dos níveis de fungicida dentro de aplicação.
+L_FinA <- head(by(X_un, pts$aplic, FUN = as.matrix), n = -2)
+L_FinA <- lapply(L_FinA,
+ FUN = "rownames<-",
+ head(levels(tb$fung), n = -2))
+lapply(L_FinA,
+ FUN = function(s) {
+ K <- all_pairwise(s, collapse = " vs ")
+ summary(glht(m1, linfct = K))
+ })
+```
+
+Baseado nas comparações múltiplas feitas, em termos práticos, não
+existem diferenças entre os fungicidas quando são aplicados antes do
+polímero (será que o polímero aplicado depois dimunui a concentração do
+fungicida?). Para o fugicida Derosal e Benlate, não existe diferença
+entre as formas de aplicação. Para o Captam, existe diferença das outras
+formas de aplicação comparadas com a aplicação do fungicida antes do
+polímero.
+
+Com o propósito didático será feita a comparação entre os 9 pontos
+experimentais. Na prática, essa abordagem pode indicar qual é tratamento
+mais promissor para o controle de *fusarium*.
+
+```{r}
+# Comparação entre todos os pontos experimentais.
+i <- !with(pts, fung %in% tail(levels(fung), n = 2))
+L <- X_un[i, ]
+rownames(L) <- with(pts[i, ], sprintf("%s:%s", fung, aplic))
+K <- all_pairwise(L, collapse = " vs ")
+
+# Número de constrastes avaliados.
+c(choose(nrow(L), 2), nrow(K))
+
+summary(glht(m1, linfct = K), test = adjusted(type = "fdr"))
+```
+
+As médias para cada ponto experimental serão determinadas agora,
+acompanhadas do respectivo intervalo de confiança para que os resultados
+sejam apresentados de forma gráfica.
+
+```{r}
+# Obtém os intervalos de confiaça para as médias ajustadas.
+ci <- confint(glht(m1, linfct = X_un),
+ calpha = univariate_calpha())
+ci <- ci$confint
+colnames(ci) <- tolower(colnames(ci))
+
+# Estimativas pontuais com intervalo de confiança de cobertura
+# individual de 95%.
+results <- cbind(pts, ci)
+results
+
+# Exibe as estimativas com IC em gráfico.
+ggplot(data = results,
+ mapping = aes(x = aplic, y = estimate)) +
+ facet_grid(facets = . ~ fung, scale = "free", space = "free") +
+ geom_point() +
+ geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
+ width = 0)
+```
+
+A função [`apmc`] do pacote [`wzRfun`] é uma função de conveniência que
+faz chamada à `glht()`, retorna o intervalo de confiança e as letras
+para indicar a significância dos contrastes. Para usar em modelos em
+que a matriz de desenho é fornecida ao contrário do convencional uso de
+fórmula, precisa adicionar `cld2 = TRUE`.
+
+```{r}
+library(wzRfun)
+
+# Para comparação entre TODOS os pontos experimentais.
+L <- X_un
+rownames(L) <- with(pts, sprintf("%s:%s", fung, aplic))
+
+# Estimativas pontuais seguido de IC e letras para os contrastes.
+results <- apmc(X = L,
+ model = m1,
+ focus = "trat",
+ test = "fdr",
+ cld2 = TRUE) %>%
+ separate(col = "trat", into = c("fung", "aplic")) %>%
+ mutate(cld = ordered_cld(cld, -fit)) %>%
+ mutate(txt = sprintf("%0.2f%s", fit, cld))
+results
+
+# Gráfico com estimativas, IC e texto.
+ggplot(data = results,
+ mapping = aes(x = aplic, y = fit)) +
+ facet_grid(facets = . ~ fung, scale = "free", space = "free") +
geom_point() +
- stat_summary(geom = "line", fun.y = "mean")
+ geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
+ width = 0) +
+ geom_text(mapping = aes(label = txt),
+ angle = 90,
+ vjust = 0,
+ nudge_x = -0.1)
```
+
+Agora que a análise chegou ao final, pode-se fazer alguns comentários
+importantes que justificam a adoção do modelo binomial. A mais saliente
+é que o limite inferior do intervalo de confiança de algumas estimativas
+de média é negativo. Isso é algo que não tem significado prático.
+
+A função $f(x) = \arcsin(\sqrt{x})$ tem limites
+$$
+\begin{aligned}
+ \lim_{x \to 0} \arcsin(\sqrt{x}) &= 0 \\
+ \lim_{x \to 1} \arcsin(\sqrt{x}) &= \pi/2
+\end{aligned}
+$$
+e a $f^{-1}(y)$ é $\sin^{2}(y)$, $y \in [0, \pi/2]$.
+
+```{r}
+curve(asin(sqrt(x)), from = 0, to = 1, asp = 1, xlim = c(0, pi/2))
+curve(sin(x)^2, from = 0, to = pi/2, add = TRUE, col = 2)
+abline(a = 0, b = 1, lty = 2)
+legend("bottomright", col = c(1, 2), lty = 1, bty = "n",
+ legend = c(expression(asin(sqrt(x))),
+ expression(sin^2 * (x))))
+```
+
+Podemos aplicar a tranformação inversa para interpretar as estimativas
+na escala original de proporção de sementes infectadas. Mas como os
+valores limites inferiores são negativos, ao fazer-se $\sin^2(x)$, eles
+se tornarão positivos e podem, inclusive deixar a estimativa pontual
+fora do intervalo de confiança, o que seria bastante estranho.
+
+```{r}
+# Apenas para provar o ponto mencionado.
+results %>%
+ mutate_at(c("fit", "lwr", "upr"), function(x) sin(x)^2) %>%
+ ggplot(mapping = aes(x = aplic, y = fit)) +
+ facet_grid(facets = . ~ fung, scale = "free", space = "free") +
+ geom_point() +
+ geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
+ width = 0)
+```
+
+A seguir faremos a análise baseada na distribuição binomial.
+
+## Análise dos dados de *fusarium* com binomial
+
+Por se tratar de um manual de análise, não farei apresentação de todos
+dos aspectos teóricos relacionados à classe de modelos lineares
+generalizados. Recomendo que procure a literatura especializada para
+isso. Surigo começar pelo material do Prof. Gauss & Profa. Clarice[^2] e
+o materiial do Prof. Gilberto de Paula[^3].
+
+O que tenho que comentar já foi brevemente mencionado: estimativa de
+proporção amostral igual a 0. Como os fatores experimentais são todos
+qualitativos, o modelo para porporção é baseado em estatísticas
+amostrais dos pontos experimentais. No entanto, sabemos da distribuição
+binomial que o espaço paramétrico de $p$ é o intervalo aberto $(0,
+1)$. Ou seja, teremos problema ao ajustar o modelo aos dados por causa
+de um ponto experimental com estimativa 0. Além do mais, o erro padrão
+da estimativa assume um valor muito grande.
+
+TODO: o problema está no erro padrão da estimativa.
+
+```{r}
+# Dados simulados para mostrar o problema da proporção amostral 0 em um
+# ponto experimental.
+y <- rep(0, 5)
+m0 <- glm(cbind(y, 40 - y) ~ 1, family = binomial)
+summary(m0)
+confint.default(m0)
+round(m0$family$linkinv(confint.default(m0)), 1)
+
+# Proporção amostral de sementes infectadas em cada ponto experimental.
+with(tb, tapply(nsi/40, list(fung, aplic), FUN = mean))
+
+tb$nsim <- tb$nsi + 0.2
+m0 <- glm(cbind(nsim, 40 - nsim) ~ fung * aplic,
+ data = tb,
+ family = quasibinomial)
+anova(m0, test = "F")
+
+m1 <- update(m0, . ~ 0 + X_ok)
+summary(m1)
+
+cbind(m0 = coef(m0)[sub("X_ok", "", names(coef(m1)))],
+ m1 = coef(m1))
+
+# Diagnóstico não será discutivo para não afastar do foco.
+# par(mfrow = c(2, 2))
+# plot(m0)
+# layout(1)
+
+# Teste de Wald para hipóteses globais.
+names(ef)
+cftest(model = m1, parm = ef[[1]], test = Ftest())
+cftest(model = m1, parm = ef[[2]], test = Ftest())
+cftest(model = m1, parm = ef[[3]], test = Ftest())
+cftest(model = m1, parm = ef[[4]], test = Ftest())
+cftest(model = m1, parm = ef[[5]], test = Ftest())
+
+#
+```
+
+Existem alguns importantes aspectos para se comentar aqui.
+
+ * Sequencial vs Marginal.
+ * Diferença de deviance vs Wald.
+
+```{r}
+# Para comparação entre TODOS os pontos experimentais.
+L <- X_un
+rownames(L) <- with(pts, sprintf("%s:%s", fung, aplic))
+
+# Estimativas pontuais seguido de IC e letras para os contrastes.
+results <- apmc(X = L,
+ model = m1,
+ focus = "trat",
+ test = "fdr",
+ cld2 = TRUE) %>%
+ separate(col = "trat", into = c("fung", "aplic")) %>%
+ mutate_at(c("fit", "lwr", "upr"), m1$family$linkinv) %>%
+ mutate(cld = ordered_cld(cld, -fit)) %>%
+ mutate(txt = sprintf("%0.2f%s", fit, cld))
+results
+
+# Gráfico com estimativas, IC e texto.
+ggplot(data = results,
+ mapping = aes(x = aplic, y = fit)) +
+ facet_grid(facets = . ~ fung, scale = "free", space = "free") +
+ geom_point() +
+ geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
+ width = 0) +
+ geom_text(mapping = aes(label = txt),
+ angle = 90,
+ vjust = 0,
+ nudge_x = -0.1)
+
+#
+```
+
+<!------------------------------------------- -->
+[^1]: https://www.google.com/search?q=stabilizing+variance+functions
+[^2]: http://www.esalq.usp.br/departamentos/lce/arquivos/aulas/2013/LCE5868/livro.pdf
+[^3]: https://www.ime.usp.br/~giapaula/texto_2013.pdf
+
+[`apmc`]: https://github.com/walmes/wzRfun/blob/master/R/pairwise.R
+[`wzRfun`]: https://github.com/walmes/wzRfun