# https://github.com/walmes/wzRfun # devtools::install_github("walmes/wzRfun") # .libPaths("/usr/lib/R/site-library") library(lattice) library(latticeExtra) library(plyr) library(rootSolve) # gradient(). library(wzRfun) # panel.cbH().
library(RDASC)
source("config/setup.R")
O experimento consistiu da observação da severidade da mancha foliar de Glomrella em ramos marcados de macieiras em dois pomares por 11 semanas. Em cada pomar, 30 plantas ao acaso tiveram um ramo marcado com 10 folhas. Aproximadamente a cada 7 dias, os ramos eram observados para a determinação da severidade de Glomerella em cada uma das folhas. Um total de 11 avaliações dos ramos foi feito produzindo 2 pomares $\times$ 30 ramos $\times$ 10 folhas $\times$ 11 avaliações $=$ 6600 observações de severidade.
As avaliações foram feitas nas mesmas datas em dois pomares indepedentes. Cada folha foi observada repetidamente nas 11 avaliações, exceto quando a folha caia do ramo. Sendo assim, um número de menor de folhas por ramo permanecia com o passar do tempo. Essas observações perdidas (missings) dificilmente foram perdidas ao acaso (missing at random), haja visto que o progresso da doença sobre as folhas é um fator que provoca a queda.
# Estrutura dos dados. str(leaf_spot) # Tabela de frequencia. ftable(xtabs(~pomar + dia, data = leaf_spot)) # Tabela de frequência de folhas presas ao ramo (sem folhas perdidas). ftable(xtabs(~pomar + dia, data = na.omit(leaf_spot))) # Acesse a documentação para mais detalhes. # help(leaf_spot, help_type = "html") # Convertendo variáveis para fator. leaf_spot <- within(leaf_spot, { pomar <- factor(pomar, labels = c("I", "II")) ramo <- factor(ramo) folha <- interaction(ramo, folha, drop = TRUE) }) xyplot(sever ~ dia | ramo, groups = folha, data = subset(leaf_spot, pomar == "I"), type = "o", xlab = "Dia de avaliação", ylab = "Severidade da mancha foliar (%)", main = "Pomar I", as.table = TRUE) xyplot(sever ~ dia | ramo, groups = folha, data = subset(leaf_spot, pomar == "II"), type = "o", xlab = "Dia de avaliação", ylab = "Severidade da mancha foliar (%)", main = "Pomar II", as.table = TRUE)
Pelos diagramas de dipersão, verifica-se que existe tanto variabilidade entre folhas de um mesmo ramo quanto entre ramos. O número de observações de cada folha também varia com o ramos. No pomar I, o ramo 13 teve poucas observações ao passo que o ramo 6 teve praticamente todas. Isso sugere que a forma como a doença se manifesta nos ramos depende de características locais não registradas, como a nutrição da planta, as condições de solo, a exposição do ramo ao sol, etc.
da <- subset(leaf_spot, pomar == "I") # Calibrando o chute inicial. start <- list(A = 80, I = 80, S = 20) xyplot(sever ~ dia, data = da) + layer(panel.curve(A/(1 + exp(-(x - I)/S)), col = 2), data = start) n0 <- nls(sever ~ A/(1 + exp(-(dia - I)/S)), data = da, start = start) # Diagnóstico m0 <- as.lm(n0) par(mfrow = c(2, 2)) plot(m0) layout(1) # Estimativas. summary(n0) # Taxa no ponto de inflexão. coef(n0)["A"]/(4 * coef(n0)["S"]) # Resultado do ajuste. xyplot(sever ~ dia, data = da) + layer(panel.curve(A/(1 + exp(-(x - I)/S)), col = 2), data = as.list(coef(n0)))
library(nlme) da <- da[complete.cases(da), ] da <- groupedData(sever ~ dia | ramo/folha, data = da, order.groups = FALSE) n1 <- nlme(sever ~ A/(1 + exp(-(dia - I)/S)), fixed = A + I + S ~ 1, random = I + S ~ 1 | ramo/folha, data = da, start = coef(n0)) logLik(n1) # n2 <- nlme(sever ~ A/(1 + exp(-(dia - I)/S)), # fixed = A + I + S ~ 1, # random = A + S ~ 1 | folha, # data = da, # start = coef(n0), # control = list(maxIter = 100)) # logLik(n2) # I + S ~ 1 | ramo/folha 'log Lik.' -4465.499 (df=10) # I + S ~ 1 | folha 'log Lik.' -4516.702 (df=7) # I ~ 1 | folha 'log Lik.' -4909.711 (df=5) # A ~ 1 | folha 'log Lik.' -4968.641 (df=5) # S ~ 1 | folha NA # A + I ~ 1 | folha NA # A + S ~ 1 | folha NA # Estimativas. summary(n1) # Taxa no ponto de inflexão. fixef(n1)["A"]/(4 * fixef(n1)["S"]) # Resultado do ajuste. plot(augPred(n1, level = 0), as.table = TRUE) # Diagnóstico. # r <- residuals(n1) # f <- fitted(n1) # xyplot(r ~ f) # qqmath(r)
# Estimates and standard error. summary(n0)$coeff summary(n1)$tTable # Confidence intervals. ci0 <- cbind(confint.default(n0), coef(n0)) ci1 <- intervals(n1)$fixed ci1 <- ci1[, c(2, 1, 3)] ci0 <- ci0[, c(3, 1, 2)] colnames(ci0) <- colnames(ci1) <- c("est", "lwr", "upr") ci <- as.data.frame(rbind(ci0, ci1)) ci$par <- factor(rownames(ci), levels = c("A", "I", "S")) rownames(ci) <- NULL ci$model <- gl(2, 3, labels = c("nls", "nlme")) ci segplot(model ~ lwr + upr | par, data = ci, centers = est, draw = FALSE, scales = list(x = "free"), layout = c(NA, 1), ylab = "Modelo", xlab = "Estimativa com IC de 95%") #----------------------------------------------------------------------- # Random effects. # a <- ranef(n1) # str(a) # # qqmath(a$ramo$I) # qqmath(a$ramo$S) # splom(a$ramo) # # qqmath(a$folha$I) # qqmath(a$folha$S) # splom(a$folha) #----------------------------------------------------------------------- # Predição. # Domínio para a predição. pred <- expand.grid(dia = 0:85) # Valores preditos. pred$y0 <- predict(n0, newdata = pred) pred$y1 <- predict(n1, newdata = pred, level = 0) # Modelo escrito como função dos parâmetros (theta). f <- function(theta, xx) { with(as.list(theta), A/(1 + exp(-(xx - I)/S))) } # Matriz com as derivadas parciais de theta no mle de theta. F0 <- gradient(f, x = coef(n0), xx = pred$dia) F1 <- gradient(f, x = fixef(n1), xx = pred$dia) # Fatoração da matriz de covariância de theta. U0 <- chol(vcov(n0)) U1 <- chol(vcov(n1)) pred$se0 <- sqrt(apply(F0 %*% t(U0), 1, function(x) sum(x^2))) pred$se1 <- sqrt(apply(F1 %*% t(U1), 1, function(x) sum(x^2))) zval <- qnorm(p = c(lwr = 0.025, fit = 0.5, upr = 0.975)) me <- outer(pred$se0, zval, "*") b <- sweep(me, 1, pred$y0, "+") colnames(b) <- paste(colnames(b), "0", sep = "") pred <- cbind(pred, b) me <- outer(pred$se1, zval, "*") b <- sweep(me, 1, pred$y1, "+") colnames(b) <- paste(colnames(b), "1", sep = "") pred <- cbind(pred, b) #----------------------------------------------------------------------- # Predição para o nível de folha. predue <- unique(subset(da, select = c(ramo, folha))) dia <- seq(0, 85, by = 2) predue <- predue[rep(1:nrow(predue), each = length(dia)), ] predue$dia <- dia str(predue) a <- predict(n1, newdata = predue, level = 2) predue$y <- unlist(a) #----------------------------------------------------------------------- xyplot(sever ~ dia, data = da, jitter.x = TRUE, pch = 19, ylab = "Severidade da mancha foliar (%)", xlab = "Dia da avaliação") + as.layer(xyplot(y ~ dia, data = predue, col = "gray50", type = "l", groups = folha), under = TRUE) + as.layer(xyplot(y0 ~ dia, data = pred, type = "l", lty = 2, lwd = 2, prepanel = prepanel.cbH, cty = "bands", ly = pred$lwr0, uy = pred$upr0, fill = "red", alpha = 0.6, panel = panel.cbH)) + as.layer(xyplot(y1 ~ dia, data = pred, type = "l", lty = 1, lwd = 2, prepanel = prepanel.cbH, cty = "bands", ly = pred$lwr1, uy = pred$upr1, fill = "blue", alpha = 0.6, panel = panel.cbH))
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.