# https://github.com/walmes/wzRfun # devtools::install_github("walmes/wzRfun") library(wzRfun) pkg <- c("lattice", "latticeExtra", "nlme", "reshape", "plyr", "doBy", "multcomp") sapply(pkg, library, character.only = TRUE, logical.return = TRUE)
library(RDASC)
source("config/setup.R")
#----------------------------------------------------------------------- # Lendo arquivos de dados. data(pineapple_swrc) str(pineapple_swrc) # Converte variáveis para fator. pineapple_swrc <- transform(pineapple_swrc, gesso = factor(gesso), ue = interaction(cober, gesso, prof, varie, bloc, drop = TRUE, sep = ":")) # Passa a tensão 0 para 0.01 para não ter problemas com o log(). pineapple_swrc$tens[pineapple_swrc$tens == 0] <- 0.01 pineapple_swrc$ltens <- log10(pineapple_swrc$tens) # Mostra a estrutura do objeto. str(pineapple_swrc)
#----------------------------------------------------------------------- # Análise exploratória. # Identifica as unidades com a combinação tripla dos fatores. pineapple_swrc$gescobpro <- with(pineapple_swrc, interaction(as.integer(gesso), as.integer(cober), as.integer(prof))) p1 <- xyplot(umid ~ tens | gescobpro + varie, data = pineapple_swrc, groups = bloc, 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(pineapple_swrc$ue) #----------------------------------------------------------------------- # Identificar os pontos atípicos. # i <- 1:nrow(pineapple_swrc) # 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) pineapple_swrc <- pineapple_swrc[-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 = pineapple_swrc, start = start) summary(cra.fit)
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 = pineapple_swrc) with(start, curve(Ur + (Us - Ur)/(1 + exp(n * (alp + x)))^(1 - 1/n), add = TRUE)) n00 <- nls(model, data = pineapple_swrc, start = start) coef(n00) #----------------------------------------------------------------------- # Ajudar para cada unidade experimental. # nlevels(pineapple_swrc$ue) db <- groupedData(umid ~ ltens | ue, data = pineapple_swrc, order.groups = FALSE) str(db) 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", "bloc") params <- cbind(equallevels(as.data.frame(aux), pineapple_swrc), c0) rownames(params) <- NULL str(params) #----------------------------------------------------------------------- # Manova. m0 <- aov(cbind(Ur, Us, n, alp) ~ cober * gesso * prof * varie, data = params) anova(m0)
#----------------------------------------------------------------------- # 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) # summary(im) del <- im$is.inf[, "dffit"] m0 <- update(m0, data = params[!del, ]) anova(m0) #----------------------------------------------------------------------- # Desdobramento. # Estimativas para as profundidades. Xm <- LE_matrix(m0, effect = "prof") rownames(Xm) <- levels(db$prof) g1 <- apmc(X = Xm, model = m0, focus = "prof") g1 # Desdobrar gesso dentro de cobertura. LSmeans(m0, effect = c("gesso", "cober")) Xm <- LE_matrix(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 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) # summary(im) del <- im$is.inf[, "dffit"] m0 <- update(m0, data = params[!del,]) anova(m0) #----------------------------------------------------------------------- # Desdobramento. # Desdobrar prof dentro de varie. LSmeans(m0, effect = c("prof", "varie")) Xm <- LE_matrix(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 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) # summary(im) del <- im$is.inf[, "dffit"] m0 <- update(m0, data = params[!del,]) anova(m0) # 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")) LSmeans(m0, effect = c("gesso", "cober")) LSmeans(m0, effect = c("cober", "varie"))
#----------------------------------------------------------------------- # 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) # summary(im) del <- im$is.inf[, "dffit"] m0 <- update(m0, data = params[!del,]) anova(m0) # Desdobramento não realizado por argumentos iguais aos do parâmetro # alpha. LSmeans(m0, effect = c("gesso", "cober"))
#----------------------------------------------------------------------- # 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) # summary(im) del <- im$is.inf[,"dffit"] m0 <- update(m0, data = params[!del, ]) anova(m0) # Não existe efeito dos fatores experimentais sobre o S. LSmeans(m0)
#----------------------------------------------------------------------- # 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) # summary(im) del <- im$is.inf[, "dffit"] m0 <- update(m0, data = params[!del, ]) anova(m0) #----------------------------------------------------------------------- # Desdobramentos. # Desdobrar prof dentro de varie. LSmeans(m0, effect = c("prof", "varie")) Xm <- LE_matrix(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 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")) Xm <- LE_matrix(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 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) # summary(im) del <- im$is.inf[, "dffit"] m0 <- update(m0, data = params[!del,]) anova(m0) #-------------------------------------------- # Desdobramentos. # Desdobrar prof dentro de varie. LSmeans(m0, effect = c("prof", "varie")) Xm <- LE_matrix(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 g1 <- ldply(g1) names(g1)[1] <- "varie" g1 <- equallevels(g1, db) # Desdobrar gesso dentro de cobertura. LSmeans(m0, effect = c("gesso", "cober")) Xm <- LE_matrix(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 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)
#----------------------------------------------------------------------- # Modelo de efeito aleatório. n1 <- nlme(model, data = pineapple_swrc, fixed = Ur + Us + n + alp ~ 1, random = Us ~ 1 | ue, start = coef(cra.fit)) summary(n1) #-------------------------------------------- dput(round(coef(cra.fit), 4)) n2 <- nlme(model, data = pineapple_swrc, fixed = Ur + Us + n + alp ~ prof, random = Us ~ 1 | ue, start = c( 0.0925, 0, 0.3672, 0, 4.0623, 0, -0.6803, 0)) summary(n2) anova(n2, n1)
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"), "----------------------------------------", sep = "\n") sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.