Múltiples comparaciones en un modelo de efectos mixtos

Question

Estoy tratando de analizar algunos datos usando un modelo de efectos mixtos. Los datos que he recogido representan el peso de algunos animales jóvenes de diferente genotipo a lo largo del tiempo.

Estoy usando el enfoque propuesto aquí: https://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/

En particular estoy usando la solución #2

Así que tengo algo como

require(nlme)
model <- lme(weight ~ time * Genotype, random = ~1|Animal/time, 
         data=weights)    
av <- anova(model)

Ahora, me gustaría tener algunas comparaciones múltiples. Usando multcomp Puedo hacerlo:

require(multcomp)
comp.geno <- glht(model, linfct=mcp(Genotype="Tukey"))
print(summary(comp.geno))

Y, por supuesto, podría hacer lo mismo con el tiempo.

Tengo dos preguntas:

1. ¿Cómo utilizo mcp para ver la interacción entre el Tiempo y el Genotipo?

2. Cuando corro glht Recibo esta advertencia:

   covariate interactions found -- default contrast might be inappropriate

   ¿Qué significa? ¿Puedo ignorarlo con seguridad? ¿O qué debo hacer para evitarlo?

EDITAR: Encontré este PDF que dice:

Debido a que en este caso es imposible determinar automáticamente los parámetros de interés, mcp() en multcomp generará por defecto comparaciones sólo para los efectos principales, ignorando las covariables y las interacciones . Desde la versión 1.1-2, se puede especificar que se promedien los términos de interacción y las covariables utilizando los argumentos interaction_average = TRUE y covariate_average = TRUE respectivamente, mientras que las versiones anteriores a la 1.0-0 promediaban automáticamente los términos de interacción. Sin embargo, sugerimos a los usuarios que escriban, manualmente, el conjunto de contrastes que deseen. Se debe hacer esto siempre que haya dudas sobre la medida de los contrastes por defecto, lo que típicamente ocurre en los modelos con términos de interacción de orden superior. Nos remitimos a Hsu (1996), Capítulo~7, y a Searle (1971), Capítulo~7.3, para más discusiones y ejemplos sobre este tema.

No tengo acceso a esos libros, pero tal vez alguien de aquí sí.

Answer 1

Si time y Genotype son ambos predictores categóricos como parecen ser, y estás interesado en comparar todos los pares tiempo/genotipo entre sí, entonces puedes crear una sola variable de interacción, y usar los contrastes de Tukey en ella:

weights$TimeGeno <- interaction(weigths$Time, weights$Geno)
model <- lme(weight ~ TimeGeno, random = ~1|Animal/time, data=weights) 
comp.timegeno <- glht(model, linfct=mcp(TimeGeno="Tukey")) 

Si está interesado en otros contrastes, entonces puede usar el hecho de que el linfct puede tomar una matriz de coeficientes para los contrastes - de esta manera puede establecer exactamente las comparaciones que desee.

EDITAR

Parece que hay cierta preocupación en los comentarios de que el modelo equipado con el TimeGeno predictor es diferente del modelo original equipado con el Time * Genotype predictor. No es así. los modelos son equivalentes. La única diferencia está en la parametrización de los efectos fijos, que se establece para facilitar el uso de la glht función.

He utilizado uno de los conjuntos de datos incorporados (tiene Dieta en lugar de Genotipo) para demostrar que los dos enfoques tienen la misma probabilidad, valores predichos, etc:

> # extract a subset of a built-in dataset for the example
> data(BodyWeight)
> ex <- as.data.frame(subset(BodyWeight, Time %in% c(1, 22, 44)))
> ex$Time <- factor(ex$Time)
> 
> #create interaction variable
> ex$TimeDiet <- interaction(ex$Time, ex$Diet)
    > 
    > model1 <- lme(weight ~ Time * Diet, random = ~1|Rat/Time,  data=ex)    
    > model2 <- lme(weight ~ TimeDiet, random = ~1|Rat/Time, data=ex)    
    > 
    > # the degrees of freedom, AIC, BIC, log-likelihood are all the same 
    > anova(model1, model2)
           Model df      AIC      BIC    logLik
    model1     1 12 367.4266 387.3893 -171.7133
    model2     2 12 367.4266 387.3893 -171.7133
    Warning message:
    In anova.lme(model1, model2) :
      fitted objects with different fixed effects. REML comparisons are not meaningful.
    > 
    > # the second model collapses the main and interaction effects of the first model
    > anova(model1)
                numDF denDF   F-value p-value
    (Intercept)     1    26 1719.5059  <.0001
    Time            2    26   28.9986  <.0001
    Diet            2    13   85.3659  <.0001
    Time:Diet       4    26    1.7610  0.1671
    > anova(model2)
                numDF denDF   F-value p-value
    (Intercept)     1    24 1719.5059  <.0001
    TimeDiet        8    24   29.4716  <.0001
    > 
    > # they give the same predicted values
    > newdata <- expand.grid(Time=levels(ex$Time), Diet=levels(ex$Diet))
    > newdata$TimeDiet <- interaction(newdata$Time, newdata$Diet)
> newdata$pred1 <- predict(model1, newdata=newdata, level=0)
    > newdata$pred2 <- predict(model2, newdata=newdata, level=0)
> newdata
  Time Diet TimeDiet   pred1   pred2
1    1    1      1.1 250.625 250.625
2   22    1     22.1 261.875 261.875
3   44    1     44.1 267.250 267.250
4    1    2      1.2 453.750 453.750
5   22    2     22.2 475.000 475.000
6   44    2     44.2 488.750 488.750
7    1    3      1.3 508.750 508.750
8   22    3     22.3 518.250 518.250
9   44    3     44.3 530.000 530.000

La única diferencia es que las hipótesis son fáciles de probar. Por ejemplo, en el primer modelo es fácil probar si los dos predictores interactúan, en el segundo modelo no hay una prueba explícita para ello. Por otra parte, el efecto conjunto de los dos predictores es fácil de probar en el segundo modelo, pero no en el primero. Las otras hipótesis son comprobables, sólo que es más difícil establecerlas.