##heights (cm) of plants grown under two fertilizers, Ex. 9.5 from
##Zar (1984): Biostatistical Analysis. Prentice Hall: New Jersey.
heights <- data.frame(Height = c(48.2, 54.6, 58.3, 47.8, 51.4, 52.0,
55.2, 49.1, 49.9, 52.6, 52.3, 57.4, 55.6, 53.2,
61.3, 58.0, 59.8, 54.8),
Fertilizer = c(rep("old", 10), rep("new", 8)))
##run linear model hypothesizing an effect of fertilizer
m1 <- lm(Height ~ Fertilizer, data = heights)
##run null model (no effect of fertilizer)
m0 <- lm(Height ~ 1, data = heights)
##assemble models in list
Cands <- list(m1, m0)
Modnames <- c("Fert", "null")
##compute model selection table to compare
##both hypotheses
aictab(cand.set = Cands, modnames = Modnames)
##note that model with fertilizer effect is much better supported
##than the null
##compute model-averaged effect sizes: one model hypothesizes a
##difference of 0, whereas the other assumes a difference
##prepare newdata object from which differences between groups
##will be computed
##the first row of the newdata data.frame relates to the first group,
##whereas the second row corresponds to the second group
pred.data <- data.frame(Fertilizer = c("new", "old"))
##compute best estimate of effect size accounting for model selection
##uncertainty
modavgEffect(cand.set = Cands, modnames = Modnames,
newdata = pred.data)
##classical one-way ANOVA type-design
if (FALSE) {
##generate data for two groups and control
set.seed(seed = 15)
y <- round(c(rnorm(n = 15, mean = 10, sd = 5),
rnorm(n = 15, mean = 15, sd = 5),
rnorm(n = 15, mean = 12, sd = 5)), digits = 2)
##groups
group <- c(rep("cont", 15), rep("trt1", 15), rep("trt2", 15))
##combine in data set
aov.data <- data.frame(Y = y, Group = group)
rm(y, group)
##run model with group effect
lm.eff <- lm(Y ~ Group, data = aov.data)
##null model
lm.0 <- lm(Y ~ 1, data = aov.data)
##compare both models
Cands <- list(lm.eff, lm.0)
Mods <- c("group effect", "no group effect")
aictab(cand.set = Cands, modnames = Mods)
##model with group effect has most of the weight
##compute model-averaged effect sizes
##trt1 - control
modavgEffect(cand.set = Cands, modnames = Modnames,
newdata = data.frame(Group = c("trt1", "cont")))
##trt1 differs from cont
##trt2 - control
modavgEffect(cand.set = Cands, modnames = Modnames,
newdata = data.frame(Group = c("trt2", "cont")))
##trt2 does not differ from cont
}
##two-way ANOVA type design, Ex. 13.1 (Zar 1984) of plasma calcium
##concentration (mg/100 ml) in birds as a function of sex and hormone
##treatment
if (FALSE) {
birds <- data.frame(Ca = c(16.87, 16.18, 17.12, 16.83, 17.19, 15.86,
14.92, 15.63, 15.24, 14.8, 19.07, 18.77, 17.63,
16.99, 18.04, 17.2, 17.64, 17.89, 16.78, 16.92,
32.45, 28.71, 34.65, 28.79, 24.46, 30.54, 32.41,
28.97, 28.46, 29.65),
Sex = c("M", "M", "M", "M", "M", "F", "F", "F", "F",
"F", "M", "M", "M", "M", "M", "F", "F", "F", "F",
"F", "M", "M", "M", "M", "M", "F", "F", "F", "F",
"F"),
Hormone = as.factor(c(1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3)))
##candidate models
##interactive effects
m.inter <- lm(Ca ~ Sex + Hormone + Sex:Hormone, data = birds)
##additive effects
m.add <- lm(Ca ~ Sex + Hormone, data = birds)
##Sex only
m.sex <- lm(Ca ~ Sex, data = birds)
##Hormone only
m.horm <- lm(Ca ~ Hormone, data = birds)
##null
m.0 <- lm(Ca ~ 1, data = birds)
##model selection
Cands <- list(m.inter, m.add, m.sex, m.horm, m.0)
Mods <- c("interaction", "additive", "sex only", "horm only", "null")
aictab(Cands, Mods)
##there is some support for a hormone only treatment, but also for
##additive effects
##compute model-averaged effects of sex, and set the other variable
##to a constant value
##M - F
sex.data <- data.frame(Sex = c("M", "F"), Hormone = c("1", "1"))
modavgEffect(Cands, Mods, newdata = sex.data)
##no support for a sex main effect
##hormone 1 - 3, but set Sex to a constant value
horm1.data <- data.frame(Sex = c("M", "M"), Hormone = c("1", "3"))
modavgEffect(Cands, Mods, newdata = horm1.data)
##hormone 2 - 3, but set Sex to a constant value
horm2.data <- data.frame(Sex = c("M", "M"), Hormone = c("2", "3"))
modavgEffect(Cands, Mods, newdata = horm2.data)
}
##Poisson regression with anuran larvae example from Mazerolle (2006)
if (FALSE) {
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,
family = poisson, offset = log(Effort),
data = min.trap)
Cand.mod[[2]] <- glm(Num_anura ~ log.Perimeter, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[3]] <- glm(Num_anura ~ Type, family = poisson,
offset = log(Effort), data = min.trap)
Cand.mod[[4]] <- glm(Num_anura ~ 1, family = poisson,
offset = log(Effort), data = min.trap)
##check c-hat for global model
vif.hat <- c_hat(Cand.mod[[1]]) #uses Pearson's chi-square/df
##assign names to each model
Modnames <- c("type + logperim", "type", "logperim", "intercept only")
##compute model-averaged estimate of difference between abundance at bog
##pond and upland pond
##create newdata object to make predictions
pred.data <- data.frame(Type = c("BOG", "UPLAND"),
log.Perimeter = mean(min.trap$log.Perimeter),
Effort = mean(min.trap$Effort))
modavgEffect(Cand.mod, Modnames, newdata = pred.data, c.hat = vif.hat,
type = "response")
##little suport for a pond type effect
}
##mixed linear model example from ?nlme
if (FALSE) {
library(nlme)
Cand.models <- list( )
Cand.models[[1]] <- lme(distance ~ age, data = Orthodont, method="ML")
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")
Modnames <- c("age", "age + sex", "null", "sex")
data.other <- data.frame(age = mean(Orthodont$age),
Sex = factor(c("Male", "Female")))
modavgEffect(cand.set = Cand.models, modnames = Modnames,
newdata = data.other, conf.level = 0.95, second.ord = TRUE,
nobs = NULL, uncond.se = "revised")
detach(package:nlme)
}
##site occupancy analysis example
if (FALSE) {
library(unmarked)
##single season model
data(frogs)
pferUMF <- unmarkedFrameOccu(pfer.bin)
##create a bogus site group
site.group <- c(rep(1, times = nrow(pfer.bin)/2), rep(0, nrow(pfer.bin)/2))
## add some fake covariates for illustration
siteCovs(pferUMF) <- data.frame(site.group, sitevar1 =
rnorm(numSites(pferUMF)),
sitevar2 = runif(numSites(pferUMF)))
## observation covariates are in site-major, observation-minor order
obsCovs(pferUMF) <- data.frame(obsvar1 =
rnorm(numSites(pferUMF) * obsNum(pferUMF)))
fm1 <- occu(~ obsvar1 ~ site.group, pferUMF)
fm2 <- occu(~ obsvar1 ~ 1, pferUMF)
Cand.mods <- list(fm1, fm2)
Modnames <- c("fm1", "fm2")
##model selection table
aictab(cand.set = Cand.mods, modnames = Modnames, second.ord = TRUE)
##model-averaged effect sizes comparing site.group 1 - site.group 0
newer.dat <- data.frame(site.group = c(0, 1))
modavgEffect(cand.set = Cand.mods, modnames = Modnames, type = "response",
second.ord = TRUE, newdata = newer.dat, parm.type = "psi")
##no support for an effect of site group
}
##single season N-mixture models
if (FALSE) {
data(mallard)
##this variable was created to illustrate the use of modavgEffect
##with detection variables
mallard.site$site.group <- c(rep(1, 119), rep(0, 120))
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
obsCovs = mallard.obs)
siteCovs(mallardUMF)
tmp.covs <- obsCovs(mallardUMF)
obsCovs(mallardUMF)$date2 <- tmp.covs$date^2
(fm.mall <- pcount(~ site.group ~ length + elev + forest, mallardUMF, K=30))
(fm.mallb <- pcount(~ 1 ~ length + elev + forest, mallardUMF, K=30))
Cands <- list(fm.mall, fm.mallb)
Modnames <- c("one", "null")
##model averaged effect size of site.group 1 - site.group 0 on response
##scale (point estimate)
modavgEffect(Cands, Modnames, newdata = data.frame(site.group = c(0, 1)),
parm.type = "detect", type = "response")
##model averaged effect size of site.group 1 - site.group 0 on link
##scale (here, logit link)
modavgEffect(Cands, Modnames, newdata = data.frame(site.group = c(0, 1)),
parm.type = "detect", type = "link")
detach(package:unmarked)
}
Run the code above in your browser using DataLab