Biocarvão e fósforo na Curva de Retenção de Água
Walmes Zeviani & Milson Serafim
Descrição do experimento
Este experimento estudou o efeito do carvão enriquecido com fósforo no desenvolvimento de mudas de eucalipto. Foram estudados 3 fontes do biocarvão (biocar):
- Biocarvão enriquecido com fósforo adicionado ao solo (S+(B+P) = S_BP).
- Biocarvão não enriquecido com fósforo mas com o fósforo adicionado ao solo separado ((S+P)+B = SP_B).
- Solo sem biocarvão mas com o fósforo (S+P = SP).
A quantidade de fósforo usada também foi um fator, com 10 níveis, praticamente igualmente espaçados na escala log base 2.
O experimento foi feito em delineamento de blocos casualizados completos em casa de vegetação. As unidades experimentais foram um vaso com uma muda de eucalipto em cada.
Para cada unidade experimental, amostras de solo indeformadas foram levadas ao laboratório para determinação da curva de retenção de água do solo. Foram aplicadas as tensões 0, 1, 2, 4, 6, 8, 10, 33, 66, 100, 300, 1500 kPa. Para cada unidade experimental ajustou-se o modelo van Genuchten para retenção de água do solo. A partir das estimativas dos parâmetros do modelo, foram calculados os parâmetros S (inclinação na inflexão) e I (tensão do ponto de inflexão). As estimativas dos parâmetros Ur (umidade residual), Us (umidade de satuação), S e I foram submetidos a análise de variância. O estudo das médias preconizou o desdobramento do efeito de biocarvão para doses fixas de P ou o estudo dos efeitos principais, caso não houvesse interação.
Por os dados serem desbalanceados, por parcelas onde não determinou-se a curva ou onde o modelo não ajustou, as comparações múltiplas não foram feitas pelo teste de Tukey mas por contrastes com p-valor corrigido pelo método single-step.
Os parâmetros a e n da curva de retenção não foram analisados por serem empíricos, com pouco significado, e porque os parâmetros S e I são dependentes destes.
#-----------------------------------------------------------------------
# Pacotes.
library(gdata)
library(lattice)
library(latticeExtra)
# library(splines)
library(doBy)
library(multcomp)
library(agricolae)
library(plyr)
source("~/Dropbox/templates/_setup.R")
opts_chunk$set(fig.width = 7,
fig.height = 5)
# "https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.cbH.R",
# "https://raw.githubusercontent.com/walmes/wzRfun/master/R/prepanel.cbH.R",
# "https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.segplotBy.R",
url <-
c("https://raw.githubusercontent.com/walmes/wzRfun/master/R/apmc.R",
"https://raw.githubusercontent.com/walmes/wzRfun/master/R/apc.R")
sapply(url, source)
#-----------------------------------------------------------------------
panel.segplotBy <- function(x, y, z,
centers, subscripts,
groups, f, letters, digits, ...){
if (!missing(groups)) {
d <- 2 * ((as.numeric(groups) - 1)/(nlevels(groups) - 1)) - 1
z <- as.numeric(z) + f * d
}
panel.segplot(x, y, z, centers = centers,
subscripts = subscripts, ...)
panel.text(x = z, y = centers,
labels = sprintf(paste0("%0.", digits, "f %s"),
centers, letters),
cex = 0.8, srt = 90, adj = c(0.5, -0.5))
}
#-----------------------------------------------------------------------
# system("file -b biocarvao.xlsx")
da <- read.xls("biocarvao.xlsx", sheet = 1,
encoding = "latin1")
da <- da[, !sapply(da, FUN = is.logical)]
da <- subset(da, select = 1:5)
str(da)
# dput(gsub(x = iconv(tolower(names(da)), to = "ASCII//TRANSLIT"),
# pattern = "\\..*$", replacement = ""))
names(da) <-
c("biocar", "dosep", "tensao", "bloco", "umid")
levels(da$biocar) <- c("S_BP", "SP_B", "SP")
da <- da[, c(1, 2, 4, 3, 5)]
str(da)
db <- droplevels(da[!is.na(da$umid) & da$umid > 0, ])
str(db)Curvas de Retenção de Água
#-----------------------------------------------------------------------
# Tabela de frequência das combinações.
xtabs(~dosep + biocar, data = db)## biocar
## dosep SP_B SP
## 0 36 36
## 1.5 36 36
## 3 36 36
## 6 24 36
## 12 36 36
## 24 36 36
## 48 36 36
## 96 36 36
## 192 36 36
## 348 36 36xtabs(~tensao + biocar, data = db)## biocar
## tensao SP_B SP
## 0 29 30
## 1 29 30
## 2 29 30
## 4 29 30
## 6 29 30
## 8 29 30
## 10 29 30
## 33 29 30
## 66 29 30
## 100 29 30
## 300 29 30
## 1500 29 30xtabs(~bloco + biocar, data = db)## biocar
## bloco SP_B SP
## 1 108 120
## 2 120 120
## 3 120 120# Cria a unidade experimental da onde são feitas as curvas.
db$ue <- with(db, interaction(bloco, biocar, dosep,
drop = TRUE, sep = ":"))
nlevels(db$ue)## [1] 59# Troca a tensão 0 por 0.1 para não virar -Inf ao tirar o log.
db$tensao[db$tensao == 0] <- 0.1
#-----------------------------------------------------------------------
xyplot(umid ~ tensao | ue, data = db,
as.table = TRUE, strip = FALSE,
ylab = "Umidade do solo",
xlab = "Tensão matricial",
scales = list(x = list(log = 10)),
xscale.components = xscale.components.log10ticks)
xyplot(umid ~ tensao | as.factor(dosep),
groups = biocar, data = db,
as.table = TRUE,
auto.key = list(title = "Biocarvão", cex.title = 1,
columns = 2),
ylab = "Umidade do solo",
xlab = "Tensão matricial",
scales = list(x = list(log = 10)),
xscale.components = xscale.components.log10ticks)
db$ltens <- log10(db$tensao)#-----------------------------------------------------------------------
# Ajuste da curva por unidade experimental.
library(rpanel)
source(paste0("https://raw.githubusercontent.com/walmes/wzRfun",
"/master/R/rp.nls.R"))
model <- umid ~ Ur + (Us - Ur)/
(1 + exp(n * (al + ltens)))^(1 - 1/n)
start <- list(Ur = c(init = 0.2, from = 0, to = 0.5),
Us = c(init = 0.6, from = 0.4, to = 0.8),
al = c(1, -2, 4),
n = c(1.5, 1, 4))
rp.nls(model = model,
data = db,
start = start,
subset = "ue",
assignTo = "cra.fit")
is.null(cra.fit)
params <- t(sapply(cra.fit[!sapply(cra.fit, is.null)], FUN = coef))
# rownames(params) <- gsub(x = gsub(x = rownames(params), "\\.", ":"),
# "1:5", "1.5")
str(db)
grid <- as.data.frame(do.call(rbind, strsplit(rownames(params), ":")))
names(grid) <- c("bloco", "biocar", "dosep")
grid <- transform(grid,
bloco = as.integer(bloco),
biocar = factor(biocar, levels = levels(db$biocar)),
dosep = as.numeric(dosep))
str(grid)
str(params)
params <- cbind(grid, as.data.frame(params))
str(params)xtabs(~biocar + dosep, data = params)## dosep
## biocar 1 2 3 4 5 6 7 8 9 10
## SP_B 3 3 3 3 2 3 3 3 2 2
## SP 3 3 3 3 3 3 3 3 3 3# 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)
})
# Calcular o I a partir das estimativas dos parâmetros da CRA.
params$I <- with(params, {
m <- 1 - 1/n
-al - log(m)/n
})
params$dosep <- as.numeric(gsub(x = rownames(params),
"^(.*:)(\\d.*)$", "\\2"))
params$catep <- factor(params$dosep)
params$bloco <- factor(params$bloco)
params$dosept <- log(params$dosep + 1, base = 2)
str(params)## 'data.frame': 57 obs. of 11 variables:
## $ bloco : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
## $ biocar: Factor w/ 2 levels "SP_B","SP": 1 1 1 2 2 2 1 1 1 2 ...
## $ dosep : num 0 0 0 0 0 0 1.5 1.5 1.5 1.5 ...
## $ Ur : num 0.118 0.103 0.109 0.102 0.111 ...
## $ Us : num 0.431 0.469 0.455 0.447 0.484 ...
## $ n : num 4.51 2.88 3.53 4.25 4.95 ...
## $ al : num -0.662 -0.208 -0.565 -0.562 -0.637 ...
## $ S : num -0.325 -0.227 -0.273 -0.335 -0.427 ...
## $ I : num 0.717 0.357 0.66 0.625 0.683 ...
## $ catep : Factor w/ 10 levels "0","1.5","3",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ dosept: num 0 0 0 0 0 ...Umidade residual (Ur)
xyplot(Ur ~ dosept, groups = biocar, data = params,
type = c("p", "a"))
m0 <- lm(Ur ~ bloco + biocar * catep, data = params)
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)
## bloco 2 0.0003039 0.00015193 1.2133 0.3094
## biocar 1 0.0000827 0.00008275 0.6608 0.4218
## catep 9 0.0078541 0.00087267 6.9692 1.086e-05 ***
## biocar:catep 9 0.0070368 0.00078186 6.2440 3.198e-05 ***
## Residuals 35 0.0043826 0.00012522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(params$biocar))
grid$catep <- factor(grid$catep, levels = levels(params$catep))
Xs <- by(L, INDICES = grid$catep, FUN = as.matrix)
Xs <- lapply(Xs, "rownames<-", levels(grid$biocar))
grid <- ldply(.id = "catep",
lapply(Xs, FUN = apmc, model = m0,
focus = "biocar", test = "single-step"))
xlab <- expression("Dose de P"~(mg~dm^{-3}))
ylab <- expression("Umidade de satuação"~(m^3~m^{-3}))
key <- list(type = "o", divide = 1, columns = 2,
text = list(levels(grid$biocar)),
lines = list(pch = c(1, 5)))
str(grid)## 'data.frame': 20 obs. of 6 variables:
## $ catep : Factor w/ 10 levels "0","1.5","3",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ biocar : Factor w/ 2 levels "SP","SP_B": 2 1 2 1 2 1 2 1 2 1 ...
## $ estimate: num 0.1099 0.1064 0.1037 0.1148 0.0981 ...
## $ lwr : num 0.0968 0.0933 0.0906 0.1017 0.085 ...
## $ upr : num 0.123 0.12 0.117 0.128 0.111 ...
## $ cld : chr "a" "a" "a" "a" ...segplot(catep ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
groups = biocar,
pch = key$lines$pch,
key = key, xlab = xlab, ylab = ylab,
f = 0.15, panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid## catep biocar estimate lwr upr cld
## 1 0 SP_B 0.10990614 0.09679044 0.1230218 a
## 2 0 SP 0.10640591 0.09329022 0.1195216 a
## 3 1.5 SP_B 0.10367659 0.09056090 0.1167923 a
## 4 1.5 SP 0.11477381 0.10165811 0.1278895 a
## 5 3 SP_B 0.09811445 0.08499875 0.1112301 b
## 6 3 SP 0.11860788 0.10549218 0.1317236 a
## 7 6 SP_B 0.09865734 0.08244374 0.1148709 b
## 8 6 SP 0.12485327 0.11173757 0.1379690 a
## 9 12 SP_B 0.10499165 0.09187595 0.1181073 a
## 10 12 SP 0.12120299 0.10808729 0.1343187 a
## 11 24 SP_B 0.09580477 0.07959117 0.1120184 b
## 12 24 SP 0.11936748 0.10625178 0.1324832 a
## 13 48 SP_B 0.10861255 0.09549686 0.1217283 a
## 14 48 SP 0.11876153 0.10564584 0.1318772 a
## 15 96 SP_B 0.11992613 0.10371836 0.1361339 a
## 16 96 SP 0.10540930 0.09229361 0.1185250 a
## 17 192 SP_B 0.12651433 0.11339863 0.1396300 a
## 18 192 SP 0.12395460 0.11083890 0.1370703 a
## 19 348 SP_B 0.17348752 0.16037182 0.1866032 b
## 20 348 SP 0.12141413 0.10829843 0.1345298 aUmidade de satuação (Us)
xyplot(Us ~ dosept, groups = biocar, data = params,
type = c("p", "a"))
m0 <- lm(Us ~ bloco + biocar * catep, data = params)
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)
## bloco 2 0.0021049 0.00105245 1.8587 0.17094
## biocar 1 0.0000133 0.00001333 0.0235 0.87893
## catep 9 0.0124833 0.00138703 2.4496 0.02789 *
## biocar:catep 9 0.0106903 0.00118781 2.0977 0.05681 .
## Residuals 35 0.0198183 0.00056624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1L <- LSmatrix(m0, effect = c("catep"))
grid <- attr(L, "grid")
grid$catep <- factor(grid$catep, levels = levels(params$catep))
rownames(L) <- levels(grid$catep)
grid <- apmc(L, model = m0, focus = "catep", test = "fdr")
grid$catep <- factor(grid$catep, levels = levels(params$catep))
xlab <- expression("Dose de P"~(mg~dm^{-3}))
ylab <- expression("Umidade residual"~(m^3~m^{-3}))
segplot(catep ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
xlab = xlab, ylab = ylab,
panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid## catep estimate lwr upr cld
## 1 0 0.4593794 0.4396578 0.4791010 b
## 2 1.5 0.4509781 0.4312565 0.4706997 b
## 3 3 0.4760507 0.4563291 0.4957724 ab
## 4 6 0.4671370 0.4449636 0.4893103 ab
## 5 12 0.4716107 0.4518891 0.4913324 ab
## 6 24 0.4540342 0.4318608 0.4762076 b
## 7 48 0.4673954 0.4476738 0.4871170 ab
## 8 96 0.4561583 0.4339898 0.4783269 b
## 9 192 0.4739578 0.4542362 0.4936794 ab
## 10 348 0.5058694 0.4861478 0.5255910 aInclinação máxima (na inflexão) da CRA (S)
xyplot(S ~ dosept, groups = biocar, data = params,
type = c("p", "a"))
m0 <- lm(S ~ bloco + biocar * catep, data = params)
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)
## bloco 2 0.014175 0.0070875 1.3785 0.26530
## biocar 1 0.000866 0.0008664 0.1685 0.68394
## catep 9 0.121015 0.0134461 2.6152 0.01999 *
## biocar:catep 9 0.088828 0.0098698 1.9196 0.08136 .
## Residuals 35 0.179952 0.0051415
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1L <- LSmatrix(m0, effect = c("catep"))
grid <- attr(L, "grid")
grid$catep <- factor(grid$catep, levels = levels(params$catep))
rownames(L) <- levels(grid$catep)
grid <- apmc(L, model = m0, focus = "catep", test = "fdr")
grid$catep <- factor(grid$catep, levels = levels(params$catep))
xlab <- expression("Dose de P"~(mg~dm^{-3}))
ylab <- expression("Inclinação na inflexão (S)")
segplot(catep ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
xlab = xlab, ylab = ylab,
panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid## catep estimate lwr upr cld
## 1 0 -0.3301326 -0.3895600 -0.2707051 a
## 2 1.5 -0.3557243 -0.4151518 -0.2962968 ab
## 3 3 -0.3530516 -0.4124790 -0.2936241 ab
## 4 6 -0.3656731 -0.4324885 -0.2988577 ab
## 5 12 -0.4527879 -0.5122154 -0.3933605 b
## 6 24 -0.3682415 -0.4350569 -0.3014261 ab
## 7 48 -0.2900491 -0.3494766 -0.2306217 a
## 8 96 -0.3217276 -0.3885285 -0.2549267 a
## 9 192 -0.2996302 -0.3590576 -0.2402027 a
## 10 348 -0.3024056 -0.3618331 -0.2429781 aPonto de inflexão da CRA (I)
xyplot(I ~ dosept, groups = biocar, data = params,
type = c("p", "a"))
m0 <- lm(I ~ bloco + biocar * catep, data = params)
par(mfrow = c(2, 2))
plot(m0); layout(1)
anova(m0)## Analysis of Variance Table
##
## Response: I
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 2 0.013738 0.006869 0.9655 0.390719
## biocar 1 0.032797 0.032797 4.6098 0.038797 *
## catep 9 0.049246 0.005472 0.7691 0.645140
## biocar:catep 9 0.218699 0.024300 3.4154 0.004146 **
## Residuals 35 0.249015 0.007115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(params$biocar))
grid$catep <- factor(grid$catep, levels = levels(params$catep))
Xs <- by(L, INDICES = grid$catep, FUN = as.matrix)
Xs <- lapply(Xs, "rownames<-", levels(grid$biocar))
grid <- ldply(.id = "catep",
lapply(Xs, FUN = apmc, model = m0,
focus = "biocar", test = "single-step"))
xlab <- expression("Dose de P"~(mg~dm^{-3}))
ylab <- expression("Tensão matricial da inflexão")
key <- list(type = "o", divide = 1, columns = 2,
text = list(levels(grid$biocar)),
lines = list(pch = c(1, 5)))
str(grid)## 'data.frame': 20 obs. of 6 variables:
## $ catep : Factor w/ 10 levels "0","1.5","3",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ biocar : Factor w/ 2 levels "SP","SP_B": 2 1 2 1 2 1 2 1 2 1 ...
## $ estimate: num 0.578 0.655 0.593 0.681 0.635 ...
## $ lwr : num 0.479 0.556 0.495 0.582 0.536 ...
## $ upr : num 0.677 0.754 0.692 0.78 0.733 ...
## $ cld : chr "a" "a" "a" "a" ...segplot(catep ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
groups = biocar,
pch = key$lines$pch,
key = key, xlab = xlab, ylab = ylab,
f = 0.15, panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid## catep biocar estimate lwr upr cld
## 1 0 SP_B 0.5779927 0.4791290 0.6768565 a
## 2 0 SP 0.6550152 0.5561515 0.7538790 a
## 3 1.5 SP_B 0.5934232 0.4945594 0.6922869 a
## 4 1.5 SP 0.6810723 0.5822086 0.7799361 a
## 5 3 SP_B 0.6345388 0.5356750 0.7334025 a
## 6 3 SP 0.5700534 0.4711897 0.6689172 a
## 7 6 SP_B 0.5865148 0.4642996 0.7087299 a
## 8 6 SP 0.6139029 0.5150392 0.7127667 a
## 9 12 SP_B 0.6041973 0.5053335 0.7030610 a
## 10 12 SP 0.6635983 0.5647346 0.7624621 a
## 11 24 SP_B 0.6534737 0.5312586 0.7756889 a
## 12 24 SP 0.6160795 0.5172157 0.7149432 a
## 13 48 SP_B 0.7057423 0.6068786 0.8046061 a
## 14 48 SP 0.3775175 0.2786537 0.4763812 b
## 15 96 SP_B 0.7246982 0.6025270 0.8468694 a
## 16 96 SP 0.5563911 0.4575274 0.6552549 b
## 17 192 SP_B 0.6611922 0.5623284 0.7600559 a
## 18 192 SP 0.6314839 0.5326201 0.7303476 a
## 19 348 SP_B 0.6585512 0.5596875 0.7574150 a
## 20 348 SP 0.5385793 0.4397156 0.6374431 a