diff --git a/40-analise-conjunta.Rmd b/40-analise-conjunta.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..a43f80e02b456ba45afdbc2f77d08262d711a499
--- /dev/null
+++ b/40-analise-conjunta.Rmd
@@ -0,0 +1,501 @@
+# Análise conjunta de experimentos
+
+**Atenção**: as análises nesse capítulo estão considerando apenas
+modelos de efeitos fixos. Posteriormente serão incluídas análises com
+modelos de efeitos mistos.
+
+```{r, include = FALSE, eval = FALSE}
+library(labestData)
+
+labestDataView()
+
+str(RamalhoEx8.1)
+str(BanzattoQd8.4.1)
+
+str(RamalhoEg8.1)
+str(RamalhoEg8.8)
+str(RamalhoEx8.2)
+str(RamalhoTb8.12)
+
+str(BanzattoQd8.2.1)
+str(BanzattoQd8.3.1)
+str(BanzattoQd8.3.1)
+str(BanzattoQd8.4.3)
+
+# DIC.
+rbind(ZimmermannTb12.1,
+ ZimmermannTb12.2)
+
+# DBC.
+rbind(ZimmermannTb12.7,
+ ZimmermannTb12.8)
+
+# DBC com alguns tratamentos comuns.
+rbind(ZimmermannTb12.19,
+ ZimmermannTb12.20)
+
+# DQL com um tratamento comum.
+rbind(ZimmermannTb12.26,
+ ZimmermannTb12.27)
+```
+
+```{r, message = FALSE}
+rm(list = objects())
+
+library(agricolae) # HSD.test().
+library(emmeans) # emmeans().
+library(multcomp) # glht() e métodos.
+library(tidyverse) # Manipulação e visualização.
+
+# Funções.
+source("mpaer_functions.R")
+```
+
+## Experimentos em delineamento inteiramente casualizado
+
+Os dados em `RamalhoEx8.1` são da produção de grãos (kg/ha) obtidos na
+avaliação de cultivares de arroz conduzidos em três locais no Estado de
+Minas Gerais. O experimento avaliou 10 cultivares em delineamento
+inteiramente casualizado com 3 repetições em cada local.
+
+### Análise exploratória
+
+```{r}
+# github.com/pet-estatistica/labestData/blob/devel/data-raw/RamalhoEx8.1.txt
+# help(RamalhoEx8.1, package = "labestData")
+da <- labestData::RamalhoEx8.1
+da$prod <- da$prod/1000
+da$rept <- factor(da$rept)
+str(da)
+
+xtabs(~local + cult, data = da)
+
+ggplot(data = da,
+ mapping = aes(x = cult, y = prod)) +
+ facet_wrap(facets = ~local, ncol = 1) +
+ geom_point() +
+ stat_summary(mapping = aes(group = 1),
+ geom = "line",
+ fun.y = "mean")
+```
+
+### Análises individuais
+
+```{r}
+da_local <- split(da, f = da$local)
+
+# Análise para cada local.
+m0_local <- lapply(da_local,
+ FUN = function(data) {
+ lm(prod ~ cult, data = data)
+ })
+
+# Gráficos de normalidade.
+par(mfrow = c(1, 3))
+sapply(m0_local, FUN = plot, which = 2)
+layout(1)
+
+# Quadrados médios.
+sapply(m0_local, FUN = function(m0) deviance(m0)/df.residual(m0))
+
+# Quadro de análise de variância.
+lapply(m0_local, FUN = anova)
+
+# Teste de Tukey.
+tk <- lapply(m0_local,
+ FUN = function(m0) {
+ tb <- HSD.test(m0, trt = "cult")$groups
+ tb$groups <- trimws(as.character(tb$groups))
+ cbind(cult = rownames(tb), tb)
+ })
+tk
+```
+
+### Análise conjunta
+
+A análise conjunta está considerando efeito fixo para local e cultivar.
+
+```{r, fig.height = 7}
+# Modelo com apenas um estrato para verificar pressupostos.
+m0 <- lm(terms(prod ~ local + local:rept + cult + local:cult,
+ keep.order = TRUE),
+ data = da)
+
+par(mfrow = c(2, 2))
+plot(m0)
+layout(1)
+
+# Modelo com dois estratos.
+m1 <- aov(prod ~ local + Error(local:rept) + cult + local:cult,
+ data = da)
+
+# Quadro de anova.
+summary(m1)
+
+# Médias marginais ajustadas por local.
+emm <- emmeans(m1, specs = ~cult | local)
+emm
+
+# Testes de Tukey por local da análise conjunta.
+# ATTENTION: assume balanceamento.
+tk <-
+ lapply(da_local,
+ FUN = function(data) {
+ rdf <- df.residual(m0)
+ rms <- deviance(m0)/rdf
+ tb <- HSD.test(y = data$prod,
+ trt = data$cult,
+ DFerror = rdf,
+ MSerror = rms)
+ tb <- tb$groups
+ tb$groups <- trimws(as.character(tb$groups))
+ names(tb)[1] <- "prod"
+ cbind(cult = rownames(tb), tb,
+ stringsAsFactors = FALSE)
+ })
+# tk
+
+# Junta para ter os intervalos de confiança.
+tk <- bind_rows(tk, .id = "local")
+tk <- inner_join(as.data.frame(emm), tk,
+ by = c("cult", "local"))
+tk[, c(1, 2, 8, 6:7, 9)]
+```
+
+```{r}
+# Gráfico de segmentos para as estimativas intervalares.
+ggplot(data = tk,
+ mapping = aes(x = prod, y = cult)) +
+ facet_wrap(facets = ~local) +
+ geom_point() +
+ geom_errorbarh(mapping = aes(xmin = lower.CL, xmax = upper.CL),
+ height = 0) +
+ geom_label(mapping = aes(label = sprintf("%0.2f%s", prod, groups)),
+ label.padding = unit(0.15, "lines"),
+ size = 3,
+ vjust = -0.25) +
+ expand_limits(y = c(NA, nlevels(da$cult) + 1)) +
+ labs(x = "Produção",
+ y = "Cultivares")
+```
+
+## Experimentos em delineamento de blocos completos
+
+Os dados em `BanzattoQd8.2.1` são referentes a um grupo de 4 ensaios de
+competição de 10 cultivares de batata que foram instalados no
+delineamento de blocos casualizados com 4 repetições. A variável
+resposta é a produção de tubérculos.
+
+### Análise exploratória
+
+```{r}
+# github.com/pet-estatistica/labestData/blob/devel/data-raw/BanzattoQd8.2.1.txt
+# help(BanzattoQd8.2.1, package = "labestData")
+da <- labestData::BanzattoQd8.2.1
+str(da)
+
+xtabs(~exper + cult, data = da)
+
+ggplot(data = da,
+ mapping = aes(x = cult, y = prod, color = bloco)) +
+ facet_wrap(facets = ~exper, ncol = 1) +
+ geom_point() +
+ geom_line(mapping = aes(group = bloco),
+ color = "gray") +
+ stat_summary(mapping = aes(group = 1),
+ geom = "line",
+ fun.y = "mean")
+```
+
+### Análises individuais
+
+```{r, fig.height = 7}
+da_exper <- split(da, f = da$exper)
+
+# Análise para cada exper.
+m0_exper <- lapply(da_exper,
+ FUN = function(data) {
+ lm(prod ~ bloco + cult, data = data)
+ })
+
+# Gráficos de normalidade.
+par(mfrow = c(2, 2))
+sapply(m0_exper, FUN = plot, which = 2)
+layout(1)
+
+# Quadrados médios.
+sapply(m0_exper, FUN = function(m0) deviance(m0)/df.residual(m0))
+
+# Quadro de análise de variância.
+lapply(m0_exper, FUN = anova)
+
+# Teste de Tukey.
+tk <- lapply(m0_exper,
+ FUN = function(m0) {
+ tb <- HSD.test(m0, trt = "cult")$groups
+ tb$groups <- trimws(as.character(tb$groups))
+ tb <- cbind(cult = rownames(tb), tb)
+ rownames(tb) <- NULL
+ return(tb)
+ })
+tk
+```
+
+### Análise conjunta
+
+A análise conjunta está considerando efeito fixo para experimento, bloco
+e cultivar.
+
+**Atenção**: esteja ciente de que a maior variância residual
+(experimento 4) é quase 5 vezes maior que a menor variância (experimento
+1).
+
+```{r, fig.height = 7}
+# Modelo com apenas um estrato para verificar pressupostos.
+m0 <- lm(terms(prod ~ exper + exper/bloco + cult + exper:cult,
+ keep.order = TRUE),
+ data = da)
+
+par(mfrow = c(2, 2))
+plot(m0)
+layout(1)
+
+# Modelo com dois estratos.
+m1 <- aov(prod ~ exper + Error(exper:bloco) + cult + exper:cult,
+ data = da)
+
+# Quadro de anova.
+summary(m1)
+
+# Médias marginais ajustadas por exper.
+emm <- emmeans(m1, specs = ~cult | exper)
+emm
+
+# Testes de Tukey por exper da análise conjunta.
+# ATTENTION: assume balanceamento.
+tk <-
+ lapply(da_exper,
+ FUN = function(data) {
+ rdf <- df.residual(m0)
+ rms <- deviance(m0)/rdf
+ tb <- HSD.test(y = data$prod,
+ trt = data$cult,
+ DFerror = rdf,
+ MSerror = rms)
+ tb <- tb$groups
+ tb$groups <- trimws(as.character(tb$groups))
+ names(tb)[1] <- "prod"
+ cbind(cult = rownames(tb), tb,
+ stringsAsFactors = FALSE)
+ })
+# tk
+
+# Junta para ter os intervalos de confiança.
+tk <- bind_rows(tk, .id = "exper")
+tk <- inner_join(as.data.frame(emm), tk,
+ by = c("cult", "exper"))
+tk[, c(1, 2, 8, 6:7, 9)] %>%
+ arrange(exper, desc(prod))
+```
+
+```{r}
+# NOTE: para ordenar os níveis de cultivar dentro de cada experimento.
+tk <- tk %>%
+ mutate(combo = {
+ f <- interaction(exper, cult, drop = TRUE)
+ reorder(f, prod)
+ })
+
+# Gráfico de segmentos para as estimativas intervalares.
+ggplot(data = tk,
+ mapping = aes(x = combo, y = prod)) +
+ facet_wrap(facets = ~exper,
+ nrow = 1,
+ scales = "free_x") +
+ geom_point() +
+ geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL),
+ width = 0) +
+ scale_x_discrete(breaks = levels(tk$combo),
+ labels = gsub("^.*\\.", "", levels(tk$combo))) +
+ geom_text(mapping = aes(label = sprintf("%0.1f %s", prod, groups)),
+ size = 3,
+ angle = 90,
+ vjust = -0.5) +
+ labs(y = "Produção",
+ x = "Cultivares") +
+ expand_limits(x = c(-0.5, NA)) +
+ theme(axis.text.x = element_text(angle = 90,
+ vjust = 0.5,
+ hjust = 1))
+```
+
+## Experimentos com tratamentos comuns
+
+```{r}
+# github.com/pet-estatistica/labestData/blob/devel/data-raw/ZimmermannTb12.19.txt
+# help(ZimmermannTb12.19, package = "labestData")
+
+da <- rbind(cbind(labestData::ZimmermannTb12.19, exper = 1),
+ cbind(labestData::ZimmermannTb12.20, exper = 2))
+da$exper <- factor(da$exper)
+da$prod <- da$prod/1000
+str(da)
+
+xtabs(~cult + exper, data = da)
+
+ggplot(data = da,
+ mapping = aes(x = cult, y = prod, color = bloco)) +
+ facet_wrap(facets = ~exper, ncol = 1) +
+ geom_point() +
+ geom_line(mapping = aes(group = bloco),
+ color = "gray") +
+ stat_summary(mapping = aes(group = 1),
+ geom = "line",
+ fun.y = "mean") +
+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
+```
+
+### Análises individuais
+
+```{r, fig.height = 7}
+da_exper <- split(da, f = da$exper)
+
+# Análise para cada exper.
+m0_exper <- lapply(da_exper,
+ FUN = function(data) {
+ lm(prod ~ bloco + cult, data = data)
+ })
+
+# Gráficos de normalidade.
+par(mfrow = c(1, 2))
+sapply(m0_exper, FUN = plot, which = 2)
+layout(1)
+
+# Quadrados médios.
+sapply(m0_exper, FUN = function(m0) deviance(m0)/df.residual(m0))
+
+# Quadro de análise de variância.
+lapply(m0_exper, FUN = anova)
+
+# Teste de Tukey.
+tk <- lapply(m0_exper,
+ FUN = function(m0) {
+ tb <- HSD.test(m0, trt = "cult")$groups
+ tb$groups <- trimws(as.character(tb$groups))
+ tb <- cbind(cult = rownames(tb), tb)
+ rownames(tb) <- NULL
+ return(tb)
+ })
+tk
+```
+
+### Análise conjunta
+
+A análise conjunta está considerando efeito fixo para experimento, bloco
+e cultivar.
+
+```{r, fig.height = 7}
+# Modelo com apenas um estrato para verificar pressupostos.
+m0 <- lm(terms(prod ~ exper + exper/bloco + cult + exper:cult,
+ keep.order = TRUE),
+ data = da)
+
+par(mfrow = c(2, 2))
+plot(m0)
+layout(1)
+
+# Modelo com dois estratos.
+m1 <- aov(prod ~ exper + Error(exper:bloco) + cult + exper:cult,
+ data = da)
+
+# Quadro de anova.
+summary(m1)
+
+# Abandona a interação não significativa para ter todos os efeitos
+# estimados.
+m1 <- aov(prod ~ exper + Error(exper:bloco) + cult,
+ data = da)
+
+# Médias marginais ajustadas do modelo de 2 estratos.
+# emm <- emmeans(m1, specs = ~cult)
+# as.data.frame(emm)
+
+# Volta para o modelo de um estrato.
+m0 <- lm(prod ~ exper + exper:bloco + cult,
+ data = da)
+
+emm <- emmeans(m0, specs = ~cult)
+as.data.frame(emm)
+
+grid <- attr(emm, "grid")
+L <- attr(emm, "linfct")
+rownames(L) <- grid[[1]]
+
+# glht(m0, linfct = mcp(cult = "Tukey"))
+results <- wzRfun::apmc(L,
+ model = m0,
+ focus = "cult",
+ test = "fdr") %>%
+ mutate(cld = wzRfun::ordered_cld(let = cld, means = fit))
+results
+```
+
+```{r}
+# Gráfico de segmentos para as estimativas intervalares.
+ggplot(data = results,
+ mapping = aes(x = reorder(cult, fit), y = fit)) +
+ geom_point() +
+ geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
+ width = 0) +
+ geom_text(mapping = aes(label = sprintf("%0.2f %s", fit, cld)),
+ angle = 90,
+ vjust = -0.5) +
+ labs(y = "Produção",
+ x = "Cultivares") +
+ expand_limits(x = c(-0.5, NA)) +
+ theme(axis.text.x = element_text(angle = 90,
+ vjust = 0.5,
+ hjust = 1))
+```
+
+## Experimentos ao longo do tempo
+
+TODO as ser feito: experimento em várias safras. E em vários locais e
+safras.
+
+ATTENTION: essa análise precisa ser refeita nos moldes de uma análise de
+experimento de medidas repetidas pois, visto que a cultura é perene, a
+mesma unidade experimental é medida ao longo do tempo. Se fosse uma
+cultura anual, cada safra seria um novo cultivo onde as unidades
+poderiam ser todas diferentes e o experimento ser em outras
+instalações/condições. Poderia ser considerado medida repetida apenas se
+as mesmas parcelas recebessem os mesmos tratamentos ao longo do tempo na
+mesma área experimental.
+
+```{r}
+# github.com/pet-estatistica/labestData/blob/devel/data-raw/RamalhoTb8.12.txt
+# help(RamalhoTb8.12, package = "labestData")
+da <- labestData::RamalhoTb8.12
+str(da)
+
+xtabs(~ano + prog, data = da)
+
+ggplot(data = da,
+ mapping = aes(x = prog, y = prod, color = bloc)) +
+ facet_wrap(facets = ~ano, ncol = 1) +
+ geom_point() +
+ geom_line(mapping = aes(group = bloc),
+ color = "gray") +
+ stat_summary(mapping = aes(group = 1),
+ geom = "line",
+ fun.y = "mean")
+
+# NOTE: o comportamento bianual de produção de café.
+ggplot(data = da,
+ mapping = aes(x = ano, y = prod)) +
+ facet_wrap(facets = ~prog, nrow = 1) +
+ geom_point() +
+ stat_summary(mapping = aes(group = 1),
+ geom = "line",
+ fun.y = "mean")
+```