https://github.com/cran/MuMIn
Tip revision: 1834bb90bade3912317a15c9b7c19771f19cf6dc authored by Kamil BartoĊ on 22 June 2024, 14:10:02 UTC
version 1.48.4
version 1.48.4
Tip revision: 1834bb9
data-Beetle.Rd
\name{Beetle}
\alias{Beetle}
\encoding{utf-8}
\docType{data}
\title{Flour beetle mortality data}
\description{
Mortality of flour beetles (\emph{Tribolium confusum}) due to exposure to gaseous
carbon disulfide CS\eqn{_{2}}{2}, from Bliss (1935).
}
\usage{
Beetle
}
\format{
\code{Beetle} is a data frame with 5 elements.
\describe{
\item{Prop}{a matrix with two columns named \bold{nkilled} and \bold{nsurvived}}
\item{mortality}{observed mortality rate}
\item{dose}{the dose of CS\eqn{_{2}}{2} in mg/L}
\item{n.tested}{number of beetles tested}
\item{n.killed}{number of beetles killed.}
}
}
\source{
Bliss, C. I. 1935 The calculation of the dosage-mortality curve.
\emph{Annals of Applied Biology} \bold{22}, 134--167.
}
\references{
Burnham, K. P. and Anderson, D. R. 2002 \emph{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)
\dontshow{ if(require(graphics)) \{ }
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))
\dontshow{ \} }
}
\keyword{datasets}
%% <eof>