data-Beetle: Flour beetle mortality data

BeetleR Documentation

Flour beetle mortality data

Description

Mortality of flour beetles (Tribolium confusum) due to exposure to gaseous carbon disulfide CS2, from Bliss (1935).

Usage

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 CS2 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


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


MuMIn documentation built on March 18, 2022, 5:28 p.m.