From 19996b023d9c31fd2718fd4d798df9f4062303d3 Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Tue, 21 Jan 2020 18:53:03 -0300
Subject: [PATCH] Adds a dataset and a preliminary analysis.

---
 50-parcela-subdividida.Rmd | 230 +++++++++++++++++++++++++++++++++++++
 1 file changed, 230 insertions(+)

diff --git a/50-parcela-subdividida.Rmd b/50-parcela-subdividida.Rmd
index 2fc743c..914d7ec 100644
--- a/50-parcela-subdividida.Rmd
+++ b/50-parcela-subdividida.Rmd
@@ -14,3 +14,233 @@ knitr::include_graphics("./img/delineamento-parcela-subdividida-plantio.png")
 str(PimentelTb9.3.1)
 help(PimentelTb9.3.1, h = "html")
 ```
+
+<!------------------------------------------- -->
+
+```{r}
+# Dados fornecidos pelo Eng. Agronômo Jenilton Gomes da Cunha, ex-aluno
+# de Mestrado em Produção Vegetal da Universidade Federal do Vale do São
+# Francisco em Petrolina.
+
+txt <- "
+DOSES  	DATAS	BLOCO	CST	Amido
+0	0	1	187,46	6,144
+0	0	2	230,43	3,137
+0	0	3	125,67	4,095
+0	0	4	198,65	4,459
+0,7	0	1	161,99	7,191
+0,7	0	2	163,13	1,738
+0,7	0	3	137,34	1,226
+0,7	0	4	200,59	3,385
+1	0	1	188,92	4,287
+1	0	2	194,75	4,543
+1	0	3	220,86	2,654
+1	0	4	203,19	4,945
+1,3	0	1	170,10	3,905
+1,3	0	2	160,86	1,978
+1,3	0	3	166,86	5,570
+1,3	0	4	236,86	1,403
+PBZ_AD	0	1	182,75	3,589
+PBZ_AD	0	2	190,21	3,179
+PBZ_AD	0	3	220,97	4,702
+PBZ_AD	0	4	277,46	5,369
+0	30	1	229,86	4,817
+0	30	2	219,83	2,244
+0	30	3	219,64	3,699
+0	30	4	113,88	3,992
+0,7	30	1	228,53	3,165
+0,7	30	2	316,97	4,238
+0,7	30	3	224,94	2,111
+0,7	30	4	272,24	3,559
+1	30	1	118,80	2,897
+1	30	2	211,51	2,837
+1	30	3	274,70	3,983
+1	30	4	149,26	3,528
+1,3	30	1	185,40	3,197
+1,3	30	2	197,51	2,505
+1,3	30	3	189,94	2,870
+1,3	30	4	197,13	2,138
+PBZ_AD	30	1	139,04	2,963
+PBZ_AD	30	2	221,91	2,745
+PBZ_AD	30	3	285,11	2,065
+PBZ_AD	30	4	224,37	4,488
+0	60	1	120,36	3,514
+0	60	2	170,17	2,816
+0	60	3	198,48	2,668
+0	60	4	177,38	1,976
+0,7	60	1	180,75	3,320
+0,7	60	2	267,32	2,289
+0,7	60	3	193,80	4,884
+0,7	60	4	262,94	2,625
+1	60	1	140,49	3,779
+1	60	2	231,89	4,618
+1	60	3	220,99	3,763
+1	60	4	246,86	2,956
+1,3	60	1	283,37	1,895
+1,3	60	2	200,89	3,500
+1,3	60	3	315,89	5,350
+1,3	60	4	253,31	3,064
+PBZ_AD	60	1	245,54	9,155
+PBZ_AD	60	2	231,58	1,769
+PBZ_AD	60	3	357,19	3,310
+PBZ_AD	60	4	191,34	2,280
+0	90	1	228,53	4,711
+0	90	2	189,17	4,420
+0	90	3	131,14	5,843
+0	90	4	211,95	3,870
+0,7	90	1	195,90	4,830
+0,7	90	2	166,70	3,396
+0,7	90	3	121,64	5,933
+0,7	90	4	164,91	4,404
+1	90	1	157,74	6,220
+1	90	2	160,38	3,436
+1	90	3	196,67	4,400
+1	90	4	181,36	4,018
+1,3	90	1	143,96	3,474
+1,3	90	2	172,00	4,819
+1,3	90	3	196,72	2,852
+1,3	90	4	79,47	2,844
+PBZ_AD	90	1	190,74	4,803
+PBZ_AD	90	2	154,68	4,488
+PBZ_AD	90	3	148,95	2,914
+PBZ_AD	90	4	174,88	2,817"
+da <- read.table(textConnection(txt),
+                 header = TRUE,
+                 sep = "\t",
+                 dec = ",",
+                 stringsAsFactors = FALSE)
+names(da) <- str_to_lower(names(da))
+
+da <- da %>%
+    mutate(bloco = factor(bloco),
+           trat = str_replace(doses, ",", "."),
+           testem = as.integer(doses == "PBZ_AD"),
+           doses = ifelse(as.logical(testem), "0", doses),
+           doses = parse_double(doses, locale = locale(decimal_mark = ","))) %>%
+    select(bloco, trat, datas, cst, amido, testem, doses)
+str(da)
+```
+
+```{r}
+# Análise exploratória.
+ggplot(data = da,
+       mapping = aes(x = datas, y = cst)) +
+    facet_wrap(facets = ~ testem + doses, nrow = 1) +
+    geom_point(alpha = 0.5) +
+    stat_summary(mapping = aes(group = 1), geom = "line", fun.y = "mean")
+
+ggplot(data = da,
+       mapping = aes(x = doses, y = cst)) +
+    facet_grid(facets = datas ~ testem, scale = "free_x", space = "free") +
+    geom_point(alpha = 0.5) +
+    stat_summary(mapping = aes(group = 1), geom = "line", fun.y = "mean")
+```
+
+```{r}
+# Pacotes para a análise dos dados.
+library(lme4)
+library(lmerTest)
+library(emmeans)
+
+# Converte e cria variáveis para usar no modelo.
+da <- da %>%
+    mutate(Datas = factor(datas),
+           Trat = factor(trat),
+           Test = factor(testem),
+           ue = interaction(trat, bloco, sep = "_", drop = TRUE))
+
+#-----------------------------------------------------------------------
+# Modelo de parcela subdividida. Todos os fatores são qualitativos.
+
+m0 <- lmer(cst ~ bloco + Trat * Datas + (1 | ue), data = da)
+# par(mfrow = c(2, 2)); plot(m0); layout(1)
+VarCorr(m0)
+anova(m0)
+
+# Comparações múltiplas entre datas.
+cld(emmeans(m0, ~Datas))
+
+#-----------------------------------------------------------------------
+# Usando efeito de segundo grau para as doses.
+
+m1 <- update(m0, . ~ bloco +
+                     (Test + poly(doses, degree = 2)) *
+                     Datas +
+                     (1 | ue))
+anova(m1, m0)
+anova(m1)
+
+# ATTENTION: O modelo mostra que apenas o efeito de datas é relevante.
+
+# Grid de valores para fazer a predição e gráfico.
+grid <- with(da,
+             crossing(doses_u = c(-0.25, seq(0, 1.5, length.out = 13)),
+                      Datas = unique(Datas))) %>%
+    mutate(Test = factor(as.integer(doses_u < 0)),
+           doses = ifelse(doses_u < 0, 0, doses_u),
+           bloco = unique(da$bloco)[1])
+str(grid)
+
+# ftable(xtabs(~Test + doses + Datas, data = grid))
+
+# Matrix do modelo para a parte fixa.
+X <- model.matrix(nobars(formula(m1))[-2],
+                  data = grid)
+dim(X)
+
+# Nivela o efeito de blocos.
+X[, attr(X, "assign") == 1] <- 1/nlevels(da$bloco)
+X[1:5, 1:5]
+
+# Faz a predição.
+grid$y <- c(X %*% fixef(m1))
+
+# Exibe os resultados.
+ggplot(data = grid,
+       mapping = aes(x = doses_u, y = y, color = Test)) +
+    facet_grid(facets = . ~ Datas) +
+    geom_point() +
+    geom_line() +
+    geom_point(data = da,
+               mapping = aes(x = doses - 0.25 * (Test == "1"),
+                             y = cst),
+               pch = 1)
+
+#-----------------------------------------------------------------------
+# Modelo com efeitos aditivos.
+
+m2 <- update(m0, . ~ bloco +
+                     Test + poly(doses, degree = 2) +
+                     Datas +
+                     (1 | ue))
+anova(m2, m1)
+anova(m2)
+
+# ATTENTION: novamente apenas há o efeito de datas.
+
+# Matrix do modelo para a parte fixa.
+X <- model.matrix(nobars(formula(m2))[-2],
+                  data = grid)
+colnames(X)
+
+# Nivela o efeito de blocos.
+X[, attr(X, "assign") == 1] <- 1/nlevels(da$bloco)
+X[1:5, 1:5]
+
+# Faz a predição.
+grid$y <- c(X %*% fixef(m2))
+
+# Exibe os resultados.
+ggplot(data = grid,
+       mapping = aes(x = doses_u, y = y, color = Test)) +
+    facet_grid(facets = . ~ Datas) +
+    geom_point() +
+    geom_line() +
+    geom_point(data = da,
+               mapping = aes(x = doses - 0.25 * (Test == "1"),
+                             y = cst),
+               pch = 1)
+
+# Comparação múltiplas entre as datas.
+cld(emmeans(m2, specs = ~Datas))
+```
-- 
GitLab