How to run a two-way ANOVA with a random variable followed by pairwise comparisons?

Question

I have a dataset with which I would like to compare the effect of species and habitat on movement rate, with pairwise comparisons. I would also like to include the effect of individual (as a random factor?) - this random factor is the part I don't know how to do, at least not in the framework of Anova().

Here's a subset of the data:

species  <- c("a", "b", "c", "a", "a", "b", "c", "a", "a", "b", "c", "a", "a", "b", "c", 
              "a", "a", "b", "c", "a")
habitat  <- c("x", "x", "x", "y", "y", "y", "x", "x", "y", "z", "y", "y", "z", "z", "x", 
              "x", "y", "y", "z", "z")
mvt.rate <- c(6, 5, 7, 8, 9, 4, 3, 5, 6, 9, 3, 6, 6, 7, 8, 9, 5, 6, 7, 8)
ind      <- as.factor(c(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4))
data1    <- data.frame(species, habitat, mvt.rate, ind)

Currently, I'm simply running a two-way ANOVA with pairwise comparisons, without considering the effect of individual, like so:

fit <- lm(mvt.rate ~ habitat + species, data=data1)
require(car)
Anova(fit, type="III")
require(agricolae)
    #pairwise comparison of habitats
comparison.hab <- HSD.test(fit, "habitat", group=TRUE)                     
    #pairwise comparison of species
comparison.sp <- HSD.test(fit, "species", group=TRUE)

In the dataset, each row represents a movement, and in many cases, individuals make several (non-independent) movements - I am currently not considering this non-independence of mvt.rate and individual. I believe the correct way to do this is to consider individual as a random variable, but I'm not entirely sure.

Answer 1

I can show you how to model the individuals as a random effect. You can simply use use the lme function from the nlme package. The syntax will be very similar to lm:

require(nlme)
fit<-lme(mvt.rate ~ habitat +species, random=~1|ind, data=data1)

As for the pairwise comparisons, for lme you can use the glht function:

require(multcomp)
summary(glht(fit,  linfct=mcp(habitat = "Tukey")))

Answer 2

Just remove the variability that is related with individuals: pre.fit=lm(mvt.rate ~ind,data=data1) and then fit<-lm(pre.fit$residuals ~ data1$habitat +data1$species)