Material complementar ao artigo com título ??? submetido à revista ???, volume ???, número ???. A citação será incluída assim que for determinado o volume e paginação.

%% Citação LaTex.
@article{SerafimAbacaxi2014,
    author    = {},
    title     = {Curva de retenção de Água na Cultura do Abacaxizeiro em
                 Função da Cobertura do Solo e Aplicação de Gesso},
    doi       = {},
    journal   = {},
    month     = {},
    pages     = {},
    publisher = {},
    url       = {},
    year      = {2014}
}

Donwnload:

Definições da sessão

#-----------------------------------------------------------------------
# Definições da sessão, pacotes a serem usados.

# Instruções para instalação do wzRfun em:
# browseURL("https://github.com/walmes/wzRfun")

pkg <- c("lattice",
         "latticeExtra",
         "nlme",
         "reshape",
         "plyr",
         "doBy",
         "multcomp",
         "wzRfun")
sapply(pkg, library, character.only = TRUE, logical.return = TRUE)
##      lattice latticeExtra         nlme      reshape         plyr         doBy 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##     multcomp       wzRfun 
##         TRUE         TRUE
#-----------------------------------------------------------------------
# Configurações de cores da lattice.

colp <- brewer.pal(11, "Spectral")
colp <- colorRampPalette(colp, space = "rgb")
mycol <- brewer.pal(8, "Dark2")
ps <- list(
    box.dot = list(col = "black", pch = "|"),
    box.rectangle = list(col = "black", fill = c("gray70")),
    box.umbrella = list(col = "black", lty = 1),
    dot.symbol = list(col = "black"),
    dot.line = list(col = "gray90", lty = 1),
    plot.symbol = list(col = "black", cex = 0.8),
    plot.line = list(col = "black", lty = 1.2),
    plot.polygon = list(col = "gray80"),
    superpose.line = list(col = mycol),
    superpose.symbol = list(col = mycol),
    superpose.polygon = list(col = mycol),
    regions = list(col = colp(100)),
    strip.background = list(col = c("gray90", "gray50")),
    strip.shingle = list(col = c("gray50", "gray90")))
trellis.par.set(ps)
dev.off()
## null device 
##           1
rm(list = c("colp", "mycol", "ps"))

Importação dos dados

#-----------------------------------------------------------------------
# Lendo arquivos de dados.

# Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/abacaxicra.txt"

path <- ifelse(Sys.info()["user"] == "walmes",
               yes = basename(url), no = url)

# Importa os dados a partir do arquivo na web.
da <- read.table(file = path, header = TRUE, sep = "\t")
str(da)
## 'data.frame':    1664 obs. of  7 variables:
##  $ cober: Factor w/ 2 levels "Com milheto",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gesso: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ prof : Factor w/ 2 levels "0-0.05","0.05-0.20": 1 1 1 1 1 1 1 1 1 1 ...
##  $ varie: Factor w/ 4 levels "IAC Fantástico",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ rept : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ tens : int  0 0 0 0 1 1 1 1 2 2 ...
##  $ umid : num  0.488 0.433 0.343 0.415 0.489 ...
da$cober <- factor(da$cober, labels = c("Millet", "No millet"))
levels(da$varie)[3] <- "Pérola"

# write.table(da,
#             file = "~/repos/wzCoop/data-raw/pineapple_swrc.csv",
#             row.names = FALSE, quote = FALSE, sep = ";")

# Converte variáveis para fator.
da <- transform(da,
                gesso = factor(gesso),
                ue = interaction(cober, gesso, prof, varie, rept,
                                 drop = TRUE, sep = ":"))

# Passa a tensão 0 para 0.01 para não ter problemas com o log().
da$tens[da$tens == 0] <- 0.01
da$ltens <- log10(da$tens)

# Mostra a estrutura do objeto.
str(da)
## 'data.frame':    1664 obs. of  9 variables:
##  $ cober: Factor w/ 2 levels "Millet","No millet": 1 1 1 1 1 1 1 1 1 1 ...
##  $ gesso: Factor w/ 2 levels "0","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ prof : Factor w/ 2 levels "0-0.05","0.05-0.20": 1 1 1 1 1 1 1 1 1 1 ...
##  $ varie: Factor w/ 4 levels "IAC Fantástico",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ rept : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ tens : num  0.01 0.01 0.01 0.01 1 1 1 1 2 2 ...
##  $ umid : num  0.488 0.433 0.343 0.415 0.489 ...
##  $ ue   : Factor w/ 128 levels "Millet:0:0-0.05:IAC Fantástico:1",..: 17 49 81 113 17 49 81 113 17 49 ...
##  $ ltens: num  -2 -2 -2 -2 0 ...

Análise exploratória

#-----------------------------------------------------------------------
# Análise exploratória.

# Remover valores fora do (0,1) na umidade (erro de digitação).
da <- droplevels(subset(da, umid < 1 & umid > 0))

# Identifica as unidades com a combinação tripla dos fatores.
da$gescobpro <- with(da, interaction(as.integer(gesso),
                                     as.integer(cober),
                                     as.integer(prof)))

p1 <- xyplot(umid ~ tens | gescobpro + varie,
             groups = rept,
             data = da,
             type = "b",
             ylab = "Moisture",
             xlab = "log 10 of matric tension",
             scales = list(x = list(log = 10)),
             xscale.components = xscale.components.log10ticks)
useOuterStrips(p1)

# Número de curvas.
nlevels(da$ue)
## [1] 128
#-----------------------------------------------------------------------
# Identificar os pontos atípicos.

# i <- 1:nrow(da)
# useOuterStrips(p1)
# panel.locs <- trellis.currentLayout()
# outliers <- vector(mode = "list",
#                    length = prod(dim(panel.locs)))
# for (row in 1:nrow(panel.locs)) {
#     for (column in 1:ncol(panel.locs)) {
#         if (panel.locs[row, column] > 0) {
#             trellis.focus("panel",
#                           row = row,
#                           column = column,
#                           highlight = TRUE)
#             wp <- panel.locs[row, column]
#             obsv <- panel.identify()
#             # print(obsv)
#             j <- column + (row - 1) * nrow(panel.locs)
#             # print(j)
#             outliers[[j]] <- i[obsv]
#             trellis.unfocus()
#         }
#     }
# }

# Clicar com o botão direito em cada cela (de baixo para cima, esquerda
# para direita) para identificar o ponto outliear. Clicar com o botão
# direito para passara para a próxima cela.

# outliers
# r <- unlist(outliers)
# dput(r)

# Outiliers identificados visualmente.
out <- c(106L, 151L, 1397L, 513L, 717L, 725L)
da <- da[-out, ]

Ajuste de forma interativa

#-----------------------------------------------------------------------
# Ajuste com rp.nls.

model <- umid ~ Ur + (Us - Ur)/(1 + exp(n * (alp + ltens)))^(1 - 1/n)
start <- list(Ur = c(0.1,   0, 0.4),
              Us = c(0.4, 0.1, 0.8),
              alp =c(1,    -5,   6),
              n =  c(1.5,   1,   4))

library(rpanel)

cra.fit <- rp.nls(model = model,
                  data = da,
                  start = start)

coef(cra.fit)
sapply(cra.fit, coef)

Ajuste por unidade experimental

O ajuste foi feito considerando a seguinte parametrização do modelo van Genuchten

\[ U(x) = U_r + \frac{U_s - U_r}{(1 + \exp\{n(\alpha + x)\})^{1 - 1/n}} \]

em que \(U\) é umidade (m3 m-3) do solo, \(x\) é o log na base 10 da tensão matricial aplicada (kPa), \(U_r\) é a umidade residual (assíntota inferior), \(U_s\) é a umidade de satuação (assíntota superior), \(\alpha\) e \(n\) são parâmetros empíricos de forma da curva de retenção de água. Uma vez conhecido valores para as quatidades mencionadas, são obtidos

\[ \begin{align*} S &= -n\cdot \frac{U_s - U_r}{(1 + 1/m)^{m + 1}} \newline I &= -\alpha - \log(m)/n \newline U_I &= U(x = I) \end{align*} \]

em que \(S\) é a taxa no ponto de inflexão, parâmetro que é tido como central para avaliação da qualidade física do solo, bem como \(I\) que corresponde ao log da tensão no ponto de inflexão da curva de retenção de água do solo. A umidade correspondente a tensão no ponto de inflexão é representada por \(U_I\).

#-----------------------------------------------------------------------
# Ajustar uma única curva.

model <- umid ~ Ur + (Us - Ur)/(1 + exp(n * (alp + ltens)))^(1 - 1/n)
start <- list(Ur = 0.1, Us = 0.4, alp = -0.5, n = 4)

plot(umid ~ ltens, data = da)
with(start, curve(Ur + (Us - Ur)/(1 + exp(n * (alp + x)))^(1 - 1/n),
                  add = TRUE))

n00 <- nls(model, data = da, start = start)
coef(n00)
##          Ur          Us         alp           n 
##  0.09186589  0.36666033 -0.68298513  4.05581675
#-----------------------------------------------------------------------
# Ajudar para cada unidade experimental.

# nlevels(da$ue)
db <- groupedData(umid ~ ltens | ue, data = da, order.groups = FALSE)
str(db)
## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   1656 obs. of  10 variables:
##  $ cober    : Factor w/ 2 levels "Millet","No millet": 1 1 1 1 1 1 1 1 1 1 ...
##  $ gesso    : Factor w/ 2 levels "0","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ prof     : Factor w/ 2 levels "0-0.05","0.05-0.20": 1 1 1 1 1 1 1 1 1 1 ...
##  $ varie    : Factor w/ 4 levels "IAC Fantástico",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ rept     : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ tens     : num  0.01 0.01 0.01 0.01 1 1 1 1 2 2 ...
##  $ umid     : num  0.488 0.433 0.343 0.415 0.489 ...
##  $ ue       : Factor w/ 128 levels "Millet:0:0-0.05:IAC Fantástico:1",..: 17 49 81 113 17 49 81 113 17 49 ...
##  $ ltens    : num  -2 -2 -2 -2 0 ...
##  $ gescobpro: Factor w/ 8 levels "1.1.1","2.1.1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "formula")=Class 'formula'  language umid ~ ltens | ue
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  - attr(*, "FUN")=function (x)  
##  - attr(*, "order.groups")= logi FALSE
n0 <- nlsList(model = model,
              data = db,
              start = as.list(coef(n00)))
c0 <- coef(n0)

# Matriz de diagramas de dispersão.
pairs(c0)

#-----------------------------------------------------------------------
# Correr análises de variância considerando parâmetros como variáveis
# resposta.

aux <- do.call(rbind, strsplit(rownames(c0), split = ":"))
colnames(aux) <- c("cober", "gesso", "prof", "varie", "rept")
params <- cbind(equallevels(as.data.frame(aux), da), c0)
rownames(params) <- NULL
str(params)
## 'data.frame':    128 obs. of  9 variables:
##  $ cober: Factor w/ 2 levels "Millet","No millet": 1 2 1 2 1 2 1 2 1 2 ...
##  $ gesso: Factor w/ 2 levels "0","4": 1 1 2 2 1 1 2 2 1 1 ...
##  $ prof : Factor w/ 2 levels "0-0.05","0.05-0.20": 1 1 1 1 2 2 2 2 1 1 ...
##  $ varie: Factor w/ 4 levels "IAC Fantástico",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ rept : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Ur   : num  0.104 0.068 0.094 0.103 0.115 ...
##  $ Us   : num  0.412 0.342 0.256 0.382 0.381 ...
##  $ alp  : num  -0.709 -0.66 -0.906 -0.648 -0.727 ...
##  $ n    : num  5.17 4.33 7.33 4.29 5.61 ...
#-----------------------------------------------------------------------
# Manova.

m0 <- aov(cbind(Ur, Us, n, alp) ~ cober * gesso * prof * varie,
          data = params)
anova(m0)
## Analysis of Variance Table
## 
##                        Df  Pillai approx F num Df den Df    Pr(>F)    
## (Intercept)             1 0.99725   8436.8      4     93 < 2.2e-16 ***
## cober                   1 0.06270      1.6      4     93  0.192794    
## gesso                   1 0.01480      0.3      4     93  0.844016    
## prof                    1 0.30140     10.0      4     93 8.575e-07 ***
## varie                   3 0.07473      0.6     12    285  0.836124    
## cober:gesso             1 0.17813      5.0      4     93  0.001014 ** 
## cober:prof              1 0.02961      0.7      4     93  0.587557    
## gesso:prof              1 0.01470      0.3      4     93  0.845534    
## cober:varie             3 0.14064      1.2     12    285  0.305700    
## gesso:varie             3 0.07338      0.6     12    285  0.845461    
## prof:varie              3 0.17447      1.5     12    285  0.136097    
## cober:gesso:prof        1 0.01719      0.4      4     93  0.803348    
## cober:gesso:varie       3 0.07002      0.6     12    285  0.867475    
## cober:prof:varie        3 0.09755      0.8     12    285  0.652301    
## gesso:prof:varie        3 0.10261      0.8     12    285  0.607928    
## cober:gesso:prof:varie  3 0.07431      0.6     12    285  0.839035    
## Residuals              96                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Curva média de acordo com a manova.

# Será obtida a média dos parâmentros de acordo com termos que foram
# signficativos na manova. No entanto, deve ser lembrado que a curva
# obtida por meio de média amostral de estimativas não corresponde a
# curva "média" ajustada a partir dos dados.

pm <- aggregate(cbind(Us = Us, Ur = Ur, alp = alp, n = n) ~
                    cober + gesso + prof,
                data = params, FUN = mean)

range(da$ltens)
## [1] -2.000000  3.176091
pred <- with(da, data.frame(ltens = seq(min(ltens), max(ltens),
                                        length.out = 30)))

preds <- ddply(pm,
               .variables = .(cober, gesso, prof),
               .fun = summarise,
               umid = Ur + (Us - Ur)/(1 + exp(n * (alp + pred$ltens)))^(1 - 1/n),
               ltens = pred$ltens)
str(preds)
## 'data.frame':    240 obs. of  5 variables:
##  $ cober: Factor w/ 2 levels "Millet","No millet": 1 1 1 1 1 1 1 1 1 1 ...
##  $ gesso: Factor w/ 2 levels "0","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ prof : Factor w/ 2 levels "0-0.05","0.05-0.20": 1 1 1 1 1 1 1 1 1 1 ...
##  $ umid : num  0.389 0.389 0.389 0.389 0.389 ...
##  $ ltens: num  -2 -1.82 -1.64 -1.46 -1.29 ...
p1 <- xyplot(umid ~ ltens | prof + cober, groups = gesso,
             data = preds, type = "l",
             xlab = expression("log of tension"~(kPa)),
             ylab = expression(Moisture~(m^3~m^{-3})),
             auto.key = list(columns = 2, lines = TRUE, points = FALSE,
                             title = expression("Gypsum"~(ton~ha^{-1})),
                             cex.title = 1.1))
useOuterStrips(p1)

p2 <- xyplot(umid ~ ltens | prof + gesso, groups = cober,
             data = preds, type = "l",
             xlab = expression("log of tension"~(kPa)),
             ylab = expression("Moisture"~(m^3~m^{-3})),
             auto.key = list(columns = 2, lines = TRUE, points = FALSE,
                             title = "Crop cover", cex.title = 1.1))
useOuterStrips(p2)

p3 <- xyplot(umid ~ ltens | cober + gesso, groups = prof,
             data = preds, type = "l",
             xlab = expression("log of tension"~(kPa)),
             ylab = expression("Moisture"~(m^3~m^{-3})),
             auto.key = list(columns = 2, lines = TRUE, points = FALSE,
                             title = expression("Soil layer"~(m)),
                             cex.title = 1.1))
useOuterStrips(p3)

Ur: umidade residual

#-----------------------------------------------------------------------
# Análise do Ur.

m0 <- lm(Ur ~ (cober + gesso + prof + varie)^2, data = params)
par(mfrow = c(2,2)); plot(m0); layout(1)

# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[,"dffit"]

m0 <- update(m0, data = params[!del,])
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: Ur
##              Df    Sum Sq   Mean Sq F value    Pr(>F)    
## cober         1 0.0000014 0.0000014  0.0065   0.93587    
## gesso         1 0.0000746 0.0000746  0.3481   0.55642    
## prof          1 0.0014229 0.0014229  6.6376   0.01134 *  
## varie         3 0.0007300 0.0002433  1.1351   0.33826    
## cober:gesso   1 0.0044550 0.0044550 20.7815 1.363e-05 ***
## cober:prof    1 0.0000297 0.0000297  0.1384   0.71056    
## cober:varie   3 0.0015884 0.0005295  2.4699   0.06580 .  
## gesso:prof    1 0.0002396 0.0002396  1.1176   0.29280    
## gesso:varie   3 0.0002934 0.0000978  0.4563   0.71340    
## prof:varie    3 0.0001364 0.0000455  0.2121   0.88785    
## Residuals   108 0.0231522 0.0002144                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Desdobramento.

# Estimativas para as profundidades.
Xm <- LSmatrix(m0, effect = "prof")
rownames(Xm) <- levels(db$prof)
g1 <- apmc(X = Xm, model = m0, focus = "prof")
g1
##        prof        fit        lwr        upr cld
## 1    0-0.05 0.09509044 0.09146270 0.09871818   a
## 2 0.05-0.20 0.08861510 0.08495424 0.09227597   b
# Desdobrar gesso dentro de cobertura.
LSmeans(m0, effect = c("gesso", "cober"))
##     estimate          se  df   t.stat      p.value gesso     cober
## 1 0.09713565 0.002635331 108 36.85899 5.306475e-63     0    Millet
## 2 0.08651437 0.002588268 108 33.42559 8.795253e-59     4    Millet
## 3 0.08520292 0.002588268 108 32.91889 3.955880e-58     0 No millet
## 4 0.09855814 0.002588268 108 38.07880 2.022286e-64     4 No millet
Xm <- LSmatrix(m0, effect = c("gesso", "cober"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES = grid$cober, FUN = as.matrix)
L <- lapply(L, "rownames<-", levels(grid$gesso))
g2 <- lapply(L, apmc, model = m0, focus = "gesso")
g2
## $Millet
##   gesso        fit        lwr        upr cld
## 1     0 0.09713565 0.09191197 0.10235934   b
## 2     4 0.08651437 0.08138398 0.09164477   a
## 
## $`No millet`
##   gesso        fit        lwr        upr cld
## 1     0 0.08520292 0.08007252 0.09033331   b
## 2     4 0.09855814 0.09342775 0.10368854   a
g2 <- ldply(g2)
names(g2)[1] <- "cober"
g2 <- equallevels(g2, db)
#--------------------------------------------
# Gráfico.

p1 <- segplot(prof ~ lwr + upr, centers = fit, data = g1, draw = FALSE,
              xlab = expression("Residual moisture"~(m^3~m^{-3})),
              ylab = "Soil layer (m)") +
    layer(panel.text(x = centers, y = z, labels = g1$cld, pos = 3))

p2 <- segplot(cober ~ lwr + upr, centers = fit, data = g2, draw = FALSE,
              groups = g2$gesso, pch = g2$gesso, gap = 0.1,
              panel = panel.groups.segplot,
              xlab = expression("Residual moisture"~(m^3~m^{-3})),
              ylab = "Cover crop",
              key = list(columns = 2, type = "o", divide = 1,
                         title = expression("Gypsum"~(ton~ha^{-1})),
                         cex.title = 1.1,
                         lines = list(pch = 1:2),
                         text = list(levels(g2$gesso)))) +
    layer(panel.text(x = centers,
                     y = as.numeric(z) +
                         gap * (2 * ((as.numeric(groups) - 1)/
                                     (nlevels(groups) - 1)) - 1),
                     labels = cld, pos = 3), data = g2)
p2

# x11(width = 6, height = 5)
d <- 0.6
plot(p1, position = c(0, d, 1, 1), more = TRUE)
plot(p2, position = c(0, 0, 1, d), more = FALSE)

Us: umidade de saturação

#-----------------------------------------------------------------------
# Análise do Us.

m0 <- lm(Us ~ (cober + gesso + prof + varie)^2, data  =  params)
par(mfrow = c(2,2)); plot(m0); layout(1)

# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[, "dffit"]

m0 <- update(m0, data = params[!del,])
par(mfrow = c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: Us
##              Df   Sum Sq  Mean Sq F value    Pr(>F)    
## cober         1 0.004879 0.004879  1.8945   0.17154    
## gesso         1 0.001788 0.001788  0.6943   0.40656    
## prof          1 0.075827 0.075827 29.4400 3.566e-07 ***
## varie         3 0.001157 0.000386  0.1497   0.92967    
## cober:gesso   1 0.006721 0.006721  2.6095   0.10914    
## cober:prof    1 0.003100 0.003100  1.2037   0.27503    
## cober:varie   3 0.009740 0.003247  1.2606   0.29165    
## gesso:prof    1 0.003934 0.003934  1.5275   0.21916    
## gesso:varie   3 0.003474 0.001158  0.4496   0.71806    
## prof:varie    3 0.030415 0.010138  3.9362   0.01041 *  
## Residuals   108 0.278168 0.002576                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#--------------------------------------------
# Desdobramento.

# Desdobrar prof dentro de varie.
LSmeans(m0, effect = c("prof", "varie"))
##    estimate         se  df   t.stat      p.value      prof          varie
## 1 0.3742672 0.01268767 108 29.49850 1.705504e-53    0-0.05 IAC Fantástico
## 2 0.3723291 0.01268767 108 29.34574 2.809382e-53 0.05-0.20 IAC Fantástico
## 3 0.3935427 0.01268767 108 31.01773 1.326966e-55    0-0.05       Imperial
## 4 0.3532441 0.01314503 108 26.87283 1.214396e-49 0.05-0.20       Imperial
## 5 0.4098936 0.01268767 108 32.30645 2.498116e-57    0-0.05         Pérola
## 6 0.3270291 0.01268767 108 25.77534 6.006926e-48 0.05-0.20         Pérola
## 7 0.3999213 0.01268767 108 31.52047 2.773782e-56    0-0.05 Smooth cayenne
## 8 0.3317214 0.01268767 108 26.14518 1.591868e-48 0.05-0.20 Smooth cayenne
Xm <- LSmatrix(m0, effect = c("prof", "varie"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES = grid$varie, FUN = as.matrix)
L <- lapply(L, "rownames<-", levels(grid$prof))
g2 <- lapply(L, apmc, model = m0, focus = "prof")
g2
## $`IAC Fantástico`
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.3742672 0.3491181 0.3994164   a
## 2 0.05-0.20 0.3723291 0.3471799 0.3974782   a
## 
## $Imperial
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.3935427 0.3683936 0.4186919   a
## 2 0.05-0.20 0.3532441 0.3271884 0.3792998   b
## 
## $Pérola
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.4098936 0.3847445 0.4350428   a
## 2 0.05-0.20 0.3270291 0.3018799 0.3521782   b
## 
## $`Smooth cayenne`
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.3999213 0.3747722 0.4250705   a
## 2 0.05-0.20 0.3317214 0.3065722 0.3568705   b
g2 <- ldply(g2); names(g2)[1] <- "varie"
g2 <- equallevels(g2, db)
#--------------------------------------------
# Gráfico.

p2 <- segplot(varie ~ lwr + upr, centers = fit, data = g2, draw = FALSE,
              groups = g2$prof, pch = g2$prof, gap = 0.1,
              panel = panel.groups.segplot,
              xlab = expression("Saturation moisture"~(m^3~m^{-3})),
              ylab = expression(Gypsum~(ton~ha^{-1})),
              key = list(columns = 2, type = "o", divide = 1,
                  title = "Soil layer (m)", cex.title = 1.1,
                  lines = list(pch = 1:2),
                  text = list(levels(g2$prof)))) +
    layer(panel.text(x = centers,
                     y = as.numeric(z) +
                         gap * (2 * ((as.numeric(groups) - 1)/
                                     (nlevels(groups) - 1)) - 1),
                     labels = cld, pos = 3), data = g2)
p2

alpha: parâmetro empírico da curva de retenção

#-----------------------------------------------------------------------
# Análise do alpha.

# min(params$alp)
# Para atender os pressupostos, somou-se 1.5 para não ter valores
# negativos e elevou-se ao quadrado.

m0 <- lm(c(alp + 1.5)^2 ~ (cober + gesso + prof + varie)^2,
         data = params)
par(mfrow = c(2,2)); plot(m0); layout(1)
# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[, "dffit"]

m0 <- update(m0, data = params[!del,])
par(mfrow = c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: c(alp + 1.5)^2
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## cober         1 0.0170 0.01700  0.5250 0.4702859    
## gesso         1 0.0032 0.00319  0.0986 0.7541183    
## prof          1 0.9220 0.92203 28.4690 5.222e-07 ***
## varie         3 0.0348 0.01161  0.3584 0.7831649    
## cober:gesso   1 0.2338 0.23375  7.2174 0.0083508 ** 
## cober:prof    1 0.0018 0.00181  0.0559 0.8136096    
## cober:varie   3 0.2634 0.08779  2.7107 0.0485873 *  
## gesso:prof    1 0.0394 0.03940  1.2165 0.2724814    
## gesso:varie   3 0.1655 0.05518  1.7038 0.1704988    
## prof:varie    3 0.5929 0.19765  6.1026 0.0007069 ***
## Residuals   109 3.5302 0.03239                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Como o alpha é um parâmetro empírico, de quase nenhum significado
# aplicado, não será feito o desdobramento da interação. As diferenças
# em alpha podem implicar em diferenças no parâmetro S que é mais
# interpretável. Abaixo segem as estimativas.

# Os valores abaixo são para simples conferência. Lembrar que estão na
# escala transformada y_t = (y + 1.5)^2.

LSmeans(m0, effect = c("prof", "varie"))
##    estimate         se  df   t.stat      p.value      prof          varie
## 1 0.6336620 0.04499117 109 14.08414 2.741151e-26    0-0.05 IAC Fantástico
## 2 0.6922503 0.04499117 109 15.38636 4.370944e-29 0.05-0.20 IAC Fantástico
## 3 0.7508621 0.04499117 109 16.68910 8.507013e-32    0-0.05       Imperial
## 4 0.5567925 0.04499117 109 12.37559 1.676863e-22 0.05-0.20       Imperial
## 5 0.8280539 0.04499117 109 18.40481 3.166319e-35    0-0.05         Pérola
## 6 0.5392821 0.04499117 109 11.98640 1.265087e-21 0.05-0.20         Pérola
## 7 0.7654165 0.04499117 109 17.01259 1.865658e-32    0-0.05 Smooth cayenne
## 8 0.5106871 0.04499117 109 11.35083 3.500925e-20 0.05-0.20 Smooth cayenne
LSmeans(m0, effect = c("gesso", "cober"))
##    estimate         se  df   t.stat      p.value gesso     cober
## 1 0.6334119 0.03181356 109 19.91012 4.216282e-38     0    Millet
## 2 0.7088900 0.03181356 109 22.28263 2.179613e-42     4    Millet
## 3 0.6958293 0.03181356 109 21.87210 1.146094e-41     0 No millet
## 4 0.6003721 0.03181356 109 18.87158 3.941312e-36     4 No millet
LSmeans(m0, effect = c("cober", "varie"))
##    estimate         se  df   t.stat      p.value     cober          varie
## 1 0.6650487 0.04499117 109 14.78176 8.489400e-28    Millet IAC Fantástico
## 2 0.6608636 0.04499117 109 14.68874 1.345018e-27 No millet IAC Fantástico
## 3 0.5951633 0.04499117 109 13.22845 2.087136e-24    Millet       Imperial
## 4 0.7124913 0.04499117 109 15.83625 4.945390e-30 No millet       Imperial
## 5 0.7397585 0.04499117 109 16.44230 2.731014e-31    Millet         Pérola
## 6 0.6275775 0.04499117 109 13.94890 5.409373e-26 No millet         Pérola
## 7 0.6846332 0.04499117 109 15.21706 9.987280e-29    Millet Smooth cayenne
## 8 0.5914704 0.04499117 109 13.14637 3.174541e-24 No millet Smooth cayenne

n: parâmetro empírico da curva de retenção

#-----------------------------------------------------------------------
# Análise do n (na escala log para atender pressupostos).

m0 <- lm(log(n) ~ (cober + gesso + prof + varie)^2, data = params)
par(mfrow = c(2,2)); plot(m0); layout(1)

# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[, "dffit"]

m0 <- update(m0, data = params[!del,])
par(mfrow = c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: log(n)
##              Df Sum Sq  Mean Sq F value  Pr(>F)  
## cober         1 0.0000 0.000044  0.0010 0.97504  
## gesso         1 0.0185 0.018477  0.4112 0.52277  
## prof          1 0.1573 0.157337  3.5010 0.06409 .
## varie         3 0.0586 0.019536  0.4347 0.72860  
## cober:gesso   1 0.1855 0.185518  4.1281 0.04468 *
## cober:prof    1 0.0232 0.023209  0.5164 0.47394  
## cober:varie   3 0.2020 0.067343  1.4985 0.21923  
## gesso:prof    1 0.0864 0.086430  1.9232 0.16841  
## gesso:varie   3 0.1008 0.033607  0.7478 0.52597  
## prof:varie    3 0.1503 0.050088  1.1145 0.34661  
## Residuals   106 4.7637 0.044941                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Desdobramento não realizado por argumentos iguais aos do parâmetro
# alpha.

LSmeans(m0, effect = c("gesso", "cober"))
##   estimate         se  df   t.stat      p.value gesso     cober
## 1 1.505749 0.03815674 106 39.46220 3.408275e-65     0    Millet
## 2 1.400516 0.03815674 106 36.70429 4.388155e-62     4    Millet
## 3 1.425495 0.03747520 106 38.03835 1.299252e-63     0 No millet
## 4 1.476608 0.03815674 106 38.69850 2.367706e-64     4 No millet

Parâmetro S: taxa da curva de retenção no ponto de inflexão

#-----------------------------------------------------------------------
# Parâmetro S: slope da curva de retenção no ponto de inflexão.

# Calcular o S a partir das estimativas dos parâmetros da CRA.
params$S <- with(params, {
    m <- 1 - 1/n
    d <- Us - Ur
    -d * n * (1 + 1/m)^(-m - 1)
})

m0 <- lm(S ~ (cober + gesso + prof + varie)^2, data = params)
par(mfrow = c(2,2)); plot(m0); layout(1)

# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[,"dffit"]

m0 <- update(m0, data = params[!del, ])
par(mfrow = c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: S
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## cober         1 0.00161 0.0016143  0.3199 0.5728
## gesso         1 0.00008 0.0000810  0.0160 0.8994
## prof          1 0.01025 0.0102480  2.0310 0.1570
## varie         3 0.00700 0.0023324  0.4622 0.7092
## cober:gesso   1 0.00001 0.0000115  0.0023 0.9621
## cober:prof    1 0.00463 0.0046337  0.9183 0.3401
## cober:varie   3 0.00400 0.0013327  0.2641 0.8511
## gesso:prof    1 0.00183 0.0018293  0.3625 0.5484
## gesso:varie   3 0.02074 0.0069133  1.3701 0.2558
## prof:varie    3 0.01092 0.0036408  0.7215 0.5412
## Residuals   108 0.54495 0.0050458
# Não existe efeito dos fatores experimentais sobre o S.
LSmeans(m0)
##            estimate          se  df    t.stat      p.value
## estimate -0.2733963 0.006307294 108 -43.34606 3.978678e-70

I: tensão correspondente ao ponto de inflexão

#-----------------------------------------------------------------------
# Parâmetro I: tensão correspondente ao ponto de inflexão.

# Calcular o I a partir das estimativas dos parâmetros da CRA.
params$I <- with(params, {
    m <- 1 - 1/n
    -alp - log(m)/n
})

# Análise sobre o recíproco de I para atender os pressupostos.
m0 <- lm(1/I ~ (cober + gesso + prof + varie)^2, data = params)
par(mfrow = c(2,2)); plot(m0); layout(1)
# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[, "dffit"]

m0 <- update(m0, data = params[!del, ])
par(mfrow = c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: 1/I
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## cober         1 0.0076 0.00764  0.2179  0.641552    
## gesso         1 0.0014 0.00142  0.0406  0.840757    
## prof          1 0.9900 0.99000 28.2402 5.735e-07 ***
## varie         3 0.0419 0.01396  0.3982  0.754568    
## cober:gesso   1 0.1745 0.17452  4.9781  0.027716 *  
## cober:prof    1 0.0202 0.02022  0.5768  0.449221    
## cober:varie   3 0.2606 0.08687  2.4780  0.065085 .  
## gesso:prof    1 0.0536 0.05364  1.5302  0.218736    
## gesso:varie   3 0.1980 0.06601  1.8829  0.136728    
## prof:varie    3 0.4783 0.15943  4.5477  0.004829 ** 
## Residuals   109 3.8211 0.03506                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#--------------------------------------------
# Desdobramentos.

# Desdobrar prof dentro de varie.
LSmeans(m0, effect = c("prof", "varie"))
##   estimate         se  df   t.stat      p.value      prof          varie
## 1 1.319915 0.04680837 109 28.19826 6.830238e-52    0-0.05 IAC Fantástico
## 2 1.339459 0.04680837 109 28.61579 1.662303e-52 0.05-0.20 IAC Fantástico
## 3 1.415107 0.04680837 109 30.23191 8.099971e-55    0-0.05       Imperial
## 4 1.238376 0.04680837 109 26.45629 2.961341e-49 0.05-0.20       Imperial
## 5 1.520414 0.04680837 109 32.48166 7.035516e-58    0-0.05         Pérola
## 6 1.210782 0.04680837 109 25.86678 2.471636e-48 0.05-0.20         Pérola
## 7 1.436758 0.04680837 109 30.69445 1.838995e-55    0-0.05 Smooth cayenne
## 8 1.200014 0.04680837 109 25.63673 5.710853e-48 0.05-0.20 Smooth cayenne
Xm <- LSmatrix(m0, effect = c("prof", "varie"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES = grid$varie, FUN = as.matrix)
L <- lapply(L, "rownames<-", levels(grid$prof))
g1 <- lapply(L, apmc, model = m0, focus = "prof")
g1
## $`IAC Fantástico`
##        prof      fit      lwr      upr cld
## 1    0-0.05 1.319915 1.227142 1.412688   a
## 2 0.05-0.20 1.339459 1.246686 1.432231   a
## 
## $Imperial
##        prof      fit      lwr      upr cld
## 1    0-0.05 1.415107 1.322334 1.507879   a
## 2 0.05-0.20 1.238376 1.145603 1.331149   b
## 
## $Pérola
##        prof      fit      lwr      upr cld
## 1    0-0.05 1.520414 1.427641 1.613187   a
## 2 0.05-0.20 1.210782 1.118009 1.303554   b
## 
## $`Smooth cayenne`
##        prof      fit      lwr      upr cld
## 1    0-0.05 1.436758 1.343985 1.529530   a
## 2 0.05-0.20 1.200014 1.107241 1.292787   b
g1 <- ldply(g1); names(g1)[1] <- "varie"
g1 <- equallevels(g1, db)
g1[, c("fit","lwr","upr")] <- exp(1/g1[, c("fit","lwr","upr")])

# Desdobrar gesso dentro de cobertura.
LSmeans(m0, effect = c("gesso", "cober"))
##   estimate         se  df   t.stat      p.value gesso     cober
## 1 1.309237 0.03309852 109 39.55577 1.763911e-66     0    Millet
## 2 1.376420 0.03309852 109 41.58554 1.055777e-68     4    Millet
## 3 1.367635 0.03309852 109 41.32011 2.035939e-68     0 No millet
## 4 1.287120 0.03309852 109 38.88753 1.000287e-65     4 No millet
Xm <- LSmatrix(m0, effect = c("gesso", "cober"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES = grid$cober, FUN = as.matrix)
L <- lapply(L, "rownames<-", levels(grid$gesso))
g3 <- lapply(L, apmc, model = m0, focus = "gesso")
g3
## $Millet
##   gesso      fit      lwr      upr cld
## 1     0 1.309237 1.243637 1.374838   a
## 2     4 1.376420 1.310820 1.442020   a
## 
## $`No millet`
##   gesso      fit      lwr      upr cld
## 1     0 1.367635 1.302034 1.433235   a
## 2     4 1.287120 1.221520 1.352720   a
g3 <- ldply(g3); names(g3)[1] <- "cober"
g3 <- equallevels(g3, db)
g3[, c("fit","lwr","upr")] <- exp(1/g3[, c("fit","lwr","upr")])
#--------------------------------------------
# Gráfico.

p1 <- segplot(varie ~ lwr + upr, centers = fit, data = g1, draw = FALSE,
              groups = g1$prof, pch = g1$prof, gap = 0.1,
              panel = panel.groups.segplot,
              # xlab = expression(log ~ da ~ tensão ~ (kPa)),
              xlab = expression(Tension ~ (kPa)),
              ylab = "Variety",
              key = list(columns = 2, type = "o", divide = 1,
                  title = "Soil layer (m)", cex.title = 1.1,
                  lines = list(pch = 1:2),
                  text = list(levels(g1$prof)))) +
    layer(panel.text(x = centers,
                     y = as.numeric(z) +
                         gap * (2 * ((as.numeric(groups) - 1)/
                                     (nlevels(groups) - 1)) - 1),
                     labels = cld, pos = 3), data = g1)

p3 <- segplot(cober ~ lwr + upr, centers = fit, data = g3, draw = FALSE,
              groups = g3$gesso, pch = g3$gesso, gap = 0.1,
              panel = panel.groups.segplot,
              # xlab = expression(log ~ da ~ tensão ~ (kPa)),
              xlab = expression(Tension ~ (kPa)),
              ylab = "Crop cover",
              key = list(columns = 2, type = "o", divide = 1,
                         title = expression(Gypsum ~ (ton ~ ha^{-1})),
                         cex.title = 1.1,
                         lines = list(pch = 1:2),
                         text = list(levels(g3$gesso)))) +
    layer(panel.text(x = centers,
                     y = as.numeric(z) +
                         gap * (2 * ((as.numeric(groups) - 1)/
                                     (nlevels(groups) - 1)) - 1),
                     labels = cld, pos = 3), data = g3)

plot(p3, position = c(0,0.6,1,1), more = TRUE)
plot(p1, position = c(0,0,1,0.6), more = FALSE)

Parâmetro CAD: Capacidade de água disponível

#-----------------------------------------------------------------------
# Parâmetro CAD: Conteúdo de água disponível (CAD  =  UI-Ur).

# Calcular a umidade na tensão e a diferência de umidade com relação à
# residual.
params$UI <- with(params, {
    UI <- Ur + (Us - Ur)/(1 + exp(n*(alp + I)))^(1 - 1/n)
    UI-Ur
})

m0 <- lm(UI ~ (cober + gesso + prof + varie)^2, data = params)
par(mfrow = c(2,2)); plot(m0); layout(1)

# MASS::boxcox(m0)

im <- influence.measures(m0)
## Warning in abbreviate(vn): abbreviate used with non-ASCII chars

## Warning in abbreviate(vn): abbreviate used with non-ASCII chars
# summary(im)
del <- im$is.inf[, "dffit"]

m0 <- update(m0, data = params[!del,])
par(mfrow = c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: UI
##              Df   Sum Sq   Mean Sq F value    Pr(>F)    
## cober         1 0.000885 0.0008848  1.2490  0.266220    
## gesso         1 0.000937 0.0009368  1.3225  0.252689    
## prof          1 0.015630 0.0156296 22.0631 7.807e-06 ***
## varie         3 0.001306 0.0004352  0.6144  0.607131    
## cober:gesso   1 0.005961 0.0059612  8.4150  0.004510 ** 
## cober:prof    1 0.000659 0.0006589  0.9301  0.336995    
## cober:varie   3 0.002605 0.0008682  1.2256  0.304028    
## gesso:prof    1 0.000801 0.0008011  1.1308  0.289966    
## gesso:varie   3 0.000626 0.0002086  0.2945  0.829282    
## prof:varie    3 0.011551 0.0038503  5.4351  0.001611 ** 
## Residuals   108 0.076507 0.0007084                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#--------------------------------------------
# Desdobramentos.

# Desdobrar prof dentro de varie.
LSmeans(m0, effect = c("prof", "varie"))
##    estimate          se  df   t.stat      p.value      prof          varie
## 1 0.1467745 0.006653957 108 22.05822 8.413805e-42    0-0.05 IAC Fantástico
## 2 0.1552068 0.006653957 108 23.32548 5.694644e-44 0.05-0.20 IAC Fantástico
## 3 0.1607525 0.006653957 108 24.15893 2.350935e-45    0-0.05       Imperial
## 4 0.1405351 0.006893816 108 20.38568 8.191351e-39 0.05-0.20       Imperial
## 5 0.1654711 0.006653957 108 24.86807 1.656559e-46    0-0.05         Pérola
## 6 0.1250955 0.006653957 108 18.80016 7.643154e-36 0.05-0.20         Pérola
## 7 0.1609968 0.006653957 108 24.19565 2.046577e-45    0-0.05 Smooth cayenne
## 8 0.1256906 0.006653957 108 18.88960 5.153559e-36 0.05-0.20 Smooth cayenne
Xm <- LSmatrix(m0, effect = c("prof", "varie"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES = grid$varie, FUN = as.matrix)
L <- lapply(L, "rownames<-", levels(grid$prof))
g1 <- lapply(L, apmc, model = m0, focus = "prof")
g1
## $`IAC Fantástico`
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.1467745 0.1335852 0.1599638   a
## 2 0.05-0.20 0.1552068 0.1420175 0.1683961   a
## 
## $Imperial
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.1607525 0.1475632 0.1739418   a
## 2 0.05-0.20 0.1405351 0.1268704 0.1541998   b
## 
## $Pérola
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.1654711 0.1522818 0.1786603   a
## 2 0.05-0.20 0.1250955 0.1119062 0.1382848   b
## 
## $`Smooth cayenne`
##        prof       fit       lwr       upr cld
## 1    0-0.05 0.1609968 0.1478075 0.1741861   a
## 2 0.05-0.20 0.1256906 0.1125013 0.1388799   b
g1 <- ldply(g1)
names(g1)[1] <- "varie"
g1 <- equallevels(g1, db)

# Desdobrar gesso dentro de cobertura.
LSmeans(m0, effect = c("gesso", "cober"))
##    estimate          se  df   t.stat      p.value gesso     cober
## 1 0.1406845 0.004705058 108 29.90068 4.627503e-54     0    Millet
## 2 0.1489789 0.004705058 108 31.66357 1.783083e-56     4    Millet
## 3 0.1600236 0.004790612 108 33.40358 9.384998e-59     0 No millet
## 4 0.1405744 0.004705058 108 29.87729 4.990420e-54     4 No millet
Xm <- LSmatrix(m0, effect = c("gesso", "cober"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES = grid$cober, FUN = as.matrix)
L <- lapply(L, "rownames<-", levels(grid$gesso))
g3 <- lapply(L, apmc, model = m0, focus = "gesso")
g3
## $Millet
##   gesso       fit       lwr       upr cld
## 1     0 0.1406845 0.1313582 0.1500107   a
## 2     4 0.1489789 0.1396527 0.1583052   a
## 
## $`No millet`
##   gesso       fit       lwr       upr cld
## 1     0 0.1600236 0.1505278 0.1695194   a
## 2     4 0.1405744 0.1312481 0.1499006   b
g3 <- ldply(g3)
names(g3)[1] <- "cober"
g3 <- equallevels(g3, db)
#--------------------------------------------
# Gráfico.

p1 <- segplot(varie ~ lwr + upr, centers = fit, data = g1, draw = FALSE,
              groups = g1$prof, pch = g1$prof, gap = 0.1,
              panel = panel.groups.segplot,
              xlab = expression(
                  "Available water content"~ (m^3~m^{-3})),
              ylab = "Variety",
              key = list(columns = 2, type = "o", divide = 1,
                  title = "Soil layer (m)", cex.title = 1.1,
                  lines = list(pch = 1:2),
                  text = list(levels(g1$prof)))) +
    layer(panel.text(x = centers,
                     y = as.numeric(z) +
                         gap * (2 * ((as.numeric(groups) - 1)/
                                     (nlevels(groups) - 1)) - 1),
                     labels = cld, pos = 3), data = g1)

p3 <- segplot(cober ~ lwr + upr, centers = fit, data = g3, draw = FALSE,
              groups = g3$gesso, pch = g3$gesso, gap = 0.1,
              panel = panel.groups.segplot,
              xlab = expression(
                  "Available water content" ~ (m^3 ~ m^{-3})),
              ylab = "Crop cover",
              key = list(columns = 2, type = "o", divide = 1,
                         title = expression(Gypsum ~ (ton ~ ha^{-1})),
                         cex.title = 1.1,
                  lines = list(pch = 1:2),
                  text = list(levels(g3$gesso)))) +
    layer(panel.text(x = centers,
                     y = as.numeric(z) +
                         gap * (2 * ((as.numeric(groups) - 1)/
                                     (nlevels(groups) - 1)) - 1),
                     labels = cld, pos = 3), data = g3)

plot(p3, position = c(0,0.6,1,1), more = TRUE)
plot(p1, position = c(0,0,1,0.6), more = FALSE)

#-----------------------------------------------------------------------
# Informações sobre as versões dos pacotes.

sessionInfo()
## R version 3.3.0 (2016-05-03)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C               LC_TIME=pt_BR.UTF-8       
##  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] methods   stats     graphics  grDevices utils     datasets  base     
## 
## other attached packages:
##  [1] wzRfun_0.65         multcomp_1.4-5      TH.data_1.0-7       MASS_7.3-45        
##  [5] survival_2.39-4     mvtnorm_1.0-5       doBy_4.5-15         plyr_1.8.4         
##  [9] reshape_0.8.5       nlme_3.1-128        latticeExtra_0.6-28 RColorBrewer_1.1-2 
## [13] lattice_0.20-33     rmarkdown_0.9.6     knitr_1.13         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.3      magrittr_1.5     splines_3.3.0    stringr_1.0.0    tools_3.3.0     
##  [6] rpanel_1.1-3     grid_3.3.0       htmltools_0.3.5  yaml_2.1.13      digest_0.6.9    
## [11] Matrix_1.2-6     formatR_1.4      codetools_0.2-14 evaluate_0.9     sandwich_2.3-4  
## [16] stringi_1.1.1    zoo_1.7-13