#-----------------------------------------------------------------------
library(lattice)
library(latticeExtra)
library(nlme)
library(doBy)
library(multcomp)
library(reshape2)
source("~/Dropbox/templates/_setup.R")
opts_chunk$set(fig.width = 7,
fig.height = 7)
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))
}
#-----------------------------------------------------------------------
# Lendo os dados.
url <- "Planilha_curva_retenção_dados_produtividade_(18-03-16).csv"
cra <- read.table(url, header = TRUE, sep = ",",
fileEncoding = "ISO-8859-1")
da <- subset(cra, Tensão..kPa. == 0)[, c(1, 3, 6:10)]
names(da) <- c("loc", "cam", "ds", "vol", "alt", "dap", "prod")
cra <- cra[, c(1, 3:5)]
names(cra) <- c("loc", "cam", "tens", "umid")
str(cra)
## 'data.frame': 1650 obs. of 4 variables:
## $ loc : int 1 1 1 1 1 1 1 1 1 1 ...
## $ cam : Factor w/ 3 levels "0 - 0,05","0,05 - 0,40",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ tens: int 0 2 4 6 8 10 33 66 100 300 ...
## $ umid: num 0.648 0.584 0.532 0.458 0.417 ...
cra$tens[cra$tens == 0] <- 0.1
#-----------------------------------------------------------------------
xyplot(umid ~ tens | factor(loc), data = cra,
groups = cam, type = c("o"), as.table = TRUE,
strip = TRUE, layout = c(NA, 5),
scales = list(x = list(log = 10)),
xscale.components = xscale.components.log10ticks,
auto.key = list(title = "Camada (cm)", cex.title = 1.1),
ylab = expression("Umidade do solo"~(m^{3}~m^{-3})),
xlab = expression(log[10]~"Tensão"~(kPa)))
# Remover curvas sem variação (impróprias).
# loc 4 cam 3
# loc 27 cam 1
# loc 37 cam 2
# loc 47 cam 2 3
del <- with(cra, {
(loc == 4 & cam == levels(cam)[3]) |
(loc == 27 & cam == levels(cam)[1]) |
(loc == 37 & cam == levels(cam)[2]) |
(loc == 47 & cam == levels(cam)[2]) |
(loc == 47 & cam == levels(cam)[3])
})
cra <- droplevels(cra[!del, ])
Ajuste da Curva de Retenção de Água do Solo
#-----------------------------------------------------------------------
# Ajuste da CRA.
xyplot(umid ~ tens , data = cra,
scales = list(x = list(log = 10)),
xscale.components = xscale.components.log10ticks,
ylab = expression("Umidade do solo"~(m^{3}~m^{-3})),
xlab = expression(log[10]~"Tensão"~(kPa)))
cra$ltens <- log10(cra$tens)
model <- umid ~ Ur + (Us - Ur)/(1 + exp(n * (alp + ltens)))^(1 - 1/n)
start <- list(Ur = 0.3, Us = 0.6, alp = -0.5, n = 4)
n00 <- nls(model, data = cra, start = start)
coef(n00)
## Ur Us alp n
## 0.2390291 0.6519457 -0.4708593 2.4380958
#-----------------------------------------------------------------------
# Ajudar para cada unidade experimental.
cra$ue <- with(cra, interaction(loc, cam, drop = TRUE))
db <- groupedData(umid ~ ltens | ue, data = cra, order.groups = FALSE)
n0 <- nlsList(model = model, data = db, start = as.list(coef(n00)))
## Warning: 1 error caught in nls(model, data = data, control =
## controlvals, start = start): singular gradient
c0 <- coef(n0)
pairs(c0)
plot(augPred(n0), strip = FALSE, as.table = TRUE,
ylab = expression("Umidade do solo"~(m^{3}~m^{-3})),
xlab = expression(log[10]~"Tensão"~(kPa)))
#-----------------------------------------------------------------------
params <- as.data.frame(do.call(rbind, strsplit(rownames(c0), "\\.")))
names(params) <- c("loc", "cam")
params <- transform(params,
loc = factor(loc),
cam = factor(cam, levels = levels(cra$cam)))
params <- na.omit(cbind(params, c0))
params$S <- with(params, {
m <- 1 - 1/n
d <- Us - Ur
-d * n * (1 + 1/m)^(-m - 1)
})
params$I <- with(params, {
m <- 1 - 1/n
-alp - log(m)/n
})
params$Ui <- with(params, {
Ur + (Us - Ur)/(1 + exp(n * (alp + I)))^(1 - 1/n)
})
params$cad <- with(params, Ui - Ur)
# with(as.list(params[1, ]), {
# curve(Ur + (Us - Ur)/(1 + exp(n * (alp + x)))^(1 - 1/n),
# from = -3, to = 6)
# abline(v = I, h = Ui, lty = 2)
# })
str(params)
## 'data.frame': 144 obs. of 10 variables:
## $ loc: Factor w/ 50 levels "1","10","11",..: 1 12 23 34 45 47 48 49 50 2 ...
## $ cam: Factor w/ 3 levels "0 - 0,05","0,05 - 0,40",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Ur : num 0.223 0.169 0.214 0.184 0.251 ...
## $ Us : num 0.648 0.622 0.667 0.66 0.576 ...
## $ alp: num -0.724 -0.578 -0.608 -0.72 -0.809 ...
## $ n : num 3.59 4.24 3.09 4.87 2.38 ...
## $ S : num -0.341 -0.439 -0.306 -0.537 -0.159 ...
## $ I : num 0.815 0.642 0.734 0.767 1.037 ...
## $ Ui : num 0.45 0.408 0.459 0.433 0.433 ...
## $ cad: num 0.227 0.239 0.245 0.249 0.182 ...
## - attr(*, "na.action")=Class 'omit' Named int 88
## .. ..- attr(*, "names")= chr "40.0,05 - 0,40"
Umidade Residual (Ur)
#-----------------------------------------------------------------------
# Ur.
m0 <- lm(Ur ~ loc + cam, data = params)
par(mfrow = c(2, 2)); plot(m0); layout(1)
## Warning: not plotting observations with leverage one:
## 46
## Warning: not plotting observations with leverage one:
## 46
anova(m0)
## Analysis of Variance Table
##
## Response: Ur
## Df Sum Sq Mean Sq F value Pr(>F)
## loc 49 0.2384 0.0048654 1.4210 0.073945 .
## cam 2 0.0406 0.0202998 5.9288 0.003785 **
## Residuals 92 0.3150 0.0034239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = "cam")
grid <- attr(L, "grid")
grid$cam <- factor(grid$cam, levels = levels(params$cam))
rownames(L) <- levels(params$cam)
grid <- apmc(X = L, model = m0, test = "single-step", focus = "cam")
xlab <- expression("Camada do solo"~(m))
ylab <- expression("Umidade residual"~(m^3~m^{-3}))
segplot(cam ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
xlab = xlab,
ylab = ylab,
f = 0.15, panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid
## cam estimate lwr upr cld
## 1 0 - 0,05 0.2186003 0.2019152 0.2352854 b
## 2 0,05 - 0,40 0.2200550 0.2027705 0.2373395 b
## 3 0,40 - 0,80 0.2551951 0.2381726 0.2722176 a
Umidade de Satuação (Us)
#-----------------------------------------------------------------------
# Us.
m0 <- lm(Us ~ loc + cam, data = params)
par(mfrow = c(2, 2)); plot(m0); layout(1)
## Warning: not plotting observations with leverage one:
## 46
## Warning: not plotting observations with leverage one:
## 46
anova(m0)
## Analysis of Variance Table
##
## Response: Us
## Df Sum Sq Mean Sq F value Pr(>F)
## loc 49 0.306960 0.0062645 1.8431 0.005864 **
## cam 2 0.039273 0.0196367 5.7774 0.004329 **
## Residuals 92 0.312698 0.0033989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = "cam")
grid <- attr(L, "grid")
grid$cam <- factor(grid$cam, levels = levels(params$cam))
rownames(L) <- levels(params$cam)
grid <- apmc(X = L, model = m0, test = "single-step", focus = "cam")
xlab <- expression("Camada do solo"~(m))
ylab <- expression("Umidade de saturação"~(m^3~m^{-3}))
segplot(cam ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
xlab = xlab,
ylab = ylab,
f = 0.15, panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid
## cam estimate lwr upr cld
## 1 0 - 0,05 0.6584775 0.6418536 0.6751014 ab
## 2 0,05 - 0,40 0.6391646 0.6219435 0.6563858 b
## 3 0,40 - 0,80 0.6801194 0.6631593 0.6970796 a
Inclinação no Ponto de Inflexão
#-----------------------------------------------------------------------
# S.
m0 <- lm(S ~ loc + cam, data = params)
par(mfrow = c(2, 2)); plot(m0); layout(1)
## Warning: not plotting observations with leverage one:
## 46
## Warning: not plotting observations with leverage one:
## 46
anova(m0)
## Analysis of Variance Table
##
## Response: S
## Df Sum Sq Mean Sq F value Pr(>F)
## loc 49 0.78671 0.016055 2.816 9.254e-06 ***
## cam 2 0.14517 0.072586 12.731 1.315e-05 ***
## Residuals 92 0.52453 0.005701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = "cam")
grid <- attr(L, "grid")
grid$cam <- factor(grid$cam, levels = levels(params$cam))
rownames(L) <- levels(params$cam)
grid <- apmc(X = L, model = m0, test = "single-step", focus = "cam")
xlab <- expression("Camada do solo"~(m))
ylab <- expression("Inclinação na inflexão (S)")
segplot(cam ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
xlab = xlab,
ylab = ylab,
f = 0.15, panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid
## cam estimate lwr upr cld
## 1 0 - 0,05 -0.2746800 -0.2962106 -0.2531494 b
## 2 0,05 - 0,40 -0.2026826 -0.2249867 -0.1803785 a
## 3 0,40 - 0,80 -0.2116415 -0.2336075 -0.1896755 a
Tensão Matricial da Inflexão da CRA (I)
#-----------------------------------------------------------------------
# I.
m0 <- lm(I ~ loc + cam, data = params)
par(mfrow = c(2, 2)); plot(m0); layout(1)
## Warning: not plotting observations with leverage one:
## 46
## Warning: not plotting observations with leverage one:
## 46
anova(m0)
## Analysis of Variance Table
##
## Response: I
## Df Sum Sq Mean Sq F value Pr(>F)
## loc 49 8.3405 0.170215 2.2645 0.0003679 ***
## cam 2 0.2176 0.108779 1.4471 0.2405446
## Residuals 92 6.9154 0.075168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = "cam")
grid <- attr(L, "grid")
grid$cam <- factor(grid$cam, levels = levels(params$cam))
rownames(L) <- levels(params$cam)
grid <- apmc(X = L, model = m0, test = "single-step", focus = "cam")
xlab <- expression("Camada do solo"~(m))
ylab <- expression("Tensão matricial da inflexão (I)")
segplot(cam ~ lwr + upr,
centers = estimate,
data = grid,
draw = FALSE, horizontal = FALSE,
xlab = xlab,
ylab = ylab,
f = 0.15, panel = panel.segplotBy,
letters = tolower(grid$cld), digits = 3)
grid
## cam estimate lwr upr cld
## 1 0 - 0,05 0.6200250 0.5418477 0.6982024 a
## 2 0,05 - 0,40 0.6479023 0.5669164 0.7288881 a
## 3 0,40 - 0,80 0.7133033 0.6335450 0.7930616 a
Análise de Componentes Principais com as Variáveis Físico Hídricas
#-----------------------------------------------------------------------
dd <- merge(da, params)
solo <- subset(dd, select = c(loc, cam, ds, Ur, Us, S, I, cad))
str(solo)
## 'data.frame': 144 obs. of 8 variables:
## $ loc: int 10 1 10 10 1 1 11 11 11 12 ...
## $ cam: Factor w/ 3 levels "0 - 0,05","0,05 - 0,40",..: 1 1 2 3 2 3 1 2 3 1 ...
## $ ds : num 1.65 1.45 1.69 1.45 1.67 ...
## $ Ur : num 0.131 0.223 0.165 0.172 0.178 ...
## $ Us : num 0.697 0.648 0.573 0.653 0.58 ...
## $ S : num -0.515 -0.341 -0.257 -0.202 -0.273 ...
## $ I : num 0.62 0.815 0.691 0.77 0.957 ...
## $ cad: num 0.3 0.227 0.222 0.274 0.217 ...
plan <- subset(dd, cam == levels(cam)[1],
select = c(loc, vol, alt, dap, prod))
str(plan)
## 'data.frame': 49 obs. of 5 variables:
## $ loc : int 10 1 11 12 13 14 15 16 17 18 ...
## $ vol : num 0.176 0.271 0.433 0.344 0.372 ...
## $ alt : num 15 12.9 19.8 21.2 19.6 ...
## $ dap : num 0.207 0.194 0.293 0.319 0.291 ...
## $ prod: num 44.4 68.2 109 86.7 93.8 ...
splom(plan)
# Já que volume é função do DAP e altura, e produção é representar em 1
# ha, então só será analisaro a produção.
#-----------------------------------------------------------------------
# Análise das variáveis de solo.
L <- split(solo, f = solo$cam)
L <- lapply(L,
FUN = function(d) {
names(d)[-c(1:2)] <- paste(names(d)[-c(1:2)],
as.integer(d$cam[1]),
sep = ".")
d[, -2]
})
L <- Reduce(f = function(x, y) { merge(x, y, by = "loc", all = TRUE) },
x = L,
accumulate = FALSE)
str(L)
## 'data.frame': 50 obs. of 19 variables:
## $ loc : int 1 2 3 4 5 6 7 8 9 10 ...
## $ ds.1 : num 1.45 1.61 1.61 1.49 1.67 ...
## $ Ur.1 : num 0.223 0.169 0.214 0.184 0.251 ...
## $ Us.1 : num 0.648 0.622 0.667 0.66 0.576 ...
## $ S.1 : num -0.341 -0.439 -0.306 -0.537 -0.159 ...
## $ I.1 : num 0.815 0.642 0.734 0.767 1.037 ...
## $ cad.1: num 0.227 0.239 0.245 0.249 0.182 ...
## $ ds.2 : num 1.67 1.65 1.71 1.62 1.68 ...
## $ Ur.2 : num 0.178 0.165 0.187 0.165 0.246 ...
## $ Us.2 : num 0.58 0.654 0.596 0.608 0.562 ...
## $ S.2 : num -0.273 -0.487 -0.219 -0.324 -0.147 ...
## $ I.2 : num 0.957 0.678 0.869 0.662 0.965 ...
## $ cad.2: num 0.217 0.257 0.226 0.238 0.178 ...
## $ ds.3 : num 1.56 1.56 1.51 NA 1.43 ...
## $ Ur.3 : num 0.188 0.179 0.253 NA 0.225 ...
## $ Us.3 : num 0.647 0.663 0.642 NA 0.628 ...
## $ S.3 : num -0.329 -0.395 -0.186 NA -0.201 ...
## $ I.3 : num 0.835 0.706 0.673 NA 0.969 ...
## $ cad.3: num 0.247 0.258 0.219 NA 0.225 ...
X <- as.matrix(L[, -1])
rownames(X) <- L$loc
X <- na.omit(X)
str(X)
## num [1:45, 1:18] 1.45 1.61 1.61 1.67 1.55 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:45] "1" "2" "3" "5" ...
## ..$ : chr [1:18] "ds.1" "Ur.1" "Us.1" "S.1" ...
## - attr(*, "na.action")=Class 'omit' Named num [1:5] 27 37 40 47 4
## .. ..- attr(*, "names")= chr [1:5] "27" "37" "40" "47" ...
pca <- princomp(x = X, cor = TRUE)
summary(pca)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 1.9429619 1.7541229 1.6387782 1.4566317
## Proportion of Variance 0.2097278 0.1709415 0.1491997 0.1178764
## Cumulative Proportion 0.2097278 0.3806694 0.5298690 0.6477455
## Comp.5 Comp.6 Comp.7 Comp.8
## Standard deviation 1.26489798 1.07400625 0.91577549 0.82418153
## Proportion of Variance 0.08888705 0.06408275 0.04659137 0.03773751
## Cumulative Proportion 0.73663251 0.80071525 0.84730663 0.88504414
## Comp.9 Comp.10 Comp.11 Comp.12
## Standard deviation 0.7895074 0.65139835 0.59404380 0.48956859
## Proportion of Variance 0.0346290 0.02357332 0.01960489 0.01331541
## Cumulative Proportion 0.9196731 0.94324646 0.96285135 0.97616676
## Comp.13 Comp.14 Comp.15
## Standard deviation 0.45683340 0.34149938 0.308434814
## Proportion of Variance 0.01159426 0.00647899 0.005285113
## Cumulative Proportion 0.98776103 0.99424002 0.999525131
## Comp.16 Comp.17 Comp.18
## Standard deviation 0.0744403128 0.0503790103 2.163896e-02
## Proportion of Variance 0.0003078533 0.0001410025 2.601358e-05
## Cumulative Proportion 0.9998329839 0.9999739864 1.000000e+00
# screeplot(pca, type = "lines", main = NULL)
## Proporção de variância acumulada.
plot(cumsum(pca$sdev^2)/sum(pca$sdev^2), type = "o",
xlab = "Componente", ylab = "Proporção de variância acumulada")
abline(h = 0.75, lty = 2)
#----------------------------------------------------------------------
# Gráficos biplot.
biplot(pca, choices = c(1, 2))
biplot(pca, choices = c(1, 3))
biplot(pca, choices = c(2, 3))
Ajuste de Modelos de Regressão para Explicar a Produção de Madeira
#-----------------------------------------------------------------------
dd <- merge(plan, L)
i <- c(5, grep(x = names(dd), pattern = "\\."))
dd <- na.omit(dd[, i])
splom(dd)
m0 <- lm(prod ~ ., data = dd)
par(mfrow = c(2, 2))
plot(m0); layout(1)
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## ds.1 1 6980.2 6980.2 8.8209 0.006329 **
## Ur.1 1 8107.7 8107.7 10.2457 0.003595 **
## Us.1 1 5976.7 5976.7 7.5527 0.010745 *
## S.1 1 49.3 49.3 0.0624 0.804772
## I.1 1 182.6 182.6 0.2307 0.635023
## cad.1 1 1664.2 1664.2 2.1030 0.158965
## ds.2 1 736.0 736.0 0.9300 0.343740
## Ur.2 1 217.3 217.3 0.2747 0.604665
## Us.2 1 444.3 444.3 0.5614 0.460416
## S.2 1 1641.6 1641.6 2.0745 0.161715
## I.2 1 3910.7 3910.7 4.9420 0.035116 *
## cad.2 1 3610.1 3610.1 4.5621 0.042272 *
## ds.3 1 262.9 262.9 0.3322 0.569300
## Ur.3 1 3065.5 3065.5 3.8739 0.059792 .
## Us.3 1 57.6 57.6 0.0728 0.789427
## S.3 1 5.9 5.9 0.0074 0.931916
## I.3 1 1825.9 1825.9 2.3074 0.140827
## cad.3 1 647.0 647.0 0.8176 0.374191
## Residuals 26 20574.5 791.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# m1 <- step(m0, k = log(nrow(dd)))
# m1 <- step(m0)
# summary(m1)
#-----------------------------------------------------------------------
# Armazenamento do solo considerando as 3 camadas.
i <- grep(x = names(dd), pattern = "cad\\.")
dd$arm <- apply(dd[, i], MARGIN = 1,
FUN = function(x) {
sum(x * c(0.05, 0.35, 0.4))
})
#-----------------------------------------------------------------------
# Ajuste de um modelo para a produção de madeira.
xyplot(prod ~ arm + ds.2 + Us.2 + Ur.2 + S.2 + I.2,
data = dd, outer = TRUE, type = c("p", "smooth"),
scales = list(x = list(relation = "free")),
xlab = "Valor da variável físico-hídrica do solo",
ylab = expression("Produção de madeira"~(m^3~ha^{-1})))
m0 <- lm(prod ~ S.2, data = dd)
par(mfrow = c(2, 2))
plot(m0); layout(1)
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## S.2 1 7869 7868.6 6.4953 0.01447 *
## Residuals 43 52091 1211.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ S.2, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -63.406 -28.400 -1.851 19.790 79.684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 115.72 14.42 8.023 4.39e-10 ***
## S.2 167.03 65.54 2.549 0.0145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34.81 on 43 degrees of freedom
## Multiple R-squared: 0.1312, Adjusted R-squared: 0.111
## F-statistic: 6.495 on 1 and 43 DF, p-value: 0.01447