Produtividade como função do teor de nutrientes do solo
Walmes Zeviani & Fábio Ono
2019-04-30
1 Importação e normalização dos dados
#-----------------------------------------------------------------------
# Definições da sessão e pacotes usados.
rm(list = objects())
library(latticeExtra)
library(tidyverse)
# Escolhe o tema para ggplot2.
# theme_set(theme_bw())
# theme_set(theme_classic())
theme_set(theme_light())
# Para não mostrar a mensagem de parse das colunas ao importar.
options(readr.num_columns = 0)
library(lme4)
library(lmerTest)
library(doBy)
library(wzRfun)
sessionInfo()## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods stats graphics grDevices datasets utils base
##
## other attached packages:
## [1] wzRfun_0.81 doBy_4.6-2 lmerTest_3.0-1
## [4] lme4_1.1-18-1 Matrix_1.2-14 forcats_0.3.0
## [7] stringr_1.4.0 dplyr_0.8.0.1 purrr_0.3.2
## [10] readr_1.1.1 tidyr_0.8.3 tibble_2.1.1
## [13] ggplot2_3.1.1 tidyverse_1.2.1 latticeExtra_0.6-28
## [16] RColorBrewer_1.1-2 lattice_0.20-38 rmarkdown_1.10
## [19] knitr_1.20.12
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 lubridate_1.7.4 mvtnorm_1.0-8
## [4] zoo_1.8-4 assertthat_0.2.1 rprojroot_1.3-2
## [7] digest_0.6.18 R6_2.4.0 cellranger_1.1.0
## [10] plyr_1.8.4 backports_1.1.2 evaluate_0.12
## [13] httr_1.3.1 pillar_1.3.1 rlang_0.3.4
## [16] lazyeval_0.2.2 multcomp_1.4-8 readxl_1.1.0
## [19] rstudioapi_0.8 minqa_1.2.4 nloptr_1.2.1
## [22] splines_3.4.4 munsell_0.5.0 broom_0.5.0
## [25] compiler_3.4.4 numDeriv_2016.8-1 modelr_0.1.2
## [28] xfun_0.4 pkgconfig_2.0.2 tcltk_3.4.4
## [31] htmltools_0.3.6 tidyselect_0.2.5 bookdown_0.9
## [34] codetools_0.2-15 crayon_1.3.4 withr_2.1.2
## [37] MASS_7.3-51.1 grid_3.4.4 nlme_3.1-137
## [40] jsonlite_1.5 gtable_0.3.0 magrittr_1.5
## [43] scales_1.0.0 cli_1.1.0 stringi_1.4.3
## [46] xml2_1.2.0 rpanel_1.1-4 sandwich_2.5-0
## [49] TH.data_1.0-9 tools_3.4.4 glue_1.3.1
## [52] hms_0.4.2 survival_2.43-1 yaml_2.2.0
## [55] colorspace_1.4-1 rvest_0.3.2 haven_1.1.2#-----------------------------------------------------------------------
# Importação dos dados.
# O melhor jeito de ler é sempre por CSV.
da <- read_csv2(file = "dados.csv",
na = c("-", " -", "----"),
skip = 1,
locale = locale(decimal_mark = ","),
col_names = TRUE)
# str(da, give.attr = FALSE)
#-----------------------------------------------------------------------
# Normalização da tabela de dados.
# Padroniza nomes das variáveis.
names(da) <- iconv(x = substr(tolower(x = names(da)),
start = 1,
stop = 4),
to = "ASCII//TRANSLIT")
# str(da, give.attr = FALSE)
# Obtém as coordenadas da unidades experimentais e cria fatores.
da <- da %>%
mutate(coorx = as.numeric(gsub("^.([0-9.]+),.*$", "\\1", coor)),
coory = as.numeric(gsub("^.*,([0-9.]+).$", "\\1", coor)),
cult = sub("^[^ ]* ", "", cult),
corr = factor(corr),
bloc = factor(bloc))
# str(da, give.attr = FALSE)
# Criando fator que representa a aplicação com fósforo.
fosf <- data.frame(coorx = sort(unique(da$coorx)),
Padic = rep(c(0, 50, 100, 50, 100),
times = 3),
Papli = rep(c("sem", "sulc", "sulc", "lanç", "lanç"),
times = 3))
fosf$Pman <- with(fosf, paste(Padic, Papli, sep = "*"))
# Incorporando a tabela principal (left join).
da <- da %>%
left_join(y = fosf) %>%
mutate(Papli = factor(Papli, levels = unique(fosf$Papli)[c(1, 3, 2)]))
# str(da)
# Decomponto rótulo de correção do solo em mais fatores (fatorial
# incompleto).
comp <- data.frame(corr = levels(da$corr),
fos = c(200, 200, 200, 200, 200, 600, 0),
cam = c(20, 20, 40, 0, 0, 40, 0),
ges = c(200, 3250, 200, 200, 3250, 200, 200))
comp <- comp %>%
arrange(fos, cam, ges)
comp$trat <- with(comp, paste(fos, cam, ges, sep = "*"))
# str(comp)
# Incorporando na tabela principal (left join).
da <- da %>%
left_join(y = comp) %>%
mutate(trat = factor(trat, levels = comp$trat))
# Nomenclatura longa.
relev <- list(c("0*0*200" = "200 gesso",
"200*0*200" = "200 gesso + 200 P 0cm",
"200*0*3250" = "3250 gesso + 200 P 0cm",
"200*20*200" = "200 gesso + 200 P 20cm",
"200*20*3250" = "3250 gesso + 200 P 20cm",
"200*40*200" = "200 gesso + 200 P 40cm",
"600*40*200" = "200 gesso + 600 P 40cm"),
# c("0*sem" = "0 ~ kg/ha ~ P[2] * O[5]",
# "50*lanç" = "50 ~ kg/ha ~ P[2] * O[5] ~ lanço",
# "50*sulc" = "50 ~ kg/ha ~ P[2] * O[5] ~ sulco",
# "100*lanç" = "100 ~ kg/ha ~ P[2] * O[5] ~ lanço",
# "100*sulc" = "100 ~ kg/ha ~ P[2] * O[5] ~ sulco")
c("0*sem" = "0",
"50*lanç" = "50 - lanço",
"50*sulc" = "50 - sulco",
"100*lanç" = "100 - lanço",
"100*sulc" = "100 - sulco"))
# Nomenclatura curta.
# relev[[1]][] <- sprintf("F%d", 1:7)
relev[[1]][] <- sprintf("N%d", 1:7)
# Renomeia fatores para algo mais descritivo.
da <- da %>%
rename("correc" = "trat", "semead" = "Pman") %>%
mutate(correc = factor(correc,
levels = names(relev[[1]]),
labels = relev[[1]]),
semead = factor(semead,
levels = names(relev[[2]]),
labels = relev[[2]])) %>%
select(safr, cult, bloc, coorx, coory, prod, semead, correc, cam, fos, ges)
str(da)## Classes 'tbl_df', 'tbl' and 'data.frame': 1155 obs. of 11 variables:
## $ safr : chr "2010/2011" "2010/2011" "2010/2011" "2010/2011" ...
## $ cult : chr "Soja" "Soja" "Soja" "Soja" ...
## $ bloc : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
## $ coorx : num 0 31.5 63 18.9 50.4 81.9 25.2 56.7 88.2 6.3 ...
## $ coory : num 69 69 69 69 69 69 69 69 69 69 ...
## $ prod : num 1959 2093 1893 2607 2455 ...
## $ semead: Factor w/ 5 levels "0","50 - lanço",..: 1 1 1 2 2 2 4 4 4 3 ...
## $ correc: Factor w/ 7 levels "N1","N2","N3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ cam : num 0 0 0 0 0 0 0 0 0 0 ...
## $ fos : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ges : num 200 200 200 200 200 200 200 200 200 200 ...2 Análise exploratória
#-----------------------------------------------------------------------
# Análise exploratória.
# Safras e cultura usada.
da %>% select(safr, cult) %>% unique()## # A tibble: 11 x 2
## safr cult
## <chr> <chr>
## 1 2010/2011 Soja
## 2 2011/2012 Soja
## 3 2012/2013 Soja
## 4 2013/2014 Soja
## 5 2013/2014 Milho Safrinha
## 6 2014/2015 Soja
## 7 2014/2015 Milho Safrinha
## 8 2015/2016 Soja
## 9 2015/2016 Milho Safrinha
## 10 2016/2017 Soja
## 11 2016/2017 Milho Safrinha# Análise apenas para o milho safrinha.
db <- da %>% filter(cult == "Milho Safrinha") %>% droplevels()
# Layout do experimento com a produção média de milho.
db %>%
group_by(coorx, coory) %>%
summarise(prod = mean(prod)) %>%
ggplot(mapping = aes(x = coorx, y = coory)) +
geom_rect(mapping = aes(xmin = coorx - 3.1,
xmax = coorx + 3.1,
ymin = coory - 4.95,
ymax = coory + 4.95,
fill = prod),
color = "white") +
scale_fill_gradient(low = "grey10", high = "gray90") +
geom_point() +
coord_equal()
ggplot(data = db,
mapping = aes(x = correc,
y = prod,
shape = factor(cam),
color = sprintf("%d kg + %d kg", fos, ges))) +
geom_point() +
stat_summary(mapping = aes(group = 1),
fun.y = mean,
colour = "gray30",
geom = "line",
group = 1) +
labs(color = expression(P[2] * O[5] + "Gesso"),
shape = "Profundidade (cm)") +
# xlab("Correção do solo\n fósforo (kg/ha) * camada (cm) * gesso (kg/ha)") +
xlab("Correção do solo") +
ylab(expression("Produtividade de grãos" ~ (kg ~ ha^{-1}))) +
facet_grid(facets = safr ~ semead, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8),
legend.position = "top") +
guides(color = guide_legend(nrow = 2,
title.position = "top"),
shape = guide_legend(nrow = 1,
title.position = "top"))
3 Análise separado por safra
Análise por safra assumindo delineamento de blocos em faixas completamente casualizadas. Na realidade a casualização não foi aleatória mas sim completamente sistemática. Por outro lado, esse é o delineamento que melhor pode aproximar do (quase) delineamento usado para o experimento. Talvez análise considerando a coordenada das parcelas possa amenizar a falta de casualização corrigindo para alguma eventual dependência especial.
Figure 3.1: Imagem aérea do experimento com retângulos delimitando os blocos (azul), as faixas verticais (amarelo) e as faixas horizontais (rosa).
Figure 3.2: Croqui do experimento que exibe a alocação dos tratamentos aplicados. Nas faixas horizontais foi casualizado um fator de 7 níveis. Nas faixas verticais foi casualizado um fator de 5 níveis. A casualização foi a mesma em todos os blocos mas em um típico experimento em faixas, em cada bloco é feita a casualização independente.
Os dados serão sumetidos para ajuste do modelo estatístico correspondende ao ensaio com de arranjo fatorial duplo completo em faixas conduzido no delineamento de blocos casualizados. Após o ajuste, a significância do termo de interação bem como dos termos para os efeitos principais serão avaliados pela estatística F do teste de Wald ao nível de 5%. Quando a interação for significativa, serão feitas comparações múltiplas entre os níveis de um fator fixando o nível do outro. Não havendo interação, será estudado cada fator separadamente. O nível de significância das comparações múltiplas será ajustado pelo método false discovery rate.
# Lista para armazenar as médias em cada safra.
u <- unique(db$safr)
resul_all <- setNames(vector(mode = "list", length = length(u)), nm = u)3.1 Safra 2013/2014
#-----------------------------------------------------------------------
# Modelo para o delineamento de blocos casualizados. Os blocos são de
# tamanho 7 * 5 = 35.
m0 <- lmer(prod ~
(1 | bloc) +
correc * semead,
data = subset(db, safr == "2013/2014"))
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 100613277 16768880 6 68 31.9299 < 2.2e-16 ***
## semead 256851189 64212797 4 68 122.2685 < 2.2e-16 ***
## correc:semead 49996137 2083172 24 68 3.9666 4.084e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo para o experimento em faixas casualizadas em blocos.
m1 <- update(m0,
formula = . ~
(1 | bloc) + # bloco.
(1 | bloc:correc) + # faixa horizontal dentro de bloco.
(1 | bloc:semead) + # faixa vertical dentro de bloco.
correc * semead)
anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 100613277 16768880 6 60 32.2461 < 2.2e-16 ***
## semead 236904236 59226059 4 8 113.8902 4.374e-07 ***
## correc:semead 49996137 2083172 24 60 4.0059 6.712e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m1)
#-----------------------------------------------------------------------
# Desdobramento da interação com comparações múltiplas.
# Para obter as médias ajustadas.
lsm <- LSmeans(m1, effect = c("correc", "semead"))
# lsm
# Cria nomes para as linhas da matriz.
rownames(lsm$L) <-
lsm$grid %>%
unite(row, correc, semead) %>%
pull(row)
# Matrizes para estudar o efeito da aplicação de fósforo (5 níveis)
# fixando o manejo do gesso (7 níveis).
s_semead <- by(data = lsm$L,
INDICES = lsm$grid[, "correc"],
FUN = as.matrix)
s_semead <- lapply(s_semead,
FUN = function(x) {
rownames(x) <- sub("^.*_", "", rownames(x))
return(x)
})
# Matrizes para estudar o efeito do manejo do gesso (7 níveis) fixando a
# aplicação de fósforo (5 níveis).
s_correc <- by(data = lsm$L,
INDICES = lsm$grid[, "semead"],
FUN = as.matrix)
s_correc <- lapply(s_correc,
FUN = function(x) {
rownames(x) <- sub("_.*$", "", rownames(x))
return(x)
})
# Faz as comparações múltiplas e acerta disposição das letras.
c_semead <- sapply(s_semead,
FUN = apmc,
model = m0,
focus = "semead",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_semead <- sapply(c_semead,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
c_correc <- sapply(s_correc,
FUN = apmc,
model = m0,
focus = "correc",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_correc <- sapply(c_correc,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
# Conversão de lista para tabela.
c_semead <- plyr::ldply(c_semead, .id = "correc")
c_correc <- plyr::ldply(c_correc, .id = "semead")
# Faz a junção e coloca letras altas e baixas.
resul <- c_semead %>%
mutate(cle = toupper(cld)) %>%
select(semead, correc, cle) %>%
inner_join(y = c_correc, by = c("semead", "correc")) %>%
unite(clde, cle, cld, sep = "", remove = FALSE)
resul <- equallevels(resul, db)
resul_all[[1]] <- resul
str(resul)## 'data.frame': 35 obs. of 8 variables:
## $ semead: Factor w/ 5 levels "0","50 - lanço",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ correc: Factor w/ 7 levels "N1","N2","N3",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ clde : chr "Ce" "Bb" "ACc" "ACa" ...
## $ cle : chr "C" "B" "AC" "AC" ...
## $ fit : num 473 4545 5586 7427 7510 ...
## $ lwr : num -360 3712 4753 6594 6677 ...
## $ upr : num 1305 5378 6419 8260 8343 ...
## $ cld : chr "e" "b" "c" "a" ...ggplot(data = resul,
mapping = aes(x = correc, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, cld)),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~semead, nrow = 1, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Correção inicial do solo com combinações de fósforo, camada e gesso",
# "[fósforo (kg/ha) * camada (cm) * gesso (kg/ha)]"))) +
xlab("Correção inicial do solo com combinações de fósforo, camada e gesso") +
ylab("Produtividade de grãos (kg/ha)")
ggplot(data = resul,
mapping = aes(x = semead, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, tolower(cle))),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~correc, nrow = 1) +
scale_x_discrete(labels = function(l) parse(text = l)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Aplicações de" ~ P[2] * O[5] ~ "na semeadura com combinação da quantidade e forma de aplicação",
# "[" * P[2] * O[5] ~ "(kg/ha) * forma de aplicação]"))) +
xlab("Aplicação de fósforo na semeadura") + xlab(expression("Doses de" ~ P[2] * O[5] ~ "na semeadura (kg/ha)")) +
ylab("Produtividade de grãos (kg/ha)")
#-----------------------------------------------------------------------3.2 Safra 2014/2015
#-----------------------------------------------------------------------
# Modelo para o delineamento de blocos casualizados. Os blocos são de
# tamanho 7 * 5 = 35.
m0 <- lmer(prod ~
(1 | bloc) +
correc * semead,
data = subset(db, safr == "2014/2015"))
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 14548854 2424809 6 70 17.2308 4.016e-12 ***
## semead 173336358 43334089 4 70 307.9336 < 2.2e-16 ***
## correc:semead 27462915 1144288 24 70 8.1314 2.928e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo para o experimento em faixas casualizadas em blocos.
m1 <- update(m0,
formula = . ~
(1 | bloc) + # bloco.
(1 | bloc:correc) + # faixa horizontal dentro de bloco.
(1 | bloc:semead) + # faixa vertical dentro de bloco.
correc * semead)
anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 14548854 2424809 6 60 21.338 3.108e-13 ***
## semead 64951580 16237895 4 10 142.894 8.889e-09 ***
## correc:semead 27462915 1144288 24 60 10.070 2.743e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m1)
#-----------------------------------------------------------------------
# Desdobramento da interação com comparações múltiplas.
# Para obter as médias ajustadas.
lsm <- LSmeans(m1, effect = c("correc", "semead"))
# lsm
# Cria nomes para as linhas da matriz.
rownames(lsm$L) <-
lsm$grid %>%
unite(row, correc, semead) %>%
pull(row)
# Matrizes para estudar o efeito da aplicação de fósforo (5 níveis)
# fixando o manejo do gesso (7 níveis).
s_semead <- by(data = lsm$L,
INDICES = lsm$grid[, "correc"],
FUN = as.matrix)
s_semead <- lapply(s_semead,
FUN = function(x) {
rownames(x) <- sub("^.*_", "", rownames(x))
return(x)
})
# Matrizes para estudar o efeito do manejo do gesso (7 níveis) fixando a
# aplicação de fósforo (5 níveis).
s_correc <- by(data = lsm$L,
INDICES = lsm$grid[, "semead"],
FUN = as.matrix)
s_correc <- lapply(s_correc,
FUN = function(x) {
rownames(x) <- sub("_.*$", "", rownames(x))
return(x)
})
# Faz as comparações múltiplas e acerta disposição das letras.
c_semead <- sapply(s_semead,
FUN = apmc,
model = m0,
focus = "semead",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_semead <- sapply(c_semead,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
c_correc <- sapply(s_correc,
FUN = apmc,
model = m0,
focus = "correc",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_correc <- sapply(c_correc,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
# Conversão de lista para tabela.
c_semead <- plyr::ldply(c_semead, .id = "correc")
c_correc <- plyr::ldply(c_correc, .id = "semead")
# Faz a junção e coloca letras altas e baixas.
resul <- c_semead %>%
mutate(cle = toupper(cld)) %>%
select(semead, correc, cle) %>%
inner_join(y = c_correc, by = c("semead", "correc")) %>%
unite(clde, cle, cld, sep = "", remove = FALSE)
resul <- equallevels(resul, db)
resul_all[[2]] <- resul
resul## semead correc clde cle fit lwr upr cld
## 1 0 N1 Bd B 2582.633 2158.137 3007.130 d
## 2 50 - lanço N1 Ca C 6710.300 6285.804 7134.796 a
## 3 50 - sulco N1 ABa AB 6731.033 6306.537 7155.530 a
## 4 100 - lanço N1 ABa AB 7896.900 7472.404 8321.396 a
## 5 100 - sulco N1 ACab AC 7342.067 6917.570 7766.563 ab
## 6 0 N2 Bb B 4441.900 4017.404 4866.396 b
## 7 50 - lanço N2 Aa A 7465.367 7040.870 7889.863 a
## 8 50 - sulco N2 ABa AB 7196.300 6771.804 7620.796 a
## 9 100 - lanço N2 ABa AB 7606.333 7181.837 8030.830 a
## 10 100 - sulco N2 Aab A 7420.433 6995.937 7844.930 ab
## 11 0 N3 Bb B 4821.600 4397.104 5246.096 b
## 12 50 - lanço N3 Aa A 7204.600 6780.104 7629.096 a
## 13 50 - sulco N3 ABa AB 7549.900 7125.404 7974.396 a
## 14 100 - lanço N3 ABa AB 7587.000 7162.504 8011.496 a
## 15 100 - sulco N3 Aab A 7636.633 7212.137 8061.130 ab
## 16 0 N4 Bc B 3371.067 2946.570 3795.563 c
## 17 50 - lanço N4 Aa A 7454.933 7030.437 7879.430 a
## 18 50 - sulco N4 ABa AB 7317.300 6892.804 7741.796 a
## 19 100 - lanço N4 ABa AB 7721.833 7297.337 8146.330 a
## 20 100 - sulco N4 Aab A 7529.000 7104.504 7953.496 ab
## 21 0 N5 Bb B 4624.033 4199.537 5048.530 b
## 22 50 - lanço N5 Aa A 7486.967 7062.470 7911.463 a
## 23 50 - sulco N5 ABa AB 7326.467 6901.970 7750.963 a
## 24 100 - lanço N5 ABa AB 7470.700 7046.204 7895.196 a
## 25 100 - sulco N5 Aab A 7673.000 7248.504 8097.496 ab
## 26 0 N6 Bcd B 3106.500 2682.004 3530.996 cd
## 27 50 - lanço N6 Ca C 7284.467 6859.970 7708.963 a
## 28 50 - sulco N6 Ba B 7011.667 6587.170 7436.163 a
## 29 100 - lanço N6 ABa AB 7974.833 7550.337 8399.330 a
## 30 100 - sulco N6 ACb AC 6927.533 6503.037 7352.030 b
## 31 0 N7 BDa BD 6866.433 6441.937 7290.930 a
## 32 50 - lanço N7 ADa AD 7235.733 6811.237 7660.230 a
## 33 50 - sulco N7 BCa BC 7438.333 7013.837 7862.830 a
## 34 100 - lanço N7 ABa AB 8050.100 7625.604 8474.596 a
## 35 100 - sulco N7 ACa AC 7948.067 7523.570 8372.563 aggplot(data = resul,
mapping = aes(x = correc, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, cld)),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~semead, nrow = 1, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Correção inicial do solo com combinações de fósforo, camada e gesso",
# "[fósforo (kg/ha) * camada (cm) * gesso (kg/ha)]"))) +
xlab("Correção inicial do solo com combinações de fósforo, camada e gesso") +
ylab("Produtividade de grãos (kg/ha)")
ggplot(data = resul,
mapping = aes(x = semead, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, tolower(cle))),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~correc, nrow = 1) +
scale_x_discrete(labels = function(l) parse(text = l)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Aplicações de" ~ P[2] * O[5] ~ "na semeadura com combinação da quantidade e forma de aplicação",
# "[" * P[2] * O[5] ~ "(kg/ha) * forma de aplicação]"))) +
xlab(expression("Doses de" ~ P[2] * O[5] ~ "na semeadura (kg/ha)")) +
ylab("Produtividade de grãos (kg/ha)")
#-----------------------------------------------------------------------3.3 Safra 2015/2016
#-----------------------------------------------------------------------
# Modelo para o delineamento de blocos casualizados. Os blocos são de
# tamanho 7 * 5 = 35.
m0 <- lmer(prod ~
(1 | bloc) +
correc * semead,
data = subset(db, safr == "2015/2016"))
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 1685988 280998 6 68 2.0322 0.07312 .
## semead 124652616 31163154 4 68 225.3750 < 2.2e-16 ***
## correc:semead 11031370 459640 24 68 3.3242 5.362e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo para o experimento em faixas casualizadas em blocos.
m1 <- update(m0,
formula = . ~
(1 | bloc) + # bloco.
(1 | bloc:correc) + # faixa horizontal dentro de bloco.
(1 | bloc:semead) + # faixa vertical dentro de bloco.
correc * semead)
anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 894634 149106 6 12 1.3931 0.2937
## semead 57860792 14465198 4 8 135.1510 2.235e-07 ***
## correc:semead 11031370 459640 24 48 4.2945 8.705e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m1)
#-----------------------------------------------------------------------
# Desdobramento da interação com comparações múltiplas.
# Para obter as médias ajustadas.
lsm <- LSmeans(m1, effect = c("correc", "semead"))
# lsm
# Cria nomes para as linhas da matriz.
rownames(lsm$L) <-
lsm$grid %>%
unite(row, correc, semead) %>%
pull(row)
# Matrizes para estudar o efeito da aplicação de fósforo (5 níveis)
# fixando o manejo do gesso (7 níveis).
s_semead <- by(data = lsm$L,
INDICES = lsm$grid[, "correc"],
FUN = as.matrix)
s_semead <- lapply(s_semead,
FUN = function(x) {
rownames(x) <- sub("^.*_", "", rownames(x))
return(x)
})
# Matrizes para estudar o efeito do manejo do gesso (7 níveis) fixando a
# aplicação de fósforo (5 níveis).
s_correc <- by(data = lsm$L,
INDICES = lsm$grid[, "semead"],
FUN = as.matrix)
s_correc <- lapply(s_correc,
FUN = function(x) {
rownames(x) <- sub("_.*$", "", rownames(x))
return(x)
})
# Faz as comparações múltiplas e acerta disposição das letras.
c_semead <- sapply(s_semead,
FUN = apmc,
model = m0,
focus = "semead",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_semead <- sapply(c_semead,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
c_correc <- sapply(s_correc,
FUN = apmc,
model = m0,
focus = "correc",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_correc <- sapply(c_correc,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
# Conversão de lista para tabela.
c_semead <- plyr::ldply(c_semead, .id = "correc")
c_correc <- plyr::ldply(c_correc, .id = "semead")
# Faz a junção e coloca letras altas e baixas.
resul <- c_semead %>%
mutate(cle = toupper(cld)) %>%
select(semead, correc, cle) %>%
inner_join(y = c_correc, by = c("semead", "correc")) %>%
unite(clde, cle, cld, sep = "", remove = FALSE)
resul <- equallevels(resul, db)
resul_all[[3]] <- resul
resul## semead correc clde cle fit lwr upr cld
## 1 0 N1 Bc B 347.3000 -91.46044 786.0604 c
## 2 50 - lanço N1 Aa A 3364.5333 2925.77289 3803.2938 a
## 3 50 - sulco N1 BCa BC 3880.3000 3441.53956 4319.0604 a
## 4 100 - lanço N1 ABa AB 4068.1000 3629.33956 4506.8604 a
## 5 100 - sulco N1 ACa AC 4012.9000 3574.13956 4451.6604 a
## 6 0 N2 Bbc B 1023.5667 584.80622 1462.3271 bc
## 7 50 - lanço N2 Aa A 3411.8667 2973.10622 3850.6271 a
## 8 50 - sulco N2 ABa AB 3600.7667 3162.00622 4039.5271 a
## 9 100 - lanço N2 ABa AB 4024.2667 3585.50622 4463.0271 a
## 10 100 - sulco N2 Aa A 3772.6333 3333.87289 4211.3938 a
## 11 0 N3 Bb B 1242.6000 803.83956 1681.3604 b
## 12 50 - lanço N3 Aa A 3821.6333 3382.87289 4260.3938 a
## 13 50 - sulco N3 ABa AB 4018.0667 3579.30622 4456.8271 a
## 14 100 - lanço N3 ABa AB 4038.8667 3600.10622 4477.6271 a
## 15 100 - sulco N3 Aa A 3830.2667 3391.50622 4269.0271 a
## 16 0 N4 Bbc B 648.6667 209.90622 1087.4271 bc
## 17 50 - lanço N4 Aa A 3679.9667 3241.20622 4118.7271 a
## 18 50 - sulco N4 ABa AB 3641.7333 3202.97289 4080.4938 a
## 19 100 - lanço N4 ABa AB 3935.8333 3497.07289 4374.5938 a
## 20 100 - sulco N4 Aa A 3762.8667 3324.10622 4201.6271 a
## 21 0 N5 Bbc B 950.6333 511.87289 1389.3938 bc
## 22 50 - lanço N5 Aa A 3574.3667 3135.60622 4013.1271 a
## 23 50 - sulco N5 ABa AB 3637.2000 3198.43956 4075.9604 a
## 24 100 - lanço N5 ABa AB 3777.6333 3338.87289 4216.3938 a
## 25 100 - sulco N5 Aa A 3813.1667 3374.40622 4251.9271 a
## 26 0 N6 Bbc B 535.1667 96.40622 973.9271 bc
## 27 50 - lanço N6 Aa A 3446.6000 3007.83956 3885.3604 a
## 28 50 - sulco N6 ABa AB 4075.4667 3636.70622 4514.2271 a
## 29 100 - lanço N6 ABa AB 3784.0667 3345.30622 4222.8271 a
## 30 100 - sulco N6 Aa A 3908.6667 3469.90622 4347.4271 a
## 31 0 N7 Ba B 2674.3000 2235.53956 3113.0604 a
## 32 50 - lanço N7 Aa A 3482.4333 3043.67289 3921.1938 a
## 33 50 - sulco N7 ABa AB 3799.7667 3361.00622 4238.5271 a
## 34 100 - lanço N7 ABa AB 3566.6333 3127.87289 4005.3938 a
## 35 100 - sulco N7 Aa A 3764.1667 3325.40622 4202.9271 aggplot(data = resul,
mapping = aes(x = correc, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, cld)),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~semead, nrow = 1, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Correção inicial do solo com combinações de fósforo, camada e gesso",
# "[fósforo (kg/ha) * camada (cm) * gesso (kg/ha)]"))) +
xlab("Correção inicial do solo com combinações de fósforo, camada e gesso") +
ylab("Produtividade de grãos (kg/ha)")
ggplot(data = resul,
mapping = aes(x = semead, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, tolower(cle))),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~correc, nrow = 1) +
scale_x_discrete(labels = function(l) parse(text = l)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Aplicações de" ~ P[2] * O[5] ~ "na semeadura com combinação da quantidade e forma de aplicação",
# "[" * P[2] * O[5] ~ "(kg/ha) * forma de aplicação]"))) +
xlab(expression("Doses de" ~ P[2] * O[5] ~ "na semeadura (kg/ha)")) +
ylab("Produtividade de grãos (kg/ha)")
#-----------------------------------------------------------------------3.4 Safra 2016/2017
#-----------------------------------------------------------------------
# Modelo para o delineamento de blocos casualizados. Os blocos são de
# tamanho 7 * 5 = 35.
m0 <- lmer(prod ~
(1 | bloc) +
correc * semead,
data = subset(db, safr == "2016/2017"))
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 21347338 3557890 6 68 12.9804 1.012e-09 ***
## semead 351890022 87972505 4 68 320.9532 < 2.2e-16 ***
## correc:semead 24228559 1009523 24 68 3.6831 1.254e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo para o experimento em faixas casualizadas em blocos.
m1 <- update(m0,
formula = . ~
(1 | bloc) + # bloco.
(1 | bloc:correc) + # faixa horizontal dentro de bloco.
(1 | bloc:semead) + # faixa vertical dentro de bloco.
correc * semead)
anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## correc 19480396 3246733 6 12 22.9705 6.413e-06 ***
## semead 39807949 9951987 4 8 70.4099 2.846e-06 ***
## correc:semead 24228559 1009523 24 48 7.1423 4.533e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m1)
#-----------------------------------------------------------------------
# Desdobramento da interação com comparações múltiplas.
# Para obter as médias ajustadas.
lsm <- LSmeans(m1, effect = c("correc", "semead"))
# lsm
# Cria nomes para as linhas da matriz.
rownames(lsm$L) <-
lsm$grid %>%
unite(row, correc, semead) %>%
pull(row)
# Matrizes para estudar o efeito da aplicação de fósforo (5 níveis)
# fixando o manejo do gesso (7 níveis).
s_semead <- by(data = lsm$L,
INDICES = lsm$grid[, "correc"],
FUN = as.matrix)
s_semead <- lapply(s_semead,
FUN = function(x) {
rownames(x) <- sub("^.*_", "", rownames(x))
return(x)
})
# Matrizes para estudar o efeito do manejo do gesso (7 níveis) fixando a
# aplicação de fósforo (5 níveis).
s_correc <- by(data = lsm$L,
INDICES = lsm$grid[, "semead"],
FUN = as.matrix)
s_correc <- lapply(s_correc,
FUN = function(x) {
rownames(x) <- sub("_.*$", "", rownames(x))
return(x)
})
# Faz as comparações múltiplas e acerta disposição das letras.
c_semead <- sapply(s_semead,
FUN = apmc,
model = m0,
focus = "semead",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_semead <- sapply(c_semead,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
c_correc <- sapply(s_correc,
FUN = apmc,
model = m0,
focus = "correc",
cld2 = TRUE,
test = "fdr",
simplify = FALSE)
c_correc <- sapply(c_correc,
simplify = FALSE,
FUN = function(x) {
x$cld <- with(x, c(ordered_cld(cld, fit)))
return(x)
})
# Conversão de lista para tabela.
c_semead <- plyr::ldply(c_semead, .id = "correc")
c_correc <- plyr::ldply(c_correc, .id = "semead")
# Faz a junção e coloca letras altas e baixas.
resul <- c_semead %>%
mutate(cle = toupper(cld)) %>%
select(semead, correc, cle) %>%
inner_join(y = c_correc, by = c("semead", "correc")) %>%
unite(clde, cle, cld, sep = "", remove = FALSE)
resul <- equallevels(resul, db)
resul_all[[4]] <- resul
resul## semead correc clde cle fit lwr upr cld
## 1 0 N1 Bd B 1667.500 1025.921 2309.079 d
## 2 50 - lanço N1 Aa A 6935.900 6294.321 7577.479 a
## 3 50 - sulco N1 ABa AB 7217.067 6575.488 7858.646 a
## 4 100 - lanço N1 ABb AB 7416.400 6774.821 8057.979 b
## 5 100 - sulco N1 Aa A 7200.167 6558.588 7841.746 a
## 6 0 N2 Bbc B 2781.900 2140.321 3423.479 bc
## 7 50 - lanço N2 Aa A 7410.600 6769.021 8052.179 a
## 8 50 - sulco N2 ABa AB 7576.767 6935.188 8218.346 a
## 9 100 - lanço N2 ABb AB 7084.833 6443.254 7726.412 b
## 10 100 - sulco N2 Aa A 7596.433 6954.854 8238.012 a
## 11 0 N3 Bb B 3491.533 2849.954 4133.112 b
## 12 50 - lanço N3 Aa A 8109.600 7468.021 8751.179 a
## 13 50 - sulco N3 ABa AB 8189.033 7547.454 8830.612 a
## 14 100 - lanço N3 ABa AB 8438.100 7796.521 9079.679 a
## 15 100 - sulco N3 Aa A 8275.733 7634.154 8917.312 a
## 16 0 N4 Bcd B 1978.933 1337.354 2620.512 cd
## 17 50 - lanço N4 Aa A 7166.033 6524.454 7807.612 a
## 18 50 - sulco N4 ABa AB 7367.033 6725.454 8008.612 a
## 19 100 - lanço N4 ABb AB 7247.900 6606.321 7889.479 b
## 20 100 - sulco N4 Aa A 7424.967 6783.388 8066.546 a
## 21 0 N5 Bb B 3636.200 2994.621 4277.779 b
## 22 50 - lanço N5 Ca C 7603.800 6962.221 8245.379 a
## 23 50 - sulco N5 ABa AB 7464.133 6822.554 8105.712 a
## 24 100 - lanço N5 ABa AB 8562.467 7920.888 9204.046 a
## 25 100 - sulco N5 ACa AC 8392.333 7750.754 9033.912 a
## 26 0 N6 Bcd B 2053.567 1411.988 2695.146 cd
## 27 50 - lanço N6 Aa A 7129.733 6488.154 7771.312 a
## 28 50 - sulco N6 ABa AB 7145.933 6504.354 7787.512 a
## 29 100 - lanço N6 ABb AB 7351.400 6709.821 7992.979 b
## 30 100 - sulco N6 Aa A 7739.433 7097.854 8381.012 a
## 31 0 N7 Ba B 5320.700 4679.121 5962.279 a
## 32 50 - lanço N7 Aa A 7532.500 6890.921 8174.079 a
## 33 50 - sulco N7 ABa AB 7778.567 7136.988 8420.146 a
## 34 100 - lanço N7 ABb AB 7223.733 6582.154 7865.312 b
## 35 100 - sulco N7 Aa A 7133.367 6491.788 7774.946 aggplot(data = resul,
mapping = aes(x = correc, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, cld)),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~semead, nrow = 1, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Correção inicial do solo com combinações de fósforo, camada e gesso",
# "[fósforo (kg/ha) * camada (cm) * gesso (kg/ha)]"))) +
xlab("Correção inicial do solo com combinações de fósforo, camada e gesso") +
ylab("Produtividade de grãos (kg/ha)")
ggplot(data = resul,
mapping = aes(x = semead, y = fit)) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, tolower(cle))),
vjust = -0.5,
size = 3,
angle = 90) +
expand_limits(x = 0) +
facet_wrap(facets = ~correc, nrow = 1) +
scale_x_discrete(labels = function(l) parse(text = l)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Aplicações de" ~ P[2] * O[5] ~ "na semeadura com combinação da quantidade e forma de aplicação",
# "[" * P[2] * O[5] ~ "(kg/ha) * forma de aplicação]"))) +
xlab(expression("Doses de" ~ P[2] * O[5] ~ "na semeadura (kg/ha)")) +
ylab("Produtividade de grãos (kg/ha)")
#-----------------------------------------------------------------------3.5 Exibição conjunta
resul <- plyr::ldply(resul_all, .id = "safr")
str(resul)## 'data.frame': 140 obs. of 9 variables:
## $ safr : Factor w/ 4 levels "2013/2014","2014/2015",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ semead: Factor w/ 5 levels "0","50 - lanço",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ correc: Factor w/ 7 levels "N1","N2","N3",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ clde : chr "Ce" "Bb" "ACc" "ACa" ...
## $ cle : chr "C" "B" "AC" "AC" ...
## $ fit : num 473 4545 5586 7427 7510 ...
## $ lwr : num -360 3712 4753 6594 6677 ...
## $ upr : num 1305 5378 6419 8260 8343 ...
## $ cld : chr "e" "b" "c" "a" ...ggplot(data = resul,
mapping = aes(x = correc, y = fit)) +
geom_point(size = 1) +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, cld)),
vjust = -0.5,
size = 2.5,
angle = 90) +
expand_limits(x = 0) +
facet_grid(facets = safr ~ semead, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Correção inicial do solo com combinações de fósforo, camada e gesso",
# "[fósforo (kg/ha) * camada (cm) * gesso (kg/ha)]"))) +
xlab("Correção inicial do solo com combinações de fósforo, camada e gesso") +
ylab("Produtividade de grãos (kg/ha)")
# write_tsv(resul, path = "resultados-medias.txt")
ggplot(data = resul,
mapping = aes(x = semead, y = fit)) +
geom_point(size = 1) +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.0f %s", fit, tolower(cle))),
vjust = -0.5,
size = 2.5,
angle = 90) +
expand_limits(x = 0) +
facet_grid(facets = safr ~ correc) +
scale_x_discrete(labels = function(l) parse(text = l)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
strip.text.x = element_text(size = 8)) +
# xlab(expression(atop(
# "Aplicações de" ~ P[2] * O[5] ~ "na semeadura com combinação da quantidade e forma de aplicação",
# "[" * P[2] * O[5] ~ "(kg/ha) * forma de aplicação]"))) +
xlab(expression("Doses de" ~ P[2] * O[5] ~ "na semeadura (kg/ha)")) +
ylab("Produtividade de grãos (kg/ha)")
ggplot(data = resul,
mapping = aes(x = safr, y = fit)) +
geom_point(size = 1) +
geom_line(mapping = aes(group = 1), color = "gray50") +
geom_errorbar(mapping = aes(ymin = lwr, ymax = upr),
width = 0) +
facet_grid(facets = correc ~ semead, labeller = label_parsed) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
strip.text.x = element_text(size = 8)) +
xlab("Safra") +
ylab("Produtividade de grãos (kg/ha)")
#-----------------------------------------------------------------------