#----------------------------------------------------------------------- # Carregando pacotes e funções necessárias. # https://github.com/walmes/wzRfun # devtools::install_github("walmes/wzRfun") # devtools::load_all("~/repos/wzRfun") library(wzRfun) library(lattice) library(latticeExtra) library(plyr) library(doBy) library(multcomp)
library(RDASC)
source("config/setup.R")
Estes dados são resultados de um experimento fatorial triplo, conduzido em laboratório, que estudou o efeito de 7 inseticidas no processo de pré parasitismo de duas espécies de Trichogramma em duas espécies de lagastas da soja (hospedeiros). Várias variáveis resposta foram registradas com a finalidade de descrever o efeito dos inseticidas para as espécies de parasitóide e hospedeiro, considerando por exemplo, o tempo de vida da fêmea ovopositora, o tempo para eclosão dos ovos, a razão sexual observada e a taxa de emergência dos ovos parasitados.
#----------------------------------------------------------------------- # Estrutura dos dados. str(egg_parasitoid) # levels(egg_parasitoid$inset) # Português Inglês # Clorpirifós Chlorpyrifos # Deltametrina Deltamethrin # Espinetoram Spinetoram # Flubendiamida Flubendiamide # Indoxacarbe Indoxacarb # Novalurom Novaluron # Testemunha Control l <- c("Chlorpyrifos", "Deltamethrin", "Spinetoram", "Flubendiamide", "Indoxacarb", "Novaluron", "Control") # Níveis dos fatores experimentais. summary(egg_parasitoid[, 1:3]) # Usando nomes curtos para os níveis. egg <- egg_parasitoid levels(egg$inset) <- substr(l, 0, 5) levels(egg$paras) <- c("Atopo", "Preti") levels(egg$hosp) <- c("Anti", "Chry") # Letra maiúscula para representar os fatores estudados. names(egg)[1:3] <- c("I", "P", "H") # Tabela de frequencia planejada do experimento (7 x 2 x 2 com 20 rep.). ftable(xtabs(~I + P + H, data = egg)) # Tabela de frequência só para casos completos. ftable(xtabs(~I + P + H, data = na.omit(egg)))
kable( aggregate(cbind(paras = opar/otot, emerg = oeme/otot) ~ H + P + I, data = egg, FUN = mean), digits = 3)
A sobreviência da fêmea, 24 horas após a liberação no tubo de ensaio
para parasitar os ovos, é representada por uma variável (mort
)
dicotômica onde 1 indica que a fêmea sobreviveu e 0 que não não
sobreviveu. A variável vivo
é o oposto da variável mort
.
# Desfechos de vivo: 1 = sobreviveu, 0 = não sobreviveu. egg$vivo <- 1 - egg$mort ftable(xtabs(vivo ~ I + P + H, data = egg)) # Modelo saturado. m0 <- glm(vivo ~ (I + P + H)^3, data = subset(egg), family = quasibinomial) anova(m0, test = "F") # Modelo reduzido. m1 <- update(m0, . ~ I * (P + H)) anova(m1, test = "F") anova(m1, m0, test = "F") #----------------------------------------------------------------------- # Comparações múltiplas. lsm <- LE_matrix(m1, effect = c("I", "P", "H")) grid <- equallevels(attr(lsm, "grid"), egg) comp <- vector(mode = "list", length = 2) # Hospedeiros dentro de inseticida x parasitóide. L <- by(lsm, INDICES = with(grid, interaction(I, P)), FUN = as.matrix) L <- lapply(L, "rownames<-", levels(egg$H)) cmp <- lapply(L, apmc, model = m1, focus = "H", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- names(grid)[1:2] names(pred)[ncol(pred)] <- "cldH" comp[[1]] <- pred # Inseticidas dentro de parasitóide e hospedeiro. L <- by(lsm, INDICES = with(grid, interaction(H, P)), FUN = as.matrix) L <- lapply(L, "rownames<-", levels(egg$I)) cmp <- lapply(L, apmc, model = m1, focus = "I", test = "fdr", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) # names(pred)[1:2] <- names(grid)[1:2] names(pred)[1:2] <- c("H", "P") names(pred)[ncol(pred)] <- "cldI" pred[, ncol(pred)] <- toupper(pred[, ncol(pred)]) comp[[2]] <- pred str(comp) pred <- merge(comp[[1]], comp[[2]], by = intersect(names(comp[[1]]), names(comp[[2]]))) #----------------------------------------------------------------------- # Passa para a escala de probabilidade. i <- c("fit", "lwr", "upr") pred[, i] <- sapply(pred[, i], m1$family$linkinv) # Ordena da tabela. pred <- pred[with(pred, order(P, I, H)), ] # Intervalos de confiança do tamanho do suporte terão apenas o ponto # representado. i <- pred$upr - pred$lwr > 0.99 if (any(i)) { pred[i, ]$lwr <- pred[i, ]$lwr <- NA } # Reordena os níveis pela probalidade de sobreviência. pred$I <- reorder(pred$I, pred$fit) # Legenda. key <- list(points = list(pch = c(1, 19)), text = list(levels(egg_parasitoid$hosp), font = 3), title = "Hosts", cex.title = 1.1) pred$cld <- with(pred, paste(cldI, cldH, sep = ""))
cap <- "Estimated probability of surviving at 24h for each inseticide on two parasiods and two hosts. Segment is a confidence interval for the probability of surviving. Parasitoids estimates followed by the same lower letters in a insetice and host combination are not different at 5%. Inseticides estimates followed by the same lower letters in a parasitoid and host combination are not different at 5%." cap <- fgn_("surv", cap) pred$vjust <- -0.5 pred$vjust[pred$cld == "ABa"] <- 1.5 # Gráfico de segmentos. segplot(I ~ lwr + upr | P, centers = fit, data = pred, xlab = "Insecticides", ylab = "Probability of surviving a 24h period", draw = FALSE, horizontal = FALSE, groups = H, key = key, strip = strip.custom( factor.levels = levels(egg_parasitoid$paras), par.strip.text = list(font = 3)), gap = 0.15, cld = pred$cld, panel = panel.groups.segplot, pch = key$points$pch[as.integer(pred$H)]) + layer({ a <- cld[which.max(nchar(cld))] l <- cld[subscripts] v <- pred$vjust[subscripts] x <- as.integer(z)[subscripts] + centfac(groups[subscripts], gap) y <- centers[subscripts] # Usa símbolo unicode: # http://www.alanwood.net/unicode/geometric_shapes.html grid.text("\u25AE", x = unit(x, "native"), y = unit(y, "native"), vjust = v, gp = gpar(col = "white") ) grid.text(l, x = unit(x, "native"), y = unit(y, "native"), vjust = v, gp = gpar(col = "black", fontsize = 10)) })
#----------------------------------------------------------------------- # Análise exploratória. useOuterStrips(xyplot(opar/otot ~ I | H + P, data = egg, jitter.x = TRUE, type = c("p", "a"))) #----------------------------------------------------------------------- # Ajuste do modelo. m0 <- glm(cbind(opar, otot - opar) ~ I * P * H, data = egg, family = quasibinomial) par(mfrow = c(2, 2)) plot(m0) layout(1) anova(m0, test = "F") # summary(m0) #----------------------------------------------------------------------- # Comparações múltiplas. lsm <- LE_matrix(m0, effect = c("I", "P", "H")) grid <- equallevels(attr(lsm, "grid"), egg) comp <- vector(mode = "list", length = 2) # Hospedeiros dentro de inseticida x parasitóide. L <- by(lsm, INDICES = with(grid, interaction(I, P)), FUN = as.matrix) L <- lapply(L, "rownames<-", levels(egg$H)) cmp <- lapply(L, apmc, model = m0, focus = "H", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- names(grid)[1:2] names(pred)[ncol(pred)] <- "cldH" comp[[1]] <- pred # Inseticidas dentro de parasitóide e hospedeiro. L <- by(lsm, INDICES = with(grid, interaction(H, P)), FUN = as.matrix) L <- lapply(L, "rownames<-", levels(egg$I)) cmp <- lapply(L, apmc, model = m0, focus = "I", test = "fdr", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) # names(pred)[1:2] <- names(grid)[1:2] names(pred)[1:2] <- c("H", "P") names(pred)[ncol(pred)] <- "cldI" pred[, ncol(pred)] <- toupper(pred[, ncol(pred)]) comp[[2]] <- pred str(comp) pred <- merge(comp[[1]], comp[[2]], by = intersect(names(comp[[1]]), names(comp[[2]]))) #----------------------------------------------------------------------- # Passa para a escala de probabilidade. i <- c("fit", "lwr", "upr") pred[, i] <- sapply(pred[, i], m1$family$linkinv) # Ordena da tabela. pred <- pred[with(pred, order(P, I, H)), ] # Intervalos de confiança do tamanho do suporte terão apenas o ponto # representado. i <- pred$upr - pred$lwr > 0.99 if (any(i)) { pred[i, ]$lwr <- pred[i, ]$lwr <- NA } # Reordena os níveis pela probalidade de sobreviência. pred$I <- reorder(pred$I, pred$fit) # Legenda. key <- list(points = list(pch = c(1, 19)), text = list(levels(egg_parasitoid$hosp), font = 3), title = "Hosts", cex.title = 1.1) pred$cld <- with(pred, paste(cldI, cldH, sep = "")) #----------------------------------------------------------------------- ab <- aggregate(cbind(paras = opar/otot) ~ I + H + P, data = egg, FUN = mean) kable(merge(pred, ab, by = intersect(names(pred), names(ab)))[, -c(5:9)]) #-----------------------------------------------------------------------
cap <- "Estimated proportion of of parasitated eggs for each inseticide on two parasiods and two hosts. Segment is a confidence interval for the probability of surviving. Parasitoids estimates followed by the same lower letters in a insetice and host combination are not different at 5%. Inseticides estimates followed by the same lower letters in a parasitoid and host combination are not different at 5%." cap <- fgn_("opar", cap) pred$vjust <- -0.5 pred$vjust[pred$cld == "ABCDa"] <- 1.5 # Gráfico de segmentos. segplot(I ~ lwr + upr | P, centers = fit, data = pred, xlab = "Insecticides", ylab = "Proportion of parasited eggs", draw = FALSE, horizontal = FALSE, groups = H, key = key, strip = strip.custom( factor.levels = levels(egg_parasitoid$paras), par.strip.text = list(font = 3)), gap = 0.15, cld = pred$cld, panel = panel.groups.segplot, pch = key$points$pch[as.integer(pred$H)]) + layer({ a <- cld[which.max(nchar(cld))] l <- cld[subscripts] v <- pred$vjust[subscripts] x <- as.integer(z)[subscripts] + centfac(groups[subscripts], gap) y <- centers[subscripts] # Usa símbolo unicode: # http://www.alanwood.net/unicode/geometric_shapes.html grid.text("\u25AE", x = unit(x, "native"), y = unit(y, "native"), vjust = v, gp = gpar(col = "white") ) grid.text(l, x = unit(x, "native"), y = unit(y, "native"), vjust = v, gp = gpar(col = "black", fontsize = 10)) })
#----------------------------------------------------------------------- # Análise exploratória. ftable(xtabs(!is.na(oeme) ~ I + P + H, data = egg)) # ATTENTION: Com os 7 inseticidas dá cela perdida. useOuterStrips(xyplot(oeme/opar ~ I | H + P, data = egg[!is.na(egg$oeme), ], jitter.x = TRUE, type = c("p", "a"))) #----------------------------------------------------------------------- # Ajuste do modelo. m0 <- glm(cbind(oeme, opar - oeme) ~ I * P * H, data = egg, family = quasibinomial) par(mfrow = c(2, 2)) plot(m0) layout(1) anova(m0, test = "F") # Modelo declarado com a matrix do modelo, apenas os efeitos estimáveis. X <- model.matrix(formula(m0)[-2], data = egg) b <- coef(m0) X <- X[, !is.na(b)] m0 <- update(m0, . ~ 0 + X) #----------------------------------------------------------------------- # Comparações entre hospedeiros dentro de inseticida x parasitóide. comp <- vector(mode = "list", length = 2) # Declarar um modelo não deficiente aqui apenas para pegar a matriz. lsm <- LE_matrix(lm(mort ~ I * P * H, data = egg), effect = c("I", "P", "H")) grid <- equallevels(attr(lsm, "grid"), egg) # Celas que serão mantidas pois são estimaveis. keep <- xtabs(!is.na(oeme) ~ interaction(I, P, H), data = egg) > 0 i <- with(grid, interaction(I, P, H) %in% names(keep[keep])) grid <- grid[i, ] lsm <- lsm[i, ] # Deixa apenas as colunas de efeitos estimados. lsm <- lsm[, !is.na(b)] # Fazer as comparações apenas onde é possível. # Hospedeiros dentro de inseticida x parasitóide. rownames(lsm) <- grid$H L <- by(lsm, INDICES = with(grid, interaction(I, P)), FUN = as.matrix) i <- sapply(L, is.null) L <- L[!i] cmp <- lapply(L, apmc, model = m0, focus = "H", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- c("I", "P") names(pred)[ncol(pred)] <- "cldH" comp[[1]] <- pred # Fazer as comparações apenas onde é possível. # Inseticidas dentro de parasitóide e hospedeiro. rownames(lsm) <- grid$I L <- by(lsm, INDICES = with(grid, interaction(H, P)), FUN = as.matrix) cmp <- lapply(L, apmc, model = m0, focus = "I", test = "fdr", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- c("H", "P") names(pred)[ncol(pred)] <- "cldI" pred[, ncol(pred)] <- toupper(pred[, ncol(pred)]) comp[[2]] <- pred str(comp) pred <- merge(comp[[1]], comp[[2]], by = intersect(names(comp[[1]]), names(comp[[2]]))) pred$cld <- with(pred, paste(cldI, cldH, sep = "")) #----------------------------------------------------------------------- # Passa para a escala de probabilidade. i <- c("fit", "lwr", "upr") pred[, i] <- sapply(pred[, i], m0$family$linkinv) # Ordena da tabela. pred <- pred[with(pred, order(P, I, H)), ] # Intervalos de confiança do tamanho do suporte terão apenas o ponto # representado. i <- pred$upr - pred$lwr > 0.99 if (any(i)) { pred[i, ]$lwr <- pred[i, ]$lwr <- NA } # Reordena os níveis pela probalidade de sobreviência. pred$I <- reorder(pred$I, pred$fit) ftable(xtabs(~I + P + H, data = pred)) # Legenda. key <- list(points = list(pch = c(1, 19)), text = list(levels(egg_parasitoid$hosp), font = 3), title = "Hosts", cex.title = 1.1) #----------------------------------------------------------------------- ab <- aggregate(cbind(emerg = oeme/opar) ~ I + H + P, data = egg, FUN = mean) kable(merge(pred, ab, by = intersect(names(pred), names(ab)))[, -c(5:9)])
cap <- "Estimated proportion of egg emergency for each inseticide on two parasiods and two hosts. Segment is a confidence interval for the probability of surviving. Parasitoids estimates followed by the same lower letters in a insetice and host combination are not different at 5%. Inseticides estimates followed by the same lower letters in a parasitoid and host combination are not different at 5%." cap <- fgn_("oeme", cap) # Gráfico de segmentos. segplot(I ~ lwr + upr | P, centers = fit, data = pred, xlab = "Insecticides", ylab = "Probability of egg emergency", draw = FALSE, horizontal = FALSE, groups = H, key = key, strip = strip.custom( factor.levels = levels(egg_parasitoid$paras), par.strip.text = list(font = 3)), gap = 0.15, cld = pred$cld, panel = panel.groups.segplot, pch = key$points$pch[as.integer(pred$H)]) + layer({ a <- cld[which.max(nchar(cld))] l <- cld[subscripts] x <- as.integer(z)[subscripts] + centfac(groups[subscripts], gap) y <- centers[subscripts] # Usa símbolo unicode: # http://www.alanwood.net/unicode/geometric_shapes.html grid.text("\u25AE", x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "white") ) grid.text(l, x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "black", fontsize = 10)) })
Essa variável apresenta pouca variação para algumas celas, portanto não informação suficiente para ser analisada.
#----------------------------------------------------------------------- # Análise exploratória. useOuterStrips(xyplot(femea + macho ~ I | H + P, data = egg[!is.na(egg$macho + egg$femea), ], jitter.x = TRUE, auto.key = TRUE, type = c("p", "a"))) #----------------------------------------------------------------------- # Ajuste do modelo. m0 <- glm(cbind(femea, macho) ~ I * P * H, data = egg, family = quasibinomial) par(mfrow = c(2, 2)) plot(m0) layout(1) anova(m0, test = "F") # summary(m0) m1 <- glm(cbind(femea, macho) ~ P * H, data = egg, family = quasibinomial) par(mfrow = c(2, 2)) plot(m1) layout(1) anova(m0, m1, test = "F") anova(m1, test = "F") # summary(m1) #----------------------------------------------------------------------- # Comparações entre hospedeiros dentro de inseticida x parasitóide. lsm <- LE_matrix(m1, effect = c("P", "H")) grid <- equallevels(attr(lsm, "grid"), egg) L <- by(lsm, INDICES = with(grid, P), FUN = as.matrix) L <- lapply(L, "rownames<-", levels(egg$H)) cmp <- lapply(L, apmc, model = m1, focus = "H", cld2 = TRUE) pred <- ldply(cmp) # cmp <- ldply(strsplit(pred$.id, "\\.")) # pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1] <- names(grid)[1] pred <- equallevels(pred, egg) # Passa para a escala de probabilidade. i <- c("fit", "lwr", "upr") pred[, i] <- sapply(pred[, i], m0$family$linkinv) # Ordena da tabela. pred <- pred[with(pred, order(P, H)), ] # Intervalos de confiança do tamanho do suporte terão apenas o ponto # representado. i <- pred$upr - pred$lwr > 0.99 if (any(i)) { pred[i, ]$lwr <- pred[i, ]$lwr <- NA } # Legenda. key <- list(points = list(pch = c(1, 19)), text = list(levels(egg_parasitoid$hosp), font = 3), title = "Hosts", cex.title = 1.1)
cap <- "Estimated proportion of female parasitoids for each inseticide on two parasiods and two hosts. Segment is a confidence interval for the probability of surviving. Parasitoids estimates followed by the same lower letters in a insetice and host combination are not different at 5%. Inseticides estimates followed by the same lower letters in a parasitoid and host combination are not different at 5%." cap <- fgn_("sexratio", cap) # Gráfico de segmentos. segplot(P ~ lwr + upr, centers = fit, data = pred, xlab = "Parasitoids species", ylab = "Proportion of female parasitoids", draw = FALSE, horizontal = FALSE, groups = H, key = key, gap = 0.15, cld = pred$cld, panel = panel.groups.segplot, pch = key$points$pch[as.integer(pred$H)]) + layer({ a <- cld[which.max(nchar(cld))] l <- cld[subscripts] x <- as.integer(z)[subscripts] + centfac(groups[subscripts], gap) y <- centers[subscripts] # Usa símbolo unicode: # http://www.alanwood.net/unicode/geometric_shapes.html grid.text("\u25AE", x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "white") ) grid.text(l, x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "black", fontsize = 10)) })
#----------------------------------------------------------------------- # Análise exploratória. egg$pem <- with(egg, femea + macho) useOuterStrips(xyplot(pem + pne ~ I | H + P, data = egg[!is.na(egg$pem + egg$pne), ], jitter.x = TRUE, auto.key = TRUE, type = c("p", "a"))) #----------------------------------------------------------------------- # Ajuste do modelo. m0 <- glm(cbind(pem, pne) ~ I * P * H, data = egg, family = quasibinomial) par(mfrow = c(2, 2)) plot(m0) layout(1) anova(m0, test = "F") # summary(m0) # Modelo declarado com a matrix do modelo, apenas os efeitos estimáveis. X <- model.matrix(formula(m0)[-2], data = egg) b <- coef(m0) X <- X[, !is.na(b)] m0 <- update(m0, . ~ 0 + X) #----------------------------------------------------------------------- # Comparações entre hospedeiros dentro de inseticida x parasitóide. comp <- vector(mode = "list", length = 2) # Declarar um modelo não deficiente aqui apenas para pegar a matriz. lsm <- LE_matrix(lm(mort ~ I * P * H, data = egg), effect = c("I", "P", "H")) grid <- equallevels(attr(lsm, "grid"), egg) # Celas que serão mantidas pois são estimaveis. keep <- xtabs(!is.na(oeme) ~ interaction(I, P, H), data = egg) > 0 i <- with(grid, interaction(I, P, H) %in% names(keep[keep])) grid <- grid[i, ] lsm <- lsm[i, ] # Deixa apenas as colunas de efeitos estimados. lsm <- lsm[, !is.na(b)] # Fazer as comparações apenas onde é possível. # Hospedeiros dentro de inseticida x parasitóide. rownames(lsm) <- grid$H L <- by(lsm, INDICES = with(grid, interaction(I, P)), FUN = as.matrix) i <- sapply(L, is.null) L <- L[!i] cmp <- lapply(L, apmc, model = m0, focus = "H", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- c("I", "P") names(pred)[ncol(pred)] <- "cldH" comp[[1]] <- pred # Fazer as comparações apenas onde é possível. # Inseticidas dentro de parasitóide e hospedeiro. rownames(lsm) <- grid$I L <- by(lsm, INDICES = with(grid, interaction(H, P)), FUN = as.matrix) cmp <- lapply(L, apmc, model = m0, focus = "I", test = "fdr", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- c("H", "P") names(pred)[ncol(pred)] <- "cldI" pred[, ncol(pred)] <- toupper(pred[, ncol(pred)]) comp[[2]] <- pred str(comp) pred <- merge(comp[[1]], comp[[2]], by = intersect(names(comp[[1]]), names(comp[[2]]))) pred$cld <- with(pred, paste(cldI, cldH, sep = "")) #----------------------------------------------------------------------- # Passa para a escala de probabilidade. i <- c("fit", "lwr", "upr") pred[, i] <- sapply(pred[, i], m0$family$linkinv) # Ordena da tabela. pred <- pred[with(pred, order(P, I, H)), ] # Intervalos de confiança do tamanho do suporte terão apenas o ponto # representado. i <- pred$upr - pred$lwr > 0.99 if (any(i)) { pred[i, ]$lwr <- pred[i, ]$lwr <- NA } # Reordena os níveis pela probalidade de sobreviência. pred$I <- reorder(pred$I, pred$fit) # Legenda. key <- list(points = list(pch = c(1, 19)), text = list(levels(egg_parasitoid$hosp), font = 3), title = "Hosts", cex.title = 1.1)
cap <- "Estimated parasitoid emergency proportion for each inseticide on two parasiods and two hosts. Segment is a confidence interval for the probability of surviving. Parasitoids estimates followed by the same lower letters in a insetice and host combination are not different at 5%. Inseticides estimates followed by the same lower letters in a parasitoid and host combination are not different at 5%." cap <- fgn_("pemer", cap) # Gráfico de segmentos. segplot(I ~ lwr + upr | P, centers = fit, data = pred, xlab = "Insecticides", ylab = "Proportion of parasitoid emergency", draw = FALSE, horizontal = FALSE, groups = H, key = key, strip = strip.custom( factor.levels = levels(egg_parasitoid$paras), par.strip.text = list(font = 3)), gap = 0.15, cld = pred$cld, panel = panel.groups.segplot, pch = key$points$pch[as.integer(pred$H)]) + layer({ a <- cld[which.max(nchar(cld))] l <- cld[subscripts] x <- as.integer(z)[subscripts] + centfac(groups[subscripts], gap) y <- centers[subscripts] # Usa símbolo unicode: # http://www.alanwood.net/unicode/geometric_shapes.html grid.text("\u25AE", x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "white") ) grid.text(l, x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "black", fontsize = 10)) })
#----------------------------------------------------------------------- # Análise exploratória. egg$ptot <- with(egg, femea + macho + pne) useOuterStrips(xyplot(ptot ~ I | H + P, data = egg[!is.na(egg$ptot), ], jitter.x = TRUE, auto.key = TRUE, type = c("p", "a"))) #----------------------------------------------------------------------- # Ajuste do modelo. m0 <- glm(ptot ~ I * P * H, data = egg, family = quasipoisson) par(mfrow = c(2, 2)) plot(m0) layout(1) anova(m0, test = "F") # summary(m0) # Modelo declarado com a matrix do modelo, apenas os efeitos estimáveis. X <- model.matrix(formula(m0)[-2], data = egg) b <- coef(m0) X <- X[, !is.na(b)] m0 <- update(m0, . ~ 0 + X) #----------------------------------------------------------------------- # Comparações entre hospedeiros dentro de inseticida x parasitóide. comp <- vector(mode = "list", length = 2) # Declarar um modelo não deficiente aqui apenas para pegar a matriz. lsm <- LE_matrix(lm(mort ~ I * P * H, data = egg), effect = c("I", "P", "H")) grid <- equallevels(attr(lsm, "grid"), egg) # Celas que serão mantidas pois são estimaveis. keep <- xtabs(!is.na(oeme) ~ interaction(I, P, H), data = egg) > 0 i <- with(grid, interaction(I, P, H) %in% names(keep[keep])) grid <- grid[i, ] lsm <- lsm[i, ] # Deixa apenas as colunas de efeitos estimados. lsm <- lsm[, !is.na(b)] # Fazer as comparações apenas onde é possível. # Hospedeiros dentro de inseticida x parasitóide. rownames(lsm) <- grid$H L <- by(lsm, INDICES = with(grid, interaction(I, P)), FUN = as.matrix) i <- sapply(L, is.null) L <- L[!i] cmp <- lapply(L, apmc, model = m0, focus = "H", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- c("I", "P") names(pred)[ncol(pred)] <- "cldH" comp[[1]] <- pred # Fazer as comparações apenas onde é possível. # Inseticidas dentro de parasitóide e hospedeiro. rownames(lsm) <- grid$I L <- by(lsm, INDICES = with(grid, interaction(H, P)), FUN = as.matrix) cmp <- lapply(L, apmc, model = m0, focus = "I", test = "fdr", cld2 = TRUE) pred <- ldply(cmp) cmp <- ldply(strsplit(pred$.id, "\\.")) pred <- cbind(as.data.frame(cmp), pred[, -1]) names(pred)[1:2] <- c("H", "P") names(pred)[ncol(pred)] <- "cldI" pred[, ncol(pred)] <- toupper(pred[, ncol(pred)]) comp[[2]] <- pred str(comp) pred <- merge(comp[[1]], comp[[2]], by = intersect(names(comp[[1]]), names(comp[[2]]))) pred$cld <- with(pred, paste(cldI, cldH, sep = "")) #----------------------------------------------------------------------- # Passa para a escala de probabilidade. i <- c("fit", "lwr", "upr") pred[, i] <- sapply(pred[, i], m0$family$linkinv) # Ordena da tabela. pred <- pred[with(pred, order(P, I, H)), ] # Reordena os níveis pela probalidade de sobreviência. pred$I <- reorder(pred$I, pred$fit) # Legenda. key <- list(points = list(pch = c(1, 19)), text = list(levels(egg_parasitoid$hosp), font = 3), title = "Hosts", cex.title = 1.1)
cap <- "Estimated total number of parasitoids for each inseticide on two parasiods and two hosts. Segment is a confidence interval for the probability of surviving. Parasitoids estimates followed by the same lower letters in a insetice and host combination are not different at 5%. Inseticides estimates followed by the same lower letters in a parasitoid and host combination are not different at 5%." cap <- fgn_("tot", cap) # Gráfico de segmentos. segplot(I ~ lwr + upr | P, centers = fit, data = pred, xlab = "Insecticides", ylab = "Total of parasitoids", draw = FALSE, horizontal = FALSE, groups = H, key = key, strip = strip.custom( factor.levels = levels(egg_parasitoid$paras), par.strip.text = list(font = 3)), gap = 0.15, cld = pred$cld, panel = panel.groups.segplot, pch = key$points$pch[as.integer(pred$H)]) + layer({ a <- cld[which.max(nchar(cld))] l <- cld[subscripts] x <- as.integer(z)[subscripts] + centfac(groups[subscripts], gap) y <- centers[subscripts] # Usa símbolo unicode: # http://www.alanwood.net/unicode/geometric_shapes.html grid.text("\u25AE", x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "white") ) grid.text(l, x = unit(x, "native"), y = unit(y, "native"), vjust = -0.5, gp = gpar(col = "black", fontsize = 10)) })
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"), "----------------------------------------", sep = "\n") sessionInfo()
# Log-verossimilhança para binomial. ll <- function(p, y, size) { sum(dbinom(x = y, size = size, prob = p, log = TRUE)) } # Dados observados (amostra pequena). y <- 20 size <- 20 # Grid de valores de p e sua ll correspondente. p <- seq(0.01, 0.99, by = 0.01) l <- sapply(p, FUN = ll, y = y, size = size) # Gráfico e IC baseado em corte da LL. plot(l ~ p, type = "l") abline(h = max(l) - 2 * qchisq(0.95, 1), v = 0.67, col = 2) # Agora com uma amostra bem ainda maior. y <- 200 size <- 200 # Grid de valores de p e sua ll correspondente. p <- seq(0.01, 0.99, by = 0.01) l <- sapply(p, FUN = ll, y = y, size = size) # Gráfico e IC baseado em corte da LL. plot(l ~ p, type = "l") abline(h = max(l) - 2 * qchisq(0.95, 1), v = 0.95, col = 2) #-----------------------------------------------------------------------
library(arm) # Cauchy(0, 1) prior. m0 <- bayesglm(mort ~ (I + P + H)^2, data = subset(egg, !I %in% c("Test", "Clor")), family = quasibinomial) summary(m0) anova(m0, test = "Chisq") # display(m0) # show(m0) # class(m0) # m0$data # names(m0$xlevels) # str(m0) # formula(m0) X <- unique(model.matrix(m0)) V <- vcov(m0) b <- coef(m0) str(b) str(X) X %*% b sqrt(diag(X %*% V %*% t(X))) predict(m0) methods(class = class(m0)) help(predict.bayesglm) predict(m0, newdata = pred) predict(m0) #-----------------------------------------------------------------------
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.