Modelo mixto de comparaciones múltiples para la interacción entre predictores continuos y categóricos

11

Me gustaría usar lme4para ajustar una regresión de efectos mixtos y multcompcalcular las comparaciones por pares. Tengo un conjunto de datos complejo con múltiples predictores continuos y categóricos, pero mi pregunta se puede demostrar utilizando el ChickWeightconjunto de datos incorporado como ejemplo:

m <- lmer(weight ~ Time * Diet + (1 | Chick), data=ChickWeight, REML=F)

Timees continuo y Dietes categórico (4 niveles) y hay múltiples pollitos por dieta. Todos los polluelos comenzaron con aproximadamente el mismo peso, pero sus dietas (pueden) afectar su tasa de crecimiento, por lo que las Dietintersecciones deben ser (más o menos) iguales, pero las pendientes pueden ser diferentes. Puedo obtener las comparaciones por pares para el efecto de intercepción de Dietesta manera:

summary(glht(m, linfct=mcp(Diet = "Tukey")))

y, de hecho, no son significativamente diferentes, pero ¿cómo puedo hacer una prueba análoga para el Time:Dietefecto? Simplemente poner el término de interacción en mcpproduce un error:

summary(glht(m, linfct=mcp('Time:Diet' = "Tukey")))
Error in summary(glht(m, linfct = mcp(`Time:Diet` = "Tukey"))) : 
  error in evaluating the argument 'object' in selecting a method for function
 'summary': Error in mcp2matrix(model, linfct = linfct) : 
Variable(s) ‘Time:Diet’ have been specified in ‘linfct’ but cannot be found in ‘model’!
r mixed-model multiple-comparisons Dan M.
fuente
Tiene Time*Diet, que es solo una simplificación de Time + Diet + Time:Diet. Usando anova(m)o summary(m)confirma que el término de interacción está en el modelo.
Dan M.

Respuestas:

8

Por defecto, lmertrata el nivel de referencia de un predictor categórico como la línea de base y estima los parámetros para los otros niveles. Por lo tanto, obtiene algunas comparaciones por pares en la salida predeterminada y puede obtener las demás utilizando relevelpara definir un nuevo nivel de referencia y volver a ajustar el modelo. Esto tiene la ventaja de permitirle usar comparaciones de modelos o MCMC para obtener los valores p, pero no corrige las comparaciones múltiples (aunque podría aplicar su propia corrección después).

Para usar multcomp, debe definir la matriz de contraste. Cada fila en la matriz de contraste representa los pesos de los efectos que obtiene en la salida del modelo predeterminado, comenzando con la Intercepción. Entonces, si desea un efecto que ya esté incluido en la salida básica, simplemente ponga un "1" en la posición correspondiente a ese efecto. Dado que las estimaciones de los parámetros son relativas a un nivel de referencia común, puede obtener comparaciones entre cualquiera de los otros dos niveles estableciendo el peso de uno en "-1" y del otro "1". Así es como funcionaría para los Time:Diettérminos en el ChickWeightejemplo:

contrast.matrix <- rbind("Time:Diet1 vs. Time:Diet2" =  c(0, 0, 0, 0, 0, 1, 0, 0),
                           "Time:Diet1 vs. Time:Diet3" =  c(0, 0, 0, 0, 0, 0, 1, 0),
                           "Time:Diet1 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, 0, 0, 1),
                           "Time:Diet2 vs. Time:Diet3" =  c(0, 0, 0, 0, 0, -1, 1, 0),
                           "Time:Diet2 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, -1, 0, 1),
                           "Time:Diet3 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, 0, -1, 1))
summary(glht(m, contrast.matrix))

Advertencia: este enfoque parece utilizar la aproximación normal para obtener valores de p, que es algo anticonservador, y luego aplica alguna corrección para comparaciones múltiples. El resultado es que este método le brinda acceso fácil a tantas estimaciones de parámetros por pares y errores estándar como desee, pero los valores p pueden o no ser lo que desea.

(Gracias a Scott Jackson de r-ling-lang-L por su ayuda con esto)

Dan M.
fuente