Mortality of flour beetles (Tribolium confusum) due to exposure to gaseous
carbon disulfide CS
Beetle
Beetle
is a data frame with 5 elements.
a matrix with two columns named nkilled and nsurvived
observed mortality rate
the dose of CS
number of beetles tested
number of beetles killed.
Burnham, K. P. and Anderson, D. R. (2002) Model selection and multimodel inference: a practical information-theoretic approach. 2nd ed. New York, Springer-Verlag.
# NOT RUN {
# "Logistic regression example"
# from Burnham & Anderson (2002) chapter 4.11
# Fit a global model with all the considered variables
globmod <- glm(Prop ~ dose + I(dose^2) + log(dose) + I(log(dose)^2),
data = Beetle, family = binomial, na.action = na.fail)
# A logical expression defining the subset of models to use:
# * either log(dose) or dose
# * the quadratic terms can appear only together with linear terms
msubset <- expression(xor(dose, `log(dose)`) &
dc(dose, `I(dose^2)`) &
dc(`log(dose)`, `I(log(dose)^2)`))
# Table 4.6
# Use 'varying' argument to fit models with different link functions
# Note the use of 'alist' rather than 'list' in order to keep the
# 'family' objects unevaluated
varying.link <- list(family = alist(
logit = binomial("logit"),
probit = binomial("probit"),
cloglog = binomial("cloglog")
))
(ms12 <- dredge(globmod, subset = msubset, varying = varying.link,
rank = AIC))
# Table 4.7 "models justifiable a priori"
(ms3 <- subset(ms12, has(dose, !`I(dose^2)`)))
# The same result, but would fit the models again:
# ms3 <- update(ms12, update(globmod, . ~ dose), subset =,
# fixed = ~dose)
mod3 <- get.models(ms3, 1:3)
# Table 4.8. Predicted mortality probability at dose 40.
# calculate confidence intervals on logit scale
logit.ci <- function(p, se, quantile = 2) {
C. <- exp(quantile * se / (p * (1 - p)))
p /(p + (1 - p) * c(C., 1/C.))
}
mavg3 <- model.avg(mod3, revised.var = FALSE)
# get predictions both from component and averaged models
pred <- lapply(c(component = mod3, list(averaged = mavg3)), predict,
newdata = list(dose = 40), type = "response", se.fit = TRUE)
# reshape predicted values
pred <- t(sapply(pred, function(x) unlist(x)[1:2]))
colnames(pred) <- c("fit", "se.fit")
# build the table
tab <- cbind(
c(Weights(ms3), NA),
pred,
matrix(logit.ci(pred[,"fit"], pred[,"se.fit"],
quantile = c(rep(1.96, 3), 2)), ncol = 2)
)
colnames(tab) <- c("Akaike weight", "Predicted(40)", "SE", "Lower CI",
"Upper CI")
rownames(tab) <- c(as.character(ms3$family), "model-averaged")
print(tab, digits = 3, na.print = "")
# Figure 4.3
newdata <- list(dose = seq(min(Beetle$dose), max(Beetle$dose), length.out = 25))
# add model-averaged prediction with CI, using the same method as above
avpred <- predict(mavg3, newdata, se.fit = TRUE, type = "response")
avci <- matrix(logit.ci(avpred$fit, avpred$se.fit, quantile = 2), ncol = 2)
# }
# NOT RUN {
matplot(newdata$dose, sapply(mod3, predict, newdata, type = "response"),
type = "l", xlab = quote(list("Dose of" ~ CS[2],(mg/L))),
ylab = "Mortality", col = 2:4, lty = 3, lwd = 1
)
matplot(newdata$dose, cbind(avpred$fit, avci), type = "l", add = TRUE,
lwd = 1, lty = c(1, 2, 2), col = 1)
legend("topleft", NULL, c(as.character(ms3$family), expression(`averaged`
%+-% CI)), lty = c(3, 3, 3, 1), col = c(2:4, 1))
# }
Run the code above in your browser using DataLab