Comparaciones múltiples en un modelo de efectos mixtos

Estoy tratando de analizar algunos datos usando un modelo de efectos mixtos. Los datos que recopilé representan el peso de algunos animales jóvenes de diferentes genotipos 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

Entonces 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 multcomppuedo hacer:

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 uso mcppara ver la interacción entre el tiempo y el genotipo?

2. Cuando corro glhtme sale esta advertencia:

   covariate interactions found -- default contrast might be inappropriate

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

EDITAR: Encontré este PDF que dice:

Debido a que es imposible determinar los parámetros de interés automáticamente en este caso, mcp () en multcomp generará por defecto comparaciones solo para los efectos principales, ignorando las covariables y las interacciones . Desde la versión 1.1-2, se puede especificar promediar los términos de interacción y las covariables usando los argumentos interaccion_average = TRUE y covariate_average = TRUE respectivamente, mientras que las versiones anteriores a 1.0-0 promediaron automáticamente sobre los términos de interacción. Sin embargo, sugerimos a los usuarios que escriban, manualmente, el conjunto de contrastes que desean.Se debe hacer esto siempre que haya dudas acerca de lo que miden los contrastes predeterminados, lo que generalmente ocurre en modelos con términos de interacción de orden superior. Nos referimos a Hsu (1996), Capítulo ~ 7, y 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 aquí tiene?

Si timey Genotypeson predictores categóricos como parecen ser, y usted está interesado en comparar todos los pares de tiempo / genotipo entre sí, entonces puede crear una variable de interacción y usar 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, puede usar el hecho de que el linfctargumento puede tomar una matriz de coeficientes para los contrastes, de esta manera puede configurar exactamente las comparaciones que desee.

EDITAR

Parece haber cierta preocupación en los comentarios de que el modelo equipado con el TimeGenopredictor es diferente del modelo original equipado con el Time * Genotypepredictor. Este no es el caso , los modelos son equivalentes. La única diferencia está en la parametrización de los efectos fijos, que se configura para facilitar el uso de la glhtfunción.

He utilizado uno de los conjuntos de datos integrados (tiene Dieta en lugar de Genotipo) para demostrar que los dos enfoques tienen la misma probabilidad, valores pronosticados, 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 esto. Por otro lado, 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, es solo más trabajo configurarlas.