---
title: >
  Doses de Nitrogênio e Potássio para Cana-de-açúcar em Diferentes
  Sistemas de Manejo no Estado do Paraná
author: >
  [Walmes Zeviani](http://www.leg.ufpr.br/doku.php/pessoais:walmes) &
  [Michael Jonathan Fernandes Alves](http://lattes.cnpq.br/0834978726176394)
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{Doses de Nitrogênio e Potássio para Cana-de-açúcar em Diferentes Sistemas de Manejo no Estado do Paraná}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Session Definition

```{r, message=FALSE, results="hide"}
# https://github.com/walmes/wzRfun
# devtools::install_github("walmes/wzRfun")
library(wzRfun)

pks <- c("lattice",
         "gridExtra",
         "lme4",
         "lmerTest",
         "doBy",
         "multcomp")
sapply(pks, library, character.only = TRUE, logical.return = TRUE)
```
```{r, eval=FALSE}
library(wzCoop)
```
```{r setup, include=FALSE}
source("_setup.R")
```

## Exploratory data analysis

```{r}
# Object structure.
str(sugarcane_straw)

# Frequencies.
ftable(xtabs(~palha + N + K, data = sugarcane_straw))

# Checking if is a complete cases dataset.
all(complete.cases(sugarcane_straw))

# Descriptive measures.
summary(sugarcane_straw)

# A more detailed description.
# Hmisc::describe(sugarcane_straw)

# Create factors.
sugarcane_straw <- within(sugarcane_straw, {
    # Convert integer to factor.
    palha <- factor(sugarcane_straw$palha)
    nit <- factor(sugarcane_straw$N)
    pot <- factor(sugarcane_straw$K)
    # Create a factor to represent main plots.
    ue <- interaction(bloc, palha)
})
```

****
## Número de colmos por metro

```{r}
#-----------------------------------------------------------------------
# To minimize modifications of code when analysing diferent responses.

# Define the legends for factors and variables.
leg <- list(N = expression("Nitrogênio"~(kg~ha^{-1})),
            K = expression("Potássio"~(kg~ha^{-1})),
            y = expression("Número de colmos por metro"),
            palha = c("Cobertura do solo", "Com palha", "Sem palha"))

# Use name y for the response.
sugarcane_straw$y <- sugarcane_straw$ncm

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
# plot(pk)
plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)

# Drop non relevant terms.
# m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue))
m1 <- lmer(y ~ palha * nit + pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m1, m0)

anova(m1)
# summary(m1)

# Reducing by using linear effect for N and K.
m2 <- lmer(y ~ palha * N + K + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m2, m0)

anova(m2)
summary(m2)

# Final model.
mod <- m2

#-----------------------------------------------------------------------
# Fitted values.

# Experimental values to get estimates.
grid <- with(sugarcane_straw,
             expand.grid(palha = levels(palha),
                         N = seq(min(N), max(N), length.out = 10),
                         K = seq(min(N), max(N), length.out = 3),
                         KEEP.OUT.ATTRS = FALSE))

# Matrix of fixed effects.
X <- model.matrix(terms(mod), data = cbind(grid, y = 0))

# Confidence intervals.
ci <- confint(glht(mod, linfct = X),
              calpha = univariate_calpha())$confint
colnames(ci)[1] <- "fit"

grid <- cbind(grid, ci)
str(grid)

# Sample averages.
# averages <- aggregate(nobars(formula(mod)),
#                       data = sugarcane_straw, FUN = mean)
#
# xyplot(y ~ N | factor(K), groups = palha, data = averages,
#        type = "o", as.table = TRUE)

xyplot(fit ~ N | factor(K), groups = palha, data = grid,
       type = "l", as.table = TRUE, layout = c(NA, 1),
       xlab = leg$N, ylab = leg$y,
       auto.key = list(title = leg$palha[1], cex.title = 1.1,
                       text = c(leg$palha[-1]), columns = 2,
                       points = FALSE, lines = TRUE),
       uy = grid$upr, ly = grid$lwr,
       cty = "bands", alpha = 0.25, fill = "gray50",
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose)

```

****
## Peso médio de colmo

```{r}
#-----------------------------------------------------------------------

leg$y <- "Peso médio de colmo (kg)"
sugarcane_straw$y <- sugarcane_straw$pmc

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
# plot(pk)
plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)

# Drop non relevant terms.
# m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue))
m1 <- lmer(y ~ palha + nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m1, m0)

anova(m1)
# summary(m1)

# Reducing by using linear effect for N and K.
m2 <- lmer(y ~ palha + N * K + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m2, m0)

anova(m2)
summary(m2)

# Final model.
mod <- m2

#-----------------------------------------------------------------------
# Fitted values.

# Experimental values to get estimates.
grid <- with(sugarcane_straw,
             expand.grid(palha = levels(palha),
                         N = seq(min(N), max(N), length.out = 10),
                         K = seq(min(N), max(N), length.out = 3),
                         KEEP.OUT.ATTRS = FALSE))

# Matrix of fixed effects.
X <- model.matrix(terms(mod), data = cbind(grid, y = 0))

# Confidence intervals.
ci <- confint(glht(mod, linfct = X),
              calpha = univariate_calpha())$confint
colnames(ci)[1] <- "fit"

grid <- cbind(grid, ci)
str(grid)

# Sample averages.
# averages <- aggregate(nobars(formula(mod)),
#                       data = sugarcane_straw, FUN = mean)
#
# xyplot(y ~ N | factor(K), groups = palha, data = averages,
#        type = "o", as.table = TRUE)

xyplot(fit ~ N | factor(K), groups = palha, data = grid,
       type = "l", as.table = TRUE, layout = c(NA, 1),
       xlab = leg$N, ylab = leg$y,
       auto.key = list(title = leg$palha[1], cex.title = 1.1,
                       text = c(leg$palha[-1]), columns = 2,
                       points = FALSE, lines = TRUE),
       uy = grid$upr, ly = grid$lwr,
       cty = "bands", alpha = 0.25, fill = "gray50",
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose)

```

****
## Produção de cana-de-açúcar

```{r}
#-----------------------------------------------------------------------

leg$y <- expression("Produção de cana"~(ton~ha^{-1}))
sugarcane_straw$y <- sugarcane_straw$tch

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
# plot(pk)
plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)

# Drop non relevant terms.
# m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue))
m1 <- lmer(y ~ palha + nit + pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m1, m0)

anova(m1)
# summary(m1)

# Reducing by using linear effect for N and K.
m2 <- lmer(y ~ palha + N + K + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m2, m0)

anova(m2)
summary(m2)

# Final model.
mod <- m2

#-----------------------------------------------------------------------
# Fitted values.

# Experimental values to get estimates.
grid <- with(sugarcane_straw,
             expand.grid(palha = levels(palha),
                         N = seq(min(N), max(N), length.out = 10),
                         K = seq(min(N), max(N), length.out = 3),
                         KEEP.OUT.ATTRS = FALSE))

# Matrix of fixed effects.
X <- model.matrix(terms(mod), data = cbind(grid, y = 0))

# Confidence intervals.
ci <- confint(glht(mod, linfct = X),
              calpha = univariate_calpha())$confint
colnames(ci)[1] <- "fit"

grid <- cbind(grid, ci)
str(grid)

# Sample averages.
# averages <- aggregate(nobars(formula(mod)),
#                       data = sugarcane_straw, FUN = mean)
#
# xyplot(y ~ N | factor(K), groups = palha, data = averages,
#        type = "o", as.table = TRUE)

xyplot(fit ~ N | factor(K), groups = palha, data = grid,
       type = "l", as.table = TRUE, layout = c(NA, 1),
       xlab = leg$N, ylab = leg$y,
       auto.key = list(title = leg$palha[1], cex.title = 1.1,
                       text = c(leg$palha[-1]), columns = 2,
                       points = FALSE, lines = TRUE),
       uy = grid$upr, ly = grid$lwr,
       cty = "bands", alpha = 0.25, fill = "gray50",
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose)
```

****
## Teor de sacarose aparente

```{r}
#-----------------------------------------------------------------------

leg$y <- "Teor de sacarose aparente"
sugarcane_straw$y <- sugarcane_straw$pol

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
# plot(pk)
plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)
```

****
## Produção de sacarose

```{r}
#-----------------------------------------------------------------------

leg$y <- expression("Sacarose"~(ton~ha^{-1}))
sugarcane_straw$y <- sugarcane_straw$tsh

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
# plot(pk)
plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)

# Drop non relevant terms.
# m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue))
m1 <- lmer(y ~ palha + nit + pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m1, m0)

anova(m1)
# summary(m1)

# Reducing by using linear effect for N and K.
m2 <- lmer(y ~ palha + N + K + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Test the reduced model.
anova(m2, m0)

anova(m2)
summary(m2)

# Final model.
mod <- m2

#-----------------------------------------------------------------------
# Fitted values.

# Experimental values to get estimates.
grid <- with(sugarcane_straw,
             expand.grid(palha = levels(palha),
                         N = seq(min(N), max(N), length.out = 10),
                         K = seq(min(N), max(N), length.out = 3),
                         KEEP.OUT.ATTRS = FALSE))

# Matrix of fixed effects.
X <- model.matrix(terms(mod), data = cbind(grid, y = 0))

# Confidence intervals.
ci <- confint(glht(mod, linfct = X),
              calpha = univariate_calpha())$confint
colnames(ci)[1] <- "fit"

grid <- cbind(grid, ci)
str(grid)

# Sample averages.
# averages <- aggregate(nobars(formula(mod)),
#                       data = sugarcane_straw, FUN = mean)
#
# xyplot(y ~ N | factor(K), groups = palha, data = averages,
#        type = "o", as.table = TRUE)

xyplot(fit ~ N | factor(K), groups = palha, data = grid,
       type = "l", as.table = TRUE, layout = c(NA, 1),
       xlab = leg$N, ylab = leg$y,
       auto.key = list(title = leg$palha[1], cex.title = 1.1,
                       text = c(leg$palha[-1]), columns = 2,
                       points = FALSE, lines = TRUE),
       uy = grid$upr, ly = grid$lwr,
       cty = "bands", alpha = 0.25, fill = "gray50",
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose)
```

****
## Teor de nitrogênio nas folhas

```{r}
#-----------------------------------------------------------------------

leg$y <- "Teor de nitrogênio nas folhas"
sugarcane_straw$y <- sugarcane_straw$tfn

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
# plot(pk)
plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)

# Final model.
mod <- m0

#-----------------------------------------------------------------------
# Fitted values.

# Experimental values to get estimates.
grid <- with(sugarcane_straw,
             expand.grid(palha = levels(palha),
                         # N = seq(min(N), max(N), length.out = 10),
                         # K = seq(min(N), max(N), length.out = 3),
                         nit = levels(nit),
                         pot = levels(pot),
                         KEEP.OUT.ATTRS = FALSE))

# Matrix of fixed effects.
X <- model.matrix(terms(mod), data = cbind(grid, y = 0))

# Confidence intervals.
ci <- confint(glht(mod, linfct = X),
              calpha = univariate_calpha())$confint
colnames(ci)[1] <- "fit"

grid <- cbind(grid, ci)
str(grid)

# Sample averages.
# averages <- aggregate(nobars(formula(mod)),
#                       data = sugarcane_straw, FUN = mean)
#
# xyplot(y ~ N | factor(K), groups = palha, data = averages,
#        type = "o", as.table = TRUE)

xyplot(fit ~ pot | nit, groups = palha, data = grid,
       as.table = TRUE, layout = c(NA, 1), type = "o",
       xlab = leg$N, ylab = leg$y,
       auto.key = list(title = leg$palha[1], cex.title = 1.1,
                       text = c(leg$palha[-1]), columns = 2,
                       points = FALSE, lines = TRUE),
       uy = grid$upr, ly = grid$lwr,
       desloc =  0.2 * scale(as.integer(grid$palha), scale = FALSE),
       cty = "bars", length = 0.02,
       alpha = 0.25, fill = "gray50",
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose)
```

High order interactions present but no smooth effect of numeric factors.

****
## Teor de potássio nas folhas

```{r}
#-----------------------------------------------------------------------

leg$y <- "Teor de potássio nas folhas"
sugarcane_straw$y <- sugarcane_straw$tfk

#-----------------------------------------------------------------------
# Scatter plots.

pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$K, ylab = leg$y,
             auto.key = list(title = leg$N,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw,
             type = c("p", "a"),
             xlab = leg$N, ylab = leg$y,
             auto.key = list(title = leg$K,
                             cex.title = 1.1, columns = 5),
             strip = strip.custom(
                 strip.names = FALSE, strip.levels = TRUE,
                 factor.levels = c("Com palha", "Sem palha")))

# grid.arrange(pn, pk)
plot(pk)
# plot(pn)
```

```{r}
#-----------------------------------------------------------------------
# Model fitting.

# Saturated model.
m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue),
           data = sugarcane_straw, REML = FALSE)

# Simple diagnostic.
plot(m0)
# qqnorm(residuals(m0, type = "pearson"))

# Estimates of the variance components.
VarCorr(m0)

# Tests for the fixed effects.
anova(m0)

# Estimates of the effects and fitting measures.
# summary(m0)

# Final model.
mod <- m0

#-----------------------------------------------------------------------
# Fitted values.

# Experimental values to get estimates.
grid <- with(sugarcane_straw,
             expand.grid(palha = levels(palha),
                         # N = seq(min(N), max(N), length.out = 10),
                         # K = seq(min(N), max(N), length.out = 3),
                         nit = levels(nit),
                         pot = levels(pot),
                         KEEP.OUT.ATTRS = FALSE))

# Matrix of fixed effects.
X <- model.matrix(terms(mod), data = cbind(grid, y = 0))

# Confidence intervals.
ci <- confint(glht(mod, linfct = X),
              calpha = univariate_calpha())$confint
colnames(ci)[1] <- "fit"

grid <- cbind(grid, ci)
str(grid)

# Sample averages.
# averages <- aggregate(nobars(formula(mod)),
#                       data = sugarcane_straw, FUN = mean)
#
# xyplot(y ~ N | factor(K), groups = palha, data = averages,
#        type = "o", as.table = TRUE)

xyplot(fit ~ pot | nit, groups = palha, data = grid,
       as.table = TRUE, layout = c(NA, 1), type = "o",
       xlab = leg$N, ylab = leg$y,
       auto.key = list(title = leg$palha[1], cex.title = 1.1,
                       text = c(leg$palha[-1]), columns = 2,
                       points = FALSE, lines = TRUE),
       uy = grid$upr, ly = grid$lwr,
       desloc =  0.2 * scale(as.integer(grid$palha), scale = FALSE),
       cty = "bars", length = 0.02,
       alpha = 0.25, fill = "gray50",
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose)
```

****
## Session information

```{r, echo=FALSE, results="hold"}
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"),
    "----------------------------------------", sep = "\n")
sessionInfo()
```
