##anuran larvae example from Mazerolle (2006)
data(min.trap)
##assign "UPLAND" as the reference level as in Mazerolle (2006)
min.trap$Type <- relevel(min.trap$Type, ref = "UPLAND")
##set up candidate models
Cand.mod <- list( )
##global model
Cand.mod[[1]] <- glm(Num_anura ~ Type + log.Perimeter + Num_ranatra,
family = poisson, offset = log(Effort),
data = min.trap)
Cand.mod[[2]] <- glm(Num_anura ~ Type + log.Perimeter, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[3]] <- glm(Num_anura ~ Type + Num_ranatra, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[4]] <- glm(Num_anura ~ Type, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[5]] <- glm(Num_anura ~ log.Perimeter + Num_ranatra,
family = poisson, offset = log(Effort),
data = min.trap)
Cand.mod[[6]] <- glm(Num_anura ~ log.Perimeter, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[7]] <- glm(Num_anura ~ Num_ranatra, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[8]] <- glm(Num_anura ~ 1, family = poisson,
offset = log(Effort), data = min.trap)
##check c-hat for global model
c_hat(Cand.mod[[1]]) #uses Pearson's chi-square/df
##note the very low overdispersion: in this case, the analysis could be
##conducted without correcting for c-hat as its value is reasonably close
##to 1
##assign names to each model
Modnames <- c("type + logperim + invertpred", "type + logperim",
"type + invertpred", "type", "logperim + invertpred",
"logperim", "invertpred", "intercept only")
##compute model-averaged estimate of TypeBOG
modavg(parm = "TypeBOG", cand.set = Cand.mod, modnames = Modnames)
##round to 4 digits after decimal point
print(modavg(parm = "TypeBOG", cand.set = Cand.mod,
modnames = Modnames), digits = 4)
##compute c-hat estimate based on residual deviance as in Mazerolle
##(2006)
Cand.mod[[1]]$deviance/Cand.mod[[1]]$df.residual
##compute model-averaged estimate of TypeBOG as in Table 4 of
##Mazerolle (2006)
modavg(parm = "TypeBOG", cand.set = Cand.mod, modnames = Modnames,
c.hat = 1.11)
##example with similarly-named variables and interaction terms
set.seed(seed = 4)
resp <- rnorm(n = 40, mean = 3, sd = 1)
size <- rep(c("small", "medsmall", "high", "medhigh"), times = 10)
set.seed(seed = 4)
mass <- rnorm(n = 40, mean = 2, sd = 0.1)
mass2 <- mass^2
age <- rpois(n = 40, lambda = 3.2)
agecorr <- rpois(n = 40, lambda = 2)
sizecat <- rep(c("a", "ab"), times = 20)
data1 <- data.frame(resp = resp, size = size, sizecat = sizecat,
mass = mass, mass2 = mass2, age = age,
agecorr = agecorr)
##set up models in list
Cand <- list( )
Cand[[1]] <- lm(resp ~ size + agecorr, data = data1)
Cand[[2]] <- lm(resp ~ size + mass + agecorr, data = data1)
Cand[[3]] <- lm(resp ~ age + mass, data = data1)
Cand[[4]] <- lm(resp ~ age + mass + mass2, data = data1)
Cand[[5]] <- lm(resp ~ mass + mass2 + size, data = data1)
Cand[[6]] <- lm(resp ~ mass + mass2 + sizecat, data = data1)
Cand[[7]] <- lm(resp ~ sizecat, data = data1)
Cand[[8]] <- lm(resp ~ sizecat + mass + sizecat:mass, data = data1)
Cand[[9]] <- lm(resp ~ agecorr + sizecat + mass + sizecat:mass,
data = data1)
##create vector of model names
Modnames <- paste("mod", 1:length(Cand), sep = "")
aictab(cand.set = Cand, modnames = Modnames, sort = TRUE) #correct
##as expected, issues warning as mass occurs sometimes with "mass2" or
##"sizecatab:mass" in some of the models
modavg(cand.set = Cand, parm = "mass", modnames = Modnames)
##no warning issued, because "age" and "agecorr" never appear in same model
modavg(cand.set = Cand, parm = "age", modnames = Modnames)
##as expected, issues warning because warn=FALSE, but it is a very bad
##idea in this example since "mass" occurs with "mass2" and "sizecat:mass"
##in some of the models - results are INCORRECT
modavg(cand.set = Cand, parm = "mass", modnames = Modnames,
warn = FALSE)
##correctly excludes models with quadratic term and interaction term
##results are CORRECT
modavg(cand.set = Cand, parm = "mass", modnames = Modnames,
exclude = list("mass2", "sizecat:mass"))
##correctly computes model-averaged estimate because no other parameter
##occurs simultaneously in any of the models
modavg(cand.set = Cand, parm = "sizesmall", modnames = Modnames) #correct
##as expected, issues a warning because "sizecatab" occurs sometimes in
##an interaction in some models
modavg(cand.set = Cand, parm = "sizecatab",
modnames = Modnames)
##exclude models with "sizecat:mass" interaction - results are CORRECT
modavg(cand.set = Cand, parm = "sizecatab", modnames = Modnames,
exclude = list("sizecat:mass"))
##example with generalized linear mixed model
##modified example from ?glmer
require(lme4)
##create proportion of incidence
cbpp$prop <- cbpp$incidence/cbpp$size
##add bogus variable
cbpp$randx <- rnorm(n = nrow(cbpp), mean = 12, sd = 3)
##model with period
gm1 <- glmer(prop ~ period + (1 | herd), family = binomial,
weights = size, data = cbpp)
gm2 <- glmer(prop ~ period + randx + (1 | herd), family = binomial,
weights = size, data = cbpp)
##model without period
gm3 <- glmer(prop ~ 1 + (1 | herd), family = binomial,
weights = size, data = cbpp)
Cands <- list(gm1, gm2, gm3)
Modnames <- c("period", "period + randx", "null")
##model selection
aictab(cand.set = Cands, modnames = Modnames)
##model average for difference between period 1 vs period 4
modavg(cand.set = Cands, modnames = Modnames, parm = "period4")
detach(package:lme4)
##modified example of Cox regression from ?coxph
require(survival)
##Create a simple data set for a time-dependent model
test2 <- list(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
x =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0),
height = c(12.3, 10.5, 9.2, 5.6, 8.9,
11.0, 16.1, 10.2, 9.9, 14.8))
m.cox <- coxph(Surv(start, stop, event) ~ x, test2)
m.cox2 <- coxph(Surv(start, stop, event) ~ x + height, test2)
Cands <- list(m.cox, m.cox2)
Mods <- c("x", "additive")
aictab(Cands, Mods)
modavg(Cands, Mods, parm = "x")
##example with multiple-season occupancy model modified from ?colext
##this is a bit longer
require(unmarked)
data(frogs)
umf <- formatMult(masspcru)
obsCovs(umf) <- scale(obsCovs(umf))
siteCovs(umf) <- rnorm(numSites(umf))
yearlySiteCovs(umf) <- data.frame(year = factor(rep(1:7,
numSites(umf))))
##set up model with constant transition rates
fm <- colext(psiformula = ~ 1, gammaformula = ~ 1, epsilonformula = ~ 1,
pformula = ~ JulianDate + I(JulianDate^2), data = umf,
control = list(trace=1, maxit=1e4))
##model with with year-dependent transition rates
fm.yearly <- colext(psiformula = ~ 1, gammaformula = ~ year,
epsilonformula = ~ year,
pformula = ~ JulianDate + I(JulianDate^2),
data = umf)
##store in list and assign model names
Cand.mods <- list(fm, fm.yearly)
Modnames <- c("psi1(.)gam(.)eps(.)p(Date + Date2)",
"psi1(.)gam(Year)eps(Year)p(Date + Date2)")
##compute model-averaged estimate of occupancy in the first year
modavg(cand.set = Cand.mods, modnames = Modnames, parm = "(Intercept)",
parm.type = "psi")
##compute model-averaged estimate of Julian Day squared on detectability
modavg(cand.set = Cand.mods, modnames = Modnames,
parm = "I(JulianDate^2)", parm.type = "detect")
##example of model-averaged estimate of area from distance model
##this is a bit longer
data(linetran) #example modified 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 detection.
fm2 <- distsamp(~ area + habitat ~ area + habitat, ltUMF)
## Hazard function. Covariates affecting both density and detection.
fm3 <- distsamp(~ habitat ~ area + habitat, ltUMF, keyfun="hazard")
##assemble model list
Cands <- list(fm1, fm2, fm3)
Modnames <- paste("mod", 1:length(Cands), sep = "")
##model-average estimate of area on abundance
modavg(cand.set = Cands, modnames = Modnames, parm = "area", parm.type = "lambda")
detach(package:unmarked)Run the code above in your browser using DataLab