In walmes/RDASC: Reproducible Data Analysis of Cientific Cooperations

Definições da Sessão

#-----------------------------------------------------------------------
# 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")

Análise Exploratória

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)

Sobreviência da Fêmea

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))
        })

Ovos Parasitados

#-----------------------------------------------------------------------
# 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))
        })

Ovos Emergidos

#-----------------------------------------------------------------------
# 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))
        })

Tempo de Incubação dos Parasitóides

Essa variável apresenta pouca variação para algumas celas, portanto não informação suficiente para ser analisada.

Razão Sexual

#-----------------------------------------------------------------------
# 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))
        })

Razão de Emergência dos Parasitóides

#-----------------------------------------------------------------------
# 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))
        })

Total de Parasitóides

#-----------------------------------------------------------------------
# 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))
        })

Session information

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)

#-----------------------------------------------------------------------