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

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.

-
You forgot to include the values of homerange in your question. –  AlefSin Sep 15 '12 at 15:42
@Luke I edited the example slightly as it contained an error. Feel free to revert it. –  Henrik Sep 15 '12 at 20:13
Some individuals are belonging to more than one species. Does this make any sense? (I don't think so.) If you can clarify that, I post another answer in addition to AlefSin's nice answer. –  Henrik Sep 15 '12 at 21:48
You're right Henrik - I didn't think of that when putting together this mock dataset - individuals should only be one species! –  Luke Sep 16 '12 at 0:58

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")))

-
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)