Curva de retenção de Água na Cultura do Abacaxizeiro em Função da Cobertura do Solo e Aplicação de Gesso
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, 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"))
#-----------------------------------------------------------------------
# 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.
# 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 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)
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)
#-----------------------------------------------------------------------
# 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)
#-----------------------------------------------------------------------
# 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
#-----------------------------------------------------------------------
# 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
#-----------------------------------------------------------------------
# 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: 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
#-----------------------------------------------------------------------
# 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: 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