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)