##cement example in Burnham and Anderson 2002
data(cement)
##setup same model set as in Table 3.2, p. 102
Cand.models <- list( )
Cand.models[[1]] <- lm(y ~ x1 + x2, data = cement)
Cand.models[[2]] <- lm(y ~ x1 + x2 + x4, data = cement)
Cand.models[[3]] <- lm(y ~ x1 + x2 + x3, data = cement)
Cand.models[[4]] <- lm(y ~ x1 + x4, data = cement)
Cand.models[[5]] <- lm(y ~ x1 + x3 + x4, data = cement)
Cand.models[[6]] <- lm(y ~ x2 + x3 + x4, data = cement)
Cand.models[[7]] <- lm(y ~ x1 + x2 + x3 + x4, data = cement)
Cand.models[[8]] <- lm(y ~ x3 + x4, data = cement)
Cand.models[[9]] <- lm(y ~ x2 + x3, data = cement)
Cand.models[[10]] <- lm(y ~ x4, data = cement)
Cand.models[[11]] <- lm(y ~ x2, data = cement)
Cand.models[[12]] <- lm(y ~ x2 + x4, data = cement)
Cand.models[[13]] <- lm(y ~ x1, data = cement)
Cand.models[[14]] <- lm(y ~ x1 + x3, data = cement)
Cand.models[[15]] <- lm(y ~ x3, data = cement)
##vector of model names
Modnames <- paste("mod", 1:15, sep="")
##AICc
aictab(cand.set = Cand.models, modnames = Modnames)
##compute model-averaged estimate with shrinkage - each parameter
##appears 8 times in the models
modavgShrink(cand.set = Cand.models, modnames = Modnames, parm = "x1")
##compare against classic model-averaging
modavg(cand.set = Cand.models, modnames = Modnames, parm = "x1")
##note that model-averaged estimate with shrinkage is closer to 0 than
##with the classic version
##remove a few models from the set and run again
Cand.unbalanced <- Cand.models[-c(3, 14, 15)]
##set up model names
Modnames <- paste("mod", 1:length(Cand.unbalanced), sep="")
##issues an error because some parameters appear more often than others
modavgShrink(cand.set = Cand.unbalanced,
modnames = Modnames, parm = "x1")
##example on Orthodont data set in nlme
require(nlme)
##set up candidate model list
##age and sex parameters appear in the same number of models
##same number of models with and without these parameters
Cand.models <- list( )
Cand.models[[1]] <- lme(distance ~ age, data = Orthodont, method = "ML")
##random is ~ age | Subject as it is a grouped data frame
Cand.models[[2]] <- lme(distance ~ age + Sex, data = Orthodont,
random = ~ 1, method = "ML")
Cand.models[[3]] <- lme(distance ~ 1, data = Orthodont, random = ~ 1,
method = "ML")
Cand.models[[4]] <- lme(distance ~ Sex, data = Orthodont, random = ~ 1,
method = "ML")
##create a vector of model names
Modnames <- paste("mod", 1:length(Cand.models), sep = "")
##compute importance values for age
imp.age <- importance(cand.set = Cand.models, parm = "age",
modnames = Modnames, second.ord = TRUE,
nobs = NULL)
##compute shrinkage version of model averaging on age
mod.avg.age.shrink <- modavgShrink(cand.set = Cand.models,
parm = "age", modnames = Modnames,
second.ord = TRUE, nobs = NULL)
##compute classic version of model averaging on age
mod.avg.age.classic <- modavg(cand.set = Cand.models, parm = "age",
modnames = Modnames, second.ord = TRUE,
nobs = NULL)
##correspondence between shrinkage version and classic version of
##model averaging
mod.avg.age.shrink$Mod.avg.beta/imp.age$w.plus
mod.avg.age.classic$Mod.avg.beta
detach(package:nlme)
##example of N-mixture model modified from ?pcount
if(require(unmarked)) {
data(mallard)
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
obsCovs = mallard.obs)
##set up models so that each variable on abundance appears twice
fm.mall.one <- pcount(~ ivel + date ~ length + forest, mallardUMF,
K = 30)
fm.mall.two <- pcount(~ ivel + date ~ elev + forest, mallardUMF,
K = 30)
fm.mall.three <- pcount(~ ivel + date ~ length + elev, mallardUMF,
K = 30)
##model list and names
Cands <- list(fm.mall.one, fm.mall.two, fm.mall.three)
Modnames <- c("length + forest", "elev + forest", "length + elev")
##compute model-averaged estimate with shrinkage for elev on abundance
modavgShrink(cand.set = Cands, modnames = Modnames, parm = "elev",
parm.type = "lambda")
detach(package:unmarked)
}Run the code above in your browser using DataLab