##example from subset of models in Table 1 in Mazerolle (2006)
data(dry.frog)
Cand.models <- list( )
Cand.models[[1]] <- lm(log_Mass_lost ~ Shade + Substrate +
cent_Initial_mass + Initial_mass2,
data = dry.frog)
Cand.models[[2]] <- lm(log_Mass_lost ~ Shade + Substrate +
cent_Initial_mass + Initial_mass2 +
Shade:Substrate, data = dry.frog)
Cand.models[[3]] <- lm(log_Mass_lost ~ cent_Initial_mass +
Initial_mass2, data = dry.frog)
Cand.models[[4]] <- lm(log_Mass_lost ~ Shade + cent_Initial_mass +
Initial_mass2, data = dry.frog)
Cand.models[[4]] <- lm(log_Mass_lost ~ Shade + cent_Initial_mass +
Initial_mass2, data = dry.frog)
Cand.models[[5]] <- lm(log_Mass_lost ~ Substrate + cent_Initial_mass +
Initial_mass2, data = dry.frog)
##setup model names
Modnames <- paste("mod", 1:length(Cand.models), sep = "")
##compute model-averaged value and unconditional SE of predicted log of
##mass lost for frogs of average mass in shade for each substrate type
##first create data set to use for predictions
new.dat <- data.frame(Shade = c(1, 1, 1),
cent_Initial_mass = c(0, 0, 0),
Initial_mass2 = c(0, 0, 0),
Substrate = c("SOIL", "SPHAGNUM", "PEAT"))
##compare unconditional SE's using both methods
modavgPred(cand.set = Cand.models, modnames = Modnames,
newdata = new.dat, type = "response", uncond.se = "old")
modavgPred(cand.set = Cand.models, modnames = Modnames,
newdata = new.dat, type = "response", uncond.se = "revised")
##round to 4 digits after decimal point
print(modavgPred(cand.set = Cand.models, modnames = Modnames,
newdata = new.dat, type = "response",
uncond.se = "revised"), digits = 4)
##Gamma glm
##clotting data example from 'gamma.shape' in MASS package of
##Venables and Ripley (2002, Modern applied statistics with
##S. Springer-Verlag: New York.)
clotting <- data.frame(u = c(5, 10, 15, 20, 30, 40, 60, 80, 100),
lot1 = c(118, 58, 42, 35, 27, 25, 21, 19, 18),
lot2 = c(69, 35, 26, 21, 18, 16, 13, 12, 12))
clot1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
library(MASS)
gamma.dispersion(clot1) #dispersion parameter
gamma.shape(clot1) #reciprocal of dispersion parameter ==
##shape parameter
summary(clot1, dispersion = gamma.dispersion(clot1)) #better
##create list with models
Cand <- list( )
Cand[[1]] <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
Cand[[2]] <- glm(lot1 ~ 1, data = clotting, family = Gamma)
##create vector of model names
Modnames <- paste("mod", 1:length(Cand), sep = "")
##compute model-averaged predictions on scale of response variable for
##all observations
modavgPred(cand.set = Cand, modnames = Modnames, newdata = clotting,
gamdisp = gamma.dispersion(clot1), type = "response")
##compute model-averaged predictions on scale of linear predictor
modavgPred(cand.set = Cand, modnames = Modnames, newdata = clotting,
gamdisp = gamma.dispersion(clot1), type = "link")
##compute model-averaged predictions on scale of linear predictor
modavgPred(cand.set = Cand, modnames = Modnames, newdata = clotting,
gamdisp = gamma.dispersion(clot1), type = "terms") #returns an error
##because type = "terms" is not defined for 'modavgPred'
modavgPred(cand.set = Cand, modnames = Modnames, newdata = clotting,
type = "terms") #returns an error because
##no gamma dispersion parameter was specified (i.e., 'gamdisp' missing)
##example of model-averaged predictions from N-mixture model
##each variable appears twice in the models - this is a bit longer
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)
fm.mall.four <- pcount(~ ivel + date ~ 1, mallardUMF, K = 30)
##model list
Cands <- list(fm.mall.one, fm.mall.two, fm.mall.three, fm.mall.four)
Modnames <- c("length + forest", "elev + forest", "length + elev",
"null")
##compute model-averaged predictions of abundance for values of elev
modavgPred(cand.set = Cands, modnames = Modnames, newdata =
data.frame(elev = seq(from = -1.4, to = 2.4, by = 0.1),
length = 0, forest = 0), parm.type = "lambda",
type = "response")
##compute model-averaged predictions of detection for values of ivel
modavgPred(cand.set = Cands, modnames = Modnames, newdata =
data.frame(ivel = seq(from = -1.75, to = 5.9, by = 0.5),
date = 0), parm.type = "detect",
type = "response")
}
##example of model-averaged abundance from distance model
##this is a bit longer
\dontrun{
data(linetran) #example from ?distsamp
ltUMF <- with(linetran, {
unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4),
siteCovs = data.frame(Length, area, habitat),
dist.breaks = c(0, 5, 10, 15, 20),
tlength = linetran$Length * 1000, survey = "line", unitsIn = "m")
})
## Half-normal detection function. Density output (log scale). No covariates.
fm1 <- distsamp(~ 1 ~ 1, ltUMF)
## Halfnormal. Covariates affecting both density and and detection.
fm2 <- distsamp(~area + habitat ~ habitat, ltUMF)
## Hazard function. Covariates affecting both density and and detection.
fm3 <- distsamp(~area + habitat ~ habitat, ltUMF, keyfun="hazard")
##assemble model list
Cands <- list(fm1, fm2, fm3)
Modnames <- paste("mod", 1:length(Cands), sep = "")
##model-average predictions on abundance
modavgPred(cand.set = Cands, modnames = Modnames, parm.type = "lambda", type = "link",
newdata = data.frame(area = mean(linetran$area), habitat = c("A", "B")))
detach(package:unmarked)
}
##example using Orthodont data set from Pinheiro and Bates (2000)
require(nlme)
##set up candidate models
m1 <- gls(distance ~ age, correlation = corCompSymm(value = 0.5, form = ~ 1 | Subject),
data = Orthodont, method = "ML")
m2 <- gls(distance ~ 1, correlation = corCompSymm(value = 0.5, form = ~ 1 | Subject),
data = Orthodont, method = "ML")
##assemble in list
Cand.models <- list(m1, m2)
##model names
Modnames <- c("age effect", "null model")
##model selection table
aictab(cand.set = Cand.models, modnames = Modnames)
##model-averaged predictions
modavgPred(cand.set = Cand.models, modnames = Modnames, newdata =
data.frame(age = c(8, 10, 12, 14)))
detach(package:nlme)Run the code above in your browser using DataLab