Hi guys,
This Friday, 29th August at 2pm, there will be in an R-lab on mixed models in ecology, a topic which many of you are interested in digging deeper into. I've set up too links to dropbox for the slides as well as a csv file which we will be playing with (dataset courtesy of Sylvia Hay =D)
Thanks.
PDF slides
Example 1: Bird counts
Example 2: Nested design dataset
(Corrected) code for bird analysis
Yours non-significantly,
FH
(Not sure if previous comment went through)
ReplyDeleteNice slides!
* recommend glmmADMB (neg binom and other non-standard distributions) and MCMCglmm (Bayesian analysis, faster and simpler than writing your own JAGS/BUGS model) packages
* p. 10: suggest centering as well as scaling continuous predictor; suggest including continuous predictor (days) as a fixed as well as a random effect, fitting *overall* slope + random variation in slope among levels of the RE (see e.g. example(lmer)
* more worked examples at http://rpubs.com/bbolker/glmmchapter
cheers
###########
ReplyDelete## Analysis of Bird Counts
############
rm(list = ls())
library(lme4)
birds <- read.csv("birdcounts.csv",header=T,sep="\t")
## Try read.delim()?
birds$year.fac <- factor(birds$year);
birds$log.counts <- log(1+birds$counts)
fixed.fit <- lm(log.counts ~ site + year.fac, data = birds)
anova(fixed.fit)
## Check out anova()
mixed.fit <- lmer(log.counts ~ (1|site) + year.fac, data = birds)
summary(mixed.fit)
anova(mixed.fit)
## Check out summary(); anova();
## Check out plot(); ranef();
mixed.fit <- lmer(log.counts ~ (1|site) + year.fac, data = birds, REML = F)
mixed.fit2 <- lmer(log.counts ~ (1|site), data = birds, REML = F)
anova(mixed.fit, mixed.fit2)
## Do a bootstrap test of the random effect for site ##
null.fit <- lm(log.counts ~ year.fac, data = birds)
lrstat <- numeric(1000)
for(i in 1:1000) {
sim.y <- unlist(simulate(null.fit))
sim.null.fit <- lm(sim.y ~ year.fac, data = birds)
sim.mixed.fit <- lmer(sim.y ~ (1|site) + year.fac, data = birds)
lrstat[i] <- as.numeric(2*(logLik(sim.mixed.fit)-logLik(sim.null.fit)))
}
null.fit <- lm(log.counts ~ year.fac, data = birds)
mixed.fit <- lmer(log.counts ~ (1|site) + year.fac, data = birds)
lrobs <- as.numeric(2*(logLik(mixed.fit)-logLik(null.fit)))
(sum(lrstat>=lrobs)+1)/1000
############
ReplyDelete## Analysis of invertebrate (Sylvia's) dataset
################
rm(list = ls())
dat <- read.csv("Inverts-adapt.csv", header = T)
dat$creek <- factor(dat$creek)
summary(dat)
hist(dat$taxa)
plot(density(dat$taxa)) ## Looks Poisson to me!
dat$dayssinceflow.cen <- scale(dat$dayssinceflow, scale = T) ## Standardize dayssinceflow as lme4 works best when things all on roughly the same scale
library(lattice)
xyplot(taxa ~ dayssinceflow | region + creek, data = dat)
## Potential models (very much open to suggestions)!
fit1 <- glmer(taxa ~ dayssinceflow.cen*region + (1|creek), family = "poisson", data = dat)
fit2 <- glmer(taxa ~ region + (dayssinceflow.cen|creek), family = "poisson", data = dat)
fit3 <- glmer(taxa ~ region + (dayssinceflow.cen|creek/site), family = "poisson", data = dat)