## Wednesday, 27 August 2014

### R-lab on Mixed models in Ecology

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

1. (Not sure if previous comment went through)

Nice 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

2. ###########
## Analysis of Bird Counts
############

rm(list = ls())
library(lme4)

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

3. ############
## Analysis of invertebrate (Sylvia's) dataset
################

rm(list = ls())
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)