Adds fatorial with additional experimental points.

10-fatorial-duplo-1-adicional.Rmd

 # Fatorial duplo com um tratamento adicional
 
-```{r, eval = FALSE}
-library(labestData)
+Para começarmos esse tutorial serão necessários os pacote `car` e
+`multcomp`. O primeiro contém a função `glht()` para realização de
+testes de hipóteses lineares gerais, úteis para aplicar procedimentos de
+comparação múltipla de hipóteses e contrastes planejados. Do segundo
+será usada a função `linearHypothesis()` que aplica o teste de Wald.  O
+pacote `tidyverse` utilizaremos para operações tabulares.
 
-ls("package:labestData")
+```{r, message = FALSE}
+# ATTENTION: Limpa espaço de trabalho.
+rm(list = objects())
 
-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))
+# Carrega pacotes.
+library(car)       # linearHypothesis()
+library(multcomp)  # glht()
+library(tidyverse)
 ```
 
-```{r, message = FALSE}
-library(tidyverse)
+Faremos uso de um conjunto de dados artificial para demonstrar o
+procedimento dessa análise. Uma parte importantíssima é fazer com que a
+testemunha seja último nível dos fatores. Dessa forma, com o uso de
+contrastes de Helmert teremos a hipótese da média do tratamento
+adicional igual a média dos pontos fatoriais especificada.  Como
+variável resposta colocaremos apenas um ruído branco gaussiano.
+Conforme a tabela de ocorrência cruzada dos pontos experimentais, temos
+com esses dados simulados, um experimento com estrutura fatorial de 2
+$\times$ 4 mais um tratamento adicional.
+
+```{r}
+# Criando um conjunto de dados artificial.
+da <- rbind(
+    # Parte fatorial.
+    crossing(A = 1:2,
+             B = 1:4,
+             r = 1:3),
+    # Parte adicional.
+    tibble(A = 0, B = 0, r = 1:3))
+
+# Converte fatores experimentais qualitativos para classe
+# correspondente. IMPORTANT: para o procedimento que será adotado, é
+# preciso que o tratamento adicional seja o último nível.
+da <- mutate(da,
+             A = fct_relevel(factor(A), "0", after = Inf),
+             B = fct_relevel(factor(B), "0", after = Inf))
+
+# Adiciona uma resposta que é apenas ruído branco.
+set.seed(18900217) # Nascimento de Fisher.
+da$y <- rnorm(nrow(da))
+str(da)
+
+# Estrutura de um fatorial 2 x 4 + 1 com 3 repetições.
+xtabs(~A + B, data = da)
+```
+
+Quando ajustarmos um modelo utilizando a função `lm()` temos efeitos não
+estimados uma vez que não existem todos os pontos experimentais do
+experimento fatorial completo.  A `model.matrix()` cria a matriz de
+delineamento com todos os termos de interação (abordagem maximal) e eles
+só serão estimados se houverem observações para todos os pontos
+experimentais do fatorial completo.
+
+É exatamente este aspecto da análise que dificulta: ser um fatorial
+incompleto.  O que será mostrado aqui são formas de contornar esta
+complicação para que a análise possa ser feita da forma mais
+reproduzível possível.
+
+```{r}
+# ATTENTION: alguns parâmetros não serão estimados porque são baseados
+# em pontos experimentais que não existem.
+m0 <- lm(y ~ A * B, data = da)
+cbind(coef(m0))
+```
+
+Os parâmetros com estimativas `NA` não foram estimados. Em geral existem
+duas principais razões para isso:
+
+  1. Não existem pontos experimentais observados para estimar tais
+     parâmetros. Esse é o caso nessa análise.
+  2. Existem variáveis completamente redundantes sendo declaradas,
+     portanto, presença de colinearidade perfeita, o que faz que efeitos
+     não sejam estimados. Esse é um caso que ocorre em regressão.
+
+O uso dos contrastes de Helmert é um facilitador haja visto que a última
+hipótese corresponde a diferença de média do ponto adicional contra
+média dos pontos fatoriais (hipótese FvsA).  Apesar da comodidade,
+outros contrastes poderiam ser empregados, incluindo os definidos pelo
+usuário.
+
+É muito comum que a hipótese FvsA seja testada. No entanto existem
+situações práticas em que ela pode não fazer sentido. Do ponto de vista
+estatístico, por exemplo, ela só terá significado quando efeito dos
+fatores `A` e `B` forem nulos, caso contrário a média dos pontos
+fatoriais é uma quantidade sem significado.
+
+A seguir cria-se a matriz do delineamento usando os contrastes de
+Helmert e ajusta-se o modelo com a função `aov()`.  Internamente, a
+`aov()` faz uma chamada da `lm()`. O que difere são métodos adicionais
+para a classe `aov` e a possibilidade de declarar termos de erro na
+fórmula (`Error()`), necessários para análise de experimento em parcelas
+subdivididas, por exemplo.
+
+```{r}
+# Codificação usada no contraste de Helmert (i.e. back difference).  A
+# última hipótese é se há diferença entre tratamento adicional contra a
+# média dos níveis anteriores.
+C(gl(4, 1), contr = contr.helmert)
+
+# Criação da matriz de delineamento.
+X_full <- model.matrix(~A * B,
+                       data = da,
+                       contrasts = list(A = contr.helmert,
+                                        B = contr.helmert))
+# colnames(X_full)
+# unique(X_full)
+attr(X_full, "contrasts")
+
+# Ajuste do modelo com aov() para usar `split =`.
+m0 <- aov(y ~ 0 + X_full, data = da)
+cbind(coef(m0))
+```
+
+A função a `aov()`, diferente da função `lm()` abandona os termos de
+efeito não estimável.  O objeto retornado por ela tem classe `aov` e
+possui métodos adicionais/modificados.  O método `summary()`, que
+retorna o quadro de análise de variância, tem um argumento `split` que
+permite fazer partições das somas de quadrados.  Apenas é necessário
+informar a posição dos coeficientes para indicar como a partição será
+feita.
+
+O segundo coeficiente corresponde ao efeito do fator A. O fator A que
+possui 2 níveis já que o seu terceiro nível é na realidade o ponto
+adicional, também presente como último nível de B.  Os coefientes de 4 a
+6 correspondem ao efeito do fator B que possui quatro níveis, portanto,
+3 efeitos estimados (o quinto nível não estimado é o ponto adicional).
+Os coeficientes de 7 a 9 correspondem ao efeito da interação entre A e B
+(1 $\times$ 3 = 3 graus de liberdade).  Finalmente, o coeficiente na
+posição 3 é o da hipótese FvsA que estabelece que a diferença entre o
+ponto adicional com relação à média dos pontos fatoriais é nula.
+
+Como a variável resposta foi simulada de um ruído branco, não se espera
+observar rejeição para nenhuma das hipóteses. Caso ocorra, é meramente
+circunstâncial.  Depois serão analisados dados reais para que possamos
+exercitar a intepretação.
+
+```{r}
+# Coeficientes e seus efeitos.
+ef <- list("A" = 2,
+           "B" = 4:6,
+           "A:B" = 7:9,
+           "fat_vs_ad" = 3)
+summary(m0, split = list(X_full = ef))
+```
+
+A função `lm()` assume distribuição gaussiana para a variável resposta e
+a `aov()` permite fazer o particionamento das somas de quadrados para as
+hipóteses relacionas aos efeitos experimentais em questão. No entanto,
+para outras classes de modelo (GLM, sobrevivência, modelos de efeitos
+aleatórios), abordagem equivalente a `lm()` e `aov()` não se aplicam.
+Então, o que se pode usar são recursos como `glht()` e
+`linearHypothesis()` que serão demontrados a seguir.  Com essas funções
+pode-se tanto obter uma estatística para a hipótese global de efeito de
+um termo do modelo como aplicar testes para outras hipóteses.
+
+O primeiro passo para usar tais funções é ajustar um modelo de classe
+`lm()` que não tenha efeitos não estimáveis. Então passamos uma matriz
+com colunas de efeito estimáveis eliminando as colunas que sabidamente
+correspondem aos efeitos não estimáveis.
+
+```{r}
+# A `aov()` já abandona parâmetros não estimados do vetor. Mas para
+# modelos de outras classes, terá que ser feito manualmente.
+X_ok <- X_full[, str_replace(names(coef(m0)), "X_full", "")]
+
+# Modelo linear com matriz onde todos efeitos são estimáveis.
+m1 <- lm(y ~ 0 + X_ok, data = da)
+cbind(lm = coef(m1), aov = coef(m0))
+```
+
+A função `multcomp::cftest()` testa hipóteses sobre vetores de
+parâmetros pelo teste de Wald com estatística F ou $\chi$-quadrado.
+Para isso, é só indicar o nome ou posição dos parâmetros sob hipótese.
+Note que a estatística F é a mesma daquela vista no `summary()` do
+modelo `m0`.
+
+Em modelos de outras classes (`glm`, `survreg`, `lmeMod`, etc) pode haver
+diferença. Por exemplo, um teste de Wald (que é uma aproximação
+quadrática da verossimilhança) poderá diferir de uma razão de deviances
+ou razão de verossimilhanças.
+
+```{r}
+# `multcomp::cftest()` é aplicável a uma ampla classe de modelos.
+cftest(m1, parm = ef[[1]], test = Ftest()) # A.
+cftest(m1, parm = ef[[2]], test = Ftest()) # B.
+cftest(m1, parm = ef[[3]], test = Ftest()) # A:B.
+cftest(m1, parm = ef[[4]], test = Ftest()) # fat_vs_ad.
+```
+
+Com a `linearHypothesis()` podemos obter os mesmos resultados.  A
+diferença é que ela precisa de uma matriz que espeficique as
+hipóteses. Por isso, pode-se especificar hipóteses mais complexas que
+não seriam possíveis de avaliar com a `cftest()`.  Um exemplo é quando
+fazemos análise de covariância, caso que será mostrado em outro
+capítulo.
+
+```{r}
+# Matriz de hipóteses lineares (linear hypothesis matrix).
+m <- diag(length(coef(m1)))
+lhm <- ef
+for (i in 1:length(ef)) {
+    lhm[[i]] <- m[ef[[i]], , drop = FALSE]
+}
+lhm
+
+# `car::linearHypothesis()` é aplicável a uma ampla classe de modelos.
+linearHypothesis(m1, hypothesis.matrix = lhm[[1]], test = "F") # A.
+linearHypothesis(m1, hypothesis.matrix = lhm[[2]], test = "F") # B.
+linearHypothesis(m1, hypothesis.matrix = lhm[[3]], test = "F") # A:B.
+linearHypothesis(m1, hypothesis.matrix = lhm[[4]], test = "F") # fat_vs_ad.
 ```
+
+Nas instruções acima testamos as hipóteses globais de feito dos fatores
+experimentais: efeito de A, efeito de B, efeito da interação
+A:B. Consideramos também a hipótese FvsA.  Agora será demonstrado como
+testar as hipóteses de comparações múltiplas entre médias para os níveis
+de cada fator.
+
+Vamos testar as hipóteses sobre as diferenças entre médias para os
+níveis do fator A, depois entre médias para os níveis do fator B. Ambos
+assumem que não houve interação entre os fatores A e B.  Quando há
+interação, recomenda-se fazer o estudo de um fator fixando-se os níveis
+dos demais.  Então serão obtidas as comparações múltiplas de médias dos
+dois sentidos: para níveis de B fixando o nível de A e para níveis de A
+fixando o nível de B. Por último, será obtido o teste que envolve a
+média do ponto adicional contra a média dos pontos fatoriais.  Essa
+hipótese já foi testada, porém não foi obtido a estimativa da diferença
+entre tais médias.
+
 ```{r}
-# Importação dados dados.
-tb <- read_tsv("./data-raw/fat2x3+2adic.txt",
-               comment = "#",
-               col_types = cols())
-attr(tb, "spec") <- NULL
-str(tb)
+# Pontos experimentais únicos de A * B.
+i <- da %>% select(A, B) %>% duplicated()
+X_un <- X_ok[!i, ]         # Estimativas.
+pts <- da[!i, c("A", "B")] # Pontos experimentais.
+cbind(pts, X_un)
 
-# Criando o fator NR.
-tb$nr <- "Baixo"
-tb$nr[with(tb, tratamento %in% c(1, 5))] <- "Convencional"
+all_pairwise <- function(L, collapse = "vs", row_names = rownames(L)) {
+    ij <- combn(x = 1:nrow(L), m = 2)
+    K <- apply(X = ij,
+               MARGIN = 2,
+               FUN = function(i) { L[i[1], ] - L[i[2], ] })
+    K <- t(K)
+    if (is.null(row_names)) {
+        nm <- 1:nrow(L)
+    } else {
+        nm <- row_names
+    }
+    rownames(K) <- apply(ij,
+                         MARGIN = 2,
+                         FUN = function(i) {
+                             paste(nm[i], collapse = collapse)
+                         })
+    return(K)
+}
 
-# Contagem das ocorrências dos pontos experimentais.
-ftable(xtabs(~nr + fitato + fitase, data = tb))
+# Contrastes dois a dois (pairwise) entre níveis de A.
+L_A <- do.call(rbind, by(X_un, pts$A, colMeans))
+K_A <- all_pairwise(head(L_A, n = -1))
+summary(glht(m1, linfct = K_A))
 
-# Análise exploratória.
-ggplot(data = tb,
-       mapping = aes(x = fitase, y = cdaeb, color = nr)) +
-    facet_wrap(facets = ~fitato) +
-    geom_point() +
-    stat_summary(geom = "line", fun.y = "mean")
+# Contrastes dois a dois (pairwise) entre níveis de B.
+L_B <- do.call(rbind, by(X_un, pts$B, colMeans))
+K_B <- all_pairwise(head(L_B, n = -1))
+summary(glht(m1, linfct = K_B))
+
+# Estudo da interação: entre níveis de B dentro de A.
+L_BinA <- by(X_un, pts$A, FUN = as.matrix)
+L_BinA <- lapply(L_BinA, "rownames<-", NULL)
+lapply(head(L_BinA, n = -1),
+       FUN = function(s) {
+           summary(glht(m1, linfct = all_pairwise(s)))
+       })
+
+# Estudo da interação: entre níveis de A dentro de B.
+L_AinB <- by(X_un, pts$B, FUN = as.matrix)
+L_AinB <- lapply(L_AinB, "rownames<-", NULL)
+lapply(head(L_AinB, n = -1),
+       FUN = function(s) {
+           summary(glht(m1, linfct = all_pairwise(s)))
+       })
+
+# Contraste do ponto adicional contra a média dos pontos fatoriais.
+L_T <- by(X_un, pts$A == "0", colMeans)
+summary(glht(m1, linfct = rbind(L_T[[1]] - L_T[[2]])))
+```
+
+A título de verificação, vamos calcular as médias amostrais.  Como
+tem-se um experimento balanceado e de efeitos ortogonais, as médias
+amostrais serão iguais as médias ajustadas.  Dessa forma, caso hajam
+diferenças é sinal de que o código precisa ser revisto.  É importante
+enfatizar que em modelos mais gerais (`glm`, `survreg`, `lmerMod`, etc)
+ou para experimentos não balanceados ou de efeitos não ortogonais, as
+médias amostrais não serão iguais as médias ajustadas.
+
+```{r}
+# Diferenças entre médias amostrais apenas para verificação.
+da %>% group_by(A) %>% summarise(m = -mean(y)) %>%
+    pull(m) %>% diff() %>% head(n = -1) %>% setNames(c("1vs2"))
+da %>% group_by(B) %>% summarise(m = -mean(y)) %>%
+    pull(m) %>% diff() %>% head(n = -1) %>%
+    setNames(c("1vs2", "2vs3", "3vs4"))
+da %>% group_by(A != "0") %>% summarise(m = -mean(y)) %>%
+    pull(m) %>% diff() %>% setNames(c("FvsA"))
 ```

11-fatorial-duplo-2-adicionais.Rmd

+# Fatorial duplo com dois tratamentos adicionais
+
+## Prototipação com dados simulados
+
+Usaremos os pacotes `car`, `multcomp` e `tidyverse` pelas mesmas razões
+descritas no capítulo anterior a esse no qual a análise de um
+experimento fatorial com um tratamento adicional foi ilustrada.
+Recomenda-se, para melhor entendimento desse capítulo, que o anterior
+tenha sido lido.
+
+```{r, message = FALSE}
+# ATTENTION: Limpa espaço de trabalho.
+rm(list = objects())
+
+# Carrega pacotes.
+library(car)       # linearHypothesis()
+library(multcomp)  # glht()
+library(tidyverse)
+```
+
+Faremos uso de um conjunto de dados artificial para demonstrar o
+procedimento dessa análise.  Uma parte importantíssima é fazer com que
+os pontos adicionais sejam os últimos níveis dos fatores.  Como variável
+resposta colocaremos apenas um ruído branco gaussiano.  Conforme a
+tabela de ocorrência cruzada dos pontos experimentais, temos com esses
+dados simulados, um experimento com estrutura fatorial de 2 $\times$ 4
+mais 2 tratamentos adicionais.
+
+```{r}
+# Criando um conjunto de dados artificial.
+da <- rbind(
+    # Parte fatorial.
+    crossing(A = as.character(1:2),
+             B = as.character(1:4),
+             r = 1:3),
+    # Parte adicional. P: positivo, N: negativo.
+    tibble(A = "P", B = "P", r = 1:3),
+    tibble(A = "N", B = "N", r = 1:3))
+
+# Converte fatores experimentais qualitativos para classe
+# correspondente. IMPORTANT: para o procedimento que será adotado, é
+# preciso que o tratamento adicional seja o último nível.
+da <- mutate(da,
+             A = fct_relevel(factor(A), c("N", "P"), after = Inf),
+             B = fct_relevel(factor(B), c("N", "P"), after = Inf))
+
+# Adiciona uma resposta que é apenas ruído branco.
+set.seed(18900217) # Nascimento de Fisher.
+da$y <- rnorm(nrow(da))
+str(da)
+
+# Estrutura de um fatorial 2 x 4 + 1 com 3 repetições.
+xtabs(~A + B, data = da)
+```
+
+Quando ajustarmos um modelo utilizando a função `lm()` temos efeitos não
+estimados uma vez que não existem todos os pontos experimentais do
+experimento fatorial completo.  A `model.matrix()` cria a matriz de
+delineamento com todos os termos de interação (abordagem maximal) e eles
+só serão estimados se houverem observações para todos os pontos
+experimentais do fatorial completo.
+
+É exatamente este aspecto da análise que dificulta: ser um fatorial
+incompleto.  O que será mostrado aqui são formas de contornar esta
+complicação para que a análise possa ser feita da forma mais
+reproduzível possível.
+
+```{r}
+# ATTENTION: alguns parâmetros não serão estimados porque são baseados
+# em pontos experimentais que não existem.
+m0 <- lm(y ~ A * B, data = da)
+cbind(coef(m0))
+```
+
+O uso dos contrastes de Helmert foi um facilitador para a análise com um
+ponto adicional.  Com dois pontos adicionais existem outras
+possibilidades das quais a mais adequada depende dos objetivos do
+experimento.  Podemos usar contrastes pra cada uma das seguintes
+hipóteses:
+
+  1. Positivo vs (Fatorial + Negativo), Negativo vs Fatorial, entre
+     pontos fatoriais.
+  2. Negativo vs (Fatorial + Positivo), Positivo vs Fatorial, entre
+     pontos fatoriais.
+  3. Positivo vs Negativo, Fatorial vs (Positivo + Negativo), entre
+     pontos fatoriais.
+  4. Positivo vs Fatorial, Negativo vs Fatorial, entre pontos fatoriais.
+
+Nas 4 opções existe as hipóteses chamadas de "entre pontos fatoriais"
+que são relacionadas ao efeito de A, de B e da interação A:B.  A
+diferença está nas duas sentenças anteriores.  As cláusulas 1 e 2 são
+simétricas, trocando apenas Positivo e Negativo de lugar.  A cláusula 3
+estabelece o contraste entre os pontos adicionais e eles juntos contra
+os pontos fatoriais.  O último caso contrasta cada ponto adicional com
+os pontos fatoriais.
+
+É importante mencionar que destes 4 casos, apenas o último não produz um
+conjunto de contrastes ortogonal.  Uma vez não sendo ortogonal, a
+estatística F do quando de anova (hipóteses sequenciais) poderá ser
+diferente da estatística F do teste de Wald (hipóteses marginais).
+
+A escolha depende das hipóteses do experimento.  Vamos considerar que o
+experimento esteja estudando dois herbicidas (fator A) e combinados com
+4 formas de aplicação no cultivo de soja estudando a produtividade.  O
+ponto adicional positivo ou testemunha positiva é a capina mecânica que
+não provoca intereferencias fisiológicos na planta.  Portanto, espera-se
+que seja o tratamento mais produtivo.  O ponto adicional negativo ou
+testemunha negativa é não fazer a capina, portanto, haverá competição
+das plantas daninhas com a cultura.  Com isso, espera-se que seja o
+tratamento menos produtivo.  Os pontos fatoriais, que são aplicações de
+herbicida de alguma forma, devem estabeler-se, em produtividade, entre
+as duas testemunhas.
+
+Suponha que não existe diferença entre pontos fatoriais, ou seja, todas
+as combinações de herbicida com a forma de aplicação apresentaram a
+mesma produtividade (hipótese nula não rejeitada).  Então o pesquisador
+deseja saber na sequência: i) o uso de herbicida difere de não fazer o
+controle (Negativo)? e ii) a capina manual (Positivo) é tão boa quando o
+uso dos herbicidas?  Para isso, as hipóteses da cláusula 4 são as
+adequadas.
+
+Agora que foi escolhido o conjunto de hipóteses, deve-se usar os
+contrastes correspondentes.  Faremos modificação do constraste de
+Helmert para adequá-lo a cláusula 4 tanto para o fator A quanto para B.
+
+```{r}
+# Atribuindo os constrastes via modificação do Helmert.
+ctr_A <- C(da$A, contr = contr.helmert)
+attr(ctr_A, "contrasts")[, nlevels(da$A) - 1] <- c(-1, -1, 0, 2)
+ctr_B <- C(da$B, contr = contr.helmert)
+attr(ctr_B, "contrasts")[, nlevels(da$B) - 1] <- c(-1, -1, -1, -1, 0, 4)
+
+# Usa contrastes definidos pelo usuário para A e B.
+contrasts(da$A) <- attr(ctr_A, "contrasts")
+contrasts(da$B) <- attr(ctr_B, "contrasts")
+
+# Criação da matriz de delineamento.
+X_full <- model.matrix(~A * B, data = da)
+attr(X_full, "contrasts")
+
+# Ajuste do modelo com aov() para usar `split =`.
+m0 <- aov(y ~ 0 + X_full, data = da)
+cbind(coef(m0))
+```
+
+O modelo `m0` foi ajustado usando-se os contrastes definidos pelo
+usuário motivados pela cláusula 4.  Para particionar as somas de
+quadrados, usamos novamente o método `summary()` com o parâmetro
+`split`.  A lista indicando os termos é passada.
+
+```{r}
+# Coeficientes e seus efeitos.
+ef <- list("fat_vs_N" = 3,
+           "fat_vs_P" = 4,
+           "A" = 2,
+           "B" = 5:7,
+           "A:B" = 8:10)
+summary(m0, split = list(X_full = ef))
+```
+
+NOTE: é importante enfatizar que a escolha que fizemos de hipóteses
+(cláusula 4) não corresponde a um conjunto de hipóteses lineares
+ortogonais.  Portanto, a estatística F da partição do quadro de análise
+de variância não será igual a estatística F do teste de Wald para alguma
+das hipóteses.
+
+```{r}
+# Produto interno nulo entre colunas -> hipóteses lineares ortogonais.
+C(da$A, contr = contr.helmert)
+C(da$A, contr = contr.treatment)
+C(da$A, contr = contr.SAS)
+C(da$A, contr = contr.poly)
+
+# Produto interno não nulo entre colunas.
+contrasts(da$A)
+```
+
+Conforme descrito no capítulo anterior, para modelos mais gerais não
+tem-se uma função que se comporta como a `aov()`.  Então reproduziremos
+os resultados usando o objeto de classe `lm()` empregado na
+`multcomp::glht()` e na `car::linearHypothesis()`.  O primeiro passo é
+ajustar o modelo com a `lm()`.  Lembre-se que os contrastes já foram
+definidos para os fatores A e B.
+
+```{r}
+# A `aov()` já abandona parâmetros não estimados do vetor. Mas para
+# modelos de outras classes, terá que ser feito manualmente.
+X_ok <- X_full[, str_replace(names(coef(m0)), "X_full", "")]
+
+# Modelo linear com matriz onde todos efeitos são estimáveis.
+m1 <- lm(y ~ 0 + X_ok, data = da)
+cbind(lm = coef(m1), aov = coef(m0))
+```
+
+A seguir, a estatística F do teste de Wald é obtida com a
+`multcomp::glht()` e com a `car::linearHypothesis()`.
+
+```{r}
+# `multcomp::cftest()` é aplicável a uma ampla classe de modelos.
+cftest(m1, parm = ef[[1]], test = Ftest()) # Fat. vs Neg.
+cftest(m1, parm = ef[[2]], test = Ftest()) # Fat. vs Pos.
+cftest(m1, parm = ef[[3]], test = Ftest()) # A.
+cftest(m1, parm = ef[[4]], test = Ftest()) # B.
+cftest(m1, parm = ef[[5]], test = Ftest()) # A:B.
+```
+
+```{r}
+# Matriz de hipóteses lineares (linear hypothesis matrix).
+m <- diag(length(coef(m1)))
+lhm <- ef
+for (i in 1:length(ef)) {
+    lhm[[i]] <- m[ef[[i]], , drop = FALSE]
+}
+lhm
+
+# `car::linearHypothesis()` é aplicável a uma ampla classe de modelos.
+linearHypothesis(m1, hypothesis.matrix = lhm[[1]], test = "F")
+linearHypothesis(m1, hypothesis.matrix = lhm[[2]], test = "F")
+linearHypothesis(m1, hypothesis.matrix = lhm[[3]], test = "F")
+linearHypothesis(m1, hypothesis.matrix = lhm[[4]], test = "F")
+linearHypothesis(m1, hypothesis.matrix = lhm[[5]], test = "F")
+```
+
+As comparações múltiplas de médias dentro de cada fator são feitas com
+os códigos a seguir.
+
+```{r}
+# Pontos experimentais únicos de A * B.
+i <- da %>% select(A, B) %>% duplicated()
+X_un <- X_ok[!i, ]         # Estimativas.
+pts <- da[!i, c("A", "B")] # Pontos experimentais.
+cbind(pts, X_un)
+
+all_pairwise <- function(L, collapse = "vs", row_names = rownames(L)) {
+    ij <- combn(x = 1:nrow(L), m = 2)
+    K <- apply(X = ij,
+               MARGIN = 2,
+               FUN = function(i) { L[i[1], ] - L[i[2], ] })
+    K <- t(K)
+    if (is.null(row_names)) {
+        nm <- 1:nrow(L)
+    } else {
+        nm <- row_names
+    }
+    rownames(K) <- apply(ij,
+                         MARGIN = 2,
+                         FUN = function(i) {
+                             paste(nm[i], collapse = collapse)
+                         })
+    return(K)
+}
+
+# Contrastes dois a dois (pairwise) entre níveis de A.
+L_A <- do.call(rbind, by(X_un, pts$A, colMeans))
+K_A <- all_pairwise(head(L_A, n = -2))
+summary(glht(m1, linfct = K_A))
+
+# Contrastes dois a dois (pairwise) entre níveis de B.
+L_B <- do.call(rbind, by(X_un, pts$B, colMeans))
+K_B <- all_pairwise(head(L_B, n = -2))
+summary(glht(m1, linfct = K_B))
+
+# Estudo da interação: entre níveis de B dentro de A.
+L_BinA <- by(X_un, pts$A, FUN = as.matrix)
+L_BinA <- lapply(L_BinA, "rownames<-", NULL)
+lapply(head(L_BinA, n = -2),
+       FUN = function(s) {
+           summary(glht(m1, linfct = all_pairwise(s)))
+       })
+
+# Estudo da interação: entre níveis de A dentro de B.
+L_AinB <- by(X_un, pts$B, FUN = as.matrix)
+L_AinB <- lapply(L_AinB, "rownames<-", NULL)
+lapply(head(L_AinB, n = -2),
+       FUN = function(s) {
+           summary(glht(m1, linfct = all_pairwise(s)))
+       })
+
+# Contraste do ponto adicional Negativo contra os pontos fatoriais.
+K <- rbind(colMeans(L_A[c("1", "2"), ]) - L_A["N", ])
+summary(glht(m1, linfct = K))
+
+# Contraste do ponto adicional Positivo contra os pontos fatoriais.
+K <- rbind(colMeans(L_A[c("1", "2"), ]) - L_A["P", ])
+summary(glht(m1, linfct = K))
+```
+
+Apenas para verificação, vamos obter as médias amostrais.
+
+```{r}
+# Diferenças entre médias amostrais apenas para verificação.
+da %>% group_by(A) %>% summarise(m = -mean(y)) %>%
+    pull(m) %>% head(n = -2) %>% outer(., ., "-") %>% as.dist()
+
+da %>% group_by(B) %>% summarise(m = -mean(y)) %>%
+    pull(m) %>% head(n = -2) %>% outer(., ., "-") %>% as.dist()
+
+x <- ifelse(da$A %in% head(levels(da$A), n = -2),
+            "fat",
+            as.character(da$A))
+da %>% group_by(!!x) %>%
+    summarise(m = -mean(y)) %>%
+    pull(m) %>% outer(., ., "-") %>% as.dist()
+```
+
+<!---------------------------------------------------------------------- -->
+
+## Análise de dados de *fusarium* em Zimmermann
+
+```{r}
+library(labestData)
+
+if (interactive()) {
+    help(ZimmermannTb9.22, help_type = "html")
+}
+
+# Cria a variável transformada pelo arco seno da raiz quadrada da
+# proporção amostral.
+tb <- ZimmermannTb9.22
+tb$y <- asin(sqrt(tb$nsi/40))
+str(tb)
+
+# Contagem dos pontos experimentais.
+xtabs(~fung + aplic, data = tb)
+
+# 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)
+```
+
+<!---------------------------------------------------------------------- -->
+
+
+```{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))
+
+# Importação dados dados.
+tb <- read_tsv("./data-raw/fat2x3+2adic.txt",
+               comment = "#",
+               col_types = cols())
+attr(tb, "spec") <- NULL
+str(tb)
+
+# Criando o fator NR.
+tb$nr <- "Baixo"
+tb$nr[with(tb, tratamento %in% c(1, 5))] <- "Convencional"
+
+# Contagem das ocorrências dos pontos experimentais.
+ftable(xtabs(~nr + fitato + fitase, data = tb))
+
+# Análise exploratória.
+ggplot(data = tb,
+       mapping = aes(x = fitase, y = cdaeb, color = nr)) +
+    facet_wrap(facets = ~fitato) +
+    geom_point() +
+    stat_summary(geom = "line", fun.y = "mean")
+```