# Flour beetle mortality data Mortality of flour beetles (**Tribolium confusum**) due to exposure to gaseous carbon disulfide CS$2$, from Bliss (1935). utf-8 data ```r Beetle ``` ## Format `Beetle` is a data frame with 5 elements. - **Prop**: a matrix with two columns named nkilled and nsurvived - **mortality**: observed mortality rate - **dose**: the dose of CS$2$ in mg/L - **n.tested**: number of beetles tested - **n.killed**: number of beetles killed. ## Source Bliss, C. I. 1935 The calculation of the dosage-mortality curve. **Annals of Applied Biology** 22 , 134--167. ## References Burnham, K. P. and Anderson, D. R. 2002 **Model selection and multimodel inference: a practical information-theoretic approach**. 2nd ed. New York, Springer-Verlag. ## Examples ```r # "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) 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)) ```

data-Beetle function