¿Cómo ejecutar un ANOVA de dos vías con una variable aleatoria seguido de comparaciones por pares?

Question

Tengo un conjunto de datos con el que me gustaría comparar el efecto de la especie y el hábitat en la tasa de movimiento, con comparaciones de pares. También me gustaría incluir el efecto del individuo (¿como un factor aleatorio?) - este factor aleatorio es la parte que no sé cómo hacer, al menos no en el marco de Anova().

Aquí hay un subconjunto de los datos:

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)

Actualmente, simplemente estoy ejecutando un ANOVA bidireccional con comparaciones de pares, sin considerar el efecto del individuo, de la siguiente manera:

fit <- lm(mvt.rate ~ habitat + species, data=data1)
require(car)
Anova(fit, type="III")
require(agricolae)
    #comparación de hábitats
comparison.hab <- HSD.test(fit, "habitat", group=TRUE)                     
    #comparación de especies
comparison.sp <- HSD.test(fit, "species", group=TRUE)

En el conjunto de datos, cada fila representa un movimiento, y en muchos casos, los individuos realizan varios movimientos (no independientes) - actualmente no estoy considerando esta falta de independencia de mvt.rate e individuo. Creo que la forma correcta de hacer esto es considerar al individuo como una variable aleatoria, pero no estoy completamente seguro.

Answer 1

Puedo mostrarte cómo modelar a los individuos como un efecto aleatorio. Simplemente puedes usar la función lme del paquete nlme. La sintaxis será muy similar a lm:

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

En cuanto a las comparaciones por pares, para lme puedes usar la función glht:

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

Answer 2

Solo elimine la variabilidad que está relacionada con los individuos: pre.fit=lm(mvt.rate ~ind,data=data1) y luego fit<-lm(pre.fit$residuals ~ data1$habitat +data1$species)