source("config/setup.R") # A revista pede que as figuras sejam em Serif. opts_chunk$set(#dpi = 96, dev = "png", dev.args = list(family = "Palatino", pointsize = 18)) # width in pixels = inches * dpi * 2. # 7 * 96 * 2
#----------------------------------------------------------------------- # Carrega os pacotes necessários. rm(list = ls()) library(lattice) library(latticeExtra) library(captioner) library(reshape2) library(plyr) library(nlme) library(splines) library(multcomp) library(wzRfun) library(EACS) #----------------------------------------------------------------------- # Funções. tb_nums <- captioner(prefix = "Tabela") fg_nums <- captioner(prefix = "Figura") tb_lab <- function(label) tb_nums(label, display = "cite") fg_lab <- function(label) fg_nums(label, display = "cite") eseq <- function(x, length.out = 25, f = 0.05) { er <- extendrange(x, f = f) seq(er[1], er[2], length.out = length.out) } # Obtém intervalo de confiança de predição para modelo de classe lme. icpred <- function(model, data) { qn <- qnorm(0.975) X <- model.matrix(formula(model)[-2], data = data) V <- model$varFix se <- sqrt(diag(X %*% V %*% t(X))) fit <- X %*% fixef(model) ic <- cbind(fit, fit - qn * se, fit + qn * se) colnames(ic) <- c("fit", "lwr", "upr") return(ic) } # Sem correlação entre estimativas de efeito fixo na nlme. # http://stackoverflow.com/questions/16518438/how-to-suppress-correlation-table-in-lme assignInNamespace("print.correlation", function(x, title) return(), ns = "nlme")
devtools::load_all()
Os dados no data frame rp_eucal
são resultados de um experimento em
delineamento de blocos casualizados onde estudou-se o efeito de 4 tipos
de manejo sobre as propriedades físicas do solo sob o cultivo de
eucalipto: densidade do solo, resistência à penetração, umidade em 0 e 6
kPa. Em cada parcela, amostras de solo indeformadas foram extraídas de 9
camadas de 5 cm de espessura a partir da superfície no mesmo ponto
(medidas repetidas), sendo assim, explorados 45 cm de profundidade.
# Nomes curtos são mais fáceis de manuzear. rp <- rp_eucal rp$dumid <- rp$umid0 - rp$umid6 rp$cam <- rp$cam/100 str(rp) # Visualização do desenho experimental. xtabs(~manejo + cam, data = rp)
des <- data.frame( Nomenclatura = c("BS", "BL", "BSL", "S"), Descrição = c( "Manejo conservacionista do solo com cultivo de braquiária nas entrelinhas e adubação total (100%) no sulco de plantio.", "Manejo conservacionista do solo com cultivo de braquiária nas entrelinhas e adubação total (100%) a lanço.", "Manejo conservacionista do solo com cultivo de braquiária nas entrelinhas e adubação metade (50%) no sulco de plantio e metade (50%) a lanço (BSL = BS + BL).", "Manejo convencional do solo sem o cultivo de braquiária na entrelinhas e adubação total (100%) no sulco de plantio (testemunha)." ), Sist = c("Conserv.", "Convenc.")[c(1, 1, 1, 2)], Braq = c("Sim", "Não")[c(1, 1, 1, 2)], Sulc = c(100, 0, 50, 100), Lanç = c(0, 100, 50, 0) ) cap <- "Nomenclatura usada para os níveis de manejo de solo, descrição, sistema (Sist), uso de braquiária nas entrelinhas (Braq), porcentagem de adubação no sulco (Sulc) e a lanço (Lanç)." cap <- tb_nums("descr", cap) kable(des, caption = cap, row.names = FALSE)
TODO modelo: considerar erros correlacionados entre camadas de uma parcela. Blocos podem ser de efeito aleatório.
txt <- "Diagrama de dispersão da umidade (m$^3$ m$^{-3}$) em 0 e 6 kPa em função da profundidade no perfil do solo (m) para cada manejo do solo. A linhas suaves foram adicionadas para marcar a tendência sobre a média dos dados." cap <- fg_nums(name = "umids", cap = txt) # Como as umidades estão na mesma escala, pode-se usar o mesmo eixo para # as duas. xyplot(umid0 + umid6 ~ cam | manejo, data = rp, xlab = "Soil depth (m)", ylab = expression("Soil water content" ~ (m^3 ~ m^{-3})), type = c("p", "smooth"), auto.key = TRUE)
Pela r fg_lab("umids")
, a umidade do solo em 0 kPa parece não depender
da profundidade do solo e ser a mesma entre os manejos. Visualmente, o
valor médio seria de 0.37 g g$^{-1}$. A umidade em 6 kPa tem o mesmo
comportamento que a em 0 kPa. A diferença de umidade mantém-se
praticamente constante em um valor de 0.15 ao longo da profundidade do
solo e é a mesma entre os manejos.
Para o manejo BSL, na profundidade 45 cm, tem-se uma observação de umidade em 0 kPa e uma em 6 kPa atípicas, com valores acima de 0.45 de umidade. Essa observação pertence ao bloco II.
# Cela dos possíveis outliers. subset(rp, manejo == "BSL" & cam == 45)
txt <- "Diagrama de dispersão de cada uma das respostas em função da profundidade no perfil do solo (m) para cada manejo do solo. A linhas suaves foram adicionadas para marcar a tendência sobre a média dos dados. Pontos com as mesmas cores pertencem ao mesmo bloco." cap <- fg_nums(name = "todas", cap = txt) # A linha preta grossa é a linha média e as coloridas são as dentro de # cada bloco (individuais). xy <- xyplot(umid0 + umid6 + dumid + dens + rp ~ cam | manejo, outer = TRUE, groups = bloc, data = rp, xlab = "Soil depth (m)", ylab = NULL, scales = list(y = list(relation = "free")), type = c("p", "smooth")) + layer(panel.smoother(..., se = FALSE, lwd = 2, span = 1.5)) combineLimits(useOuterStrips(xy))
A r fg_lab("todas")
exibe cada uma das respostas, incluindo a
diferença de umidade entre 0 e 6 kPa (dumid
), em função da
profundidade para cada manejo. Visualmente, o maior efeito da
profundudade e diferença entre manejos está na densidade do solo e na
resistência à penetração, sendo o manejo BSL o mais distinto dos demais.
Conforme o delineamento considerado para realização do experimento e aquisição dos dados, o modelo proposto é
$$ \mu_{ijk} = \mu + \text{BLC}i + \text{MAN}_j \times s(\text{CAM}_k) + u{ij} + e_{ijk} $$
em que
O efeito de bloco ($BLC$) foi considerado como aleatório. Para o efeito de camada considerou-se base splines cúbico de 3 graus de liberdade.
Para o erro em nível de amostra, considerou-se inicialmente uma estrutura de correlação exponencial $$ \rho(d) = \exp\left{\frac{-d}{\lambda}\right} $$ em que $d$ é a distância entre duas camadas e $\lambda$ é o parâmetro relacionado ao alcance. Na prática, o alcançe é a distância na qual a correlação é 5% que corresponde à $d = 3\rho \approx 0.05$.
Para acomodar o efeito de parcela, foi criada a variável ue
concatenando
os rótulos de bloco e manejo. Para fazer análise de modelos mistos com a
função lme()
é recomendado criar um data frame de classe
groupedData
.
Para cada resposta, três variações do modelo acima descrito foram ajustados:
m0
: modelo em que $e_{ijk}$ é considerado independente dentro da
parcela $ij$;m1
: modelo em que $e_{ijk}$ é considerado autocorrelacionado
dentro da parcela $ij$ descrito pela função exponencial. O modelo
m0
está encaixado no modelo m1
;m2
: modelo igual ao m1
porém com efeito de camada descrito pelo
polinômio de grau 1 (reta). O modelo m2
está encaixado no modelo
m2
.O modelo foi ajustado aos dados pelo método da máxima verossimilhança (ML). Testes entre modelos encaixados foram feitos pelo teste da razão de verossimilhanças considerando um nível de significância de 5%. O resultado de ajuste é apresentado pelo gráfico com as curvas ajustadas acompanhadas das bandas de confiança.
# Ordena e cria UE para representar parcelas (unidades experimentais). rp <- rp[with(rp, order(bloc, manejo, cam)), ] rp$ue <- with(rp, interaction(bloc, manejo, drop = TRUE)) rownames(rp) <- NULL # nlevels(rp$ue) # Talvez um outlier na umidade. Obs 63. # subset(rp, umid0 > 0.5 & umid6 > 0.4) # Cria o `groupedData` para ter `augPred` do ajuste. rpg <- groupedData(rp ~ cam | ue, data = rp, order.groups = FALSE)
# Modelo que assume independência entre camadas. m0 <- lme(umid0 ~ manejo * ns(cam), data = rpg[-63, ], random = ~1 | bloc/ue, method = "ML") # Modelo que expecifica correlação exponencial. m1 <- update(m0, correlation = corExp(value = 0.05, form = ~cam)) anova(m0, m1) # Comparação visual dos ajustes. # plot(comparePred(m0, m1), layout = c(4, NA)) # Modelo reduzido em termos fixos. m2 <- update(m1, fixed = . ~ manejo * cam) anova(m1, m2)
Para a resposta umid0
, foi verificado que o modelo que assume
correlação nula entre camadas (m0
) foi rejeitado em favor do modelo
que especifica uma correlação do tipo exponencial (m1
). O modelo com
efeito linear para a camada (m2
) mostrou-se tão adequado quanto o
modelo com base splines (m1
). Neste modelo final (m2
) não foi
verificado, pela tabela de testes de Wald[^2], efeito de manejos e
camada do solo.
# Resultados pelo modelo final. anova(m2) summary(m2) # Estimativa de alcance. r <- coef(m2$modelStruct$corStruct, unconstrained = FALSE) # Um alcance de 95% (correlação 5%) está na distância de 3 * r. 3 * r
txt <- "Correlação entre camadas em função da distância (m) entre elas no perfil do solo." cap <- fg_nums(name = "corr-umid0", cap = txt) curve(exp(-x/r), from = 0, to = max(rp$cam), ylab = "Correlação entre observações", xlab = "Distância entre camadas do solo (m)") abline(v = 3 * r, h = exp(-3), lty = 2)
A estimativa de alcance[^1] foi de r 3*r
m
(r fg_lab("corr-umid0")
). O gráfico mostra o decaimento da correlação
entre camadas em função da distância entre elas segundo a função
exponencial considerada no ajuste.
txt <- "Valores observados de umidade do solo (m$^3$ m$^{-3}$) em 0 kPa em função da profundidade do solo (m) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%)." cap <- fg_nums(name = "pred-umid0", cap = txt) pred <- with(rp, expand.grid(manejo = levels(manejo), cam = eseq(cam, f = 0))) pred <- cbind(pred, icpred(model = m2, data = pred)) xyplot(umid0 ~ cam | manejo, data = rp[-63, ], ylab = expression("Soil water content at 0 kPa" ~ (m^3 ~ m^{-3})), xlab = "Soil depth (m)") + as.layer(xyplot(fit ~ cam | manejo, data = pred, type = "l", ly = pred$lwr, uy = pred$upr, cty = "bands", prepanel = prepanel.cbH, panel = panel.cbH))
Na r fg_lab("pred-umid0")
estão os valores observados e ajustados pelo
modelo. As bandas de confiança podem conter uma linha horizontal o que
uma confirmalção visual de que o modelo nulo (reta horizontal ou apenas
intercepto) não é um modelo rejeitado para descrever o comportamento do
a umidade do solo em relação à profundidade.
# Modelo que assume independência entre camadas. m0 <- lme(umid6 ~ manejo * bs(cam), data = rpg[-63, ], random = ~1 | bloc/ue, method = "ML") # Modelo que expecifica correlação exponencial. m1 <- update(m0, correlation = corExp(value = 0.02, form = ~cam)) anova(m0, m1) # Comparação visual dos ajustes. # plot(comparePred(m0, m1), layout = c(4, NA)) # Modelo reduzido em termos fixos. m2 <- update(m1, fixed = . ~ manejo * cam) anova(m1, m2)
Para a resposta umid6
, assim como ocorre para umid0
, foi verificado
que o modelo que assume correlação nula entre camadas (m0
) foi
rejeitado em favor do modelo que especifica uma correlação do tipo
exponencial (m1
). O modelo com efeito linear para a camada (m2
)
mostrou-se tão adequado quanto o modelo com base splines (m1
). Neste
modelo final (m2
) não foi verificado à 5% (mas foi à 10%), pela
tabela de testes de Wald, efeito de manejos e camada do solo.
# Resultados pelo modelo final. anova(m2) summary(m2) # Estimativa de alcance. r <- coef(m2$modelStruct$corStruct, unconstrained = FALSE) # Um alcance de 95% (correlação 5%) está na distância de 3 * r. 3 * r
txt <- "Correlação entre camadas em função da distância (m) entre elas no perfil do solo." cap <- fg_nums(name = "corr-umid6", cap = txt) curve(exp(-x/r), from = 0, to = max(rp$cam), ylab = "Correlação entre observações", xlab = "Distância entre camadas do solo (m)") abline(v = 3 * r, h = exp(-3), lty = 2)
A estimativa de alcance foi de r 3*r
m
(r fg_lab("corr-umid6")
). O gráfico mostra o decaimento da correlação
entre camadas em função da distância entre elas segundo a função
exponencial considerada no ajuste.
txt <- "Valores observados de umidade do solo (m$^3$ m$^{-3}$) em 6 kPa em função da profundidade do solo (m) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%)." cap <- fg_nums(name = "pred-umid6", cap = txt) pred <- with(rp, expand.grid(manejo = levels(manejo), cam = eseq(cam, f = 0))) pred <- cbind(pred, icpred(model = m2, data = pred)) xyplot(umid6 ~ cam | manejo, data = rp[-63, ], ylab = expression("Soil water content at 6 kPa" ~ (m^3 ~ m^{-3})), xlab = "Soil depth (m)") + as.layer(xyplot(fit ~ cam | manejo, data = pred, type = "l", ly = pred$lwr, uy = pred$upr, cty = "bands", prepanel = prepanel.cbH, panel = panel.cbH))
Na r fg_lab("pred-umid6")
estão os valores observados e ajustados pelo
modelo. Exceto para o manejo BSL, as bandas de confiança podem conter
uma linha horizontal. Para o BSL verifica-se um efeito significativo de
camada.
# Modelo que assume independência entre camadas. m0 <- lme(dumid ~ manejo * bs(cam), data = rpg, random = ~1 | bloc/ue, method = "ML") # Modelo que expecifica correlação exponencial. m1 <- update(m0, correlation = corExp(value = 0.05, form = ~cam)) anova(m0, m1) # Comparação visual dos ajustes. # plot(comparePred(m0, m1), layout = c(4, NA)) # Modelo reduzido em termos fixos. m2 <- update(m0, fixed = . ~ manejo * cam) anova(m1, m2)
Para a resposta dumid
, ao contrário de umid0
e umid6
, foi
verificado que o modelo que assume correlação nula entre camadas (m0
)
não foi rejeitado em favor do modelo que especifica uma correlação do
tipo exponencial (m1
). O modelo com efeito linear para a camada (m2
)
não se mostrou superior ao modelo com base splines (m1
). Neste
modelo final (m0
), foi detectado interação entre manejo e camada do
solo.
# Resultados pelo modelo final. anova(m0) summary(m0) L <- rbind(c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0), c(0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0), c(0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1)) rownames(L) <- levels(rp$manejo) # Teste para o efeito de profundidade dentro de cada camada. summary(glht(m0, linfct = L))
txt <- "Valores observados para a diferença entre umidade do solo em 0 e 6 kPa (m$^3$ m$^{-3}$) em função da profundidade do solo (m) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%)." cap <- fg_nums(name = "pred-dumid", cap = txt) pred <- with(rp, expand.grid(manejo = levels(manejo), cam = eseq(cam, f = 0))) pred <- cbind(pred, icpred(model = m0, data = pred)) xyplot(dumid ~ cam | manejo, data = rp, ylab = expression("Soil water content - 0 to 6 kPa" ~ (m^3 ~ m^{-3})), xlab = "Soil depth (m)") + as.layer(xyplot(fit ~ cam | manejo, data = pred, type = "l", ly = pred$lwr, uy = pred$upr, cty = "bands", prepanel = prepanel.cbH, panel = panel.cbH))
Na r fg_lab("pred-dumid")
estão os valores observados e ajustados pelo
modelo. Exceto para o manejo BSL, as bandas de confiança podem conter
uma linha horizontal. Para o BSL verifica-se um efeito significativo de
camada descrito pelo base splines.
# Modelo que assume independência entre camadas. m0 <- lme(dens ~ manejo * bs(cam), data = rpg, random = ~1 | bloc/ue, method = "ML") # Modelo que expecifica correlação exponencial. m1 <- update(m0, correlation = corExp(value = 0.05, form = ~cam)) anova(m0, m1) # Comparação visual dos ajustes. # plot(comparePred(m0, m1), layout = c(4, NA)) # Modelo reduzido em termos fixos. m2 <- update(m0, fixed = . ~ manejo * cam) anova(m1, m2)
Para a resposta dens
, assim como ocorre com dumid
, foi verificado
que o modelo que assume correlação nula entre camadas (m0
) não foi
rejeitado em favor do modelo que especifica uma correlação do tipo
exponencial (m1
). O modelo com efeito linear para a camada (m2
) não
se mostrou superior ao modelo com base splines (m1
). Neste modelo
final (m0
), foi detectado interação entre manejo e camada do solo.
# Resultados pelo modelo final. anova(m0) summary(m0) L <- rbind(c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0), c(0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0), c(0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1)) rownames(L) <- levels(rp$manejo) # Teste para o efeito de profundidade dentro de cada camada. summary(glht(m0, linfct = L))
txt <- "Valores observados para a densidade do solo (Mg m$^{-3}$) em função da profundidade do solo (m) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%)." cap <- fg_nums(name = "pred-dens", cap = txt) pred <- with(rp, expand.grid(manejo = levels(manejo), cam = eseq(cam, f = 0))) pred <- cbind(pred, icpred(model = m0, data = pred)) xyplot(dens ~ cam | manejo, data = rp, ylab = expression("Bulk density of soil" ~ (Mg ~ m^{-3})), xlab = "Soil depth (m)") + as.layer(xyplot(fit ~ cam | manejo, data = pred, type = "l", ly = pred$lwr, uy = pred$upr, cty = "bands", prepanel = prepanel.cbH, panel = panel.cbH))
Na r fg_lab("pred-dens")
estão os valores observados e ajustados pelo
modelo. Exceto para o manejo BSL, as bandas de confiança podem conter
uma linha horizontal. Para o BSL verifica-se um efeito significativo de
camada descrito pelo base splines.
# Modelo que assume independência entre camadas. m0 <- lme(rp ~ manejo * bs(cam), data = rpg, random = ~1 | bloc/ue, method = "ML") # Modelo que expecifica correlação exponencial. m1 <- update(m0, correlation = corExp(value = 0.05, form = ~cam)) anova(m0, m1) # Comparação visual dos ajustes. # plot(comparePred(m0, m1), layout = c(4, NA)) # Modelo reduzido em termos fixos. m2 <- update(m0, fixed = . ~ manejo * cam) anova(m1, m2)
Para a resposta rp
, foi verificado que o modelo que assume correlação
nula entre camadas (m0
) foi rejeitado em favor do modelo que
especifica uma correlação do tipo exponencial (m1
). O modelo com
efeito linear para a camada (m2
) não se mostrou superior ao modelo com
base splines (m1
). Neste modelo final (m0
), foi detectado
interação entre manejo e camada do solo.
# Resultados pelo modelo final. anova(m0) summary(m0) L <- rbind(c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0), c(0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0), c(0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1)) rownames(L) <- levels(rp$manejo) # Teste para o efeito de profundidade dentro de cada camada. summary(glht(m0, linfct = L))
txt <- "Valores observados para a resistência a penetração do solo (MPa) em função da profundidade do solo (m) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%)." cap <- fg_nums(name = "pred-rp", cap = txt) pred <- with(rp, expand.grid(manejo = levels(manejo), cam = eseq(cam, f = 0))) pred <- cbind(pred, icpred(model = m0, data = pred)) xyplot(rp ~ cam | manejo, data = rp, ylab = expression("Soi penetration resistance" ~ (MPa)), xlab = "Soil depth (m)") + as.layer(xyplot(fit ~ cam | manejo, data = pred, type = "l", ly = pred$lwr, uy = pred$upr, cty = "bands", prepanel = prepanel.cbH, panel = panel.cbH))
Na r fg_lab("pred-rp")
estão os valores observados e ajustados pelo
modelo. Apenas para o manejo BSL, as bandas de confiança pode conter uma
linha horizontal. Para os demais manejos verifica-se um efeito
significativo de camada descrito pelo base splines.
# Resíduos dentro de cada UE. r <- cbind(res = unlist(residuals(m0)), rpg[, c("ue", "cam")]) r <- dcast(data = r, formula = ue ~ cam, value.var = "res") # splom(r[, -1], type = c("p", "r")) r <- cor(r[, -1]) colnames(r) <- rownames(r) <- unique(sort(rp$cam)) levelplot(r, at = seq(-1, 1, by = 0.1), main = list("Correlação entre camadas", font = 1), xlab = "Camada (m)", ylab = "Camada (m)") # Predomínio e correlação positiva entre camadas adjacentes. A # correlação tende a se tornar negativa com o aumento da separação entre # camadas. # Autocorrelação. plot(ACF(m0)) # Extraindo as estimativas e erros padrões. L <- list(m0 = m0, m1 = m1) L <- lapply(L, FUN = function(m) { est <- cbind(data.frame(coef = names(coef(m))), summary(m)$tTable[, 1:2]) return(est) }) L <- ldply(L, .id = "model") # Gráfico. xyplot(m0 ~ m1, aspect = "iso", data = dcast(L, coef ~ model, value.var = "Std.Error")) + layer(panel.abline(a = 0, b = 1)) # O modelo com correlação nos erros apresenta erros padrões # maiores. Isso está ligado ao fato de que por haver correlação entre as # camadas existe informação redundante e isso não diminui a incerteza.
cat(format(Sys.time(), format = "Atualizado em %d de %B de %Y.\n\n")) sessionInfo()
[^1]: Alcance é definida como sendo a distância para uma correlação de 5%. [^2]: Em modelos mistos, o teste para os efeitos fixos é o teste F de Wald e não o teste F de um quadro de ANOVA. A interpretação permanece a mesma.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.