[R-br] dose-resposta

[Walmes Zeviani]

Ajuste um modelo com as três espécies como se fosse um fatorial. Para obter
as DLs é mais fácil usar a parametrização de "estimativas separadas por
nível". Veja o código abaixo.

rm(list = ls())
da <- data.frame(
    dose = rep(c(0, 0.15625, 0.3125, 0.625, 1.25, 2.5, 5, 10), each = 4),
    n = rep(10, 32),
    m1 = c(1, 3, 4, 0, 5, 2, 5, 5, 4, 4, 2, 2, 0, 3, 3, 5, 10,
           10, 7, 0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
           10),
    m2 = c(0, 3, 4, 4, 1, 0, 4, 0, 0, 1, 0, 0, 3, 2, 0, 3, 3, 2,
           0, 3, 3, 3, 3, 6, 9, 4, 7, 2, 10, 10, 10, 10),
    m3 = c(4, 2, 2, 3, 4, 6, 8, 6, 4, 6, 6, 5, 5, 9, 5, 8, 10,
            8, 5, 5, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
            10))
da$Dose <- as.factor(da$dose)

da <- reshape2::melt(da,
                     id.vars = c("dose", "Dose", "n"),
                     variable.name = "esp",
                     value.name = "mort")
str(da)

library(lattice)

bwplot(mort/n ~ Dose | esp,
       data = da,
       pch = "|",
       as.table = TRUE,
       layout = c(NA, 1))

xyplot(mort/n ~ Dose | esp,
       data = da,
       type = c("p", "a"),
       as.table = TRUE,
       layout = c(NA, 1))

m0 <- glm(cbind(mort, n - mort) ~ esp * dose,
          data = da,
          family = quasibinomial)

par(mfrow = c(2, 2))
plot(m0)
layout(1)

summary(m0)

anova(m0, test = "F")

# Em um modelo binomial do tipo ~ \beta_0 + \beta_1 * x, a DL_{50} é
# -\beta_0/\beta_1.

# Ajuste do modelo com estimativas separadas para cada nível.
m1 <- glm(cbind(mort, n - mort) ~ 0 + esp/dose,
          data = da,
          family = quasibinomial)

# Deviances iguais porque são modelos iguais.
deviance(m0)
deviance(m1)

summary(m1)

beta <- matrix(coef(m1), ncol = 2)
dl <- -beta[, 1]/beta[, 2]

pred <- with(da,
             expand.grid(esp = levels(esp),
                         dose = seq(min(dose), max(dose),
                                    length.out = 30)))

pred$p <- predict(m1, newdata = pred, type = "response")

library(latticeExtra)

xyplot(mort/n ~ dose | esp,
       data = da,
       as.table = TRUE,
       layout = c(NA, 1)) +
    as.layer(xyplot(p ~ dose | esp,
                    data = pred,
                    type = "l",
                    col = 2,
                    lwd = 2)) +
    layer(panel.abline(v = dl[which.packet()],
                       h = 0.5,
                       lty = 2))

À disposição.
Walmes.
​
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20161222/f03d87aa/attachment.html>