AICcmodavg.R
In AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c)

### R code from vignette source 'AICcmodavg.Rnw'

###################################################
### code chunk number 1: AICcmodavg.Rnw:37-38
###################################################
options(width=70, continue = "  ")


###################################################
### code chunk number 2: import
###################################################
library(AICcmodavg)
data(dry.frog)


###################################################
### code chunk number 3: subData
###################################################
##extract only first 7 columns
frog <- dry.frog[, 1:7]
##first lines
head(frog)
##structure of data frame
str(frog)


###################################################
### code chunk number 4: na
###################################################
any(is.na(frog))


###################################################
### code chunk number 5: centInitialMass
###################################################
##center initial mass
frog$InitMass_cent <- frog$Initial_mass - mean(frog$Initial_mass)


###################################################
### code chunk number 6: InitialMass2
###################################################
frog$InitMass2 <- frog$InitMass_cent^2


###################################################
### code chunk number 7: checkDiag
###################################################
##run global model
global <- lm(Mass_lost ~ InitMass_cent + InitMass2 + Substrate + Shade, 
             data = frog)
par(mfrow = c(2, 2))
plot(global)


###################################################
### code chunk number 8: logMass
###################################################
frog$logMass_lost <- log(frog$Mass_lost + 1) #adding 1 due to presence of 0's


###################################################
### code chunk number 9: checkDiag2
###################################################
##run global model
global.log <- lm(logMass_lost ~ InitMass_cent + InitMass2 + Substrate + Shade, 
                 data = frog)
par(mfrow = c(2, 2))
plot(global.log)


###################################################
### code chunk number 10: fitCands
###################################################
m.null <- lm(logMass_lost ~ 1,
             data = frog)
m.shade <- lm(logMass_lost ~ Shade, 
              data = frog)
m.substrate <- lm(logMass_lost ~ Substrate, 
                  data = frog)
m.shade.substrate <- lm(logMass_lost ~ Shade + Substrate, 
                        data = frog)
m.null.mass <- lm(logMass_lost ~ InitMass_cent + InitMass2,
                  data = frog)
m.shade.mass <- lm(logMass_lost ~ InitMass_cent + InitMass2 + Shade, 
                   data = frog)
m.substrate.mass <- lm(logMass_lost ~ InitMass_cent + InitMass2 + Substrate, 
                       data = frog)
m.global.mass <- global.log


###################################################
### code chunk number 11: storeList
###################################################
##store models in named list
Cand.models <- list("null" = m.null, "shade" = m.shade, 
                    "substrate" = m.substrate, 
                    "shade + substrate" = m.shade.substrate, 
                    "mass" = m.null.mass, "mass + shade" = m.shade.mass, 
                    "mass + substrate" = m.substrate.mass, 
                    "global" = m.global.mass)


###################################################
### code chunk number 12: modTableAICc
###################################################
selectionTable <- aictab(cand.set = Cand.models)
selectionTable


###################################################
### code chunk number 13: modTableAIC
###################################################
aictab(Cand.models, second.ord = FALSE)


###################################################
### code chunk number 14: exportTable3 (eval = FALSE)
###################################################
## library(xtable)
## print(xtable(selectionTable, caption = "Model selection table on frog mass lost.", 
##              label = "tab:selection"),
##       include.rownames = FALSE, caption.placement = "top")


###################################################
### code chunk number 15: exportTable4
###################################################
library(xtable)
print(xtable(selectionTable, caption = "Model selection table on frog mass lost.", 
             label = "tab:selection"),
      include.rownames = FALSE, caption.placement = "top", )


###################################################
### code chunk number 16: confSet
###################################################
##confidence set of models
confset(cand.set = Cand.models)


###################################################
### code chunk number 17: evidence
###################################################
##evidence ratios
evidence(aic.table = selectionTable)


###################################################
### code chunk number 18: evidenceSilent
###################################################
evRatio <- evidence(selectionTable)


###################################################
### code chunk number 19: evidence2
###################################################
##compare "substrate" vs "shade"
evidence(selectionTable, model.high = "substrate",
         model.low = "shade")


###################################################
### code chunk number 20: evidence2Silent
###################################################
##compare "substrate" vs "shade"
evRatio2 <- evidence(selectionTable, model.high = "substrate",
                     model.low = "shade")


###################################################
### code chunk number 21: evidenceNull
###################################################
evidence(selectionTable, model.high = "global",
         model.low = "null")


###################################################
### code chunk number 22: confint
###################################################
confint(m.global.mass)


###################################################
### code chunk number 23: modavg
###################################################
modavg(cand.set = Cand.models, parm = "Shade")


###################################################
### code chunk number 24: modavg2
###################################################
modavgShade <- modavg(cand.set = Cand.models, parm = "Shade")


###################################################
### code chunk number 25: coef
###################################################
coef(m.global.mass)


###################################################
### code chunk number 26: substrateSPHAG
###################################################
modavg(Cand.models, parm = "SubstrateSPHAGNUM")


###################################################
### code chunk number 27: substrateSPHAG2
###################################################
modavgSphag <- modavg(Cand.models, parm = "SubstrateSPHAGNUM")


###################################################
### code chunk number 28: modavg
###################################################
modavgShrink(cand.set = Cand.models, parm = "Shade")


###################################################
### code chunk number 29: substrateSPHAGShrink
###################################################
modavgShrink(Cand.models, parm = "SubstrateSPHAGNUM")


###################################################
### code chunk number 30: shadePred
###################################################
##data frame to make predictions
##all variables are held constant, except Shade
predData <- data.frame(InitMass_cent = c(0, 0),
                       InitMass2 = c(0, 0),
                       Substrate = factor("SOIL", 
                                          levels = levels(frog$Substrate)),
                       Shade = c(0, 1))
##predictions from global model
predict(m.global.mass, newdata = predData, se.fit = TRUE)
##predictions from null model
predict(m.null, newdata = predData, se.fit = TRUE)


###################################################
### code chunk number 31: extractX
###################################################
extractX(cand.set = Cand.models)


###################################################
### code chunk number 32: modavgPred
###################################################
modavgPred(cand.set = Cand.models, newdata = predData)


###################################################
### code chunk number 33: modavgPredSub
###################################################
##data frame holding all variables constant, except Substrate
predSub <- data.frame(InitMass_cent = c(0, 0, 0),
                      InitMass2 = c(0, 0, 0),
                      Substrate = factor(c("PEAT", "SOIL", "SPHAGNUM"),
                                         levels = levels(frog$Substrate)),
                      Shade = c(1, 1, 1))
##model-average predictions
predsMod <- modavgPred(Cand.models, newdata = predSub)
predsMod


###################################################
### code chunk number 34: checkContent
###################################################
##check content of object
str(predsMod)


###################################################
### code chunk number 35: savePreds
###################################################
##add predictions, lower CL, and upper CL
predSub$fit <- predsMod$mod.avg.pred
predSub$low95 <- predsMod$lower.CL
predSub$upp95 <- predsMod$upper.CL


###################################################
### code chunk number 36: plotPreds
###################################################
##create vector for X axis
predSub$xvals <- c(0.25, 0.5, 0.75)
##create empty box
plot(fit ~ xvals, 
     data = predSub,
     xlim = c(0, 1),
     ylim = range(low95, upp95),
     xlab = "Substrate type",
     ylab = "Predicted mass lost (log of mass in g)",
     xaxt = "n",
     cex.axis = 1.2, 
     cex.lab = 1.2)
##add x axis
axis(side = 1, at = predSub$xvals,
     labels = c("Peat", "Soil", "Sphagnum"),
     cex.axis = 1.2)
##add CI's
segments(x0 = predSub$xvals, x1 = predSub$xvals,
         y0 = predSub$low95, predSub$upp95)


###################################################
### code chunk number 37: compGroups
###################################################
predComp <- data.frame(InitMass_cent = c(0, 0),
                       InitMass2 = c(0, 0),
                       Substrate = factor(c("PEAT", "SPHAGNUM"),
                                          levels = levels(frog$Substrate)),
                       Shade = c(1, 1))
##model-average predictions
modavgEffect(Cand.models, newdata = predComp)


###################################################
### code chunk number 38: customAICc
###################################################
##log-likelihoods
modL <- c(-225.4180, -224.0697, -225.4161)
##number of parameters
modK <- c(2, 3, 3)
##model selection
outTab <- aictabCustom(logL = modL,
                       K = modK,
                       modnames = c("null", "phi(SVL)p(.)",
                                    "phi(Road)p(.)"),
                       nobs = 621)


###################################################
### code chunk number 39: evRatioCustom (eval = FALSE)
###################################################
## evidence(outTab, model.high = "phi(SVL)p(.)",
##          model.low = "phi(Road)p(.)")


###################################################
### code chunk number 40: evRatioCustom
###################################################
evRatioCust <- evidence(outTab, model.high = "phi(SVL)p(.)",
                        model.low = "phi(Road)p(.)")


###################################################
### code chunk number 41: estSE
###################################################
##survival estimates with road mitigation
modEst <- c(0.1384450, 0.1266030, 0.1378745)
##SE's of survival estimates with road mitigation
modSE <- c(0.03670327, 0.03347475, 0.03862634)


###################################################
### code chunk number 42: customModavg
###################################################
##model-averaged survival with road mitigation
modavgCustom(logL = modL,
             K = modK,
             modnames = c("null", "phi(SVL)p(.)",
                          "phi(Road)p(.)"),
             estimate = modEst,
             se = modSE,
             nobs = 621)


###################################################
### code chunk number 43: customModavg2
###################################################
##survival estimates without road mitigation
modEst2 <- c(0.1384450, 0.1266030, 0.1399727)

##SE's of survival estimates without road mitigation
modSE2 <- c(0.03670327, 0.03347475, 0.04981635)

##model-averaged survival
modavgCustom(logL = modL,
             K = modK,
             modnames = c("null", "phi(SVL)p(.)",
                          "phi(Road)p(.)"),
             estimate = modEst2,
             se = modSE2,
             nobs = 621)